exaFOAM/foam-extend4.1-coherent-io
Archived
1
1
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

Updated for dynamic meshes; improved numerics

Browse source
This commit is contained in:
Hrvoje Jasak 2014-11-27 19:09:11 +00:00
parent 5305c8f4ae
commit 953b9dea77
11 changed files with 439 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

3
applications/solvers/compressible/sonicDyMFoam/Make/files Normal file
View file

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

21
applications/solvers/compressible/sonicDyMFoam/Make/options Normal file
View file

@ -0,0 +1,21 @@
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
-I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \
-I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude
EXE_LIBS = \
-lbasicThermophysicalModels \
-lspecie \
-lcompressibleTurbulenceModel \
-lcompressibleRASModels \
-lcompressibleLESModels \
-lfiniteVolume \
-ldynamicFvMesh \
-ltopoChangerFvMesh \
-ldynamicMesh \
-lmeshTools \
-llduSolvers

24
applications/solvers/compressible/sonicDyMFoam/UEqn.H Normal file
View file

@ -0,0 +1,24 @@
fvVectorMatrix UEqn
(
fvm::ddt(rho, U)
+ fvm::div(phi, U)
+ turbulence->divDevRhoReff(U)
);
if (oCorr == nOuterCorr - 1)
{
if (mesh.solutionDict().relax("UFinal"))
{
UEqn.relax(mesh.solutionDict().relaxationFactor("UFinal"));
}
else
{
UEqn.relax(1);
}
}
else
{
UEqn.relax();
}
solve(UEqn == -fvc::grad(p));

20
applications/solvers/compressible/sonicDyMFoam/compressibleContinuityErrs.H Normal file
View file

@ -0,0 +1,20 @@
{
# include "rhoEqn.H"
}
{
scalar sumLocalContErr =
sum
(
mag(rho.internalField() - (psi*p)().internalField())
)/sum(rho.internalField());
scalar globalContErr =
sum(rho.internalField() - (psi*p)().internalField())
/sum(rho.internalField());
cumulativeContErr += globalContErr;
Info<< "time step continuity errors : sum local = " << sumLocalContErr
<< ", global = " << globalContErr
<< ", cumulative = " << cumulativeContErr << endl;
}

64
applications/solvers/compressible/sonicDyMFoam/createFields.H Normal file
View file

@ -0,0 +1,64 @@
Info<< "Reading thermophysical properties\n" << endl;
autoPtr<basicPsiThermo> pThermo
(
basicPsiThermo::New(mesh)
);
basicPsiThermo& thermo = pThermo();
volScalarField& p = thermo.p();
volScalarField& e = thermo.e();
const volScalarField& psi = thermo.psi();
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh
),
thermo.rho()
);
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
# include "compressibleCreatePhi.H"
// Store energy source term for under-relaxation
volScalarField pDivU
(
IOobject
(
"pDivU",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
p*fvc::div(phi/fvc::interpolate(rho))
);
Info<< "Creating turbulence model\n" << endl;
autoPtr<compressible::turbulenceModel> turbulence
(
compressible::turbulenceModel::New
(
rho,
U,
phi,
thermo
)
);

44
applications/solvers/compressible/sonicDyMFoam/eEqn.H Normal file
View file

@ -0,0 +1,44 @@
{
fvScalarMatrix eEqn
(
fvm::ddt(rho, e)
+ fvm::div(phi, e)
- fvm::laplacian(turbulence->alphaEff(), e)
==
// - p*fvc::div(phi/fvc::interpolate(rho) + fvc::meshPhi(rho, U))
- fvm::SuSp(pDivU/e, e)
// viscous heating?
);
if (oCorr == nOuterCorr - 1)
{
if (mesh.solutionDict().relax("eFinal"))
{
eEqn.relax(mesh.solutionDict().relaxationFactor("eFinal"));
}
else
{
eEqn.relax(1);
}
}
else
{
eEqn.relax();
}
eEqn.solve();
// Bound the energy using TMin and TMax
{
dimensionedScalar Tstd("Tstd", dimTemperature, specie::Tstd);
volScalarField Cv = thermo.Cv();
volScalarField R = thermo.Cp() - Cv;
e = Foam::min(e, TMax*Cv + R*Tstd);
e = Foam::max(e, TMin*Cv + R*Tstd);
e.correctBoundaryConditions();
}
thermo.correct();
}

17
applications/solvers/compressible/sonicDyMFoam/limitU.H Normal file
View file

@ -0,0 +1,17 @@
{
// Bound the velocity
volScalarField magU = mag(U);
if (max(magU) > UMax)
{
volScalarField Ulimiter =
pos(magU - UMax)*UMax/(magU + smallU)
+ neg(magU - UMax);
Ulimiter.max(scalar(0));
Ulimiter.min(scalar(1));
U *= Ulimiter;
U.correctBoundaryConditions();
}
}

14
applications/solvers/compressible/sonicDyMFoam/readControls.H Normal file
View file

@ -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"));
}

