ENH: Added semi-implicit treatment for particle sources

src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C

@@ -126,6 +126,16 @@ Foam::scalar Foam::KinematicCloud<ParcelType>::cloudSolution::relaxCoeff
 }
 
 
+template<class ParcelType>
+bool Foam::KinematicCloud<ParcelType>::cloudSolution::semiImplicit
+(
+    const word& fieldName
+) const
+{
+    return readBool(sourceTermDict().subDict(fieldName).lookup("semiImplicit"));
+}
+
+
 template<class ParcelType>
 bool Foam::KinematicCloud<ParcelType>::cloudSolution::sourceActive() const
 {
@@ -495,6 +505,22 @@ Foam::KinematicCloud<ParcelType>::KinematicCloud
             mesh_,
             dimensionedVector("zero", dimMass*dimVelocity, vector::zero)
         )
+    ),
+    UCoeff_
+    (
+        new DimensionedField<scalar, volMesh>
+        (
+            IOobject
+            (
+                this->name() + "UCoeff",
+                this->db().time().timeName(),
+                this->db(),
+                IOobject::READ_IF_PRESENT,
+                IOobject::AUTO_WRITE
+            ),
+            mesh_,
+            dimensionedScalar("zero",  dimMass/dimTime, 0.0)
+        )
     )
 {
     if (readFields)
@@ -556,7 +582,24 @@ Foam::KinematicCloud<ParcelType>::KinematicCloud
             c.UTrans_().dimensions(),
             c.UTrans_().field()
         )
+    ),
+    UCoeff_
+    (
+        new DimensionedField<scalar, volMesh>
+        (
+            IOobject
+            (
+                name + "UCoeff",
+                this->db().time().timeName(),
+                this->db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            c.UCoeff_()
         )
+    )
+
 {}
 
 
@@ -602,7 +645,8 @@ Foam::KinematicCloud<ParcelType>::KinematicCloud
     postProcessingModel_(NULL),
     surfaceFilmModel_(NULL),
     UIntegrator_(NULL),
-    UTrans_(NULL)
+    UTrans_(NULL),
+    UCoeff_(NULL)
 {}
 
 
@@ -658,6 +702,7 @@ template<class ParcelType>
 void Foam::KinematicCloud<ParcelType>::resetSourceTerms()
 {
     UTrans().field() = vector::zero;
+    UCoeff().field() = 0.0;
 }
 
 

src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H

@@ -187,6 +187,9 @@ public:
                 //- Return relaxation coefficient for field
                 scalar relaxCoeff(const word& fieldName) const;
 
+                //- Return semi-implicit flag coefficient for field
+                bool semiImplicit(const word& fieldName) const;
+
                 //- Return reference to the mesh
                 inline const fvMesh& mesh() const;
 
@@ -323,6 +326,9 @@ protected:
             //- Momentum
             autoPtr<DimensionedField<vector, volMesh> > UTrans_;
 
+            //- Coefficient for carrier phase U equation
+            autoPtr<DimensionedField<scalar, volMesh> > UCoeff_;
+
 
         // Cloud evolution functions
 
@@ -536,8 +542,15 @@ public:
                     inline const DimensionedField<vector, volMesh>&
                         UTrans() const;
 
-                    //- Return tmp momentum source term - fully explicit
+                     //- Return coefficient for carrier phase U equation
-                    inline tmp<DimensionedField<vector, volMesh> > SU() const;
+                    inline DimensionedField<scalar, volMesh>& UCoeff();
+
+                    //- Return const coefficient for carrier phase U equation
+                    inline const DimensionedField<scalar, volMesh>&
+                        UCoeff() const;
+
+                    //- Return tmp momentum source term
+                    inline tmp<fvVectorMatrix> SU(volVectorField& U) const;
 
 
         // Check

src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H

@@ -410,38 +410,54 @@ Foam::KinematicCloud<ParcelType>::UTrans() const
 
 
 template<class ParcelType>
-inline Foam::tmp<Foam::DimensionedField<Foam::vector, Foam::volMesh> >
+inline Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
-Foam::KinematicCloud<ParcelType>::SU() const
+Foam::KinematicCloud<ParcelType>::UCoeff()
 {
-    tmp<DimensionedField<vector, volMesh> > tSU
+    return UCoeff_();
-    (
+}
-        new DimensionedField<vector, volMesh>
+
-        (
+
-            IOobject
+template<class ParcelType>
-            (
+inline const Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
-                this->name() + "SU",
+Foam::KinematicCloud<ParcelType>::UCoeff() const
-                this->db().time().timeName(),
+{
-                this->mesh(),
+    return UCoeff_();
-                IOobject::NO_READ,
+}
-                IOobject::NO_WRITE
+
-            ),
+
-            this->mesh(),
+template<class ParcelType>
-            dimensionedVector
+inline Foam::tmp<Foam::fvVectorMatrix>
-            (
+Foam::KinematicCloud<ParcelType>::SU(volVectorField& U) const
-                "zero",
+{
-                dimDensity*dimVelocity/dimTime,
+    if (debug)
-                vector::zero
+    {
-            )
+        Info<< "UTrans min/max = " << min(UTrans()).value() << ", "
-        )
+            << max(UTrans()).value() << nl
-    );
+            << "UCoeff min/max = " << min(UCoeff()).value() << ", "
+            << max(UCoeff()).value() << endl;
+    }
 
     if (solution_.sourceActive())
     {
-        vectorField& SU = tSU().field();
+        if (solution_.semiImplicit("UTrans"))
-        SU = UTrans()/(mesh_.V()*this->db().time().deltaT());
+        {
+            return
+                UTrans()/(mesh_.V()*this->db().time().deltaT())
+              + fvm::SuSp(-UCoeff().dimensionedInternalField()/mesh_.V(), U)
+              + UCoeff()/mesh_.V()*U;
+        }
+        else
+        {
+            tmp<fvVectorMatrix> tfvm(new fvVectorMatrix(U, dimForce));
+            fvVectorMatrix& fvm = tfvm();
+
+            fvm.source() = -UTrans()/(this->db().time().deltaT());
+
+            return tfvm;
+        }
     }
 
-    return tSU;
+    return tmp<fvVectorMatrix>(new fvVectorMatrix(U, dimForce));
 }
 
 

src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C

@@ -103,7 +103,25 @@ Foam::ThermoCloud<ParcelType>::ThermoCloud
             this->mesh(),
             dimensionedScalar("zero", dimEnergy, 0.0)
         )
