Commit 92d2bc1b authored by sergio's avatar sergio
Browse files

ENH: Aligning solid thermo mixing with gas mixing and solid reaction

parent 3083282e
......@@ -7,6 +7,7 @@ EXE_INC = \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/solidThermo/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/solidChemistryModel/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/radiationModels/lnInclude \
-I$(LIB_SRC)/turbulenceModels \
-I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel/lnInclude \
......@@ -28,4 +29,5 @@ LIB_LIBS = \
-lcompressibleLESModels \
-lLESdeltas \
-lregionModels \
-lradiationModels
-lradiationModels \
-lreactionThermophysicalModels
EXE_INC = \
-I$(LIB_SRC)/regionModels/regionModel/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/solidThermo/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/solid/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude \
......
......@@ -4,6 +4,7 @@ makeType=${1:-libso}
set -x
wmake $makeType specie
wmake $makeType solidSpecie
wmake $makeType thermophysicalFunctions
./properties/Allwmake $*
......@@ -15,7 +16,7 @@ wmake $makeType barotropicCompressibilityModel
wmake $makeType thermalPorousZone
wmake $makeType SLGThermo
wmake $makeType solidSpecie
wmake $makeType solidThermo
wmake $makeType solidChemistryModel
......
......@@ -87,6 +87,27 @@ public:
return mixture_;
}
const ThermoType& cellVolMixture
(
const scalar,
const scalar,
const label
) const
{
return mixture_;
}
const ThermoType& patchFaceVolMixture
(
const scalar,
const scalar,
const label,
const label
) const
{
return mixture_;
}
//- Read dictionary
void read(const dictionary&);
};
......
......@@ -27,7 +27,7 @@ License
#include "addToRunTimeSelectionTable.H"
#include "unitConversion.H"
#include "zeroGradientFvPatchFields.H"
#include "basicSolidMixture.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
......@@ -71,7 +71,7 @@ greyMeanSolidAbsorptionEmission::X(const word specie) const
}
}
const scalarField& Yj = mixture_.Y(specie);
const label mySpecieI = mixture_.components()[specie];
const label mySpecieI = mixture_.species()[specie];
forAll(Xj, iCell)
{
Xj[iCell] = Yj[iCell]/mixture_.rho(mySpecieI, p[iCell], T[iCell]);
......@@ -93,10 +93,10 @@ greyMeanSolidAbsorptionEmission
coeffsDict_((dict.subDict(typeName + "Coeffs"))),
thermo_(mesh.lookupObject<solidThermo>("thermophysicalProperties")),
speciesNames_(0),
mixture_(dynamic_cast<const basicSolidMixture&>(thermo_)),
mixture_(dynamic_cast<const basicMultiComponentMixture&>(thermo_)),
solidData_(mixture_.Y().size())
{
if (!isA<basicSolidMixture>(thermo_))
if (!isA<basicMultiComponentMixture>(thermo_))
{
FatalErrorIn
(
......
......@@ -90,7 +90,7 @@ private:
HashTable<label> speciesNames_;
//- Basic multicomponent mixture
const basicSolidMixture& mixture_;
const basicMultiComponentMixture& mixture_;
//- List of solid species data
List<FixedList<scalar, 2> > solidData_;
......
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/specie/lnInclude
-I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/solidSpecie/lnInclude
LIB_LIBS = \
-lfiniteVolume
......@@ -25,6 +25,7 @@ License
#include "makeReactionThermo.H"
#include "thermoPhysicsTypes.H"
#include "solidThermoPhysicsTypes.H"
#include "chemistryReader.H"
#include "foamChemistryReader.H"
......@@ -41,6 +42,8 @@ makeChemistryReader(gasThermoPhysics);
makeChemistryReader(constIncompressibleGasThermoPhysics);
makeChemistryReader(incompressibleGasThermoPhysics);
makeChemistryReader(icoPoly8ThermoPhysics);
makeChemistryReader(hConstSolidThermoPhysics);
makeChemistryReader(hExponentialSolidThermoPhysics);
makeChemistryReaderType(foamChemistryReader, constGasThermoPhysics);
makeChemistryReaderType(foamChemistryReader, gasThermoPhysics);
......@@ -51,6 +54,8 @@ makeChemistryReaderType
);
makeChemistryReaderType(foamChemistryReader, incompressibleGasThermoPhysics);
makeChemistryReaderType(foamChemistryReader, icoPoly8ThermoPhysics);
makeChemistryReaderType(foamChemistryReader, hConstSolidThermoPhysics);
makeChemistryReaderType(foamChemistryReader, hExponentialSolidThermoPhysics);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......
......@@ -41,6 +41,7 @@ License
#include "powerSeriesReactionRate.H"
#include "addToRunTimeSelectionTable.H"
/* * * * * * * * * * * * * * * * * Static data * * * * * * * * * * * * * * * */
namespace Foam
......@@ -178,7 +179,8 @@ void Foam::chemkinReader::addReactionType
{
reactions_.append
(
new IrreversibleReaction<gasThermoPhysics, ReactionRateType>
new IrreversibleReaction
<Reaction, gasThermoPhysics, ReactionRateType>
(
Reaction<gasThermoPhysics>
(
......@@ -197,7 +199,8 @@ void Foam::chemkinReader::addReactionType
{
reactions_.append
(
new ReversibleReaction<gasThermoPhysics, ReactionRateType>
new ReversibleReaction
<Reaction, gasThermoPhysics, ReactionRateType>
(
Reaction<gasThermoPhysics>
(
......@@ -496,7 +499,7 @@ void Foam::chemkinReader::addReaction
reactions_.append
(
new NonEquilibriumReversibleReaction
<gasThermoPhysics, ArrheniusReactionRate>
<Reaction, gasThermoPhysics, ArrheniusReactionRate>
(
Reaction<gasThermoPhysics>
(
......@@ -549,7 +552,11 @@ void Foam::chemkinReader::addReaction
reactions_.append
(
new NonEquilibriumReversibleReaction
<gasThermoPhysics, thirdBodyArrheniusReactionRate>
<
Reaction,
gasThermoPhysics,
thirdBodyArrheniusReactionRate
>
(
Reaction<gasThermoPhysics>
(
......@@ -654,7 +661,7 @@ void Foam::chemkinReader::addReaction
reactions_.append
(
new NonEquilibriumReversibleReaction
<gasThermoPhysics, LandauTellerReactionRate>
<Reaction, gasThermoPhysics, LandauTellerReactionRate>
(
Reaction<gasThermoPhysics>
(
......
......@@ -77,7 +77,8 @@ Foam::multiComponentMixture<ThermoType>::multiComponentMixture
:
basicMultiComponentMixture(thermoDict, specieNames, mesh),
speciesData_(species_.size()),
mixture_("mixture", *thermoData[specieNames[0]])
mixture_("mixture", *thermoData[specieNames[0]]),
mixtureVol_("volMixture", *thermoData[specieNames[0]])
{
forAll(species_, i)
{
......@@ -101,7 +102,8 @@ Foam::multiComponentMixture<ThermoType>::multiComponentMixture
:
basicMultiComponentMixture(thermoDict, thermoDict.lookup("species"), mesh),
speciesData_(species_.size()),
mixture_("mixture", constructSpeciesData(thermoDict))
mixture_("mixture", constructSpeciesData(thermoDict)),
mixtureVol_("volMixture", speciesData_[0])
{
correctMassFractions();
}
......@@ -148,6 +150,78 @@ const ThermoType& Foam::multiComponentMixture<ThermoType>::patchFaceMixture
}
template<class ThermoType>
const ThermoType& Foam::multiComponentMixture<ThermoType>::cellVolMixture
(
const scalar p,
const scalar T,
const label celli
) const
{
scalar rhoInv = 0.0;
forAll(speciesData_, i)
{
rhoInv += Y_[i][celli]/speciesData_[i].rho(p, T);
}
mixtureVol_ =
Y_[0][celli]
/ speciesData_[0].rho(p, T)
/ rhoInv
/ speciesData_[0].W()
* speciesData_[0];
for (label n=1; n<Y_.size(); n++)
{
mixtureVol_ +=
Y_[n][celli]
/ speciesData_[n].rho(p, T)
/ rhoInv
/ speciesData_[n].W()
* speciesData_[n];
}
return mixtureVol_;
}
template<class ThermoType>
const ThermoType& Foam::multiComponentMixture<ThermoType>::
patchFaceVolMixture
(
const scalar p,
const scalar T,
const label patchi,
const label facei
) const
{
scalar rhoInv = 0.0;
forAll(speciesData_, i)
{
rhoInv += Y_[i].boundaryField()[patchi][facei]/speciesData_[i].rho(p, T);
}
mixtureVol_ =
Y_[0].boundaryField()[patchi][facei]
/ speciesData_[0].rho(p, T)
/ rhoInv
/ speciesData_[0].W()
* speciesData_[0];
for (label n=1; n<Y_.size(); n++)
{
mixtureVol_ +=
Y_[n].boundaryField()[patchi][facei]
/ speciesData_[n].rho(p,T)
/ rhoInv
/ speciesData_[n].W()
* speciesData_[n];
}
return mixtureVol_;
}
template<class ThermoType>
void Foam::multiComponentMixture<ThermoType>::read
(
......
......@@ -60,6 +60,10 @@ class multiComponentMixture
//- Temporary storage for the cell/face mixture thermo data
mutable ThermoType mixture_;
//- Temporary storage for the volume weighted
// cell/face mixture thermo data
mutable ThermoType mixtureVol_;
// Private Member Functions
......@@ -110,6 +114,21 @@ public:
const label facei
) const;
const ThermoType& cellVolMixture
(
const scalar p,
const scalar T,
const label celli
) const;
const ThermoType& patchFaceVolMixture
(
const scalar p,
const scalar T,
const label patchi,
const label facei
) const;
//- Return the raw specie thermodynamic data
const PtrList<ThermoType>& speciesData() const
{
......
......@@ -12,4 +12,5 @@ EXE_INC = \
LIB_LIBS = \
-lchemistryModel \
-lfiniteVolume \
-lODE
-lODE\
-lreactionThermophysicalModels
......@@ -24,7 +24,8 @@ License
\*---------------------------------------------------------------------------*/
#include "ODESolidChemistryModel.H"
#include "reactingSolidMixture.H"
#include "reactingMixture.H"
#include "solidReaction.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
......@@ -38,24 +39,20 @@ ODESolidChemistryModel
CompType(mesh),
ODE(),
Ys_(this->solidThermo().composition().Y()),
pyrolisisGases_
(
mesh.lookupObject<dictionary>
("chemistryProperties").lookup("species")
),
reactions_
(
dynamic_cast<const reactingSolidMixture<SolidThermo>& >
dynamic_cast<const reactingMixture<SolidThermo>& >
(
this->solidThermo()
)
),
pyrolisisGases_(reactions_[0].gasSpecies()),
solidThermo_
(
dynamic_cast<const reactingSolidMixture<SolidThermo>& >
dynamic_cast<const reactingMixture<SolidThermo>& >
(
this->solidThermo()
).solidData()
).speciesData()
),
gasThermo_(pyrolisisGases_.size()),
nGases_(pyrolisisGases_.size()),
......@@ -183,7 +180,7 @@ ODESolidChemistryModel
dictionary thermoDict =
mesh.lookupObject<dictionary>
(
"chemistryProperties"
"thermophysicalProperties"
).subDict(pyrolisisGases_[gasI]);
gasThermo_.set
......@@ -193,14 +190,13 @@ ODESolidChemistryModel
);
}
Info<< "ODESolidChemistryModel: Number of solids = " << nSolids_
<< " and reactions = " << nReaction_ << endl;
Info<< "solidChemistryModel: Number of solids = " << nSolids_ << endl;
Info<< "Number of gases from pyrolysis = " << nGases_ << endl;
Info<< "solidChemistryModel: Number of gases = " << nGases_ << endl;
forAll(reactions_, i)
{
Info<< indent << "Reaction " << i << nl << reactions_[i] << nl;
Info<< indent << reactions_[i] << nl;
}
}
......@@ -235,23 +231,23 @@ ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::omega
forAll(reactions_, i)
{
const solidReaction& R = reactions_[i];
const Reaction<SolidThermo>& R = reactions_[i];
scalar omegai = omega
(
R, c, T, p, pf, cf, lRef, pr, cr, rRef
);
scalar rhoL = 0.0;
forAll(R.slhs(), s)
forAll(R.lhs(), s)
{
label si = R.slhs()[s];
label si = R.lhs()[s].index;
om[si] -= omegai;
rhoL = solidThermo_[si].rho(p, T);
}
scalar sr = 0.0;
forAll(R.srhs(), s)
forAll(R.rhs(), s)
{
label si = R.srhs()[s];
label si = R.rhs()[s].index;
scalar rhoR = solidThermo_[si].rho(p, T);
sr = rhoR/rhoL;
om[si] += sr*omegai;
......@@ -263,7 +259,7 @@ ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::omega
}
forAll(R.grhs(), g)
{
label gi = R.grhs()[g];
label gi = R.grhs()[g].index;
om[gi + nSolids_] += (1.0 - sr)*omegai;
}
}
......@@ -276,7 +272,7 @@ template<class CompType, class SolidThermo, class GasThermo>
Foam::scalar
Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::omega
(
const solidReaction& R,
const Reaction<SolidThermo>& R,
const scalarField& c,
const scalar T,
const scalar p,
......@@ -299,16 +295,17 @@ Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::omega
scalar kf = R.kf(p, T, c1);
const scalar exponent = R.nReact();
//const scalar exponent = R.nReact();
const label Nl = R.slhs().size();
const label Nl = R.lhs().size();
for (label s=0; s<Nl; s++)
{
label si = R.slhs()[s];
label si = R.lhs()[s].index;
const scalar exp = R.lhs()[si].exponent;
kf *=
pow(c1[si]/Ys0_[si][cellI], exponent)
pow(c1[si]/Ys0_[si][cellI], exp)
*(Ys0_[si][cellI]);
}
......@@ -390,18 +387,18 @@ void Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::jacobian
for (label ri=0; ri<reactions_.size(); ri++)
{
const solidReaction& R = reactions_[ri];
const Reaction<SolidThermo>& R = reactions_[ri];
scalar kf0 = R.kf(p, T, c2);
forAll(R.slhs(), j)
forAll(R.lhs(), j)
{
label sj = R.slhs()[j];
label sj = R.lhs()[j].index;
scalar kf = kf0;
forAll(R.slhs(), i)
forAll(R.lhs(), i)
{
label si = R.slhs()[i];
scalar exp = R.nReact();
label si = R.lhs()[i].index;
scalar exp = R.lhs()[i].exponent;
if (i == j)
{
if (exp < 1.0)
......@@ -428,14 +425,14 @@ void Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::jacobian
}
}
forAll(R.slhs(), i)
forAll(R.lhs(), i)
{
label si = R.slhs()[i];
label si = R.lhs()[i].index;
dfdc[si][sj] -= kf;
}
forAll(R.srhs(), i)
forAll(R.rhs(), i)
{
label si = R.srhs()[i];
label si = R.rhs()[i].index;
dfdc[si][sj] += kf;
}
}
......@@ -551,7 +548,8 @@ ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::nEqns() const
template<class CompType, class SolidThermo, class GasThermo>
void Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::calculate()
void Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::
calculate()
{
if (!this->chemistry_)
{
......
......@@ -38,7 +38,7 @@ SourceFiles
#ifndef ODESolidChemistryModel_H
#define ODESolidChemistryModel_H
#include "solidReaction.H"
#include "Reaction.H"
#include "ODE.H"
#include "volFieldsFwd.H"
#include "DimensionedField.H"
......@@ -72,12 +72,12 @@ protected:
//- Reference to solid mass fractions
PtrList<volScalarField>& Ys_;
//- Reactions
const PtrList<Reaction<SolidThermo> >& reactions_;
//- List of gas species present in reaction system
speciesTable pyrolisisGases_;
//- Reactions
const PtrList<solidReaction>& reactions_;
//- Thermodynamic data of solids
const PtrList<SolidThermo>& solidThermo_;