research/OpenFOAM-4.x
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

remove CMC methods and pure virtual methods added to chemistryModels

Browse source
This commit is contained in:
ignis 2018-01-23 21:35:04 +09:00
parent 9626afd9f8
commit a8044e1f53
3 changed files with 8 additions and 408 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

34
src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H
View file

@ -178,45 +178,23 @@ public:
//- Return the mass fraction of the *specieName
//- 28.Nov.2017 Jinwoo Park
virtual tmp<volScalarField> Yspecie(const word *specieName) const = 0;
virtual tmp<volScalarField> Yspecie(const word *specieName) {}
//- Return the mass fraction of the *specieName
//- 28.Nov.2017 Jinwoo Park
virtual tmp<volScalarField> MoleFracSpecie(const word *specieName) const = 0;
virtual tmp<volScalarField> MoleFracSpecie(const word *specieName) {}
//- Return the epsilon
virtual tmp<volScalarField> epsilon() const = 0;
virtual tmp<volScalarField> epsilon() {}
//- Return the k
virtual tmp<volScalarField> k() const = 0;
virtual tmp<volScalarField> k() {}
//- Return the U field
virtual tmp<volVectorField> U() const = 0;
virtual tmp<volVectorField> U() {}
//- Return the turbulent kinematic diffusivity of chemical species (Le=1 is assumed)
virtual tmp<volScalarField> nut() const = 0;
virtual scalar solveCMCchem
(
scalar deltaT, scalar& rho, scalar& T,
scalar& p, scalarField& Qi, scalar& Qh,
scalarField& RRCMC, scalar& ChemDeltaT
) = 0; //CMC
virtual scalar calculateTCMC(scalar& Qh, scalarField& Qi, scalar& Told, scalar& rho, scalar& p) = 0; //CMC
virtual scalar calculateRHOCMC(scalarField& Qi, scalar& T, scalar& p) = 0; //CMC
virtual scalar calculateHCMC(scalarField& Qi, scalar& T, scalar& rho, scalar& p) = 0; //CMC
virtual scalar calculateShCMC(scalarField& RRCMC, label i) = 0; //CMC
virtual void correction
(
scalarField& b, scalarField& Ta, scalarField& Wa, scalarField& Wb,
scalarField& Gab, scalarField& Gat, scalarField& Gbt,
scalarField& Gaa, scalarField& Gbb, scalar& Gtt, scalar& rhorho
) = 0;
//- Calculate the reaction rate (EDM or finite-rate/EDM)
//- 11.Apr.2017 Karam Han (updated by Jinwoo Park)
virtual void calculateEDM(bool finiteRate, const scalar A, const scalar B) = 0;
virtual tmp<volScalarField> nut() {}
};

353
src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C
View file