+    ),
+    hsCoeff_
+    (
+        new DimensionedField<scalar, volMesh>
+        (
+            IOobject
+            (
+                this->name() + "hsCoeff",
+                this->db().time().timeName(),
+                this->db(),
+                IOobject::READ_IF_PRESENT,
+                IOobject::AUTO_WRITE
+            ),
+            this->mesh(),
+            dimensionedScalar("zero", dimEnergy/dimTime/dimTemperature, 0.0),
+            zeroGradientFvPatchScalarField::typeName
         )
+    )
+
 {
     if (readFields)
     {
@@ -149,6 +167,22 @@ Foam::ThermoCloud<ParcelType>::ThermoCloud
             ),
             c.hsTrans()
         )
+    ),
+    hsCoeff_
+    (
+        new DimensionedField<scalar, volMesh>
+        (
+            IOobject
+            (
+                this->name() + "hsCoeff",
+                this->db().time().timeName(),
+                this->db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            c.hsCoeff()
+        )
     )
 {}
 
@@ -171,7 +205,8 @@ Foam::ThermoCloud<ParcelType>::ThermoCloud
     heatTransferModel_(NULL),
     TIntegrator_(NULL),
     radiation_(false),
-    hsTrans_(NULL)
+    hsTrans_(NULL),
+    hsCoeff_(NULL)
 {}
 
 
@@ -233,6 +268,7 @@ void Foam::ThermoCloud<ParcelType>::resetSourceTerms()
 {
     KinematicCloud<ParcelType>::resetSourceTerms();
     hsTrans_->field() = 0.0;
+    hsCoeff_->field() = 0.0;
 }
 
 

src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H

@@ -122,6 +122,9 @@ protected:
             //- Sensible enthalpy transfer [J/kg]
             autoPtr<DimensionedField<scalar, volMesh> > hsTrans_;
 
+            //- Coefficient for carrier phase hs equation [W/K]
+            autoPtr<DimensionedField<scalar, volMesh> > hsCoeff_;
+
 
     // Protected Member Functions
 
@@ -232,8 +235,15 @@ public:
                     inline const DimensionedField<scalar, volMesh>&
                         hsTrans() const;
 
-                    //- Return enthalpy source [J/kg/m3/s]
+                    //- Return coefficient for carrier phase hs equation
-                    inline tmp<DimensionedField<scalar, volMesh> > Sh() const;
+                    inline DimensionedField<scalar, volMesh>& hsCoeff();
+
+                    //- Return const coefficient for carrier phase hs equation
+                    inline const DimensionedField<scalar, volMesh>&
+                        hsCoeff() const;
+
+                    //- Return sensible enthalpy source term [J/kg/m3/s]
+                    inline tmp<fvScalarMatrix> Sh(volScalarField& hs) const;
 
 
                 // Radiation - overrides thermoCloud virtual abstract members

src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H

@@ -98,34 +98,60 @@ Foam::ThermoCloud<ParcelType>::hsTrans() const
 
 
 template<class ParcelType>
-inline Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
+inline Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
-Foam::ThermoCloud<ParcelType>::Sh() const
+Foam::ThermoCloud<ParcelType>::hsCoeff()
 {
-    tmp<DimensionedField<scalar, volMesh> > tSh
+    return hsCoeff_();
-    (
+}
-        new DimensionedField<scalar, volMesh>
+
-        (
+
-            IOobject
+template<class ParcelType>
-            (
+inline const Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
-                this->name() + "Sh",
+Foam::ThermoCloud<ParcelType>::hsCoeff() const
-                this->db().time().timeName(),
+{
-                this->mesh(),
+    return hsCoeff_();
-                IOobject::NO_READ,
+}
-                IOobject::NO_WRITE,
+
-                false
+
-            ),
+template<class ParcelType>
-            this->mesh(),
+inline Foam::tmp<Foam::fvScalarMatrix>
-            dimensionedScalar("zero", dimEnergy/dimVolume/dimTime, 0.0)
+Foam::ThermoCloud<ParcelType>::Sh(volScalarField& hs) const
-        )
+{
-    );
+    if (debug)
+    {
+        Info<< "hsTrans min/max = " << min(hsTrans()).value() << ", "
+            << max(hsTrans()).value() << nl
+            << "hsCoeff min/max = " << min(hsCoeff()).value() << ", "
+            << max(hsCoeff()).value() << endl;
+    }
 
     if (this->solution().sourceActive())
     {
-        scalarField& Sh = tSh();
+        if (this->solution().semiImplicit("hsTrans"))
-        Sh = hsTrans_()/(this->mesh().V()*this->db().time().deltaT());
+        {
+            const volScalarField Cp = thermo_.thermo().Cp();
+
+            return
+                hsTrans()/(this->mesh().V()*this->db().time().deltaT())
+              + fvm::SuSp
+                (
+                   -hsCoeff().dimensionedInternalField()/(Cp*this->mesh().V()),
+                    hs
+                )
+              + hsCoeff()/(Cp*this->mesh().V())*hs;
+        }
+        else
+        {
+            tmp<fvScalarMatrix> tfvm(new fvScalarMatrix(hs, dimEnergy/dimTime));
+            fvScalarMatrix& fvm = tfvm();
+
+            fvm.source() = -hsTrans()/(this->db().time().deltaT());
+
+            return tfvm;
+        }
     }
 
-    return tSh;
+    return tmp<fvScalarMatrix>(new fvScalarMatrix(hs, dimEnergy/dimTime));
 }
 
 

src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C

@@ -121,8 +121,23 @@ void Foam::KinematicParcel<ParcelType>::calc
     // ~~~~~~
 
     // Calculate new particle velocity
+    scalar Cud = 0.0;
     vector U1 =
-        calcVelocity(td, dt, cellI, Re, muc_, d0, U0, rho0, mass0, Su, dUTrans);
+        calcVelocity
+        (
+            td,
+            dt,
+            cellI,
+            Re,
+            muc_,
+            d0,
+            U0,
+            rho0,
+            mass0,
+            Su,
+            dUTrans,
+            Cud
+        );
 
 
     // Accumulate carrier phase source terms
@@ -131,6 +146,9 @@ void Foam::KinematicParcel<ParcelType>::calc
     {
         // Update momentum transfer
         td.cloud().UTrans()[cellI] += np0*dUTrans;
+
+        // Update momentum transfer coefficient
+        td.cloud().UCoeff()[cellI] += np0*mass0*Cud;
     }
 
 
@@ -154,7 +172,8 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
     const scalar rho,
     const scalar mass,
     const vector& Su,
-    vector& dUTrans
+    vector& dUTrans,
+    scalar& Cud
 ) const
 {
     const polyMesh& mesh = this->cloud().pMesh();
@@ -165,7 +184,7 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
     tetIndices tetIs = this->currentTetIndices();
 
     // Momentum source due to particle forces
-    const vector FCoupled = mass*td.cloud().forces().calcCoupled
+    const vector Fcp = mass*td.cloud().forces().calcCoupled
     (
         this->position(),
         tetIs,
@@ -177,7 +196,7 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
         d
     );
 
-    const vector FNonCoupled = mass*td.cloud().forces().calcNonCoupled
+    const vector Fncp = mass*td.cloud().forces().calcNonCoupled
     (
         this->position(),
         tetIs,
@@ -195,15 +214,17 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
 
     // Update velocity - treat as 3-D
     const scalar As = this->areaS(d);
-    const vector ap = Uc_ + (FCoupled + FNonCoupled + Su)/(utc*As);
+    const vector ap = Uc_ + (Fcp + Fncp + Su)/(utc*As);
     const scalar bp = 6.0*utc/(rho*d);
 
+    Cud = bp;
+
     IntegrationScheme<vector>::integrationResult Ures =
         td.cloud().UIntegrator().integrate(U, dt, ap, bp);
 
     vector Unew = Ures.value();
 
-    dUTrans += dt*(utc*As*(Ures.average() - Uc_) - FCoupled);
+    dUTrans += dt*(utc*As*(Ures.average() - Uc_) - Fcp);
 
     // Apply correction to velocity and dUTrans for reduced-D cases
     meshTools::constrainDirection(mesh, mesh.solutionD(), Unew);

src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H

@@ -308,7 +308,8 @@ protected:
             const scalar rho,          // density
             const scalar mass,         // mass
             const vector& Su,          // explicit particle momentum source
-            vector& dUTrans            // momentum transfer to carrier
+            vector& dUTrans,           // momentum transfer to carrier
+            scalar& Cud                // linearised drag coefficient
         ) const;
 
 

src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C

@@ -353,6 +353,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     // ~~~~~~~~~~~~~
 
     // Calculate new particle temperature
+    scalar Cuh = 0.0;
     scalar T1 =
         calcHeatTransfer
         (
@@ -368,7 +369,8 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
             Cp0,
             NCpW,
             Sh,
-            dhsTrans
+            dhsTrans,
+            Cuh
         );
 
 
@@ -376,8 +378,23 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     // ~~~~~~
 
     // Calculate new particle velocity
+    scalar Cud = 0;
     vector U1 =
-        calcVelocity(td, dt, cellI, Re, mus, d0, U0, rho0, mass0, Su, dUTrans);
+        calcVelocity
+        (
+            td,
+            dt,
+            cellI,
+            Re,
+            mus,
+            d0,
+            U0,
+            rho0,
+            mass0,
+            Su,
+            dUTrans,
+            Cud
+        );
 
     dUTrans += 0.5*(mass0 - mass1)*(U0 + U1);
 
@@ -414,8 +431,14 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         // Update momentum transfer
         td.cloud().UTrans()[cellI] += np0*dUTrans;
 
+        // Update momentum transfer coefficient
+        td.cloud().UCoeff()[cellI] += np0*0.5*(mass0 + mass1)*Cud;
+
         // Update sensible enthalpy transfer
         td.cloud().hsTrans()[cellI] += np0*dhsTrans;
+
+        // Update sensible enthalpy coefficient
+        td.cloud().hsCoeff()[cellI] += np0*Cuh*this->areaS();
     }
 
 

