Loading...
Searching...
No Matches
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 "faMatrices.H"
30
31// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
32
33namespace Foam
34{
35
36// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
37
38namespace fa
39{
40
41// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
42
43template<class Type>
44scalar EulerFaD2dt2Scheme<Type>::deltaT_() const
45{
46 return mesh().time().deltaTValue();
47}
48
49
50template<class Type>
51scalar EulerFaD2dt2Scheme<Type>::deltaT0_() const
52{
53 return mesh().time().deltaT0Value();
54}
55
56
57template<class Type>
58tmp<GeometricField<Type, faPatchField, areaMesh>>
59EulerFaD2dt2Scheme<Type>::facD2dt2
60(
61 const dimensioned<Type> dt
62)
63{
64
65 scalar deltaT = deltaT_();
66 scalar deltaT0 = deltaT0_();
67
68 dimensionedScalar rDeltaT2 =
69 4.0/sqr(mesh().time().deltaT() + mesh().time().deltaT0());
70
71 scalar coefft = (deltaT + deltaT0)/(2*deltaT);
72 scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
73
74 const IOobject d2dt2IOobject
75 (
76 mesh().thisDb().newIOobject
77 (
78 "d2dt2("+dt.name()+')',
79 { IOobject::REGISTER }
80 )
81 );
82
83 if (mesh().moving())
84 {
85 scalar halfRdeltaT2 = 0.5*rDeltaT2.value();
86
87 scalarField SS0(mesh().S() + mesh().S0());
88 scalarField S0S00(mesh().S0() + mesh().S00());
89
90 tmp<GeometricField<Type, faPatchField, areaMesh>> tdt2dt2
91 (
92 new GeometricField<Type, faPatchField, areaMesh>
93 (
94 d2dt2IOobject,
95 mesh(),
96 dimensioned<Type>(dt.dimensions()/dimTime/dimTime, Zero)
97 )
98 );
99
100 tdt2dt2.ref().primitiveFieldRef() =
101 halfRdeltaT2*dt.value()
102 *(coefft*SS0 - (coefft*SS0 + coefft00*S0S00) + coefft00*S0S00)
103 /mesh().S();
104
105 return tdt2dt2;
106 }
107 else
108 {
109 return tmp<GeometricField<Type, faPatchField, areaMesh>>
110 (
111 new GeometricField<Type, faPatchField, areaMesh>
112 (
113 d2dt2IOobject,
114 mesh(),
115 dimensioned<Type>(dt.dimensions()/dimTime/dimTime, Zero),
116 faPatchFieldBase::calculatedType()
117 )
118 );
119 }
120}
121
122
123template<class Type>
124tmp<GeometricField<Type, faPatchField, areaMesh>>
125EulerFaD2dt2Scheme<Type>::facD2dt2
126(
127 const GeometricField<Type, faPatchField, areaMesh>& vf
128)
129{
130 dimensionedScalar rDeltaT2 =
131 4.0/sqr(mesh().time().deltaT() + mesh().time().deltaT0());
132
133 const IOobject d2dt2IOobject
134 (
135 mesh().thisDb().newIOobject
136 (
137 "d2dt2("+vf.name()+')',
138 { IOobject::REGISTER }
139 )
140 );
141
142 scalar deltaT = deltaT_();
143 scalar deltaT0 = deltaT0_();
144
145 scalar coefft = (deltaT + deltaT0)/(2*deltaT);
146 scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
147 scalar coefft0 = coefft + coefft00;
148
149 if (mesh().moving())
150 {
151 scalar halfRdeltaT2 = 0.5*rDeltaT2.value();
152
153 scalarField SS0(mesh().S() + mesh().S0());
154 scalarField S0S00(mesh().S0() + mesh().S00());
155
156 return tmp<GeometricField<Type, faPatchField, areaMesh>>
157 (
158 new GeometricField<Type, faPatchField, areaMesh>
159 (
160 d2dt2IOobject,
161 mesh(),
162 rDeltaT2.dimensions()*vf.dimensions(),
163 halfRdeltaT2*
164 (
165 coefft*SS0*vf.primitiveField()
166
167 - (coefft*SS0 + coefft00*S0S00)
168 *vf.oldTime().primitiveField()
169
170 + (coefft00*S0S00)*vf.oldTime().oldTime().primitiveField()
171 )/mesh().S(),
172 rDeltaT2.value()*
173 (
174 coefft*vf.boundaryField()
175 - coefft0*vf.oldTime().boundaryField()
176 + coefft00*vf.oldTime().oldTime().boundaryField()
177 )
178 )
179 );
180 }
181 else
182 {
183 return tmp<GeometricField<Type, faPatchField, areaMesh>>
184 (
185 new GeometricField<Type, faPatchField, areaMesh>
186 (
187 d2dt2IOobject,
188 rDeltaT2*
189 (
190 coefft*vf
191 - coefft0*vf.oldTime()
192 + coefft00*vf.oldTime().oldTime()
193 )
194 )
195 );
196 }
197}
198
199
200template<class Type>
201tmp<GeometricField<Type, faPatchField, areaMesh>>
202EulerFaD2dt2Scheme<Type>::facD2dt2
203(
204 const dimensionedScalar& rho,
205 const GeometricField<Type, faPatchField, areaMesh>& vf
206)
207{
208 dimensionedScalar rDeltaT2 =
209 4.0/sqr(mesh().time().deltaT() + mesh().time().deltaT0());
210
211 const IOobject d2dt2IOobject
212 (
213 mesh().thisDb().newIOobject
214 (
215 "d2dt2("+rho.name()+','+vf.name()+')',
216 { IOobject::REGISTER }
217 )
218 );
219
220 scalar deltaT = deltaT_();
221 scalar deltaT0 = deltaT0_();
222
223 scalar coefft = (deltaT + deltaT0)/(2*deltaT);
224 scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
225 scalar coefft0 = coefft + coefft00;
226
227 if (mesh().moving())
228 {
229 scalar halfRdeltaT2 = 0.5*rDeltaT2.value();
230
231 scalarField SS0rhoRho0
232 (
233 (mesh().S() + mesh().S0())*rho.value()
234 );
235
236 scalarField S0S00rho0Rho00
237 (
238 (mesh().S0() + mesh().S00())*rho.value()
239 );
240
241 return tmp<GeometricField<Type, faPatchField, areaMesh>>
242 (
243 new GeometricField<Type, faPatchField, areaMesh>
244 (
245 d2dt2IOobject,
246 mesh(),
247 rDeltaT2.dimensions()*rho.dimensions()*vf.dimensions(),
248 halfRdeltaT2*
249 (
250 coefft*SS0rhoRho0*vf.primitiveField()
251
252 - (coefft*SS0rhoRho0 + coefft00*S0S00rho0Rho00)
253 *vf.oldTime().primitiveField()
254
255 + (coefft00*S0S00rho0Rho00)
256 *vf.oldTime().oldTime().primitiveField()
257 )/mesh().S(),
258 rDeltaT2.value()*rho.value()*
259 (
260 coefft*vf.boundaryField()
261 - coefft0*vf.oldTime().boundaryField()
262 + coefft00*vf.oldTime().oldTime().boundaryField()
263 )
264 )
265 );
266 }
267 else
268 {
269
270 return tmp<GeometricField<Type, faPatchField, areaMesh>>
271 (
272 new GeometricField<Type, faPatchField, areaMesh>
273 (
274 d2dt2IOobject,
275 rDeltaT2 * rho.value() *
276 (
277 coefft*vf
278 - coefft0*vf.oldTime()
279 + coefft00*vf.oldTime().oldTime()
280 )
281 )
282 );
283 }
284}
285
286
287template<class Type>
288tmp<GeometricField<Type, faPatchField, areaMesh>>
289EulerFaD2dt2Scheme<Type>::facD2dt2
290(
291 const areaScalarField& rho,
292 const GeometricField<Type, faPatchField, areaMesh>& vf
293)
294{
295 dimensionedScalar rDeltaT2 =
296 4.0/sqr(mesh().time().deltaT() + mesh().time().deltaT0());
297
298 const IOobject d2dt2IOobject
299 (
300 mesh().thisDb().newIOobject
301 (
302 "d2dt2("+rho.name()+','+vf.name()+')',
303 { IOobject::REGISTER }
304 )
305 );
306
307 scalar deltaT = deltaT_();
308 scalar deltaT0 = deltaT0_();
309
310 scalar coefft = (deltaT + deltaT0)/(2*deltaT);
311 scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
312
313 if (mesh().moving())
314 {
315 scalar halfRdeltaT2 = 0.5*rDeltaT2.value();
316 scalar quarterRdeltaT2 = 0.25*rDeltaT2.value();
317
318 scalarField SS0rhoRho0
319 (
320 (mesh().S() + mesh().S0())
321 *(rho.primitiveField() + rho.oldTime().primitiveField())
322 );
323
324 scalarField S0S00rho0Rho00
325 (
326 (mesh().S0() + mesh().S00())
327 *(
328 rho.oldTime().primitiveField()
329 + rho.oldTime().oldTime().primitiveField()
330 )
331 );
332
333 return tmp<GeometricField<Type, faPatchField, areaMesh>>
334 (
335 new GeometricField<Type, faPatchField, areaMesh>
336 (
337 d2dt2IOobject,
338 mesh(),
339 rDeltaT2.dimensions()*rho.dimensions()*vf.dimensions(),
340 quarterRdeltaT2*
341 (
342 coefft*SS0rhoRho0*vf.primitiveField()
343
344 - (coefft*SS0rhoRho0 + coefft00*S0S00rho0Rho00)
345 *vf.oldTime().primitiveField()
346
347 + (coefft00*S0S00rho0Rho00)
348 *vf.oldTime().oldTime().primitiveField()
349 )/mesh().S(),
350 halfRdeltaT2*
351 (
352 coefft
353 *(rho.boundaryField() + rho.oldTime().boundaryField())
354 *vf.boundaryField()
355
356 - (
357 coefft
358 *(
359 rho.boundaryField()
360 + rho.oldTime().boundaryField()
361 )
362 + coefft00
363 *(
364 rho.oldTime().boundaryField()
365 + rho.oldTime().oldTime().boundaryField()
366 )
367 )*vf.oldTime().boundaryField()
368
369 + coefft00
370 *(
371 rho.oldTime().boundaryField()
372 + rho.oldTime().oldTime().boundaryField()
373 )*vf.oldTime().oldTime().boundaryField()
374 )
375 )
376 );
377 }
378 else
379 {
380 dimensionedScalar halfRdeltaT2 = 0.5*rDeltaT2;
381
382 areaScalarField rhoRho0(rho + rho.oldTime());
383 areaScalarField rho0Rho00(rho.oldTime() + rho.oldTime().oldTime());
384
385 return tmp<GeometricField<Type, faPatchField, areaMesh>>
386 (
387 new GeometricField<Type, faPatchField, areaMesh>
388 (
389 d2dt2IOobject,
390 halfRdeltaT2*
391 (
392 coefft*rhoRho0*vf
393 - (coefft*rhoRho0 + coefft00*rho0Rho00)*vf.oldTime()
394 + coefft00*rho0Rho00*vf.oldTime().oldTime()
395 )
396 )
397 );
398 }
399}
400
401
402template<class Type>
403tmp<faMatrix<Type>>
404EulerFaD2dt2Scheme<Type>::famD2dt2
405(
406 const GeometricField<Type, faPatchField, areaMesh>& vf
407)
408{
409 tmp<faMatrix<Type>> tfam
410 (
411 new faMatrix<Type>
412 (
413 vf,
414 vf.dimensions()*dimArea/dimTime/dimTime
415 )
416 );
417
418 faMatrix<Type>& fam = tfam.ref();
419
420 scalar deltaT = deltaT_();
421 scalar deltaT0 = deltaT0_();
422
423 scalar coefft = (deltaT + deltaT0)/(2*deltaT);
424 scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
425 scalar coefft0 = coefft + coefft00;
426
427 scalar rDeltaT2 = 4.0/sqr(deltaT + deltaT0);
428
429 if (mesh().moving())
430 {
431 scalar halfRdeltaT2 = rDeltaT2/2.0;
432
433 scalarField SS0(mesh().S() + mesh().S0());
434 scalarField S0S00(mesh().S0() + mesh().S00());
435
436 fam.diag() = (coefft*halfRdeltaT2)*SS0;
437
438 fam.source() = halfRdeltaT2*
439 (
440 (coefft*SS0 + coefft00*S0S00)
441 *vf.oldTime().primitiveField()
442
443 - (coefft00*S0S00)*vf.oldTime().oldTime().primitiveField()
444 );
445 }
446 else
447 {
448 fam.diag() = (coefft*rDeltaT2)*mesh().S();
449
450 fam.source() = rDeltaT2*mesh().S()*
451 (
452 coefft0*vf.oldTime().primitiveField()
453 - coefft00*vf.oldTime().oldTime().primitiveField()
454 );
455 }
456
457 return tfam;
458}
459
460
461template<class Type>
462tmp<faMatrix<Type>>
463EulerFaD2dt2Scheme<Type>::famD2dt2
464(
465 const dimensionedScalar& rho,
466 const GeometricField<Type, faPatchField, areaMesh>& vf
467)
468{
469 tmp<faMatrix<Type>> tfam
470 (
471 new faMatrix<Type>
472 (
473 vf,
474 rho.dimensions()*vf.dimensions()*dimArea
475 /dimTime/dimTime
476 )
477 );
478
479 faMatrix<Type>& fam = tfam.ref();
480
481 scalar deltaT = deltaT_();
482 scalar deltaT0 = deltaT0_();
483
484 scalar coefft = (deltaT + deltaT0)/(2*deltaT);
485 scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
486
487 scalar rDeltaT2 = 4.0/sqr(deltaT + deltaT0);
488
489 if (mesh().moving())
490 {
491 scalar halfRdeltaT2 = 0.5*rDeltaT2;
492
493 scalarField SS0(mesh().S() + mesh().S0());
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 noexcept
Return dimensions.
Definition DimensionedField.H:693
Foam::GeometricField
Generic GeometricField class.
Definition GeometricField.H:75
Foam::GeometricField::oldTime
const GeometricField< Type, PatchField, GeoMesh > & oldTime() const
Return old time field.
Definition GeometricField.C:961
Foam::GeometricField::primitiveField
const Internal::FieldType & primitiveField() const noexcept
Return a const-reference to the internal field values.
Definition GeometricFieldI.H:36
Foam::GeometricField::boundaryField
const Boundary & boundaryField() const noexcept
Return const-reference to the boundary field.
Definition GeometricFieldI.H:45
Foam::IOobject
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition IOobject.H:191
Foam::IOobject::name
const word & name() const noexcept
Return the object name.
Definition IOobjectI.H:205
Foam::dimensioned
Generic dimensioned Type class.
Definition dimensionedType.H:70
Foam::dimensioned::dimensions
const dimensionSet & dimensions() const noexcept
Return const reference to dimensions.
Definition dimensionedType.H:362
Foam::dimensioned::name
const word & name() const noexcept
Return const reference to name.
Definition dimensionedType.H:352
Foam::dimensioned::value
const Type & value() const noexcept
Return const reference to value.
Definition dimensionedType.H:372
Foam::faMatrix
A special matrix type and solver, designed for finite area solutions of scalar equations....
Definition faMatrix.H:108
Foam::faPatchFieldBase::calculatedType
static const word & calculatedType() noexcept
The type name for calculated patch fields.
Definition faPatchField.H:174
Foam::fa::EulerFaD2dt2Scheme::famD2dt2
tmp< faMatrix< Type > > famD2dt2(const GeometricField< Type, faPatchField, areaMesh > &)
Definition EulerFaD2dt2Scheme.C:398
Foam::fa::EulerFaD2dt2Scheme::mesh
const faMesh & mesh() const
Return mesh reference.
Definition EulerFaD2dt2Scheme.H:115
Foam::fa::EulerFaD2dt2Scheme::facD2dt2
tmp< GeometricField< Type, faPatchField, areaMesh > > facD2dt2(const dimensioned< Type >)
Definition EulerFaD2dt2Scheme.C:53
Foam::tmp
A class for managing temporary objects.
Definition tmp.H:75
Foam::tmp::ref
T & ref() const
Return non-const reference to the contents of a non-null managed pointer.
Definition tmpI.H:235
mesh
dynamicFvMesh & mesh
Definition createDynamicFvMesh.H:6
Foam::fa
Namespace for finite-area.
Definition limitHeight.C:30
Foam
Namespace for OpenFOAM.
Definition atmBoundaryLayer.C:27
Foam::sqr
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Definition dimensionedSymmTensor.C:44
Foam::dimTime
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition dimensionSets.H:51
Foam::dimArea
const dimensionSet dimArea(sqr(dimLength))
Definition dimensionSets.H:57
Foam::scalarField
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Definition primitiveFieldsFwd.H:46
Foam::areaScalarField
GeometricField< scalar, faPatchField, areaMesh > areaScalarField
Definition areaFieldsFwd.H:72
Foam::Zero
static constexpr const zero Zero
Global zero (0).
Definition zero.H:127
Foam::dimensionedScalar
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
Definition dimensionedScalarFwd.H:35
fam
Calculate the matrix for the second temporal derivative.