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

GIT: Resolve conflict associated with cherry-pick of Foundation commit 79ff91350

Browse Source 
79ff91350 - rhoPimpleFoam: Improved support for compressible liquids
(2017-05-17 17:05:43 +0100) <Henry Weller>
This commit is contained in:
Henry Weller
2017-05-17 17:05:43 +01:00
committed by Andrew Heather
parent b83af3b085
commit 39476bde1c
44 changed files with 1194 additions and 607 deletions
Download Patch File Download Diff File Expand all files Collapse all files

22
applications/solvers/combustion/reactingFoam/createFields.H
View File

@ -52,27 +52,7 @@ volScalarField& p = thermo.p();
#include "compressibleCreatePhi.H"
dimensionedScalar rhoMax
(
dimensionedScalar::lookupOrDefault
(
"rhoMax",
pimple.dict(),
dimDensity,
GREAT
)
);
dimensionedScalar rhoMin
(
dimensionedScalar::lookupOrDefault
(
"rhoMin",
pimple.dict(),
dimDensity,
0
)
);
pressureControl pressureControl(p, rho, pimple.dict(), false);
mesh.setFluxRequired(p.name());

17
applications/solvers/combustion/reactingFoam/pEqn.H
View File

@ -1,7 +1,4 @@
rho = thermo.rho();
rho = max(rho, rhoMin);
rho = min(rho, rhoMax);
rho.relax();
volScalarField rAU(1.0/UEqn.A());
surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
@ -87,19 +84,17 @@ else
// Explicitly relax pressure for momentum corrector
p.relax();
// Recalculate density from the relaxed pressure
rho = thermo.rho();
rho = max(rho, rhoMin);
rho = min(rho, rhoMax);
rho.relax();
Info<< "rho max/min : " << max(rho).value()
<< " " << min(rho).value() << endl;
U = HbyA - rAU*fvc::grad(p);
U.correctBoundaryConditions();
fvOptions.correct(U);
K = 0.5*magSqr(U);
if (pressureControl.limit(p))
{
p.correctBoundaryConditions();
rho = thermo.rho();
}
if (thermo.dpdt())
{
dpdt = fvc::ddt(p);

21
applications/solvers/combustion/reactingFoam/pcEqn.H
View File

@ -1,7 +1,4 @@
rho = thermo.rho();
rho = max(rho, rhoMin);
rho = min(rho, rhoMax);
rho.relax();
volScalarField rAU(1.0/UEqn.A());
volScalarField rAtU(1.0/(1.0/rAU - UEqn.H1()));
@ -109,19 +106,13 @@ U.correctBoundaryConditions();
fvOptions.correct(U);
K = 0.5*magSqr(U);
if (pressureControl.limit(p))
{
p.correctBoundaryConditions();
rho = thermo.rho();
}
if (thermo.dpdt())
{
dpdt = fvc::ddt(p);
}
// Recalculate density from the relaxed pressure
rho = thermo.rho();
rho = max(rho, rhoMin);
rho = min(rho, rhoMax);
if (!pimple.transonic())
{
rho.relax();
}
Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;

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

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -37,6 +37,7 @@ Description
#include "psiCombustionModel.H"
#include "multivariateScheme.H"
#include "pimpleControl.H"
#include "pressureControl.H"
#include "fvOptions.H"
#include "localEulerDdtScheme.H"
#include "fvcSmooth.H"
@ -114,6 +115,8 @@ int main(int argc, char *argv[])
}
}
rho = thermo.rho();
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

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

@ -1,76 +1,74 @@
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::flux(rho*HbyA)
+ 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())
{
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;
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
fvScalarMatrix p_rghEqn
(
"phiHbyA",
(
fvc::flux(rho*HbyA)
+ rhorAUf*fvc::ddtCorr(rho, U, phi)
)
+ phig
p_rghDDtEqn
- fvm::laplacian(rhorAUf, p_rgh)
);
MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
// 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())
if (pimple.finalNonOrthogonalIter())
{
fvScalarMatrix p_rghEqn
(
p_rghDDtEqn
- fvm::laplacian(rhorAUf, p_rgh)
);
// Calculate the conservative fluxes
phi = phiHbyA + p_rghEqn.flux();
p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
// Explicitly relax pressure for momentum corrector
p_rgh.relax();
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);
}
// 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;
// Second part of thermodynamic density update
thermo.rho() += psi*p;
if (thermo.dpdt())
{
dpdt = fvc::ddt(p);
}
#include "rhoEqn.H"
#include "compressibleContinuityErrs.H"
}
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"