20
applications/solvers/compressible/sonicDyMFoam/readFieldBounds.H Normal file
View file

@ -0,0 +1,20 @@
// Read field bounds
dictionary fieldBounds = mesh.solutionDict().subDict("fieldBounds");
// Pressure bounds
dimensionedScalar pMin("pMin", dimPressure, 0);
dimensionedScalar pMax("pMax", dimPressure, GREAT);
fieldBounds.lookup("p") >> pMin.value() >> pMax.value();
// Temperature bounds
dimensionedScalar TMin("TMin", dimTemperature, 0);
dimensionedScalar TMax("TMax", dimTemperature, GREAT);
fieldBounds.lookup("T") >> TMin.value() >> TMax.value();
// Velocity bound
dimensionedScalar UMax("UMax", dimVelocity, GREAT);
fieldBounds.lookup(U.name()) >> UMax.value();
dimensionedScalar smallU("smallU", dimVelocity, 1e-10);

18
applications/solvers/compressible/sonicDyMFoam/readTransportProperties.H Normal file
View file

@ -0,0 +1,18 @@
Info<< "Reading transportProperties\n" << endl;
IOdictionary transportProperties
(
IOobject
(
"transportProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
)
);
dimensionedScalar mu
(
transportProperties.lookup("mu")
);

194
applications/solvers/compressible/sonicDyMFoam/sonicDyMFoam.C Normal file
View file

@ -0,0 +1,194 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
-------------------------------------------------------------------------------
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
sonicDyMFoam
Description
Transient solver for trans-sonic/supersonic, laminar flow of a
compressible gas with support for mesh motion and topological changes
Updated from sonicFoamAutoMotion by Hrvoje Jasak
Author
Hrvoje Jasak, Wikki Ltd. All rights reserved.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "dynamicFvMesh.H"
#include "specie.H"
#include "basicPsiThermo.H"
#include "turbulenceModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
# include "setRootCase.H"
# include "createTime.H"
# include "createDynamicFvMesh.H"
# include "createFields.H"
# include "initContinuityErrs.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.run())
{
# include "readControls.H"
# include "readFieldBounds.H"
# include "compressibleCourantNo.H"
# include "setDeltaT.H"
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)
// {
// # include "correctPhi.H"
// }
if (meshChanged)
{
# include "CourantNo.H"
}
// --- PIMPLE loop
label oCorr = 0;
do
{
// Under-relax pDivU term
pDivU.storePrevIter();
pDivU = p*fvc::div(phi/fvc::interpolate(rho));
pDivU.relax();
# include "rhoEqn.H"
# include "eEqn.H"
# include "UEqn.H"
// --- PISO loop
for (int corr = 0; corr < nCorr; corr++)
{
U = UEqn.H()/UEqn.A();
# include "limitU.H"
surfaceScalarField phid
(
"phid",
fvc::interpolate(psi)*
(
(fvc::interpolate(U) & mesh.Sf())
- fvc::meshPhi(rho, U)
)
);
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
// Store pressure for under-relaxation
p.storePrevIter();
fvScalarMatrix pEqn
(
fvm::ddt(psi, p)
+ fvm::div(phid, p)
- fvm::laplacian(rho/UEqn.A(), p)
);
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();
}
// Bound the pressure
if (min(p) < pMin || max(p) > pMax)
{
p.max(pMin);
p.min(pMax);
p.correctBoundaryConditions();
}
// Relax the pressure
p.relax();
}
# include "compressibleContinuityErrs.H"
U -= fvc::grad(p)/UEqn.A();
U.correctBoundaryConditions();
# include "limitU.H"
}
// Recalculate density
rho = thermo.rho();
turbulence->correct();
} while (++oCorr < nOuterCorr);
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return(0);
}
// ************************************************************************* //