Updated enthalpy equation

Conflicts:
	applications/solvers/compressible/steadyCompressibleFoam/pEqn.H
This commit is contained in:
Hrvoje Jasak 2015-12-04 19:27:05 +00:00
parent 6b6a9cbb13
commit 9daa8e28f4
10 changed files with 42 additions and 51 deletions

View file

@ -14,9 +14,10 @@
fvm::ddt(rho, h)
+ fvm::div(phi, h)
- fvm::laplacian(turbulence->alphaEff(), h)
+ fvm::SuSp((fvc::div(faceU, p, "div(U,p)") - p*fvc::div(faceU))/h, h)
==
// ddt(p) term removed: steady-state. HJ, 27/Apr/2010
fvc::div(faceU, p, "div(U,p)")
- p*fvc::div(faceU)
// Viscous heating: note sign (devRhoReff has a minus in it)
- (turbulence->devRhoReff() && fvc::grad(U))
);
@ -26,12 +27,6 @@
eqnResidual = hEqn.solve().initialResidual();
maxResidual = max(eqnResidual, maxResidual);
// Bound the enthalpy using TMin and TMax
volScalarField Cp = thermo.Cp();
h = Foam::min(h, TMax*Cp);
h = Foam::max(h, TMin*Cp);
h.correctBoundaryConditions();
// Bounding of enthalpy taken out
thermo.correct();
}

View file

@ -1,6 +1,15 @@
{
volScalarField rUA = 1.0/UEqn.A();
<<<<<<< HEAD
=======
surfaceScalarField psisf = fvc::interpolate(psis);
surfaceScalarField rhof = fvc::interpolate(rho);
// Needs to be outside of loop since p is changing, but psi and rho are not
surfaceScalarField rhoReff = rhof - psisf*fvc::interpolate(p);
>>>>>>> 07efd81... Updated enthalpy equation
for (int corr = 0; corr < nCorr; corr++)
{
U = rUA*UEqn.H();

View file

@ -2,19 +2,19 @@
dictionary fieldBounds = mesh.solutionDict().subDict("fieldBounds");
// Pressure bounds
dimensionedScalar pMin("pMin", dimPressure, 0);
dimensionedScalar pMax("pMax", dimPressure, GREAT);
dimensionedScalar pMin("pMin", p.dimensions(), 0);
dimensionedScalar pMax("pMax", p.dimensions(), GREAT);
fieldBounds.lookup("p") >> pMin.value() >> pMax.value();
fieldBounds.lookup(p.name()) >> pMin.value() >> pMax.value();
// Temperature bounds
dimensionedScalar TMin("TMin", dimTemperature, 0);
dimensionedScalar TMax("TMax", dimTemperature, GREAT);
dimensionedScalar TMin("TMin", T.dimensions(), 0);
dimensionedScalar TMax("TMax", T.dimensions(), GREAT);
fieldBounds.lookup("T") >> TMin.value() >> TMax.value();
fieldBounds.lookup(T.name()) >> TMin.value() >> TMax.value();
// Velocity bound
dimensionedScalar UMax("UMax", dimVelocity, GREAT);
dimensionedScalar UMax("UMax", U.dimensions(), GREAT);
fieldBounds.lookup(U.name()) >> UMax.value();
dimensionedScalar smallU("smallU", dimVelocity, 1e-10);

View file

@ -1,7 +1,7 @@
{
// Calculate density from pressure
rho.storePrevIter();
rho = thermo.rho();
rho = thermo.rho()();
// Bound rho
volScalarField R = thermo.Cp() - thermo.Cv();

View file

@ -10,10 +10,9 @@
+ fvm::div(phi, h)
- fvm::laplacian(turbulence->alphaEff(), h)
==
// Note: potential issue with reconstructed relative velocity.
// HJ, 12/Dec/2009
// fvc::div(faceU, p, "div(U,p)")
fvc::div(faceU + mrfZones.fluxCorrection(), p, "div(U,p)")
// Note: div flux correction is zero so there is no need to
// carry it. HJ, 4/Dec/2015
- p*fvc::div(faceU)
// Viscous heating: note sign (devRhoReff has a minus in it)
- (turbulence->devRhoReff() && fvc::grad(U))
@ -24,12 +23,6 @@
eqnResidual = hEqn.solve().initialResidual();
maxResidual = max(eqnResidual, maxResidual);
// Bound the enthalpy using TMin and TMax
volScalarField Cp = thermo.Cp();
h = Foam::min(h, TMax*Cp);
h = Foam::max(h, TMin*Cp);
h.correctBoundaryConditions();
// Bounding of enthalpy taken out
thermo.correct();
}

View file

@ -2,19 +2,19 @@
dictionary fieldBounds = mesh.solutionDict().subDict("fieldBounds");
// Pressure bounds
dimensionedScalar pMin("pMin", dimPressure, 0);
dimensionedScalar pMax("pMax", dimPressure, GREAT);
dimensionedScalar pMin("pMin", p.dimensions(), 0);
dimensionedScalar pMax("pMax", p.dimensions(), GREAT);
fieldBounds.lookup("p") >> pMin.value() >> pMax.value();
fieldBounds.lookup(p.name()) >> pMin.value() >> pMax.value();
// Temperature bounds
dimensionedScalar TMin("TMin", dimTemperature, 0);
dimensionedScalar TMax("TMax", dimTemperature, GREAT);
dimensionedScalar TMin("TMin", T.dimensions(), 0);
dimensionedScalar TMax("TMax", T.dimensions(), GREAT);
fieldBounds.lookup("T") >> TMin.value() >> TMax.value();
fieldBounds.lookup(T.name()) >> TMin.value() >> TMax.value();
// Velocity bound
dimensionedScalar UMax("UMax", dimVelocity, GREAT);
dimensionedScalar UMax("UMax", U.dimensions(), GREAT);
fieldBounds.lookup(U.name()) >> UMax.value();
dimensionedScalar smallU("smallU", dimVelocity, 1e-10);

View file

@ -1,7 +1,7 @@
{
// Calculate density from pressure
rho.storePrevIter();
rho = thermo.rho();
rho = thermo.rho()();
// Bound rho
volScalarField R = thermo.Cp() - thermo.Cv();

View file

@ -14,10 +14,10 @@
fvm::ddt(rho, i)
+ fvm::div(phi, i)
- fvm::laplacian(turbulence->alphaEff(), i)
// u & gradP term (steady-state formulation)
+ fvm::SuSp((fvc::div(faceU, p, "div(U,p)") - p*fvc::div(faceU))/i, i)
==
// ddt(p) term removed: steady-state. HJ, 27/Apr/2010
fvc::div(faceU, p, "div(U,p)")
- p*fvc::div(faceU))
// Viscous heating: note sign (devRhoReff has a minus in it)
- (turbulence->devRhoReff() && fvc::grad(Urel))
);
@ -33,15 +33,9 @@
h = i + 0.5*magSqr(Urot);
h.correctBoundaryConditions();
// Bound the enthalpy using TMin and TMax
volScalarField Cp = thermo.Cp();
h = Foam::min(h, TMax*Cp);
h = Foam::max(h, TMin*Cp);
h.correctBoundaryConditions();
// Re-initialise rothalpy based on limited enthalpy
i = h - 0.5*magSqr(Urot);
// Bounding of enthalpy taken out
thermo.correct();
}

