85 lines
2 KiB
C
85 lines
2 KiB
C
{
|
|
tmp<fv::convectionScheme<scalar>> hConvection
|
|
(
|
|
fv::convectionScheme<scalar>::New
|
|
(
|
|
mesh,
|
|
phi,
|
|
mesh.divScheme("div(phi,Yi_h)")
|
|
)
|
|
);
|
|
|
|
typedef multiComponentMixture<gasHThermoPhysics> MMix;
|
|
|
|
MMix &janafComposition = dynamic_cast<MMix&>(composition);
|
|
|
|
volScalarField& he = thermo.he();
|
|
|
|
volScalarField hsi(he);
|
|
|
|
volVectorField hDiffusionSrc("rho*Sum(hYV)", he * Vc * 0.0);
|
|
|
|
forAll(Y, i)
|
|
{
|
|
const MMix::thermoType &tti = janafComposition.getLocalThermo(i);
|
|
|
|
forAll(hsi, cellI)
|
|
{
|
|
const scalar pI = p[cellI];
|
|
const scalar TI = T[cellI];
|
|
|
|
hsi[cellI] = tti.HE(pI,TI);
|
|
}
|
|
|
|
forAll(hsi.boundaryFieldRef(), patchI)
|
|
{
|
|
volScalarField::Patch &hsiP = hsi.boundaryFieldRef()[patchI];
|
|
const volScalarField::Patch &pP = p.boundaryField()[patchI];
|
|
const volScalarField::Patch &TP = T.boundaryField()[patchI];
|
|
|
|
forAll(hsiP, faceI)
|
|
{
|
|
const scalar pI = pP[faceI];
|
|
const scalar TI = TP[faceI];
|
|
|
|
hsiP[faceI] = tti.HE(pI,TI);
|
|
}
|
|
}
|
|
|
|
hDiffusionSrc += hsi * (- diff.D(i) * fvc::grad(Y[i]) + Y[i] * Vc);
|
|
}
|
|
|
|
tmp<volVectorField> tTauU(diff.mu()*(U&(Foam::dev2(Foam::T(gradU)) + gradU)));
|
|
|
|
hDiffusionSrc *= rho;
|
|
|
|
volVectorField hWorkSrc1(rho * U);
|
|
|
|
volVectorField hWorkSrc2(rho * Vc);
|
|
|
|
fvScalarMatrix EEqn
|
|
(
|
|
fvm::ddt(rho, he) + hConvection->fvmDiv(phi, he)
|
|
+ fvc::ddt(rho, K) + fvc::div(phi, K)
|
|
- dpdt
|
|
- fvc::laplacian(diff.k(), T)
|
|
+ fvc::div(hDiffusionSrc)
|
|
==
|
|
reaction->Sh()
|
|
+ fvc::div(tTauU)
|
|
+ fvOptions(rho, he)
|
|
);
|
|
|
|
EEqn.relax();
|
|
|
|
fvOptions.constrain(EEqn);
|
|
|
|
EEqn.solve();
|
|
|
|
fvOptions.correct(he);
|
|
|
|
thermo.correct();
|
|
|
|
Info<< "min/max(T) = "
|
|
<< min(T).value() << ", " << max(T).value() << endl;
|
|
}
|