zhyang-dev/OpenFOAM-12
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

buoyantPimpleFoam, chtMultiRegionFoam, rhoReactingBuoyantFoam: Share pEqn.H

Browse Source
This commit is contained in:
Will Bainbridge
2020-03-05 09:22:05 +00:00
parent f933861c48
commit 27ab3edc5e
11 changed files with 145 additions and 267 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/Make/options
View File

@ -2,6 +2,7 @@ EXE_INC = \
-I. \
-I$(FOAM_SOLVERS)/combustion/reactingFoam \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/finiteVolume/cfdTools \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \

3
applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/UEqn.H
View File

@ -1,6 +1,6 @@
MRF.correctBoundaryVelocity(U);
fvVectorMatrix UEqn
tmp<fvVectorMatrix> tUEqn
(
fvm::ddt(rho, U) + fvm::div(phi, U)
+ MRF.DDt(rho, U)
@ -8,6 +8,7 @@
==
fvOptions(rho, U)
);
fvVectorMatrix& UEqn = tUEqn.ref();
UEqn.relax();

19
applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/createFields.H
View File

@ -89,6 +89,17 @@ volScalarField p_rgh
// Force p_rgh to be consistent with p
p_rgh = p - rho*gh;
pressureControl pressureControl
(
p,
p_rgh,
rho,
pimple.dict(),
thermo.incompressible()
);
mesh.setFluxRequired(p_rgh.name());
Info<< "Creating field dpdt\n" << endl;
volScalarField dpdt
(
@ -105,9 +116,10 @@ volScalarField dpdt
Info<< "Creating field kinetic energy K\n" << endl;
volScalarField K("K", 0.5*magSqr(U));
dimensionedScalar initialMass = fvc::domainIntegrate(rho);
multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
forAll(Y, i)
{
fields.add(Y[i]);
@ -116,3 +128,8 @@ fields.add(thermo.he());
#include "createMRF.H"
#include "createFvOptions.H"
// This solver does not support moving mesh but it uses the pressure equation
// of one which does, so we need a dummy face-momentum field
autoPtr<surfaceVectorField> rhoUf(nullptr);

74
applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H
View File

@ -1,74 +0,0 @@
rho = thermo.rho();
// Thermodynamic density needs to be updated by psi*d(p) after the
// pressure solution
const volScalarField psip0(psi*p);
volScalarField rAU(1.0/UEqn.A());
surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
surfaceScalarField phiHbyA
(
"phiHbyA",
(
fvc::interpolate(rho)*fvc::flux(HbyA)
+ MRF.zeroFilter(rhorAUf*fvc::ddtCorr(rho, U, phi))
)
+ phig
);
MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
// Update the pressure BCs to ensure flux consistency
constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
fvScalarMatrix p_rghDDtEqn
(
fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
+ fvc::div(phiHbyA)
==
fvOptions(psi, p_rgh, rho.name())
);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix p_rghEqn
(
p_rghDDtEqn
- fvm::laplacian(rhorAUf, p_rgh)
);
p_rghEqn.solve();
if (pimple.finalNonOrthogonalIter())
{
// Calculate the conservative fluxes
phi = phiHbyA + p_rghEqn.flux();
// Explicitly relax pressure for momentum corrector
p_rgh.relax();
// Correct the momentum source with the pressure gradient flux
// calculated from the relaxed pressure
U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
U.correctBoundaryConditions();
fvOptions.correct(U);
K = 0.5*magSqr(U);
}
}
p = p_rgh + rho*gh;
// Thermodynamic density update
thermo.correctRho(psi*p - psip0);
if (thermo.dpdt())
{
dpdt = fvc::ddt(p);
}
#include "rhoEqn.H"
#include "compressibleContinuityErrs.H"

5
applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/rhoReactingBuoyantFoam.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2012-2019 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2012-2020 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -36,6 +36,7 @@ Description
#include "turbulentFluidThermoModel.H"
#include "multivariateScheme.H"
#include "pimpleControl.H"
#include "pressureControl.H"
#include "fvOptions.H"
#include "localEulerDdtScheme.H"
#include "fvcSmooth.H"
@ -97,7 +98,7 @@ int main(int argc, char *argv[])
// --- Pressure corrector loop
while (pimple.correct())
{
#include "pEqn.H"
#include "../../../heatTransfer/buoyantPimpleFoam/pEqn.H"
}
if (pimple.turbCorr())

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

@ -2,7 +2,7 @@
MRF.correctBoundaryVelocity(U);
fvVectorMatrix UEqn
tmp<fvVectorMatrix> tUEqn
(
fvm::ddt(rho, U) + fvm::div(phi, U)
+ MRF.DDt(rho, U)
@ -10,6 +10,7 @@
==
fvOptions(rho, U)
);
fvVectorMatrix& UEqn = tUEqn.ref();
UEqn.relax();

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

@ -102,6 +102,8 @@ volScalarField dpdt
Info<< "Creating field kinetic energy K\n" << endl;
volScalarField K("K", 0.5*magSqr(U));
dimensionedScalar initialMass = fvc::domainIntegrate(rho);
#include "createMRF.H"
#include "createRadiationModel.H"
#include "createFvOptions.H"

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

@ -1,4 +1,4 @@
if (!pimple.simpleRho())
if (!mesh.steady() && !pimple.simpleRho())
{
rho = thermo.rho();
}
@ -11,17 +11,22 @@ volScalarField rAU(1.0/UEqn.A());
surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
if (pimple.nCorrPiso() <= 1)
{
tUEqn.clear();
}
surfaceScalarField phiHbyA
(
"phiHbyA",
(
fvc::interpolate(rho)*fvc::flux(HbyA)
+ MRF.zeroFilter(rhorAUf*fvc::ddtCorr(rho, U, phi, rhoUf))
)
fvc::interpolate(rho)*fvc::flux(HbyA)
+ MRF.zeroFilter(rhorAUf*fvc::ddtCorr(rho, U, phi, rhoUf))
);
MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
bool adjustMass = false;
surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
phiHbyA += phig;
@ -31,66 +36,128 @@ constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
fvc::makeRelative(phiHbyA, rho, U);
fvScalarMatrix p_rghDDtEqn
(
fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
+ fvc::div(phiHbyA)
==
fvOptions(psi, p_rgh, rho.name())
);
fvScalarMatrix p_rghEqn(p_rgh, dimMass/dimTime);
while (pimple.correctNonOrthogonal())
if (pimple.transonic())
{
fvScalarMatrix p_rghEqn(p_rghDDtEqn - fvm::laplacian(rhorAUf, p_rgh));
p_rghEqn.setReference
surfaceScalarField phid
(
pressureControl.refCell(),
pressureControl.refValue()
"phid",
(fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
);
p_rghEqn.solve();
phiHbyA -= fvc::interpolate(psi*p_rgh)*phiHbyA/fvc::interpolate(rho);
if (pimple.finalNonOrthogonalIter())
fvScalarMatrix p_rghDDtEqn
(
fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
+ fvc::div(phiHbyA) + fvm::div(phid, p_rgh)
==
fvOptions(psi, p_rgh, rho.name())
);
while (pimple.correctNonOrthogonal())
{
// Calculate the conservative fluxes
phi = phiHbyA + p_rghEqn.flux();
p_rghEqn = p_rghDDtEqn - fvm::laplacian(rhorAUf, p_rgh);
// Explicitly relax pressure for momentum corrector
p_rgh.relax();
// Relax the pressure equation to ensure diagonal-dominance
p_rghEqn.relax();
// Correct the momentum source with the pressure gradient flux
// calculated from the relaxed pressure
U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
U.correctBoundaryConditions();
fvOptions.correct(U);
K = 0.5*magSqr(U);
p_rghEqn.setReference
(
pressureControl.refCell(),
pressureControl.refValue()
);
p_rghEqn.solve();
}
}
else
{
if (mesh.steady())
{
adjustMass = adjustPhi(phiHbyA, U, p_rgh);
}
fvScalarMatrix p_rghDDtEqn
(
fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
+ fvc::div(phiHbyA)
==
fvOptions(psi, p_rgh, rho.name())
);
while (pimple.correctNonOrthogonal())
{
p_rghEqn = p_rghDDtEqn - fvm::laplacian(rhorAUf, p_rgh);
p_rghEqn.setReference
(
pressureControl.refCell(),
pressureControl.refValue()
);
p_rghEqn.solve();
}
}
phi = phiHbyA + p_rghEqn.flux();
p = p_rgh + rho*gh;
bool limitedp = pressureControl.limit(p);
if (limitedp)
if (mesh.steady())
{
#include "incompressible/continuityErrs.H"
}
else
{
const bool limitedp = pressureControl.limit(p);
// Thermodynamic density update
thermo.correctRho(psi*p - psip0);
if (limitedp)
{
rho = thermo.rho();
}
#include "rhoEqn.H"
#include "compressibleContinuityErrs.H"
}
// Explicitly relax pressure for momentum corrector
p_rgh.relax();
// Correct the momentum source with the pressure gradient flux
// calculated from the relaxed pressure
U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
U.correctBoundaryConditions();
fvOptions.correct(U);
K = 0.5*magSqr(U);
if (mesh.steady())
{
pressureControl.limit(p);
}
// For steady closed-volume compressible cases adjust the pressure level
// to obey overall mass continuity
if (adjustMass && !thermo.incompressible())
{
p += (initialMass - fvc::domainIntegrate(thermo.rho()))
/fvc::domainIntegrate(psi);
p_rgh = p - rho*gh;
p.correctBoundaryConditions();
}
// Thermodynamic density update
thermo.correctRho(psi*p - psip0);
if (limitedp)
if (mesh.steady() || pimple.simpleRho() || adjustMass)
{
rho = thermo.rho();
}
#include "rhoEqn.H"
#include "compressibleContinuityErrs.H"
if (pimple.simpleRho())
if ((mesh.steady() || pimple.simpleRho()) && !pimple.transonic())
{
rho = thermo.rho();
rho.relax();
}
// Correct rhoUf if the mesh is moving
@ -105,3 +172,6 @@ if (thermo.dpdt())
dpdt -= fvc::div(fvc::meshPhi(rho, U), p);
}
}
Info<< "Min/max rho:" << min(rho).value() << ' '
<< max(rho).value() << endl;

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

