zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

ENH: Updated heat transfer solvers and tutorials to use p_rgh

Browse Source
This commit is contained in:
andy
2010-10-06 17:47:53 +01:00
parent aadacc5812
commit 590e2c3a52
128 changed files with 1255 additions and 1780 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/TEqn.H
View File

@ -12,7 +12,7 @@
);
TEqn.relax();
TEqn.solve();
TEqn.solve(mesh.solver(T.select(finalIter)));
rhok = 1.0 - beta*(T - TRef);
}

9
applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/UEqn.H
View File

@ -18,9 +18,10 @@
fvc::reconstruct
(
(
fvc::interpolate(rhok)*(g & mesh.Sf())
- fvc::snGrad(p)*mesh.magSf()
)
)
- ghf*fvc::snGrad(rhok)
- fvc::snGrad(p_rgh)
)*mesh.magSf()
),
mesh.solver(U.select(finalIter))
);
}

2
applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C
View File

@ -87,7 +87,7 @@ int main(int argc, char *argv[])
if (nOuterCorr != 1)
{
p.storePrevIter();
p_rgh.storePrevIter();
}
#include "UEqn.H"

56
applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/createFields.H
View File

@ -14,12 +14,12 @@
mesh
);
Info<< "Reading field p\n" << endl;
volScalarField p
Info<< "Reading field p_rgh\n" << endl;
volScalarField p_rgh
(
IOobject
(
"p",
"p_rgh",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
@ -52,6 +52,18 @@
incompressible::RASModel::New(U, phi, laminarTransport)
);
// Kinematic density for buoyancy force
volScalarField rhok
(
IOobject
(
"rhok",
runTime.timeName(),
mesh
),
1.0 - beta*(T - TRef)
);
// kinematic turbulent thermal thermal conductivity m2/s
Info<< "Reading field kappat\n" << endl;
volScalarField kappat
@ -67,25 +79,41 @@
mesh
);
Info<< "Calculating field g.h\n" << endl;
volScalarField gh("gh", g & mesh.C());
surfaceScalarField ghf("ghf", g & mesh.Cf());
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
p_rgh + rhok*gh
);
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell
(
p,
p_rgh,
mesh.solutionDict().subDict("PIMPLE"),
pRefCell,
pRefValue
);
// Kinematic density for buoyancy force
volScalarField rhok
(
IOobject
if (p_rgh.needReference())
{
p += dimensionedScalar
(
"rhok",
runTime.timeName(),
mesh
),
1.0 - beta*(T - TRef)
);
"p",
p.dimensions(),
pRefValue - getRefCellValue(p, pRefCell)
);
}

34
applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H
View File

@ -7,22 +7,23 @@
phi = (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, U, phi);
surfaceScalarField buoyancyPhi =
rUAf*fvc::interpolate(rhok)*(g & mesh.Sf());
phi += buoyancyPhi;
surfaceScalarField buoyancyPhi = rUAf*ghf*fvc::snGrad(rhok)*mesh.magSf();
phi -= buoyancyPhi;
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
fvScalarMatrix p_rghEqn
(
fvm::laplacian(rUAf, p) == fvc::div(phi)
fvm::laplacian(rUAf, p_rgh) == fvc::div(phi)
);
pEqn.solve
p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
p_rghEqn.solve
(
mesh.solver
(
p.select
p_rgh.select
(
(
finalIter
@ -36,17 +37,30 @@
if (nonOrth == nNonOrthCorr)
{
// Calculate the conservative fluxes
phi -= pEqn.flux();
phi -= p_rghEqn.flux();
// Explicitly relax pressure for momentum corrector
p.relax();
p_rgh.relax();
// Correct the momentum source with the pressure gradient flux
// calculated from the relaxed pressure
U += rUA*fvc::reconstruct((buoyancyPhi - pEqn.flux())/rUAf);
U -= rUA*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rUAf);
U.correctBoundaryConditions();
}
}
#include "continuityErrs.H"
p = p_rgh + rhok*gh;
if (p_rgh.needReference())
{
p += dimensionedScalar
(
"p",
p.dimensions(),
pRefValue - getRefCellValue(p, pRefCell)
);
p_rgh = p - rhok*gh;
}
}

18
applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H
View File

@ -96,29 +96,23 @@
p_rgh + rhok*gh
);
label p_rghRefCell = 0;
scalar p_rghRefValue = 0.0;
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell
(
p,
p_rgh,
mesh.solutionDict().subDict("SIMPLE"),
p_rghRefCell,
p_rghRefValue
pRefCell,
pRefValue
);
scalar pRefValue = 0.0;
if (p_rgh.needReference())
{
pRefValue = readScalar
(
mesh.solutionDict().subDict("SIMPLE").lookup("pRefValue")
);
p += dimensionedScalar
(
"p",
p.dimensions(),
pRefValue - getRefCellValue(p, p_rghRefCell)
pRefValue - getRefCellValue(p, pRefCell)
);
}