View file

@ -2,19 +2,19 @@
dictionary fieldBounds = mesh.solutionDict().subDict("fieldBounds");
// Pressure bounds
dimensionedScalar pMin("pMin", dimPressure, 0);
dimensionedScalar pMax("pMax", dimPressure, GREAT);
dimensionedScalar pMin("pMin", p.dimensions(), 0);
dimensionedScalar pMax("pMax", p.dimensions(), GREAT);
fieldBounds.lookup("p") >> pMin.value() >> pMax.value();
fieldBounds.lookup(p.name()) >> pMin.value() >> pMax.value();
// Temperature bounds
dimensionedScalar TMin("TMin", dimTemperature, 0);
dimensionedScalar TMax("TMax", dimTemperature, GREAT);
dimensionedScalar TMin("TMin", T.dimensions(), 0);
dimensionedScalar TMax("TMax", T.dimensions(), GREAT);
fieldBounds.lookup("T") >> TMin.value() >> TMax.value();
fieldBounds.lookup(T.name()) >> TMin.value() >> TMax.value();
// Velocity bound
dimensionedScalar UrelMax("UrelMax", dimVelocity, GREAT);
dimensionedScalar UrelMax("UrelMax", Urel.dimensions(), GREAT);
fieldBounds.lookup(Urel.name()) >> UrelMax.value();
dimensionedScalar smallUrel("smallUrel", dimVelocity, 1e-10);

View file

@ -1,7 +1,7 @@
{
// Calculate density from pressure
rho.storePrevIter();
rho = thermo.rho();
rho = thermo.rho()();
// Bound rho
volScalarField R = thermo.Cp() - thermo.Cv();