FreeFOAM The Cross-Platform CFD Toolkit
boundedBackwardDdtScheme.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 
28 #include <finiteVolume/fvcDiv.H>
29 
30 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
31 
32 namespace Foam
33 {
34 
35 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
36 
37 namespace fv
38 {
39 
40 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
41 
42 scalar boundedBackwardDdtScheme::deltaT_() const
43 {
44  return mesh().time().deltaT().value();
45 }
46 
47 
48 scalar boundedBackwardDdtScheme::deltaT0_() const
49 {
50  return mesh().time().deltaT0().value();
51 }
52 
53 
54 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
55 
56 tmp<volScalarField>
58 (
59  const dimensionedScalar& dt
60 )
61 {
62  // No change compared to backward
63 
64  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
65 
66  IOobject ddtIOobject
67  (
68  "ddt("+dt.name()+')',
69  mesh().time().timeName(),
70  mesh()
71  );
72 
73  scalar deltaT = deltaT_();
74  scalar deltaT0 = deltaT0_();
75 
76  scalar coefft = 1 + deltaT/(deltaT + deltaT0);
77  scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
78  scalar coefft0 = coefft + coefft00;
79 
80  if (mesh().moving())
81  {
83  (
84  new volScalarField
85  (
86  ddtIOobject,
87  mesh(),
89  (
90  "0",
91  dt.dimensions()/dimTime,
92  0.0
93  )
94  )
95  );
96 
97  tdtdt().internalField() = rDeltaT.value()*dt.value()*
98  (
99  coefft - (coefft0*mesh().V0() - coefft00*mesh().V00())/mesh().V()
100  );
101 
102  return tdtdt;
103  }
104  else
105  {
106  return tmp<volScalarField>
107  (
108  new volScalarField
109  (
110  ddtIOobject,
111  mesh(),
113  (
114  "0",
115  dt.dimensions()/dimTime,
116  0.0
117  ),
118  calculatedFvPatchScalarField::typeName
119  )
120  );
121  }
122 }
123 
124 
127 (
128  const volScalarField& vf
129 )
130 {
131  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
132 
133  IOobject ddtIOobject
134  (
135  "ddt("+vf.name()+')',
136  mesh().time().timeName(),
137  mesh()
138  );
139 
140  scalar deltaT = deltaT_();
141  scalar deltaT0 = deltaT0_(vf);
142 
143  // Calculate unboundedness indicator
144  // Note: all times moved by one because access to internal field
145  // copies current field into the old-time level.
146  volScalarField phict =
147  mag
148  (
149  vf.oldTime().oldTime()
150  - vf.oldTime().oldTime().oldTime()
151  )/
152  (
153  mag
154  (
155  vf.oldTime()
156  - vf.oldTime().oldTime()
157  )
158  + dimensionedScalar("small", vf.dimensions(), VSMALL)
159  );
160 
161  volScalarField limiter(pos(phict) - pos(phict - scalar(1)));
162 
163  volScalarField coefft = scalar(1) + limiter*deltaT/(deltaT + deltaT0);
164  volScalarField coefft00 = limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0));
165  volScalarField coefft0 = coefft + coefft00;
166 
167  if (mesh().moving())
168  {
169  return tmp<volScalarField>
170  (
171  new volScalarField
172  (
173  ddtIOobject,
174  mesh(),
175  rDeltaT.dimensions()*vf.dimensions(),
176  rDeltaT.value()*
177  (
178  coefft*vf.internalField() -
179  (
180  coefft0.internalField()
181  *vf.oldTime().internalField()*mesh().V0()
182  - coefft00.internalField()
183  *vf.oldTime().oldTime().internalField()
184  *mesh().V00()
185  )/mesh().V()
186  ),
187  rDeltaT.value()*
188  (
189  coefft.boundaryField()*vf.boundaryField() -
190  (
191  coefft0.boundaryField()*
192  vf.oldTime().boundaryField()
193  - coefft00.boundaryField()*
194  vf.oldTime().oldTime().boundaryField()
195  )
196  )
197  )
198  );
199  }
200  else
201  {
202  return tmp<volScalarField>
203  (
204  new volScalarField
205  (
206  ddtIOobject,
207  rDeltaT*
208  (
209  coefft*vf
210  - coefft0*vf.oldTime()
211  + coefft00*vf.oldTime().oldTime()
212  )
213  )
214  );
215  }
216 }
217 
218 
221 (
222  const dimensionedScalar& rho,
223  const volScalarField& vf
224 )
225 {
226  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
227 
228  IOobject ddtIOobject
229  (
230  "ddt("+rho.name()+','+vf.name()+')',
231  mesh().time().timeName(),
232  mesh()
233  );
234 
235  scalar deltaT = deltaT_();
236  scalar deltaT0 = deltaT0_(vf);
237 
238  // Calculate unboundedness indicator
239  // Note: all times moved by one because access to internal field
240  // copies current field into the old-time level.
241  volScalarField phict =
242  mag
243  (
244  vf.oldTime().oldTime()
245  - vf.oldTime().oldTime().oldTime()
246  )/
247  (
248  mag
249  (
250  vf.oldTime()
251  - vf.oldTime().oldTime()
252  )
253  + dimensionedScalar("small", vf.dimensions(), VSMALL)
254  );
255 
256  volScalarField limiter(pos(phict) - pos(phict - scalar(1)));
257 
258  volScalarField coefft = scalar(1) + limiter*deltaT/(deltaT + deltaT0);
259  volScalarField coefft00 = limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0));
260  volScalarField coefft0 = coefft + coefft00;
261 
262  if (mesh().moving())
263  {
264  return tmp<volScalarField>
265  (
266  new volScalarField
267  (
268  ddtIOobject,
269  mesh(),
270  rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
271  rDeltaT.value()*rho.value()*
272  (
273  coefft*vf.internalField() -
274  (
275  coefft0.internalField()*
276  vf.oldTime().internalField()*mesh().V0()
277  - coefft00.internalField()*
278  vf.oldTime().oldTime().internalField()
279  *mesh().V00()
280  )/mesh().V()
281  ),
282  rDeltaT.value()*rho.value()*
283  (
284  coefft.boundaryField()*vf.boundaryField() -
285  (
286  coefft0.boundaryField()*
287  vf.oldTime().boundaryField()
288  - coefft00.boundaryField()*
289  vf.oldTime().oldTime().boundaryField()
290  )
291  )
292  )
293  );
294  }
295  else
296  {
297  return tmp<volScalarField>
298  (
299  new volScalarField
300  (
301  ddtIOobject,
302  rDeltaT*rho*
303  (
304  coefft*vf
305  - coefft0*vf.oldTime()
306  + coefft00*vf.oldTime().oldTime()
307  )
308  )
309  );
310  }
311 }
312 
313 
316 (
317  const volScalarField& rho,
318  const volScalarField& vf
319 )
320 {
321  dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
322 
323  IOobject ddtIOobject
324  (
325  "ddt("+rho.name()+','+vf.name()+')',
326  mesh().time().timeName(),
327  mesh()
328  );
329 
330  scalar deltaT = deltaT_();
331  scalar deltaT0 = deltaT0_(vf);
332 
333  // Calculate unboundedness indicator
334  // Note: all times moved by one because access to internal field
335  // copies current field into the old-time level.
336  volScalarField phict =
337  mag
338  (
339  rho.oldTime().oldTime()*vf.oldTime().oldTime()
340  - rho.oldTime().oldTime().oldTime()*vf.oldTime().oldTime().oldTime()
341  )/
342  (
343  mag
344  (
345  rho.oldTime()*vf.oldTime()
346  - rho.oldTime().oldTime()*vf.oldTime().oldTime()
347  )
348  + dimensionedScalar("small", rho.dimensions()*vf.dimensions(), VSMALL)
349  );
350 
351  volScalarField limiter(pos(phict) - pos(phict - scalar(1)));
352 
353  volScalarField coefft = scalar(1) + limiter*deltaT/(deltaT + deltaT0);
354  volScalarField coefft00 = limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0));
355  volScalarField coefft0 = coefft + coefft00;
356 
357  if (mesh().moving())
358  {
359  return tmp<volScalarField>
360  (
361  new volScalarField
362  (
363  ddtIOobject,
364  mesh(),
365  rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
366  rDeltaT.value()*
367  (
368  coefft*rho.internalField()*vf.internalField() -
369  (
370  coefft0.internalField()*
371  rho.oldTime().internalField()*
372  vf.oldTime().internalField()*mesh().V0()
373  - coefft00.internalField()*
374  rho.oldTime().oldTime().internalField()
375  *vf.oldTime().oldTime().internalField()*mesh().V00()
376  )/mesh().V()
377  ),
378  rDeltaT.value()*
379  (
380  coefft.boundaryField()*vf.boundaryField() -
381  (
382  coefft0.boundaryField()*
383  rho.oldTime().boundaryField()*
384  vf.oldTime().boundaryField()
385  - coefft00.boundaryField()*
386  rho.oldTime().oldTime().boundaryField()*
387  vf.oldTime().oldTime().boundaryField()
388  )
389  )
390  )
391  );
392  }
393  else
394  {
395  return tmp<volScalarField>
396  (
397  new volScalarField
398  (
399  ddtIOobject,
400  rDeltaT*
401  (
402  coefft*rho*vf
403  - coefft0*rho.oldTime()*vf.oldTime()
404  + coefft00*rho.oldTime().oldTime()*vf.oldTime().oldTime()
405  )
406  )
407  );
408  }
409 }
410 
411 
414 (
415  volScalarField& vf
416 )
417 {
419  (
420  new fvScalarMatrix
421  (
422  vf,
424  )
425  );
426 
427  fvScalarMatrix& fvm = tfvm();
428 
429  scalar rDeltaT = 1.0/deltaT_();
430 
431  scalar deltaT = deltaT_();
432  scalar deltaT0 = deltaT0_(vf);
433 
434  // Calculate unboundedness indicator
435  // Note: all times moved by one because access to internal field
436  // copies current field into the old-time level.
437  scalarField phict =
438  mag
439  (
440  vf.oldTime().oldTime().internalField()
441  - vf.oldTime().oldTime().oldTime().internalField()
442  )/
443  (
444  mag
445  (
446  vf.oldTime().internalField()
447  - vf.oldTime().oldTime().internalField()
448  )
449  + VSMALL
450  );
451 
452  scalarField limiter(pos(phict) - pos(phict - 1.0));
453 
454  scalarField coefft = 1.0 + limiter*deltaT/(deltaT + deltaT0);
455  scalarField coefft00 = limiter*deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
456  scalarField coefft0 = coefft + coefft00;
457 
458  fvm.diag() = (coefft*rDeltaT)*mesh().V();
459 
460  if (mesh().moving())
461  {
462  fvm.source() = rDeltaT*
463  (
464  coefft0*vf.oldTime().internalField()*mesh().V0()
465  - coefft00*vf.oldTime().oldTime().internalField()
466  *mesh().V00()
467  );
468  }
469  else
470  {
471  fvm.source() = rDeltaT*mesh().V()*
472  (
473  coefft0*vf.oldTime().internalField()
474  - coefft00*vf.oldTime().oldTime().internalField()
475  );
476  }
477 
478  return tfvm;
479 }
480 
481 
484 (
485  const dimensionedScalar& rho,
486  volScalarField& vf
487 )
488 {
490  (
491  new fvScalarMatrix
492  (
493  vf,
495  )
496  );
497  fvScalarMatrix& fvm = tfvm();
498 
499  scalar rDeltaT = 1.0/deltaT_();
500 
501  scalar deltaT = deltaT_();
502  scalar deltaT0 = deltaT0_(vf);
503 
504  // Calculate unboundedness indicator
505  // Note: all times moved by one because access to internal field
506  // copies current field into the old-time level.
507  scalarField phict =
508  mag
509  (
510  vf.oldTime().oldTime().internalField()
511  - vf.oldTime().oldTime().oldTime().internalField()
512  )/
513  (
514  mag
515  (
516  vf.oldTime().internalField()
517  - vf.oldTime().oldTime().internalField()
518  )
519  + VSMALL
520  );
521 
522  scalarField limiter(pos(phict) - pos(phict - 1.0));
523 
524  scalarField coefft = 1.0 + limiter*deltaT/(deltaT + deltaT0);
525  scalarField coefft00 = limiter*deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
526  scalarField coefft0 = coefft + coefft00;
527 
528  fvm.diag() = (coefft*rDeltaT*rho.value())*mesh().V();
529 
530  if (mesh().moving())
531  {
532  fvm.source() = rDeltaT*rho.value()*
533  (
534  coefft0*vf.oldTime().internalField()*mesh().V0()
535  - coefft00*vf.oldTime().oldTime().internalField()
536  *mesh().V00()
537  );
538  }
539  else
540  {
541  fvm.source() = rDeltaT*mesh().V()*rho.value()*
542  (
543  coefft0*vf.oldTime().internalField()
544  - coefft00*vf.oldTime().oldTime().internalField()
545  );
546  }
547 
548  return tfvm;
549 }
550 
551 
554 (
555  const volScalarField& rho,
556  volScalarField& vf
557 )
558 {
560  (
561  new fvScalarMatrix
562  (
563  vf,
565  )
566  );
567  fvScalarMatrix& fvm = tfvm();
568 
569  scalar rDeltaT = 1.0/deltaT_();
570 
571  scalar deltaT = deltaT_();
572  scalar deltaT0 = deltaT0_(vf);
573 
574  // Calculate unboundedness indicator
575  // Note: all times moved by one because access to internal field
576  // copies current field into the old-time level.
577  scalarField phict =
578  mag
579  (
580  rho.oldTime().oldTime().internalField()*
581  vf.oldTime().oldTime().internalField()
582  - rho.oldTime().oldTime().oldTime().internalField()*
584  )/
585  (
586  mag
587  (
588  rho.oldTime().internalField()*
589  vf.oldTime().internalField()
590  - rho.oldTime().oldTime().internalField()*
591  vf.oldTime().oldTime().internalField()
592  )
593  + VSMALL
594  );
595 
596  scalarField limiter(pos(phict) - pos(phict - 1.0));
597 
598  scalarField coefft = 1.0 + limiter*deltaT/(deltaT + deltaT0);
599  scalarField coefft00 = limiter*deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
600  scalarField coefft0 = coefft + coefft00;
601 
602  fvm.diag() = (coefft*rDeltaT)*rho.internalField()*mesh().V();
603 
604  if (mesh().moving())
605  {
606  fvm.source() = rDeltaT*
607  (
608  coefft0*rho.oldTime().internalField()
609  *vf.oldTime().internalField()*mesh().V0()
610  - coefft00*rho.oldTime().oldTime().internalField()
611  *vf.oldTime().oldTime().internalField()*mesh().V00()
612  );
613  }
614  else
615  {
616  fvm.source() = rDeltaT*mesh().V()*
617  (
618  coefft0*rho.oldTime().internalField()
619  *vf.oldTime().internalField()
620  - coefft00*rho.oldTime().oldTime().internalField()
621  *vf.oldTime().oldTime().internalField()
622  );
623  }
624 
625  return tfvm;
626 }
627 
628 
630 (
631  const volScalarField& rA,
632  const volScalarField& U,
633  const surfaceScalarField& phi
634 )
635 {
637  (
638  "boundedBackwardDdtScheme::fvcDdtPhiCorr"
639  );
640 
641  return surfaceScalarField::null();
642 }
643 
644 
646 (
647  const volScalarField& rA,
648  const volScalarField& rho,
649  const volScalarField& U,
650  const surfaceScalarField& phi
651 )
652 {
654  (
655  "boundedBackwardDdtScheme::fvcDdtPhiCorr"
656  );
657 
658  return surfaceScalarField::null();
659 }
660 
661 
663 (
664  const volScalarField& vf
665 )
666 {
668  (
669  "boundedBackwardDdtScheme::meshPhi(const volScalarField& vf)"
670  );
671 
672  return surfaceScalarField::null();
673 }
674 
675 
676 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
677 
678 } // End namespace fv
679 
680 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
681 
682 } // End namespace Foam
683 
684 // ************************ vim: set sw=4 sts=4 et: ************************ //