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 | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2011-2016 OpenFOAM Foundation
9 -------------------------------------------------------------------------------
10 License
11  This file is part of OpenFOAM.
12 
13  OpenFOAM is free software: you can redistribute it and/or modify it
14  under the terms of the GNU General Public License as published by
15  the Free Software Foundation, either version 3 of the License, or
16  (at your option) any later version.
17 
18  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21  for more details.
22 
23  You should have received a copy of the GNU General Public License
24  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25 
26 \*---------------------------------------------------------------------------*/
27 
28 #include "EulerD2dt2Scheme.H"
29 #include "fvcDiv.H"
30 #include "fvMatrices.H"
31 
32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
36 
37 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
38 
39 namespace fv
40 {
41 
42 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
43 
44 template<class Type>
45 tmp<GeometricField<Type, fvPatchField, volMesh>>
46 EulerD2dt2Scheme<Type>::fvcD2dt2
47 (
48  const GeometricField<Type, fvPatchField, volMesh>& vf
49 )
50 {
51  dimensionedScalar rDeltaT2 =
52  4.0/sqr(mesh().time().deltaT() + mesh().time().deltaT0());
53 
54  IOobject d2dt2IOobject
55  (
56  "d2dt2("+vf.name()+')',
57  mesh().time().timeName(),
58  mesh(),
59  IOobject::NO_READ,
60  IOobject::NO_WRITE
61  );
62 
63  scalar deltaT = mesh().time().deltaTValue();
64  scalar deltaT0 = mesh().time().deltaT0Value();
65 
66  scalar coefft = (deltaT + deltaT0)/(2*deltaT);
67  scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
68  scalar coefft0 = coefft + coefft00;
69 
70  if (mesh().moving())
71  {
72  scalar halfRdeltaT2 = rDeltaT2.value()/2.0;
73 
74  scalarField VV0 = mesh().V() + mesh().V0();
75  scalarField V0V00 = mesh().V0() + mesh().V00();
76 
77  return tmp<GeometricField<Type, fvPatchField, volMesh>>
78  (
79  new GeometricField<Type, fvPatchField, volMesh>
80  (
81  d2dt2IOobject,
82  mesh(),
83  rDeltaT2.dimensions()*vf.dimensions(),
84  halfRdeltaT2*
85  (
86  coefft*VV0*vf.primitiveField()
87 
88  - (coefft*VV0 + coefft00*V0V00)
89  *vf.oldTime().primitiveField()
90 
91  + (coefft00*V0V00)*vf.oldTime().oldTime().primitiveField()
92  )/mesh().V(),
93  rDeltaT2.value()*
94  (
95  coefft*vf.boundaryField()
96  - coefft0*vf.oldTime().boundaryField()
97  + coefft00*vf.oldTime().oldTime().boundaryField()
98  )
99  )
100  );
101  }
102  else
103  {
104  return tmp<GeometricField<Type, fvPatchField, volMesh>>
105  (
106  new GeometricField<Type, fvPatchField, volMesh>
107  (
108  d2dt2IOobject,
109  rDeltaT2*
110  (
111  coefft*vf
112  - coefft0*vf.oldTime()
113  + coefft00*vf.oldTime().oldTime()
114  )
115  )
116  );
117  }
118 }
119 
120 
121 template<class Type>
122 tmp<GeometricField<Type, fvPatchField, volMesh>>
123 EulerD2dt2Scheme<Type>::fvcD2dt2
124 (
125  const volScalarField& rho,
126  const GeometricField<Type, fvPatchField, volMesh>& vf
127 )
128 {
129  dimensionedScalar rDeltaT2 =
130  4.0/sqr(mesh().time().deltaT() + mesh().time().deltaT0());
131 
132  IOobject d2dt2IOobject
133  (
134  "d2dt2("+rho.name()+','+vf.name()+')',
135  mesh().time().timeName(),
136  mesh(),
137  IOobject::NO_READ,
138  IOobject::NO_WRITE
139  );
140 
141  scalar deltaT = mesh().time().deltaTValue();
142  scalar deltaT0 = mesh().time().deltaT0Value();
143 
144  scalar coefft = (deltaT + deltaT0)/(2*deltaT);
145  scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
146 
147  if (mesh().moving())
148  {
149  scalar halfRdeltaT2 = 0.5*rDeltaT2.value();
150  scalar quarterRdeltaT2 = 0.25*rDeltaT2.value();
151 
152  const scalarField VV0rhoRho0
153  (
154  (mesh().V() + mesh().V0())
155  * (rho.primitiveField() + rho.oldTime().primitiveField())
156  );
157 
158  const scalarField V0V00rho0Rho00
159  (
160  (mesh().V0() + mesh().V00())
161  * (
162  rho.oldTime().primitiveField()
163  + rho.oldTime().oldTime().primitiveField()
164  )
165  );
166 
167  return tmp<GeometricField<Type, fvPatchField, volMesh>>
168  (
169  new GeometricField<Type, fvPatchField, volMesh>
170  (
171  d2dt2IOobject,
172  mesh(),
173  rDeltaT2.dimensions()*rho.dimensions()*vf.dimensions(),
174  quarterRdeltaT2*
175  (
176  coefft*VV0rhoRho0*vf.primitiveField()
177 
178  - (coefft*VV0rhoRho0 + coefft00*V0V00rho0Rho00)
179  *vf.oldTime().primitiveField()
180 
181  + (coefft00*V0V00rho0Rho00)
182  *vf.oldTime().oldTime().primitiveField()
183  )/mesh().V(),
184  halfRdeltaT2*
185  (
186  coefft
187  *(rho.boundaryField() + rho.oldTime().boundaryField())
188  *vf.boundaryField()
189 
190  - (
191  coefft
192  *(
193  rho.boundaryField()
194  + rho.oldTime().boundaryField()
195  )
196  + coefft00
197  *(
198  rho.oldTime().boundaryField()
199  + rho.oldTime().oldTime().boundaryField()
200  )
201  )*vf.oldTime().boundaryField()
202 
203  + coefft00
204  *(
205  rho.oldTime().boundaryField()
206  + rho.oldTime().oldTime().boundaryField()
207  )*vf.oldTime().oldTime().boundaryField()
208  )
209  )
210  );
211  }
212  else
213  {
214  dimensionedScalar halfRdeltaT2 = 0.5*rDeltaT2;
215 
216  const volScalarField rhoRho0(rho + rho.oldTime());
217  const volScalarField rho0Rho00(rho.oldTime() +rho.oldTime().oldTime());
218 
219  return tmp<GeometricField<Type, fvPatchField, volMesh>>
220  (
221  new GeometricField<Type, fvPatchField, volMesh>
222  (
223  d2dt2IOobject,
224  halfRdeltaT2*
225  (
226  coefft*rhoRho0*vf
227  - (coefft*rhoRho0 + coefft00*rho0Rho00)*vf.oldTime()
228  + coefft00*rho0Rho00*vf.oldTime().oldTime()
229  )
230  )
231  );
232  }
233 }
234 
235 
236 template<class Type>
237 tmp<fvMatrix<Type>>
238 EulerD2dt2Scheme<Type>::fvmD2dt2
239 (
240  const GeometricField<Type, fvPatchField, volMesh>& vf
241 )
242 {
243  tmp<fvMatrix<Type>> tfvm
244  (
245  new fvMatrix<Type>
246  (
247  vf,
248  vf.dimensions()*dimVol/dimTime/dimTime
249  )
250  );
251 
252  fvMatrix<Type>& fvm = tfvm.ref();
253 
254  scalar deltaT = mesh().time().deltaTValue();
255  scalar deltaT0 = mesh().time().deltaT0Value();
256 
257  scalar coefft = (deltaT + deltaT0)/(2*deltaT);
258  scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
259  scalar coefft0 = coefft + coefft00;
260 
261  scalar rDeltaT2 = 4.0/sqr(deltaT + deltaT0);
262 
263  if (mesh().moving())
264  {
265  scalar halfRdeltaT2 = rDeltaT2/2.0;
266 
267  const scalarField VV0(mesh().V() + mesh().V0());
268  const scalarField V0V00(mesh().V0() + mesh().V00());
269 
270  fvm.diag() = (coefft*halfRdeltaT2)*VV0;
271 
272  fvm.source() = halfRdeltaT2*
273  (
274  (coefft*VV0 + coefft00*V0V00)
275  *vf.oldTime().primitiveField()
276 
277  - (coefft00*V0V00)*vf.oldTime().oldTime().primitiveField()
278  );
279  }
280  else
281  {
282  fvm.diag() = (coefft*rDeltaT2)*mesh().V();
283 
284  fvm.source() = rDeltaT2*mesh().V()*
285  (
286  coefft0*vf.oldTime().primitiveField()
287  - coefft00*vf.oldTime().oldTime().primitiveField()
288  );
289  }
290 
291  return tfvm;
292 }
293 
294 
295 template<class Type>
296 tmp<fvMatrix<Type>>
297 EulerD2dt2Scheme<Type>::fvmD2dt2
298 (
299  const dimensionedScalar& rho,
300  const GeometricField<Type, fvPatchField, volMesh>& vf
301 )
302 {
303  tmp<fvMatrix<Type>> tfvm
304  (
305  new fvMatrix<Type>
306  (
307  vf,
308  rho.dimensions()*vf.dimensions()*dimVol
309  /dimTime/dimTime
310  )
311  );
312 
313  fvMatrix<Type>& fvm = tfvm.ref();
314 
315  scalar deltaT = mesh().time().deltaTValue();
316  scalar deltaT0 = mesh().time().deltaT0Value();
317 
318  scalar coefft = (deltaT + deltaT0)/(2*deltaT);
319  scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
320 
321  scalar rDeltaT2 = 4.0/sqr(deltaT + deltaT0);
322 
323  if (mesh().moving())
324  {
325  scalar halfRdeltaT2 = 0.5*rDeltaT2;
326 
327  const scalarField VV0(mesh().V() + mesh().V0());
328  const scalarField V0V00(mesh().V0() + mesh().V00());
329 
330  fvm.diag() = rho.value()*(coefft*halfRdeltaT2)*VV0;
331 
332  fvm.source() = halfRdeltaT2*rho.value()*
333  (
334  (coefft*VV0 + coefft00*V0V00)
335  *vf.oldTime().primitiveField()
336 
337  - (coefft00*V0V00)*vf.oldTime().oldTime().primitiveField()
338  );
339  }
340  else
341  {
342  fvm.diag() = (coefft*rDeltaT2)*mesh().V()*rho.value();
343 
344  fvm.source() = rDeltaT2*mesh().V()*rho.value()*
345  (
346  (coefft + coefft00)*vf.oldTime().primitiveField()
347  - coefft00*vf.oldTime().oldTime().primitiveField()
348  );
349  }
350 
351  return tfvm;
352 }
353 
354 
355 template<class Type>
356 tmp<fvMatrix<Type>>
357 EulerD2dt2Scheme<Type>::fvmD2dt2
358 (
359  const volScalarField& rho,
360  const GeometricField<Type, fvPatchField, volMesh>& vf
361 )
362 {
363  tmp<fvMatrix<Type>> tfvm
364  (
365  new fvMatrix<Type>
366  (
367  vf,
368  rho.dimensions()*vf.dimensions()*dimVol
369  /dimTime/dimTime
370  )
371  );
372 
373  fvMatrix<Type>& fvm = tfvm.ref();
374 
375  scalar deltaT = mesh().time().deltaTValue();
376  scalar deltaT0 = mesh().time().deltaT0Value();
377 
378  scalar coefft = (deltaT + deltaT0)/(2*deltaT);
379  scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
380 
381  scalar rDeltaT2 = 4.0/sqr(deltaT + deltaT0);
382 
383  if (mesh().moving())
384  {
385  scalar quarterRdeltaT2 = 0.25*rDeltaT2;
386 
387  const scalarField VV0rhoRho0
388  (
389  (mesh().V() + mesh().V0())
390  *(rho.primitiveField() + rho.oldTime().primitiveField())
391  );
392 
393  const scalarField V0V00rho0Rho00
394  (
395  (mesh().V0() + mesh().V00())
396  *(
397  rho.oldTime().primitiveField()
398  + rho.oldTime().oldTime().primitiveField()
399  )
400  );
401 
402  fvm.diag() = (coefft*quarterRdeltaT2)*VV0rhoRho0;
403 
404  fvm.source() = quarterRdeltaT2*
405  (
406  (coefft*VV0rhoRho0 + coefft00*V0V00rho0Rho00)
407  *vf.oldTime().primitiveField()
408 
409  - (coefft00*V0V00rho0Rho00)
410  *vf.oldTime().oldTime().primitiveField()
411  );
412  }
413  else
414  {
415  scalar halfRdeltaT2 = 0.5*rDeltaT2;
416 
417  const scalarField rhoRho0
418  (
419  rho.primitiveField()
420  + rho.oldTime().primitiveField()
421  );
422 
423  const scalarField rho0Rho00
424  (
425  rho.oldTime().primitiveField()
426  + rho.oldTime().oldTime().primitiveField()
427  );
428 
429  fvm.diag() = (coefft*halfRdeltaT2)*mesh().V()*rhoRho0;
430 
431  fvm.source() = halfRdeltaT2*mesh().V()*
432  (
433  (coefft*rhoRho0 + coefft00*rho0Rho00)
434  *vf.oldTime().primitiveField()
435 
436  - (coefft00*rho0Rho00)
437  *vf.oldTime().oldTime().primitiveField()
438  );
439  }
440 
441  return tfvm;
442 }
443 
444 
445 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
446 
447 } // End namespace fv
448 
449 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
450 
451 } // End namespace Foam
452 
453 // ************************************************************************* //
Foam::IOobject::NO_WRITE
Definition: IOobject.H:195
Foam::IOobject
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:169
Foam::tmp
A class for managing temporary objects.
Definition: PtrList.H:61
Foam::GeometricField::primitiveField
const Internal::FieldType & primitiveField() const
Return a const-reference to the internal field.
Definition: GeometricFieldI.H:53
Foam::GeometricField::oldTime
const GeometricField< Type, PatchField, GeoMesh > & oldTime() const
Return old time field.
Definition: GeometricField.C:850
fvcDiv.H
Calculate the divergence of the given field.
Foam::dimensioned::value
const Type & value() const
Return const reference to value.
Definition: dimensionedType.C:434
rho
rho
Definition: readInitialConditions.H:88
fvMatrices.H
A special matrix type and solver, designed for finite volume solutions of scalar equations.
Foam::dimTime
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:53
Foam::tmp::ref
T & ref() const
Definition: tmpI.H:227
EulerD2dt2Scheme.H
Foam::Field< scalar >
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:6
Foam::dimensioned< scalar >
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.C:33
fv
labelList fv(nPoints)
Foam::sqr
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Definition: dimensionedSymmTensor.C:51
Foam::fvMatrix::source
Field< Type > & source()
Definition: fvMatrix.H:445
Foam::fvMatrix
A special matrix type and solver, designed for finite volume solutions of scalar equations....
Definition: fvPatchField.H:68
Foam::fv::EulerD2dt2Scheme::fvcD2dt2
tmp< GeometricField< Type, fvPatchField, volMesh > > fvcD2dt2(const GeometricField< Type, fvPatchField, volMesh > &)
Definition: EulerD2dt2Scheme.C:47
Foam::dimVol
const dimensionSet dimVol(dimVolume)
Older spelling for dimVolume.
Definition: dimensionSets.H:64
Foam::dimensioned::dimensions
const dimensionSet & dimensions() const
Return const reference to dimensions.
Definition: dimensionedType.C:420
Foam::GeometricField< Type, fvPatchField, volMesh >
Foam::IOobject::NO_READ
Definition: IOobject.H:188
Foam::GeometricField::boundaryField
const Boundary & boundaryField() const
Return const-reference to the boundary field.
Definition: GeometricFieldI.H:62
Foam::fv::EulerD2dt2Scheme::fvmD2dt2
tmp< fvMatrix< Type > > fvmD2dt2(const GeometricField< Type, fvPatchField, volMesh > &)
Definition: EulerD2dt2Scheme.C:239