EulerFaD2dt2Scheme.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) 2017 Volkswagen AG
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 "EulerFaD2dt2Scheme.H"
29 #include "facDiv.H"
30 #include "faMatrices.H"
31 
32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
36 
37 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
38 
39 namespace fa
40 {
41 
42 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
43 
44 template<class Type>
45 scalar EulerFaD2dt2Scheme<Type>::deltaT_() const
46 {
47  return mesh().time().deltaT().value();
48 }
49 
50 
51 template<class Type>
52 scalar EulerFaD2dt2Scheme<Type>::deltaT0_() const
53 {
54  return mesh().time().deltaT0().value();
55 }
56 
57 
58 template<class Type>
59 tmp<GeometricField<Type, faPatchField, areaMesh>>
60 EulerFaD2dt2Scheme<Type>::facD2dt2
61 (
62  const dimensioned<Type> dt
63 )
64 {
65 
66  scalar deltaT = deltaT_();
67  scalar deltaT0 = deltaT0_();
68 
69  dimensionedScalar rDeltaT2 =
70  4.0/sqr(mesh().time().deltaT() + mesh().time().deltaT0());
71 
72  scalar coefft = (deltaT + deltaT0)/(2*deltaT);
73  scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
74 
75  IOobject d2dt2IOobject
76  (
77  "d2dt2("+dt.name()+')',
78  mesh()().time().timeName(),
79  mesh()(),
80  IOobject::NO_READ,
81  IOobject::NO_WRITE
82  );
83 
84  if (mesh().moving())
85  {
86  scalar halfRdeltaT2 = rDeltaT2.value()/2.0;
87 
88  scalarField SS0 = mesh().S() + mesh().S0();
89  scalarField S0S00 = mesh().S0() + mesh().S00();
90 
91  tmp<GeometricField<Type, faPatchField, areaMesh>> tdt2dt2
92  (
93  new GeometricField<Type, faPatchField, areaMesh>
94  (
95  d2dt2IOobject,
96  mesh(),
97  dimensioned<Type>(dt.dimensions()/dimTime/dimTime, Zero)
98  )
99  );
100 
101  tdt2dt2.ref().primitiveFieldRef() =
102  halfRdeltaT2*dt.value()
103  *(coefft*SS0 - (coefft*SS0 + coefft00*S0S00) + coefft00*S0S00)
104  /mesh().S();
105 
106  return tdt2dt2;
107  }
108  else
109  {
110  return tmp<GeometricField<Type, faPatchField, areaMesh>>
111  (
112  new GeometricField<Type, faPatchField, areaMesh>
113  (
114  d2dt2IOobject,
115  mesh(),
116  dimensioned<Type>(dt.dimensions()/dimTime/dimTime, Zero),
117  calculatedFaPatchField<Type>::typeName
118  )
119  );
120  }
121 }
122 
123 
124 template<class Type>
125 tmp<GeometricField<Type, faPatchField, areaMesh>>
126 EulerFaD2dt2Scheme<Type>::facD2dt2
127 (
128  const GeometricField<Type, faPatchField, areaMesh>& vf
129 )
130 {
131  dimensionedScalar rDeltaT2 =
132  4.0/sqr(mesh().time().deltaT() + mesh().time().deltaT0());
133 
134  IOobject d2dt2IOobject
135  (
136  "d2dt2("+vf.name()+')',
137  mesh()().time().timeName(),
138  mesh()(),
139  IOobject::NO_READ,
140  IOobject::NO_WRITE
141  );
142 
143  scalar deltaT = deltaT_();
144  scalar deltaT0 = deltaT0_();
145 
146  scalar coefft = (deltaT + deltaT0)/(2*deltaT);
147  scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
148  scalar coefft0 = coefft + coefft00;
149 
150  if (mesh().moving())
151  {
152  scalar halfRdeltaT2 = rDeltaT2.value()/2.0;
153 
154  scalarField SS0 = mesh().S() + mesh().S0();
155  scalarField S0S00 = mesh().S0() + mesh().S00();
156 
157  return tmp<GeometricField<Type, faPatchField, areaMesh>>
158  (
159  new GeometricField<Type, faPatchField, areaMesh>
160  (
161  d2dt2IOobject,
162  mesh(),
163  rDeltaT2.dimensions()*vf.dimensions(),
164  halfRdeltaT2*
165  (
166  coefft*SS0*vf.primitiveField()
167 
168  - (coefft*SS0 + coefft00*S0S00)
169  *vf.oldTime().primitiveField()
170 
171  + (coefft00*S0S00)*vf.oldTime().oldTime().primitiveField()
172  )/mesh().S(),
173  rDeltaT2.value()*
174  (
175  coefft*vf.boundaryField()
176  - coefft0*vf.oldTime().boundaryField()
177  + coefft00*vf.oldTime().oldTime().boundaryField()
178  )
179  )
180  );
181  }
182  else
183  {
184  return tmp<GeometricField<Type, faPatchField, areaMesh>>
185  (
186  new GeometricField<Type, faPatchField, areaMesh>
187  (
188  d2dt2IOobject,
189  rDeltaT2*
190  (
191  coefft*vf
192  - coefft0*vf.oldTime()
193  + coefft00*vf.oldTime().oldTime()
194  )
195  )
196  );
197  }
198 }
199 
200 
201 template<class Type>
202 tmp<GeometricField<Type, faPatchField, areaMesh>>
203 EulerFaD2dt2Scheme<Type>::facD2dt2
204 (
205  const dimensionedScalar& rho,
206  const GeometricField<Type, faPatchField, areaMesh>& vf
207 )
208 {
209  dimensionedScalar rDeltaT2 =
210  4.0/sqr(mesh().time().deltaT() + mesh().time().deltaT0());
211 
212  IOobject d2dt2IOobject
213  (
214  "d2dt2("+rho.name()+','+vf.name()+')',
215  mesh()().time().timeName(),
216  mesh()(),
217  IOobject::NO_READ,
218  IOobject::NO_WRITE
219  );
220 
221  scalar deltaT = deltaT_();
222  scalar deltaT0 = deltaT0_();
223 
224  scalar coefft = (deltaT + deltaT0)/(2*deltaT);
225  scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
226  scalar coefft0 = coefft + coefft00;
227 
228  if (mesh().moving())
229  {
230  scalar halfRdeltaT2 = 0.5*rDeltaT2.value();
231 
232  scalarField SS0rhoRho0 =
233  (mesh().S() + mesh().S0())
234  *rho.value();
235 
236  scalarField S0S00rho0Rho00 =
237  (mesh().S0() + mesh().S00())
238  *rho.value();
239 
240  return tmp<GeometricField<Type, faPatchField, areaMesh>>
241  (
242  new GeometricField<Type, faPatchField, areaMesh>
243  (
244  d2dt2IOobject,
245  mesh(),
246  rDeltaT2.dimensions()*rho.dimensions()*vf.dimensions(),
247  halfRdeltaT2*
248  (
249  coefft*SS0rhoRho0*vf.primitiveField()
250 
251  - (coefft*SS0rhoRho0 + coefft00*S0S00rho0Rho00)
252  *vf.oldTime().primitiveField()
253 
254  + (coefft00*S0S00rho0Rho00)
255  *vf.oldTime().oldTime().primitiveField()
256  )/mesh().S(),
257  rDeltaT2.value()*rho.value()*
258  (
259  coefft*vf.boundaryField()
260  - coefft0*vf.oldTime().boundaryField()
261  + coefft00*vf.oldTime().oldTime().boundaryField()
262  )
263  )
264  );
265  }
266  else
267  {
268 
269  return tmp<GeometricField<Type, faPatchField, areaMesh>>
270  (
271  new GeometricField<Type, faPatchField, areaMesh>
272  (
273  d2dt2IOobject,
274  rDeltaT2 * rho.value() *
275  (
276  coefft*vf
277  - coefft0*vf.oldTime()
278  + coefft00*vf.oldTime().oldTime()
279  )
280  )
281  );
282  }
283 }
284 
285 
286 template<class Type>
287 tmp<GeometricField<Type, faPatchField, areaMesh>>
288 EulerFaD2dt2Scheme<Type>::facD2dt2
289 (
290  const areaScalarField& rho,
291  const GeometricField<Type, faPatchField, areaMesh>& vf
292 )
293 {
294  dimensionedScalar rDeltaT2 =
295  4.0/sqr(mesh().time().deltaT() + mesh().time().deltaT0());
296 
297  IOobject d2dt2IOobject
298  (
299  "d2dt2("+rho.name()+','+vf.name()+')',
300  mesh()().time().timeName(),
301  mesh()(),
302  IOobject::NO_READ,
303  IOobject::NO_WRITE
304  );
305 
306  scalar deltaT = deltaT_();
307  scalar deltaT0 = deltaT0_();
308 
309  scalar coefft = (deltaT + deltaT0)/(2*deltaT);
310  scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
311 
312  if (mesh().moving())
313  {
314  scalar halfRdeltaT2 = 0.5*rDeltaT2.value();
315  scalar quarterRdeltaT2 = 0.25*rDeltaT2.value();
316 
317  scalarField SS0rhoRho0
318  (
319  (mesh().S() + mesh().S0())
320  *(rho.primitiveField() + rho.oldTime().primitiveField())
321  );
322 
323  scalarField S0S00rho0Rho00
324  (
325  (mesh().S0() + mesh().S00())
326  *(
327  rho.oldTime().primitiveField()
328  + rho.oldTime().oldTime().primitiveField()
329  )
330  );
331 
332  return tmp<GeometricField<Type, faPatchField, areaMesh>>
333  (
334  new GeometricField<Type, faPatchField, areaMesh>
335  (
336  d2dt2IOobject,
337  mesh(),
338  rDeltaT2.dimensions()*rho.dimensions()*vf.dimensions(),
339  quarterRdeltaT2*
340  (
341  coefft*SS0rhoRho0*vf.primitiveField()
342 
343  - (coefft*SS0rhoRho0 + coefft00*S0S00rho0Rho00)
344  *vf.oldTime().primitiveField()
345 
346  + (coefft00*S0S00rho0Rho00)
347  *vf.oldTime().oldTime().primitiveField()
348  )/mesh().S(),
349  halfRdeltaT2*
350  (
351  coefft
352  *(rho.boundaryField() + rho.oldTime().boundaryField())
353  *vf.boundaryField()
354 
355  - (
356  coefft
357  *(
358  rho.boundaryField()
359  + rho.oldTime().boundaryField()
360  )
361  + coefft00
362  *(
363  rho.oldTime().boundaryField()
364  + rho.oldTime().oldTime().boundaryField()
365  )
366  )*vf.oldTime().boundaryField()
367 
368  + coefft00
369  *(
370  rho.oldTime().boundaryField()
371  + rho.oldTime().oldTime().boundaryField()
372  )*vf.oldTime().oldTime().boundaryField()
373  )
374  )
375  );
376  }
377  else
378  {
379  dimensionedScalar halfRdeltaT2 = 0.5*rDeltaT2;
380 
381  areaScalarField rhoRho0(rho + rho.oldTime());
382  areaScalarField rho0Rho00(rho.oldTime() + rho.oldTime().oldTime());
383 
384  return tmp<GeometricField<Type, faPatchField, areaMesh>>
385  (
386  new GeometricField<Type, faPatchField, areaMesh>
387  (
388  d2dt2IOobject,
389  halfRdeltaT2*
390  (
391  coefft*rhoRho0*vf
392  - (coefft*rhoRho0 + coefft00*rho0Rho00)*vf.oldTime()
393  + coefft00*rho0Rho00*vf.oldTime().oldTime()
394  )
395  )
396  );
397  }
398 }
399 
400 
401 template<class Type>
402 tmp<faMatrix<Type>>
403 EulerFaD2dt2Scheme<Type>::famD2dt2
404 (
405  const GeometricField<Type, faPatchField, areaMesh>& vf
406 )
407 {
408  tmp<faMatrix<Type>> tfam
409  (
410  new faMatrix<Type>
411  (
412  vf,
413  vf.dimensions()*dimArea/dimTime/dimTime
414  )
415  );
416 
417  faMatrix<Type>& fam = tfam.ref();
418 
419  scalar deltaT = deltaT_();
420  scalar deltaT0 = deltaT0_();
421 
422  scalar coefft = (deltaT + deltaT0)/(2*deltaT);
423  scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
424  scalar coefft0 = coefft + coefft00;
425 
426  scalar rDeltaT2 = 4.0/sqr(deltaT + deltaT0);
427 
428  if (mesh().moving())
429  {
430  scalar halfRdeltaT2 = rDeltaT2/2.0;
431 
432  scalarField SS0(mesh().S() + mesh().S0());
433  scalarField S0S00(mesh().S0() + mesh().S00());
434 
435  fam.diag() = (coefft*halfRdeltaT2)*SS0;
436 
437  fam.source() = halfRdeltaT2*
438  (
439  (coefft*SS0 + coefft00*S0S00)
440  *vf.oldTime().primitiveField()
441 
442  - (coefft00*S0S00)*vf.oldTime().oldTime().primitiveField()
443  );
444  }
445  else
446  {
447  fam.diag() = (coefft*rDeltaT2)*mesh().S();
448 
449  fam.source() = rDeltaT2*mesh().S()*
450  (
451  coefft0*vf.oldTime().primitiveField()
452  - coefft00*vf.oldTime().oldTime().primitiveField()
453  );
454  }
455 
456  return tfam;
457 }
458 
459 
460 template<class Type>
461 tmp<faMatrix<Type>>
462 EulerFaD2dt2Scheme<Type>::famD2dt2
463 (
464  const dimensionedScalar& rho,
465  const GeometricField<Type, faPatchField, areaMesh>& vf
466 )
467 {
468  tmp<faMatrix<Type>> tfam
469  (
470  new faMatrix<Type>
471  (
472  vf,
473  rho.dimensions()*vf.dimensions()*dimArea
474  /dimTime/dimTime
475  )
476  );
477 
478  faMatrix<Type>& fam = tfam.ref();
479 
480  scalar deltaT = deltaT_();
481  scalar deltaT0 = deltaT0_();
482 
483  scalar coefft = (deltaT + deltaT0)/(2*deltaT);
484  scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
485 
486  scalar rDeltaT2 = 4.0/sqr(deltaT + deltaT0);
487 
488  if (mesh().moving())
489  {
490  scalar halfRdeltaT2 = 0.5*rDeltaT2;
491 
492  scalarField SS0(mesh().S() + mesh().S0());
493 
494  scalarField S0S00(mesh().S0() + mesh().S00());
495 
496  fam.diag() = rho.value()*(coefft*halfRdeltaT2)*SS0;
497 
498  fam.source() = halfRdeltaT2*rho.value()*
499  (
500  (coefft*SS0 + coefft00*S0S00)
501  *vf.oldTime().primitiveField()
502 
503  - (coefft00*S0S00)*vf.oldTime().oldTime().primitiveField()
504  );
505  }
506  else
507  {
508  fam.diag() = (coefft*rDeltaT2)*mesh().S()*rho.value();
509 
510  fam.source() = rDeltaT2*mesh().S()*rho.value()*
511  (
512  (coefft + coefft00)*vf.oldTime().primitiveField()
513  - coefft00*vf.oldTime().oldTime().primitiveField()
514  );
515  }
516 
517  return tfam;
518 }
519 
520 
521 template<class Type>
522 tmp<faMatrix<Type>>
523 EulerFaD2dt2Scheme<Type>::famD2dt2
524 (
525  const areaScalarField& rho,
526  const GeometricField<Type, faPatchField, areaMesh>& vf
527 )
528 {
529  tmp<faMatrix<Type>> tfam
530  (
531  new faMatrix<Type>
532  (
533  vf,
534  rho.dimensions()*vf.dimensions()*dimArea/dimTime/dimTime
535  )
536  );
537  faMatrix<Type>& fam = tfam.ref();
538 
539  scalar deltaT =deltaT_();
540  scalar deltaT0 = deltaT0_();
541 
542  scalar coefft = (deltaT + deltaT0)/(2*deltaT);
543  scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
544 
545  scalar rDeltaT2 = 4.0/sqr(deltaT + deltaT0);
546 
547  if (mesh().moving())
548  {
549  scalar quarterRdeltaT2 = 0.25*rDeltaT2;
550 
551  scalarField SS0rhoRho0
552  (
553  (mesh().S() + mesh().S0())
554  *(rho.primitiveField() + rho.oldTime().primitiveField())
555  );
556 
557  scalarField S0S00rho0Rho00
558  (
559  (mesh().S0() + mesh().S00())
560  *(
561  rho.oldTime().primitiveField()
562  + rho.oldTime().oldTime().primitiveField()
563  )
564  );
565 
566  fam.diag() = (coefft*quarterRdeltaT2)*SS0rhoRho0;
567 
568  fam.source() = quarterRdeltaT2*
569  (
570  (coefft*SS0rhoRho0 + coefft00*S0S00rho0Rho00)
571  *vf.oldTime().primitiveField()
572 
573  - (coefft00*S0S00rho0Rho00)
574  *vf.oldTime().oldTime().primitiveField()
575  );
576  }
577  else
578  {
579  scalar halfRdeltaT2 = 0.5*rDeltaT2;
580 
581  scalarField rhoRho0
582  (
583  rho.primitiveField() + rho.oldTime().primitiveField()
584  );
585 
586  scalarField rho0Rho00
587  (
588  rho.oldTime().primitiveField()
589  + rho.oldTime().oldTime().primitiveField()
590  );
591 
592  fam.diag() = (coefft*halfRdeltaT2)*mesh().S()*rhoRho0;
593 
594  fam.source() = halfRdeltaT2*mesh().S()*
595  (
596  (coefft*rhoRho0 + coefft00*rho0Rho00)
597  *vf.oldTime().primitiveField()
598 
599  - (coefft00*rho0Rho00)
600  *vf.oldTime().oldTime().primitiveField()
601  );
602  }
603 
604  return tfam;
605 }
606 
607 
608 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
609 
610 } // End namespace fa
611 
612 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
613 
614 } // End namespace Foam
615 
616 // ************************************************************************* //
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::faMatrix
A special matrix type and solver, designed for finite area solutions of scalar equations....
Definition: faMatricesFwd.H:43
Foam::tmp
A class for managing temporary objects.
Definition: PtrList.H:61
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
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
Foam::dimensioned::name
const word & name() const
Return const reference to name.
Definition: dimensionedType.C:406
faMatrices.H
facDiv.H
Calculate the divergence of the given field.
Foam::dimensioned::value
const Type & value() const
Return const reference to value.
Definition: dimensionedType.C:434
Foam::fa::EulerFaD2dt2Scheme::facD2dt2
tmp< GeometricField< Type, faPatchField, areaMesh > > facD2dt2(const dimensioned< Type >)
Definition: EulerFaD2dt2Scheme.C:61
rho
rho
Definition: readInitialConditions.H:88
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
Foam::dimArea
const dimensionSet dimArea(sqr(dimLength))
Definition: dimensionSets.H:59
Foam::Field< scalar >
Foam::calculatedFaPatchField
Author Zeljko Tukovic, FMENA Hrvoje Jasak, Wikki Ltd.
Definition: calculatedFaPatchField.H:54
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:6
Foam::dimensioned
Generic dimensioned Type class.
Definition: dimensionedScalarFwd.H:42
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.C:33
Foam::sqr
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Definition: dimensionedSymmTensor.C:51
Foam::fa::EulerFaD2dt2Scheme::famD2dt2
tmp< faMatrix< Type > > famD2dt2(const GeometricField< Type, faPatchField, areaMesh > &)
Definition: EulerFaD2dt2Scheme.C:404
EulerFaD2dt2Scheme.H
fam
Calculate the matrix for the second temporal derivative.
Foam::dimensioned::dimensions
const dimensionSet & dimensions() const
Return const reference to dimensions.
Definition: dimensionedType.C:420
Foam::GeometricField
Generic GeometricField class.
Definition: areaFieldsFwd.H:53
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