diff --git a/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C b/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C index 43e2017b02..29325f602e 100644 --- a/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C +++ b/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 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,22 +242,28 @@ void reactingOneDim::solveContinuity() Info<< "reactingOneDim::solveContinuity()" << endl; } - if (moveMesh_) - { - const scalarField mass0 = rho_*regionMesh().V(); + const scalarField mass0 = rho_*regionMesh().V(); - fvScalarMatrix rhoEqn + fvScalarMatrix rhoEqn + ( + fvm::ddt(rho_) + == + - solidChemistry_->RRg() + ); + + if (regionMesh().moving()) + { + surfaceScalarField phiRhoMesh ( - fvm::ddt(rho_) - == - - solidChemistry_->RRg() + fvc::interpolate(rho_)*regionMesh().phi() ); - rhoEqn.solve(); - - updateMesh(mass0); - + 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; } diff --git a/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.H b/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.H index d323f952f5..4a92fd6967 100644 --- a/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.H +++ b/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(); diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C index 226e19c226..ed9cf97b20 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C +++ b/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::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::updateRRInReactionI const scalar corr, const label lRef, const label rRef, + const scalar p, + const scalar T, simpleMatrix& RR ) const { diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H index b663769d2a..7b3d1c663a 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H +++ b/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& RR ) const; diff --git a/src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.C b/src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.C index f0394a1c69..349ef77f8c 100644 --- a/src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.C +++ b/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::solve } } - this->updateRRInReactionI(i, pr, pf, corr, lRef, rRef, RR); + this->updateRRInReactionI(i, pr, pf, corr, lRef, rRef, p, T, RR); } diff --git a/src/thermophysicalModels/chemistryModel/chemistrySolver/sequential/sequential.C b/src/thermophysicalModels/chemistryModel/chemistrySolver/sequential/sequential.C index defef83418..b5e7297bd8 100644 --- a/src/thermophysicalModels/chemistryModel/chemistrySolver/sequential/sequential.C +++ b/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::solve tChemInv = max(tChemInv, mag(omega)); - this->updateConcsInReactionI(i, dt, omega, c); + this->updateConcsInReactionI(i, dt, omega, p, T, c); } return cTauChem_/tChemInv; diff --git a/src/thermophysicalModels/radiationModels/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C b/src/thermophysicalModels/radiationModels/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C index 52b3192c6d..b4d433d425 100644 --- a/src/thermophysicalModels/radiationModels/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C +++ b/src/thermophysicalModels/radiationModels/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C @@ -97,15 +97,7 @@ greyDiffusiveRadiationMixedFvPatchScalarField } else { - // No value given. Restart as fixedValue b.c. - - const scalarField& Tp = - patch().lookupPatchField(TName_); - - //NOTE: Assumes emissivity = 1 as the solidThermo might - // not be constructed yet - refValue() = - 4.0*physicoChemical::sigma.value()*pow4(Tp)/pi; + refValue() = 0.0; refGrad() = 0.0; valueFraction() = 1.0; diff --git a/src/thermophysicalModels/radiationModels/submodels/absorptionEmissionModel/greyMeanSolidAbsorptionEmission/greyMeanSolidAbsorptionEmission.C b/src/thermophysicalModels/radiationModels/submodels/absorptionEmissionModel/greyMeanSolidAbsorptionEmission/greyMeanSolidAbsorptionEmission.C index 5121436412..56c34ea71d 100644 --- a/src/thermophysicalModels/radiationModels/submodels/absorptionEmissionModel/greyMeanSolidAbsorptionEmission/greyMeanSolidAbsorptionEmission.C +++ b/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]; diff --git a/src/thermophysicalModels/reactionThermo/chemistryReaders/chemistryReader/makeChemistryReaders.C b/src/thermophysicalModels/reactionThermo/chemistryReaders/chemistryReader/makeChemistryReaders.C index b90e447c0b..7b0561ee7b 100644 --- a/src/thermophysicalModels/reactionThermo/chemistryReaders/chemistryReader/makeChemistryReaders.C +++ b/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); // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/thermophysicalModels/solidChemistryModel/basicSolidChemistryModel/basicSolidChemistryModels.C b/src/thermophysicalModels/solidChemistryModel/basicSolidChemistryModel/basicSolidChemistryModels.C index 905f6b17a2..85306aaf72 100644 --- a/src/thermophysicalModels/solidChemistryModel/basicSolidChemistryModel/basicSolidChemistryModels.C +++ b/src/thermophysicalModels/solidChemistryModel/basicSolidChemistryModel/basicSolidChemistryModels.C @@ -59,6 +59,15 @@ namespace Foam hExponentialSolidThermoPhysics, gasHThermoPhysics ); + + makeSolidChemistryModel + ( + solidChemistryModel, + pyrolysisChemistryModel, + basicSolidChemistryModel, + hExpKappaConstSolidThermoPhysics, + gasHThermoPhysics + ); } // ************************************************************************* // diff --git a/src/thermophysicalModels/solidChemistryModel/makeSolidChemistryModel.H b/src/thermophysicalModels/solidChemistryModel/makeSolidChemistryModel.H index febcc0ac69..7f24af81d6 100644 --- a/src/thermophysicalModels/solidChemistryModel/makeSolidChemistryModel.H +++ b/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 \ ); diff --git a/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.C b/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.C index e036fc7a1b..e301a1b85e 100644 --- a/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.C +++ b/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.C @@ -274,6 +274,29 @@ Foam::pyrolysisChemistryModel::omega } +template +Foam::scalar Foam::pyrolysisChemistryModel:: +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& R = this->reactions_[index]; + scalar w = omega(R, c, T, p, pf, cf, lRef, pr, cr, rRef); + return(w); +} + + template void Foam::pyrolysisChemistryModel:: derivatives @@ -489,6 +512,76 @@ calculate() } +template +void Foam::pyrolysisChemistryModel:: +updateConcsInReactionI +( + const label index, + const scalar dt, + const scalar omeg, + const scalar p, + const scalar T, + scalarField& c +) const +{ + // update species + const Reaction& 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 +void Foam::pyrolysisChemistryModel:: +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& RR +) const +{ + const Reaction& 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 Foam::scalar Foam::pyrolysisChemistryModel::solve diff --git a/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.H b/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.H index 6527089c36..2d4a11cc90 100644 --- a/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.H +++ b/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& RR + ) const; + // Chemistry model functions diff --git a/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H b/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H index 6f7dbe48ce..03076d79ce 100644 --- a/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H +++ b/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& RR + ) const = 0; + //- Calculates the reaction rates virtual void calculate() = 0; diff --git a/src/thermophysicalModels/solidChemistryModel/solidChemistrySolver/makeSolidChemistrySolverType.H b/src/thermophysicalModels/solidChemistryModel/solidChemistrySolver/makeSolidChemistrySolverType.H index b1e33933f6..12f842a054 100644 --- a/src/thermophysicalModels/solidChemistryModel/solidChemistrySolver/makeSolidChemistrySolverType.H +++ b/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 diff --git a/src/thermophysicalModels/solidChemistryModel/solidChemistrySolver/makeSolidChemistrySolvers.C b/src/thermophysicalModels/solidChemistryModel/solidChemistrySolver/makeSolidChemistrySolvers.C index 1d16cb506c..f1eab789ca 100644 --- a/src/thermophysicalModels/solidChemistryModel/solidChemistrySolver/makeSolidChemistrySolvers.C +++ b/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 + ) + } diff --git a/src/thermophysicalModels/solidSpecie/include/solidThermoPhysicsTypes.H b/src/thermophysicalModels/solidSpecie/include/solidThermoPhysicsTypes.H index b7eac0fc3f..968898479a 100644 --- a/src/thermophysicalModels/solidSpecie/include/solidThermoPhysicsTypes.H +++ b/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 + >, + sensibleEnthalpy + > + > hExpKappaConstSolidThermoPhysics; } diff --git a/src/thermophysicalModels/solidSpecie/reaction/reactions/makeSolidReactions.C b/src/thermophysicalModels/solidSpecie/reaction/reactions/makeSolidReactions.C index 8fe56961ca..3241eb9ed8 100644 --- a/src/thermophysicalModels/solidSpecie/reaction/reactions/makeSolidReactions.C +++ b/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 diff --git a/src/thermophysicalModels/solidThermo/solidReactionThermo/solidReactionThermos.C b/src/thermophysicalModels/solidThermo/solidReactionThermo/solidReactionThermos.C index e3a95e079c..f8fde90475 100644 --- a/src/thermophysicalModels/solidThermo/solidReactionThermo/solidReactionThermos.C +++ b/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, diff --git a/src/thermophysicalModels/specie/thermo/hExponential/hExponentialThermo.C b/src/thermophysicalModels/specie/thermo/hExponential/hExponentialThermo.C index 552ac7bb01..562ae0d63d 100644 --- a/src/thermophysicalModels/specie/thermo/hExponential/hExponentialThermo.C +++ b/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 +Foam::hExponentialThermo::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 Foam::hExponentialThermo::hExponentialThermo ( @@ -39,7 +54,10 @@ Foam::hExponentialThermo::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 * * * * * * * * * * * * * // diff --git a/src/thermophysicalModels/specie/thermo/hExponential/hExponentialThermo.H b/src/thermophysicalModels/specie/thermo/hExponential/hExponentialThermo.H index a54dfafcda..05ee681ea6 100644 --- a/src/thermophysicalModels/specie/thermo/hExponential/hExponentialThermo.H +++ b/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 hExponentialThermo; +template +inline hExponentialThermo operator+ +( + const hExponentialThermo&, + const hExponentialThermo& +); + +template +inline hExponentialThermo operator- +( + const hExponentialThermo&, + const hExponentialThermo& +); + template inline hExponentialThermo operator* ( @@ -56,6 +70,14 @@ inline hExponentialThermo operator* ); +template +inline hExponentialThermo operator== +( + const hExponentialThermo&, + const hExponentialThermo& +); + + template 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 clone() const; + + //- Selector from Istream + inline static autoPtr New(Istream& is); + + //- Selector from dictionary + inline static autoPtr 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+ + ( + const hExponentialThermo&, + const hExponentialThermo& + ); + + friend hExponentialThermo operator- + ( + const hExponentialThermo&, + const hExponentialThermo& + ); friend hExponentialThermo operator* ( @@ -166,6 +206,13 @@ public: ); + friend hExponentialThermo operator== + ( + const hExponentialThermo&, + const hExponentialThermo& + ); + + // Ostream Operator friend Ostream& operator<< diff --git a/src/thermophysicalModels/specie/thermo/hExponential/hExponentialThermoI.H b/src/thermophysicalModels/specie/thermo/hExponential/hExponentialThermoI.H index ac4a671959..d0130eb512 100644 --- a/src/thermophysicalModels/specie/thermo/hExponential/hExponentialThermoI.H +++ b/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::integrateCp { return ( - c0_*pow(T, n0_ + 1.0) - /(pow(Tref_, n0_)*(n0_ + 1.0)) + c0_*pow(T, n0_ + 1.0)/(pow(Tref_, n0_)*(n0_ + 1.0)) ); } -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - template inline Foam::hExponentialThermo::hExponentialThermo ( @@ -78,6 +75,8 @@ inline Foam::hExponentialThermo::hExponentialThermo {} +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + template inline Foam::hExponentialThermo::hExponentialThermo ( @@ -96,6 +95,39 @@ inline Foam::hExponentialThermo::hExponentialThermo {} +template +inline Foam::autoPtr > +Foam::hExponentialThermo::clone() const +{ + return autoPtr > + ( + new hExponentialThermo(*this) + ); +} + + +template +inline Foam::autoPtr > +Foam::hExponentialThermo::New(Istream& is) +{ + return autoPtr > + ( + new hExponentialThermo(is) + ); +} + + +template +inline Foam::autoPtr > +Foam::hExponentialThermo::New(const dictionary& dict) +{ + return autoPtr > + ( + new hExponentialThermo(dict) + ); +} + + // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // template @@ -114,7 +146,7 @@ inline Foam::scalar Foam::hExponentialThermo::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::ha return ( - (integrateCp(T) + Hf_ - hOffset)*this->W() + (integrateCp(T) + Hf_ - hOffset) ); } @@ -140,14 +172,14 @@ inline Foam::scalar Foam::hExponentialThermo::hs ) const { scalar hOffset = integrateCp(specie::Tstd); - return (integrateCp(T) - hOffset)*this->W(); + return (integrateCp(T) - hOffset); } template inline Foam::scalar Foam::hExponentialThermo::hc() const { - return Hf_*this->W(); + return Hf_; } @@ -167,25 +199,6 @@ inline Foam::scalar Foam::hExponentialThermo::s // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // - -template -inline Foam::hExponentialThermo& -Foam::hExponentialThermo::operator= -( - const hExponentialThermo& ct -) -{ - equationOfState::operator=(ct); - - Hf_ = ct.Hf_; - c0_ = ct.c0_; - n0_ = ct.n0_; - Tref_ = ct.Tref_; - - return *this; -} - - template inline void Foam::hExponentialThermo::operator+= ( @@ -195,14 +208,13 @@ inline void Foam::hExponentialThermo::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_); - Tref_ = (molr1*Tref_ + molr2*ct.Tref_); + n0_ = molr1*n0_ + molr2*ct.n0_; + Tref_ = molr1*Tref_ + molr2*ct.Tref_; } @@ -228,6 +240,62 @@ inline void Foam::hExponentialThermo::operator-= // * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * // +template +inline Foam::hExponentialThermo Foam::operator+ +( + const hExponentialThermo& ct1, + const hExponentialThermo& ct2 +) +{ + equationOfState eofs + ( + static_cast(ct1) + + static_cast(ct2) + ); + + return hExponentialThermo + ( + 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 +inline Foam::hExponentialThermo Foam::operator- +( + const hExponentialThermo& ct1, + const hExponentialThermo& ct2 +) +{ + equationOfState eofs + ( + static_cast(ct1) + + static_cast(ct2) + ); + + return hExponentialThermo + ( + 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 inline Foam::hExponentialThermo Foam::operator* ( @@ -238,12 +306,23 @@ inline Foam::hExponentialThermo Foam::operator* return hExponentialThermo ( s*static_cast(ct), - ct.Hf_, ct.c0_, ct.n0_, - ct.Tref_ + ct.Tref_, + ct.Hf_ ); } +template +inline Foam::hExponentialThermo Foam::operator== +( + const hExponentialThermo& ct1, + const hExponentialThermo& ct2 +) +{ + return ct2 - ct1; +} + + // ************************************************************************* // diff --git a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/0/panelRegion/Qr b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/0/panelRegion/Qr index 6cbd558c50..6b451a9d53 100644 --- a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/0/panelRegion/Qr +++ b/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; diff --git a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/panelRegion/radiationProperties b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/panelRegion/radiationProperties index bf592a115d..18b20332bb 100644 --- a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/panelRegion/radiationProperties +++ b/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; } } diff --git a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/pyrolysisZones b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/pyrolysisZones index cef16b64a5..1c0c248345 100644 --- a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/pyrolysisZones +++ b/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;