zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Merge branch 'master' into dsmc

Browse Source
This commit is contained in:
graham
2009-11-11 13:24:48 +00:00
parent c3eac7b801 3ec8579187
commit 0870e22400
517 changed files with 18508 additions and 4296 deletions
Download Patch File Download Diff File Expand all files Collapse all files

8
applications/solvers/combustion/dieselEngineFoam/createFields.H
View File

@ -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

3
applications/solvers/incompressible/porousSimpleFoam/Make/files Normal file
View File

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

14
applications/solvers/incompressible/porousSimpleFoam/Make/options Normal file
View File

@ -0,0 +1,14 @@
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

47
applications/solvers/incompressible/porousSimpleFoam/UEqn.H Normal file
View File

@ -0,0 +1,47 @@
// 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");
}

64
applications/solvers/incompressible/porousSimpleFoam/createFields.H Normal file
View File

@ -0,0 +1,64 @@
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;
}
}

59
applications/solvers/incompressible/porousSimpleFoam/pEqn.H Normal file
View File

@ -0,0 +1,59 @@
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();

85
applications/solvers/incompressible/porousSimpleFoam/porousSimpleFoam.C Normal file
View File

@ -0,0 +1,85 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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;
}
// ************************************************************************* //

2
applications/solvers/lagrangian/coalChemistryFoam/YEqn.H
View File