@ -27,7 +27,6 @@ License
#include "reactingMixture.H"
#include "UniformField.H"
#include "extrapolatedCalculatedFvPatchFields.H"
#include "wallFvPatch.H" //jinwoo, implemented according to the of231 version of simpleCoalCombustionfoam from karam
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
@ -50,12 +49,11 @@ Foam::chemistryModel<CompType, ThermoType>::chemistryModel
dynamic_cast<const reactingMixture<ThermoType>&>
(this->thermo()).speciesData()
),
nSpecie_(Y_.size()),
nReaction_(reactions_.size()),
Treact_(CompType::template lookupOrDefault<scalar>("Treact", 0.0)),
RR_(nSpecie_),
edmRR_(reactions_.size()),
kineticRR_(reactions_.size())
RR_(nSpecie_)
{
// create the fields for the chemistry sources
forAll(RR_, fieldi)
@ -79,45 +77,6 @@ Foam::chemistryModel<CompType, ThermoType>::chemistryModel
);
}
forAll(edmRR_, reactioni)
{
edmRR_.set
(
reactioni,
new volScalarField
(
IOobject
(
"edmRR." + reactions_[reactioni].name(),
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
)
);
kineticRR_.set
(
reactioni,
new volScalarField
(
IOobject
(
"kineticRR." + reactions_[reactioni].name(),
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
)
);
}
Info<< "chemistryModel: Number of species = " << nSpecie_
<< " and reactions = " << nReaction_ << endl;
}
@ -1161,313 +1120,5 @@ void Foam::chemistryModel<CompType, ThermoType>::solve
NotImplemented;
}
template<class CompType, class ThermoType>
Foam::scalar Foam::chemistryModel<CompType, ThermoType>::solveCMCchem
(
scalar deltaT, scalar& rho, scalar& T, scalar& p, scalarField& Qi, scalar& Qh, scalarField& RRCMC, scalar& ChemDeltaT
)
{
scalar deltaTMin = GREAT;
//if(!this->chemistry_)
//{
// return deltaTMin;
//}
scalar rhoi = rho;
scalar Ti = T;
scalar pi = p;
scalarField c(nSpecie_);
scalarField c0(nSpecie_);
for(label i=0 ; i<nSpecie_ ; i++)
{
c[i] = rhoi*Qi[i]/specieThermo_[i].W();
c0[i] = c[i];
}
scalar timeLeft = deltaT;
while(timeLeft > SMALL)
{
scalar dt = timeLeft;
this->solve(c, Ti, pi, dt, ChemDeltaT);
timeLeft -= dt;
}
deltaTMin = min(ChemDeltaT, deltaTMin);
for(label i=0 ; i<nSpecie_; i++)
{
RRCMC[i] = (c[i]-c0[i])*specieThermo_[i].W()/deltaT;
}
return deltaTMin;
}
template<class CompType, class ThermoType>
Foam::scalar Foam::chemistryModel<CompType, ThermoType>::calculateTCMC
(
scalar& Qh, scalarField& Qi, scalar& Told, scalar& rho, scalar& p
)
{
scalar rhoi = rho;
scalar Ti = Told;
scalar pi = p;
scalarField c(nSpecie_);
for (label i=0; i<nSpecie_; i++)
{
c[i] = rhoi*Qi[i]/specieThermo_[i].W();
}
scalar cTot = 0.0;
// update the temperature
cTot = sum(c);
ThermoType mixture(0.0*specieThermo_[0]);
for (label i=0; i<nSpecie_; i++)
{
mixture += (c[i]/cTot)*specieThermo_[i];
}
scalar TCMC = mixture.THa(Qh, pi, Ti);
//scalar TCMCnew = mixture.THs(Qh, Ti);
return TCMC;//
}
template<class CompType, class ThermoType>
Foam::scalar Foam::chemistryModel<CompType, ThermoType>::calculateRHOCMC
(
scalarField& Qi, scalar& T, scalar& p
)
{
scalar Ti = T; //conditional temperature [K]
scalar pi = p; //pressure [Pa]
scalar MeanMW(0); //mean molecular weight [kg/kmol]
scalar InvMeanMW(0); //inverse of MeanMW [kmol/kg]
scalar RHOCMC(0); //conditional density [kg/kmol]
for (label i=0; i<nSpecie_; i++) //get mean molecular weight respect to Qi
{
InvMeanMW += (Qi[i]/specieThermo_[i].W());
}
MeanMW = 1.0/InvMeanMW ;
RHOCMC = pi/(8314.4621/MeanMW*Ti); //universial gas constant 8314.51[J/kmol K] = [kg.m2/kmol.K.sec2]
return RHOCMC;//
}
template<class CompType, class ThermoType>
Foam::scalar Foam::chemistryModel<CompType, ThermoType>::calculateHCMC
(
scalarField& Qi, scalar& T, scalar& rho, scalar& p
)
{
scalar rhoi = rho; //conditional density [kg/kmol]
scalar pi = p;
scalar Ti = T; //conditional temperature [K]
scalar HCMC(0); //conditional total enthalpy [J/kg]
scalarField c(nSpecie_);
for (label i=0; i<nSpecie_; i++)
{
c[i] = rhoi*Qi[i]/specieThermo_[i].W();
}
scalar cTot = 0.0;
// update the temperature
cTot = sum(c);
ThermoType mixture(0.0*specieThermo_[0]);
for (label i=0; i<nSpecie_; i++)
{
mixture += (c[i]/cTot)*specieThermo_[i];
}
HCMC = mixture.Ha(pi, Ti);
return HCMC;
}
template<class CompType, class ThermoType>
Foam::scalar Foam::chemistryModel<CompType, ThermoType>::calculateShCMC
(
scalarField& RRCMC, label i
)
{
scalar tSh(0);
if (this->chemistry_)
{
scalar hi = specieThermo_[i].Hc();
//Info<<"hi = "<<hi<<endl;
//Info<<"RRCMC[i] = "<<RRCMC[i]<<endl;
tSh = hi*RRCMC[i];
//Info<<"tSh = "<<tSh<<endl;
}
return tSh;
}
template<class CompType, class ThermoType>
void Foam::chemistryModel<CompType, ThermoType>::correction
(
scalarField& b, scalarField& Ta, scalarField& Wa, scalarField& Wb,
scalarField& Gab, scalarField& Gat, scalarField& Gbt,
scalarField& Gaa, scalarField& Gbb, scalar& Gtt, scalar& rhorho
)
{
for(label i=0 ; i<3 ; i++)
{
b_[i] = b[i];
Ta_[i] = Ta[i];
Wa_[i] = Wa[i];
Wb_[i] = Wb[i];
Gab_[i] = Gab[i];
Gat_[i] = Gat[i];
Gbt_[i] = Gbt[i];
Gaa_[i] = Gaa[i];
Gbb_[i] = Gbb[i];
}
Gtt_ = Gtt;
rhorho_ = rhorho;
}
/*
template<class CompType, class ThermoType>
void Foam::chemistryModel<CompType, ThermoType>::specieMoleFrac
(
char specieName
)
{
}
*/
template<class CompType, class ThermoType>
void Foam::chemistryModel<CompType, ThermoType>::calculateEDM
(
bool finiteRate, const scalar A, const scalar B
)
{
if (!this->chemistry_)
{
Info<<"Reaction rates are not calculated!!"<<endl;
Info<<"Please check constant/chemistryProperties"<<endl;
return;
}
const volScalarField& rho = this->thermo().rho();
const volScalarField& k = this->db().objectRegistry::lookupObject<volScalarField>("k");
const volScalarField& epsilon = this->db().objectRegistry::lookupObject<volScalarField>("epsilon");
const scalarField& T = this->thermo().T();
const scalarField& p = this->thermo().p();
scalar EDMrate(0.0);
forAll(rho, celli)
{
scalarField RRedm(nSpecie_, 0.0);
scalarField RRkinetic(nSpecie_, 0.0);
forAll(reactions_, reactioni)
{
const Reaction<ThermoType>& R = reactions_[reactioni];
scalar Ymin(1.0);
forAll(R.lhs(), lhs_speciei)
{
const label speciei = R.lhs()[lhs_speciei].index;
const scalar stoiCoeff = R.lhs()[lhs_speciei].stoichCoeff;
if(stoiCoeff != 0)
{
Ymin = min(Ymin, Y_[speciei][celli]/stoiCoeff/specieThermo_[speciei].W());
}
}
scalar YPmin(0.0);
scalar Num(0.0);
scalar Den(0.0);
forAll(R.rhs(), rhs_speciei)
{
const label speciei = R.rhs()[rhs_speciei].index;
const scalar stoiCoeff = R.rhs()[rhs_speciei].stoichCoeff;
if(stoiCoeff != 0)
{
Num += Y_[speciei][celli];
Den += stoiCoeff*specieThermo_[speciei].W();
}
}
YPmin = Num / Den;
Ymin = min(Ymin, B*YPmin);
EDMrate = A*rho[celli]*epsilon[celli]/k[celli]*Ymin;
scalar omegai = 0.0;
if(finiteRate == true)
{
scalar pf, cf, pr, cr;
label lRef, rRef;
scalarField c(nSpecie_, 0.0);
for(label i=0 ; i<nSpecie_ ; i++)
{
const scalar Yi = Y_[i][celli];
c[i] = rho[celli]*Yi/specieThermo_[i].W();
}
omegai = omega
(
R, c, T[celli], p[celli], pf, cf, lRef, pr, cr, rRef
);
}
forAll(R.lhs(), lhs_speciei)
{
const label speciei = R.lhs()[lhs_speciei].index;
const scalar stoiCoeff = R.lhs()[lhs_speciei].stoichCoeff;
RRedm[speciei] -= stoiCoeff*specieThermo_[speciei].W()*EDMrate;
RRkinetic[speciei] -= stoiCoeff*specieThermo_[speciei].W()*omegai;
if(lhs_speciei == 0)
{
//for post-processing
edmRR_[reactioni][celli] = stoiCoeff*specieThermo_[speciei].W()*EDMrate;
kineticRR_[reactioni][celli] = stoiCoeff*specieThermo_[speciei].W()*omegai;
}
}
forAll(R.rhs(), rhs_speciei)
{
const label speciei = R.rhs()[rhs_speciei].index;
const scalar stoiCoeff = R.rhs()[rhs_speciei].stoichCoeff;
RRedm[speciei] += stoiCoeff*specieThermo_[speciei].W()*EDMrate;
RRkinetic[speciei] += stoiCoeff*specieThermo_[speciei].W()*omegai;
}
}
for (label i=0; i<nSpecie_; i++)
{
if(finiteRate == true)
{
if(RRedm[i] == 0.0 || RRkinetic[i] == 0.0)
{
RR_[i][celli] = 0.0;
}
else
{
RR_[i][celli] = (RRedm[i]*RRkinetic[i])/(RRedm[i]+RRkinetic[i]);
}
}
else
{
RR_[i][celli] = RRedm[i];
}
}
}
}
// ************************************************************************* //

