Final iteration information now available in mesh::data (used to be mesh::fvData)

Relaxation and solution parameters for the final iteration in PIMPLE loops are
now selected according to the value of the "finalIteration" entry in the
mesh::data dictionary.

rhoPimpleFoam significantly updates and now replaces rhoPisoFoam.
This commit is contained in:
henry
2010-05-25 18:45:25 +01:00
parent 49ccf0ffaa
commit 361b153343
161 changed files with 685 additions and 859 deletions

View File

@ -7,22 +7,6 @@ fvVectorMatrix UEqn
UEqn.relax();
if (oCorr == nOuterCorr - 1)
{
solve
(
UEqn
==
fvc::reconstruct
(
fvc::interpolate(rho)*(g & mesh.Sf())
- fvc::snGrad(p)*mesh.magSf()
),
mesh.solver("UFinal")
);
}
else
{
solve
(
UEqn
@ -33,4 +17,3 @@ else
- fvc::snGrad(p)*mesh.magSf()
)
);
}

View File

@ -70,6 +70,12 @@ int main(int argc, char *argv[])
// --- Pressure-velocity PIMPLE corrector loop
for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
{
bool finalIter = oCorr == nOuterCorr-1;
if (finalIter)
{
mesh.data::add("finalIteration", true);
}
#include "UEqn.H"
#include "ftEqn.H"
@ -80,6 +86,11 @@ int main(int argc, char *argv[])
{
#include "pEqn.H"
}
if (finalIter)
{
mesh.data::remove("finalIteration");
}
}
turbulence->correct();

View File

@ -30,14 +30,20 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
closedVolume = p.needReference();
if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
{
pEqn.solve(mesh.solver(p.name() + "Final"));
}
else
{
pEqn.solve(mesh.solver(p.name()));
}
pEqn.solve
(
mesh.solver
(
p.select
(
(
finalIter
&& corr == nCorr-1
&& nonOrth == nNonOrthCorr
)
)
)
);
if (nonOrth == nNonOrthCorr)
{

View File

@ -31,14 +31,20 @@
- fvm::laplacian(rho*rUA, p)
);
if (ocorr == nOuterCorr && corr == nCorr && nonOrth == nNonOrthCorr)
{
pEqn.solve(mesh.solver(p.name() + "Final"));
}
else
{
pEqn.solve();
}
pEqn.solve
(
mesh.solver
(
p.select
(
(
finalIter
&& corr == nCorr-1
&& nonOrth == nNonOrthCorr
)
)
)
);
if (nonOrth == nNonOrthCorr)
{
@ -64,14 +70,20 @@
- fvm::laplacian(rho*rUA, p)
);
if (ocorr == nOuterCorr && corr == nCorr && nonOrth == nNonOrthCorr)
{
pEqn.solve(mesh.solver(p.name() + "Final"));
}
else
{
pEqn.solve();
}
pEqn.solve
(
mesh.solver
(
p.select
(
(
finalIter
&& corr == nCorr-1
&& nonOrth == nNonOrthCorr
)
)
)
);
if (nonOrth == nNonOrthCorr)
{

View File

@ -69,8 +69,14 @@ int main(int argc, char *argv[])
#include "chemistry.H"
#include "rhoEqn.H"
for (label ocorr=1; ocorr <= nOuterCorr; ocorr++)
for (label oCorr=1; oCorr <= nOuterCorr; oCorr++)
{
bool finalIter = oCorr == nOuterCorr-1;
if (finalIter)
{
mesh.data::add("finalIteration", true);
}
#include "UEqn.H"
#include "YEqn.H"
#include "hsEqn.H"
@ -80,6 +86,11 @@ int main(int argc, char *argv[])
{
#include "pEqn.H"
}
if (finalIter)
{
mesh.data::remove("finalIteration");
}
}
turbulence->correct();

View File

@ -7,28 +7,14 @@ tmp<fvVectorMatrix> UEqn
+ turbulence->divDevRhoReff(U)
);
if (oCorr == nOuterCorr-1)
{
UEqn().relax(1);
}
else
{
UEqn().relax();
}
volScalarField rUA = 1.0/UEqn().A();
if (momentumPredictor)
{
if (oCorr == nOuterCorr-1)
{
solve(UEqn() == -fvc::grad(p), mesh.solver("UFinal"));
}
else
{
solve(UEqn() == -fvc::grad(p));
}
}
else
{
U = rUA*(UEqn().H() - fvc::grad(p));

View File

@ -39,11 +39,6 @@
#include "compressibleCreatePhi.H"
dimensionedScalar pMin
(
mesh.solutionDict().subDict("PIMPLE").lookup("pMin")
);
Info<< "Creating turbulence model\n" << endl;
autoPtr<compressible::turbulenceModel> turbulence
(
@ -56,9 +51,6 @@
)
);
//dimensionedScalar initialMass = fvc::domainIntegrate(rho);
Info<< "Creating field DpDt\n" << endl;
volScalarField DpDt =
fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);

View File

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

View File

