backwardDdtScheme.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2024 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 "backwardDdtScheme.H"
27 #include "surfaceInterpolate.H"
28 #include "fvcDiv.H"
29 #include "fvMatrices.H"
30 
31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
32 
33 namespace Foam
34 {
35 
36 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
37 
38 namespace fv
39 {
40 
41 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
42 
43 template<class Type>
44 scalar backwardDdtScheme<Type>::deltaT_() const
45 {
46  return mesh().time().deltaTValue();
47 }
48 
49 
50 template<class Type>
51 scalar backwardDdtScheme<Type>::deltaT0_() const
52 {
53  return mesh().time().deltaT0Value();
54 }
55 
56 
57 template<class Type>
58 template<class GeoField>
59 scalar backwardDdtScheme<Type>::deltaT0_(const GeoField& vf) const
60 {
61  if (vf.nOldTimes() < 2)
62  {
63  return great;
64  }
65  else
66  {
67  return deltaT0_();
68  }
69 }
70 
71 
72 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
73 
74 template<class Type>
75 tmp<VolField<Type>>
76 backwardDdtScheme<Type>::fvcDdt
77 (
78  const dimensioned<Type>& dt
79 )
80 {
81  const dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
82 
83  const word ddtName("ddt("+dt.name()+')');
84 
85  const scalar deltaT = deltaT_();
86  const scalar deltaT0 = deltaT0_();
87 
88  const scalar coefft = 1 + deltaT/(deltaT + deltaT0);
89  const scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
90  const scalar coefft0 = coefft + coefft00;
91 
92  if (mesh().moving())
93  {
94  tmp<VolField<Type>> tdtdt
95  (
96  VolField<Type>::New
97  (
98  ddtName,
99  mesh(),
100  dimensioned<Type>
101  (
102  "0",
103  dt.dimensions()/dimTime,
104  Zero
105  )
106  )
107  );
108 
109  tdtdt.ref().primitiveFieldRef() = rDeltaT.value()*dt.value()*
110  (
111  coefft - (coefft0*mesh().V0() - coefft00*mesh().V00())/mesh().V()
112  );
113 
114  return tdtdt;
115  }
116  else
117  {
118  return tmp<VolField<Type>>
119  (
120  VolField<Type>::New
121  (
122  ddtName,
123  mesh(),
124  dimensioned<Type>
125  (
126  "0",
127  dt.dimensions()/dimTime,
128  Zero
129  ),
130  calculatedFvPatchField<Type>::typeName
131  )
132  );
133  }
134 }
135 
136 
137 template<class Type>
138 tmp<VolField<Type>>
139 backwardDdtScheme<Type>::fvcDdt
140 (
141  const VolField<Type>& vf
142 )
143 {
144  const dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
145 
146  const word ddtName("ddt("+vf.name()+')');
147 
148  const scalar deltaT = deltaT_();
149  const scalar deltaT0 = deltaT0_(vf);
150 
151  const scalar coefft = 1 + deltaT/(deltaT + deltaT0);
152  const scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
153  const scalar coefft0 = coefft + coefft00;
154 
155  if (mesh().moving())
156  {
157  return tmp<VolField<Type>>
158  (
159  VolField<Type>::New
160  (
161  ddtName,
162  rDeltaT*
163  (
164  coefft*vf() -
165  (
166  coefft0*vf.oldTime()()*mesh().V0()
167  - coefft00*vf.oldTime().oldTime()()
168  *mesh().V00()
169  )/mesh().V()
170  ),
171  rDeltaT.value()*
172  (
173  coefft*vf.boundaryField() -
174  (
175  coefft0*vf.oldTime().boundaryField()
176  - coefft00*vf.oldTime().oldTime().boundaryField()
177  )
178  )
179  )
180  );
181  }
182  else
183  {
184  return tmp<VolField<Type>>
185  (
186  VolField<Type>::New
187  (
188  ddtName,
189  rDeltaT*
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<VolField<Type>>
203 backwardDdtScheme<Type>::fvcDdt
204 (
205  const dimensionedScalar& rho,
206  const VolField<Type>& vf
207 )
208 {
209  const dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
210 
211  const word ddtName("ddt("+rho.name()+','+vf.name()+')');
212 
213  const scalar deltaT = deltaT_();
214  const scalar deltaT0 = deltaT0_(vf);
215 
216  const scalar coefft = 1 + deltaT/(deltaT + deltaT0);
217  const scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
218  const scalar coefft0 = coefft + coefft00;
219 
220  if (mesh().moving())
221  {
222  return tmp<VolField<Type>>
223  (
224  VolField<Type>::New
225  (
226  ddtName,
227  rDeltaT*rho*
228  (
229  coefft*vf() -
230  (
231  coefft0*vf.oldTime()()*mesh().V0()
232  - coefft00*vf.oldTime().oldTime()()
233  *mesh().V00()
234  )/mesh().V()
235  ),
236  rDeltaT.value()*rho.value()*
237  (
238  coefft*vf.boundaryField() -
239  (
240  coefft0*vf.oldTime().boundaryField()
241  - coefft00*vf.oldTime().oldTime().boundaryField()
242  )
243  )
244  )
245  );
246  }
247  else
248  {
249  return tmp<VolField<Type>>
250  (
251  VolField<Type>::New
252  (
253  ddtName,
254  rDeltaT*rho*
255  (
256  coefft*vf
257  - coefft0*vf.oldTime()
258  + coefft00*vf.oldTime().oldTime()
259  )
260  )
261  );
262  }
263 }
264 
265 
266 template<class Type>
267 tmp<VolField<Type>>
268 backwardDdtScheme<Type>::fvcDdt
269 (
270  const volScalarField& rho,
271  const VolField<Type>& vf
272 )
273 {
274  const dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
275 
276  const word ddtName("ddt("+rho.name()+','+vf.name()+')');
277 
278  const scalar deltaT = deltaT_();
279  const scalar deltaT0 = deltaT0_(vf);
280 
281  const scalar coefft = 1 + deltaT/(deltaT + deltaT0);
282  const scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
283  const scalar coefft0 = coefft + coefft00;
284 
285  if (mesh().moving())
286  {
287  return tmp<VolField<Type>>
288  (
289  VolField<Type>::New
290  (
291  ddtName,
292  rDeltaT*
293  (
294  coefft*rho()*vf() -
295  (
296  coefft0*rho.oldTime()()
297  *vf.oldTime()()*mesh().V0()
298  - coefft00*rho.oldTime().oldTime()()
299  *vf.oldTime().oldTime()()*mesh().V00()
300  )/mesh().V()
301  ),
302  rDeltaT.value()*
303  (
304  coefft*rho.boundaryField()*vf.boundaryField() -
305  (
306  coefft0*rho.oldTime().boundaryField()
307  *vf.oldTime().boundaryField()
308  - coefft00*rho.oldTime().oldTime().boundaryField()
309  *vf.oldTime().oldTime().boundaryField()
310  )
311  )
312  )
313  );
314  }
315  else
316  {
317  return tmp<VolField<Type>>
318  (
319  VolField<Type>::New
320  (
321  ddtName,
322  rDeltaT*
323  (
324  coefft*rho*vf
325  - coefft0*rho.oldTime()*vf.oldTime()
326  + coefft00*rho.oldTime().oldTime()*vf.oldTime().oldTime()
327  )
328  )
329  );
330  }
331 }
332 
333 
334 template<class Type>
335 tmp<VolField<Type>>
336 backwardDdtScheme<Type>::fvcDdt
337 (
338  const volScalarField& alpha,
339  const volScalarField& rho,
340  const VolField<Type>& vf
341 )
342 {
343  const dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
344 
345  const word ddtName("ddt("+alpha.name()+','+rho.name()+','+vf.name()+')');
346 
347  const scalar deltaT = deltaT_();
348  const scalar deltaT0 = deltaT0_(vf);
349 
350  const scalar coefft = 1 + deltaT/(deltaT + deltaT0);
351  const scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
352  const scalar coefft0 = coefft + coefft00;
353 
354  if (mesh().moving())
355  {
356  return tmp<VolField<Type>>
357  (
358  VolField<Type>::New
359  (
360  ddtName,
361  rDeltaT*
362  (
363  coefft*alpha()*rho()*vf()
364  - (
365  coefft0
366  *alpha.oldTime()()*rho.oldTime()()
367  *vf.oldTime()()*mesh().V0()
368  - coefft00
369  *alpha.oldTime().oldTime()()*rho.oldTime().oldTime()()
370  *vf.oldTime().oldTime()()*mesh().V00()
371  )/mesh().V()
372  ),
373  rDeltaT.value()*
374  (
375  coefft
376  *alpha.boundaryField()
377  *rho.boundaryField()
378  *vf.boundaryField() -
379  (
380  coefft0
381  *alpha.oldTime().boundaryField()
382  *rho.oldTime().boundaryField()
383  *vf.oldTime().boundaryField()
384 
385  - coefft00
386  *alpha.oldTime().oldTime().boundaryField()
387  *rho.oldTime().oldTime().boundaryField()
388  *vf.oldTime().oldTime().boundaryField()
389  )
390  )
391  )
392  );
393  }
394  else
395  {
396  return tmp<VolField<Type>>
397  (
398  VolField<Type>::New
399  (
400  ddtName,
401  rDeltaT*
402  (
403  coefft*alpha*rho*vf
404  - coefft0*alpha.oldTime()*rho.oldTime()*vf.oldTime()
405  + coefft00*alpha.oldTime().oldTime()
406  *rho.oldTime().oldTime()*vf.oldTime().oldTime()
407  )
408  )
409  );
410  }
411 }
412 
413 
414 template<class Type>
415 tmp<fvMatrix<Type>>
416 backwardDdtScheme<Type>::fvmDdt
417 (
418  const VolField<Type>& vf
419 )
420 {
421  tmp<fvMatrix<Type>> tfvm
422  (
423  new fvMatrix<Type>
424  (
425  vf,
426  vf.dimensions()*dimVolume/dimTime
427  )
428  );
429 
430  fvMatrix<Type>& fvm = tfvm.ref();
431 
432  const scalar rDeltaT = 1.0/deltaT_();
433 
434  const scalar deltaT = deltaT_();
435  const scalar deltaT0 = deltaT0_(vf);
436 
437  const scalar coefft = 1 + deltaT/(deltaT + deltaT0);
438  const scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
439  const scalar coefft0 = coefft + coefft00;
440 
441  fvm.diag() = (coefft*rDeltaT)*mesh().V();
442 
443  if (mesh().moving())
444  {
445  fvm.source() = rDeltaT*
446  (
447  coefft0*vf.oldTime().primitiveField()*mesh().V0()
448  - coefft00*vf.oldTime().oldTime().primitiveField()
449  *mesh().V00()
450  );
451  }
452  else
453  {
454  fvm.source() = rDeltaT*mesh().V()*
455  (
456  coefft0*vf.oldTime().primitiveField()
457  - coefft00*vf.oldTime().oldTime().primitiveField()
458  );
459  }
460 
461  return tfvm;
462 }
463 
464 
465 template<class Type>
466 tmp<fvMatrix<Type>>
467 backwardDdtScheme<Type>::fvmDdt
468 (
469  const dimensionedScalar& rho,
470  const VolField<Type>& vf
471 )
472 {
473  tmp<fvMatrix<Type>> tfvm
474  (
475  new fvMatrix<Type>
476  (
477  vf,
478  rho.dimensions()*vf.dimensions()*dimVolume/dimTime
479  )
480  );
481  fvMatrix<Type>& fvm = tfvm.ref();
482 
483  const scalar rDeltaT = 1.0/deltaT_();
484 
485  const scalar deltaT = deltaT_();
486  const scalar deltaT0 = deltaT0_(vf);
487 
488  const scalar coefft = 1 + deltaT/(deltaT + deltaT0);
489  const scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
490  const scalar coefft0 = coefft + coefft00;
491 
492  fvm.diag() = (coefft*rDeltaT*rho.value())*mesh().V();
493 
494  if (mesh().moving())
495  {
496  fvm.source() = rDeltaT*rho.value()*
497  (
498  coefft0*vf.oldTime().primitiveField()*mesh().V0()
499  - coefft00*vf.oldTime().oldTime().primitiveField()
500  *mesh().V00()
501  );
502  }
503  else
504  {
505  fvm.source() = rDeltaT*mesh().V()*rho.value()*
506  (
507  coefft0*vf.oldTime().primitiveField()
508  - coefft00*vf.oldTime().oldTime().primitiveField()
509  );
510  }
511 
512  return tfvm;
513 }
514 
515 
516 template<class Type>
517 tmp<fvMatrix<Type>>
518 backwardDdtScheme<Type>::fvmDdt
519 (
520  const volScalarField& rho,
521  const VolField<Type>& vf
522 )
523 {
524  tmp<fvMatrix<Type>> tfvm
525  (
526  new fvMatrix<Type>
527  (
528  vf,
529  rho.dimensions()*vf.dimensions()*dimVolume/dimTime
530  )
531  );
532  fvMatrix<Type>& fvm = tfvm.ref();
533 
534  const scalar rDeltaT = 1.0/deltaT_();
535 
536  const scalar deltaT = deltaT_();
537  const scalar deltaT0 = deltaT0_(vf);
538 
539  const scalar coefft = 1 + deltaT/(deltaT + deltaT0);
540  const scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
541  const scalar coefft0 = coefft + coefft00;
542 
543  fvm.diag() = (coefft*rDeltaT)*rho.primitiveField()*mesh().V();
544 
545  if (mesh().moving())
546  {
547  fvm.source() = rDeltaT*
548  (
549  coefft0*rho.oldTime().primitiveField()
550  *vf.oldTime().primitiveField()*mesh().V0()
551  - coefft00*rho.oldTime().oldTime().primitiveField()
552  *vf.oldTime().oldTime().primitiveField()*mesh().V00()
553  );
554  }
555  else
556  {
557  fvm.source() = rDeltaT*mesh().V()*
558  (
559  coefft0*rho.oldTime().primitiveField()
560  *vf.oldTime().primitiveField()
561  - coefft00*rho.oldTime().oldTime().primitiveField()
562  *vf.oldTime().oldTime().primitiveField()
563  );
564  }
565 
566  return tfvm;
567 }
568 
569 
570 template<class Type>
571 tmp<fvMatrix<Type>>
572 backwardDdtScheme<Type>::fvmDdt
573 (
574  const volScalarField& alpha,
575  const volScalarField& rho,
576  const VolField<Type>& vf
577 )
578 {
579  tmp<fvMatrix<Type>> tfvm
580  (
581  new fvMatrix<Type>
582  (
583  vf,
584  alpha.dimensions()
585  *rho.dimensions()
586  *vf.dimensions()
587  *dimVolume
588  /dimTime
589  )
590  );
591  fvMatrix<Type>& fvm = tfvm.ref();
592 
593  const scalar rDeltaT = 1.0/deltaT_();
594 
595  const scalar deltaT = deltaT_();
596  const scalar deltaT0 = deltaT0_(vf);
597 
598  const scalar coefft = 1 + deltaT/(deltaT + deltaT0);
599  const scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
600  const scalar coefft0 = coefft + coefft00;
601 
602  fvm.diag() =
603  (coefft*rDeltaT)*alpha.primitiveField()*rho.primitiveField()*mesh().V();
604 
605  if (mesh().moving())
606  {
607  fvm.source() = rDeltaT*
608  (
609  coefft0
610  *alpha.oldTime().primitiveField()
611  *rho.oldTime().primitiveField()
612  *vf.oldTime().primitiveField()*mesh().V0()
613 
614  - coefft00
615  *alpha.oldTime().oldTime().primitiveField()
616  *rho.oldTime().oldTime().primitiveField()
617  *vf.oldTime().oldTime().primitiveField()*mesh().V00()
618  );
619  }
620  else
621  {
622  fvm.source() = rDeltaT*mesh().V()*
623  (
624  coefft0
625  *alpha.oldTime().primitiveField()
626  *rho.oldTime().primitiveField()
627  *vf.oldTime().primitiveField()
628 
629  - coefft00
630  *alpha.oldTime().oldTime().primitiveField()
631  *rho.oldTime().oldTime().primitiveField()
632  *vf.oldTime().oldTime().primitiveField()
633  );
634  }
635 
636  return tfvm;
637 }
638 
639 
640 template<class Type>
641 tmp<typename backwardDdtScheme<Type>::fluxFieldType>
642 backwardDdtScheme<Type>::fvcDdtUfCorr
643 (
644  const VolField<Type>& U,
645  const SurfaceField<Type>& Uf
646 )
647 {
648  const dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
649 
650  const scalar deltaT = deltaT_();
651  const scalar deltaT0 = deltaT0_(U);
652 
653  const scalar coefft = 1 + deltaT/(deltaT + deltaT0);
654  const scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
655  const scalar coefft0 = coefft + coefft00;
656 
657  return fluxFieldType::New
658  (
659  "ddtCorr(" + U.name() + ',' + Uf.name() + ')',
660  this->fvcDdtPhiCoeff(U.oldTime(), (mesh().Sf() & Uf.oldTime()))
661  *rDeltaT
662  *(
663  mesh().Sf()
664  & (
665  (coefft0*Uf.oldTime() - coefft00*Uf.oldTime().oldTime())
666  - fvc::interpolate
667  (
668  coefft0*U.oldTime() - coefft00*U.oldTime().oldTime()
669  )
670  )
671  )
672  );
673 }
674 
675 
676 template<class Type>
677 tmp<typename backwardDdtScheme<Type>::fluxFieldType>
678 backwardDdtScheme<Type>::fvcDdtPhiCorr
679 (
680  const VolField<Type>& U,
681  const fluxFieldType& phi
682 )
683 {
684  const dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
685 
686  const scalar deltaT = deltaT_();
687  const scalar deltaT0 = deltaT0_(U);
688 
689  const scalar coefft = 1 + deltaT/(deltaT + deltaT0);
690  const scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
691  const scalar coefft0 = coefft + coefft00;
692 
693  return fluxFieldType::New
694  (
695  "ddtCorr(" + U.name() + ',' + phi.name() + ')',
696  this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())
697  *rDeltaT
698  *(
699  (coefft0*phi.oldTime() - coefft00*phi.oldTime().oldTime())
700  - fvc::dotInterpolate
701  (
702  mesh().Sf(),
703  coefft0*U.oldTime() - coefft00*U.oldTime().oldTime()
704  )
705  )
706  );
707 }
708 
709 
710 template<class Type>
711 tmp<typename backwardDdtScheme<Type>::fluxFieldType>
712 backwardDdtScheme<Type>::fvcDdtUfCorr
713 (
714  const volScalarField& rho,
715  const VolField<Type>& U,
716  const SurfaceField<Type>& rhoUf
717 )
718 {
719  const dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
720 
721  const scalar deltaT = deltaT_();
722  const scalar deltaT0 = deltaT0_(U);
723 
724  const scalar coefft = 1 + deltaT/(deltaT + deltaT0);
725  const scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
726  const scalar coefft0 = coefft + coefft00;
727 
728  if
729  (
730  U.dimensions() == dimVelocity
731  && rhoUf.dimensions() == rho.dimensions()*dimVelocity
732  )
733  {
734  const VolField<Type> rhoU0
735  (
736  rho.oldTime()*U.oldTime()
737  );
738 
739  const VolField<Type> rhoU00
740  (
741  rho.oldTime().oldTime()*U.oldTime().oldTime()
742  );
743 
744  return fluxFieldType::New
745  (
746  "ddtCorr(" + rho.name() + ',' + U.name() + ',' + rhoUf.name() + ')',
747  this->fvcDdtPhiCoeff
748  (
749  rhoU0,
750  mesh().Sf() & rhoUf.oldTime(),
751  rho.oldTime()
752  )
753  *rDeltaT
754  *(
755  mesh().Sf()
756  & (
757  (
758  coefft0*rhoUf.oldTime()
759  - coefft00*rhoUf.oldTime().oldTime()
760  )
761  - fvc::interpolate(coefft0*rhoU0 - coefft00*rhoU00)
762  )
763  )
764  );
765  }
766  else if
767  (
768  U.dimensions() == rho.dimensions()*dimVelocity
769  && rhoUf.dimensions() == rho.dimensions()*dimVelocity
770  )
771  {
772  return fluxFieldType::New
773  (
774  "ddtCorr(" + rho.name() + ',' + U.name() + ',' + rhoUf.name() + ')',
775  this->fvcDdtPhiCoeff
776  (
777  U.oldTime(),
778  mesh().Sf() & rhoUf.oldTime(),
779  rho.oldTime()
780  )
781  *rDeltaT
782  *(
783  mesh().Sf()
784  & (
785  (
786  coefft0*rhoUf.oldTime()
787  - coefft00*rhoUf.oldTime().oldTime()
788  )
789  - fvc::interpolate
790  (
791  coefft0*U.oldTime()
792  - coefft00*U.oldTime().oldTime()
793  )
794  )
795  )
796  );
797  }
798  else
799  {
800  FatalErrorInFunction
801  << "dimensions of phi are not correct"
802  << abort(FatalError);
803 
804  return fluxFieldType::null();
805  }
806 }
807 
808 
809 template<class Type>
810 tmp<typename backwardDdtScheme<Type>::fluxFieldType>
811 backwardDdtScheme<Type>::fvcDdtPhiCorr
812 (
813  const volScalarField& rho,
814  const VolField<Type>& U,
815  const fluxFieldType& phi
816 )
817 {
818  const dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
819 
820  const scalar deltaT = deltaT_();
821  const scalar deltaT0 = deltaT0_(U);
822 
823  const scalar coefft = 1 + deltaT/(deltaT + deltaT0);
824  const scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
825  const scalar coefft0 = coefft + coefft00;
826 
827  if
828  (
829  U.dimensions() == dimVelocity
830  && phi.dimensions() == rho.dimensions()*dimVolumetricFlux
831  )
832  {
833  const VolField<Type> rhoU0
834  (
835  rho.oldTime()*U.oldTime()
836  );
837 
838  const VolField<Type> rhoU00
839  (
840  rho.oldTime().oldTime()*U.oldTime().oldTime()
841  );
842 
843  return fluxFieldType::New
844  (
845  "ddtCorr(" + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
846  this->fvcDdtPhiCoeff(rhoU0, phi.oldTime(), rho.oldTime())
847  *rDeltaT
848  *(
849  (coefft0*phi.oldTime() - coefft00*phi.oldTime().oldTime())
850  - fvc::dotInterpolate
851  (
852  mesh().Sf(),
853  coefft0*rhoU0 - coefft00*rhoU00
854  )
855  )
856  );
857  }
858  else if
859  (
860  U.dimensions() == rho.dimensions()*dimVelocity
861  && phi.dimensions() == rho.dimensions()*dimVolumetricFlux
862  )
863  {
864  return fluxFieldType::New
865  (
866  "ddtCorr(" + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
867  this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime(), rho.oldTime())
868  *rDeltaT
869  *(
870  (coefft0*phi.oldTime() - coefft00*phi.oldTime().oldTime())
871  - fvc::dotInterpolate
872  (
873  mesh().Sf(),
874  coefft0*U.oldTime() - coefft00*U.oldTime().oldTime()
875  )
876  )
877  );
878  }
879  else
880  {
881  FatalErrorInFunction
882  << "dimensions of phi are not correct"
883  << abort(FatalError);
884 
885  return fluxFieldType::null();
886  }
887 }
888 
889 
890 template<class Type>
891 tmp<typename backwardDdtScheme<Type>::fluxFieldType>
892 backwardDdtScheme<Type>::fvcDdtUfCorr
893 (
894  const volScalarField& alpha,
895  const volScalarField& rho,
896  const VolField<Type>& U,
897  const SurfaceField<Type>& Uf
898 )
899 {
900  const dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
901 
902  const scalar deltaT = deltaT_();
903  const scalar deltaT0 = deltaT0_(U);
904 
905  const scalar coefft = 1 + deltaT/(deltaT + deltaT0);
906  const scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
907  const scalar coefft0 = coefft + coefft00;
908 
909  if (U.dimensions() == dimVelocity && Uf.dimensions() == dimVelocity)
910  {
911  const volScalarField alphaRho0(alpha.oldTime()*rho.oldTime());
912  const volScalarField alphaRho00
913  (
914  alpha.oldTime().oldTime()*rho.oldTime().oldTime()
915  );
916 
917  return fluxFieldType::New
918  (
919  "ddtCorr("
920  + alpha.name() + rho.name() + ',' + U.name() + ',' + Uf.name()
921  + ')',
922  this->fvcDdtPhiCoeff(U.oldTime(), mesh().Sf() & Uf.oldTime())
923  *rDeltaT
924  *(
925  mesh().Sf()
926  & (
927  (
928  coefft0*fvc::interpolate(alphaRho0)*Uf.oldTime()
929  - coefft00*fvc::interpolate(alphaRho00)
930  *Uf.oldTime().oldTime()
931  )
932  - fvc::interpolate
933  (
934  coefft0*alphaRho0*U.oldTime()
935  - coefft00*alphaRho00*U.oldTime().oldTime()
936  )
937  )
938  )
939  );
940  }
941  else
942  {
943  FatalErrorInFunction
944  << "dimensions of phi are not correct"
945  << abort(FatalError);
946 
947  return fluxFieldType::null();
948  }
949 }
950 
951 
952 template<class Type>
953 tmp<typename backwardDdtScheme<Type>::fluxFieldType>
954 backwardDdtScheme<Type>::fvcDdtPhiCorr
955 (
956  const volScalarField& alpha,
957  const volScalarField& rho,
958  const VolField<Type>& U,
959  const fluxFieldType& phi
960 )
961 {
962  const dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
963 
964  const scalar deltaT = deltaT_();
965  const scalar deltaT0 = deltaT0_(U);
966 
967  const scalar coefft = 1 + deltaT/(deltaT + deltaT0);
968  const scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
969  const scalar coefft0 = coefft + coefft00;
970 
971  if (U.dimensions() == dimVelocity && phi.dimensions() == dimVolumetricFlux)
972  {
973  const volScalarField alphaRho0(alpha.oldTime()*rho.oldTime());
974  const volScalarField alphaRho00
975  (
976  alpha.oldTime().oldTime()*rho.oldTime().oldTime()
977  );
978 
979  return fluxFieldType::New
980  (
981  "ddtCorr("
982  + alpha.name() + rho.name() + ',' + U.name() + ',' + phi.name()
983  + ')',
984  this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())
985  *rDeltaT
986  *(
987  (
988  coefft0*fvc::interpolate(alphaRho0)*phi.oldTime()
989  - coefft00*fvc::interpolate(alphaRho00)
990  *phi.oldTime().oldTime()
991  )
992  - fvc::dotInterpolate
993  (
994  mesh().Sf(),
995  coefft0*alphaRho0*U.oldTime()
996  - coefft00*alphaRho00*U.oldTime().oldTime()
997  )
998  )
999  );
1000  }
1001  else
1002  {
1003  FatalErrorInFunction
1004  << "dimensions of phi are not correct"
1005  << abort(FatalError);
1006 
1007  return fluxFieldType::null();
1008  }
1009 }
1010 
1011 
1012 template<class Type>
1013 tmp<surfaceScalarField> backwardDdtScheme<Type>::meshPhi
1014 (
1015  const VolField<Type>& vf
1016 )
1017 {
1018  const scalar deltaT = deltaT_();
1019  const scalar deltaT0 = deltaT0_(vf);
1020 
1021  // Coefficient for t-3/2 (between times 0 and 00)
1022  const scalar coefft0_00 = deltaT/(deltaT + deltaT0);
1023 
1024  // Coefficient for t-1/2 (between times n and 0)
1025  const scalar coefftn_0 = 1 + coefft0_00;
1026 
1027  return surfaceScalarField::New
1028  (
1029  mesh().phi().name(),
1030  coefftn_0*mesh().phi() - coefft0_00*mesh().phi().oldTime()
1031  );
1032 }
1033 
1034 
1035 template<class Type>
1036 tmp<scalarField> backwardDdtScheme<Type>::meshPhi
1037 (
1038  const VolField<Type>& vf,
1039  const label patchi
1040 )
1041 {
1042  const scalar deltaT = deltaT_();
1043  const scalar deltaT0 = deltaT0_(vf);
1044 
1045  // Coefficient for t-3/2 (between times 0 and 00)
1046  const scalar coefft0_00 = deltaT/(deltaT + deltaT0);
1047 
1048  // Coefficient for t-1/2 (between times n and 0)
1049  const scalar coefftn_0 = 1 + coefft0_00;
1050 
1051  return
1052  (
1053  coefftn_0*mesh().phi().boundaryField()[patchi]
1054  - coefft0_00*mesh().phi().oldTime().boundaryField()[patchi]
1055  );
1056 }
1057 
1058 
1059 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1060 
1061 } // End namespace fv
1062 
1063 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1064 
1065 } // End namespace Foam
1066 
1067 // ************************************************************************* //
Foam::DimensionedField::dimensions
const dimensionSet & dimensions() const
Return dimensions.
Definition: DimensionedFieldI.H:48
Foam::GeometricField
Generic GeometricField class.
Definition: GeometricField.H:81
Foam::GeometricField::boundaryField
const Boundary & boundaryField() const
Return const-reference to the boundary field.
Definition: GeometricFieldI.H:60
Foam::GeometricField::New
static tmp< GeometricField< Type, PatchField, GeoMesh > > New(const word &name, const Internal &, const PtrList< PatchField< Type >> &, const HashPtrTable< Source > &=HashPtrTable< Source >())
Return a temporary field constructed from name,.
Definition: GeometricField.C:709
Foam::IOobject::name
const word & name() const
Return name.
Definition: IOobject.H:310
Foam::OldTimeField::oldTime
const FieldType & oldTime() const
Return the old time field.
Definition: OldTimeField.C:315
Foam::calculatedFvPatchField
This boundary condition is not designed to be evaluated; it is assumed that the value is assigned via...
Definition: calculatedFvPatchField.H:64
Foam::dimensioned
Generic dimensioned Type class.
Definition: dimensionedType.H:62
Foam::dimensioned::dimensions
const dimensionSet & dimensions() const
Return const reference to dimensions.
Definition: dimensionedType.C:274
Foam::dimensioned::value
const Type & value() const
Return const reference to value.
Definition: dimensionedType.C:287
Foam::dimensioned::name
const word & name() const
Return const reference to name.
Definition: dimensionedType.C:261
Foam::fvMatrix
A special matrix type and solver, designed for finite volume solutions of scalar equations....
Definition: fvMatrix.H:118
Foam::fvMatrix::source
Field< Type > & source()
Definition: fvMatrix.H:307
Foam::fv::backwardDdtScheme::fvcDdtPhiCorr
virtual tmp< fluxFieldType > fvcDdtPhiCorr(const VolField< Type > &U, const fluxFieldType &phi)
Definition: backwardDdtScheme.C:679
Foam::fv::backwardDdtScheme::fvmDdt
virtual tmp< fvMatrix< Type > > fvmDdt(const VolField< Type > &)
Definition: backwardDdtScheme.C:417
Foam::fv::backwardDdtScheme::meshPhi
virtual tmp< surfaceScalarField > meshPhi(const VolField< Type > &)
Definition: backwardDdtScheme.C:1014
Foam::fv::backwardDdtScheme::fvcDdt
virtual tmp< VolField< Type > > fvcDdt(const dimensioned< Type > &)
Definition: backwardDdtScheme.C:77
Foam::fv::backwardDdtScheme::fvcDdtUfCorr
virtual tmp< fluxFieldType > fvcDdtUfCorr(const VolField< Type > &U, const SurfaceField< Type > &Uf)
Definition: backwardDdtScheme.C:643
Foam::lduMatrix::diag
scalarField & diag()
Definition: lduMatrix.C:186
Foam::tmp
A class for managing temporary objects.
Definition: tmp.H:55
Foam::tmp::ref
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:181
Foam::word
A class for handling words, derived from string.
Definition: word.H:62
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:334
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.
patchi
label patchi
Definition: getPatchFieldScalar.H:1
U
U
Definition: pEqn.H:72
alpha
volScalarField alpha(IOobject("alpha", runTime.name(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), lambda *max(Ua &U, zeroSensitivity))
Foam::compressible::New
autoPtr< CompressibleMomentumTransportModel > New(const volScalarField &rho, const volVectorField &U, const surfaceScalarField &phi, const viscosity &viscosity)
Definition: compressibleMomentumTransportModelTemplates.C:34
Foam::fvc::interpolate
static tmp< SurfaceField< Type > > interpolate(const VolField< Type > &tvf, const surfaceScalarField &faceFlux, Istream &schemeData)
Interpolate field onto faces using scheme given by Istream.
Foam::fvc::dotInterpolate
static tmp< SurfaceField< typename innerProduct< vector, Type >::type > > dotInterpolate(const surfaceVectorField &Sf, const VolField< Type > &tvf)
Interpolate field onto faces.
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.H:214
Foam::Zero
static const zero Zero
Definition: zero.H:97
Foam::label
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
Foam::name
word name(const bool)
Return a word representation of a bool.
Definition: boolIO.C:39
Foam::dimVolumetricFlux
const dimensionSet dimVolumetricFlux
Foam::abort
errorManip< error > abort(error &err)
Definition: errorManip.H:131
Foam::dimTime
const dimensionSet dimTime
Foam::dimVolume
const dimensionSet dimVolume
Foam::dimVelocity
const dimensionSet dimVelocity
Foam::FatalError
error FatalError
fv
labelList fv(nPoints)