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 | Copyright (C) 2011 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> >
44 EulerD2dt2Scheme<Type>::fvcD2dt2
45 (
46  const GeometricField<Type, fvPatchField, volMesh>& vf
47 )
48 {
49  dimensionedScalar rDeltaT2 =
50  4.0/sqr(mesh().time().deltaT() + mesh().time().deltaT0());
51 
52  IOobject d2dt2IOobject
53  (
54  "d2dt2("+vf.name()+')',
55  mesh().time().timeName(),
56  mesh(),
57  IOobject::NO_READ,
58  IOobject::NO_WRITE
59  );
60 
61  scalar deltaT = mesh().time().deltaTValue();
62  scalar deltaT0 = mesh().time().deltaT0Value();
63 
64  scalar coefft = (deltaT + deltaT0)/(2*deltaT);
65  scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
66  scalar coefft0 = coefft + coefft00;
67 
68  if (mesh().moving())
69  {
70  scalar halfRdeltaT2 = rDeltaT2.value()/2.0;
71 
72  scalarField VV0 = mesh().V() + mesh().V0();
73  scalarField V0V00 = mesh().V0() + mesh().V00();
74 
75  return tmp<GeometricField<Type, fvPatchField, volMesh> >
76  (
77  new GeometricField<Type, fvPatchField, volMesh>
78  (
79  d2dt2IOobject,
80  mesh(),
81  rDeltaT2.dimensions()*vf.dimensions(),
82  halfRdeltaT2*
83  (
84  coefft*VV0*vf.internalField()
85 
86  - (coefft*VV0 + coefft00*V0V00)
87  *vf.oldTime().internalField()
88 
89  + (coefft00*V0V00)*vf.oldTime().oldTime().internalField()
90  )/mesh().V(),
91  rDeltaT2.value()*
92  (
93  coefft*vf.boundaryField()
94  - coefft0*vf.oldTime().boundaryField()
95  + coefft00*vf.oldTime().oldTime().boundaryField()
96  )
97  )
98  );
99  }
100  else
101  {
102  return tmp<GeometricField<Type, fvPatchField, volMesh> >
103  (
104  new GeometricField<Type, fvPatchField, volMesh>
105  (
106  d2dt2IOobject,
107  rDeltaT2*
108  (
109  coefft*vf
110  - coefft0*vf.oldTime()
111  + coefft00*vf.oldTime().oldTime()
112  )
113  )
114  );
115  }
116 }
117 
118 
119 template<class Type>
120 tmp<GeometricField<Type, fvPatchField, volMesh> >
121 EulerD2dt2Scheme<Type>::fvcD2dt2
122 (
123  const volScalarField& rho,
124  const GeometricField<Type, fvPatchField, volMesh>& vf
125 )
126 {
127  dimensionedScalar rDeltaT2 =
128  4.0/sqr(mesh().time().deltaT() + mesh().time().deltaT0());
129 
130  IOobject d2dt2IOobject
131  (
132  "d2dt2("+rho.name()+','+vf.name()+')',
133  mesh().time().timeName(),
134  mesh(),
135  IOobject::NO_READ,
136  IOobject::NO_WRITE
137  );
138 
139  scalar deltaT = mesh().time().deltaTValue();
140  scalar deltaT0 = mesh().time().deltaT0Value();
141 
142  scalar coefft = (deltaT + deltaT0)/(2*deltaT);
143  scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
144 
145  if (mesh().moving())
146  {
147  scalar halfRdeltaT2 = 0.5*rDeltaT2.value();
148  scalar quarterRdeltaT2 = 0.25*rDeltaT2.value();
149 
150  const scalarField VV0rhoRho0
151  (
152  (mesh().V() + mesh().V0())
153  * (rho.internalField() + rho.oldTime().internalField())
154  );
155 
156  const scalarField V0V00rho0Rho00
157  (
158  (mesh().V0() + mesh().V00())
159  * (
160  rho.oldTime().internalField()
161  + rho.oldTime().oldTime().internalField()
162  )
163  );
164 
165  return tmp<GeometricField<Type, fvPatchField, volMesh> >
166  (
167  new GeometricField<Type, fvPatchField, volMesh>
168  (
169  d2dt2IOobject,
170  mesh(),
171  rDeltaT2.dimensions()*rho.dimensions()*vf.dimensions(),
172  quarterRdeltaT2*
173  (
174  coefft*VV0rhoRho0*vf.internalField()
175 
176  - (coefft*VV0rhoRho0 + coefft00*V0V00rho0Rho00)
177  *vf.oldTime().internalField()
178 
179  + (coefft00*V0V00rho0Rho00)
180  *vf.oldTime().oldTime().internalField()
181  )/mesh().V(),
182  halfRdeltaT2*
183  (
184  coefft
185  *(rho.boundaryField() + rho.oldTime().boundaryField())
186  *vf.boundaryField()
187 
188  - (
189  coefft
190  *(
191  rho.boundaryField()
192  + rho.oldTime().boundaryField()
193  )
194  + coefft00
195  *(
196  rho.oldTime().boundaryField()
197  + rho.oldTime().oldTime().boundaryField()
198  )
199  )*vf.oldTime().boundaryField()
200 
201  + coefft00
202  *(
203  rho.oldTime().boundaryField()
204  + rho.oldTime().oldTime().boundaryField()
205  )*vf.oldTime().oldTime().boundaryField()
206  )
207  )
208  );
209  }
210  else
211  {
212  dimensionedScalar halfRdeltaT2 = 0.5*rDeltaT2;
213 
214  const volScalarField rhoRho0(rho + rho.oldTime());
215  const volScalarField rho0Rho00(rho.oldTime() +rho.oldTime().oldTime());
216 
217  return tmp<GeometricField<Type, fvPatchField, volMesh> >
218  (
219  new GeometricField<Type, fvPatchField, volMesh>
220  (
221  d2dt2IOobject,
222  halfRdeltaT2*
223  (
224  coefft*rhoRho0*vf
225  - (coefft*rhoRho0 + coefft00*rho0Rho00)*vf.oldTime()
226  + coefft00*rho0Rho00*vf.oldTime().oldTime()
227  )
228  )
229  );
230  }
231 }
232 
233 
234 template<class Type>
235 tmp<fvMatrix<Type> >
236 EulerD2dt2Scheme<Type>::fvmD2dt2
237 (
238  const GeometricField<Type, fvPatchField, volMesh>& vf
239 )
240 {
241  tmp<fvMatrix<Type> > tfvm
242  (
243  new fvMatrix<Type>
244  (
245  vf,
246  vf.dimensions()*dimVol/dimTime/dimTime
247  )
248  );
249 
250  fvMatrix<Type>& fvm = tfvm();
251 
252  scalar deltaT = mesh().time().deltaTValue();
253  scalar deltaT0 = mesh().time().deltaT0Value();
254 
255  scalar coefft = (deltaT + deltaT0)/(2*deltaT);
256  scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
257  scalar coefft0 = coefft + coefft00;
258 
259  scalar rDeltaT2 = 4.0/sqr(deltaT + deltaT0);
260 
261  if (mesh().moving())
262  {
263  scalar halfRdeltaT2 = rDeltaT2/2.0;
264 
265  const scalarField VV0(mesh().V() + mesh().V0());
266  const scalarField V0V00(mesh().V0() + mesh().V00());
267 
268  fvm.diag() = (coefft*halfRdeltaT2)*VV0;
269 
270  fvm.source() = halfRdeltaT2*
271  (
272  (coefft*VV0 + coefft00*V0V00)
273  *vf.oldTime().internalField()
274 
275  - (coefft00*V0V00)*vf.oldTime().oldTime().internalField()
276  );
277  }
278  else
279  {
280  fvm.diag() = (coefft*rDeltaT2)*mesh().V();
281 
282  fvm.source() = rDeltaT2*mesh().V()*
283  (
284  coefft0*vf.oldTime().internalField()
285  - coefft00*vf.oldTime().oldTime().internalField()
286  );
287  }
288 
289  return tfvm;
290 }
291 
292 
293 template<class Type>
294 tmp<fvMatrix<Type> >
295 EulerD2dt2Scheme<Type>::fvmD2dt2
296 (
297  const dimensionedScalar& rho,
298  const GeometricField<Type, fvPatchField, volMesh>& vf
299 )
300 {
301  tmp<fvMatrix<Type> > tfvm
302  (
303  new fvMatrix<Type>
304  (
305  vf,
306  rho.dimensions()*vf.dimensions()*dimVol
307  /dimTime/dimTime
308  )
309  );
310 
311  fvMatrix<Type>& fvm = tfvm();
312 
313  scalar deltaT = mesh().time().deltaTValue();
314  scalar deltaT0 = mesh().time().deltaT0Value();
315 
316  scalar coefft = (deltaT + deltaT0)/(2*deltaT);
317  scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
318 
319  scalar rDeltaT2 = 4.0/sqr(deltaT + deltaT0);
320 
321  if (mesh().moving())
322  {
323  scalar halfRdeltaT2 = 0.5*rDeltaT2;
324 
325  const scalarField VV0(mesh().V() + mesh().V0());
326  const scalarField V0V00(mesh().V0() + mesh().V00());
327 
328  fvm.diag() = rho.value()*(coefft*halfRdeltaT2)*VV0;
329 
330  fvm.source() = halfRdeltaT2*rho.value()*
331  (
332  (coefft*VV0 + coefft00*V0V00)
333  *vf.oldTime().internalField()
334 
335  - (coefft00*V0V00)*vf.oldTime().oldTime().internalField()
336  );
337  }
338  else
339  {
340  fvm.diag() = (coefft*rDeltaT2)*mesh().V()*rho.value();
341 
342  fvm.source() = rDeltaT2*mesh().V()*rho.value()*
343  (
344  (coefft + coefft00)*vf.oldTime().internalField()
345  - coefft00*vf.oldTime().oldTime().internalField()
346  );
347  }
348 
349  return tfvm;
350 }
351 
352 
353 template<class Type>
354 tmp<fvMatrix<Type> >
355 EulerD2dt2Scheme<Type>::fvmD2dt2
356 (
357  const volScalarField& rho,
358  const GeometricField<Type, fvPatchField, volMesh>& vf
359 )
360 {
361  tmp<fvMatrix<Type> > tfvm
362  (
363  new fvMatrix<Type>
364  (
365  vf,
366  rho.dimensions()*vf.dimensions()*dimVol
367  /dimTime/dimTime
368  )
369  );
370 
371  fvMatrix<Type>& fvm = tfvm();
372 
373  scalar deltaT = mesh().time().deltaTValue();
374  scalar deltaT0 = mesh().time().deltaT0Value();
375 
376  scalar coefft = (deltaT + deltaT0)/(2*deltaT);
377  scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
378 
379  scalar rDeltaT2 = 4.0/sqr(deltaT + deltaT0);
380 
381  if (mesh().moving())
382  {
383  scalar quarterRdeltaT2 = 0.25*rDeltaT2;
384 
385  const scalarField VV0rhoRho0
386  (
387  (mesh().V() + mesh().V0())
388  *(rho.internalField() + rho.oldTime().internalField())
389  );
390 
391  const scalarField V0V00rho0Rho00
392  (
393  (mesh().V0() + mesh().V00())
394  *(
395  rho.oldTime().internalField()
396  + rho.oldTime().oldTime().internalField()
397  )
398  );
399 
400  fvm.diag() = (coefft*quarterRdeltaT2)*VV0rhoRho0;
401 
402  fvm.source() = quarterRdeltaT2*
403  (
404  (coefft*VV0rhoRho0 + coefft00*V0V00rho0Rho00)
405  *vf.oldTime().internalField()
406 
407  - (coefft00*V0V00rho0Rho00)
408  *vf.oldTime().oldTime().internalField()
409  );
410  }
411  else
412  {
413  scalar halfRdeltaT2 = 0.5*rDeltaT2;
414 
415  const scalarField rhoRho0
416  (
417  rho.internalField()
418  + rho.oldTime().internalField()
419  );
420 
421  const scalarField rho0Rho00
422  (
423  rho.oldTime().internalField()
424  + rho.oldTime().oldTime().internalField()
425  );
426 
427  fvm.diag() = (coefft*halfRdeltaT2)*mesh().V()*rhoRho0;
428 
429  fvm.source() = halfRdeltaT2*mesh().V()*
430  (
431  (coefft*rhoRho0 + coefft00*rho0Rho00)
432  *vf.oldTime().internalField()
433 
434  - (coefft00*rho0Rho00)
435  *vf.oldTime().oldTime().internalField()
436  );
437  }
438 
439  return tfvm;
440 }
441 
442 
443 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
444 
445 } // End namespace fv
446 
447 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
448 
449 } // End namespace Foam
450 
451 // ************************************************************************* //
Foam::IOobject
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:91
Foam::tmp
A class for managing temporary objects.
Definition: PtrList.H:118
Foam::GeometricField::oldTime
const GeometricField< Type, PatchField, GeoMesh > & oldTime() const
Return old time field.
Definition: GeometricField.C:804
fvcDiv.H
Calculate the divergence of the given field.
Foam::GeometricField::boundaryField
GeometricBoundaryField & boundaryField()
Return reference to GeometricBoundaryField.
Definition: GeometricField.C:735
Foam::dimensioned::value
const Type & value() const
Return const reference to value.
Definition: dimensionedType.C:261
Foam::IOobject::NO_WRITE
@ NO_WRITE
Definition: IOobject.H:118
fvMatrices.H
A special matrix type and solver, designed for finite volume solutions of scalar equations.
Foam::dimTime
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:51
EulerD2dt2Scheme.H
Foam::Field
Pre-declare SubField and related Field type.
Definition: Field.H:57
Foam::IOobject::NO_READ
@ NO_READ
Definition: IOobject.H:111
Foam::GeometricField::internalField
InternalField & internalField()
Return internal field.
Definition: GeometricField.C:724
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:18
Foam::dimensioned
Generic dimensioned Type class.
Definition: dimensionedScalarFwd.H:41
Foam
Namespace for OpenFOAM.
Definition: combustionModel.C:30
rho
rho
Definition: pEqn.H:3
fv
labelList fv(nPoints)
Foam::sqr
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Definition: dimensionedSymmTensor.C:49
Foam::fvMatrix::source
Field< Type > & source()
Definition: fvMatrix.H:291
timeName
word timeName
Definition: getTimeIndex.H:3
Foam::fvMatrix< Type >
Foam::fv::EulerD2dt2Scheme::fvcD2dt2
tmp< GeometricField< Type, fvPatchField, volMesh > > fvcD2dt2(const GeometricField< Type, fvPatchField, volMesh > &)
Definition: EulerD2dt2Scheme.C:45
Foam::dimVol
const dimensionSet dimVol(dimVolume)
Definition: dimensionSets.H:59
Foam::dimensioned::dimensions
const dimensionSet & dimensions() const
Return const reference to dimensions.
Definition: dimensionedType.C:248
Foam::GeometricField
Generic GeometricField class.
Definition: surfaceFieldsFwd.H:52
Foam::fv::EulerD2dt2Scheme::fvmD2dt2
tmp< fvMatrix< Type > > fvmD2dt2(const GeometricField< Type, fvPatchField, volMesh > &)
Definition: EulerD2dt2Scheme.C:237

Copyright © 2011-2018  OpenFOAM | OPENFOAM® is a registered trademark of OpenCFD Ltd.

Generated by   doxygen 1.8.17