2
applications/solvers/combustion/reactingFoam/rhoReactingFoam/createFields.H
View File

@ -52,6 +52,8 @@ volScalarField& p = thermo.p();
#include "compressibleCreatePhi.H"
pressureControl pressureControl(p, rho, pimple.dict(), false);
mesh.setFluxRequired(p.name());

109
applications/solvers/combustion/reactingFoam/rhoReactingFoam/pEqn.H
View File

@ -1,109 +0,0 @@
{
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;
volScalarField rAU(1.0/UEqn.A());
surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
if (pimple.transonic())
{
surfaceScalarField phiHbyA
(
"phiHbyA",
(
fvc::flux(HbyA)
+ rhorAUf*fvc::ddtCorr(rho, U, phi)/fvc::interpolate(rho)
)
);
MRF.makeRelative(phiHbyA);
surfaceScalarField phid("phid", fvc::interpolate(thermo.psi())*phiHbyA);
phiHbyA *= fvc::interpolate(rho);
fvScalarMatrix pDDtEqn
(
fvc::ddt(rho) + fvc::div(phiHbyA)
+ correction(psi*fvm::ddt(p) + fvm::div(phid, p))
);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
pDDtEqn
- fvm::laplacian(rhorAUf, p)
==
fvOptions(psi, p, rho.name())
);
pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
if (pimple.finalNonOrthogonalIter())
{
phi = phiHbyA + pEqn.flux();
}
}
}
else
{
surfaceScalarField phiHbyA
(
"phiHbyA",
(
fvc::flux(rho*HbyA)
+ rhorAUf*fvc::ddtCorr(rho, U, phi)
)
);
MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
// Update the pressure BCs to ensure flux consistency
constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
fvScalarMatrix pDDtEqn
(
fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+ fvc::div(phiHbyA)
==
fvOptions(psi, p, rho.name())
);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
pDDtEqn
- fvm::laplacian(rhorAUf, p)
);
pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
if (pimple.finalNonOrthogonalIter())
{
phi = phiHbyA + pEqn.flux();
}
}
}
// Second part of thermodynamic density update
thermo.rho() += psi*p;
#include "rhoEqn.H"
#include "compressibleContinuityErrs.H"
U = HbyA - rAU*fvc::grad(p);
U.correctBoundaryConditions();
fvOptions.correct(U);
K = 0.5*magSqr(U);
if (thermo.dpdt())
{
dpdt = fvc::ddt(p);
}
}

12
applications/solvers/combustion/reactingFoam/rhoReactingFoam/rhoReactingFoam.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -38,6 +38,7 @@ Description
#include "turbulentFluidThermoModel.H"
#include "multivariateScheme.H"
#include "pimpleControl.H"
#include "pressureControl.H"
#include "fvOptions.H"
#include "localEulerDdtScheme.H"
#include "fvcSmooth.H"
@ -100,7 +101,14 @@ int main(int argc, char *argv[])
// --- Pressure corrector loop
while (pimple.correct())
{
#include "pEqn.H"
if (pimple.consistent())
{
#include "../../../compressible/rhoPimpleFoam/pcEqn.H"
}
else
{
#include "../../../compressible/rhoPimpleFoam/pEqn.H"
}
}
if (pimple.turbCorr())

4
applications/solvers/compressible/rhoPimpleFoam/createFields.H
View File

