Convection-diffusion steady intertial ddt scheme

src/finiteVolume/finiteVolume/ddtSchemes/steadyInertialDdtScheme/steadyInertialDdtScheme.C

@@ -26,6 +26,7 @@ License
 #include "steadyInertialDdtScheme.H"
 #include "surfaceInterpolate.H"
 #include "fvMatrices.H"
+#include "autoPtr.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -42,8 +43,52 @@ namespace fv
 template<class Type>
 tmp<volScalarField> steadyInertialDdtScheme<Type>::CorDeltaT() const
 {
+    const objectRegistry& registry = this->mesh();
+
+    autoPtr<surfaceScalarField> cofrDeltaTPtr;
+
+    if (registry.foundObject<surfaceScalarField>(phiName_))
+    {
+        cofrDeltaTPtr.reset(convectionCofrDeltaT().ptr());
+        Info<< "Convection min ddt: "
+            << gMin(1/cofrDeltaTPtr().internalField())
+            << endl;
+    }
+
+    if
+    (
+        registry.foundObject<volScalarField>(nuName_)
+     || registry.foundObject<surfaceScalarField>(nuName_)
+    )
+    {
+        if (cofrDeltaTPtr.valid())
+        {
+            cofrDeltaTPtr() =
+                Foam::max(cofrDeltaTPtr(), diffusionCofrDeltaT());
+            Info<< "Combined min ddt: "
+                << gMin(1/cofrDeltaTPtr().internalField())
+                << endl;
+        }
+        else
+        {
+            cofrDeltaTPtr.reset(diffusionCofrDeltaT().ptr());
+            Info<< "Diffusion min ddt: "
+                << gMin(1/cofrDeltaTPtr().internalField())
+                << endl;
+        }
+    }
+    
+    if (cofrDeltaTPtr.empty())
+    {
+        FatalErrorIn
+        (
+            "steaddyInertialDdtScheme<Type>::CorDeltaT() const"
+        )   << "Cannot find phi or nu: " << phiName_ << " " << nuName_
+            << abort(FatalError);
+    }
+
     // Collect face delta t and pick the smallest for the cell
-    surfaceScalarField cofrDeltaT = CofrDeltaT();
+    surfaceScalarField& cofrDeltaT = cofrDeltaTPtr();
 
     tmp<volScalarField> tcorDeltaT
     (
@@ -60,7 +105,6 @@ tmp<volScalarField> steadyInertialDdtScheme<Type>::CorDeltaT() const
             zeroGradientFvPatchScalarField::typeName
         )
     );
-
     volScalarField& corDeltaT = tcorDeltaT();
 
     const unallocLabelList& owner = mesh().owner();
@@ -100,7 +144,8 @@ tmp<volScalarField> steadyInertialDdtScheme<Type>::CorDeltaT() const
 
 
 template<class Type>
