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-------------------------------------------------------------------------------
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 "EulerD2dt2Scheme.H"
29#include "fvcDiv.H"
30#include "fvMatrices.H"
31
32// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
33
34namespace Foam
35{
36
37// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
38
39namespace fv
40{
41
42// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
43
44template<class Type>
45tmp<GeometricField<Type, fvPatchField, volMesh>>
46EulerD2dt2Scheme<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
121template<class Type>
122tmp<GeometricField<Type, fvPatchField, volMesh>>
123EulerD2dt2Scheme<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
236template<class Type>
237tmp<fvMatrix<Type>>
238EulerD2dt2Scheme<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
295template<class Type>
296tmp<fvMatrix<Type>>
297EulerD2dt2Scheme<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
355template<class Type>
356tmp<fvMatrix<Type>>
357EulerD2dt2Scheme<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::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::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::fvMatrix
A special matrix type and solver, designed for finite volume solutions of scalar equations....
Definition: fvMatrix.H:121
Foam::fvMatrix::source
Field< Type > & source() noexcept
Definition: fvMatrix.H:458
Foam::fv::EulerD2dt2Scheme::fvmD2dt2
tmp< fvMatrix< Type > > fvmD2dt2(const GeometricField< Type, fvPatchField, volMesh > &)
Definition: EulerD2dt2Scheme.C:239
Foam::fv::EulerD2dt2Scheme::fvcD2dt2
tmp< GeometricField< Type, fvPatchField, volMesh > > fvcD2dt2(const GeometricField< Type, fvPatchField, volMesh > &)
Definition: EulerD2dt2Scheme.C:47
Foam::lduMatrix::diag
scalarField & diag()
Definition: lduMatrix.C:192
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
fvMatrices.H
A special matrix type and solver, designed for finite volume solutions of scalar equations.
fvcDiv.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::dimVol
const dimensionSet dimVol(dimVolume)
Older spelling for dimVolume.
Definition: dimensionSets.H:64
fv
labelList fv(nPoints)