@ -42,6 +42,8 @@ volVectorField U
pressureControl pressureControl(p, rho, pimple.dict(), false);
mesh.setFluxRequired(p.name());
Info<< "Creating turbulence model\n" << endl;
autoPtr<compressible::turbulenceModel> turbulence
(
@ -54,8 +56,6 @@ autoPtr<compressible::turbulenceModel> turbulence
)
);
mesh.setFluxRequired(p.name());
Info<< "Creating field dpdt\n" << endl;
volScalarField dpdt
(

52
applications/solvers/compressible/rhoPimpleFoam/pEqn.H
View File

@ -1,4 +1,11 @@
rho = thermo.rho();
if (!pimple.SIMPLErho())
{
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));
@ -31,17 +38,17 @@ if (pimple.transonic())
phiHbyA -= fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
fvScalarMatrix pDDtEqn
(
fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+ fvc::div(phiHbyA) + fvm::div(phid, p)
==
fvOptions(psi, p, rho.name())
);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
fvm::ddt(psi, p)
+ fvc::div(phiHbyA)
+ fvm::div(phid, p)
- fvm::laplacian(rhorAUf, p)
==
fvOptions(psi, p, rho.name())
);
fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAUf, p));
// Relax the pressure equation to ensure diagonal-dominance
pEqn.relax();
@ -56,16 +63,17 @@ if (pimple.transonic())
}
else
{
fvScalarMatrix pDDtEqn
(
fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+ fvc::div(phiHbyA)
==
fvOptions(psi, p, rho.name())
);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
fvm::ddt(psi, p)
+ fvc::div(phiHbyA)
- fvm::laplacian(rhorAUf, p)
==
fvOptions(psi, p, rho.name())
);
fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAUf, p));
pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
@ -76,6 +84,9 @@ else
}
}
// Thermodynamic density update
thermo.correctRho(psi*p - psip0);
#include "rhoEqn.H"
#include "compressibleContinuityErrs.H"
@ -92,10 +103,9 @@ if (pressureControl.limit(p))
p.correctBoundaryConditions();
rho = thermo.rho();
}
if (!pimple.transonic())
else if (pimple.SIMPLErho())
{
rho.relax();
rho = thermo.rho();
}
if (thermo.dpdt())

52
applications/solvers/compressible/rhoPimpleFoam/pcEqn.H
View File

@ -1,4 +1,11 @@
rho = thermo.rho();
if (!pimple.SIMPLErho())
{
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());
volScalarField rAtU(1.0/(1.0/rAU - UEqn.H1()));
@ -39,17 +46,17 @@ if (pimple.transonic())
HbyA -= (rAU - rAtU)*fvc::grad(p);
fvScalarMatrix pDDtEqn
(
fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+ fvc::div(phiHbyA) + fvm::div(phid, p)
==
fvOptions(psi, p, rho.name())
);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
fvm::ddt(psi, p)
+ fvc::div(phiHbyA)
+ fvm::div(phid, p)
- fvm::laplacian(rhorAtU, p)
==
fvOptions(psi, p, rho.name())
);
fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAtU, p));
// Relax the pressure equation to ensure diagonal-dominance
pEqn.relax();
@ -67,16 +74,17 @@ else
phiHbyA += fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf();
HbyA -= (rAU - rAtU)*fvc::grad(p);
fvScalarMatrix pDDtEqn
(
fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+ fvc::div(phiHbyA)
==
fvOptions(psi, p, rho.name())
);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
fvm::ddt(psi, p)
+ fvc::div(phiHbyA)
- fvm::laplacian(rhorAtU, p)
==
fvOptions(psi, p, rho.name())
);
fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAtU, p));
pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
@ -87,6 +95,9 @@ else
}
}
// Thermodynamic density update
thermo.correctRho(psi*p - psip0);
#include "rhoEqn.H"
#include "compressibleContinuityErrs.H"
@ -103,10 +114,9 @@ if (pressureControl.limit(p))
p.correctBoundaryConditions();
rho = thermo.rho();
}
if (!pimple.transonic())
else if (pimple.SIMPLErho())
{
rho.relax();
rho = thermo.rho();
}
if (thermo.dpdt())

9
applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
View File

