FreeFOAM The Cross-Platform CFD Toolkit
EulerD2dt2Scheme.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "EulerD2dt2Scheme.H"
27 #include <finiteVolume/fvcDiv.H>
29 
30 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
31 
32 namespace Foam
33 {
34 
35 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
36 
37 namespace fv
38 {
39 
40 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
41 
42 template<class Type>
43 tmp<GeometricField<Type, fvPatchField, volMesh> >
45 (
47 )
48 {
49  dimensionedScalar rDeltaT2 =
50  4.0/sqr(mesh().time().deltaT() + mesh().time().deltaT0());
51 
52  IOobject d2dt2IOobject
53  (
54  "d2dt2("+vf.name()+')',
55  mesh().time().timeName(),
56  mesh(),
59  );
60 
61  scalar deltaT = mesh().time().deltaT().value();
62  scalar deltaT0 = mesh().time().deltaT0().value();
63 
64  scalar coefft = (deltaT + deltaT0)/(2*deltaT);
65  scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
66  scalar coefft0 = coefft + coefft00;
67 
68  if (mesh().moving())
69  {
70  scalar halfRdeltaT2 = rDeltaT2.value()/2.0;
71 
72  scalarField VV0 = mesh().V() + mesh().V0();
73  scalarField V0V00 = mesh().V0() + mesh().V00();
74 
76  (
78  (
79  d2dt2IOobject,
80  mesh(),
81  rDeltaT2.dimensions()*vf.dimensions(),
82  halfRdeltaT2*
83  (
84  coefft*VV0*vf.internalField()
85 
86  - (coefft*VV0 + coefft00*V0V00)
87  *vf.oldTime().internalField()
88 
89  + (coefft00*V0V00)*vf.oldTime().oldTime().internalField()
90  )/mesh().V(),
91  rDeltaT2.value()*
92  (
93  coefft*vf.boundaryField()
94  - coefft0*vf.oldTime().boundaryField()
95  + coefft00*vf.oldTime().oldTime().boundaryField()
96  )
97  )
98  );
99  }
100  else
101  {
103  (
105  (
106  d2dt2IOobject,
107  rDeltaT2*
108  (
109  coefft*vf
110  - coefft0*vf.oldTime()
111  + coefft00*vf.oldTime().oldTime()
112  )
113  )
114  );
115  }
116 }
117 
118 
119 template<class Type>
122 (
123  const volScalarField& rho,
125 )
126 {
127  dimensionedScalar rDeltaT2 =
128  4.0/sqr(mesh().time().deltaT() + mesh().time().deltaT0());
129 
130  IOobject d2dt2IOobject
131  (
132  "d2dt2("+rho.name()+','+vf.name()+')',
133  mesh().time().timeName(),
134  mesh(),
137  );
138 
139  scalar deltaT = mesh().time().deltaT().value();
140  scalar deltaT0 = mesh().time().deltaT0().value();
141 
142  scalar coefft = (deltaT + deltaT0)/(2*deltaT);
143  scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
144 
145  if (mesh().moving())
146  {
147  scalar halfRdeltaT2 = 0.5*rDeltaT2.value();
148  scalar quarterRdeltaT2 = 0.25*rDeltaT2.value();
149 
150  scalarField VV0rhoRho0 =
151  (mesh().V() + mesh().V0())
152  *(rho.internalField() + rho.oldTime().internalField());
153 
154  scalarField V0V00rho0Rho00 =
155  (mesh().V0() + mesh().V00())
156  *(
157  rho.oldTime().internalField()
158  + rho.oldTime().oldTime().internalField()
159  );
160 
162  (
164  (
165  d2dt2IOobject,
166  mesh(),
167  rDeltaT2.dimensions()*rho.dimensions()*vf.dimensions(),
168  quarterRdeltaT2*
169  (
170  coefft*VV0rhoRho0*vf.internalField()
171 
172  - (coefft*VV0rhoRho0 + coefft00*V0V00rho0Rho00)
173  *vf.oldTime().internalField()
174 
175  + (coefft00*V0V00rho0Rho00)
176  *vf.oldTime().oldTime().internalField()
177  )/mesh().V(),
178  halfRdeltaT2*
179  (
180  coefft
181  *(rho.boundaryField() + rho.oldTime().boundaryField())
182  *vf.boundaryField()
183 
184  - (
185  coefft
186  *(
187  rho.boundaryField()
188  + rho.oldTime().boundaryField()
189  )
190  + coefft00
191  *(
192  rho.oldTime().boundaryField()
193  + rho.oldTime().oldTime().boundaryField()
194  )
195  )*vf.oldTime().boundaryField()
196 
197  + coefft00
198  *(
199  rho.oldTime().boundaryField()
200  + rho.oldTime().oldTime().boundaryField()
201  )*vf.oldTime().oldTime().boundaryField()
202  )
203  )
204  );
205  }
206  else
207  {
208  dimensionedScalar halfRdeltaT2 = 0.5*rDeltaT2;
209 
210  volScalarField rhoRho0 = rho + rho.oldTime();
211  volScalarField rho0Rho00 = rho.oldTime() +rho.oldTime().oldTime();
212 
214  (
216  (
217  d2dt2IOobject,
218  halfRdeltaT2*
219  (
220  coefft*rhoRho0*vf
221  - (coefft*rhoRho0 + coefft00*rho0Rho00)*vf.oldTime()
222  + coefft00*rho0Rho00*vf.oldTime().oldTime()
223  )
224  )
225  );
226  }
227 }
228 
229 
230 template<class Type>
233 (
235 )
236 {
237  tmp<fvMatrix<Type> > tfvm
238  (
239  new fvMatrix<Type>
240  (
241  vf,
243  )
244  );
245 
246  fvMatrix<Type>& fvm = tfvm();
247 
248  scalar deltaT = mesh().time().deltaT().value();
249  scalar deltaT0 = mesh().time().deltaT0().value();
250 
251  scalar coefft = (deltaT + deltaT0)/(2*deltaT);
252  scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
253  scalar coefft0 = coefft + coefft00;
254 
255  scalar rDeltaT2 = 4.0/sqr(deltaT + deltaT0);
256 
257  if (mesh().moving())
258  {
259  scalar halfRdeltaT2 = rDeltaT2/2.0;
260 
261  scalarField VV0 = mesh().V() + mesh().V0();
262  scalarField V0V00 = mesh().V0() + mesh().V00();
263 
264  fvm.diag() = (coefft*halfRdeltaT2)*VV0;
265 
266  fvm.source() = halfRdeltaT2*
267  (
268  (coefft*VV0 + coefft00*V0V00)
269  *vf.oldTime().internalField()
270 
271  - (coefft00*V0V00)*vf.oldTime().oldTime().internalField()
272  );
273  }
274  else
275  {
276  fvm.diag() = (coefft*rDeltaT2)*mesh().V();
277 
278  fvm.source() = rDeltaT2*mesh().V()*
279  (
280  coefft0*vf.oldTime().internalField()
281  - coefft00*vf.oldTime().oldTime().internalField()
282  );
283  }
284 
285  return tfvm;
286 }
287 
288 
289 template<class Type>
292 (
293  const dimensionedScalar& rho,
295 )
296 {
297  tmp<fvMatrix<Type> > tfvm
298  (
299  new fvMatrix<Type>
300  (
301  vf,
302  rho.dimensions()*vf.dimensions()*dimVol
304  )
305  );
306 
307  fvMatrix<Type>& fvm = tfvm();
308 
309  scalar deltaT = mesh().time().deltaT().value();
310  scalar deltaT0 = mesh().time().deltaT0().value();
311 
312  scalar coefft = (deltaT + deltaT0)/(2*deltaT);
313  scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
314 
315  scalar rDeltaT2 = 4.0/sqr(deltaT + deltaT0);
316 
317  if (mesh().moving())
318  {
319  scalar halfRdeltaT2 = 0.5*rDeltaT2;
320 
321  scalarField VV0 = mesh().V() + mesh().V0();
322 
323  scalarField V0V00 = mesh().V0() + mesh().V00();
324 
325  fvm.diag() = rho.value()*(coefft*halfRdeltaT2)*VV0;
326 
327  fvm.source() = halfRdeltaT2*rho.value()*
328  (
329  (coefft*VV0 + coefft00*V0V00)
330  *vf.oldTime().internalField()
331 
332  - (coefft00*V0V00)*vf.oldTime().oldTime().internalField()
333  );
334  }
335  else
336  {
337  fvm.diag() = (coefft*rDeltaT2)*mesh().V()*rho.value();
338 
339  fvm.source() = rDeltaT2*mesh().V()*rho.value()*
340  (
341  (coefft + coefft00)*vf.oldTime().internalField()
342  - coefft00*vf.oldTime().oldTime().internalField()
343  );
344  }
345 
346  return tfvm;
347 }
348 
349 
350 template<class Type>
353 (
354  const volScalarField& rho,
356 )
357 {
358  tmp<fvMatrix<Type> > tfvm
359  (
360  new fvMatrix<Type>
361  (
362  vf,
363  rho.dimensions()*vf.dimensions()*dimVol
365  )
366  );
367 
368  fvMatrix<Type>& fvm = tfvm();
369 
370  scalar deltaT = mesh().time().deltaT().value();
371  scalar deltaT0 = mesh().time().deltaT0().value();
372 
373  scalar coefft = (deltaT + deltaT0)/(2*deltaT);
374  scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
375 
376  scalar rDeltaT2 = 4.0/sqr(deltaT + deltaT0);
377 
378  if (mesh().moving())
379  {
380  scalar quarterRdeltaT2 = 0.25*rDeltaT2;
381 
382  scalarField VV0rhoRho0 =
383  (mesh().V() + mesh().V0())
384  *(rho.internalField() + rho.oldTime().internalField());
385 
386  scalarField V0V00rho0Rho00 =
387  (mesh().V0() + mesh().V00())
388  *(
389  rho.oldTime().internalField()
390  + rho.oldTime().oldTime().internalField()
391  );
392 
393  fvm.diag() = (coefft*quarterRdeltaT2)*VV0rhoRho0;
394 
395  fvm.source() = quarterRdeltaT2*
396  (
397  (coefft*VV0rhoRho0 + coefft00*V0V00rho0Rho00)
398  *vf.oldTime().internalField()
399 
400  - (coefft00*V0V00rho0Rho00)
401  *vf.oldTime().oldTime().internalField()
402  );
403  }
404  else
405  {
406  scalar halfRdeltaT2 = 0.5*rDeltaT2;
407 
408  scalarField rhoRho0 =
409  (rho.internalField() + rho.oldTime().internalField());
410 
411  scalarField rho0Rho00 =
412  (
413  rho.oldTime().internalField()
414  + rho.oldTime().oldTime().internalField()
415  );
416 
417  fvm.diag() = (coefft*halfRdeltaT2)*mesh().V()*rhoRho0;
418 
419  fvm.source() = halfRdeltaT2*mesh().V()*
420  (
421  (coefft*rhoRho0 + coefft00*rho0Rho00)
422  *vf.oldTime().internalField()
423 
424  - (coefft00*rho0Rho00)
425  *vf.oldTime().oldTime().internalField()
426  );
427  }
428 
429  return tfvm;
430 }
431 
432 
433 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
434 
435 } // End namespace fv
436 
437 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
438 
439 } // End namespace Foam
440 
441 // ************************ vim: set sw=4 sts=4 et: ************************ //