ENH: Update to pyrolysis model, adding solid physical types and solid reaction solver types,

fixing bug in hExponentialThermo. Updating oppositeBurningPanels tutorial.

src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C

@@ -57,8 +57,10 @@ void reactingOneDim::readReactingOneDimControls()
     solution.lookup("nNonOrthCorr") >> nNonOrthCorr_;
     time().controlDict().lookup("maxDi") >> maxDiff_;
 
-    coeffs().lookup("radFluxName") >> primaryRadFluxName_;
     coeffs().lookup("minimumDelta") >> minimumDelta_;
+
+    coeffs().lookup("gasHSource") >> gasHSource_;
+    coeffs().lookup("QrHSource") >> QrHSource_;
 }
 
 
@@ -90,6 +92,58 @@ bool reactingOneDim::read(const dictionary& dict)
 }
 
 
+void reactingOneDim::updateQr()
+{
+    // Update local Qr from coupled Qr field
+    Qr_ == dimensionedScalar("zero", Qr_.dimensions(), 0.0);
+
+    // Retrieve field from coupled region using mapped boundary conditions
+    Qr_.correctBoundaryConditions();
+
+    forAll(intCoupledPatchIDs_, i)
+    {
+        const label patchI = intCoupledPatchIDs_[i];
+
+        scalarField& Qrp = Qr_.boundaryField()[patchI];
+
+        // Qr is positive going in the solid
+        // If the surface is emitting the radiative flux is set to zero
+        Qrp = max(Qrp, scalar(0.0));
+    }
+
+    const vectorField& cellC = regionMesh().cellCentres();
+
+    tmp<volScalarField> kappa = kappaRad();
+
+    // Propagate Qr through 1-D regions
+    label localPyrolysisFaceI = 0;
+    forAll(intCoupledPatchIDs_, i)
+    {
+        const label patchI = intCoupledPatchIDs_[i];
+
+        const scalarField& Qrp = Qr_.boundaryField()[patchI];
+        const vectorField& Cf = regionMesh().Cf().boundaryField()[patchI];
+
+        forAll(Qrp, faceI)
+        {
+            const scalar Qr0 = Qrp[faceI];
+            point Cf0 = Cf[faceI];
+            const labelList& cells = boundaryFaceCells_[localPyrolysisFaceI++];
+            scalar kappaInt = 0.0;
+            forAll(cells, k)
+            {
+                const label cellI = cells[k];
+                const point& Cf1 = cellC[cellI];
+                const scalar delta = mag(Cf1 - Cf0);
+                kappaInt += kappa()[cellI]*delta;
+                Qr_[cellI] = Qr0*exp(-kappaInt);
+                Cf0 = Cf1;
+            }
+        }
+    }
+}
+
+
 void reactingOneDim::updatePhiGas()
 {
     phiHsGas_ ==  dimensionedScalar("zero", phiHsGas_.dimensions(), 0.0);
@@ -112,20 +166,21 @@ void reactingOneDim::updatePhiGas()
         {
             const label patchI = intCoupledPatchIDs_[i];
 
-            const scalarField& phiGasp = phiHsGas_.boundaryField()[patchI];
+            scalarField& phiGasp = phiGas_.boundaryField()[patchI];
+            const scalarField& cellVol = regionMesh().V();
 
             forAll(phiGasp, faceI)
             {
-                const labelList& cells = boundaryFaceCells_[totalFaceId];
+                const labelList& cells = boundaryFaceCells_[totalFaceId++];
                 scalar massInt = 0.0;
                 forAllReverse(cells, k)
                 {
                     const label cellI = cells[k];
-                    massInt += RRiGas[cellI]*regionMesh().V()[cellI];
+                    massInt += RRiGas[cellI]*cellVol[cellI];
                     phiHsGas_[cellI] += massInt*HsiGas[cellI];
                 }
 
-                phiGas_.boundaryField()[patchI][faceI] += massInt;
+                phiGasp[faceI] += massInt;
 
                 if (debug)
                 {
@@ -136,7 +191,6 @@ void reactingOneDim::updatePhiGas()
                         << " is : " << massInt
                         << " [kg/s] " << endl;
                 }
-                totalFaceId ++;
             }
         }
         tHsiGas().clear();
@@ -146,6 +200,11 @@ void reactingOneDim::updatePhiGas()
 
 void reactingOneDim::updateFields()
 {
+    if (QrHSource_)
+    {
+        updateQr();
+    }
+
     updatePhiGas();
 }
 
@@ -183,8 +242,6 @@ void reactingOneDim::solveContinuity()
         Info<< "reactingOneDim::solveContinuity()" << endl;
     }
 
-    if (moveMesh_)
-    {
     const scalarField mass0 = rho_*regionMesh().V();
 
     fvScalarMatrix rhoEqn
@@ -194,11 +251,19 @@ void reactingOneDim::solveContinuity()
       - solidChemistry_->RRg()
     );
 
