diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/UEqn.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/UEqn.H
index 2e8f979be419dba827b3cf63867d13db67c5aa92..61d805e01dffa692d0292b556d1da74f97c19d68 100644
--- a/applications/solvers/lagrangian/reactingParcelFilmFoam/UEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/UEqn.H
@@ -4,13 +4,23 @@
+ fvm::div(phi, U)
+ turbulence->divDevRhoReff(U)
==
- rho.dimensionedInternalField()*g
- + parcels.SU(U)
+ parcels.SU(U)
);
UEqn.relax();
if (momentumPredictor)
{
- solve(UEqn == -fvc::grad(p));
+ solve
+ (
+ UEqn
+ ==
+ fvc::reconstruct
+ (
+ (
+ - ghf*fvc::snGrad(rho)
+ - fvc::snGrad(p_rgh)
+ )*mesh.magSf()
+ )
+ );
}
diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/YEqn.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/YEqn.H
index 8526b590a297c6f8fd351af4deb63cdebd713f0d..5eed16aef407f694fd67b644eabba11a14c8a005 100644
--- a/applications/solvers/lagrangian/reactingParcelFilmFoam/YEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/YEqn.H
@@ -5,7 +5,7 @@ tmp<fv::convectionScheme<scalar> > mvConvection
mesh,
fields,
phi,
- mesh.divScheme("div(phi,Yi_h)")
+ mesh.divScheme("div(phi,Yi_hs)")
)
);
diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/chemistry.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/chemistry.H
index 5dd7ce9cb03c26f4bd5219bb3b7acf6ae72a925f..ce36fa0d1b5b6b472958ed9c38a0010c4f710888 100644
--- a/applications/solvers/lagrangian/reactingParcelFilmFoam/chemistry.H
+++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/chemistry.H
@@ -1,3 +1,4 @@
+if (chemistry.chemistry())
{
Info << "Solving chemistry" << endl;
diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H
index 2e3208b0ee99d2d65ca78ca9f4ec56bef8737dd7..a5a87c73efdc4b747e903373126a929bcee6901f 100644
--- a/applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H
+++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H
@@ -15,11 +15,6 @@
const word inertSpecie(thermo.lookup("inertSpecie"));
- volScalarField& p = thermo.p();
- volScalarField& hs = thermo.hs();
- const volScalarField& T = thermo.T();
- const volScalarField& psi = thermo.psi();
-
Info<< "Creating field rho\n" << endl;
volScalarField rho
(
@@ -34,6 +29,11 @@
thermo.rho()
);
+ volScalarField& p = thermo.p();
+ volScalarField& hs = thermo.hs();
+ const volScalarField& T = thermo.T();
+ const volScalarField& psi = thermo.psi();
+
Info<< "\nReading field U\n" << endl;
volVectorField U
(
@@ -84,6 +84,28 @@
fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p)
);
+
+ Info<< "Calculating field g.h\n" << endl;
+ volScalarField gh("gh", g & mesh.C());
+
+ surfaceScalarField ghf("gh", g & mesh.Cf());
+
+ 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;
+
multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
forAll(Y, i)
diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/hsEqn.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/hsEqn.H
index 0cb2318252b46fbb7734e0fd5b46b3abd9dbb7a3..feb112f652a3f9478af26963b590aa65e2e47de4 100644
--- a/applications/solvers/lagrangian/reactingParcelFilmFoam/hsEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/hsEqn.H
@@ -19,4 +19,6 @@
thermo.correct();
radiation->correct();
+
+ Info<< "min/max(T) = " << min(T).value() << ", " << max(T).value() << endl;
}
diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H
index 9d028d4502add888fb67dea12d4766c539e62e18..ff7c7169680e21fe9f46b67901dda27b5f143128 100644
--- a/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H
@@ -1,74 +1,58 @@
rho = thermo.rho();
volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField rhorAUf(rAU.name(), fvc::interpolate(rho*rAU));
U = rAU*UEqn.H();
-if (transonic)
+surfaceScalarField phiU
+(
+ fvc::interpolate(rho)
+ *(
+ (fvc::interpolate(U) & mesh.Sf())
+ + fvc::ddtPhiCorr(rAU, rho, U, phi)
+ )
+);
+
+phi = phiU - rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf();
+
+for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
- surfaceScalarField phid
+ fvScalarMatrix p_rghEqn
(
- "phid",
- fvc::interpolate(psi)
- *(
- (fvc::interpolate(U) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, rho, U, phi)
- )
+ fvc::ddt(psi, rho)*gh
+ + fvc::div(phi)
+ + fvm::ddt(psi, p_rgh)
+ - fvm::laplacian(rhorAUf, p_rgh)
+ ==
+ parcels.Srho()
+ + surfaceFilm.Srho()
);
- for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
- {
- fvScalarMatrix pEqn
+ p_rghEqn.solve
+ (
+ mesh.solver
(
- fvm::ddt(psi, p)
- + fvm::div(phid, p)
- - fvm::laplacian(rho*rAU, p)
- ==
- parcels.Srho()
- + surfaceFilm.Srho()
- );
-
- pEqn.solve();
-
- if (nonOrth == nNonOrthCorr)
- {
- phi == pEqn.flux();
- }
- }
-}
-else
-{
- phi =
- fvc::interpolate(rho)
- *(
- (fvc::interpolate(U) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, rho, U, phi)
- );
+ p_rgh.select
+ (
+ pimpleCorr.finalIter()
+ && corr == nCorr-1
+ && nonOrth == nNonOrthCorr
+ )
+ )
+ );
- for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+ if (nonOrth == nNonOrthCorr)
{
- fvScalarMatrix pEqn
- (
- fvm::ddt(psi, p)
- + fvc::div(phi)
- - fvm::laplacian(rho*rAU, p)
- ==
- parcels.Srho()
- + surfaceFilm.Srho()
- );
-
- pEqn.solve();
-
- if (nonOrth == nNonOrthCorr)
- {
- phi += pEqn.flux();
- }
+ phi += p_rghEqn.flux();
}
}
+p = p_rgh + rho*gh;
+
#include "rhoEqn.H"
#include "compressibleContinuityErrs.H"
-U -= rAU*fvc::grad(p);
+U += rAU*fvc::reconstruct((phi - phiU)/rhorAUf);
U.correctBoundaryConditions();
DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C b/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C
index f788da57e468464e44f4964b0a423ebabe60bf27..fdedebb78550ef5a0be066be06f7be306464668b 100644
--- a/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C
+++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2009-2010 OpenCFD Ltd.
+ \\ / A nd | Copyright (C) 2009-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -39,6 +39,7 @@ Description
#include "chemistrySolver.H"
#include "radiationModel.H"
#include "SLGThermo.H"
+#include "pimpleLoop.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -68,6 +69,7 @@ int main(int argc, char *argv[])
#include "readTimeControls.H"
#include "readPISOControls.H"
#include "compressibleCourantNo.H"
+ #include "setMultiRegionDeltaT.H"
#include "setDeltaT.H"
runTime++;
@@ -84,20 +86,22 @@ int main(int argc, char *argv[])
#include "rhoEqn.H"
// --- PIMPLE loop
- for (int ocorr=1; ocorr<=nOuterCorr; ocorr++)
+ for
+ (
+ pimpleLoop pimpleCorr(mesh, nOuterCorr);
+ pimpleCorr.loop();
+ pimpleCorr++
+ )
{
#include "UEqn.H"
#include "YEqn.H"
+ #include "hsEqn.H"
// --- PISO loop
for (int corr=1; corr<=nCorr; corr++)
{
- #include "hsEqn.H"
#include "pEqn.H"
}
-
- Info<< "T gas min/max = " << min(T).value() << ", "
- << max(T).value() << endl;
}
turbulence->correct();
diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/setMultiRegionDeltaT.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/setMultiRegionDeltaT.H
new file mode 100644
index 0000000000000000000000000000000000000000..ed3b584675f94ae862e5e1db06ade682fda0654a
--- /dev/null
+++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/setMultiRegionDeltaT.H
@@ -0,0 +1,60 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2010-2011 OpenCFD Ltd.
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
+
+Global
+ setMultiRegionDeltaT
+
+Description
+ Reset the timestep to maintain a constant maximum Courant numbers.
+ Reduction of time-step is immediate, but increase is damped to avoid
+ unstable oscillations.
+
+\*---------------------------------------------------------------------------*/
+
+if (adjustTimeStep)
+{
+ if (CoNum == -GREAT)
+ {
+ CoNum = SMALL;
+ }
+
+ const scalar TFactorFluid = maxCo/(CoNum + SMALL);
+ const scalar TFactorFilm = maxCo/(surfaceFilm.CourantNumber() + SMALL);
+
+ const scalar dt0 = runTime.deltaTValue();
+
+ const scalar TFactorFluidLim =
+ min(min(TFactorFluid, 1.0 + 0.1*TFactorFluid), 1.2);
+
+ runTime.setDeltaT
+ (
+ min
+ (
+ dt0*min(TFactorFluidLim, TFactorFilm),
+ maxDeltaT
+ )
+ );
+}
+
+
+// ************************************************************************* //