15
applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H
View File

@ -18,17 +18,9 @@
fvm::laplacian(rUAf, p_rgh) == fvc::div(phi)
);
p_rghEqn.setReference(p_rghRefCell, p_rghRefValue);
p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
// retain the residual from the first iteration
if (nonOrth == 0)
{
p_rghEqn.solve();
}
else
{
p_rghEqn.solve();
}
p_rghEqn.solve();
if (nonOrth == nNonOrthCorr)
{
@ -55,7 +47,8 @@
(
"p",
p.dimensions(),
pRefValue - getRefCellValue(p, p_rghRefCell)
pRefValue - getRefCellValue(p, pRefCell)
);
p_rgh = p - rhok*gh;
}
}

9
applications/solvers/heatTransfer/buoyantPimpleFoam/UEqn.H
View File

@ -17,8 +17,11 @@
==
fvc::reconstruct
(
fvc::interpolate(rho)*(g & mesh.Sf())
- fvc::snGrad(p)*mesh.magSf()
)
(
- ghf*fvc::snGrad(rho)
- fvc::snGrad(p_rgh)
)*mesh.magSf()
),
mesh.solver(U.select(finalIter))
);
}

2
applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
View File

@ -80,7 +80,7 @@ int main(int argc, char *argv[])
if (nOuterCorr != 1)
{
p.storePrevIter();
p_rgh.storePrevIter();
}
#include "UEqn.H"

27
applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H
View File

@ -53,15 +53,30 @@
)
);
Info<< "Calculating field g.h\n" << endl;
volScalarField gh("gh", g & mesh.C());
surfaceScalarField ghf("ghf", g & mesh.Cf());
Info<< "Reading field p_rgh\n" << endl;
volScalarField p_rgh
(
IOobject
(
"p_rgh",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
// Force p_rgh to be consistent with p
p_rgh = p - rho*gh;
Info<< "Creating field DpDt\n" << endl;
volScalarField DpDt
(
"DpDt",
fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p)
);
thermo.correct();
dimensionedScalar initialMass = fvc::domainIntegrate(rho);
dimensionedScalar totalVolume = sum(mesh.V());

2
applications/solvers/heatTransfer/buoyantPimpleFoam/hEqn.H
View File

@ -9,7 +9,7 @@
);
hEqn.relax();
hEqn.solve();
hEqn.solve(mesh.solver(h.select(finalIter)));
thermo.correct();
}

38
applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
View File