+    if (regionMesh().moving())
+    {
+        surfaceScalarField phiRhoMesh
+        (
+            fvc::interpolate(rho_)*regionMesh().phi()
+        );
+
+        rhoEqn += fvc::div(phiRhoMesh);
+    }
+
     rhoEqn.solve();
 
     updateMesh(mass0);
-
-    }
 }
 
 
@@ -239,6 +304,7 @@ void reactingOneDim::solveSpeciesMass()
     }
 
     Ys_[Ys_.size() - 1] = 1.0 - Yt;
+
 }
 
 
@@ -259,6 +325,18 @@ void reactingOneDim::solveEnergy()
         chemistrySh_
     );
 
+    if (gasHSource_)
+    {
+        const surfaceScalarField phiGas(fvc::interpolate(phiHsGas_));
+        hEqn += fvc::div(phiGas);
+    }
+
+    if (QrHSource_)
+    {
+        const surfaceScalarField phiQr(fvc::interpolate(Qr_)*nMagSf());
+        hEqn += fvc::div(phiQr);
+    }
+
     if (regionMesh().moving())
     {
         surfaceScalarField phihMesh
@@ -271,9 +349,6 @@ void reactingOneDim::solveEnergy()
 
     hEqn.relax();
     hEqn.solve();
-
-    Info<< "pyrolysis min/max(T) = " << min(solidThermo_.T()) << ", "
-        << max(solidThermo_.T()) << endl;
 }
 
 
@@ -366,10 +441,27 @@ reactingOneDim::reactingOneDim(const word& modelType, const fvMesh& mesh)
         dimensionedScalar("zero", dimEnergy/dimTime/dimVolume, 0.0)
     ),
 
+    Qr_
+    (
+        IOobject
+        (
+            "Qr",
+            time().timeName(),
+            regionMesh(),
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        regionMesh()
+        //dimensionedScalar("zero", dimEnergy/dimArea/dimTime, 0.0),
+        //zeroGradientFvPatchVectorField::typeName
+    ),
+
     lostSolidMass_(dimensionedScalar("zero", dimMass, 0.0)),
     addedGasMass_(dimensionedScalar("zero", dimMass, 0.0)),
     totalGasMassFlux_(0.0),
-    totalHeatRR_(dimensionedScalar("zero", dimEnergy/dimTime, 0.0))
+    totalHeatRR_(dimensionedScalar("zero", dimEnergy/dimTime, 0.0)),
+    gasHSource_(false),
+    QrHSource_(false)
 {
     if (active_)
     {
@@ -449,10 +541,25 @@ reactingOneDim::reactingOneDim
         dimensionedScalar("zero", dimEnergy/dimTime/dimVolume, 0.0)
     ),
 
+    Qr_
+    (
+        IOobject
+        (
+            "Qr",
+            time().timeName(),
+            regionMesh(),
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        regionMesh()
+    ),
+
     lostSolidMass_(dimensionedScalar("zero", dimMass, 0.0)),
     addedGasMass_(dimensionedScalar("zero", dimMass, 0.0)),
     totalGasMassFlux_(0.0),
-    totalHeatRR_(dimensionedScalar("zero", dimEnergy/dimTime, 0.0))
+    totalHeatRR_(dimensionedScalar("zero", dimEnergy/dimTime, 0.0)),
+    gasHSource_(false),
+    QrHSource_(false)
 {
     if (active_)
     {
@@ -579,8 +686,6 @@ void reactingOneDim::evolveRegion()
         time().deltaTValue()
     );
 
-    calculateMassTransfer();
-
     solveContinuity();
 
     chemistrySh_ = solidChemistry_->Sh()();
@@ -594,9 +699,15 @@ void reactingOneDim::evolveRegion()
         solveEnergy();
     }
 
+    calculateMassTransfer();
+
     solidThermo_.correct();
 
-    rho_ = solidThermo_.rho();
+    Info<< "pyrolysis min/max(T) = "
+        << min(solidThermo_.T().internalField())
+        << ", "
+        << max(solidThermo_.T().internalField())
+        << endl;
 }
 
 

src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.H

@@ -97,10 +97,6 @@ protected:
             volScalarField& h_;
 
 
-        //- Name of the radiative flux in the primary region
-        word primaryRadFluxName_;
-
-
         // Solution parameters
 
             //- Number of non-orthogonal correctors
@@ -125,6 +121,16 @@ protected:
             volScalarField chemistrySh_;
 
 
+        // Source term fields
+
+            //- Coupled region radiative heat flux [W/m2]
+            //  Requires user to input mapping info for coupled patches
+            //volScalarField QrCoupled_;
+
+            //- In depth radiative heat flux [W/m2]
+            volScalarField Qr_;
+
+
         // Checks
 
             //- Cumulative lost mass of the condensed phase [kg]
@@ -140,6 +146,15 @@ protected:
             dimensionedScalar totalHeatRR_;
 
 
+        // Options
+
+            //- Add gas enthalpy source term
+            bool  gasHSource_;
+
+            //- Add in depth radiation source term
+            bool  QrHSource_;
+
+
     // Protected member functions
 
         //- Read control parameters from dictionary
@@ -154,6 +169,9 @@ protected:
         //- Update/move mesh based on change in mass
         void updateMesh(const scalarField& mass0);
 
+        //- Update radiative flux in pyrolysis region
+        void updateQr();
+
         //- Update enthalpy flux for pyrolysis gases
         void updatePhiGas();
 

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -159,6 +159,8 @@ void Foam::chemistryModel<CompType, ThermoType>::updateConcsInReactionI
     const label index,
     const scalar dt,
     const scalar omeg,
+    const scalar p,
+    const scalar T,
     scalarField& c
 ) const
 {
@@ -191,6 +193,8 @@ void Foam::chemistryModel<CompType, ThermoType>::updateRRInReactionI
     const scalar corr,
     const label lRef,
     const label rRef,
+    const scalar p,
+    const scalar T,
     simpleMatrix<scalar>& RR
 ) const
 {

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -182,6 +182,8 @@ public:
             const label i,
             const scalar dt,
             const scalar omega,
+            const scalar p,
+            const scalar T,
             scalarField& c
         ) const;
 
@@ -194,6 +196,8 @@ public:
             const scalar corr,
             const label lRef,
             const label rRef,
+            const scalar p,
+            const scalar T,
             simpleMatrix<scalar>& RR
         ) const;
 