29
src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H
View file

@ -103,13 +103,6 @@ protected:
//- List of reaction rate per specie [kg/m3/s]
PtrList<DimensionedField<scalar, volMesh>> RR_;
scalarField b_, Ta_, Wa_, Wb_, Gab_, Gat_, Gbt_, Gaa_, Gbb_;
scalar Gtt_, rhorho_;
PtrList<DimensionedField<scalar, volMesh>> edmRR_;
PtrList<DimensionedField<scalar, volMesh>> kineticRR_;
// Protected Member Functions
@ -229,24 +222,6 @@ public:
// and return the characteristic time
virtual scalar solve(const scalarField& deltaT);
virtual scalar solveCMCchem
(
scalar deltaT, scalar& rho, scalar& T,
scalar& p, scalarField& Qi, scalar& Qh,
scalarField& RRCMC, scalar& ChemDeltaT
);
virtual scalar calculateTCMC(scalar& Qh, scalarField& Qi, scalar& Told, scalar& rho, scalar& p);
virtual scalar calculateRHOCMC(scalarField& Qi, scalar& T, scalar& p);
virtual scalar calculateHCMC(scalarField& Qi, scalar& T, scalar& rho, scalar& p);
virtual scalar calculateShCMC(scalarField& RRCMC, label i);
virtual void correction
(
scalarField& b, scalarField& Ta, scalarField& Wa, scalarField& Wb,
scalarField& Gab, scalarField& Gat, scalarField& Gbt,
scalarField& Gaa, scalarField& Gbb, scalar& Gtt, scalar& rhorho
);
//- Return the chemical time scale
virtual tmp<volScalarField> tc() const;
@ -280,10 +255,6 @@ public:
//- 28.Nov.2017 Jinwoo Park
virtual tmp<volScalarField> nut() const;
//- Calculate the reaction rates (EDM or finite-rate EDM)
//- 11.Apr.2017 Karam Han (updated by Jinwoo Park)
virtual void calculateEDM(bool finiteRate, const scalar A, const scalar B);
// ODE functions (overriding abstract functions in ODE.H)
//- Number of ODE's to solve