src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C

@@ -303,6 +303,7 @@ void Foam::ReactingParcel<ParcelType>::calc
     // ~~~~~~~~~~~~~
 
     // Calculate new particle temperature
+    scalar Cuh = 0.0;
     scalar T1 =
         calcHeatTransfer
         (
@@ -318,7 +319,8 @@ void Foam::ReactingParcel<ParcelType>::calc
             Cp0,
             NCpW,
             Sh,
-            dhsTrans
+            dhsTrans,
+            Cuh
         );
 
 
@@ -326,8 +328,23 @@ void Foam::ReactingParcel<ParcelType>::calc
     // ~~~~~~
 
     // Calculate new particle velocity
+    scalar Cud = 0.0;
     vector U1 =
-        calcVelocity(td, dt, cellI, Re, mus, d0, U0, rho0, mass0, Su, dUTrans);
+        calcVelocity
+        (
+            td,
+            dt,
+            cellI,
+            Re,
+            mus,
+            d0,
+            U0,
+            rho0,
+            mass0,
+            Su,
+            dUTrans,
+            Cud
+        );
 
     dUTrans += 0.5*(mass0 - mass1)*(U0 + U1);
 
@@ -345,8 +362,14 @@ void Foam::ReactingParcel<ParcelType>::calc
         // Update momentum transfer
         td.cloud().UTrans()[cellI] += np0*dUTrans;
 
