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) 2011-2016 OpenFOAM Foundation
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 "fvcDiv.H"
28 #include "fvMatrices.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().deltaTValue();
62  scalar deltaT0 = mesh().time().deltaT0Value();
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.primitiveField()
85 
86  - (coefft*VV0 + coefft00*V0V00)
87  *vf.oldTime().primitiveField()
88 
89  + (coefft00*V0V00)*vf.oldTime().oldTime().primitiveField()
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().deltaTValue();
140  scalar deltaT0 = mesh().time().deltaT0Value();
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  const scalarField VV0rhoRho0
151  (
152  (mesh().V() + mesh().V0())
153  * (rho.primitiveField() + rho.oldTime().primitiveField())
154  );
155 
156  const scalarField V0V00rho0Rho00
157  (
158  (mesh().V0() + mesh().V00())
159  * (
160  rho.oldTime().primitiveField()
161  + rho.oldTime().oldTime().primitiveField()
162  )
163  );
164 
166  (
168  (
169  d2dt2IOobject,
170  mesh(),
171  rDeltaT2.dimensions()*rho.dimensions()*vf.dimensions(),
172  quarterRdeltaT2*
173  (
174  coefft*VV0rhoRho0*vf.primitiveField()
175 
176  - (coefft*VV0rhoRho0 + coefft00*V0V00rho0Rho00)
177  *vf.oldTime().primitiveField()
178 
179  + (coefft00*V0V00rho0Rho00)
180  *vf.oldTime().oldTime().primitiveField()
181  )/mesh().V(),
182  halfRdeltaT2*
183  (
184  coefft
185  *(rho.boundaryField() + rho.oldTime().boundaryField())
186  *vf.boundaryField()
187 
188  - (
189  coefft
190  *(
191  rho.boundaryField()
192  + rho.oldTime().boundaryField()
193  )
194  + coefft00
195  *(
196  rho.oldTime().boundaryField()
197  + rho.oldTime().oldTime().boundaryField()
198  )
199  )*vf.oldTime().boundaryField()
200 
201  + coefft00
202  *(
203  rho.oldTime().boundaryField()
204  + rho.oldTime().oldTime().boundaryField()
205  )*vf.oldTime().oldTime().boundaryField()
206  )
207  )
208  );
209  }
210  else
211  {
212  dimensionedScalar halfRdeltaT2 = 0.5*rDeltaT2;
213 
214  const volScalarField rhoRho0(rho + rho.oldTime());
215  const volScalarField rho0Rho00(rho.oldTime() +rho.oldTime().oldTime());
216 
218  (
220  (
221  d2dt2IOobject,
222  halfRdeltaT2*
223  (
224  coefft*rhoRho0*vf
225  - (coefft*rhoRho0 + coefft00*rho0Rho00)*vf.oldTime()
226  + coefft00*rho0Rho00*vf.oldTime().oldTime()
227  )
228  )
229  );
230  }
231 }
232 
233 
234 template<class Type>
237 (
239 )
240 {
241  tmp<fvMatrix<Type>> tfvm
242  (
243  new fvMatrix<Type>
244  (
245  vf,
247  )
248  );
249 
250  fvMatrix<Type>& fvm = tfvm.ref();
251 
252  scalar deltaT = mesh().time().deltaTValue();
253  scalar deltaT0 = mesh().time().deltaT0Value();
254 
255  scalar coefft = (deltaT + deltaT0)/(2*deltaT);
256  scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
257  scalar coefft0 = coefft + coefft00;
258 
259  scalar rDeltaT2 = 4.0/sqr(deltaT + deltaT0);
260 
261  if (mesh().moving())
262  {
263  scalar halfRdeltaT2 = rDeltaT2/2.0;
264 
265  const scalarField VV0(mesh().V() + mesh().V0());
266  const scalarField V0V00(mesh().V0() + mesh().V00());
267 
268  fvm.diag() = (coefft*halfRdeltaT2)*VV0;
269 
270  fvm.source() = halfRdeltaT2*
271  (
272  (coefft*VV0 + coefft00*V0V00)
273  *vf.oldTime().primitiveField()
274 
275  - (coefft00*V0V00)*vf.oldTime().oldTime().primitiveField()
276  );
277  }
278  else
279  {
280  fvm.diag() = (coefft*rDeltaT2)*mesh().V();
281 
282  fvm.source() = rDeltaT2*mesh().V()*
283  (
284  coefft0*vf.oldTime().primitiveField()
285  - coefft00*vf.oldTime().oldTime().primitiveField()
286  );
287  }
288 
289  return tfvm;
290 }
291 
292 
293 template<class Type>
296 (
297  const dimensionedScalar& rho,
299 )
300 {
301  tmp<fvMatrix<Type>> tfvm
302  (
303  new fvMatrix<Type>
304  (
305  vf,
306  rho.dimensions()*vf.dimensions()*dimVol
308  )
309  );
310 
311  fvMatrix<Type>& fvm = tfvm.ref();
312 
313  scalar deltaT = mesh().time().deltaTValue();
314  scalar deltaT0 = mesh().time().deltaT0Value();
315 
316  scalar coefft = (deltaT + deltaT0)/(2*deltaT);
317  scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
318 
319  scalar rDeltaT2 = 4.0/sqr(deltaT + deltaT0);
320 
321  if (mesh().moving())
322  {
323  scalar halfRdeltaT2 = 0.5*rDeltaT2;
324 
325  const scalarField VV0(mesh().V() + mesh().V0());
326  const scalarField V0V00(mesh().V0() + mesh().V00());
327 
328  fvm.diag() = rho.value()*(coefft*halfRdeltaT2)*VV0;
329 
330  fvm.source() = halfRdeltaT2*rho.value()*
331  (
332  (coefft*VV0 + coefft00*V0V00)
333  *vf.oldTime().primitiveField()
334 
335  - (coefft00*V0V00)*vf.oldTime().oldTime().primitiveField()
336  );
337  }
338  else
339  {
340  fvm.diag() = (coefft*rDeltaT2)*mesh().V()*rho.value();
341 
342  fvm.source() = rDeltaT2*mesh().V()*rho.value()*
343  (
344  (coefft + coefft00)*vf.oldTime().primitiveField()
345  - coefft00*vf.oldTime().oldTime().primitiveField()
346  );
347  }
348 
349  return tfvm;
350 }
351 
352 
353 template<class Type>
356 (
357  const volScalarField& rho,
359 )
360 {
361  tmp<fvMatrix<Type>> tfvm
362  (
363  new fvMatrix<Type>
364  (
365  vf,
366  rho.dimensions()*vf.dimensions()*dimVol
368  )
369  );
370 
371  fvMatrix<Type>& fvm = tfvm.ref();
372 
373  scalar deltaT = mesh().time().deltaTValue();
374  scalar deltaT0 = mesh().time().deltaT0Value();
375 
376  scalar coefft = (deltaT + deltaT0)/(2*deltaT);
377  scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
378 
379  scalar rDeltaT2 = 4.0/sqr(deltaT + deltaT0);
380 
381  if (mesh().moving())
382  {
383  scalar quarterRdeltaT2 = 0.25*rDeltaT2;
384 
385  const scalarField VV0rhoRho0
386  (
387  (mesh().V() + mesh().V0())
388  *(rho.primitiveField() + rho.oldTime().primitiveField())
389  );
390 
391  const scalarField V0V00rho0Rho00
392  (
393  (mesh().V0() + mesh().V00())
394  *(
395  rho.oldTime().primitiveField()
396  + rho.oldTime().oldTime().primitiveField()
397  )
398  );
399 
400  fvm.diag() = (coefft*quarterRdeltaT2)*VV0rhoRho0;
401 
402  fvm.source() = quarterRdeltaT2*
403  (
404  (coefft*VV0rhoRho0 + coefft00*V0V00rho0Rho00)
405  *vf.oldTime().primitiveField()
406 
407  - (coefft00*V0V00rho0Rho00)
408  *vf.oldTime().oldTime().primitiveField()
409  );
410  }
411  else
412  {
413  scalar halfRdeltaT2 = 0.5*rDeltaT2;
414 
415  const scalarField rhoRho0
416  (
417  rho.primitiveField()
418  + rho.oldTime().primitiveField()
419  );
420 
421  const scalarField rho0Rho00
422  (
423  rho.oldTime().primitiveField()
424  + rho.oldTime().oldTime().primitiveField()
425  );
426 
427  fvm.diag() = (coefft*halfRdeltaT2)*mesh().V()*rhoRho0;
428 
429  fvm.source() = halfRdeltaT2*mesh().V()*
430  (
431  (coefft*rhoRho0 + coefft00*rho0Rho00)
432  *vf.oldTime().primitiveField()
433 
434  - (coefft00*rho0Rho00)
435  *vf.oldTime().oldTime().primitiveField()
436  );
437  }
438 
439  return tfvm;
440 }
441 
442 
443 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
444 
445 } // End namespace fv
446 
447 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
448 
449 } // End namespace Foam
450 
451 // ************************************************************************* //
const dimensionSet & dimensions() const
Return const reference to dimensions.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
const dimensionSet dimVol(dimVolume)
Definition: dimensionSets.H:59
const Internal::FieldType & primitiveField() const
Return a const-reference to the internal field.
const Type & value() const
Return const reference to value.
Generic GeometricField class.
tmp< fvMatrix< Type > > fvmD2dt2(const GeometricField< Type, fvPatchField, volMesh > &)
const GeometricField< Type, PatchField, GeoMesh > & oldTime() const
Return old time field.
dynamicFvMesh & mesh
const Boundary & boundaryField() const
Return const-reference to the boundary field.
labelList fv(nPoints)
A special matrix type and solver, designed for finite volume solutions of scalar equations. Face addressing is used to make all matrix assembly and solution loops vectorise.
Definition: fvPatchField.H:71
const dimensionSet & dimensions() const
Return dimensions.
Calculate the divergence of the given field.
Field< Type > & source()
Definition: fvMatrix.H:291
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:51
A special matrix type and solver, designed for finite volume solutions of scalar equations.
A class for managing temporary objects.
Definition: PtrList.H:54
scalarField & diag()
Definition: lduMatrix.C:183
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:174
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:91
tmp< GeometricField< Type, fvPatchField, volMesh > > fvcD2dt2(const GeometricField< Type, fvPatchField, volMesh > &)
const word & name() const
Return name.
Definition: IOobject.H:260
Namespace for OpenFOAM.