applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H

   1 {
   2     const dictionary& pimpleDict = pimple.dict();
   3
   4     scalar maxCo
   5     (
   6         pimpleDict.lookupOrDefault<scalar>("maxCo", 0.9)
   7     );
   8
   9     scalar maxAlphaCo
  10     (
  11         pimpleDict.lookupOrDefault<scalar>("maxAlphaCo", 0.2)
  12     );
  13
  14     scalar rDeltaTSmoothingCoeff
  15     (
  16         pimpleDict.lookupOrDefault<scalar>("rDeltaTSmoothingCoeff", 0.1)
  17     );
  18
  19     label nAlphaSpreadIter
  20     (
  21         pimpleDict.lookupOrDefault<label>("nAlphaSpreadIter", 1)
  22     );
  23
  24     scalar alphaSpreadDiff
  25     (
  26         pimpleDict.lookupOrDefault<scalar>("alphaSpreadDiff", 0.2)
  27     );
  28
  29     scalar alphaSpreadMax
  30     (
  31         pimpleDict.lookupOrDefault<scalar>("alphaSpreadMax", 0.99)
  32     );
  33
  34     scalar alphaSpreadMin
  35     (
  36         pimpleDict.lookupOrDefault<scalar>("alphaSpreadMin", 0.01)
  37     );
  38
  39     label nAlphaSweepIter
  40     (
  41         pimpleDict.lookupOrDefault<label>("nAlphaSweepIter", 5)
  42     );
  43
  44     scalar rDeltaTDampingCoeff
  45     (
  46         pimpleDict.lookupOrDefault<scalar>("rDeltaTDampingCoeff", 1.0)
  47     );
  48
  49     scalar maxDeltaT
  50     (
  51         pimpleDict.lookupOrDefault<scalar>("maxDeltaT", GREAT)
  52     );
  53
  54     volScalarField rDeltaT0("rDeltaT0", rDeltaT);
  55
  56     // Set the reciprocal time-step from the local Courant number
  57     rDeltaT.dimensionedInternalField() = max
  58     (
  59         1/dimensionedScalar("maxDeltaT", dimTime, maxDeltaT),
  60         fvc::surfaceSum(mag(rhoPhi))().dimensionedInternalField()
  61        /((2*maxCo)*mesh.V()*rho.dimensionedInternalField())
  62     );
  63
  64     if (maxAlphaCo < maxCo)
  65     {
  66         // Further limit the reciprocal time-step
  67         // in the vicinity of the interface
  68
  69         volScalarField alpha1Bar(fvc::average(alpha1));
  70
  71         rDeltaT.dimensionedInternalField() = max
  72         (
  73             rDeltaT.dimensionedInternalField(),
  74             pos(alpha1Bar.dimensionedInternalField() - alphaSpreadMin)
  75            *pos(alphaSpreadMax - alpha1Bar.dimensionedInternalField())
  76            *fvc::surfaceSum(mag(phi))().dimensionedInternalField()
  77            /((2*maxAlphaCo)*mesh.V())
  78         );
  79     }
  80
  81     // Update tho boundary values of the reciprocal time-step
  82     rDeltaT.correctBoundaryConditions();
  83
  84     Info<< "Flow time scale min/max = "
  85         << gMin(1/rDeltaT.internalField())
  86         << ", " << gMax(1/rDeltaT.internalField()) << endl;
  87
  88     if (rDeltaTSmoothingCoeff < 1.0)
  89     {
  90         fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff);
  91     }
  92
  93     if (nAlphaSpreadIter > 0)
  94     {
  95         fvc::spread
  96         (
  97             rDeltaT,
  98             alpha1,
  99             nAlphaSpreadIter,
 100             alphaSpreadDiff,
 101             alphaSpreadMax,
 102             alphaSpreadMin
 103         );
 104     }
 105
 106     if (nAlphaSweepIter > 0)
 107     {
 108         fvc::sweep(rDeltaT, alpha1, nAlphaSweepIter, alphaSpreadDiff);
 109     }
 110
 111     Info<< "Smoothed flow time scale min/max = "
 112         << gMin(1/rDeltaT.internalField())
 113         << ", " << gMax(1/rDeltaT.internalField()) << endl;
 114
 115     // Limit rate of change of time scale
 116     // - reduce as much as required
 117     // - only increase at a fraction of old time scale
 118     if
 119     (
 120         rDeltaTDampingCoeff < 1.0
 121      && runTime.timeIndex() > runTime.startTimeIndex() + 1
 122     )
 123     {
 124         rDeltaT =
 125              rDeltaT0
 126             *max(rDeltaT/rDeltaT0, scalar(1.0) - rDeltaTDampingCoeff);
 127
 128         Info<< "Damped flow time scale min/max = "
 129             << gMin(1/rDeltaT.internalField())
 130             << ", " << gMax(1/rDeltaT.internalField()) << endl;
 131     }
 132
 133     label nAlphaSubCycles
 134     (
 135         readLabel(pimpleDict.lookup("nAlphaSubCycles"))
 136     );
 137
 138     rSubDeltaT = rDeltaT*nAlphaSubCycles;
 139 }