FreeFOAM The Cross-Platform CFD Toolkit
backwardDdtScheme.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) 1991-2010 OpenCFD Ltd.
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"
28 #include <finiteVolume/fvcDiv.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().deltaT().value();
47 }
48 
49 
50 template<class Type>
51 scalar backwardDdtScheme<Type>::deltaT0_() const
52 {
53  return mesh().time().deltaT0().value();
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<GeometricField<Type, fvPatchField, volMesh> >
77 (
78  const dimensioned<Type>& dt
79 )
80 {
81  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
82 
83  IOobject ddtIOobject
84  (
85  "ddt("+dt.name()+')',
86  mesh().time().timeName(),
87  mesh()
88  );
89 
90  scalar deltaT = deltaT_();
91  scalar deltaT0 = deltaT0_();
92 
93  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
94  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
95  scalar coefft0 = coefft + coefft00;
96 
97  if (mesh().moving())
98  {
100  (
102  (
103  ddtIOobject,
104  mesh(),
106  (
107  "0",
108  dt.dimensions()/dimTime,
110  )
111  )
112  );
113 
114  tdtdt().internalField() = rDeltaT.value()*dt.value()*
115  (
116  coefft - (coefft0*mesh().V0() - coefft00*mesh().V00())/mesh().V()
117  );
118 
119  return tdtdt;
120  }
121  else
122  {
124  (
126  (
127  ddtIOobject,
128  mesh(),
130  (
131  "0",
132  dt.dimensions()/dimTime,
134  ),
136  )
137  );
138  }
139 }
140 
141 
142 template<class Type>
145 (
147 )
148 {
149  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
150 
151  IOobject ddtIOobject
152  (
153  "ddt("+vf.name()+')',
154  mesh().time().timeName(),
155  mesh()
156  );
157 
158  scalar deltaT = deltaT_();
159  scalar deltaT0 = deltaT0_(vf);
160 
161  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
162  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
163  scalar coefft0 = coefft + coefft00;
164 
165  if (mesh().moving())
166  {
168  (
170  (
171  ddtIOobject,
172  mesh(),
173  rDeltaT.dimensions()*vf.dimensions(),
174  rDeltaT.value()*
175  (
176  coefft*vf.internalField() -
177  (
178  coefft0*vf.oldTime().internalField()*mesh().V0()
179  - coefft00*vf.oldTime().oldTime().internalField()
180  *mesh().V00()
181  )/mesh().V()
182  ),
183  rDeltaT.value()*
184  (
185  coefft*vf.boundaryField() -
186  (
187  coefft0*vf.oldTime().boundaryField()
188  - coefft00*vf.oldTime().oldTime().boundaryField()
189  )
190  )
191  )
192  );
193  }
194  else
195  {
197  (
199  (
200  ddtIOobject,
201  rDeltaT*
202  (
203  coefft*vf
204  - coefft0*vf.oldTime()
205  + coefft00*vf.oldTime().oldTime()
206  )
207  )
208  );
209  }
210 }
211 
212 
213 template<class Type>
216 (
217  const dimensionedScalar& rho,
219 )
220 {
221  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
222 
223  IOobject ddtIOobject
224  (
225  "ddt("+rho.name()+','+vf.name()+')',
226  mesh().time().timeName(),
227  mesh()
228  );
229 
230  scalar deltaT = deltaT_();
231  scalar deltaT0 = deltaT0_(vf);
232 
233  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
234  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
235  scalar coefft0 = coefft + coefft00;
236 
237  if (mesh().moving())
238  {
240  (
242  (
243  ddtIOobject,
244  mesh(),
245  rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
246  rDeltaT.value()*rho.value()*
247  (
248  coefft*vf.internalField() -
249  (
250  coefft0*vf.oldTime().internalField()*mesh().V0()
251  - coefft00*vf.oldTime().oldTime().internalField()
252  *mesh().V00()
253  )/mesh().V()
254  ),
255  rDeltaT.value()*rho.value()*
256  (
257  coefft*vf.boundaryField() -
258  (
259  coefft0*vf.oldTime().boundaryField()
260  - coefft00*vf.oldTime().oldTime().boundaryField()
261  )
262  )
263  )
264  );
265  }
266  else
267  {
269  (
271  (
272  ddtIOobject,
273  rDeltaT*rho*
274  (
275  coefft*vf
276  - coefft0*vf.oldTime()
277  + coefft00*vf.oldTime().oldTime()
278  )
279  )
280  );
281  }
282 }
283 
284 template<class Type>
287 (
288  const volScalarField& rho,
290 )
291 {
292  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
293 
294  IOobject ddtIOobject
295  (
296  "ddt("+rho.name()+','+vf.name()+')',
297  mesh().time().timeName(),
298  mesh()
299  );
300 
301  scalar deltaT = deltaT_();
302  scalar deltaT0 = deltaT0_(vf);
303 
304  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
305  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
306  scalar coefft0 = coefft + coefft00;
307 
308  if (mesh().moving())
309  {
311  (
313  (
314  ddtIOobject,
315  mesh(),
316  rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
317  rDeltaT.value()*
318  (
319  coefft*rho.internalField()*vf.internalField() -
320  (
321  coefft0*rho.oldTime().internalField()
322  *vf.oldTime().internalField()*mesh().V0()
323  - coefft00*rho.oldTime().oldTime().internalField()
324  *vf.oldTime().oldTime().internalField()*mesh().V00()
325  )/mesh().V()
326  ),
327  rDeltaT.value()*
328  (
329  coefft*rho.boundaryField()*vf.boundaryField() -
330  (
331  coefft0*rho.oldTime().boundaryField()
332  *vf.oldTime().boundaryField()
333  - coefft00*rho.oldTime().oldTime().boundaryField()
334  *vf.oldTime().oldTime().boundaryField()
335  )
336  )
337  )
338  );
339  }
340  else
341  {
343  (
345  (
346  ddtIOobject,
347  rDeltaT*
348  (
349  coefft*rho*vf
350  - coefft0*rho.oldTime()*vf.oldTime()
351  + coefft00*rho.oldTime().oldTime()*vf.oldTime().oldTime()
352  )
353  )
354  );
355  }
356 }
357 
358 
359 template<class Type>
362 (
364 )
365 {
366  tmp<fvMatrix<Type> > tfvm
367  (
368  new fvMatrix<Type>
369  (
370  vf,
372  )
373  );
374 
375  fvMatrix<Type>& fvm = tfvm();
376 
377  scalar rDeltaT = 1.0/deltaT_();
378 
379  scalar deltaT = deltaT_();
380  scalar deltaT0 = deltaT0_(vf);
381 
382  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
383  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
384  scalar coefft0 = coefft + coefft00;
385 
386  fvm.diag() = (coefft*rDeltaT)*mesh().V();
387 
388  if (mesh().moving())
389  {
390  fvm.source() = rDeltaT*
391  (
392  coefft0*vf.oldTime().internalField()*mesh().V0()
393  - coefft00*vf.oldTime().oldTime().internalField()
394  *mesh().V00()
395  );
396  }
397  else
398  {
399  fvm.source() = rDeltaT*mesh().V()*
400  (
401  coefft0*vf.oldTime().internalField()
402  - coefft00*vf.oldTime().oldTime().internalField()
403  );
404  }
405 
406  return tfvm;
407 }
408 
409 
410 template<class Type>
413 (
414  const dimensionedScalar& rho,
416 )
417 {
418  tmp<fvMatrix<Type> > tfvm
419  (
420  new fvMatrix<Type>
421  (
422  vf,
424  )
425  );
426  fvMatrix<Type>& fvm = tfvm();
427 
428  scalar rDeltaT = 1.0/deltaT_();
429 
430  scalar deltaT = deltaT_();
431  scalar deltaT0 = deltaT0_(vf);
432 
433  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
434  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
435  scalar coefft0 = coefft + coefft00;
436 
437  fvm.diag() = (coefft*rDeltaT*rho.value())*mesh().V();
438 
439  if (mesh().moving())
440  {
441  fvm.source() = rDeltaT*rho.value()*
442  (
443  coefft0*vf.oldTime().internalField()*mesh().V0()
444  - coefft00*vf.oldTime().oldTime().internalField()
445  *mesh().V00()
446  );
447  }
448  else
449  {
450  fvm.source() = rDeltaT*mesh().V()*rho.value()*
451  (
452  coefft0*vf.oldTime().internalField()
453  - coefft00*vf.oldTime().oldTime().internalField()
454  );
455  }
456 
457  return tfvm;
458 }
459 
460 
461 template<class Type>
464 (
465  const volScalarField& rho,
467 )
468 {
469  tmp<fvMatrix<Type> > tfvm
470  (
471  new fvMatrix<Type>
472  (
473  vf,
475  )
476  );
477  fvMatrix<Type>& fvm = tfvm();
478 
479  scalar rDeltaT = 1.0/deltaT_();
480 
481  scalar deltaT = deltaT_();
482  scalar deltaT0 = deltaT0_(vf);
483 
484  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
485  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
486  scalar coefft0 = coefft + coefft00;
487 
488  fvm.diag() = (coefft*rDeltaT)*rho.internalField()*mesh().V();
489 
490  if (mesh().moving())
491  {
492  fvm.source() = rDeltaT*
493  (
494  coefft0*rho.oldTime().internalField()
495  *vf.oldTime().internalField()*mesh().V0()
496  - coefft00*rho.oldTime().oldTime().internalField()
497  *vf.oldTime().oldTime().internalField()*mesh().V00()
498  );
499  }
500  else
501  {
502  fvm.source() = rDeltaT*mesh().V()*
503  (
504  coefft0*rho.oldTime().internalField()
505  *vf.oldTime().internalField()
506  - coefft00*rho.oldTime().oldTime().internalField()
507  *vf.oldTime().oldTime().internalField()
508  );
509  }
510 
511  return tfvm;
512 }
513 
514 
515 template<class Type>
518 (
519  const volScalarField& rA,
521  const fluxFieldType& phi
522 )
523 {
524  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
525 
526  IOobject ddtIOobject
527  (
528  "ddtPhiCorr(" + rA.name() + ',' + U.name() + ',' + phi.name() + ')',
529  mesh().time().timeName(),
530  mesh()
531  );
532 
533  scalar deltaT = deltaT_();
534  scalar deltaT0 = deltaT0_(U);
535 
536  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
537  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
538  scalar coefft0 = coefft + coefft00;
539 
540  return tmp<fluxFieldType>
541  (
542  new fluxFieldType
543  (
544  ddtIOobject,
545  rDeltaT*fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())
546  *(
547  fvc::interpolate(rA)
548  *(
549  coefft0*phi.oldTime()
550  - coefft00*phi.oldTime().oldTime()
551  )
552  - (
554  (
555  rA*
556  (
557  coefft0*U.oldTime()
558  - coefft00*U.oldTime().oldTime()
559  )
560  ) & mesh().Sf()
561  )
562  )
563  )
564  );
565 }
566 
567 
568 template<class Type>
571 (
572  const volScalarField& rA,
573  const volScalarField& rho,
575  const fluxFieldType& phiAbs
576 )
577 {
578  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
579 
580  IOobject ddtIOobject
581  (
582  "ddtPhiCorr("
583  + rA.name() + ','
584  + rho.name() + ','
585  + U.name() + ','
586  + phiAbs.name() + ')',
587  mesh().time().timeName(),
588  mesh()
589  );
590 
591  scalar deltaT = deltaT_();
592  scalar deltaT0 = deltaT0_(U);
593 
594  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
595  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
596  scalar coefft0 = coefft + coefft00;
597 
598  if
599  (
600  U.dimensions() == dimVelocity
601  && phiAbs.dimensions() == dimVelocity*dimArea
602  )
603  {
604  return tmp<fluxFieldType>
605  (
606  new fluxFieldType
607  (
608  ddtIOobject,
609  rDeltaT*fvcDdtPhiCoeff(U.oldTime(), phiAbs.oldTime())
610  *(
611  coefft0*fvc::interpolate(rA*rho.oldTime())
612  *phiAbs.oldTime()
613  - coefft00*fvc::interpolate(rA*rho.oldTime().oldTime())
614  *phiAbs.oldTime().oldTime()
615  - (
617  (
618  rA*
619  (
620  coefft0*rho.oldTime()*U.oldTime()
621  - coefft00*rho.oldTime().oldTime()
622  *U.oldTime().oldTime()
623  )
624  ) & mesh().Sf()
625  )
626  )
627  )
628  );
629  }
630  else if
631  (
632  U.dimensions() == dimVelocity
633  && phiAbs.dimensions() == rho.dimensions()*dimVelocity*dimArea
634  )
635  {
636  return tmp<fluxFieldType>
637  (
638  new fluxFieldType
639  (
640  ddtIOobject,
641  rDeltaT
642  *fvcDdtPhiCoeff
643  (
644  U.oldTime(),
645  phiAbs.oldTime()/fvc::interpolate(rho.oldTime())
646  )
647  *(
648  fvc::interpolate(rA*rho.oldTime())
649  *(
650  coefft0*phiAbs.oldTime()
651  /fvc::interpolate(rho.oldTime())
652  - coefft00*phiAbs.oldTime().oldTime()
654  )
655  - (
657  (
658  rA*rho.oldTime()*
659  (
660  coefft0*U.oldTime()
661  - coefft00*U.oldTime().oldTime()
662  )
663  ) & mesh().Sf()
664  )
665  )
666  )
667  );
668  }
669  else if
670  (
671  U.dimensions() == rho.dimensions()*dimVelocity
672  && phiAbs.dimensions() == rho.dimensions()*dimVelocity*dimArea
673  )
674  {
675  return tmp<fluxFieldType>
676  (
677  new fluxFieldType
678  (
679  ddtIOobject,
680  rDeltaT
681  *fvcDdtPhiCoeff(rho.oldTime(), U.oldTime(), phiAbs.oldTime())
682  *(
683  fvc::interpolate(rA)
684  *(
685  coefft0*phiAbs.oldTime()
686  - coefft00*phiAbs.oldTime().oldTime()
687  )
688  - (
690  (
691  rA*
692  (
693  coefft0*U.oldTime()
694  - coefft00*U.oldTime().oldTime()
695  )
696  ) & mesh().Sf()
697  )
698  )
699  )
700  );
701  }
702  else
703  {
705  (
706  "backwardDdtScheme<Type>::fvcDdtPhiCorr"
707  ) << "dimensions of phiAbs are not correct"
708  << abort(FatalError);
709 
710  return fluxFieldType::null();
711  }
712 }
713 
714 
715 template<class Type>
717 (
719 )
720 {
721  scalar deltaT = deltaT_();
722  scalar deltaT0 = deltaT0_(vf);
723 
724  // Coefficient for t-3/2 (between times 0 and 00)
725  scalar coefft0_00 = deltaT/(deltaT + deltaT0);
726 
727  // Coefficient for t-1/2 (between times n and 0)
728  scalar coefftn_0 = 1 + coefft0_00;
729 
730  return coefftn_0*mesh().phi() - coefft0_00*mesh().phi().oldTime();
731 }
732 
733 
734 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
735 
736 } // End namespace fv
737 
738 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
739 
740 } // End namespace Foam
741 
742 // ************************ vim: set sw=4 sts=4 et: ************************ //