+        // Update momentum transfer coefficient
+        td.cloud().UCoeff()[cellI] += np0*0.5*(mass0 + mass1)*Cud;
+
         // Update sensible enthalpy transfer
         td.cloud().hsTrans()[cellI] += np0*dhsTrans;
+
+        // Update sensible enthalpy coefficient
+        td.cloud().hsCoeff()[cellI] += np0*Cuh*this->areaS();
     }
 
 

src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C

@@ -164,6 +164,7 @@ void Foam::ThermoParcel<ParcelType>::calc
     scalar NCpW = 0.0;
 
     // Calculate new particle velocity
+    scalar Cuh = 0.0;
     scalar T1 =
         calcHeatTransfer
         (
@@ -179,7 +180,8 @@ void Foam::ThermoParcel<ParcelType>::calc
             Cp0,
             NCpW,
             Sh,
-            dhsTrans
+            dhsTrans,
+            Cuh
         );
 
 
@@ -187,8 +189,23 @@ void Foam::ThermoParcel<ParcelType>::calc
     // ~~~~~~
 
     // Calculate new particle velocity
+    scalar Cud = 0.0;
     vector U1 =
-        calcVelocity(td, dt, cellI, Re, mus, d0, U0, rho0, mass0, Su, dUTrans);
+        calcVelocity
+        (
+            td,
+            dt,
+            cellI,
+            Re,
+            mus,
+            d0,
+            U0,
+            rho0,
+            mass0,
+            Su,
+            dUTrans,
+            Cud
+        );
 
 
     //  Accumulate carrier phase source terms