@ -36,7 +36,7 @@ if (pimple.transonic())
{
fvScalarMatrix pEqn
(
fvm::ddt(psi, p)
fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+ fvc::div(phiHbyA)
+ fvm::div(phid, p)
- fvm::laplacian(rhorAUf, p)
@ -62,7 +62,7 @@ else
// Pressure corrector
fvScalarMatrix pEqn
(
fvm::ddt(psi, p)
fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+ fvc::div(phiHbyA)
- fvm::laplacian(rhorAUf, p)
==
@ -95,11 +95,6 @@ if (pressureControl.limit(p))
rho = thermo.rho();
}
if (!pimple.transonic())
{
rho.relax();
}
{
rhoUf = fvc::interpolate(rho*U);
surfaceVectorField n(mesh.Sf()/mesh.magSf());

2
applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
View File

@ -122,6 +122,8 @@ int main(int argc, char *argv[])
}
}
rho = thermo.rho();
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

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

@ -1,115 +1,74 @@
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_rgh));
surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
surfaceScalarField phiHbyA
(
"phiHbyA",
(
fvc::flux(rho*HbyA)
+ 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())
{
bool closedVolume = p_rgh.needReference();
dimensionedScalar compressibility = fvc::domainIntegrate(psi);
bool compressible = (compressibility.value() > SMALL);
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_rgh;
volScalarField rAU(1.0/UEqn.A());
surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
surfaceScalarField phiHbyA
fvScalarMatrix p_rghEqn
(
"phiHbyA",
(
fvc::flux(rho*HbyA)
+ rhorAUf*fvc::ddtCorr(rho, U, phi)
)
+ phig
p_rghDDtEqn
- fvm::laplacian(rhorAUf, p_rgh)
);
MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
// Update the pressure BCs to ensure flux consistency
constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
tmp<fvScalarMatrix> p_rghDDtEqn
(
new fvScalarMatrix(p_rgh, dimMass/dimTime)
);
if (compressible)
if (pimple.finalNonOrthogonalIter())
{
p_rghDDtEqn =
(
fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
);
// 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);
}
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix p_rghEqn
(
p_rghDDtEqn()
+ fvc::div(phiHbyA)
- fvm::laplacian(rhorAUf, p_rgh)
==
fvOptions(psi, p_rgh, rho.name())
);
p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
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;
// Second part of thermodynamic density update
thermo.rho() += psi*p_rgh;
if (thermo.dpdt())
{
dpdt = fvc::ddt(p);
}
if (compressible)
{
#include "rhoEqn.H"
}
#include "compressibleContinuityErrs.H"
if (closedVolume)
{
if (!compressible)
{
p += dimensionedScalar
(
"p",
p.dimensions(),
pRefValue - getRefCellValue(p, pRefCell)
);
}
else
{
p += (initialMass - fvc::domainIntegrate(thermo.rho()))
/compressibility;
rho = thermo.rho();
}
p_rgh = p - rho*gh;
}
Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value()
<< endl;
}
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"

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

@ -1,114 +1,100 @@
bool closedVolume = p_rgh.needReference();
dimensionedScalar compressibility = fvc::domainIntegrate(psi);
bool compressible = (compressibility.value() > SMALL);
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));
surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
surfaceScalarField phiHbyA
(
"phiHbyA",
(
fvc::flux(rho*HbyA)
+ 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);
{
bool closedVolume = p_rgh.needReference();
dimensionedScalar compressibility = fvc::domainIntegrate(psi);
bool compressible = (compressibility.value() > SMALL);
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));
surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
surfaceScalarField phiHbyA
fvScalarMatrix p_rghDDtEqn
(
"phiHbyA",
(
fvc::flux(rho*HbyA)
+ rhorAUf*fvc::ddtCorr(rho, U, phi)
)
+ phig
fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
+ fvc::div(phiHbyA)
);
MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
// Update the pressure BCs to ensure flux consistency
constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
tmp<fvScalarMatrix> p_rghDDtEqn
(
new fvScalarMatrix(p_rgh, dimMass/dimTime)
);
// Thermodynamic density needs to be updated by psi*d(p) after the
// pressure solution
const volScalarField psip0(psi*p);
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
if (compressible)
{
p_rghDDtEqn =
(
fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
==
fvOptions(psi, p_rgh, rho.name())
);
}
fvScalarMatrix p_rghEqn
(
p_rghDDtEqn
- fvm::laplacian(rhorAUf, p_rgh)
);
// Thermodynamic density needs to be updated by psi*d(p) after the
// pressure solution - done in 2 parts. Part 1:
thermo.rho() -= psi*p_rgh;
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix p_rghEqn
solvPerfp_rgh = p_rghEqn.solve
(
mesh.solver
(
p_rghDDtEqn()
+ fvc::div(phiHbyA)
- fvm::laplacian(rhorAUf, p_rgh)
);
p_rghEqn.solve
(
mesh.solver
p_rgh.select
(
p_rgh.select
(
(
oCorr == nOuterCorr-1
&& corr == nCorr-1
&& nonOrth == nNonOrthCorr
)
oCorr == nOuterCorr-1
&& corr == nCorr-1
&& nonOrth == nNonOrthCorr
)
)
);
)
);
if (nonOrth == nNonOrthCorr)
{
phi = phiHbyA + p_rghEqn.flux();
U = HbyA
+ rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
U.correctBoundaryConditions();
fvOptions.correct(U);
K = 0.5*magSqr(U);
}
if (nonOrth == nNonOrthCorr)
{
phi = phiHbyA + p_rghEqn.flux();
U = HbyA
+ rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
U.correctBoundaryConditions();
fvOptions.correct(U);
K = 0.5*magSqr(U);
}
// Second part of thermodynamic density update
thermo.rho() += psi*p_rgh;
}
p = p_rgh + rho*gh;
// Update pressure time derivative if needed
if (thermo.dpdt())
{
dpdt = fvc::ddt(p);
}
if (compressible)
{
// Solve continuity
#include "rhoEqn.H"
}
// Update continuity errors
#include "compressibleContinuityErrors.H"
// For closed-volume cases adjust the pressure and density levels
// to obey overall mass continuity
if (closedVolume && compressible)
{
p += (initialMass - fvc::domainIntegrate(thermo.rho()))
/compressibility;
rho = thermo.rho();
p_rgh = p - rho*gh;
}
// Thermodynamic density update
thermo.correctRho(psi*p - psip0);
}
// Update pressure time derivative if needed
if (thermo.dpdt())
{
dpdt = fvc::ddt(p);
}
// Solve continuity
#include "rhoEqn.H"
// Update continuity errors
#include "compressibleContinuityErrors.H"
// For closed-volume cases adjust the pressure and density levels
// to obey overall mass continuity
if (closedVolume && compressible)
{
p += (initialMass - fvc::domainIntegrate(thermo.rho()))
/compressibility;
rho = thermo.rho();
p_rgh = p - rho*gh;
}

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

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -42,6 +42,7 @@ Description
#include "radiationModel.H"
#include "SLGThermo.H"
#include "pimpleControl.H"
#include "pressureControl.H"
#include "localEulerDdtScheme.H"
#include "fvcSmooth.H"

