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

Improved the flux-velocity correspondence for cases where hydrostatic balance is

Browse Source 
important e.g. in atmospheric flows.
This commit is contained in:
henry
2009-02-03 16:51:07 +00:00
parent 6c36054c43
commit 338b72b1eb
11 changed files with 151 additions and 180 deletions
Download Patch File Download Diff File Expand all files Collapse all files

16
applications/solvers/heatTransfer/buoyantFoam/UEqn.H
View File

@ -9,4 +9,18 @@
UEqn.relax();
solve(UEqn == -fvc::grad(pd) - fvc::grad(rho)*gh);
if (momentumPredictor)
{
solve
(
UEqn
==
-fvc::reconstruct
(
(
fvc::snGrad(pd)
+ ghf*fvc::snGrad(rho)
) * mesh.magSf()
)
);
}

4
applications/solvers/heatTransfer/buoyantFoam/buoyantFoam.C
View File

@ -41,7 +41,6 @@ Description
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
@ -55,9 +54,10 @@ int main(int argc, char *argv[])
while (runTime.run())
{
#include "readTimeControls.H"
#include "readPISOControls.H"
#include "compressibleCourantNo.H"
//# include "setDeltaT.H"
#include "setDeltaT.H"
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;

1
applications/solvers/heatTransfer/buoyantFoam/createFields.H
View File

@ -57,6 +57,7 @@
Info<< "Calculating field g.h\n" << endl;
volScalarField gh("gh", g & mesh.C());
surfaceScalarField ghf("gh", g & mesh.Cf());
dimensionedScalar pRef("pRef", p.dimensions(), thermo->lookup("pRef"));

23
applications/solvers/heatTransfer/buoyantFoam/pEqn.H
View File

@ -1,17 +1,23 @@
{
bool closedVolume = pd.needReference();
rho = thermo->rho();
volScalarField rUA = 1.0/UEqn.A();
surfaceScalarField rhorUAf = fvc::interpolate(rho*rUA);
U = rUA*UEqn.H();
phi =
surfaceScalarField phiU
(
fvc::interpolate(rho)
*(
(fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi)
)
- fvc::interpolate(rho*rUA*gh)*fvc::snGrad(rho)*mesh.magSf();
);
phi = phiU - ghf*fvc::snGrad(rho)*rhorUAf*mesh.magSf();
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
@ -21,7 +27,7 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+ fvc::ddt(psi)*pRef
+ fvc::ddt(psi, rho)*gh
+ fvc::div(phi)
- fvm::laplacian(rho*rUA, pd)
- fvm::laplacian(rhorUAf, pd)
);
if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
@ -39,22 +45,23 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
}
}
U += rUA*fvc::reconstruct((phi - phiU)/rhorUAf);
U.correctBoundaryConditions();
p == pd + rho*gh + pRef;
dpdt = fvc::ddt(p);
#include "rhoEqn.H"
#include "compressibleContinuityErrs.H"
U -= rUA*(fvc::grad(pd) + fvc::grad(rho)*gh);
U.correctBoundaryConditions();
// For closed-volume cases adjust the pressure and density levels
// to obey overall mass continuity
if (closedVolume)
{
p += (initialMass - fvc::domainIntegrate(thermo->psi()*p))
/fvc::domainIntegrate(thermo->psi());
pd == p - (rho*gh + pRef);
rho = thermo->rho();
}
pd == p - (rho*gh + pRef);
}

10
applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H
View File

@ -11,7 +11,15 @@
eqnResidual = solve
(
UEqn() == -fvc::grad(pd) - fvc::grad(rho)*gh
UEqn()
==
-fvc::reconstruct
(
(
fvc::snGrad(pd)
+ ghf*fvc::snGrad(rho)
) * mesh.magSf()
)
).initialResidual();
maxResidual = max(eqnResidual, maxResidual);

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

@ -53,6 +53,7 @@
Info<< "Calculating field g.h\n" << endl;
volScalarField gh("gh", g & mesh.C());
surfaceScalarField ghf("gh", g & mesh.Cf());
dimensionedScalar pRef("pRef", p.dimensions(), thermo->lookup("pRef"));

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