src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -94,7 +94,7 @@ Foam::scalar Foam::EulerImplicit<ChemistryModel>::solve
             }
         }
 
-        this->updateRRInReactionI(i, pr, pf, corr, lRef, rRef, RR);
+        this->updateRRInReactionI(i, pr, pf, corr, lRef, rRef, p, T, RR);
     }
 
 

src/thermophysicalModels/chemistryModel/chemistrySolver/sequential/sequential.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -83,7 +83,7 @@ Foam::scalar Foam::sequential<ChemistryModel>::solve
 
         tChemInv = max(tChemInv, mag(omega));
 
-        this->updateConcsInReactionI(i, dt, omega, c);
+        this->updateConcsInReactionI(i, dt, omega, p, T, c);
     }
 
     return cTauChem_/tChemInv;

src/thermophysicalModels/radiationModels/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C

@@ -97,15 +97,7 @@ greyDiffusiveRadiationMixedFvPatchScalarField
     }
     else
     {
-        // No value given. Restart as fixedValue b.c.
+        refValue() = 0.0;
-
-        const scalarField& Tp =
-            patch().lookupPatchField<volScalarField, scalar>(TName_);
-
-        //NOTE: Assumes emissivity = 1 as the solidThermo might
-        // not be constructed yet
-        refValue() =
-            4.0*physicoChemical::sigma.value()*pow4(Tp)/pi;
         refGrad() = 0.0;
         valueFraction() = 1.0;
 

src/thermophysicalModels/radiationModels/submodels/absorptionEmissionModel/greyMeanSolidAbsorptionEmission/greyMeanSolidAbsorptionEmission.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -121,6 +121,21 @@ greyMeanSolidAbsorptionEmission
             continue;
         }
         const word& key = iter().keyword();
+        if (!mixture_.contains(key))
+        {
+            WarningIn
+            (
+                "greyMeanSolidAbsorptionEmission::"
+                "greyMeanSolidAbsorptionEmission "
+                "("
+                "   const dictionary& dict,"
+                "   const fvMesh& mesh"
+                ")"
+            )   << " specie: " << key << " is not found in the solid mixture"
+                << nl
+                << " specie is the mixture are:" << mixture_.species() << nl
+                << nl << endl;
+        }
         speciesNames_.insert(key, nFunc);
         const dictionary& dict = iter().dict();
         dict.lookup("absorptivity") >> solidData_[nFunc][absorptivity];

src/thermophysicalModels/reactionThermo/chemistryReaders/chemistryReader/makeChemistryReaders.C

@@ -80,9 +80,11 @@ makeChemistryReaderType(foamChemistryReader, icoPoly8EThermoPhysics);
 
 makeChemistryReader(hConstSolidThermoPhysics);
 makeChemistryReader(hExponentialSolidThermoPhysics);
+makeChemistryReader(hExpKappaConstSolidThermoPhysics);
 
 makeChemistryReaderType(foamChemistryReader, hConstSolidThermoPhysics);
 makeChemistryReaderType(foamChemistryReader, hExponentialSolidThermoPhysics);
+makeChemistryReaderType(foamChemistryReader, hExpKappaConstSolidThermoPhysics);
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/thermophysicalModels/solidChemistryModel/basicSolidChemistryModel/basicSolidChemistryModels.C

@@ -59,6 +59,15 @@ namespace Foam
         hExponentialSolidThermoPhysics,
         gasHThermoPhysics
     );
+
+    makeSolidChemistryModel
+    (
+        solidChemistryModel,
+        pyrolysisChemistryModel,
+        basicSolidChemistryModel,
+        hExpKappaConstSolidThermoPhysics,
+        gasHThermoPhysics
+    );
 }
 
 // ************************************************************************* //

src/thermophysicalModels/solidChemistryModel/makeSolidChemistryModel.H

