EulerD2dt2Scheme.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2019 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  const dimensionedScalar rDeltaT2
50  (
51  4.0/sqr(mesh().time().deltaT() + mesh().time().deltaT0())
52  );
53 
54  const word d2dt2name("d2dt2("+vf.name()+')');
55 
56  const scalar deltaT = mesh().time().deltaTValue();
57  const scalar deltaT0 = mesh().time().deltaT0Value();
58 
59  const scalar coefft = (deltaT + deltaT0)/(2*deltaT);
60  const scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
61  const scalar coefft0 = coefft + coefft00;
62 
63  if (mesh().moving())
64  {
65  const dimensionedScalar halfRdeltaT2(rDeltaT2/2.0);
66 
67  const volScalarField::Internal VV0(mesh().V() + mesh().V0());
68  const volScalarField::Internal V0V00(mesh().V0() + mesh().V00());
69 
71  (
72  d2dt2name,
73  halfRdeltaT2*
74  (
75  coefft*VV0*vf()
76 
77  - (coefft*VV0 + coefft00*V0V00)
78  *vf.oldTime()()
79 
80  + (coefft00*V0V00)*vf.oldTime().oldTime()()
81  )/mesh().V(),
82  rDeltaT2.value()*
83  (
84  coefft*vf.boundaryField()
85  - coefft0*vf.oldTime().boundaryField()
86  + coefft00*vf.oldTime().oldTime().boundaryField()
87  )
88  );
89  }
90  else
91  {
93  (
94  d2dt2name,
95  rDeltaT2*
96  (
97  coefft*vf
98  - coefft0*vf.oldTime()
99  + coefft00*vf.oldTime().oldTime()
100  )
101  );
102  }
103 }
104 
105 
106 template<class Type>
109 (
110  const volScalarField& rho,
112 )
113 {
114  const dimensionedScalar rDeltaT2
115  (
116  4.0/sqr(mesh().time().deltaT() + mesh().time().deltaT0())
117  );
118 
119  const word d2dt2name("d2dt2("+rho.name()+','+vf.name()+')');
120 
121  const scalar deltaT = mesh().time().deltaTValue();
122  const scalar deltaT0 = mesh().time().deltaT0Value();
123 
124  const scalar coefft = (deltaT + deltaT0)/(2*deltaT);
125  const scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
126 
127  if (mesh().moving())
128  {
129  const dimensionedScalar halfRdeltaT2(0.5*rDeltaT2);
130  const dimensionedScalar quarterRdeltaT2(0.25*rDeltaT2);
131 
132  const volScalarField::Internal VV0rhoRho0
133  (
134  (mesh().V() + mesh().V0())*(rho() + rho.oldTime()())
135  );
136 
137  const volScalarField::Internal V0V00rho0Rho00
138  (
139  (mesh().V0() + mesh().V00())
140  *(
141  rho.oldTime()()
142  + rho.oldTime().oldTime()()
143  )
144  );
145 
147  (
148  d2dt2name,
149  quarterRdeltaT2*
150  (
151  coefft*VV0rhoRho0*vf()
152 
153  - (coefft*VV0rhoRho0 + coefft00*V0V00rho0Rho00)
154  *vf.oldTime()()
155 
156  + (coefft00*V0V00rho0Rho00)
157  *vf.oldTime().oldTime()()
158  )/mesh().V(),
159  halfRdeltaT2.value()*
160  (
161  coefft
162  *(rho.boundaryField() + rho.oldTime().boundaryField())
163  *vf.boundaryField()
164 
165  - (
166  coefft
167  *(
168  rho.boundaryField()
169  + rho.oldTime().boundaryField()
170  )
171  + coefft00
172  *(
173  rho.oldTime().boundaryField()
174  + rho.oldTime().oldTime().boundaryField()
175  )
176  )*vf.oldTime().boundaryField()
177 
178  + coefft00
179  *(
180  rho.oldTime().boundaryField()
181  + rho.oldTime().oldTime().boundaryField()
182  )*vf.oldTime().oldTime().boundaryField()
183  )
184  );
185  }
186  else
187  {
188  const dimensionedScalar halfRdeltaT2 = 0.5*rDeltaT2;
189 
190  const volScalarField rhoRho0(rho + rho.oldTime());
191  const volScalarField rho0Rho00(rho.oldTime() +rho.oldTime().oldTime());
192 
194  (
195  d2dt2name,
196  halfRdeltaT2*
197  (
198  coefft*rhoRho0*vf
199  - (coefft*rhoRho0 + coefft00*rho0Rho00)*vf.oldTime()
200  + coefft00*rho0Rho00*vf.oldTime().oldTime()
201  )
202  );
203  }
204 }
205 
206 
207 template<class Type>
210 (
212 )
213 {
214  tmp<fvMatrix<Type>> tfvm
215  (
216  new fvMatrix<Type>
217  (
218  vf,
220  )
221  );
222 
223  fvMatrix<Type>& fvm = tfvm.ref();
224 
225  const scalar deltaT = mesh().time().deltaTValue();
226  const scalar deltaT0 = mesh().time().deltaT0Value();
227 
228  const scalar coefft = (deltaT + deltaT0)/(2*deltaT);
229  const scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
230  const scalar coefft0 = coefft + coefft00;
231 
232  const scalar rDeltaT2 = 4.0/sqr(deltaT + deltaT0);
233 
234  if (mesh().moving())
235  {
236  const scalar halfRdeltaT2 = rDeltaT2/2.0;
237 
238  const scalarField VV0(mesh().V() + mesh().V0());
239  const scalarField V0V00(mesh().V0() + mesh().V00());
240 
241  fvm.diag() = (coefft*halfRdeltaT2)*VV0;
242 
243  fvm.source() = halfRdeltaT2*
244  (
245  (coefft*VV0 + coefft00*V0V00)
246  *vf.oldTime().primitiveField()
247 
248  - (coefft00*V0V00)*vf.oldTime().oldTime().primitiveField()
249  );
250  }
251  else
252  {
253  fvm.diag() = (coefft*rDeltaT2)*mesh().V();
254 
255  fvm.source() = rDeltaT2*mesh().V()*
256  (
257  coefft0*vf.oldTime().primitiveField()
258  - coefft00*vf.oldTime().oldTime().primitiveField()
259  );
260  }
261 
262  return tfvm;
263 }
264 
265 
266 template<class Type>
269 (
270  const dimensionedScalar& rho,
272 )
273 {
274  tmp<fvMatrix<Type>> tfvm
275  (
276  new fvMatrix<Type>
277  (
278  vf,
279  rho.dimensions()*vf.dimensions()*dimVol
281  )
282  );
283 
284  fvMatrix<Type>& fvm = tfvm.ref();
285 
286  const scalar deltaT = mesh().time().deltaTValue();
287  const scalar deltaT0 = mesh().time().deltaT0Value();
288 
289  const scalar coefft = (deltaT + deltaT0)/(2*deltaT);
290  const scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
291 
292  const scalar rDeltaT2 = 4.0/sqr(deltaT + deltaT0);
293 
294  if (mesh().moving())
295  {
296  const scalar halfRdeltaT2 = 0.5*rDeltaT2;
297 
298  const scalarField VV0(mesh().V() + mesh().V0());
299  const scalarField V0V00(mesh().V0() + mesh().V00());
300 
301  fvm.diag() = rho.value()*(coefft*halfRdeltaT2)*VV0;
302 
303  fvm.source() = halfRdeltaT2*rho.value()*
304  (
305  (coefft*VV0 + coefft00*V0V00)
306  *vf.oldTime().primitiveField()
307 
308  - (coefft00*V0V00)*vf.oldTime().oldTime().primitiveField()
309  );
310  }
311  else
312  {
313  fvm.diag() = (coefft*rDeltaT2)*mesh().V()*rho.value();
314 
315  fvm.source() = rDeltaT2*mesh().V()*rho.value()*
316  (
317  (coefft + coefft00)*vf.oldTime().primitiveField()
318  - coefft00*vf.oldTime().oldTime().primitiveField()
319  );
320  }
321 
322  return tfvm;
323 }
324 
325 
326 template<class Type>
329 (
330  const volScalarField& rho,
332 )
333 {
334  tmp<fvMatrix<Type>> tfvm
335  (
336  new fvMatrix<Type>
337  (
338  vf,
339  rho.dimensions()*vf.dimensions()*dimVol
341  )
342  );
343 
344  fvMatrix<Type>& fvm = tfvm.ref();
345 
346  const scalar deltaT = mesh().time().deltaTValue();
347  const scalar deltaT0 = mesh().time().deltaT0Value();
348 
349  const scalar coefft = (deltaT + deltaT0)/(2*deltaT);
350  const scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
351 
352  const scalar rDeltaT2 = 4.0/sqr(deltaT + deltaT0);
353 
354  if (mesh().moving())
355  {
356  const scalar quarterRdeltaT2 = 0.25*rDeltaT2;
357 
358  const scalarField VV0rhoRho0
359  (
360  (mesh().V() + mesh().V0())
361  *(rho.primitiveField() + rho.oldTime().primitiveField())
362  );
363 
364  const scalarField V0V00rho0Rho00
365  (
366  (mesh().V0() + mesh().V00())
367  *(
368  rho.oldTime().primitiveField()
369  + rho.oldTime().oldTime().primitiveField()
370  )
371  );
372 
373  fvm.diag() = (coefft*quarterRdeltaT2)*VV0rhoRho0;
374 
375  fvm.source() = quarterRdeltaT2*
376  (
377  (coefft*VV0rhoRho0 + coefft00*V0V00rho0Rho00)
378  *vf.oldTime().primitiveField()
379 
380  - (coefft00*V0V00rho0Rho00)
381  *vf.oldTime().oldTime().primitiveField()
382  );
383  }
384  else
385  {
386  const scalar halfRdeltaT2 = 0.5*rDeltaT2;
387 
388  const scalarField rhoRho0
389  (
390  rho.primitiveField()
391  + rho.oldTime().primitiveField()
392  );
393 
394  const scalarField rho0Rho00
395  (
396  rho.oldTime().primitiveField()
397  + rho.oldTime().oldTime().primitiveField()
398  );
399 
400  fvm.diag() = (coefft*halfRdeltaT2)*mesh().V()*rhoRho0;
401 
402  fvm.source() = halfRdeltaT2*mesh().V()*
403  (
404  (coefft*rhoRho0 + coefft00*rho0Rho00)
405  *vf.oldTime().primitiveField()
406 
407  - (coefft00*rho0Rho00)
408  *vf.oldTime().oldTime().primitiveField()
409  );
410  }
411 
412  return tfvm;
413 }
414 
415 
416 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
417 
418 } // End namespace fv
419 
420 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
421 
422 } // End namespace Foam
423 
424 // ************************************************************************* //
const GeometricField< Type, PatchField, GeoMesh > & oldTime() const
Return old time field.
const word & name() const
Return name.
Definition: IOobject.H:315
const Boundary & boundaryField() const
Return const-reference to the boundary field.
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:181
dimensionedSymmTensor sqr(const dimensionedVector &dv)
static tmp< GeometricField< Type, PatchField, GeoMesh > > New(const word &name, const Internal &, const PtrList< PatchField< Type >> &)
Return a temporary field constructed from name,.
const Internal::FieldType & primitiveField() const
Return a const-reference to the internal field.
Generic GeometricField class.
fvMesh & mesh
tmp< fvMatrix< Type > > fvmD2dt2(const GeometricField< Type, fvPatchField, volMesh > &)
const dimensionSet dimTime
const dimensionSet & dimensions() const
Return dimensions.
A class for handling words, derived from string.
Definition: word.H:59
labelList fv(nPoints)
const dimensionSet dimVol
const Type & value() const
Return const reference to value.
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:72
Calculate the divergence of the given field.
Field< Type > & source()
Definition: fvMatrix.H:300
Field with dimensions and associated with geometry type GeoMesh which is used to size the field and a...
A special matrix type and solver, designed for finite volume solutions of scalar equations.
const dimensionSet & dimensions() const
Return const reference to dimensions.
A class for managing temporary objects.
Definition: PtrList.H:53
scalarField & diag()
Definition: lduMatrix.C:186
tmp< GeometricField< Type, fvPatchField, volMesh > > fvcD2dt2(const GeometricField< Type, fvPatchField, volMesh > &)
Namespace for OpenFOAM.