@ -1,6 +1,5 @@
rho = thermo.rho();
volScalarField rUA = 1.0/UEqn().A();
U = rUA*UEqn().H();
if (nCorr <= 1)
@ -29,19 +28,20 @@ if (transonic)
- fvm::laplacian(rho*rUA, p)
);
if
pEqn.solve
(
oCorr == nOuterCorr-1
mesh.solver
(
p.select
(
(
finalIter
&& corr == nCorr-1
&& nonOrth == nNonOrthCorr
)
{
pEqn.solve(mesh.solver("pFinal"));
}
else
{
pEqn.solve();
}
)
)
);
if (nonOrth == nNonOrthCorr)
{
@ -55,7 +55,7 @@ else
fvc::interpolate(rho)*
(
(fvc::interpolate(U) & mesh.Sf())
//+ fvc::ddtPhiCorr(rUA, rho, U, phi)
+ fvc::ddtPhiCorr(rUA, rho, U, phi)
);
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
@ -68,19 +68,20 @@ else
- fvm::laplacian(rho*rUA, p)
);
if
pEqn.solve
(
oCorr == nOuterCorr-1
mesh.solver
(
p.select
(
(
finalIter
&& corr == nCorr-1
&& nonOrth == nNonOrthCorr
)
{
pEqn.solve(mesh.solver("pFinal"));
}
else
{
pEqn.solve();
}
)
)
);
if (nonOrth == nNonOrthCorr)
{
@ -92,30 +93,15 @@ else
#include "rhoEqn.H"
#include "compressibleContinuityErrs.H"
//if (oCorr != nOuterCorr-1)
{
// Explicitly relax pressure for momentum corrector
p.relax();
// Recalculate density from the relaxed pressure
rho = thermo.rho();
rho.relax();
Info<< "rho max/min : " << max(rho).value()
<< " " << min(rho).value() << endl;
}
U -= rUA*fvc::grad(p);
U.correctBoundaryConditions();
DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
bound(p, pMin);
// 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);
}
*/

View File

@ -61,15 +61,22 @@ int main(int argc, char *argv[])
Info<< "Time = " << runTime.timeName() << nl << endl;
if (nOuterCorr != 1)
{
p.storePrevIter();
rho.storePrevIter();
}
#include "rhoEqn.H"
// --- Pressure-velocity PIMPLE corrector loop
for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
{
bool finalIter = oCorr == nOuterCorr-1;
if (finalIter)
{
mesh.data::add("finalIteration", true);
}
if (nOuterCorr != 1)
{
p.storePrevIter();
}
#include "UEqn.H"
#include "hEqn.H"
@ -80,6 +87,11 @@ int main(int argc, char *argv[])
}
turbulence->correct();
if (finalIter)
{
mesh.data::remove("finalIteration");
}
}
runTime.write();

View File

@ -1,3 +0,0 @@
rhoPisoFoam.C
EXE = $(FOAM_APPBIN)/rhoPisoFoam

View File

@ -1,13 +0,0 @@
EXE_INC = \
-I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
-I$(LIB_SRC)/finiteVolume/lnInclude
EXE_LIBS = \
-lfiniteVolume \
-lmeshTools \
-lbasicThermophysicalModels \
-lspecie \
-lcompressibleRASModels \
-lcompressibleLESModels

View File

@ -1,11 +0,0 @@
fvVectorMatrix UEqn
(
fvm::ddt(rho, U)
+ fvm::div(phi, U)
+ turbulence->divDevRhoReff(U)
);
if (momentumPredictor)
{
solve(UEqn == -fvc::grad(p));
}

View File

@ -1,58 +0,0 @@
Info<< "Reading thermophysical properties\n" << endl;
autoPtr<basicPsiThermo> pThermo
(
basicPsiThermo::New(mesh)
);
basicPsiThermo& thermo = pThermo();
volScalarField& p = thermo.p();
volScalarField& h = thermo.h();
const volScalarField& psi = thermo.psi();
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
thermo.rho()
);
Info<< "\nReading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
# include "compressibleCreatePhi.H"
Info<< "Creating turbulence model\n" << endl;
autoPtr<compressible::turbulenceModel> turbulence
(
compressible::turbulenceModel::New
(
rho,
U,
phi,
thermo
)
);
Info<< "Creating field DpDt\n" << endl;
volScalarField DpDt =
fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);

View File

@ -1,12 +0,0 @@
{
solve
(
fvm::ddt(rho, h)
+ fvm::div(phi, h)
- fvm::laplacian(turbulence->alphaEff(), h)
==
DpDt
);
thermo.correct();
}

View File

@ -1,68 +0,0 @@
rho = thermo.rho();
volScalarField rUA = 1.0/UEqn.A();
U = rUA*UEqn.H();
if (transonic)
{
surfaceScalarField phid
(
"phid",
fvc::interpolate(psi)
*(
(fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi)
)
);
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
(
fvm::ddt(psi, p)
+ fvm::div(phid, p)
- fvm::laplacian(rho*rUA, p)
);
pEqn.solve();
if (nonOrth == nNonOrthCorr)
{
phi == pEqn.flux();
}
}
}
else
{
phi =
fvc::interpolate(rho)*
(
(fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi)
);
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
(
fvm::ddt(psi, p)
+ fvc::div(phi)
- fvm::laplacian(rho*rUA, p)
);
pEqn.solve();
if (nonOrth == nNonOrthCorr)
{
phi += pEqn.flux();
}
}
}
#include "rhoEqn.H"
#include "compressibleContinuityErrs.H"
U -= rUA*fvc::grad(p);
U.correctBoundaryConditions();
DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);