@@ -56,7 +56,8 @@ namespace Foam
     defineTemplateTypeNameAndDebugWithName                                    \
     (                                                                         \
         SS##Comp##SThermo##GThermo,                                           \
-        (#SS"<"#Comp"," + SThermo::typeName() + "," + GThermo::typeName() +   \
+        (word(SS::typeName_()) + "<"#Comp"," + SThermo::typeName() + "," +    \
+        GThermo::typeName() +                                                 \
         ">").c_str(),                                                         \
         0                                                                     \
     );

src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.C

@@ -274,6 +274,29 @@ Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::omega
 }
 
 
+template<class CompType, class SolidThermo, class GasThermo>
+Foam::scalar Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::
+omegaI
+(
+    const label index,
+    const scalarField& c,
+    const scalar T,
+    const scalar p,
+    scalar& pf,
+    scalar& cf,
+    label& lRef,
+    scalar& pr,
+    scalar& cr,
+    label& rRef
+) const
+{
+
+    const Reaction<SolidThermo>& R = this->reactions_[index];
+    scalar w = omega(R, c, T, p, pf, cf, lRef, pr, cr, rRef);
+    return(w);
+}
+
+
 template<class CompType, class SolidThermo, class GasThermo>
 void Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::
 derivatives
@@ -489,6 +512,76 @@ calculate()
 }
 
 
+template<class CompType, class SolidThermo, class GasThermo>
+void Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::
+updateConcsInReactionI
+(
+    const label index,
+    const scalar dt,
+    const scalar omeg,
+    const scalar p,
+    const scalar T,
+    scalarField& c
+) const
+{
+     // update species
+    const Reaction<SolidThermo>& R = this->reactions_[index];
+    scalar rhoL = 0.0;
+    forAll(R.lhs(), s)
+    {
+        label si = R.lhs()[s].index;
+        rhoL = this->solidThermo_[si].rho(p, T);
+        c[si] -= dt*omeg;
+        c[si] = max(0.0, c[si]);
+    }
+
+    forAll(R.rhs(), s)
+    {
+        label si = R.rhs()[s].index;
+        const scalar rhoR = this->solidThermo_[si].rho(p, T);
+        const scalar sr = rhoR/rhoL;
+        c[si] += dt*sr*omeg;
+        c[si] = max(0.0, c[si]);
+    }
+}
+
+
+template<class CompType, class SolidThermo, class GasThermo>
+void Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::
+updateRRInReactionI
+(
+    const label index,
+    const scalar pr,
+    const scalar pf,
+    const scalar corr,
+    const label lRef,
+    const label rRef,
+    const scalar p,
+    const scalar T,
+    simpleMatrix<scalar>& RR
+) const
+{
+    const Reaction<SolidThermo>& R = this->reactions_[index];
+    scalar rhoL = 0.0;
+    forAll(R.lhs(), s)
+    {
+        label si = R.lhs()[s].index;
+        rhoL = this->solidThermo_[si].rho(p, T);
+        RR[si][rRef] -= pr*corr;
+        RR[si][lRef] += pf*corr;
+    }
+
+    forAll(R.rhs(), s)
+    {
+        label si = R.rhs()[s].index;
+        const scalar rhoR = this->solidThermo_[si].rho(p, T);
+        const scalar sr = rhoR/rhoL;
+        RR[si][lRef] -= sr*pf*corr;
+        RR[si][rRef] += sr*pr*corr;
+    }
+}
+
+
 template<class CompType, class SolidThermo, class GasThermo>
 Foam::scalar
 Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::solve

src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.H

@@ -153,9 +153,52 @@ public:
             label& rRef
         ) const;
 
+        //- Return the reaction rate for iReaction and the reference
+        //  species and charateristic times
+        virtual scalar omegaI
+        (
+            label iReaction,
+            const scalarField& c,
+            const scalar T,
+            const scalar p,
+            scalar& pf,
+            scalar& cf,
+            label& lRef,
+            scalar& pr,
+            scalar& cr,
+            label& rRef
+        ) const;
+
         //- Calculates the reaction rates
         virtual void calculate();
 
+        //- Update concentrations in reaction i given dt and reaction rate
+        // omega used by sequential solver
+        virtual void updateConcsInReactionI
+        (
+            const label i,
+            const scalar dt,
+            const scalar omega,
+            const scalar p,
+            const scalar T,
+            scalarField& c
+        ) const;
+
+
+        //- Update matrix RR for reaction i. Used by EulerImplicit
+        virtual void updateRRInReactionI
+        (
+            const label i,
+            const scalar pr,
+            const scalar pf,
+            const scalar corr,
+            const label lRef,
+            const label rRef,
+            const scalar p,
+            const scalar T,
+            simpleMatrix<scalar>& RR
+        ) const;
+
 
         // Chemistry model functions
 

src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H

@@ -43,6 +43,7 @@ SourceFiles
 #include "ODE.H"
 #include "volFieldsFwd.H"
 #include "DimensionedField.H"
+#include "simpleMatrix.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -150,6 +151,49 @@ public:
             label& rRef
         ) const = 0;
 