@ -1,11 +1,9 @@
{
bool closedVolume = p.needReference();
rho = thermo.rho();
// Thermodynamic density needs to be updated by psi*d(p) after the
// pressure solution - done in 2 parts. Part 1:
thermo.rho() -= psi*p;
thermo.rho() -= psi*p_rgh;
volScalarField rUA = 1.0/UEqn.A();
surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA));
@ -18,24 +16,23 @@
+ fvc::ddtPhiCorr(rUA, rho, U, phi)
);
surfaceScalarField buoyancyPhi =
rhorUAf*fvc::interpolate(rho)*(g & mesh.Sf());
surfaceScalarField buoyancyPhi = -rhorUAf*ghf*fvc::snGrad(rho)*mesh.magSf();
phi += buoyancyPhi;
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
fvScalarMatrix p_rghEqn
(
fvc::ddt(rho) + psi*correction(fvm::ddt(p))
fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
+ fvc::div(phi)
- fvm::laplacian(rhorUAf, p)
- fvm::laplacian(rhorUAf, p_rgh)
);
pEqn.solve
p_rghEqn.solve
(
mesh.solver
(
p.select
p_rgh.select
(
(
finalIter
@ -49,34 +46,25 @@
if (nonOrth == nNonOrthCorr)
{
// Calculate the conservative fluxes
phi += pEqn.flux();
phi += p_rghEqn.flux();
// Explicitly relax pressure for momentum corrector
p.relax();
p_rgh.relax();
// Correct the momentum source with the pressure gradient flux
// calculated from the relaxed pressure
U += rUA*fvc::reconstruct((buoyancyPhi + pEqn.flux())/rhorUAf);
U += rUA*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rhorUAf);
U.correctBoundaryConditions();
}
}
p = p_rgh + rho*gh;
// Second part of thermodynamic density update
thermo.rho() += psi*p;
thermo.rho() += psi*p_rgh;
DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
#include "rhoEqn.H"
#include "compressibleContinuityErrs.H"
// For closed-volume cases adjust the pressure and density levels
// to obey overall mass continuity
if (closedVolume)
{
p +=
(initialMass - fvc::domainIntegrate(psi*p))
/fvc::domainIntegrate(psi);
thermo.rho() = psi*p;
rho += (initialMass - fvc::domainIntegrate(rho))/totalVolume;
}
}

3
applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C
View File

@ -62,10 +62,7 @@ int main(int argc, char *argv[])
{
#include "UEqn.H"
#include "hEqn.H"
for (int i=0; i<3; i++)
{
#include "pEqn.H"
}
}
turbulence->correct();

68
applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H
View File

@ -23,20 +23,6 @@
volScalarField& h = thermo.h();
const volScalarField& psi = thermo.psi();
Info<< "Reading field p_rgh\n" << endl;
volScalarField p_rgh
(
IOobject
(
"p_rgh",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading field U\n" << endl;
volVectorField U
(
@ -53,7 +39,6 @@
#include "compressibleCreatePhi.H"
Info<< "Creating turbulence model\n" << endl;
autoPtr<compressible::RASModel> turbulence
(
@ -66,40 +51,39 @@
)
);
Info<< "Calculating field g.h\n" << endl;
volScalarField gh("gh", g & mesh.C());
surfaceScalarField ghf("ghf", g & mesh.Cf());
p = p_rgh + rho*gh;
thermo.correct();
rho = thermo.rho();
p_rgh = p - rho*gh;
label p_rghRefCell = 0;
scalar p_rghRefValue = 0.0;
setRefCell
Info<< "Reading field p_rgh\n" << endl;
volScalarField p_rgh
(
p_rgh,
mesh.solutionDict().subDict("SIMPLE"),
p_rghRefCell,
p_rghRefValue
IOobject
(
"p_rgh",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
// Force p_rgh to be consistent with p
p_rgh = p - rho*gh;
label pRefCell = 0;
scalar pRefValue = 0.0;
if (p_rgh.needReference())
{
pRefValue = readScalar
(
mesh.solutionDict().subDict("SIMPLE").lookup("pRefValue")
);
p += dimensionedScalar
(
"p",
p.dimensions(),
pRefValue - getRefCellValue(p, p_rghRefCell)
);
}
setRefCell
(
p,
p_rgh,
mesh.solutionDict().subDict("SIMPLE"),
pRefCell,
pRefValue
);
dimensionedScalar initialMass = fvc::domainIntegrate(rho);
dimensionedScalar totalVolume = sum(mesh.V());

9
applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
View File

@ -1,11 +1,12 @@
{
rho = thermo.rho();
rho.relax();
volScalarField rUA = 1.0/UEqn().A();
surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA));
U = rUA*UEqn().H();
//UEqn.clear();
UEqn.clear();
phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
bool closedVolume = adjustPhi(phi, U, p_rgh);
@ -20,7 +21,7 @@
fvm::laplacian(rhorUAf, p_rgh) == fvc::div(phi)
);
p_rghEqn.setReference(p_rghRefCell, p_rghRefValue);
p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
p_rghEqn.solve();
if (nonOrth == nNonOrthCorr)
@ -42,13 +43,13 @@
p = p_rgh + rho*gh;
// For closed-volume cases adjust the pressure and density levels
// For closed-volume cases adjust the pressure level
// to obey overall mass continuity
if (closedVolume)
{
p += (initialMass - fvc::domainIntegrate(psi*p))
/fvc::domainIntegrate(psi);
p_rgh == p - rho*gh;
p_rgh = p - rho*gh;
}
rho = thermo.rho();

2
applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C
View File

@ -58,7 +58,7 @@ int main(int argc, char *argv[])
#include "readSIMPLEControls.H"
p.storePrevIter();
p_rgh.storePrevIter();
rho.storePrevIter();
// Pressure-velocity SIMPLE corrector

2
applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
View File

@ -93,6 +93,8 @@ int main(int argc, char *argv[])
// --- PIMPLE loop
for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
{
bool finalIter = oCorr == nOuterCorr-1;
forAll(fluidRegions, i)
{
Info<< "\nSolving for fluid region "

6
applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/UEqn.H
View File

@ -13,7 +13,9 @@
==
fvc::reconstruct
(
fvc::interpolate(rho)*(g & mesh.Sf())
- fvc::snGrad(p)*mesh.magSf()
(
- ghf*fvc::snGrad(rho)
- fvc::snGrad(p_rgh)
)*mesh.magSf()
)
);

59
applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidFields.H
View File

@ -6,12 +6,16 @@
PtrList<surfaceScalarField> phiFluid(fluidRegions.size());
PtrList<uniformDimensionedVectorField> gFluid(fluidRegions.size());
PtrList<compressible::turbulenceModel> turbulence(fluidRegions.size());
PtrList<volScalarField> DpDtf(fluidRegions.size());
PtrList<volScalarField> p_rghFluid(fluidRegions.size());
PtrList<volScalarField> ghFluid(fluidRegions.size());
PtrList<surfaceScalarField> ghfFluid(fluidRegions.size());
List<scalar> initialMassFluid(fluidRegions.size());
List<label> pRefCellFluid(fluidRegions.size(),0);
List<scalar> pRefValueFluid(fluidRegions.size(),0.0);
PtrList<dimensionedScalar> rhoMax(fluidRegions.size());
PtrList<dimensionedScalar> rhoMin(fluidRegions.size());
// Populate fluid field pointer lists
forAll(fluidRegions, i)
@ -130,15 +134,68 @@
).ptr()
);
Info<< " Adding to ghFluid\n" << endl;
ghFluid.set
(
i,
new volScalarField("gh", gFluid[i] & fluidRegions[i].C())
);
Info<< " Adding to ghfFluid\n" << endl;
ghfFluid.set
(
i,
new surfaceScalarField("ghf", gFluid[i] & fluidRegions[i].Cf())
);
p_rghFluid.set
(
i,
new volScalarField
(
IOobject
(
"p_rgh",
runTime.timeName(),
fluidRegions[i],
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
fluidRegions[i]
)
);
// Force p_rgh to be consistent with p
p_rghFluid[i] = thermoFluid[i].p() - rhoFluid[i]*ghFluid[i];
initialMassFluid[i] = fvc::domainIntegrate(rhoFluid[i]).value();
setRefCell
(
thermoFluid[i].p(),
p_rghFluid[i],
fluidRegions[i].solutionDict().subDict("SIMPLE"),
pRefCellFluid[i],
pRefValueFluid[i]
);
rhoMax.set
(
i,
new dimensionedScalar
(
fluidRegions[i].solutionDict().subDict("SIMPLE").lookup("rhoMax")
)
);
rhoMin.set
(
i,
new dimensionedScalar
(
fluidRegions[i].solutionDict().subDict("SIMPLE").lookup("rhoMin")
)
);
}

4
applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/hEqn.H
View File

@ -5,8 +5,8 @@
- fvm::Sp(fvc::div(phi), h)
- fvm::laplacian(turb.alphaEff(), h)
==
fvc::div(phi/fvc::interpolate(rho)*fvc::interpolate(p))
- p*fvc::div(phi/fvc::interpolate(rho))
fvc::div(phi/fvc::interpolate(rho), rho/psi, "div(U,p)")
- (rho/psi)*fvc::div(phi/fvc::interpolate(rho))
);
hEqn.relax();

53
applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
View File

@ -1,7 +1,9 @@
{
// From buoyantSimpleFoam
rho = thermo.rho();
rho = max(rho, rhoMin[i]);
rho = min(rho, rhoMax[i]);
rho.relax();
volScalarField rUA = 1.0/UEqn().A();
surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA));
@ -10,59 +12,54 @@
UEqn.clear();
phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
bool closedVolume = adjustPhi(phi, U, p);
bool closedVolume = adjustPhi(phi, U, p_rgh);
surfaceScalarField buoyancyPhi =
rhorUAf*fvc::interpolate(rho)*(g & mesh.Sf());
phi += buoyancyPhi;
surfaceScalarField buoyancyPhi = rhorUAf*ghf*fvc::snGrad(rho)*mesh.magSf();
phi -= buoyancyPhi;
// Solve pressure
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
fvScalarMatrix p_rghEqn
(
fvm::laplacian(rhorUAf, p) == fvc::div(phi)
fvm::laplacian(rhorUAf, p_rgh) == fvc::div(phi)
);
pEqn.setReference(pRefCell, pRefValue);
p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
// retain the residual from the first iteration
if (nonOrth == 0)
{
pEqn.solve();
}
else
{
pEqn.solve();
}
p_rghEqn.solve();
if (nonOrth == nNonOrthCorr)
{
// For closed-volume cases adjust the pressure and density levels
// to obey overall mass continuity
if (closedVolume)
{
p += (initialMass - fvc::domainIntegrate(psi*p))
/fvc::domainIntegrate(psi);
}
// Calculate the conservative fluxes
phi -= pEqn.flux();
phi -= p_rghEqn.flux();
// Explicitly relax pressure for momentum corrector
p.relax();
p_rgh.relax();
// Correct the momentum source with the pressure gradient flux
// calculated from the relaxed pressure
U += rUA*fvc::reconstruct((buoyancyPhi - pEqn.flux())/rhorUAf);
U -= rUA*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rhorUAf);
U.correctBoundaryConditions();
}
}
p = p_rgh + rho*gh;
#include "continuityErrs.H"
// For closed-volume cases adjust the pressure level
// to obey overall mass continuity
if (closedVolume)
{
p += (initialMass - fvc::domainIntegrate(psi*p))
/fvc::domainIntegrate(psi);
p_rgh = p - rho*gh;
}
rho = thermo.rho();
rho = max(rho, rhoMin[i]);
rho = min(rho, rhoMax[i]);
rho.relax();
Info<< "Min/max rho:" << min(rho).value() << ' '