@ -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)
{

8
applications/solvers/lagrangian/coalChemistryFoam/createFields.H
View File

@ -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();

2
applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/YEqn.H
View File

@ -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)
{

8
applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createFields.H
View File

@ -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();

2
applications/solvers/lagrangian/reactingParcelFoam/YEqn.H
View File

@ -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)
{

8
applications/solvers/lagrangian/reactingParcelFoam/createFields.H
View File

@ -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();

6
applications/solvers/multiphase/interMixingFoam/Make/files Normal file
View File

@ -0,0 +1,6 @@
incompressibleThreePhaseMixture/threePhaseMixture.C
threePhaseInterfaceProperties/threePhaseInterfaceProperties.C
interMixingFoam.C
EXE = $(FOAM_APPBIN)/interMixingFoam

17
applications/solvers/multiphase/interMixingFoam/Make/options Normal file
View File

@ -0,0 +1,17 @@
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

164
applications/solvers/multiphase/interMixingFoam/alphaEqns.H Normal file
View File

@ -0,0 +1,164 @@
{
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;
}

43
applications/solvers/multiphase/interMixingFoam/alphaEqnsSubCycle.H Normal file
View File

@ -0,0 +1,43 @@
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
{
# include "alphaEqns.H"
}
interface.correct();
{
volScalarField rhoNew = alpha1*rho1 + alpha2*rho2 + alpha3*rho3;
//solve(fvm::ddt(rho) + fvc::div(rhoPhi));
//Info<< "density error = "
// << max((mag(rho - rhoNew)/mag(rhoNew))().internalField()) << endl;
rho == rhoNew;
}

132
applications/solvers/multiphase/interMixingFoam/createFields.H Normal file
View File

@ -0,0 +1,132 @@
Info<< "Reading field p\n" << endl;
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading field alpha1\n" << endl;
volScalarField alpha1
(
IOobject
(
"alpha1",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading field alpha2\n" << endl;
volScalarField alpha2
(
IOobject
(
"alpha2",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading field alpha3\n" << endl;
volScalarField alpha3
(
IOobject
(
"alpha3",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
alpha3 == 1.0 - alpha1 - alpha2;
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
# include "createPhi.H"
threePhaseMixture threePhaseProperties(U, phi);
const dimensionedScalar& rho1 = threePhaseProperties.rho1();
const dimensionedScalar& rho2 = threePhaseProperties.rho2();
const dimensionedScalar& rho3 = threePhaseProperties.rho3();
dimensionedScalar D23(threePhaseProperties.lookup("D23"));
// Need to store rho for ddt(rho, U)
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT
),
alpha1*rho1 + alpha2*rho2 + alpha3*rho3,
alpha1.boundaryField().types()
);
rho.oldTime();
// Mass flux
// Initialisation does not matter because rhoPhi is reset after the
// alpha solution before it is used in the U equation.
surfaceScalarField rhoPhi
(
IOobject
(
"rho*phi",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
rho1*phi
);
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
// Construct interface from alpha distribution
threePhaseInterfaceProperties interface(threePhaseProperties);
// Construct incompressible turbulence model
autoPtr<incompressible::turbulenceModel> turbulence
(
incompressible::turbulenceModel::New(U, phi, threePhaseProperties)
);

204
applications/solvers/multiphase/interMixingFoam/incompressibleThreePhaseMixture/threePhaseMixture.C Normal file
View File

@ -0,0 +1,204 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
Class
threePhaseMixture
\*---------------------------------------------------------------------------*/
#include "threePhaseMixture.H"
#include "addToRunTimeSelectionTable.H"
#include "surfaceFields.H"
#include "fvc.H"
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
//- Calculate and return the laminar viscosity
void Foam::threePhaseMixture::calcNu()
{
// Average kinematic viscosity calculated from dynamic viscosity
nu_ = mu()/(alpha1_*rho1_ + alpha2_*rho2_ + alpha3_*rho3_);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::threePhaseMixture::threePhaseMixture
(
const volVectorField& U,
const surfaceScalarField& phi
)
:
transportModel(U, phi),
phase1Name_("phase1"),
phase2Name_("phase2"),
phase3Name_("phase3"),
nuModel1_
(
viscosityModel::New
(
"nu1",
subDict(phase1Name_),
U,
phi
)
),
nuModel2_
(
viscosityModel::New
(
"nu2",
subDict(phase2Name_),
U,
phi
)
),
nuModel3_
(
viscosityModel::New
(
"nu3",
subDict(phase2Name_),
U,
phi
)
),
rho1_(nuModel1_->viscosityProperties().lookup("rho")),
rho2_(nuModel2_->viscosityProperties().lookup("rho")),
rho3_(nuModel3_->viscosityProperties().lookup("rho")),
U_(U),
phi_(phi),
alpha1_(U_.db().lookupObject<const volScalarField> ("alpha1")),
alpha2_(U_.db().lookupObject<const volScalarField> ("alpha2")),
alpha3_(U_.db().lookupObject<const volScalarField> ("alpha3")),
nu_
(
IOobject
(
"nu",
U_.time().timeName(),
U_.db()
),
U_.mesh(),
dimensionedScalar("nu", dimensionSet(0, 2, -1, 0, 0), 0),
calculatedFvPatchScalarField::typeName
)
{
calcNu();
}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
Foam::tmp<Foam::volScalarField> Foam::threePhaseMixture::mu() const
{
return tmp<volScalarField>
(
new volScalarField
(
"mu",
alpha1_*rho1_*nuModel1_->nu()
+ alpha2_*rho2_*nuModel2_->nu()
+ alpha3_*rho3_*nuModel3_->nu()
)
);
}
Foam::tmp<Foam::surfaceScalarField> Foam::threePhaseMixture::muf() const
{
surfaceScalarField alpha1f = fvc::interpolate(alpha1_);
surfaceScalarField alpha2f = fvc::interpolate(alpha2_);
surfaceScalarField alpha3f = fvc::interpolate(alpha3_);
return tmp<surfaceScalarField>
(
new surfaceScalarField
(
"mu",
alpha1f*rho1_*fvc::interpolate(nuModel1_->nu())
+ alpha2f*rho2_*fvc::interpolate(nuModel2_->nu())
+ alpha3f*rho3_*fvc::interpolate(nuModel3_->nu())
)
);
}
Foam::tmp<Foam::surfaceScalarField> Foam::threePhaseMixture::nuf() const
{
surfaceScalarField alpha1f = fvc::interpolate(alpha1_);
surfaceScalarField alpha2f = fvc::interpolate(alpha2_);
surfaceScalarField alpha3f = fvc::interpolate(alpha3_);
return tmp<surfaceScalarField>
(
new surfaceScalarField
(
"nu",
(
alpha1f*rho1_*fvc::interpolate(nuModel1_->nu())
+ alpha2f*rho2_*fvc::interpolate(nuModel2_->nu())
+ alpha3f*rho3_*fvc::interpolate(nuModel3_->nu())
)/(alpha1f*rho1_ + alpha2f*rho2_ + alpha3f*rho3_)
)
);
}
bool Foam::threePhaseMixture::read()
{
if (transportModel::read())
{
if
(
nuModel1_().read(*this)
&& nuModel2_().read(*this)
&& nuModel3_().read(*this)
)
{
nuModel1_->viscosityProperties().lookup("rho") >> rho1_;
nuModel2_->viscosityProperties().lookup("rho") >> rho2_;
nuModel3_->viscosityProperties().lookup("rho") >> rho3_;
return true;
}
else
{
return false;
}
}
else
{
return false;
}
}
// ************************************************************************* //

197
applications/solvers/multiphase/interMixingFoam/incompressibleThreePhaseMixture/threePhaseMixture.H Normal file
View File

@ -0,0 +1,197 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
Class
threePhaseMixture
Description
SourceFiles
threePhaseMixture.C
\*---------------------------------------------------------------------------*/
#ifndef threePhaseMixture_H
#define threePhaseMixture_H
#include "incompressible/transportModel/transportModel.H"
#include "incompressible/viscosityModels/viscosityModel/viscosityModel.H"
#include "dimensionedScalar.H"
#include "volFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class threePhaseMixture Declaration
\*---------------------------------------------------------------------------*/
class threePhaseMixture
:
public transportModel
{
// Private data
word phase1Name_;
word phase2Name_;
word phase3Name_;
autoPtr<viscosityModel> nuModel1_;
autoPtr<viscosityModel> nuModel2_;
autoPtr<viscosityModel> nuModel3_;
dimensionedScalar rho1_;
dimensionedScalar rho2_;
dimensionedScalar rho3_;
const volVectorField& U_;
const surfaceScalarField& phi_;
const volScalarField& alpha1_;
const volScalarField& alpha2_;
const volScalarField& alpha3_;
volScalarField nu_;
// Private Member Functions
//- Calculate and return the laminar viscosity
void calcNu();
public:
// Constructors
//- Construct from components
threePhaseMixture
(
const volVectorField& U,
const surfaceScalarField& phi
);
// Destructor
~threePhaseMixture()
{}
// Member Functions
//- Return const-access to phase1 viscosityModel
const viscosityModel& nuModel1() const
{
return nuModel1_();
}
//- Return const-access to phase2 viscosityModel
const viscosityModel& nuModel2() const
{
return nuModel2_();
}
//- Return const-access to phase3 viscosityModel
const viscosityModel& nuModel3() const
{
return nuModel3_();
}
//- Return const-access to phase1 density
const dimensionedScalar& rho1() const
{
return rho1_;
}
//- Return const-access to phase2 density
const dimensionedScalar& rho2() const
{
return rho2_;
};
//- Return const-access to phase3 density
const dimensionedScalar& rho3() const
{
return rho3_;
};
const volScalarField& alpha1() const
{
return alpha1_;
}
const volScalarField& alpha2() const
{
return alpha2_;
}
const volScalarField& alpha3() const
{
return alpha3_;
}
//- Return the velocity
const volVectorField& U() const
{
return U_;
}
//- Return the dynamic laminar viscosity
tmp<volScalarField> mu() const;
//- Return the face-interpolated dynamic laminar viscosity
tmp<surfaceScalarField> muf() const;
//- Return the kinematic laminar viscosity
tmp<volScalarField> nu() const
{
return nu_;
}
//- Return the face-interpolated dynamic laminar viscosity
tmp<surfaceScalarField> nuf() const;
//- Correct the laminar viscosity
void correct()
{
calcNu();
}
//- Read base transportProperties dictionary
bool read();
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

102
applications/solvers/multiphase/interMixingFoam/interMixingFoam.C Normal file
View File

@ -0,0 +1,102 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
Application
interMixingFoam
Description
Solver for 3 incompressible fluids, two of which are miscible,
using a VOF method to capture the interface.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "MULES.H"
#include "subCycle.H"
#include "threePhaseMixture.H"
#include "threePhaseInterfaceProperties.H"
#include "turbulenceModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "readGravitationalAcceleration.H"
#include "readPISOControls.H"
#include "initContinuityErrs.H"
#include "createFields.H"
#include "readTimeControls.H"
#include "CourantNo.H"
#include "setInitialDeltaT.H"
#include "correctPhi.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.run())
{
#include "readPISOControls.H"
#include "readTimeControls.H"
#include "CourantNo.H"
#include "setDeltaT.H"
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
threePhaseProperties.correct();
#include "alphaEqnsSubCycle.H"
#define twoPhaseProperties threePhaseProperties
#include "UEqn.H"
// --- PISO loop
for (int corr=0; corr<nCorr; corr++)
{
#include "pEqn.H"
}
#include "continuityErrs.H"
turbulence->correct();
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "\n end \n";
return(0);
}
// ************************************************************************* //

209
applications/solvers/multiphase/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.C Normal file
View File

@ -0,0 +1,209 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
Application
threePhaseInterfaceProperties
Description
Properties to aid interFoam :
1. Correct the alpha boundary condition for dynamic contact angle.
2. Calculate interface curvature.
\*---------------------------------------------------------------------------*/
#include "threePhaseInterfaceProperties.H"
#include "alphaContactAngleFvPatchScalarField.H"
#include "mathematicalConstants.H"
#include "surfaceInterpolate.H"
#include "fvcDiv.H"
#include "fvcGrad.H"
// * * * * * * * * * * * * * * * Static Member Data * * * * * * * * * * * * //
const Foam::scalar Foam::threePhaseInterfaceProperties::convertToRad =
Foam::constant::mathematical::pi/180.0;
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// Correction for the boundary condition on the unit normal nHat on
// walls to produce the correct contact angle.
// The dynamic contact angle is calculated from the component of the
// velocity on the direction of the interface, parallel to the wall.
void Foam::threePhaseInterfaceProperties::correctContactAngle
(
surfaceVectorField::GeometricBoundaryField& nHatb
) const
{
const volScalarField::GeometricBoundaryField& alpha1 =
mixture_.alpha1().boundaryField();
const volScalarField::GeometricBoundaryField& alpha2 =
mixture_.alpha2().boundaryField();
const volScalarField::GeometricBoundaryField& alpha3 =
mixture_.alpha3().boundaryField();
const volVectorField::GeometricBoundaryField& U =
mixture_.U().boundaryField();
const fvMesh& mesh = mixture_.U().mesh();
const fvBoundaryMesh& boundary = mesh.boundary();
forAll(boundary, patchi)
{
if (isA<alphaContactAngleFvPatchScalarField>(alpha1[patchi]))
{
const alphaContactAngleFvPatchScalarField& a2cap =
refCast<const alphaContactAngleFvPatchScalarField>
(alpha2[patchi]);
const alphaContactAngleFvPatchScalarField& a3cap =
refCast<const alphaContactAngleFvPatchScalarField>
(alpha3[patchi]);
scalarField twoPhaseAlpha2 = max(a2cap, scalar(0));
scalarField twoPhaseAlpha3 = max(a3cap, scalar(0));
scalarField sumTwoPhaseAlpha =
twoPhaseAlpha2 + twoPhaseAlpha3 + SMALL;
twoPhaseAlpha2 /= sumTwoPhaseAlpha;
twoPhaseAlpha3 /= sumTwoPhaseAlpha;
fvsPatchVectorField& nHatp = nHatb[patchi];
scalarField theta =
convertToRad
*(
twoPhaseAlpha2*(180 - a2cap.theta(U[patchi], nHatp))
+ twoPhaseAlpha3*(180 - a3cap.theta(U[patchi], nHatp))
);
vectorField nf = boundary[patchi].nf();
// Reset nHatPatch to correspond to the contact angle
scalarField a12 = nHatp & nf;
scalarField b1 = cos(theta);
scalarField b2(nHatp.size());
forAll(b2, facei)
{
b2[facei] = cos(acos(a12[facei]) - theta[facei]);
}
scalarField det = 1.0 - a12*a12;
scalarField a = (b1 - a12*b2)/det;
scalarField b = (b2 - a12*b1)/det;
nHatp = a*nf + b*nHatp;
nHatp /= (mag(nHatp) + deltaN_.value());
}
}
}
void Foam::threePhaseInterfaceProperties::calculateK()
{
const volScalarField& alpha1 = mixture_.alpha1();
const fvMesh& mesh = alpha1.mesh();
const surfaceVectorField& Sf = mesh.Sf();
// Cell gradient of alpha
volVectorField gradAlpha = fvc::grad(alpha1);
// Interpolated face-gradient of alpha
surfaceVectorField gradAlphaf = fvc::interpolate(gradAlpha);
// Face unit interface normal
surfaceVectorField nHatfv = gradAlphaf/(mag(gradAlphaf) + deltaN_);
correctContactAngle(nHatfv.boundaryField());
// Face unit interface normal flux
nHatf_ = nHatfv & Sf;
// Simple expression for curvature
K_ = -fvc::div(nHatf_);
// Complex expression for curvature.
// Correction is formally zero but numerically non-zero.
//volVectorField nHat = gradAlpha/(mag(gradAlpha) + deltaN_);
//nHat.boundaryField() = nHatfv.boundaryField();
//K_ = -fvc::div(nHatf_) + (nHat & fvc::grad(nHatfv) & nHat);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::threePhaseInterfaceProperties::threePhaseInterfaceProperties
(
const threePhaseMixture& mixture
)
:
mixture_(mixture),
cAlpha_
(
readScalar
(
mixture.U().mesh().solutionDict().subDict("PISO").lookup("cAlpha")
)
),
sigma12_(mixture.lookup("sigma12")),
sigma13_(mixture.lookup("sigma13")),
deltaN_
(
"deltaN",
1e-8/pow(average(mixture.U().mesh().V()), 1.0/3.0)
),
nHatf_
(
(
fvc::interpolate(fvc::grad(mixture.alpha1()))
/(mag(fvc::interpolate(fvc::grad(mixture.alpha1()))) + deltaN_)
) & mixture.alpha1().mesh().Sf()
),
K_
(
IOobject
(
"K",
mixture.alpha1().time().timeName(),
mixture.alpha1().mesh()
),
-fvc::div(nHatf_)
)
{
calculateK();
}
// ************************************************************************* //

157
applications/solvers/multiphase/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.H Normal file
View File

@ -0,0 +1,157 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
Class
threePhaseInterfaceProperties
Description
Properties to aid interFoam :
1. Correct the alpha boundary condition for dynamic contact angle.
2. Calculate interface curvature.
SourceFiles
threePhaseInterfaceProperties.C
\*---------------------------------------------------------------------------*/
#ifndef threePhaseInterfaceProperties_H
#define threePhaseInterfaceProperties_H
#include "threePhaseMixture.H"
#include "surfaceFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class threePhaseInterfaceProperties Declaration
\*---------------------------------------------------------------------------*/
class threePhaseInterfaceProperties
{
// Private data
const threePhaseMixture& mixture_;
//- Compression coefficient
scalar cAlpha_;
//- Surface tension 1-2
dimensionedScalar sigma12_;
//- Surface tension 1-3
dimensionedScalar sigma13_;
//- Stabilisation for normalisation of the interface normal
const dimensionedScalar deltaN_;
surfaceScalarField nHatf_;
volScalarField K_;
// Private Member Functions
//- Disallow default bitwise copy construct and assignment
threePhaseInterfaceProperties(const threePhaseInterfaceProperties&);
void operator=(const threePhaseInterfaceProperties&);
//- Correction for the boundary condition on the unit normal nHat on
// walls to produce the correct contact dynamic angle
// calculated from the component of U parallel to the wall
void correctContactAngle
(
surfaceVectorField::GeometricBoundaryField& nHat
) const;
//- Re-calculate the interface curvature
void calculateK();
public:
//- Conversion factor for degrees into radians
static const scalar convertToRad;
// Constructors
//- Construct from volume fraction field alpha and IOdictionary
threePhaseInterfaceProperties(const threePhaseMixture& mixture);
// Member Functions
scalar cAlpha() const
{
return cAlpha_;
}
const dimensionedScalar& deltaN() const
{
return deltaN_;
}
const surfaceScalarField& nHatf() const
{
return nHatf_;
}
const volScalarField& K() const
{
return K_;
}
tmp<volScalarField> sigma() const
{
volScalarField limitedAlpha2 = max(mixture_.alpha2(), scalar(0));
volScalarField limitedAlpha3 = max(mixture_.alpha3(), scalar(0));
return
(limitedAlpha2*sigma12_ + limitedAlpha3*sigma13_)
/(limitedAlpha2 + limitedAlpha3 + SMALL);
}
tmp<volScalarField> sigmaK() const
{
return sigma()*K_;
}
void correct()
{
calculateK();
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //