chemistryModel/laminar/PaSR: Add support for LTS to laminar model and to PaSR by derivation

src/combustionModels/PaSR/PaSR.C

@@ -35,7 +35,7 @@ Foam::combustionModels::PaSR<Type>::PaSR
     const fvMesh& mesh
 )
 :
-    Type(modelType, mesh),
+    laminar<Type>(modelType, mesh),
     Cmix_(readScalar(this->coeffs().lookup("Cmix"))),
     turbulentReaction_(this->coeffs().lookup("turbulentReaction")),
     kappa_
@@ -50,21 +50,8 @@ Foam::combustionModels::PaSR<Type>::PaSR
         ),
         mesh,
         dimensionedScalar("kappa", dimless, 0.0)
-    ),
-    integrateReactionRate_
-    (
-        this->coeffs().lookupOrDefault("integrateReactionRate", true)
     )
-{
+{}
-    if (integrateReactionRate_)
-    {
-        Info<< "    using integrated reaction rate" << endl;
-    }
-    else
-    {
-        Info<< "    using instantaneous reaction rate" << endl;
-    }
-}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
@@ -76,28 +63,12 @@ Foam::combustionModels::PaSR<Type>::~PaSR()
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-template<class Type>
-Foam::tmp<Foam::volScalarField> Foam::combustionModels::PaSR<Type>::tc() const
-{
-    return this->chemistryPtr_->tc();
-}
-
-
 template<class Type>
 void Foam::combustionModels::PaSR<Type>::correct()
 {
     if (this->active())
     {
-        const scalar dt = this->mesh().time().deltaTValue();
+        laminar<Type>::correct();
-
-        if (integrateReactionRate_)
-        {
-            this->chemistryPtr_->solve(dt);
-        }
-        else
-        {
-            this->chemistryPtr_->calculate();
-        }
 
         if (turbulentReaction_)
         {
@@ -108,7 +79,7 @@ void Foam::combustionModels::PaSR<Type>::correct()
             tmp<volScalarField> tmuEff(this->turbulence().muEff());
             const volScalarField& muEff = tmuEff();
 
-            tmp<volScalarField> ttc(tc());
+            tmp<volScalarField> ttc(this->tc());
             const volScalarField& tc = ttc();
 
             forAll(epsilon, i)
@@ -138,18 +109,7 @@ template<class Type>
 Foam::tmp<Foam::fvScalarMatrix>
 Foam::combustionModels::PaSR<Type>::R(volScalarField& Y) const
 {
-    tmp<fvScalarMatrix> tSu(new fvScalarMatrix(Y, dimMass/dimTime));
+    return kappa_*laminar<Type>::R(Y);
-
-    fvScalarMatrix& Su = tSu();
-
-    if (this->active())
-    {
-        const label specieI = this->thermo().composition().species()[Y.name()];
-
-        Su += kappa_*this->chemistryPtr_->RR(specieI);
-    }
-
-    return tSu;
 }
 
 
@@ -157,32 +117,7 @@ template<class Type>
 Foam::tmp<Foam::volScalarField>
 Foam::combustionModels::PaSR<Type>::dQ() const
 {
-    tmp<volScalarField> tdQ
+    return kappa_*laminar<Type>::dQ();
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                typeName + ":dQ",
-                this->mesh().time().timeName(),
-                this->mesh(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            this->mesh(),
-            dimensionedScalar("dQ", dimEnergy/dimTime, 0.0),
-            zeroGradientFvPatchScalarField::typeName
-        )
-    );
-
-    if (this->active())
-    {
-        volScalarField& dQ = tdQ();
-        dQ = kappa_*this->chemistryPtr_->dQ();
-    }
-
-    return tdQ;
 }
 
 
@@ -190,44 +125,17 @@ template<class Type>
 Foam::tmp<Foam::volScalarField>
 Foam::combustionModels::PaSR<Type>::Sh() const
 {
-    tmp<volScalarField> tSh
+    return kappa_*laminar<Type>::Sh();
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                typeName + ":Sh",
-                this->mesh().time().timeName(),
-                this->mesh(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            this->mesh(),
-            dimensionedScalar("zero", dimEnergy/dimTime/dimVolume, 0.0),
-            zeroGradientFvPatchScalarField::typeName
-        )
-    );
-
-    if (this->active())
-    {
-        scalarField& Sh = tSh();
-        Sh = kappa_*this->chemistryPtr_->Sh();
-    }
-
-    return tSh;
 }
 
 
 template<class Type>
 bool Foam::combustionModels::PaSR<Type>::read()
 {
-    if (Type::read())
+    if (laminar<Type>::read())
     {
         this->coeffs().lookup("Cmix") >> Cmix_;
         this->coeffs().lookup("turbulentReaction") >> turbulentReaction_;
-        this->coeffs().lookup("integrateReactionRate")
-            >> integrateReactionRate_;
         return true;
     }
     else

src/combustionModels/PaSR/PaSR.H

@@ -38,6 +38,8 @@ SourceFiles
 #ifndef PaSR_H
 #define PaSR_H
 
+#include "laminar.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -52,7 +54,7 @@ namespace combustionModels
 template<class Type>
 class PaSR
 :
-    public Type
+    public laminar<Type>
 {
     // Private data
 
@@ -65,16 +67,9 @@ class PaSR
         //- Mixing parameter
         volScalarField kappa_;
 
-        //- Integrate reaction rate over the time-step
-        //  using the selected ODE solver
-        bool integrateReactionRate_;
-
 
     // Private Member Functions
 
-        //- Return the chemical time scale
-        tmp<volScalarField> tc() const;
-
         //- Disallow copy construct
         PaSR(const PaSR&);
 

src/combustionModels/laminar/laminar.C

@@ -25,6 +25,7 @@ License
 
 #include "laminar.H"
 #include "fvmSup.H"
+#include "localEulerDdtScheme.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -76,7 +77,37 @@ void Foam::combustionModels::laminar<Type>::correct()
     {
         if (integrateReactionRate_)
         {
-            this->chemistryPtr_->solve(this->mesh().time().deltaTValue());
+            word ddtScheme(this->mesh().ddtScheme("Yi"));
+
+            if (ddtScheme == fv::localEulerDdtScheme<scalar>::typeName)
+            {
+                const scalarField& rDeltaT =
+                this->mesh().objectRegistry::lookupObject<volScalarField>
+                (
+                    "rDeltaT"
+                );
+
+                if (this->coeffs().found("maxIntegrationTime"))
+                {
+                    scalar maxIntegrationTime
+                    (
+                        readScalar(this->coeffs().lookup("maxIntegrationTime"))
+                    );
+
+                    this->chemistryPtr_->solve
+                    (
+                        min(1.0/rDeltaT, maxIntegrationTime)()
+                    );
+                }
+                else
+                {
+                    this->chemistryPtr_->solve((1.0/rDeltaT)());
+                }
+            }
+            else
+            {
+                this->chemistryPtr_->solve(this->mesh().time().deltaTValue());
+            }
         }
         else
         {

src/combustionModels/laminar/laminar.H

@@ -57,12 +57,17 @@ class laminar
         //  using the selected ODE solver
         bool integrateReactionRate_;
 
+protected:
 
-    // Private Member Functions
+    // Protected Member Functions
 
         //- Return the chemical time scale
         tmp<volScalarField> tc() const;
 
+private:
+
+    // Private Member Functions
+
         //- Disallow copy construct
         laminar(const laminar&);
 

src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H

@@ -156,6 +156,10 @@ public:
                 //  and return the characteristic time
                 virtual scalar solve(const scalar deltaT) = 0;
 
+                //- Solve the reaction system for the given time step
+                //  and return the characteristic time
+                virtual scalar solve(const scalarField& deltaT) = 0;
+
                 //- Return the chemical time scale
                 virtual tmp<volScalarField> tc() const = 0;
 

src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C

@@ -25,6 +25,7 @@ License
 
 #include "chemistryModel.H"
 #include "reactingMixture.H"
+#include "UniformField.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -742,9 +743,10 @@ void Foam::chemistryModel<CompType, ThermoType>::calculate()
 
 
 template<class CompType, class ThermoType>
+template<class DeltaTType>
 Foam::scalar Foam::chemistryModel<CompType, ThermoType>::solve
 (
-    const scalar deltaT
+    const DeltaTType& deltaT
 )
 {
     CompType::correct();
@@ -795,9 +797,9 @@ Foam::scalar Foam::chemistryModel<CompType, ThermoType>::solve
 
         // initialise timing parameters
         scalar t = 0;
+        scalar timeLeft = deltaT[celli];
         scalar tauC = this->deltaTChem_[celli];
-        scalar dt = min(deltaT, tauC);
+        scalar dt = min(timeLeft, tauC);
-        scalar timeLeft = deltaT;
 
         // calculate the chemical source terms
         while (timeLeft > SMALL)
@@ -823,17 +825,39 @@ Foam::scalar Foam::chemistryModel<CompType, ThermoType>::solve
         dc = c - c0;
         for (label i=0; i<nSpecie_; i++)
         {
-            RR_[i][celli] = dc[i]*specieThermo_[i].W()/deltaT;
+            RR_[i][celli] = dc[i]*specieThermo_[i].W()/deltaT[celli];
         }
     }
 
-    // Don't allow the time-step to change more than a factor of 2
-    deltaTMin = min(deltaTMin, 2*deltaT);
-
     return deltaTMin;
 }
 
 
+template<class CompType, class ThermoType>
+Foam::scalar Foam::chemistryModel<CompType, ThermoType>::solve
+(
+    const scalar deltaT
+)
+{
+    // Don't allow the time-step to change more than a factor of 2
+    return min
+    (
+        this->solve<UniformField<scalar> >(UniformField<scalar>(deltaT)),
+        2*deltaT
+    );
+}
+
+
+template<class CompType, class ThermoType>
+Foam::scalar Foam::chemistryModel<CompType, ThermoType>::solve
+(
+    const scalarField& deltaT
+)
+{
+    return this->solve<scalarField>(deltaT);
+}
+
+
 template<class CompType, class ThermoType>
 Foam::scalar Foam::chemistryModel<CompType, ThermoType>::solve
 (

src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H

@@ -70,6 +70,11 @@ class chemistryModel
         //- Disallow default bitwise assignment
         void operator=(const chemistryModel&);
 
+        //- Solve the reaction system for the given time step
+        //  of given type and return the characteristic time
+        template<class DeltaTType>
+        scalar solve(const DeltaTType& deltaT);
+
 
 protected:
 
@@ -221,6 +226,10 @@ public:
             //  and return the characteristic time
             virtual scalar solve(const scalar deltaT);
 
+            //- Solve the reaction system for the given time step
+            //  and return the characteristic time
+            virtual scalar solve(const scalarField& deltaT);
+
             //- Return the chemical time scale
             virtual tmp<volScalarField> tc() const;