View File

@ -1,93 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 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/>.
Application
rhoPisoFoam
Description
Transient PISO solver for compressible, laminar or turbulent flow.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "basicPsiThermo.H"
#include "turbulenceModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "createFields.H"
#include "initContinuityErrs.H"
#include "readTimeControls.H"
#include "compressibleCourantNo.H"
#include "setInitialDeltaT.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.run())
{
#include "readTimeControls.H"
#include "readPISOControls.H"
#include "compressibleCourantNo.H"
#include "setDeltaT.H"
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
#include "rhoEqn.H"
#include "UEqn.H"
// --- PISO loop
for (int corr=1; corr<=nCorr; corr++)
{
#include "hEqn.H"
#include "pEqn.H"
}
turbulence->correct();
rho = thermo.rho();
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //

View File

@ -7,14 +7,7 @@ tmp<fvVectorMatrix> UEqn
+ turbulence->divDevRhoReff(U)
);
if (oCorr == nOuterCorr-1)
{
UEqn().relax(1);
}
else
{
UEqn().relax();
}
mrfZones.addCoriolis(rho, UEqn());
pZones.addResistance(UEqn());
@ -22,16 +15,9 @@ pZones.addResistance(UEqn());
volScalarField rUA = 1.0/UEqn().A();
if (momentumPredictor)
{
if (oCorr == nOuterCorr-1)
{
solve(UEqn() == -fvc::grad(p), mesh.solver("UFinal"));
}
else
{
solve(UEqn() == -fvc::grad(p));
}
}
else
{
U = rUA*(UEqn().H() - fvc::grad(p));

View File

@ -30,19 +30,20 @@ if (transonic)
- fvm::laplacian(rho*rUA, p)
);
if
pEqn.solve
(
oCorr == nOuterCorr-1
mesh.solver
(
p.select
(
(
finalIter
&& corr == nCorr-1
&& nonOrth == nNonOrthCorr
)
{
pEqn.solve(mesh.solver("pFinal"));
}
else
{
pEqn.solve();
}
)
)
);
if (nonOrth == nNonOrthCorr)
{
@ -70,19 +71,20 @@ else
- fvm::laplacian(rho*rUA, p)
);
if
pEqn.solve
(
oCorr == nOuterCorr-1
mesh.solver
(
p.select
(
(
finalIter
&& corr == nCorr-1
&& nonOrth == nNonOrthCorr
)
{
pEqn.solve(mesh.solver("pFinal"));
}
else
{
pEqn.solve();
}
)
)
);
if (nonOrth == nNonOrthCorr)
{

View File

@ -65,17 +65,23 @@ int main(int argc, char *argv[])
Info<< "Time = " << runTime.timeName() << nl << endl;
#include "rhoEqn.H"
// --- Pressure-velocity PIMPLE corrector loop
for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
{
bool finalIter = oCorr == nOuterCorr-1;
if (finalIter)
{
mesh.data::add("finalIteration", true);
}
if (nOuterCorr != 1)
{
p.storePrevIter();
rho.storePrevIter();
}
#include "rhoEqn.H"
// --- Pressure-velocity PIMPLE corrector loop
for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
{
#include "UEqn.H"
#include "hEqn.H"
@ -86,6 +92,11 @@ int main(int argc, char *argv[])
}
turbulence->correct();
if (finalIter)
{
mesh.data::remove("finalIteration");
}
}
runTime.write();

View File

@ -45,9 +45,14 @@
scalar pRefValue = 0.0;
setRefCell(p, mesh.solutionDict().subDict("SIMPLE"), pRefCell, pRefValue);
dimensionedScalar pMin
dimensionedScalar rhoMax
(
mesh.solutionDict().subDict("SIMPLE").lookup("pMin")
mesh.solutionDict().subDict("SIMPLE").lookup("rhoMax")
);
dimensionedScalar rhoMin
(
mesh.solutionDict().subDict("SIMPLE").lookup("rhoMin")
);
Info<< "Creating turbulence model\n" << endl;

View File

@ -5,8 +5,8 @@
- fvm::Sp(fvc::div(phi), h)
- fvm::laplacian(turbulence->alphaEff(), h)
==
fvc::div(phi/fvc::interpolate(rho)*fvc::interpolate(p, "div(U,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))
);
pZones.addEnthalpySource(thermo, rho, hEqn);

View File

@ -1,3 +1,8 @@
rho = thermo.rho();
rho = max(rho, rhoMin);
rho = min(rho, rhoMax);
rho.relax();
if (pressureImplicitPorosity)
{
U = trTU()&UEqn().H();
@ -58,8 +63,6 @@ else
U.correctBoundaryConditions();
bound(p, pMin);
// For closed-volume cases adjust the pressure and density levels
// to obey overall mass continuity
if (closedVolume)
@ -69,5 +72,7 @@ if (closedVolume)
}
rho = thermo.rho();
rho = max(rho, rhoMin);
rho = min(rho, rhoMax);
rho.relax();
Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;

View File

@ -0,0 +1,3 @@
buoyantBoussinesqPimpleFoam.C
EXE = $(FOAM_APPBIN)/buoyantBoussinesqPimpleFoam

View File

@ -13,7 +13,6 @@
);
TEqn.relax();
TEqn.solve();
rhok = 1.0 - beta*(T - TRef);

View File

@ -22,7 +22,7 @@ License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Application
buoyantBoussinesqPisoFoam
buoyantBoussinesqPimpleFoam
Description
Transient solver for buoyant, turbulent flow of incompressible fluids
@ -72,10 +72,24 @@ int main(int argc, char *argv[])
Info<< "Time = " << runTime.timeName() << nl << endl;
#include "readTimeControls.H"
#include "readPISOControls.H"
#include "readPIMPLEControls.H"
#include "CourantNo.H"
#include "setDeltaT.H"
// --- Pressure-velocity PIMPLE corrector loop
for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
{
bool finalIter = oCorr == nOuterCorr-1;
if (finalIter)
{
mesh.data::add("finalIteration", true);
}
if (nOuterCorr != 1)
{
p.storePrevIter();
}
#include "UEqn.H"
#include "TEqn.H"
@ -87,6 +101,12 @@ int main(int argc, char *argv[])
turbulence->correct();
if (finalIter)
{
mesh.data::remove("finalIteration");
}
}
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

View File

@ -57,7 +57,7 @@
setRefCell
(
p,
mesh.solutionDict().subDict("PISO"),
mesh.solutionDict().subDict("PIMPLE"),
pRefCell,
pRefValue
);

View File

@ -0,0 +1,52 @@
{
volScalarField rUA("rUA", 1.0/UEqn.A());
surfaceScalarField rUAf("(1|A(U))", fvc::interpolate(rUA));
U = rUA*UEqn.H();
phi = (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, U, phi);
surfaceScalarField buoyancyPhi =
rUAf*fvc::interpolate(rhok)*(g & mesh.Sf());
phi += buoyancyPhi;
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
(
fvm::laplacian(rUAf, p) == fvc::div(phi)
);
pEqn.solve
(
mesh.solver
(
p.select
(
(
finalIter
&& corr == nCorr-1
&& nonOrth == nNonOrthCorr
)
)
)
);
if (nonOrth == nNonOrthCorr)
{
// Calculate the conservative fluxes
phi -= pEqn.flux();
// Explicitly relax pressure for momentum corrector
p.relax();
// Correct the momentum source with the pressure gradient flux
// calculated from the relaxed pressure
U += rUA*fvc::reconstruct((buoyancyPhi - pEqn.flux())/rUAf);
U.correctBoundaryConditions();
}
}
#include "continuityErrs.H"
}