@ -1,19 +1,24 @@
{
volScalarField rUA = 1.0/UEqn().A();
surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA));
U = rUA*UEqn().H();
UEqn.clear();
phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
bool closedVolume = adjustPhi(phi, U, p);
phi -= fvc::interpolate(rho*gh*rUA)*fvc::snGrad(rho)*mesh.magSf();
surfaceScalarField buoyancyPhi = ghf*fvc::snGrad(rho)*rhorUAf*mesh.magSf();
phi -= buoyancyPhi;
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pdEqn
(
fvm::laplacian(rho*rUA, pd) == fvc::div(phi)
fvm::laplacian(rhorUAf, pd) == fvc::div(phi)
);
pdEqn.setReference(pdRefCell, pdRefValue);
// retain the residual from the first iteration
if (nonOrth == 0)
{
@ -27,19 +32,22 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
if (nonOrth == nNonOrthCorr)
{
// Calculate the conservative fluxes
phi -= pdEqn.flux();
// Explicitly relax pressure for momentum corrector
pd.relax();
// Correct the momentum source with the pressure gradient flux
// calculated from the relaxed pressure
U -= rUA*fvc::reconstruct((buoyancyPhi + pdEqn.flux())/rhorUAf);
U.correctBoundaryConditions();
}
}
#include "continuityErrs.H"
// Explicitly relax pressure for momentum corrector
pd.relax();
p = pd + rho*gh + pRef;
U -= rUA*(fvc::grad(pd) + fvc::grad(rho)*gh);
U.correctBoundaryConditions();
p == pd + rho*gh + pRef;
// For closed-volume cases adjust the pressure and density levels
// to obey overall mass continuity
@ -47,9 +55,11 @@ if (closedVolume)
{
p += (initialMass - fvc::domainIntegrate(thermo->psi()*p))
/fvc::domainIntegrate(thermo->psi());
pd == p - (rho*gh + pRef);
}
rho = thermo->rho();
rho.relax();
Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;
pd == p - (rho*gh + pRef);
}

1
applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/Make/options
View File

@ -1,4 +1,5 @@
EXE_INC = \
-I../buoyantSimpleFoam \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
-I$(LIB_SRC)/turbulenceModels \

18
applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/UEqn.H
View File

@ -1,18 +0,0 @@
// Solve the Momentum equation
tmp<fvVectorMatrix> UEqn
(
fvm::div(phi, U)
- fvm::Sp(fvc::div(phi), U)
+ turbulence->divDevRhoReff(U)
);
UEqn().relax();
eqnResidual = solve
(
UEqn() == -fvc::grad(pd) - fvc::grad(rho)*gh
).initialResidual();
maxResidual = max(eqnResidual, maxResidual);

1
applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/createFields.H
View File

@ -54,6 +54,7 @@
Info<< "Calculating field g.h\n" << endl;
volScalarField gh("gh", g & mesh.C());
surfaceScalarField ghf("gh", g & mesh.Cf());
dimensionedScalar pRef("pRef", p.dimensions(), thermo->lookup("pRef"));

54
applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/pEqn.H
View File

@ -1,54 +0,0 @@
volScalarField rUA = 1.0/UEqn().A();
U = rUA*UEqn().H();
UEqn.clear();
phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
bool closedVolume = adjustPhi(phi, U, p);
phi -= fvc::interpolate(rho*gh*rUA)*fvc::snGrad(rho)*mesh.magSf();
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pdEqn
(
fvm::laplacian(rho*rUA, pd) == fvc::div(phi)
);
pdEqn.setReference(pdRefCell, pdRefValue);
// retain the residual from the first iteration
if (nonOrth == 0)
{
eqnResidual = pdEqn.solve().initialResidual();
maxResidual = max(eqnResidual, maxResidual);
}
else
{
pdEqn.solve();
}
if (nonOrth == nNonOrthCorr)
{
phi -= pdEqn.flux();
}
}
#include "continuityErrs.H"
// Explicitly relax pressure for momentum corrector
pd.relax();
p = pd + rho*gh + pRef;
U -= rUA*(fvc::grad(pd) + fvc::grad(rho)*gh);
U.correctBoundaryConditions();
// For closed-volume cases adjust the pressure and density levels
// to obey overall mass continuity
if (closedVolume)
{
p += (initialMass - fvc::domainIntegrate(thermo->psi()*p))
/fvc::domainIntegrate(thermo->psi());
pd == p - (rho*gh + pRef);
}
rho = thermo->rho();
rho.relax();
Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;