@ -1,145 +0,0 @@
if (!mesh.steady() && !pimple.simpleRho())
{
rho = thermo.rho();
}
volScalarField rAU("rAU", 1.0/UEqn.A());
surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
if (pimple.nCorrPiso() <= 1)
{
tUEqn.clear();
}
surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
surfaceScalarField phiHbyA
(
"phiHbyA",
fvc::interpolate(rho)*fvc::flux(HbyA)
+ MRF.zeroFilter(rhorAUf*fvc::ddtCorr(rho, U, phi))
);
MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
const bool closedVolume = mesh.steady() && adjustPhi(phiHbyA, U, p_rgh);
const bool adjustMass = closedVolume && !thermo.incompressible();
phiHbyA += phig;
// Update the pressure BCs to ensure flux consistency
constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
{
fvScalarMatrix p_rghEqnComp
(
fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
==
fvOptions(psi, p_rgh, rho.name())
);
if (pimple.transonic())
{
surfaceScalarField phid
(
"phid",
(fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
);
phiHbyA -= fvc::interpolate(psi*p_rgh)*phiHbyA/fvc::interpolate(rho);
p_rghEqnComp += fvm::div(phid, p_rgh);
}
// Thermodynamic density needs to be updated by psi*d(p) after the
// pressure solution
tmp<volScalarField> psip0(mesh.steady() ? tmp<volScalarField>() : psi*p);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix p_rghEqnIncomp
(
fvc::div(phiHbyA)
- fvm::laplacian(rhorAUf, p_rgh)
);
fvScalarMatrix p_rghEqn(p_rghEqnComp + p_rghEqnIncomp);
p_rghEqn.setReference
(
pressureControl.refCell(),
pressureControl.refValue()
);
p_rghEqn.solve();
if (pimple.finalNonOrthogonalIter())
{
// Calculate the conservative fluxes
phi = phiHbyA + p_rghEqn.flux();
// Explicitly relax pressure for momentum corrector
p_rgh.relax();
// Correct the momentum source with the pressure gradient flux
// calculated from the relaxed pressure
U = HbyA
+ rAU*fvc::reconstruct((phig + p_rghEqnIncomp.flux())/rhorAUf);
U.correctBoundaryConditions();
fvOptions.correct(U);
K = 0.5*magSqr(U);
}
}
p = p_rgh + rho*gh;
pressureControl.limit(p);
// Thermodynamic density update
if (!mesh.steady())
{
thermo.correctRho(psi*p - psip0);
}
}
// Update pressure time derivative if needed
if (thermo.dpdt())
{
dpdt = fvc::ddt(p);
}
// Solve continuity
if (!mesh.steady())
{
#include "rhoEqn.H"
#include "compressibleContinuityErrs.H"
}
else
{
#include "incompressible/continuityErrs.H"
}
// For closed-volume compressible cases adjust the pressure level
// to obey overall mass continuity
if (adjustMass)
{
p += (initialMass - fvc::domainIntegrate(thermo.rho()))
/fvc::domainIntegrate(psi);
p_rgh = p - rho*gh;
p.correctBoundaryConditions();
}
// Density updates
if (adjustMass || mesh.steady() || pimple.simpleRho())
{
rho = thermo.rho();
}
if (mesh.steady() && !pimple.transonic())
{
rho.relax();
}
Info<< "Min/max rho:" << min(rho).value() << ' '
<< max(rho).value() << endl;

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

@ -60,3 +60,7 @@
pressureControl& pressureControl = pressureControlFluid[i];
scalar cumulativeContErr = cumulativeContErrs[i];
// This solver does not support moving mesh but it uses the pressure
// equation of one which does, so we need a dummy face-momentum field
autoPtr<surfaceVectorField> rhoUf(nullptr);

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

@ -16,7 +16,7 @@ else
// --- PISO loop
while (pimple.correct())
{
#include "pEqn.H"
#include "../../buoyantPimpleFoam/pEqn.H"
}
if (pimples.pimpleTurbCorr(i))