+        //- Return the reaction rate for iReaction and the reference
+        //  species and charateristic times
+        virtual scalar omegaI
+        (
+            label iReaction,
+            const scalarField& c,
+            const scalar T,
+            const scalar p,
+            scalar& pf,
+            scalar& cf,
+            label& lRef,
+            scalar& pr,
+            scalar& cr,
+            label& rRef
+        ) const = 0;
+
+        //- Update concentrations in reaction i given dt and reaction rate
+        // omega used by sequential solver
+        virtual void updateConcsInReactionI
+        (
+            const label i,
+            const scalar dt,
+            const scalar omega,
+            const scalar p,
+            const scalar T,
+            scalarField& c
+        ) const = 0;
+
+
+        //- Update matrix RR for reaction i. Used by EulerImplicit
+        virtual void updateRRInReactionI
+        (
+            const label i,
+            const scalar pr,
+            const scalar pf,
+            const scalar corr,
+            const label lRef,
+            const label rRef,
+            const scalar p,
+            const scalar T,
+            simpleMatrix<scalar>& RR
+        ) const = 0;
+
         //- Calculates the reaction rates
         virtual void calculate() = 0;
 

src/thermophysicalModels/solidChemistryModel/solidChemistrySolver/makeSolidChemistrySolverType.H

@@ -32,6 +32,11 @@ Description
 
 #include "addToRunTimeSelectionTable.H"
 
+#include "noChemistrySolver.H"
+#include "EulerImplicit.H"
+#include "ode.H"
+#include "sequential.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -47,7 +52,8 @@ namespace Foam
     defineTemplateTypeNameAndDebugWithName                                    \
     (                                                                         \
         SS##Schem##Comp##SThermo##GThermo,                                    \
-        (#SS"<" + word(Schem::typeName_()) + "<"#Comp"," + SThermo::typeName()\
+        (#SS"<" + word(Schem::typeName_())                                    \
+      + "<"#Comp"," + SThermo::typeName()                                     \
       + "," + GThermo::typeName() + ">>").c_str(),                            \
         0                                                                     \
     );                                                                        \
@@ -60,6 +66,45 @@ namespace Foam
     );
 
 
+#define makeSolidChemistrySolverTypes(SolidChem, Comp, SThermo, GThermo)      \
+                                                                              \
+    makeSolidChemistrySolverType                                              \
+    (                                                                         \
+        noChemistrySolver,                                                    \
+        SolidChem,                                                            \
+        Comp,                                                                 \
+        SThermo,                                                              \
+        GThermo                                                               \
+    );                                                                        \
+                                                                              \
+    makeSolidChemistrySolverType                                              \
+    (                                                                         \
+        EulerImplicit,                                                        \
+        SolidChem,                                                            \
+        Comp,                                                                 \
+        SThermo,                                                              \
+        GThermo                                                               \
+    );                                                                        \
+                                                                              \
+    makeSolidChemistrySolverType                                              \
+    (                                                                         \
+        ode,                                                                  \
+        SolidChem,                                                            \
+        Comp,                                                                 \
+        SThermo,                                                              \
+        GThermo                                                               \
+    );                                                                        \
+                                                                              \
+    makeSolidChemistrySolverType                                              \
+    (                                                                         \
+        sequential,                                                           \
+        SolidChem,                                                            \
+        Comp,                                                                 \
+        SThermo,                                                              \
+        GThermo                                                               \
+    );
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/thermophysicalModels/solidChemistryModel/solidChemistrySolver/makeSolidChemistrySolvers.C

@@ -30,30 +30,34 @@ License
 #include "pyrolysisChemistryModel.H"
 #include "basicSolidChemistryModel.H"
 
-#include "ode.H"
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    makeSolidChemistrySolverType
+    makeSolidChemistrySolverTypes
     (
-        ode,
         pyrolysisChemistryModel,
         basicSolidChemistryModel,
         hConstSolidThermoPhysics,
         gasHThermoPhysics
     )
 
-    makeSolidChemistrySolverType
+    makeSolidChemistrySolverTypes
     (
-        ode,
         pyrolysisChemistryModel,
         basicSolidChemistryModel,
         hExponentialSolidThermoPhysics,
         gasHThermoPhysics
     )
 
+    makeSolidChemistrySolverTypes
+    (
+        pyrolysisChemistryModel,
+        basicSolidChemistryModel,
+        hExpKappaConstSolidThermoPhysics,
+        gasHThermoPhysics
+    )
+
 }
 
 

src/thermophysicalModels/solidSpecie/include/solidThermoPhysicsTypes.H

@@ -76,7 +76,6 @@ namespace Foam
             >
         > hExponentialSolidThermoPhysics;
 
-
     typedef
         polynomialSolidTransport
         <
@@ -91,6 +90,19 @@ namespace Foam
             >,
             8
         > hTransportThermoPoly8SolidThermoPhysics;
+
+     typedef
+        constIsoSolidTransport
+        <
+            species::thermo
+            <
+                hExponentialThermo
+                <
+                    rhoConst<specie>
+                >,
+                sensibleEnthalpy
+            >
+        > hExpKappaConstSolidThermoPhysics;
 }
 
 