5
applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H
View File

@ -5,7 +5,6 @@
volScalarField& K = KFluid[i];
volVectorField& U = UFluid[i];
surfaceScalarField& phi = phiFluid[i];
const dimensionedVector& g = gFluid[i];
compressible::turbulenceModel& turb = turbulence[i];
@ -22,3 +21,7 @@
const label pRefCell = pRefCellFluid[i];
const scalar pRefValue = pRefValueFluid[i];
volScalarField& p_rgh = p_rghFluid[i];
const volScalarField& gh = ghFluid[i];
const surfaceScalarField& ghf = ghfFluid[i];

2
applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/solveFluid.H
View File

@ -1,6 +1,6 @@
// Pressure-velocity SIMPLE corrector
p.storePrevIter();
p_rgh.storePrevIter();
rho.storePrevIter();
{
#include "UEqn.H"

9
applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H
View File

@ -16,8 +16,11 @@
==
fvc::reconstruct
(
fvc::interpolate(rho)*(g & mesh.Sf())
- fvc::snGrad(p)*mesh.magSf()
)
(
- ghf*fvc::snGrad(rho)
- fvc::snGrad(p_rgh)
)*mesh.magSf()
),
mesh.solver(U.select(finalIter))
);
}

41
applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
View File

@ -6,6 +6,9 @@
PtrList<surfaceScalarField> phiFluid(fluidRegions.size());
PtrList<uniformDimensionedVectorField> gFluid(fluidRegions.size());
PtrList<compressible::turbulenceModel> turbulence(fluidRegions.size());
PtrList<volScalarField> p_rghFluid(fluidRegions.size());
PtrList<volScalarField> ghFluid(fluidRegions.size());
PtrList<surfaceScalarField> ghfFluid(fluidRegions.size());
PtrList<volScalarField> DpDtFluid(fluidRegions.size());
List<scalar> initialMassFluid(fluidRegions.size());
@ -129,6 +132,42 @@
).ptr()
);
Info<< " Adding to ghFluid\n" << endl;
ghFluid.set
(
i,
new volScalarField("gh", gFluid[i] & fluidRegions[i].C())
);
Info<< " Adding to ghfFluid\n" << endl;
ghfFluid.set
(
i,
new surfaceScalarField("ghf", gFluid[i] & fluidRegions[i].Cf())
);
p_rghFluid.set
(
i,
new volScalarField
(
IOobject
(
"p_rgh",
runTime.timeName(),
fluidRegions[i],
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
fluidRegions[i]
)
);
// Force p_rgh to be consistent with p
p_rghFluid[i] = thermoFluid[i].p() - rhoFluid[i]*ghFluid[i];
initialMassFluid[i] = fvc::domainIntegrate(rhoFluid[i]).value();
Info<< " Adding to DpDtFluid\n" << endl;
DpDtFluid.set
(
@ -147,6 +186,4 @@
)
)
);
initialMassFluid[i] = fvc::domainIntegrate(rhoFluid[i]).value();
}