128
applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
View File

@ -1,73 +1,71 @@
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 phiHbyA
(
"phiHbyA",
(
fvc::flux(rho*HbyA)
+ rhorAUf*fvc::ddtCorr(rho, U, phi)
)
);
MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
// Update the pressure BCs to ensure flux consistency
constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
fvScalarMatrix pDDtEqn
(
fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+ fvc::div(phiHbyA)
==
parcels.Srho()
+ fvOptions(psi, p, rho.name())
);
while (pimple.correctNonOrthogonal())
{
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;
volScalarField rAU(1.0/UEqn.A());
surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
surfaceScalarField phiHbyA
fvScalarMatrix pEqn
(
"phiHbyA",
(
fvc::flux(rho*HbyA)
+ rhorAUf*fvc::ddtCorr(rho, U, phi)
)
pDDtEqn
- fvm::laplacian(rhorAUf, p)
);
MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
// Update the pressure BCs to ensure flux consistency
constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
fvScalarMatrix pDDtEqn
(
fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+ fvc::div(phiHbyA)
==
parcels.Srho()
+ fvOptions(psi, p, rho.name())
);
while (pimple.correctNonOrthogonal())
if (pimple.finalNonOrthogonalIter())
{
fvScalarMatrix pEqn
(
pDDtEqn
- fvm::laplacian(rhorAUf, p)
);
pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
if (pimple.finalNonOrthogonalIter())
{
phi = phiHbyA + pEqn.flux();
}
phi = phiHbyA + pEqn.flux();
}
p.relax();
// Second part of thermodynamic density update
thermo.rho() += psi*p;
#include "rhoEqn.H" // NOTE: flux and time scales now inconsistent
#include "compressibleContinuityErrs.H"
U = HbyA - rAU*fvc::grad(p);
U.correctBoundaryConditions();
fvOptions.correct(U);
K = 0.5*magSqr(U);
if (thermo.dpdt())
{
dpdt = fvc::ddt(p);
}
rho = thermo.rho();
rho = max(rho, rhoMin);
rho = min(rho, rhoMax);
Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl;
}
p.relax();
// Thermodynamic density update
thermo.correctRho(psi*p - psip0);
#include "rhoEqn.H" // NOTE: flux and time scales now inconsistent
#include "compressibleContinuityErrs.H"
U = HbyA - rAU*fvc::grad(p);
U.correctBoundaryConditions();
fvOptions.correct(U);
K = 0.5*magSqr(U);
if (thermo.dpdt())
{
dpdt = fvc::ddt(p);
}
rho = thermo.rho();
rho = max(rho, rhoMin);
rho = min(rho, rhoMax);
Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl;