src/thermophysicalModels/solidSpecie/reaction/reactions/makeSolidReactions.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -42,6 +42,12 @@ makeSolidIRReactions
     solidArrheniusReactionRate
 )
 
+makeSolidIRReactions
+(
+    hExpKappaConstSolidThermoPhysics,
+    solidArrheniusReactionRate
+)
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/thermophysicalModels/solidThermo/solidReactionThermo/solidReactionThermos.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -62,6 +62,19 @@ makeReactingSolidThermo
 );
 
 
+makeReactingSolidThermo
+(
+    solidReactionThermo,
+    heSolidThermo,
+    reactingMixture,
+    constIsoSolidTransport,
+    sensibleEnthalpy,
+    hExponentialThermo,
+    rhoConst,
+    specie
+);
+
+
 makeReactingSolidThermo
 (
     solidThermo,

src/thermophysicalModels/specie/thermo/hExponential/hExponentialThermo.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -28,6 +28,21 @@ License
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
+template<class equationOfState>
+Foam::hExponentialThermo<equationOfState>::hExponentialThermo(Istream& is)
+:
+    equationOfState(is),
+    n0_(readScalar(is)),
+    Tref_(readScalar(is)),
+    Hf_(readScalar(is))
+{
+    is.check("hExponentialThermo::hExponentialThermo(Istream& is)");
+
+    c0_ *= this->W();
+    Hf_ *= this->W();
+}
+
+
 template<class equationOfState>
 Foam::hExponentialThermo<equationOfState>::hExponentialThermo
 (
@@ -39,7 +54,10 @@ Foam::hExponentialThermo<equationOfState>::hExponentialThermo
     n0_(readScalar(dict.subDict("thermodynamics").lookup("n0"))),
     Tref_(readScalar(dict.subDict("thermodynamics").lookup("Tref"))),
     Hf_(readScalar(dict.subDict("thermodynamics").lookup("Hf")))
-{}
+{
+    c0_ *= this->W();
+    Hf_ *= this->W();
+}
 
 
 // * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //

src/thermophysicalModels/specie/thermo/hExponential/hExponentialThermo.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -48,6 +48,20 @@ namespace Foam
 
 template<class EquationOfState> class hExponentialThermo;
 
+template<class EquationOfState>
+inline hExponentialThermo<EquationOfState> operator+
+(
+    const hExponentialThermo<EquationOfState>&,
+    const hExponentialThermo<EquationOfState>&
+);
+
+template<class EquationOfState>
+inline hExponentialThermo<EquationOfState> operator-
+(
+    const hExponentialThermo<EquationOfState>&,
+    const hExponentialThermo<EquationOfState>&
+);
+
 template<class EquationOfState>
 inline hExponentialThermo<EquationOfState> operator*
 (
@@ -56,6 +70,14 @@ inline hExponentialThermo<EquationOfState> operator*
 );
 
 
+template<class EquationOfState>
+inline hExponentialThermo<EquationOfState> operator==
+(
+    const hExponentialThermo<EquationOfState>&,
+    const hExponentialThermo<EquationOfState>&
+);
+
+
 template<class EquationOfState>
 Ostream& operator<<
 (
@@ -90,11 +112,6 @@ class hExponentialThermo
         //- Integrate Cp expression
         inline scalar integrateCp(const scalar T) const;
 
-
-public:
-
-    // Constructors
-
         //- Construct from components
         inline hExponentialThermo
         (
@@ -105,6 +122,13 @@ public:
             const scalar Hf
         );
 
+public:
+
+    // Constructors
+
+        //- Construct from Istream
+        hExponentialThermo(Istream&);
+
         //- Construct from dictionary
         hExponentialThermo(const dictionary&);
 
@@ -115,6 +139,15 @@ public:
             const hExponentialThermo&
         );
 
+         //- Construct and return a clone
+        inline autoPtr<hExponentialThermo> clone() const;
+
+        //- Selector from Istream
+        inline static autoPtr<hExponentialThermo> New(Istream& is);
+
+        //- Selector from dictionary
+        inline static autoPtr<hExponentialThermo> New(const dictionary& dict);
+
 
     // Member Functions
 
@@ -148,16 +181,23 @@ public:
 
     // Member operators
 
-        inline hExponentialThermo& operator=
-        (
-            const hExponentialThermo&
-        );
         inline void operator+=(const hExponentialThermo&);
         inline void operator-=(const hExponentialThermo&);
 
 
     // Friend operators
 
+        friend hExponentialThermo operator+ <EquationOfState>
+        (
+            const hExponentialThermo&,
+            const hExponentialThermo&
+        );
+
+        friend hExponentialThermo operator- <EquationOfState>
+        (
+            const hExponentialThermo&,
+            const hExponentialThermo&
+        );
 
         friend hExponentialThermo operator* <EquationOfState>
         (
@@ -166,6 +206,13 @@ public:
         );
 
 
+        friend hExponentialThermo operator== <EquationOfState>
+        (
+            const hExponentialThermo&,
+            const hExponentialThermo&
+        );
+
+
     // Ostream Operator
 
         friend Ostream& operator<< <EquationOfState>

src/thermophysicalModels/specie/thermo/hExponential/hExponentialThermoI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -55,14 +55,11 @@ inline Foam::scalar Foam::hExponentialThermo<equationOfState>::integrateCp
 {
     return
     (
-        c0_*pow(T, n0_ + 1.0)
+        c0_*pow(T, n0_ + 1.0)/(pow(Tref_, n0_)*(n0_ + 1.0))
-       /(pow(Tref_, n0_)*(n0_ + 1.0))
     );
 }
 
 
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
 template<class equationOfState>
 inline Foam::hExponentialThermo<equationOfState>::hExponentialThermo
 (
@@ -78,6 +75,8 @@ inline Foam::hExponentialThermo<equationOfState>::hExponentialThermo
 {}
 
 
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
 template<class equationOfState>
 inline Foam::hExponentialThermo<equationOfState>::hExponentialThermo
 (
@@ -96,6 +95,39 @@ inline Foam::hExponentialThermo<equationOfState>::hExponentialThermo
 {}
 
 
+template<class equationOfState>
+inline Foam::autoPtr<Foam::hExponentialThermo<equationOfState> >
+Foam::hExponentialThermo<equationOfState>::clone() const
+{
+    return autoPtr<hExponentialThermo<equationOfState> >
+    (
+        new hExponentialThermo<equationOfState>(*this)
+    );
+}
+
+
+template<class equationOfState>
+inline Foam::autoPtr<Foam::hExponentialThermo<equationOfState> >
+Foam::hExponentialThermo<equationOfState>::New(Istream& is)
+{
+    return autoPtr<hExponentialThermo<equationOfState> >
+    (
+        new hExponentialThermo<equationOfState>(is)
+    );
+}
+
+
+template<class equationOfState>
+inline Foam::autoPtr<Foam::hExponentialThermo<equationOfState> >
+Foam::hExponentialThermo<equationOfState>::New(const dictionary& dict)
+{
+    return autoPtr<hExponentialThermo<equationOfState> >
+    (
+        new hExponentialThermo<equationOfState>(dict)
+    );
+}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class equationOfState>
@@ -114,7 +146,7 @@ inline Foam::scalar Foam::hExponentialThermo<equationOfState>::cp
     const scalar p, const scalar T
 ) const
 {
-    return c0_*pow(T/Tref_, n0_)*this->W();
+    return c0_*pow(T/Tref_, n0_);
 }
 
 
@@ -128,7 +160,7 @@ inline Foam::scalar Foam::hExponentialThermo<equationOfState>::ha
 
     return
     (
-        (integrateCp(T) + Hf_ - hOffset)*this->W()
+        (integrateCp(T) + Hf_ - hOffset)
     );
 }
 
@@ -140,14 +172,14 @@ inline Foam::scalar Foam::hExponentialThermo<equationOfState>::hs
 ) const
 {
     scalar hOffset = integrateCp(specie::Tstd);
-    return (integrateCp(T) - hOffset)*this->W();
+    return (integrateCp(T) - hOffset);
 }
 
 
 template<class equationOfState>
 inline Foam::scalar Foam::hExponentialThermo<equationOfState>::hc() const
 {
-    return Hf_*this->W();
+    return Hf_;
 }
 
 
@@ -167,25 +199,6 @@ inline Foam::scalar Foam::hExponentialThermo<equationOfState>::s
 
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
-
-template<class equationOfState>
-inline Foam::hExponentialThermo<equationOfState>&
-Foam::hExponentialThermo<equationOfState>::operator=
-(
-    const hExponentialThermo<equationOfState>& ct
-)
-{
-    equationOfState::operator=(ct);
-
-    Hf_ = ct.Hf_;
-    c0_ = ct.c0_;
-    n0_ = ct.n0_;
-    Tref_ = ct.Tref_;
-
-    return *this;
-}
-
-
 template<class equationOfState>
 inline void Foam::hExponentialThermo<equationOfState>::operator+=
 (
@@ -195,14 +208,13 @@ inline void Foam::hExponentialThermo<equationOfState>::operator+=
     scalar molr1 = this->nMoles();
 
     equationOfState::operator+=(ct);
-
     molr1 /= this->nMoles();
     scalar molr2 = ct.nMoles()/this->nMoles();
 
     Hf_ = molr1*Hf_ + molr2*ct.Hf_;
     c0_ = molr1*c0_ + molr2*ct.c0_;
-    n0_ = (molr1*n0_ + molr2*ct.n0_);
+    n0_ = molr1*n0_ + molr2*ct.n0_;
-    Tref_ = (molr1*Tref_ + molr2*ct.Tref_);
+    Tref_ = molr1*Tref_ + molr2*ct.Tref_;
 }
 
 
@@ -228,6 +240,62 @@ inline void Foam::hExponentialThermo<equationOfState>::operator-=
 
 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 
+template<class equationOfState>
+inline Foam::hExponentialThermo<equationOfState> Foam::operator+
+(
+    const hExponentialThermo<equationOfState>& ct1,
+    const hExponentialThermo<equationOfState>& ct2
+)
+{
+    equationOfState eofs
+    (
+        static_cast<const equationOfState&>(ct1)
+      + static_cast<const equationOfState&>(ct2)
+    );
+
+    return hExponentialThermo<equationOfState>
+    (
+        eofs,
+        ct1.nMoles()/eofs.nMoles()*ct1.c0_
+      + ct2.nMoles()/eofs.nMoles()*ct2.c0_,
+        ct1.nMoles()/eofs.nMoles()*ct1.n0_
+      + ct2.nMoles()/eofs.nMoles()*ct2.n0_,
+        ct1.nMoles()/eofs.nMoles()*ct1.Tref_
+      + ct2.nMoles()/eofs.nMoles()*ct2.Tref_,
+        ct1.nMoles()/eofs.nMoles()*ct1.Hf_
+      + ct2.nMoles()/eofs.nMoles()*ct2.Hf_
+    );
+}
+
+
+template<class equationOfState>
+inline Foam::hExponentialThermo<equationOfState> Foam::operator-
+(
+    const hExponentialThermo<equationOfState>& ct1,
+    const hExponentialThermo<equationOfState>& ct2
+)
+{
+    equationOfState eofs
+    (
+        static_cast<const equationOfState&>(ct1)
+      + static_cast<const equationOfState&>(ct2)
+    );
+
+    return hExponentialThermo<equationOfState>
+    (
+        eofs,
+        ct1.nMoles()/eofs.nMoles()*ct1.c0_
+      - ct2.nMoles()/eofs.nMoles()*ct2.c0_,
+        ct1.nMoles()/eofs.nMoles()*ct1.n0_
+      - ct2.nMoles()/eofs.nMoles()*ct2.n0_,
+        ct1.nMoles()/eofs.nMoles()*ct1.Tref_
+      - ct2.nMoles()/eofs.nMoles()*ct2.Tref_,
+        ct1.nMoles()/eofs.nMoles()*ct1.Hf_
+      - ct2.nMoles()/eofs.nMoles()*ct2.Hf_
+    );
+}
+
+
 template<class equationOfState>
 inline Foam::hExponentialThermo<equationOfState> Foam::operator*
 (
@@ -238,12 +306,23 @@ inline Foam::hExponentialThermo<equationOfState> Foam::operator*
     return hExponentialThermo<equationOfState>
     (
         s*static_cast<const equationOfState&>(ct),
-        ct.Hf_,
         ct.c0_,
         ct.n0_,
-        ct.Tref_
+        ct.Tref_,
+        ct.Hf_
     );
 }
 
 
+template<class equationOfState>
+inline Foam::hExponentialThermo<equationOfState> Foam::operator==
+(
+    const hExponentialThermo<equationOfState>& ct1,
+    const hExponentialThermo<equationOfState>& ct2
+)
+{
+    return ct2 - ct1;
+}
+
+
 // ************************************************************************* //

tutorials/combustion/fireFoam/les/oppositeBurningPanels/0/panelRegion/Qr

@@ -36,12 +36,12 @@ boundaryField
         type            wedge;
     }
 