13
applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/hEqn.H
View File

@ -7,16 +7,9 @@
==
DpDt
);
if (oCorr == nOuterCorr-1)
{
hEqn.relax();
hEqn.solve(mesh.solver("hFinal"));
}
else
{
hEqn.relax();
hEqn.solve();
}
hEqn.relax();
hEqn.solve(mesh.solver(h.select(finalIter)));
thermo.correct();

42
applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
View File

@ -1,5 +1,5 @@
{
bool closedVolume = p.needReference();
bool closedVolume = p_rgh.needReference();
rho = thermo.rho();
@ -17,34 +17,35 @@
)
);
phi = phiU + fvc::interpolate(rho)*(g & mesh.Sf())*rhorUAf;
phi = phiU - rhorUAf*ghf*fvc::snGrad(rho)*mesh.magSf();
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
fvScalarMatrix p_rghEqn
(
fvm::ddt(psi, p)
fvm::ddt(psi, p_rgh) + fvc::ddt(psi, rho)*gh
+ fvc::div(phi)
- fvm::laplacian(rhorUAf, p)
- fvm::laplacian(rhorUAf, p_rgh)
);
if
p_rghEqn.solve
(
oCorr == nOuterCorr-1
&& corr == nCorr-1
&& nonOrth == nNonOrthCorr
)
{
pEqn.solve(mesh.solver(p.name() + "Final"));
}
else
{
pEqn.solve(mesh.solver(p.name()));
}
mesh.solver
(
p_rgh.select
(
(
oCorr == nOuterCorr-1
&& corr == nCorr-1
&& nonOrth == nNonOrthCorr
)
)
)
);
if (nonOrth == nNonOrthCorr)
{
phi += pEqn.flux();
phi += p_rghEqn.flux();
}
}
@ -52,6 +53,8 @@
U += rUA*fvc::reconstruct((phi - phiU)/rhorUAf);
U.correctBoundaryConditions();
p = p_rgh + rho*gh;
// Update pressure substantive derivative
DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
@ -65,9 +68,10 @@
// to obey overall mass continuity
if (closedVolume)
{
p += (massIni - fvc::domainIntegrate(psi*p))
p += (initialMass - fvc::domainIntegrate(psi*p))
/fvc::domainIntegrate(psi);
rho = thermo.rho();
p_rgh = p - rho*gh;
}
// Update thermal conductivity

