Files
openfoam/applications/solvers/combustion/chemFoam/readInitialConditions.H

109 lines
2.7 KiB
C

word constProp(initialConditions.get<word>("constantProperty"));
if (constProp != "pressure" && constProp != "volume")
{
FatalError << "in initialConditions, unknown constantProperty type "
<< constProp << nl << " Valid types are: pressure volume."
<< abort(FatalError);
}
word fractionBasis(initialConditions.get<word>("fractionBasis"));
if (fractionBasis != "mass" && fractionBasis != "mole")
{
FatalError << "in initialConditions, unknown fractionBasis type " << nl
<< "Valid types are: mass or mole."
<< fractionBasis << abort(FatalError);
}
label nSpecie = Y.size();
PtrList<gasHThermoPhysics> specieData(Y.size());
forAll(specieData, i)
{
specieData.set
(
i,
new gasHThermoPhysics
(
dynamic_cast<const reactingMixture<gasHThermoPhysics>&>
(thermo).speciesData()[i]
)
);
}
scalarList Y0(nSpecie, 0.0);
scalarList X0(nSpecie, 0.0);
dictionary fractions(initialConditions.subDict("fractions"));
if (fractionBasis == "mole")
{
forAll(Y, i)
{
const word& name = Y[i].name();
if (fractions.found(name))
{
X0[i] = fractions.get<scalar>(name);
}
}
scalar mw = 0.0;
const scalar mTot = sum(X0);
forAll(Y, i)
{
X0[i] /= mTot;
mw += specieData[i].W()*X0[i];
}
forAll(Y, i)
{
Y0[i] = X0[i]*specieData[i].W()/mw;
}
}
else // mass fraction
{
forAll(Y, i)
{
const word& name = Y[i].name();
if (fractions.found(name))
{
Y0[i] = fractions.get<scalar>(name);
}
}
scalar invW = 0.0;
const scalar mTot = sum(Y0);
forAll(Y, i)
{
Y0[i] /= mTot;
invW += Y0[i]/specieData[i].W();
}
const scalar mw = 1.0/invW;
forAll(Y, i)
{
X0[i] = Y0[i]*mw/specieData[i].W();
}
}
scalar h0 = 0.0;
forAll(Y, i)
{
Y[i] = Y0[i];
h0 += Y0[i]*specieData[i].Hs(p[0], T0);
}
thermo.he() = dimensionedScalar("h", dimEnergy/dimMass, h0);
thermo.correct();
rho = thermo.rho();
scalar rho0 = rho[0];
scalar u0 = h0 - p0/rho0;
scalar R0 = p0/(rho0*T0);
Rspecific[0] = R0;
scalar integratedHeat = 0.0;
Info << constProp << " will be held constant." << nl
<< " p = " << p[0] << " [Pa]" << nl
<< " T = " << thermo.T()[0] << " [K] " << nl
<< " rho = " << rho[0] << " [kg/m3]" << nl
<< endl;