Temporary: consistent solver examples
This commit is contained in:
Hrvoje Jasak
parent
5766de6e54
commit
2d5052d714
20 changed files with 813 additions and 0 deletions
|
|
@ -0,0 +1,3 @@
|
|||
consistentIcoFoam.C
|
||||
|
||||
EXE = $(FOAM_APPBIN)/consistentIcoFoam
|
||||
|
|
@ -0,0 +1,6 @@
|
|||
EXE_INC = \
|
||||
-I$(LIB_SRC)/finiteVolume/lnInclude
|
||||
|
||||
EXE_LIBS = \
|
||||
-lfiniteVolume \
|
||||
-llduSolvers
|
||||
|
|
@ -0,0 +1,125 @@
|
|||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | foam-extend: Open Source CFD
|
||||
\\ / O peration | Version: 3.2
|
||||
\\ / A nd | Web: http://www.foam-extend.org
|
||||
\\/ M anipulation | For copyright notice see file Copyright
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
This file is part of foam-extend.
|
||||
|
||||
foam-extend 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.
|
||||
|
||||
foam-extend 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 foam-extend. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
Application
|
||||
consistentIcoFoam
|
||||
|
||||
Description
|
||||
Transient solver for incompressible, laminar flow of Newtonian fluids.
|
||||
Consistent formulation without time-step and relaxation dependence by Jasak
|
||||
|
||||
Author
|
||||
Hrvoje Jasak, Wikki Ltd. All rights reserved
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#include "fvCFD.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
# include "setRootCase.H"
|
||||
|
||||
# include "createTime.H"
|
||||
# include "createMesh.H"
|
||||
# include "createFields.H"
|
||||
# include "initContinuityErrs.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
Info<< "\nStarting time loop\n" << endl;
|
||||
|
||||
while (runTime.loop())
|
||||
{
|
||||
Info<< "Time = " << runTime.timeName() << nl << endl;
|
||||
|
||||
# include "readPISOControls.H"
|
||||
# include "CourantNo.H"
|
||||
|
||||
// Convection-diffusion matrix
|
||||
fvVectorMatrix HUEqn
|
||||
(
|
||||
fvm::div(phi, U)
|
||||
- fvm::laplacian(nu, U)
|
||||
);
|
||||
|
||||
// ddt matrix
|
||||
fvVectorMatrix ddtUEqn(fvm::ddt(U));
|
||||
|
||||
solve(ddtUEqn + HUEqn == -fvc::grad(p));
|
||||
|
||||
// Prepare clean Ap without time derivative contribution
|
||||
// HJ, 26/Oct/2015
|
||||
volScalarField aU = HUEqn.A();
|
||||
|
||||
// --- PISO loop
|
||||
|
||||
for (int corr = 0; corr < nCorr; corr++)
|
||||
{
|
||||
U = HUEqn.H()/aU;
|
||||
phi = (fvc::interpolate(U) & mesh.Sf());
|
||||
|
||||
adjustPhi(phi, U, p);
|
||||
|
||||
for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++)
|
||||
{
|
||||
fvScalarMatrix pEqn
|
||||
(
|
||||
fvm::laplacian(1/aU, p) == fvc::div(phi)
|
||||
);
|
||||
|
||||
pEqn.setReference(pRefCell, pRefValue);
|
||||
pEqn.solve();
|
||||
|
||||
if (nonOrth == nNonOrthCorr)
|
||||
{
|
||||
phi -= pEqn.flux();
|
||||
}
|
||||
}
|
||||
|
||||
# include "continuityErrs.H"
|
||||
|
||||
// Note: cannot call H(U) here because the velocity is not complete
|
||||
// HJ, 22/Jan/2016
|
||||
U = 1.0/(aU + ddtUEqn.A())*
|
||||
(
|
||||
U*aU - fvc::grad(p) + ddtUEqn.H()
|
||||
);
|
||||
U.correctBoundaryConditions();
|
||||
}
|
||||
|
||||
runTime.write();
|
||||
|
||||
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
|
||||
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
|
||||
<< nl << endl;
|
||||
}
|
||||
|
||||
Info<< "End\n" << endl;
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
|
|
@ -0,0 +1,55 @@
|
|||
Info<< "Reading transportProperties\n" << endl;
|
||||
|
||||
IOdictionary transportProperties
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"transportProperties",
|
||||
runTime.constant(),
|
||||
mesh,
|
||||
IOobject::MUST_READ,
|
||||
IOobject::NO_WRITE
|
||||
)
|
||||
);
|
||||
|
||||
dimensionedScalar nu
|
||||
(
|
||||
transportProperties.lookup("nu")
|
||||
);
|
||||
|
||||
Info<< "Reading field p\n" << endl;
|
||||
volScalarField p
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"p",
|
||||
runTime.timeName(),
|
||||
mesh,
|
||||
IOobject::MUST_READ,
|
||||
IOobject::AUTO_WRITE
|
||||
),
|
||||
mesh
|
||||
);
|
||||
|
||||
|
||||
Info<< "Reading field U\n" << endl;
|
||||
volVectorField U
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"U",
|
||||
runTime.timeName(),
|
||||
mesh,
|
||||
IOobject::MUST_READ,
|
||||
IOobject::AUTO_WRITE
|
||||
),
|
||||
mesh
|
||||
);
|
||||
|
||||
|
||||
# include "createPhi.H"
|
||||
|
||||
|
||||
label pRefCell = 0;
|
||||
scalar pRefValue = 0.0;
|
||||
setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
|
||||
|
|
@ -0,0 +1,3 @@
|
|||
consistentPimpleDyMFoam.C
|
||||
|
||||
EXE = $(FOAM_APPBIN)/consistentPimpleDyMFoam
|
||||
|
|
@ -0,0 +1,21 @@
|
|||
EXE_INC = \
|
||||
-I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \
|
||||
-I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
|
||||
-I$(LIB_SRC)/meshTools/lnInclude \
|
||||
-I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
|
||||
-I$(LIB_SRC)/transportModels \
|
||||
-I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
|
||||
-I$(LIB_SRC)/finiteVolume/lnInclude
|
||||
|
||||
EXE_LIBS = \
|
||||
-ldynamicFvMesh \
|
||||
-ltopoChangerFvMesh \
|
||||
-ldynamicMesh \
|
||||
-lmeshTools \
|
||||
-lincompressibleTransportModels \
|
||||
-lincompressibleTurbulenceModel \
|
||||
-lincompressibleRASModels \
|
||||
-lincompressibleLESModels \
|
||||
-lfiniteVolume \
|
||||
-llduSolvers \
|
||||
-L$(MESQUITE_LIB_DIR) -lmesquite
|
||||
|
|
@ -0,0 +1,33 @@
|
|||
// Convection-diffusion matrix
|
||||
fvVectorMatrix HUEqn
|
||||
(
|
||||
fvm::div(phi, U)
|
||||
+ turbulence->divDevReff(U)
|
||||
);
|
||||
|
||||
// ddt matrix
|
||||
fvVectorMatrix ddtUEqn(fvm::ddt(U));
|
||||
|
||||
// Get under-relaxation factor and under-relax the diagonal
|
||||
scalar UUrf = mesh.solutionDict().relaxationFactor(U.name());
|
||||
|
||||
if (oCorr == nOuterCorr - 1)
|
||||
{
|
||||
if (mesh.solutionDict().relax("UFinal"))
|
||||
{
|
||||
UUrf = mesh.solutionDict().relaxationFactor("UFinal");
|
||||
}
|
||||
else
|
||||
{
|
||||
UUrf = 1;
|
||||
}
|
||||
}
|
||||
|
||||
// Solve momentum predictor
|
||||
solve
|
||||
(
|
||||
ddtUEqn
|
||||
+ relax(HUEqn, UUrf)
|
||||
==
|
||||
-fvc::grad(p)
|
||||
);
|
||||
|
|
@ -0,0 +1,128 @@
|
|||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | foam-extend: Open Source CFD
|
||||
\\ / O peration | Version: 3.2
|
||||
\\ / A nd | Web: http://www.foam-extend.org
|
||||
\\/ M anipulation | For copyright notice see file Copyright
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
This file is part of foam-extend.
|
||||
|
||||
foam-extend 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.
|
||||
|
||||
foam-extend 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 foam-extend. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
Application
|
||||
pimpleDyMFoam.C
|
||||
|
||||
Description
|
||||
Transient solver for incompressible, flow of Newtonian fluids
|
||||
with dynamic mesh using the PIMPLE (merged PISO-SIMPLE) algorithm.
|
||||
|
||||
Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#include "fvCFD.H"
|
||||
#include "singlePhaseTransportModel.H"
|
||||
#include "turbulenceModel.H"
|
||||
#include "dynamicFvMesh.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
# include "setRootCase.H"
|
||||
|
||||
# include "createTime.H"
|
||||
# include "createDynamicFvMesh.H"
|
||||
# include "readPIMPLEControls.H"
|
||||
# include "initContinuityErrs.H"
|
||||
# include "createFields.H"
|
||||
# include "readTimeControls.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
Info<< "\nStarting time loop\n" << endl;
|
||||
|
||||
while (runTime.run())
|
||||
{
|
||||
# include "readControls.H"
|
||||
# include "CourantNo.H"
|
||||
# include "setDeltaT.H"
|
||||
|
||||
// Make the fluxes absolute
|
||||
fvc::makeAbsolute(phi, U);
|
||||
|
||||
runTime++;
|
||||
|
||||
Info<< "Time = " << runTime.timeName() << nl << endl;
|
||||
|
||||
bool meshChanged = mesh.update();
|
||||
|
||||
# include "volContinuity.H"
|
||||
|
||||
if (checkMeshCourantNo)
|
||||
{
|
||||
# include "meshCourantNo.H"
|
||||
}
|
||||
|
||||
// Mesh motion update
|
||||
if (correctPhi && meshChanged)
|
||||
{
|
||||
// Fluxes will be corrected to absolute velocity
|
||||
// HJ, 6/Feb/2009
|
||||
# include "correctPhi.H"
|
||||
}
|
||||
|
||||
if (meshChanged)
|
||||
{
|
||||
# include "CourantNo.H"
|
||||
}
|
||||
|
||||
// Make the fluxes relative to the mesh motion
|
||||
fvc::makeRelative(phi, U);
|
||||
|
||||
// --- PIMPLE loop
|
||||
label oCorr = 0;
|
||||
do
|
||||
{
|
||||
if (nOuterCorr != 1)
|
||||
{
|
||||
p.storePrevIter();
|
||||
}
|
||||
|
||||
# include "UEqn.H"
|
||||
|
||||
// --- PISO loop
|
||||
for (int corr = 0; corr < nCorr; corr++)
|
||||
{
|
||||
# include "pEqn.H"
|
||||
}
|
||||
|
||||
turbulence->correct();
|
||||
} while (++oCorr < nOuterCorr);
|
||||
|
||||
runTime.write();
|
||||
|
||||
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
|
||||
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
|
||||
<< nl << endl;
|
||||
}
|
||||
|
||||
Info<< "End\n" << endl;
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
|
|
@ -0,0 +1,58 @@
|
|||
{
|
||||
# include "continuityErrs.H"
|
||||
|
||||
wordList pcorrTypes
|
||||
(
|
||||
p.boundaryField().size(),
|
||||
zeroGradientFvPatchScalarField::typeName
|
||||
);
|
||||
|
||||
for (label i = 0; i < p.boundaryField().size(); i++)
|
||||
{
|
||||
if (p.boundaryField()[i].fixesValue())
|
||||
{
|
||||
pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
|
||||
}
|
||||
}
|
||||
|
||||
volScalarField pcorr
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"pcorr",
|
||||
runTime.timeName(),
|
||||
mesh,
|
||||
IOobject::NO_READ,
|
||||
IOobject::NO_WRITE
|
||||
),
|
||||
mesh,
|
||||
dimensionedScalar("pcorr", p.dimensions(), 0),
|
||||
pcorrTypes
|
||||
);
|
||||
|
||||
// Initialise flux with interpolated velocity
|
||||
phi = fvc::interpolate(U) & mesh.Sf();
|
||||
|
||||
adjustPhi(phi, U, pcorr);
|
||||
|
||||
for(int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++)
|
||||
{
|
||||
fvScalarMatrix pcorrEqn
|
||||
(
|
||||
fvm::laplacian(1/aU, pcorr) == fvc::div(phi)
|
||||
);
|
||||
|
||||
pcorrEqn.setReference(pRefCell, pRefValue);
|
||||
pcorrEqn.solve();
|
||||
|
||||
if (nonOrth == nNonOrthCorr)
|
||||
{
|
||||
phi -= pcorrEqn.flux();
|
||||
}
|
||||
|
||||
// Fluxes are corrected to absolute velocity and further corrected
|
||||
// later. HJ, 6/Feb/2009
|
||||
}
|
||||
|
||||
# include "continuityErrs.H"
|
||||
}
|
||||
|
|
@ -0,0 +1,58 @@
|
|||
Info<< "Reading field p\n" << endl;
|
||||
volScalarField p
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"p",
|
||||
runTime.timeName(),
|
||||
mesh,
|
||||
IOobject::MUST_READ,
|
||||
IOobject::AUTO_WRITE
|
||||
),
|
||||
mesh
|
||||
);
|
||||
|
||||
|
||||
Info<< "Reading field U\n" << endl;
|
||||
volVectorField U
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"U",
|
||||
runTime.timeName(),
|
||||
mesh,
|
||||
IOobject::MUST_READ,
|
||||
IOobject::AUTO_WRITE
|
||||
),
|
||||
mesh
|
||||
);
|
||||
|
||||
# include "createPhi.H"
|
||||
|
||||
|
||||
label pRefCell = 0;
|
||||
scalar pRefValue = 0.0;
|
||||
setRefCell(p, mesh.solutionDict().subDict("PIMPLE"), pRefCell, pRefValue);
|
||||
|
||||
singlePhaseTransportModel laminarTransport(U, phi);
|
||||
|
||||
autoPtr<incompressible::turbulenceModel> turbulence
|
||||
(
|
||||
incompressible::turbulenceModel::New(U, phi, laminarTransport)
|
||||
);
|
||||
|
||||
Info<< "Reading field aU if present\n" << endl;
|
||||
volScalarField aU
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"aU",
|
||||
runTime.timeName(),
|
||||
mesh,
|
||||
IOobject::READ_IF_PRESENT,
|
||||
IOobject::AUTO_WRITE
|
||||
),
|
||||
mesh,
|
||||
1/runTime.deltaT(),
|
||||
zeroGradientFvPatchScalarField::typeName
|
||||
);
|
||||
|
|
@ -0,0 +1,71 @@
|
|||
{
|
||||
p.boundaryField().updateCoeffs();
|
||||
|
||||
// Prepare clean Ap without time derivative contribution and
|
||||
// without contribution from under-relaxation
|
||||
// HJ, 26/Oct/2015
|
||||
aU = HUEqn.A();
|
||||
|
||||
// Store velocity under-relaxation point before using U for
|
||||
// the flux precursor
|
||||
U.storePrevIter();
|
||||
|
||||
U = HUEqn.H()/aU;
|
||||
phi = (fvc::interpolate(U) & mesh.Sf());
|
||||
|
||||
// Global flux balance
|
||||
adjustPhi(phi, U, p);
|
||||
|
||||
for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++)
|
||||
{
|
||||
fvScalarMatrix pEqn
|
||||
(
|
||||
fvm::laplacian(1/aU, p) == fvc::div(phi)
|
||||
);
|
||||
|
||||
pEqn.setReference(pRefCell, pRefValue);
|
||||
|
||||
if
|
||||
(
|
||||
// oCorr == nOuterCorr - 1
|
||||
corr == nCorr - 1
|
||||
&& nonOrth == nNonOrthCorr
|
||||
)
|
||||
{
|
||||
pEqn.solve
|
||||
(
|
||||
mesh.solutionDict().solver(p.name() + "Final")
|
||||
);
|
||||
}
|
||||
else
|
||||
{
|
||||
pEqn.solve(mesh.solutionDict().solver(p.name()));
|
||||
}
|
||||
|
||||
if (nonOrth == nNonOrthCorr)
|
||||
{
|
||||
phi -= pEqn.flux();
|
||||
}
|
||||
}
|
||||
|
||||
// Explicitly relax pressure for momentum corrector
|
||||
if (oCorr != nOuterCorr - 1)
|
||||
{
|
||||
p.relax();
|
||||
}
|
||||
|
||||
// Make the fluxes relative to the mesh motion
|
||||
fvc::makeRelative(phi, U);
|
||||
|
||||
# include "movingMeshContinuityErrs.H"
|
||||
|
||||
U = UUrf*
|
||||
(
|
||||
1.0/(aU + ddtUEqn.A())*
|
||||
(
|
||||
U*aU - fvc::grad(p) + ddtUEqn.H()
|
||||
)
|
||||
)
|
||||
+ (1 - UUrf)*U.prevIter();
|
||||
U.correctBoundaryConditions();
|
||||
}
|
||||
|
|
@ -0,0 +1,14 @@
|
|||
# include "readTimeControls.H"
|
||||
# include "readPIMPLEControls.H"
|
||||
|
||||
bool correctPhi = false;
|
||||
if (pimple.found("correctPhi"))
|
||||
{
|
||||
correctPhi = Switch(pimple.lookup("correctPhi"));
|
||||
}
|
||||
|
||||
bool checkMeshCourantNo = false;
|
||||
if (pimple.found("checkMeshCourantNo"))
|
||||
{
|
||||
checkMeshCourantNo = Switch(pimple.lookup("checkMeshCourantNo"));
|
||||
}
|
||||
|
|
@ -0,0 +1,3 @@
|
|||
consistentSimpleFoam.C
|
||||
|
||||
EXE = $(FOAM_APPBIN)/consistentSimpleFoam
|
||||
|
|
@ -0,0 +1,12 @@
|
|||
EXE_INC = -g \
|
||||
-I$(LIB_SRC)/turbulenceModels \
|
||||
-I$(LIB_SRC)/turbulenceModels/incompressible/RAS/RASModel \
|
||||
-I$(LIB_SRC)/transportModels \
|
||||
-I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
|
||||
-I$(LIB_SRC)/finiteVolume/lnInclude
|
||||
|
||||
EXE_LIBS = \
|
||||
-lincompressibleRASModels \
|
||||
-lincompressibleTransportModels \
|
||||
-lfiniteVolume \
|
||||
-llduSolvers
|
||||
|
|
@ -0,0 +1,20 @@
|
|||
// Solve the momentum equation
|
||||
|
||||
tmp<fvVectorMatrix> HUEqn
|
||||
(
|
||||
fvm::div(phi, U)
|
||||
+ turbulence->divDevReff(U)
|
||||
);
|
||||
|
||||
// Get under-relaxation factor and under-relax the diagonal
|
||||
const scalar UUrf = mesh.solutionDict().relaxationFactor(U.name());
|
||||
|
||||
// Momentum solution
|
||||
eqnResidual = solve
|
||||
(
|
||||
relax(HUEqn(), UUrf)
|
||||
==
|
||||
-fvc::grad(p)
|
||||
).initialResidual();
|
||||
|
||||
maxResidual = max(eqnResidual, maxResidual);
|
||||
|
|
@ -0,0 +1,86 @@
|
|||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | foam-extend: Open Source CFD
|
||||
\\ / O peration | Version: 3.2
|
||||
\\ / A nd | Web: http://www.foam-extend.org
|
||||
\\/ M anipulation | For copyright notice see file Copyright
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
This file is part of foam-extend.
|
||||
|
||||
foam-extend 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.
|
||||
|
||||
foam-extend 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 foam-extend. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
Application
|
||||
consistentSimpleFoam
|
||||
|
||||
Description
|
||||
Steady-state solver for incompressible, turbulent flow
|
||||
Consistent formulation without time-step and relaxation dependence by Jasak
|
||||
|
||||
Author
|
||||
Hrvoje Jasak, Wikki Ltd. All rights reserved
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#include "fvCFD.H"
|
||||
#include "singlePhaseTransportModel.H"
|
||||
#include "RASModel.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
# include "setRootCase.H"
|
||||
# include "createTime.H"
|
||||
# include "createMesh.H"
|
||||
# include "createFields.H"
|
||||
# include "initContinuityErrs.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
Info<< "\nStarting time loop\n" << endl;
|
||||
|
||||
while (runTime.loop())
|
||||
{
|
||||
Info<< "Time = " << runTime.timeName() << nl << endl;
|
||||
|
||||
# include "readSIMPLEControls.H"
|
||||
# include "initConvergenceCheck.H"
|
||||
|
||||
p.storePrevIter();
|
||||
|
||||
// Pressure-velocity SIMPLE corrector
|
||||
{
|
||||
# include "UEqn.H"
|
||||
# include "pEqn.H"
|
||||
}
|
||||
|
||||
turbulence->correct();
|
||||
|
||||
runTime.write();
|
||||
|
||||
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
|
||||
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
|
||||
<< nl << endl;
|
||||
|
||||
# include "convergenceCheck.H"
|
||||
}
|
||||
|
||||
Info<< "End\n" << endl;
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
|
|
@ -0,0 +1,9 @@
|
|||
// check convergence
|
||||
|
||||
if (maxResidual < convergenceCriterion)
|
||||
{
|
||||
Info<< "reached convergence criterion: " << convergenceCriterion << endl;
|
||||
runTime.writeAndEnd();
|
||||
Info<< "latestTime = " << runTime.timeName() << endl;
|
||||
}
|
||||
|
||||
|
|
@ -0,0 +1,42 @@
|
|||
Info << "Reading field p\n" << endl;
|
||||
volScalarField p
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"p",
|
||||
runTime.timeName(),
|
||||
mesh,
|
||||
IOobject::MUST_READ,
|
||||
IOobject::AUTO_WRITE
|
||||
),
|
||||
mesh
|
||||
);
|
||||
|
||||
Info << "Reading field U\n" << endl;
|
||||
volVectorField U
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"U",
|
||||
runTime.timeName(),
|
||||
mesh,
|
||||
IOobject::MUST_READ,
|
||||
IOobject::AUTO_WRITE
|
||||
),
|
||||
mesh
|
||||
);
|
||||
|
||||
# include "createPhi.H"
|
||||
|
||||
|
||||
label pRefCell = 0;
|
||||
scalar pRefValue = 0.0;
|
||||
setRefCell(p, mesh.solutionDict().subDict("SIMPLE"), pRefCell, pRefValue);
|
||||
|
||||
|
||||
singlePhaseTransportModel laminarTransport(U, phi);
|
||||
|
||||
autoPtr<incompressible::RASModel> turbulence
|
||||
(
|
||||
incompressible::RASModel::New(U, phi, laminarTransport)
|
||||
);
|
||||
|
|
@ -0,0 +1,7 @@
|
|||
// initialize values for convergence checks
|
||||
|
||||
scalar eqnResidual = 1, maxResidual = 0;
|
||||
scalar convergenceCriterion = 0;
|
||||
|
||||
simple.readIfPresent("convergence", convergenceCriterion);
|
||||
|
||||
|
|
@ -0,0 +1,59 @@
|
|||
p.boundaryField().updateCoeffs();
|
||||
|
||||
// Prepare clean 1/Ap without contribution from under-relaxation
|
||||
// HJ, 26/Oct/2015
|
||||
volScalarField rUA
|
||||
(
|
||||
"(1|A(U))",
|
||||
1/HUEqn().A()
|
||||
);
|
||||
|
||||
// Store velocity under-relaxation point before using U for
|
||||
// the flux precursor
|
||||
U.storePrevIter();
|
||||
|
||||
U = rUA*HUEqn().H();
|
||||
HUEqn.clear();
|
||||
phi = fvc::interpolate(U) & mesh.Sf();
|
||||
|
||||
// Global flux balance
|
||||
adjustPhi(phi, U, p);
|
||||
|
||||
// Non-orthogonal pressure corrector loop
|
||||
for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++)
|
||||
{
|
||||
fvScalarMatrix pEqn
|
||||
(
|
||||
fvm::laplacian(rUA, p) == fvc::div(phi)
|
||||
);
|
||||
|
||||
pEqn.setReference(pRefCell, pRefValue);
|
||||
|
||||
// Retain the residual from the first iteration
|
||||
if (nonOrth == 0)
|
||||
{
|
||||
eqnResidual = pEqn.solve().initialResidual();
|
||||
maxResidual = max(eqnResidual, maxResidual);
|
||||
}
|
||||
else
|
||||
{
|
||||
pEqn.solve();
|
||||
}
|
||||
|
||||
if (nonOrth == nNonOrthCorr)
|
||||
{
|
||||
phi -= pEqn.flux();
|
||||
}
|
||||
}
|
||||
|
||||
# include "continuityErrs.H"
|
||||
|
||||
// Explicitly relax pressure for momentum corrector
|
||||
p.relax();
|
||||
|
||||
// Momentum corrector
|
||||
// Note: since under-relaxation does not change aU, H/a in U can be
|
||||
// re-used. HJ, 22/Jan/2016
|
||||
U = UUrf*(U - rUA*fvc::grad(p)) + (1 - UUrf)*U.prevIter();
|
||||
U.correctBoundaryConditions();
|
||||
|
||||
Reference in a new issue