-tmp<surfaceScalarField> steadyInertialDdtScheme<Type>::CofrDeltaT() const
+tmp<surfaceScalarField>
+steadyInertialDdtScheme<Type>::convectionCofrDeltaT() const
 {
     const objectRegistry& registry = this->mesh();
 
@@ -143,8 +188,76 @@ tmp<surfaceScalarField> steadyInertialDdtScheme<Type>::CofrDeltaT() const
     }
     else
     {
-        FatalErrorIn("steaddyInertialDdtScheme<Type>::CofrDeltaT() const")
-            << "Incorrect dimensions of phi: " << phi.dimensions()
+        FatalErrorIn
+        (
+            "steaddyInertialDdtScheme<Type>::convectionCofrDeltaT() const"
+        )   << "Incorrect dimensions of phi: " << phi.dimensions()
+            << abort(FatalError);
+
+        return tmp<surfaceScalarField>(NULL);
+    }
+}
+
+
+template<class Type>
+tmp<surfaceScalarField>
+steadyInertialDdtScheme<Type>::diffusionCofrDeltaT() const
+{
+    const objectRegistry& registry = this->mesh();
+
+    if (registry.foundObject<volScalarField>(nuName_))
+    {
+        const volScalarField& nu =
+            registry.lookupObject<volScalarField>(nuName_);
+
+        return diffusionCofrDeltaT(fvc::interpolate(nu)());
+    }
+    else if (registry.foundObject<volScalarField>(nuName_))
+    {
+        const surfaceScalarField& nuf =
+            registry.lookupObject<surfaceScalarField>(nuName_);
+
+        return diffusionCofrDeltaT(nuf);
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "steaddyInertialDdtScheme<Type>::diffusionCofrDeltaT() const"
+        )   << "Cannot find nu"
+            << abort(FatalError);
+
+        return tmp<surfaceScalarField>(NULL);
+    }
+}
+
+template<class Type>
+tmp<surfaceScalarField>
+steadyInertialDdtScheme<Type>::diffusionCofrDeltaT
+(
+    const surfaceScalarField& nuf
+) const
+{
+    const objectRegistry& registry = this->mesh();
+
+    if (nuf.dimensions() == dimensionSet(0, 2, -1, 0, 0))
+    {
+        return nuf*sqr(mesh().surfaceInterpolation::deltaCoeffs())/maxCo_;
+    }
+    else if (nuf.dimensions() == dimensionSet(1, -1, -1, 0, 0))
+    {
+        const volScalarField& rho =
+            registry.lookupObject<volScalarField>(rhoName_);
+
+        return nuf*sqr(mesh().surfaceInterpolation::deltaCoeffs())/
+            (fvc::interpolate(rho)*maxCo_);
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "steaddyInertialDdtScheme<Type>::diffusionCofrDeltaT() const"
+        )   << "Incorrect dimensions of nu: " << nuf.dimensions()
             << abort(FatalError);
 
         return tmp<surfaceScalarField>(NULL);

src/finiteVolume/finiteVolume/ddtSchemes/steadyInertialDdtScheme/steadyInertialDdtScheme.H

@@ -25,7 +25,8 @@ Class
     Foam::fv::steadyInertialDdtScheme
 
 Description
-    Stabilised local time-step first-order Euler implicit/explicit ddt.
+    Stabilised local time-step first-order Euler implicit/explicit ddt
+    for intertially relaxed steady-state simulations.
     The time-step is adjusted locally so that an advective equations remains
     diagonally dominant.
 
@@ -35,6 +36,9 @@ Description
 
     See also CoEulerDdtScheme.
 
+Author
+    Hrvoje Jasak, Wikki Ltd.  All rights reserved.
+
 SourceFiles
     steadyInertialDdtScheme.C
 
@@ -66,8 +70,13 @@ class steadyInertialDdtScheme
 {
     // Private Data
 
-        //- Name of the flux field used to calculate the local time-step
-        word phiName_;
+        //- Name of the flux field used to calculate the local
+        //  convective time-step
+        word phiName_; 
+
+        //- Name of the diffusivity field used to calculate the local
+        //  diffusive time-step
+        word nuName_;
 
         //- Name of the density field used to obtain the volumetric flux
         //  from the mass flux if required
@@ -89,7 +98,19 @@ class steadyInertialDdtScheme
         tmp<volScalarField> CorDeltaT() const;
 
         //- Return the reciprocal of the face-Courant-number limited time-step
-        tmp<surfaceScalarField> CofrDeltaT() const;
+        //  for convective transport
+        tmp<surfaceScalarField> convectionCofrDeltaT() const;
+
+        //- Return the reciprocal of the face-Courant-number limited time-step
+        //  for diffusive transport
+        tmp<surfaceScalarField> diffusionCofrDeltaT() const;
+
+        //- Return the reciprocal of the face-Courant-number limited time-step
+        //  for diffusive transport ginev face diffusivity
+        tmp<surfaceScalarField> diffusionCofrDeltaT
+        (
+            const surfaceScalarField& nuf
+        ) const;
 
 
 public:
@@ -105,6 +126,7 @@ public:
         :
             ddtScheme<Type>(mesh, is),
             phiName_(is),
+            nuName_(is),
             rhoName_(is),
             maxCo_(readScalar(is))
         {}