17
applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/readFluidMultiRegionPISOControls.H
View File

@ -1,17 +0,0 @@
const dictionary& piso = fluidRegions[i].solutionDict().subDict("PISO");
const int nOuterCorr =
piso.lookupOrDefault<int>("nOuterCorrectors", 1);
const int nCorr =
piso.lookupOrDefault<int>("nCorrectors", 1);
const int nNonOrthCorr =
piso.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0);
const bool momentumPredictor =
piso.lookupOrDefault("momentumPredictor", true);
const bool transonic =
piso.lookupOrDefault("transonic", false);

14
applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
View File

@ -1,11 +1,10 @@
const fvMesh& mesh = fluidRegions[i];
fvMesh& mesh = fluidRegions[i];
basicPsiThermo& thermo = thermoFluid[i];
volScalarField& rho = rhoFluid[i];
volScalarField& K = KFluid[i];
volVectorField& U = UFluid[i];
surfaceScalarField& phi = phiFluid[i];
const dimensionedVector& g = gFluid[i];
compressible::turbulenceModel& turb = turbulence[i];
volScalarField& DpDt = DpDtFluid[i];
@ -14,4 +13,13 @@
const volScalarField& psi = thermo.psi();
volScalarField& h = thermo.h();
const dimensionedScalar massIni("massIni", dimMass, initialMassFluid[i]);
volScalarField& p_rgh = p_rghFluid[i];
const volScalarField& gh = ghFluid[i];
const surfaceScalarField& ghf = ghfFluid[i];
const dimensionedScalar initialMass
(
"initialMass",
dimMass,
initialMassFluid[i]
);