-    region0_to_panelRegion_left_face
+    region0_to_panelRegion_fLeft_zone
     {
         type            mappedField;
         sampleRegion    region0;
         sampleMode      nearestPatchFace;
-        samplePatch     region0_to_panelRegion_left_face;
+        samplePatch     region0_to_panelRegion_fLeft_zone;
         offset          (0 0 0);
         fieldName       Qr;
         setAverage      no;
@@ -49,12 +49,12 @@ boundaryField
         value           uniform 0;
     }
 
-    region0_to_panelRegion_right_face
+    region0_to_panelRegion_fRight_zone
     {
         type            mappedField;
         sampleRegion    region0;
         sampleMode      nearestPatchFace;
-        samplePatch     region0_to_panelRegion_right_face;
+        samplePatch     region0_to_panelRegion_fRight_zone;
         offset          (0 0 0);
         fieldName       Qr;
         setAverage      no;

tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/panelRegion/radiationProperties

@@ -23,15 +23,15 @@ absorptionEmissionModel greyMeanSolidAbsorptionEmission;
 
 greyMeanSolidAbsorptionEmissionCoeffs
 {
-    v
+      wood
     {
-        absorptivity    0.0;  //opaque
+        absorptivity    0.17;
         emissivity      0.17;
     }
 
     char
     {
-        absorptivity    0.0;
+        absorptivity    0.85;
         emissivity      0.85;
     }
 }

tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/pyrolysisZones

@@ -27,8 +27,8 @@ FoamFile
 
         reactingOneDimCoeffs
         {
-            radFluxName     Qr;
+            gasHSource      false; //Energy source term due to pyrolysis gas
-
+            QrHSource       false; //Energy source term due in depht radiation
             minimumDelta    1e-12;
 
             reactionDeltaMin 1e-12;