View File

@ -1,3 +0,0 @@
buoyantBoussinesqPisoFoam.C
EXE = $(FOAM_APPBIN)/buoyantBoussinesqPisoFoam

View File

@ -1,41 +0,0 @@
{
volScalarField rUA("rUA", 1.0/UEqn.A());
surfaceScalarField rUAf("(1|A(U))", fvc::interpolate(rUA));
U = rUA*UEqn.H();
surfaceScalarField phiU
(
(fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, U, phi)
);
phi = phiU + rUAf*fvc::interpolate(rhok)*(g & mesh.Sf());
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
(
fvm::laplacian(rUAf, p) == fvc::div(phi)
);
if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
{
pEqn.solve(mesh.solver(p.name() + "Final"));
}
else
{
pEqn.solve(mesh.solver(p.name()));
}
if (nonOrth == nNonOrthCorr)
{
phi -= pEqn.flux();
}
}
U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
U.correctBoundaryConditions();
#include "continuityErrs.H"
}

View File

@ -0,0 +1,3 @@
buoyantPimpleFoam.C
EXE = $(FOAM_APPBIN)/buoyantPimpleFoam

View File

@ -22,7 +22,7 @@ License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Application
buoyantPisoFoam
buoyantPimpleFoam
Description
Transient solver for buoyant, turbulent flow of compressible fluids for
@ -59,7 +59,7 @@ int main(int argc, char *argv[])
while (runTime.run())
{
#include "readTimeControls.H"
#include "readPISOControls.H"
#include "readPIMPLEControls.H"
#include "compressibleCourantNo.H"
#include "setDeltaT.H"
@ -69,12 +69,24 @@ int main(int argc, char *argv[])
#include "rhoEqn.H"
#include "UEqn.H"
// --- Pressure-velocity PIMPLE corrector loop
for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
{
bool finalIter = oCorr == nOuterCorr-1;
if (finalIter)
{
mesh.data::add("finalIteration", true);
}
if (nOuterCorr != 1)
{
p.storePrevIter();
}
#include "UEqn.H"
#include "hEqn.H"
// --- PISO loop
for (int corr=0; corr<nCorr; corr++)
{
#include "pEqn.H"
@ -84,6 +96,12 @@ int main(int argc, char *argv[])
rho = thermo.rho();
if (finalIter)
{
mesh.data::remove("finalIteration");
}
}
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

View File

@ -12,16 +12,15 @@
U = rUA*UEqn.H();
surfaceScalarField phiU
phi = fvc::interpolate(rho)*
(
fvc::interpolate(rho)
*(
(fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi)
)
);
phi = phiU + rhorUAf*fvc::interpolate(rho)*(g & mesh.Sf());
surfaceScalarField buoyancyPhi =
rhorUAf*fvc::interpolate(rho)*(g & mesh.Sf());
phi += buoyancyPhi;
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
@ -32,27 +31,39 @@
- fvm::laplacian(rhorUAf, p)
);
if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
{
pEqn.solve(mesh.solver(p.name() + "Final"));
}
else
{
pEqn.solve(mesh.solver(p.name()));
}
pEqn.solve
(
mesh.solver
(
p.select
(
(
finalIter
&& corr == nCorr-1
&& nonOrth == nNonOrthCorr
)
)
)
);
if (nonOrth == nNonOrthCorr)
{
// Calculate the conservative fluxes
phi += pEqn.flux();
// Explicitly relax pressure for momentum corrector
p.relax();
// Correct the momentum source with the pressure gradient flux
// calculated from the relaxed pressure
U += rUA*fvc::reconstruct((buoyancyPhi + pEqn.flux())/rhorUAf);
U.correctBoundaryConditions();
}
}
// Second part of thermodynamic density update
thermo.rho() += psi*p;
U += rUA*fvc::reconstruct((phi - phiU)/rhorUAf);
U.correctBoundaryConditions();
DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
#include "rhoEqn.H"