96
applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/pEqn.H
View File

@ -1,57 +1,55 @@
{
// Thermodynamic density needs to be updated by psi*d(p) after the
// pressure solution - done in 2 parts. Part 1:
thermo.rho() -= psi*p;
// 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));
tUEqn.clear();
surfaceScalarField phiHbyA
volScalarField rAU(1.0/UEqn.A());
surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
tUEqn.clear();
surfaceScalarField phiHbyA
(
"phiHbyA",
fvc::interpolate(rho)*fvc::flux(HbyA)
);
MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
// Update the pressure BCs to ensure flux consistency
constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
while (simple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
"phiHbyA",
fvc::interpolate(rho)*fvc::flux(HbyA)
fvc::div(phiHbyA)
- fvm::laplacian(rhorAUf, p)
==
parcels.Srho()
+ fvOptions(psi, p, rho.name())
);
MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
pEqn.solve();
// Update the pressure BCs to ensure flux consistency
constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
while (simple.correctNonOrthogonal())
if (simple.finalNonOrthogonalIter())
{
fvScalarMatrix pEqn
(
fvc::div(phiHbyA)
- fvm::laplacian(rhorAUf, p)
==
parcels.Srho()
+ fvOptions(psi, p, rho.name())
);
pEqn.solve();
if (simple.finalNonOrthogonalIter())
{
phi = phiHbyA + pEqn.flux();
}
phi = phiHbyA + pEqn.flux();
}
p.relax();
// Second part of thermodynamic density update
thermo.rho() += psi*p;
#include "compressibleContinuityErrs.H"
U = HbyA - rAU*fvc::grad(p);
U.correctBoundaryConditions();
fvOptions.correct(U);
rho = thermo.rho();
rho = max(rho, rhoMin);
rho = min(rho, rhoMax);
rho.relax();
Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl;
}
p.relax();
// Thermodynamic density update
thermo.correctRho(psi*p - psip0);
#include "compressibleContinuityErrs.H"
U = HbyA - rAU*fvc::grad(p);
U.correctBoundaryConditions();
fvOptions.correct(U);
rho = thermo.rho();
rho = max(rho, rhoMin);
rho = min(rho, rhoMax);
rho.relax();
Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl;

4
applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
View File

@ -104,8 +104,8 @@
}
// Update densities from change in p_rgh
rho1 += psi1*(p_rgh - p_rgh_0);
rho2 += psi2*(p_rgh - p_rgh_0);
mixture.thermo1().correctRho(psi1*(p_rgh - p_rgh_0));
mixture.thermo2().correctRho(psi2*(p_rgh - p_rgh_0));
rho = alpha1*rho1 + alpha2*rho2;

4
applications/solvers/multiphase/compressibleInterFoam/createFields.H
View File

@ -38,8 +38,8 @@ volScalarField& alpha2(mixture.alpha2());
Info<< "Reading thermophysical properties\n" << endl;
volScalarField& rho1 = mixture.thermo1().rho();
volScalarField& rho2 = mixture.thermo2().rho();
const volScalarField& rho1 = mixture.thermo1().rho();
const volScalarField& rho2 = mixture.thermo2().rho();
volScalarField rho
(

4
applications/solvers/multiphase/compressibleInterFoam/pEqn.H
View File

@ -104,8 +104,8 @@
}
// Update densities from change in p_rgh
rho1 += psi1*(p_rgh - p_rgh_0);
rho2 += psi2*(p_rgh - p_rgh_0);
mixture.thermo1().correctRho(psi1*(p_rgh - p_rgh_0));
mixture.thermo2().correctRho(psi2*(p_rgh - p_rgh_0));
rho = alpha1*rho1 + alpha2*rho2;