Commit 6930293d authored by mattijs's avatar mattijs
Browse files

Merge branch 'master' of /home/noisy3/OpenFOAM/OpenFOAM-dev

parents afeb6029 51bb18ce
......@@ -13,6 +13,14 @@ PtrList<volScalarField>& Y = composition.Y();
word inertSpecie(thermo.lookup("inertSpecie"));
if (!composition.contains(inertSpecie))
{
FatalErrorIn(args.executable())
<< "Specified inert specie '" << inertSpecie << "' not found in "
<< "species list. Available species:" << composition.species()
<< exit(FatalError);
}
volScalarField rho
(
IOobject
......
porousSimpleFoam.C
EXE = $(FOAM_APPBIN)/porousSimpleFoam
EXE_INC = \
-I../simpleFoam \
-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 \
-I$(LIB_SRC)/meshTools/lnInclude
EXE_LIBS = \
-lincompressibleRASModels \
-lincompressibleTransportModels \
-lfiniteVolume \
-lmeshTools
// Construct the Momentum equation
tmp<fvVectorMatrix> UEqn
(
fvm::div(phi, U)
- fvm::Sp(fvc::div(phi), U)
+ turbulence->divDevReff(U)
);
UEqn().relax();
// Include the porous media resistance and solve the momentum equation
// either implicit in the tensorial resistance or transport using by
// including the spherical part of the resistance in the momentum diagonal
tmp<volScalarField> trAU;
tmp<volTensorField> trTU;
if (pressureImplicitPorosity)
{
tmp<volTensorField> tTU = tensor(I)*UEqn().A();
pZones.addResistance(UEqn(), tTU());
trTU = inv(tTU());
trTU().rename("rAU");
volVectorField gradp = fvc::grad(p);
for (int UCorr=0; UCorr<nUCorr; UCorr++)
{
U = trTU() & (UEqn().H() - gradp);
}
U.correctBoundaryConditions();
}
else
{
pZones.addResistance(UEqn());
eqnResidual = solve
(
UEqn() == -fvc::grad(p)
). initialResidual();
maxResidual = max(eqnResidual, maxResidual);
trAU = 1.0/UEqn().A();
trAU().rename("rAU");
}
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)
);
porousZones pZones(mesh);
Switch pressureImplicitPorosity(false);
int nUCorr = 0;
if (pZones.size())
{
// nUCorrectors for pressureImplicitPorosity
if (mesh.solutionDict().subDict("SIMPLE").found("nUCorrectors"))
{
nUCorr = readInt
(
mesh.solutionDict().subDict("SIMPLE").lookup("nUCorrectors")
);
}
if (nUCorr > 0)
{
pressureImplicitPorosity = true;
}
}
if (pressureImplicitPorosity)
{
U = trTU()&UEqn().H();
}
else
{
U = trAU()*UEqn().H();
}
UEqn.clear();
phi = fvc::interpolate(U) & mesh.Sf();
adjustPhi(phi, U, p);
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
tmp<fvScalarMatrix> tpEqn;
if (pressureImplicitPorosity)
{
tpEqn = (fvm::laplacian(trTU(), p) == fvc::div(phi));
}
else
{
tpEqn = (fvm::laplacian(trAU(), p) == fvc::div(phi));
}
tpEqn().setReference(pRefCell, pRefValue);
// retain the residual from the first iteration
if (nonOrth == 0)
{
eqnResidual = tpEqn().solve().initialResidual();
maxResidual = max(eqnResidual, maxResidual);
}
else
{
tpEqn().solve();
}
if (nonOrth == nNonOrthCorr)
{
phi -= tpEqn().flux();
}
}
#include "continuityErrs.H"
// Explicitly relax pressure for momentum corrector
p.relax();
if (pressureImplicitPorosity)
{
U -= trTU()&fvc::grad(p);
}
else
{
U -= trAU()*fvc::grad(p);
}
U.correctBoundaryConditions();
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Application
porousSimpleFoam
Description
Steady-state solver for incompressible, turbulent flow with
implicit or explicit porosity treatment
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "singlePhaseTransportModel.H"
#include "RASModel.H"
#include "porousZones.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;
}
// ************************************************************************* //
......@@ -14,7 +14,7 @@ tmp<fv::convectionScheme<scalar> > mvConvection
label inertIndex = -1;
volScalarField Yt = 0.0*Y[0];
for (label i=0; i<Y.size(); i++)
forAll(Y, i)
{
if (Y[i].name() != inertSpecie)
{
......
......@@ -13,6 +13,14 @@
word inertSpecie(thermo.lookup("inertSpecie"));
if (!composition.contains(inertSpecie))
{
FatalErrorIn(args.executable())
<< "Specified inert specie '" << inertSpecie << "' not found in "
<< "species list. Available species:" << composition.species()
<< exit(FatalError);
}
volScalarField& p = thermo.p();
volScalarField& h = thermo.h();
const volScalarField& T = thermo.T();
......
......@@ -15,7 +15,7 @@ tmp<fv::convectionScheme<scalar> > mvConvection
label inertIndex = -1;
volScalarField Yt = 0.0*Y[0];
for (label i=0; i<Y.size(); i++)
forAll(Y, i)
{
if (Y[i].name() != inertSpecie)
{
......
......@@ -13,6 +13,14 @@
word inertSpecie(thermo.lookup("inertSpecie"));
if (!composition.contains(inertSpecie))
{
FatalErrorIn(args.executable())
<< "Specified inert specie '" << inertSpecie << "' not found in "
<< "species list. Available species:" << composition.species()
<< exit(FatalError);
}
volScalarField& p = thermo.p();
volScalarField& h = thermo.h();
const volScalarField& T = thermo.T();
......
......@@ -14,7 +14,7 @@ tmp<fv::convectionScheme<scalar> > mvConvection
label inertIndex = -1;
volScalarField Yt = 0.0*Y[0];
for (label i=0; i<Y.size(); i++)
forAll(Y, i)
{
if (Y[i].name() != inertSpecie)
{
......
......@@ -13,6 +13,14 @@
word inertSpecie(thermo.lookup("inertSpecie"));
if (!composition.contains(inertSpecie))
{
FatalErrorIn(args.executable())
<< "Specified inert specie '" << inertSpecie << "' not found in "
<< "species list. Available species:" << composition.species()
<< exit(FatalError);
}
volScalarField& p = thermo.p();
volScalarField& h = thermo.h();
const volScalarField& T = thermo.T();
......
incompressibleThreePhaseMixture/threePhaseMixture.C
threePhaseInterfaceProperties/threePhaseInterfaceProperties.C
interMixingFoam.C
EXE = $(FOAM_APPBIN)/interMixingFoam
INTERFOAM = $(FOAM_SOLVERS)/multiphase/interFoam
EXE_INC = \
-I$(INTERFOAM) \
-IincompressibleThreePhaseMixture \
-IthreePhaseInterfaceProperties \
-I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
-I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/transportModels
EXE_LIBS = \
-linterfaceProperties \
-lincompressibleTransportModels \
-lincompressibleRASModels \
-lincompressibleLESModels \
-lfiniteVolume
{
word alphaScheme("div(phi,alpha)");
word alpharScheme("div(phirb,alpha)");
surfaceScalarField phir
(
IOobject
(
"phir",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
interface.cAlpha()*mag(phi/mesh.magSf())*interface.nHatf()
);
for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)
{
// Create the limiter to be used for all phase-fractions
scalarField allLambda(mesh.nFaces(), 1.0);
// Split the limiter into a surfaceScalarField
slicedSurfaceScalarField lambda
(
IOobject
(
"lambda",
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
mesh,
dimless,
allLambda
);
// Create the complete convection flux for alpha1
surfaceScalarField phiAlpha1 =
fvc::flux
(
phi,
alpha1,
alphaScheme
)
+ fvc::flux
(
-fvc::flux(-phir, alpha2, alpharScheme),
alpha1,
alpharScheme
)
+ fvc::flux
(
-fvc::flux(-phir, alpha3, alpharScheme),
alpha1,
alpharScheme
);
// Create the bounded (upwind) flux for alpha1
surfaceScalarField phiAlpha1BD =
upwind<scalar>(mesh, phi).flux(alpha1);
// Calculate the flux correction for alpha1
phiAlpha1 -= phiAlpha1BD;
// Calculate the limiter for alpha1
MULES::limiter
(
allLambda,
oneField(),
alpha1,
phiAlpha1BD,
phiAlpha1,
zeroField(),
zeroField(),
1,
0,
3
);
// Create the complete flux for alpha2
surfaceScalarField phiAlpha2 =
fvc::flux
(
phi,
alpha2,
alphaScheme
)
+ fvc::flux
(
-fvc::flux(phir, alpha1, alpharScheme),
alpha2,
alpharScheme
);
// Create the bounded (upwind) flux for alpha2
surfaceScalarField phiAlpha2BD =
upwind<scalar>(mesh, phi).flux(alpha2);
// Calculate the flux correction for alpha2
phiAlpha2 -= phiAlpha2BD;
// Further limit the limiter for alpha2
MULES::limiter
(
allLambda,
oneField(),
alpha2,
phiAlpha2BD,
phiAlpha2,
zeroField(),
zeroField(),
1,
0,
3
);
// Construct the limited fluxes
phiAlpha1 = phiAlpha1BD + lambda*phiAlpha1;
phiAlpha2 = phiAlpha2BD + lambda*phiAlpha2;
// Solve for alpha1
solve(fvm::ddt(alpha1) + fvc::div(phiAlpha1));
// Create the diffusion coefficients for alpha2<->alpha3
volScalarField Dc23 = D23*max(alpha3, scalar(0))*pos(alpha2);
volScalarField Dc32 = D23*max(alpha2, scalar(0))*pos(alpha3);
// Add the diffusive flux for alpha3->alpha2
phiAlpha2 -= fvc::interpolate(Dc32)*mesh.magSf()*fvc::snGrad(alpha1);
// Solve for alpha2
fvScalarMatrix alpha2Eqn
(
fvm::ddt(alpha2)
+ fvc::div(phiAlpha2)
- fvm::laplacian(Dc23 + Dc32, alpha2)
);
alpha2Eqn.solve();
// Construct the complete mass flux
rhoPhi =
phiAlpha1*(rho1 - rho3)
+ (phiAlpha2 + alpha2Eqn.flux())*(rho2 - rho3)
+ phi*rho3;
alpha3 = 1.0 - alpha1 - alpha2;
}
Info<< "Air phase volume fraction = "
<< alpha1.weightedAverage(mesh.V()).value()
<< " Min(alpha1) = " << min(alpha1).value()
<< " Max(alpha1) = " << max(alpha1).value()
<< endl;
Info<< "Liquid phase volume fraction = "
<< alpha2.weightedAverage(mesh.V()).value()
<< " Min(alpha2) = " << min(alpha2).value()
<< " Max(alpha2) = " << max(alpha2).value()
<< endl;
}
label nAlphaCorr
(
readLabel(piso.lookup("nAlphaCorr"))
);
label nAlphaSubCycles
(
readLabel(piso.lookup("nAlphaSubCycles"))
);
if (nAlphaSubCycles > 1)
{
surfaceScalarField rhoPhiSum = 0.0*rhoPhi;
dimensionedScalar totalDeltaT = runTime.deltaT();
for
(
subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
!(++alphaSubCycle).end();
)
{
# include "alphaEqns.H"
rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
}
rhoPhi = rhoPhiSum;
}
else