View File

@ -1,3 +0,0 @@
buoyantPisoFoam.C
EXE = $(FOAM_APPBIN)/buoyantPisoFoam

View File

@ -22,16 +22,7 @@
);
pEqn.setReference(pRefCell, pRefValue);
// retain the residual from the first iteration
if (nonOrth == 0)
{
pEqn.solve();
}
else
{
pEqn.solve();
}
if (nonOrth == nNonOrthCorr)
{

View File

@ -1,16 +0,0 @@
fvVectorMatrix UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
+ turbulence->divDevReff(U)
);
if (ocorr != nOuterCorr-1)
{
UEqn.relax();
}
if (momentumPredictor)
{
solve(UEqn == -fvc::grad(p));
}

View File

@ -0,0 +1,8 @@
#!/bin/sh
cd ${0%/*} || exit 1 # run from this directory
set -x
wmake
wmake pimpleDyMFoam
# ----------------------------------------------------------------- end-of-file

View File

@ -7,28 +7,14 @@ tmp<fvVectorMatrix> UEqn
+ turbulence->divDevReff(U)
);
if (oCorr == nOuterCorr-1)
{
UEqn().relax(1);
}
else
{
UEqn().relax();
}
volScalarField rUA = 1.0/UEqn().A();
if (momentumPredictor)
{
if (oCorr == nOuterCorr-1)
{
solve(UEqn() == -fvc::grad(p), mesh.solver("UFinal"));
}
else
{
solve(UEqn() == -fvc::grad(p));
}
}
else
{
U = rUA*(UEqn().H() - fvc::grad(p));

View File

@ -21,19 +21,20 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
pEqn.setReference(pRefCell, pRefValue);
if
pEqn.solve
(
oCorr == nOuterCorr-1
mesh.solver
(
p.select
(
(
finalIter
&& corr == nCorr-1
&& nonOrth == nNonOrthCorr
)
{
pEqn.solve(mesh.solver("pFinal"));
}
else
{
pEqn.solve();
}
)
)
);
if (nonOrth == nNonOrthCorr)
{
@ -43,11 +44,8 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
#include "continuityErrs.H"
// Explicitly relax pressure for momentum corrector except for last corrector
if (oCorr != nOuterCorr-1)
{
// Explicitly relax pressure for momentum corrector
p.relax();
}
U -= rUA*fvc::grad(p);
U.correctBoundaryConditions();

View File

@ -1,4 +1,5 @@
EXE_INC = \
-I.. \
-I$(LIB_SRC)/dynamicFvMesh/lnInclude \
-I$(LIB_SRC)/dynamicMesh/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \

View File

@ -0,0 +1,55 @@
U = rAU*UEqn().H();
if (nCorr <= 1)
{
UEqn.clear();
}
phi = (fvc::interpolate(U) & mesh.Sf());
if (p.needReference())
{
fvc::makeRelative(phi, U);
adjustPhi(phi, U, p);
fvc::makeAbsolute(phi, U);
}
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
(
fvm::laplacian(rAU, p) == fvc::div(phi)
);
pEqn.setReference(pRefCell, pRefValue);
pEqn.solve
(
mesh.solver
(
p.select
(
(
finalIter
&& corr == nCorr-1
&& nonOrth == nNonOrthCorr
)
)
)
);
if (nonOrth == nNonOrthCorr)
{
phi -= pEqn.flux();
}
}
#include "continuityErrs.H"
p.relax();
// Make the fluxes relative to the mesh motion
fvc::makeRelative(phi, U);
U -= rAU*fvc::grad(p);
U.correctBoundaryConditions();

View File

@ -84,8 +84,14 @@ int main(int argc, char *argv[])
}
// --- PIMPLE loop
for (int ocorr=0; ocorr<nOuterCorr; ocorr++)
for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
{
bool finalIter = oCorr == nOuterCorr-1;
if (finalIter)
{
mesh.data::add("finalIteration", true);
}
if (nOuterCorr != 1)
{
p.storePrevIter();
@ -96,62 +102,15 @@ int main(int argc, char *argv[])
// --- PISO loop
for (int corr=0; corr<nCorr; corr++)
{
rAU = 1.0/UEqn.A();
U = rAU*UEqn.H();
phi = (fvc::interpolate(U) & mesh.Sf());
if (p.needReference())
{
fvc::makeRelative(phi, U);
adjustPhi(phi, U, p);
fvc::makeAbsolute(phi, U);
}
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
(
fvm::laplacian(rAU, p) == fvc::div(phi)
);
pEqn.setReference(pRefCell, pRefValue);
if
(
ocorr == nOuterCorr-1
&& corr == nCorr-1
&& nonOrth == nNonOrthCorr)
{
pEqn.solve(mesh.solver(p.name() + "Final"));
}
else
{
pEqn.solve(mesh.solver(p.name()));
}
if (nonOrth == nNonOrthCorr)
{
phi -= pEqn.flux();
}
}
#include "continuityErrs.H"
// Explicitly relax pressure for momentum corrector
if (ocorr != nOuterCorr-1)
{
p.relax();
}
// Make the fluxes relative to the mesh motion
fvc::makeRelative(phi, U);
U -= rAU*fvc::grad(p);
U.correctBoundaryConditions();
#include "pEqn.H"
}
turbulence->correct();
if (finalIter)
{
mesh.data::remove("finalIteration");
}
}
runTime.write();

View File

@ -62,6 +62,12 @@ int main(int argc, char *argv[])
// --- Pressure-velocity PIMPLE corrector loop
for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
{
bool finalIter = oCorr == nOuterCorr-1;
if (finalIter)
{
mesh.data::add("finalIteration", true);
}
if (nOuterCorr != 1)
{
p.storePrevIter();
@ -76,6 +82,11 @@ int main(int argc, char *argv[])
}
turbulence->correct();
if (finalIter)
{
mesh.data::remove("finalIteration");
}
}
runTime.write();

View File

@ -15,14 +15,7 @@
//- fvc::div(rho*turbulence->nuEff()*dev(fvc::grad(U)().T()))
);
if (oCorr == nOuterCorr-1)
{
UEqn.relax(1);
}
else
{
UEqn.relax();
}
if (momentumPredictor)
{
@ -34,7 +27,6 @@
(
fvc::interpolate(rho)*(g & mesh.Sf())
- mesh.magSf()*fvc::snGrad(p)
),
mesh.solver(oCorr == nOuterCorr-1 ? "UFinal" : "U")
)
);
}

View File

@ -22,18 +22,20 @@
pEqn.setReference(pRefCell, pRefValue);
if
pEqn.solve
(
corr == nCorr-1
mesh.solver
(
p.select
(
(
finalIter
&& corr == nCorr-1
&& nonOrth == nNonOrthCorr
)
{
pEqn.solve(mesh.solver(p.name() + "Final"));
}
else
{
pEqn.solve(mesh.solver(p.name()));
}
)
)
);
if (nonOrth == nNonOrthCorr)
{

View File

@ -68,6 +68,12 @@ int main(int argc, char *argv[])
// --- Pressure-velocity PIMPLE corrector loop
for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
{
bool finalIter = oCorr == nOuterCorr-1;
if (finalIter)
{
mesh.data::add("finalIteration", true);
}
twoPhaseProperties.correct();
#include "alphaEqn.H"
@ -83,6 +89,11 @@ int main(int argc, char *argv[])
#include "continuityErrs.H"
turbulence->correct();
if (finalIter)
{
mesh.data::remove("finalIteration");
}
}
runTime.write();

View File

@ -543,4 +543,6 @@ $(writers)/gnuplotGraph/gnuplotGraph.C
$(writers)/xmgrGraph/xmgrGraph.C
$(writers)/jplotGraph/jplotGraph.C
meshes/data/data.C
LIB = $(FOAM_LIBBIN)/libOpenFOAM

View File

@ -27,6 +27,7 @@ License
#include "Time.H"
#include "demandDrivenData.H"
#include "dictionary.H"
#include "data.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -918,6 +919,15 @@ bool Foam::GeometricField<Type, PatchField, GeoMesh>::needReference() const
template<class Type, template<class> class PatchField, class GeoMesh>
void Foam::GeometricField<Type, PatchField, GeoMesh>::relax(const scalar alpha)
{
if (debug)
{
InfoIn
(
"GeometricField<Type, PatchField, GeoMesh>::relax"
"(const scalar alpha)"
) << "Relaxing" << endl << this->info() << " by " << alpha << endl;
}
operator==(prevIter() + alpha*(*this - prevIter()));
}
@ -925,16 +935,33 @@ void Foam::GeometricField<Type, PatchField, GeoMesh>::relax(const scalar alpha)
template<class Type, template<class> class PatchField, class GeoMesh>
void Foam::GeometricField<Type, PatchField, GeoMesh>::relax()
{
scalar alpha = 0;
word name = this->name();
if (this->mesh().relax(this->name()))
if (this->mesh().data::lookupOrDefault<bool>("finalIteration", false))
{
alpha = this->mesh().relaxationFactor(this->name());
name += "Final";
}
if (alpha > 0)
if (this->mesh().relax(name))
{
relax(alpha);
relax(this->mesh().relaxationFactor(name));
}
}
template<class Type, template<class> class PatchField, class GeoMesh>
Foam::word Foam::GeometricField<Type, PatchField, GeoMesh>::select
(
bool final
) const
{
if (final)
{
return this->name() + "Final";
}
else
{
return this->name();
}
}

View File

@ -476,6 +476,11 @@ public:
// alpha is read from controlDict
void relax();
//- Select the final iteration parameters if `final' is true
// by returning the field name + "Final"
// otherwise the standard parameters by returning the field name
word select(bool final) const;
// Member function *this operators

View File

@ -23,23 +23,23 @@ License
\*---------------------------------------------------------------------------*/
#include "fvData.H"
#include "data.H"
#include "Time.H"
#include "lduMatrix.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
int Foam::fvData::debug(Foam::debug::debugSwitch("fvData", false));
int Foam::data::debug(Foam::debug::debugSwitch("data", false));
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::fvData::fvData(const objectRegistry& obr)
Foam::data::data(const objectRegistry& obr)
:
IOdictionary
(
IOobject
(
"fvData",
"data",
obr.time().system(),
obr,
IOobject::NO_READ,
@ -53,13 +53,13 @@ Foam::fvData::fvData(const objectRegistry& obr)
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
const Foam::dictionary& Foam::fvData::solverPerformanceDict() const
const Foam::dictionary& Foam::data::solverPerformanceDict() const
{
return subDict("solverPerformance");
}
void Foam::fvData::setSolverPerformance
void Foam::data::setSolverPerformance
(
const word& name,
const lduMatrix::solverPerformance& sp

View File

@ -22,20 +22,21 @@ License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::fvData
Foam::data
Description
Database for finite volume solution data, solver performance and
other reduced data. fvMesh is derived from fvData so that all fields have
access to the fvData from the mesh reference they hold.
Database for solution data, solver performance and other reduced data.
fvMesh is derived from data so that all fields have access to the data from
the mesh reference they hold.
SourceFiles
fvData.C
data.C
\*---------------------------------------------------------------------------*/
#ifndef fvData_H
#define fvData_H
#ifndef data_H
#define data_H
#include "IOdictionary.H"
#include "lduMatrix.H"
@ -46,20 +47,20 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
Class fvData Declaration
Class data Declaration
\*---------------------------------------------------------------------------*/
class fvData
class data
:
public IOdictionary
{
// Private Member Functions
//- Disallow default bitwise copy construct
fvData(const fvData&);
data(const data&);
//- Disallow default bitwise assignment
void operator=(const fvData&);
void operator=(const data&);
public:
@ -71,7 +72,7 @@ public:
// Constructors
//- Construct for objectRegistry
fvData(const objectRegistry& obr);
data(const objectRegistry& obr);
// Member Functions

View File

@ -270,7 +270,6 @@ $(multivariateSchemes)/limitedCubic/multivariateLimitedCubic.C
finiteVolume/fv/fv.C
finiteVolume/fvSchemes/fvSchemes.C
finiteVolume/fvData/fvData.C
ddtSchemes = finiteVolume/ddtSchemes
$(ddtSchemes)/ddtScheme/ddtSchemes.C

View File

@ -506,6 +506,13 @@ void Foam::fvMatrix<Type>::relax(const scalar alpha)
return;
}
if (debug)
{
InfoIn("fvMatrix<Type>::relax(const scalar alpha)")
<< "Relaxing " << psi_.name() << " by " << alpha
<< endl;
}
Field<Type>& S = source();
scalarField& D = diag();
@ -591,9 +598,14 @@ void Foam::fvMatrix<Type>::relax(const scalar alpha)
template<class Type>
void Foam::fvMatrix<Type>::relax()
{
if (psi_.mesh().relax(psi_.name()))
word name = psi_.select
(
psi_.mesh().data::lookupOrDefault<bool>("finalIteration", false)
);
if (psi_.mesh().relax(name))
{
relax(psi_.mesh().relaxationFactor(psi_.name()));
relax(psi_.mesh().relaxationFactor(name));
}
}

View File

@ -167,21 +167,48 @@ template<class Type>
Foam::autoPtr<typename Foam::fvMatrix<Type>::fvSolver>
Foam::fvMatrix<Type>::solver()
{
return solver(psi_.mesh().solverDict(psi_.name()));
return solver
(
psi_.mesh().solverDict
(
psi_.select
(
psi_.mesh().data::lookupOrDefault<bool>("finalIteration", false)
)
)
);
}
template<class Type>
Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::fvSolver::solve()
{
return solve(psi_.mesh().solverDict(psi_.name()));
return solve
(
psi_.mesh().solverDict
(
psi_.select
(
psi_.mesh().data::lookupOrDefault<bool>("finalIteration", false)
)
)
);
}
template<class Type>
Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::solve()
{
return solve(psi_.mesh().solverDict(psi_.name()));
return solve
(
psi_.mesh().solverDict
(
psi_.select
(
psi_.mesh().data::lookupOrDefault<bool>("finalIteration", false)
)
)
);
}

View File

@ -153,7 +153,7 @@ Foam::fvMesh::fvMesh(const IOobject& io)
surfaceInterpolation(*this),
fvSchemes(static_cast<const objectRegistry&>(*this)),
fvSolution(static_cast<const objectRegistry&>(*this)),
fvData(static_cast<const objectRegistry&>(*this)),
data(static_cast<const objectRegistry&>(*this)),
boundary_(*this, boundaryMesh()),
lduPtr_(NULL),
curTimeIndex_(time().timeIndex()),
@ -248,7 +248,7 @@ Foam::fvMesh::fvMesh
surfaceInterpolation(*this),
fvSchemes(static_cast<const objectRegistry&>(*this)),
fvSolution(static_cast<const objectRegistry&>(*this)),
fvData(static_cast<const objectRegistry&>(*this)),
data(static_cast<const objectRegistry&>(*this)),
boundary_(*this),
lduPtr_(NULL),
curTimeIndex_(time().timeIndex()),
@ -281,7 +281,7 @@ Foam::fvMesh::fvMesh
surfaceInterpolation(*this),
fvSchemes(static_cast<const objectRegistry&>(*this)),
fvSolution(static_cast<const objectRegistry&>(*this)),
fvData(static_cast<const objectRegistry&>(*this)),
data(static_cast<const objectRegistry&>(*this)),
boundary_(*this),
lduPtr_(NULL),
curTimeIndex_(time().timeIndex()),

View File

@ -54,7 +54,7 @@ SourceFiles
#include "surfaceInterpolation.H"
#include "fvSchemes.H"
#include "fvSolution.H"
#include "fvData.H"
#include "data.H"
#include "DimensionedField.H"
#include "volFieldsFwd.H"
#include "surfaceFieldsFwd.H"
@ -83,7 +83,7 @@ class fvMesh
public surfaceInterpolation,
public fvSchemes,
public fvSolution,
public fvData
public data
{
// Private data

View File

@ -1,7 +1,7 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.6 |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/

View File

@ -15,7 +15,7 @@ FoamFile
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application rhoPisoFoam;
application rhoPimpleFoam;
startFrom startTime;

View File

@ -23,15 +23,13 @@ ddtSchemes
gradSchemes
{
default Gauss linear;
grad(p) Gauss linear;
grad(U) Gauss linear;
}
divSchemes
{
default none;
div(phi,U) Gauss linear;
div(phi,h) Gauss linear;
div(phi,U) Gauss filteredLinear2V 0.2 0;
div(phi,h) Gauss filteredLinear2 0.2 0;
div(phiU,p) Gauss linear;
div(phi,k) Gauss limitedLinear 1;
div(phi,B) Gauss limitedLinear 1;
@ -54,7 +52,6 @@ laplacianSchemes
interpolationSchemes
{
default linear;
interpolate(HbyA) linear;
}
snGradSchemes

View File

@ -21,31 +21,42 @@ solvers
{
solver PCG;
preconditioner DIC;
tolerance 1e-06;
tolerance 1e-6;
relTol 0;
}
h
"(p|rho)Final"
{
$p;
relTol 0;
}
"(U|h|k|nuTilda)"
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-06;
tolerance 1e-6;
relTol 0;
}
"(U|k|B|nuTilda)"
"(U|h|k|nuTilda)Final"
{
$h;
tolerance 1e-05;
$U;
relTol 0;
}
}
PISO
PIMPLE
{
nCorrectors 2;
momentumPredictor yes;
nOuterCorrectors 2;
nCorrectors 1;
nNonOrthogonalCorrectors 0;
}
relaxationFactors
{
"(U|h|k|epsilon|omega).*" 1;
}
// ************************************************************************* //

View File

@ -32,14 +32,7 @@ solvers
relTol 0;
}
rho
{
$p;
tolerance 1e-05;
relTol 0;
}
U
"(rho|U|h|k|epsilon|omega)"
{
solver PBiCG;
preconditioner DILU;
@ -47,9 +40,7 @@ solvers
relTol 0.1;
}
h { $U; }
"(UFinal|hFinal|R|k|epsilon|omega)"
"(rho|U|h|k|epsilon|omega)Final"
{
$U;
tolerance 1e-05;
@ -60,20 +51,18 @@ solvers
PIMPLE
{
momentumPredictor yes;
nOuterCorrectors 50;
nCorrectors 1;
nNonOrthogonalCorrectors 0;
momentumPredictor yes;
pMin pMin [ 1 -1 -2 0 0 0 0 ] 1000;
}
relaxationFactors
{
p 0.3;
rho 0.05;
U 0.7;
h 0.7;
"(k|epsilon|omega)" 0.7;
"p.*" 0.3;
"rho.*" 1;
"(U|h|k|epsilon|omega)" 0.7;
"(U|h|k|epsilon|omega)Final" 1;
}

Some files were not shown because too many files have changed in this diff Show More