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-------------------------------------------------------------------------------
10License
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
34namespace Foam
35{
36
37// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
38
39namespace fa
40{
41
42// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
43
44template<class Type>
45scalar EulerFaD2dt2Scheme<Type>::deltaT_() const
46{
47 return mesh().time().deltaT().value();
48}
49
50
51template<class Type>
52scalar EulerFaD2dt2Scheme<Type>::deltaT0_() const
53{
54 return mesh().time().deltaT0().value();
55}
56
57
58template<class Type>
59tmp<GeometricField<Type, faPatchField, areaMesh>>
60EulerFaD2dt2Scheme<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
124template<class Type>
125tmp<GeometricField<Type, faPatchField, areaMesh>>
126EulerFaD2dt2Scheme<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
201template<class Type>
202tmp<GeometricField<Type, faPatchField, areaMesh>>
203EulerFaD2dt2Scheme<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
286template<class Type>
287tmp<GeometricField<Type, faPatchField, areaMesh>>
288EulerFaD2dt2Scheme<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
401template<class Type>
402tmp<faMatrix<Type>>
403EulerFaD2dt2Scheme<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
460template<class Type>
461tmp<faMatrix<Type>>
462EulerFaD2dt2Scheme<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
521template<class Type>
522tmp<faMatrix<Type>>
523EulerFaD2dt2Scheme<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::DimensionedField::dimensions
const dimensionSet & dimensions() const
Return dimensions.
Definition: DimensionedFieldI.H:49
Foam::GeometricField
Generic GeometricField class.
Definition: GeometricField.H:80
Foam::GeometricField::oldTime
const GeometricField< Type, PatchField, GeoMesh > & oldTime() const
Return old time field.
Definition: GeometricField.C:978
Foam::GeometricField::primitiveField
const Internal::FieldType & primitiveField() const
Return a const-reference to the internal field.
Definition: GeometricFieldI.H:53
Foam::GeometricField::boundaryField
const Boundary & boundaryField() const
Return const-reference to the boundary field.
Definition: GeometricFieldI.H:62
Foam::IOobject
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:170
Foam::IOobject::name
const word & name() const noexcept
Return the object name.
Definition: IOobjectI.H:65
Foam::calculatedFaPatchField
Author Zeljko Tukovic, FMENA Hrvoje Jasak, Wikki Ltd.
Definition: calculatedFaPatchField.H:57
Foam::dimensioned
Generic dimensioned Type class.
Definition: dimensionedType.H:67
Foam::dimensioned::dimensions
const dimensionSet & dimensions() const
Return const reference to dimensions.
Definition: dimensionedType.C:420
Foam::dimensioned::value
const Type & value() const
Return const reference to value.
Definition: dimensionedType.C:434
Foam::dimensioned::name
const word & name() const
Return const reference to name.
Definition: dimensionedType.C:406
Foam::faMatrix
A special matrix type and solver, designed for finite area solutions of scalar equations....
Definition: faMatrix.H:76
Foam::fa::EulerFaD2dt2Scheme::famD2dt2
tmp< faMatrix< Type > > famD2dt2(const GeometricField< Type, faPatchField, areaMesh > &)
Definition: EulerFaD2dt2Scheme.C:404
Foam::fa::EulerFaD2dt2Scheme::facD2dt2
tmp< GeometricField< Type, faPatchField, areaMesh > > facD2dt2(const dimensioned< Type >)
Definition: EulerFaD2dt2Scheme.C:61
Foam::tmp
A class for managing temporary objects.
Definition: tmp.H:65
Foam::tmp::ref
T & ref() const
Definition: tmpI.H:227
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:6
facDiv.H
Calculate the divergence of the given field.
timeName
word timeName
Definition: getTimeIndex.H:3
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.C:34
Foam::sqr
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Definition: dimensionedSymmTensor.C:51
Foam::dimTime
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:53
Foam::dimArea
const dimensionSet dimArea(sqr(dimLength))
Definition: dimensionSets.H:59
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
fam
Calculate the matrix for the second temporal derivative.