@@ -198,8 +215,14 @@ void Foam::ThermoParcel<ParcelType>::calc
         // Update momentum transfer
         td.cloud().UTrans()[cellI] += np0*dUTrans;
 
+        // Update momentum transfer coefficient
+        td.cloud().UCoeff()[cellI] += np0*mass0*Cud;
+
         // Update sensible enthalpy transfer
         td.cloud().hsTrans()[cellI] += np0*dhsTrans;
+
+        // Update sensible enthalpy coefficient
+        td.cloud().hsCoeff()[cellI] += np0*Cuh*this->areaS();
     }
 
     // Set new particle properties
@@ -225,7 +248,8 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
     const scalar Cp,
     const scalar NCpW,
     const scalar Sh,
-    scalar& dhsTrans
+    scalar& dhsTrans,
+    scalar& Cuh
 )
 {
     if (!td.cloud().heatTransfer().active())
@@ -270,6 +294,8 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
 
     dhsTrans += dt*htc*As*(0.5*(T + Tnew) - Tc_);
 
+    Cuh = bp;
+
     return Tnew;
 }
 

src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H

@@ -223,7 +223,8 @@ protected:
             const scalar Cp,           // specific heat capacity
             const scalar NCpW,         // Sum of N*Cp*W of emission species
             const scalar Sh,           // explicit particle enthalpy source
-            scalar& dhsTrans           // sensible enthalpy transfer to carrier
+            scalar& dhsTrans,          // sensible enthalpy transfer to carrier
+            scalar& Cuh                // linearised heat transfer coefficient
         );