10
applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H
View File

@ -1,3 +1,8 @@
if (finalIter)
{
mesh.data::add("finalIteration", true);
}
if (oCorr == 0)
{
#include "rhoEqn.H"
@ -16,3 +21,8 @@ for (int corr=0; corr<nCorr; corr++)
turb.correct();
rho = thermo.rho();
if (finalIter)
{
mesh.data::remove("finalIteration");
}

2
applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/storeOldFluidFields.H
View File

@ -1,2 +1,2 @@
p.storePrevIter();
p_rgh.storePrevIter();
rho.storePrevIter();

6
applications/solvers/heatTransfer/chtMultiRegionFoam/solid/createSolidFields.H
View File

@ -8,9 +8,5 @@
<< solidRegions[i].name() << nl << endl;
Info<< " Adding to thermos\n" << endl;
thermos.set
(
i,
basicSolidThermo::New(solidRegions[i])
);
thermos.set(i, basicSolidThermo::New(solidRegions[i]));
}

12
applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H
View File

@ -1,3 +1,8 @@
if (finalIter)
{
mesh.data::add("finalIteration", true);
}
{
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
@ -7,8 +12,13 @@
- fvm::laplacian(K, T)
);
TEqn().relax();
TEqn().solve();
TEqn().solve(mesh.solver(T.select(finalIter)));
}
Info<< "Min/max T:" << min(T) << ' ' << max(T) << endl;
}
if (finalIter)
{
mesh.data::remove("finalIteration");
}

2
applications/solvers/lagrangian/coalChemistryFoam/UEqn.H
View File

@ -13,5 +13,5 @@
if (momentumPredictor)
{
solve(UEqn == -fvc::grad(p));
solve(UEqn == -fvc::grad(p), mesh.solver(U.select(finalIter)));
}

15
applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C
View File

@ -90,17 +90,28 @@ int main(int argc, char *argv[])
#include "rhoEqn.H"
// --- PIMPLE loop
for (int ocorr=1; ocorr<=nOuterCorr; ocorr++)
for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
{
bool finalIter = oCorr == nOuterCorr - 1;
if (finalIter)
{
mesh.data::add("finalIteration", true);
}
#include "UEqn.H"
#include "YEqn.H"
#include "hsEqn.H"
// --- PISO loop
for (int corr=1; corr<=nCorr; corr++)
for (int corr=0; corr<nCorr; corr++)
{
#include "pEqn.H"
}
if (finalIter)
{
mesh.data::remove("finalIteration");
}
}
turbulence->correct();

2
applications/solvers/lagrangian/coalChemistryFoam/hsEqn.H
View File

@ -15,7 +15,7 @@
hsEqn.relax();
hsEqn.solve();
hsEqn.solve(mesh.solver(hs.select(finalIter)));
thermo.correct();

30
applications/solvers/lagrangian/coalChemistryFoam/pEqn.H
View File

@ -26,7 +26,20 @@ if (transonic)
coalParcels.Srho()
);
pEqn.solve();
pEqn.solve
(
mesh.solver
(
p.select
(
(
finalIter
&& corr == nCorr-1
&& nonOrth == nNonOrthCorr
)
)
)
);
if (nonOrth == nNonOrthCorr)
{
@ -54,7 +67,20 @@ else
coalParcels.Srho()
);
pEqn.solve();
pEqn.solve
(
mesh.solver
(
p.select
(
(
finalIter
&& corr == nCorr-1
&& nonOrth == nNonOrthCorr
)
)
)
);
if (nonOrth == nNonOrthCorr)
{