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

Added generic turbulenceModel base class to incompressible turbulence models.

Browse Source
This commit is contained in:
henry
2008-11-20 20:33:06 +00:00
parent 326b86ec2d
commit 0479165024
307 changed files with 64527 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

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

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

@ -0,0 +1,14 @@
EXE_INC = \
-I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
-I$(LIB_SRC)/turbulenceModels/incompressible/LES/LESModel \
-I$(LIB_SRC)/turbulenceModels/LES/LESdeltas/lnInclude \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude
EXE_LIBS = \
-lincompressibleLESModels \
-lincompressibleTransportModels \
-lfiniteVolume \
-lmeshTools

155
applications/solvers/incompressible/channelFoam/channelFoam.C Normal file
View File

@ -0,0 +1,155 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 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
channelFoam
Description
Incompressible LES solver for flow in a channel.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "singlePhaseTransportModel.H"
#include "LESModel.H"
#include "IFstream.H"
#include "OFstream.H"
#include "Random.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "readTransportProperties.H"
#include "createFields.H"
#include "initContinuityErrs.H"
#include "createGradP.H"
Info<< "\nStarting time loop\n" << endl;
for(runTime++; !runTime.end(); runTime++)
{
Info<< "Time = " << runTime.timeName() << nl << endl;
#include "readPISOControls.H"
#include "CourantNo.H"
sgsModel->correct();
fvVectorMatrix UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
+ sgsModel->divDevBeff(U)
==
flowDirection*gradP
);
if (momentumPredictor)
{
solve(UEqn == -fvc::grad(p));
}
// --- PISO loop
volScalarField rUA = 1.0/UEqn.A();
for (int corr=0; corr<nCorr; corr++)
{
U = rUA*UEqn.H();
phi = (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, U, phi);
adjustPhi(phi, U, p);
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
(
fvm::laplacian(rUA, p) == fvc::div(phi)
);
pEqn.setReference(pRefCell, pRefValue);
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();
}
}
#include "continuityErrs.H"
U -= rUA*fvc::grad(p);
U.correctBoundaryConditions();
}
// Correct driving force for a constant mass flow rate
// Extract the velocity in the flow direction
dimensionedScalar magUbarStar =
(flowDirection & U)().weightedAverage(mesh.V());
// Calculate the pressure gradient increment needed to
// adjust the average flow-rate to the correct value
dimensionedScalar gragPplus =
(magUbar - magUbarStar)/rUA.weightedAverage(mesh.V());
U += flowDirection*rUA*gragPplus;
gradP += gragPplus;
Info<< "Uncorrected Ubar = " << magUbarStar.value() << tab
<< "pressure gradient = " << gradP.value() << endl;
runTime.write();
#include "writeGradP.H"
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return(0);
}
// ************************************************************************* //

43
applications/solvers/incompressible/channelFoam/createFields.H Normal file
View File

@ -0,0 +1,43 @@
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);
singlePhaseTransportModel laminarTransport(U, phi);
autoPtr<incompressible::LESModel> sgsModel
(
incompressible::LESModel::New(U, phi, laminarTransport)
);

24
applications/solvers/incompressible/channelFoam/createGradP.H Normal file
View File

@ -0,0 +1,24 @@
dimensionedScalar gradP
(
"gradP",
dimensionSet(0, 1, -2, 0, 0),
0.0
);
IFstream gradPFile
(
runTime.path()/runTime.timeName()/"uniform"/"gradP.raw"
);
if(gradPFile.good())
{
gradPFile >> gradP;
Info<< "Reading average pressure gradient" <<endl
<< endl;
}
else
{
Info<< "Initializing with 0 pressure gradient" <<endl
<< endl;
};

28
applications/solvers/incompressible/channelFoam/readTransportProperties.H Normal file
View File

@ -0,0 +1,28 @@
Info<< "\nReading transportProperties\n" << endl;
IOdictionary transportProperties
(
IOobject
(
"transportProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
)
);
dimensionedScalar nu
(
transportProperties.lookup("nu")
);
// Read centerline velocity for channel simulations
dimensionedVector Ubar
(
transportProperties.lookup("Ubar")
);
dimensionedScalar magUbar = mag(Ubar);
vector flowDirection = (Ubar/magUbar).value();

19
applications/solvers/incompressible/channelFoam/writeGradP.H Normal file
View File

@ -0,0 +1,19 @@
if (runTime.outputTime())
{
OFstream gradPFile
(
runTime.path()/runTime.timeName()/"uniform"/"gradP.raw"
);
if(gradPFile.good())
{
gradPFile << gradP << endl;
}
else
{
FatalErrorIn(args.executable())
<< "Cannot open file "
<< runTime.path()/runTime.timeName()/"uniform"/"gradP.raw"
<< exit(FatalError);
};
};

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

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

17
applications/solvers/incompressible/pimpleDyMFoam/Make/options Normal file
View File

@ -0,0 +1,17 @@
EXE_INC = \
-I$(LIB_SRC)/dynamicFvMesh/lnInclude \
-I$(LIB_SRC)/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 \
-ldynamicMesh \
-lmeshTools \
-lincompressibleTransportModels \
-lincompressibleRASModels \
-lincompressibleLESModels \
-lfiniteVolume

16
applications/solvers/incompressible/pimpleDyMFoam/UEqn.H Normal file
View File

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

44
applications/solvers/incompressible/pimpleDyMFoam/correctPhi.H Normal file
View File

@ -0,0 +1,44 @@
{
wordList pcorrTypes(p.boundaryField().types());
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.0),
pcorrTypes
);
for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pcorrEqn
(
fvm::laplacian(rAU, pcorr) == fvc::div(phi)
);
pcorrEqn.setReference(pRefCell, pRefValue);
pcorrEqn.solve();
if (nonOrth == nNonOrthCorr)
{
phi -= pcorrEqn.flux();
}
}
}
#include "continuityErrs.H"

59
applications/solvers/incompressible/pimpleDyMFoam/createFields.H Normal file
View File

@ -0,0 +1,59 @@
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);
singlePhaseTransportModel laminarTransport(U, phi);
autoPtr<incompressible::turbulenceModel> turbulence
(
incompressible::turbulenceModel::New(U, phi, laminarTransport)
);
Info<< "Reading field rAU if present\n" << endl;
volScalarField rAU
(
IOobject
(
"rAU",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mesh,
runTime.deltaT(),
zeroGradientFvPatchScalarField::typeName
);

172
applications/solvers/incompressible/pimpleDyMFoam/pimpleDyMFoam.C Normal file
View File

@ -0,0 +1,172 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 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
turbDyMFoam
Description
Transient solver for incompressible, flow of Newtonian fluids
on a moving 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 "readPISOControls.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"
// Make the fluxes absolute
fvc::makeAbsolute(phi, U);
# include "setDeltaT.H"
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
mesh.update();
if (mesh.changing() && correctPhi)
{
# include "correctPhi.H"
}
// Make the fluxes relative to the mesh motion
fvc::makeRelative(phi, U);
if (mesh.changing() && checkMeshCourantNo)
{
# include "meshCourantNo.H"
}
// --- PIMPLE loop
for (int ocorr=0; ocorr<nOuterCorr; ocorr++)
{
if (nOuterCorr != 1)
{
p.storePrevIter();
}
# include "UEqn.H"
// --- 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();
}
}
turbulence->correct();
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return(0);
}
// ************************************************************************* //

14
applications/solvers/incompressible/pimpleDyMFoam/readControls.H Normal file
View File

@ -0,0 +1,14 @@
# include "readTimeControls.H"
# include "readPISOControls.H"
bool correctPhi = false;
if (piso.found("correctPhi"))
{
correctPhi = Switch(piso.lookup("correctPhi"));
}
bool checkMeshCourantNo = false;
if (piso.found("checkMeshCourantNo"))
{
checkMeshCourantNo = Switch(piso.lookup("checkMeshCourantNo"));
}

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

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

12
applications/solvers/incompressible/pisoFoam/Make/options Normal file
View File

@ -0,0 +1,12 @@
EXE_INC = \
-I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
-I$(LIB_SRC)/finiteVolume/lnInclude
EXE_LIBS = \
-lincompressibleRASModels \
-lincompressibleLESModels \
-lincompressibleTransportModels \
-lfiniteVolume \
-lmeshTools

42
applications/solvers/incompressible/pisoFoam/createFields.H Normal file
View File

@ -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("PISO"), pRefCell, pRefValue);
singlePhaseTransportModel laminarTransport(U, phi);
autoPtr<incompressible::turbulenceModel> turbulence
(
incompressible::turbulenceModel::New(U, phi, laminarTransport)
);

131
applications/solvers/incompressible/pisoFoam/pisoFoam.C Normal file
View File

@ -0,0 +1,131 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 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
turbFoam
Description
Transient solver for incompressible flow.
Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "singlePhaseTransportModel.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"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
for (runTime++; !runTime.end(); runTime++)
{
Info<< "Time = " << runTime.timeName() << nl << endl;
# include "readPISOControls.H"
# include "CourantNo.H"
// Pressure-velocity PISO corrector
{
// Momentum predictor
fvVectorMatrix UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
+ turbulence->divDevReff(U)
);
if (momentumPredictor)
{
solve(UEqn == -fvc::grad(p));
}
// --- PISO loop
for (int corr=0; corr<nCorr; corr++)
{
volScalarField rUA = 1.0/UEqn.A();
U = rUA*UEqn.H();
phi = (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, U, phi);
adjustPhi(phi, U, p);
// Non-orthogonal pressure corrector loop
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
// Pressure corrector
fvScalarMatrix pEqn
(
fvm::laplacian(rUA, p) == fvc::div(phi)
);
pEqn.setReference(pRefCell, pRefValue);
pEqn.solve();
if (nonOrth == nNonOrthCorr)
{
phi -= pEqn.flux();
}
}
# include "continuityErrs.H"
U -= rUA*fvc::grad(p);
U.correctBoundaryConditions();
}
}
turbulence->correct();
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return(0);
}
// ************************************************************************* //

59
applications/solvers/multiphase/cavitatingFoam/CourantNo.H Normal file
View File

@ -0,0 +1,59 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 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
Global
CourantNo
Description
Calculates and outputs the mean and maximum Courant Numbers.
\*---------------------------------------------------------------------------*/
scalar CoNum = 0.0;
scalar meanCoNum = 0.0;
scalar acousticCoNum = 0.0;
if (mesh.nInternalFaces())
{
surfaceScalarField SfUfbyDelta =
mesh.surfaceInterpolation::deltaCoeffs()*mag(phiv);
CoNum = max(SfUfbyDelta/mesh.magSf())
.value()*runTime.deltaT().value();
meanCoNum = (sum(SfUfbyDelta)/sum(mesh.magSf()))
.value()*runTime.deltaT().value();
acousticCoNum = max
(
mesh.surfaceInterpolation::deltaCoeffs()/sqrt(fvc::interpolate(psi))
).value()*runTime.deltaT().value();
}
Info<< "phiv Courant Number mean: " << meanCoNum
<< " max: " << CoNum
<< " acoustic max: " << acousticCoNum
<< endl;
// ************************************************************************* //

3
applications/solvers/multiphase/cavitatingFoam/Make/files Normal file
View File

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

14
applications/solvers/multiphase/cavitatingFoam/Make/options Normal file
View File

@ -0,0 +1,14 @@
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/lnInclude \
-I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
-I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
-I$(LIB_SRC)/thermophysicalModels/barotropicCompressibilityModel/lnInclude
EXE_LIBS = \
-lincompressibleTransportModels \
-lincompressibleRASModels \
-lincompressibleLESModels \
-lfiniteVolume \
-lbarotropicCompressibilityModel

22
applications/solvers/multiphase/cavitatingFoam/UEqn.H Normal file
View File

@ -0,0 +1,22 @@
surfaceScalarField muEff
(
"muEff",
twoPhaseProperties.muf()
+ fvc::interpolate(rho*turbulence->nut())
);
fvVectorMatrix UEqn
(
fvm::ddt(rho, U)
+ fvm::div(phi, U)
- fvm::laplacian(muEff, U)
//- (fvc::grad(U) & fvc::grad(muf))
- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
);
UEqn.relax();
if (momentumPredictor)
{
solve(UEqn == -fvc::grad(p));
}

96
applications/solvers/multiphase/cavitatingFoam/cavitatingFoam.C Normal file
View File

@ -0,0 +1,96 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 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
cavitatingFoam
Description
Transient cavitation code based on the barotropic equation of state.
Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "barotropicCompressibilityModel.H"
#include "twoPhaseMixture.H"
#include "turbulenceModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
# include "setRootCase.H"
# include "createTime.H"
# include "createMesh.H"
# include "readThermodynamicProperties.H"
# include "readControls.H"
# include "createFields.H"
# include "initContinuityErrs.H"
# include "compressibleCourantNo.H"
# include "setInitialDeltaT.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.run())
{
# include "readControls.H"
# include "CourantNo.H"
# include "setDeltaT.H"
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
for (int outerCorr=0; outerCorr<nOuterCorr; outerCorr++)
{
# include "rhoEqn.H"
# include "gammaPsi.H"
# include "UEqn.H"
for (int corr=0; corr<nCorr; corr++)
{
# include "pEqn.H"
}
}
turbulence->correct();
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "\n end \n";
return(0);
}
// ************************************************************************* //

22
applications/solvers/multiphase/cavitatingFoam/continuityErrs.H Normal file
View File

@ -0,0 +1,22 @@
{
volScalarField thermoRho = psi*p + (1.0 - gamma)*rhol0;
dimensionedScalar totalMass = fvc::domainIntegrate(rho);
scalar sumLocalContErr =
(
fvc::domainIntegrate(mag(rho - thermoRho))/totalMass
).value();
scalar globalContErr =
(
fvc::domainIntegrate(rho - thermoRho)/totalMass
).value();
cumulativeContErr += globalContErr;
Info<< "time step continuity errors : sum local = " << sumLocalContErr
<< ", global = " << globalContErr
<< ", cumulative = " << cumulativeContErr
<< endl;
}

85
applications/solvers/multiphase/cavitatingFoam/createFields.H Normal file
View File

@ -0,0 +1,85 @@
Info<< "Reading field p\n" << endl;
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
volScalarField gamma
(
IOobject
(
"gamma",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
max(min((rho - rholSat)/(rhovSat - rholSat), scalar(1)), scalar(0))
);
gamma.oldTime();
Info<< "Creating compressibilityModel\n" << endl;
autoPtr<barotropicCompressibilityModel> psiModel =
barotropicCompressibilityModel::New
(
thermodynamicProperties,
gamma
);
const volScalarField& psi = psiModel->psi();
rho == max
(
psi*p
+ (1.0 - gamma)*rhol0
+ ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat,
rhoMin
);
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
# include "createPhiv.H"
# include "compressibleCreatePhi.H"
Info<< "Reading transportProperties\n" << endl;
twoPhaseMixture twoPhaseProperties(U, phiv, "gamma");
// Create incompressible turbulence model
autoPtr<incompressible::turbulenceModel> turbulence
(
incompressible::turbulenceModel::New(U, phiv, twoPhaseProperties)
);

10
applications/solvers/multiphase/cavitatingFoam/gammaPsi.H Normal file
View File

@ -0,0 +1,10 @@
{
gamma = max(min((rho - rholSat)/(rhovSat - rholSat), scalar(1)), scalar(0));
Info<< "max-min gamma: " << max(gamma).value()
<< " " << min(gamma).value() << endl;
psiModel->correct();
//Info<< "min a: " << 1.0/sqrt(max(psi)).value() << endl;
}

80
applications/solvers/multiphase/cavitatingFoam/pEqn.H Normal file
View File

@ -0,0 +1,80 @@
{
if (nOuterCorr == 1)
{
p =
(
rho
- (1.0 - gamma)*rhol0
- ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat
)/psi;
}
surfaceScalarField rhof = fvc::interpolate(rho, "rhof");
volScalarField rUA = 1.0/UEqn.A();
surfaceScalarField rUAf("rUAf", rhof*fvc::interpolate(rUA));
volVectorField HbyA = rUA*UEqn.H();
phiv = (fvc::interpolate(HbyA) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phiv);
p.boundaryField().updateCoeffs();
surfaceScalarField phiGradp = rUAf*mesh.magSf()*fvc::snGrad(p);
phiv -= phiGradp/rhof;
# include "resetPhivPatches.H"
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
(
fvm::ddt(psi, p)
- (rhol0 + (psil - psiv)*pSat)*fvc::ddt(gamma) - pSat*fvc::ddt(psi)
+ fvc::div(phiv, rho)
+ fvc::div(phiGradp)
- fvm::laplacian(rUAf, p)
);
pEqn.solve();
if (nonOrth == nNonOrthCorr)
{
phiv += (phiGradp + pEqn.flux())/rhof;
}
}
Info<< "max-min p: " << max(p).value()
<< " " << min(p).value() << endl;
U = HbyA - rUA*fvc::grad(p);
// Remove the swirl component of velocity for "wedge" cases
if (piso.found("removeSwirl"))
{
label swirlCmpt(readLabel(piso.lookup("removeSwirl")));
Info<< "Removing swirl component-" << swirlCmpt << " of U" << endl;
U.field().replace(swirlCmpt, 0.0);
}
U.correctBoundaryConditions();
Info<< "max(U) " << max(mag(U)).value() << endl;
rho == max
(
psi*p
+ (1.0 - gamma)*rhol0
+ ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat,
rhoMin
);
Info<< "max-min rho: " << max(rho).value()
<< " " << min(rho).value() << endl;
# include "gammaPsi.H"
}

9
applications/solvers/multiphase/cavitatingFoam/readControls.H Normal file
View File

@ -0,0 +1,9 @@
#include "readTimeControls.H"
scalar maxAcousticCo
(
readScalar(runTime.controlDict().lookup("maxAcousticCo"))
);
#include "readPISOControls.H"

27
applications/solvers/multiphase/cavitatingFoam/readThermodynamicProperties.H Normal file
View File

@ -0,0 +1,27 @@
Info<< "Reading thermodynamicProperties\n" << endl;
IOdictionary thermodynamicProperties
(
IOobject
(
"thermodynamicProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
)
);
dimensionedScalar psil(thermodynamicProperties.lookup("psil"));
dimensionedScalar rholSat(thermodynamicProperties.lookup("rholSat"));
dimensionedScalar psiv(thermodynamicProperties.lookup("psiv"));
dimensionedScalar pSat(thermodynamicProperties.lookup("pSat"));
dimensionedScalar rhovSat("rhovSat", psiv*pSat);
dimensionedScalar rhol0("rhol0", rholSat - pSat*psil);
dimensionedScalar rhoMin(thermodynamicProperties.lookup("rhoMin"));

15
applications/solvers/multiphase/cavitatingFoam/resetPhiPatches.H Normal file
View File

@ -0,0 +1,15 @@
fvsPatchScalarFieldField& phiPatches = phi.boundaryField();
const fvPatchScalarFieldField& rhoPatches = rho.boundaryField();
const fvPatchVectorFieldField& Upatches = U.boundaryField();
const fvsPatchVectorFieldField& SfPatches = mesh.Sf().boundaryField();
forAll(phiPatches, patchI)
{
if (phi.boundaryField().types()[patchI] == "calculated")
{
calculatedFvsPatchScalarField& phiPatch =
refCast<calculatedFvsPatchScalarField>(phiPatches[patchI]);
phiPatch == ((rhoPatches[patchI]*Upatches[patchI]) & SfPatches[patchI]);
}
}

14
applications/solvers/multiphase/cavitatingFoam/resetPhivPatches.H Normal file
View File

@ -0,0 +1,14 @@
surfaceScalarField::GeometricBoundaryField& phivPatches = phiv.boundaryField();
const volVectorField::GeometricBoundaryField& Upatches = U.boundaryField();
const surfaceVectorField::GeometricBoundaryField& SfPatches = mesh.Sf().boundaryField();
forAll(phivPatches, patchI)
{
if (phiv.boundaryField().types()[patchI] == "calculated")
{
calculatedFvsPatchScalarField& phivPatch =
refCast<calculatedFvsPatchScalarField>(phivPatches[patchI]);
phivPatch == (Upatches[patchI] & SfPatches[patchI]);
}
}

16
applications/solvers/multiphase/cavitatingFoam/rhoEqn.H Normal file
View File

@ -0,0 +1,16 @@
{
fvScalarMatrix rhoEqn
(
fvm::ddt(rho)
+ fvm::div(phiv, rho)
);
rhoEqn.solve();
phi = rhoEqn.flux();
Info<< "max-min rho: " << max(rho).value()
<< " " << min(rho).value() << endl;
rho == max(rho, rhoMin);
}

54
applications/solvers/multiphase/cavitatingFoam/setDeltaT.H Normal file
View File

@ -0,0 +1,54 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 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
Global
setDeltaT
Description
Reset the timestep to maintain a constant maximum courant Number.
Reduction of time-step is imediate but increase is damped to avoid
unstable oscillations.
\*---------------------------------------------------------------------------*/
if (adjustTimeStep)
{
scalar maxDeltaTFact =
min(maxCo/(CoNum + SMALL), maxAcousticCo/(acousticCoNum + SMALL));
scalar deltaTFact = min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
runTime.setDeltaT
(
min
(
deltaTFact*runTime.deltaT().value(),
maxDeltaT
)
);
Info<< "deltaT = " << runTime.deltaT().value() << endl;
}
// ************************************************************************* //

54
applications/solvers/multiphase/cavitatingFoam/setInitialDeltaT.H Normal file
View File

@ -0,0 +1,54 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 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
Global
setInitialDeltaT
Description
Set the initial timestep corresponding to the timestep adjustment
algorithm in setDeltaT
\*---------------------------------------------------------------------------*/
if (adjustTimeStep)
{
# include "CourantNo.H"
if (CoNum > SMALL)
{
scalar maxDeltaTFact =
min(maxCo/(CoNum + SMALL), maxAcousticCo/(acousticCoNum + SMALL));
runTime.setDeltaT
(
min
(
maxDeltaTFact*runTime.deltaT().value(),
maxDeltaT
)
);
}
}
// ************************************************************************* //

3
applications/solvers/multiphase/compressibleInterDyMFoam/Make/files Normal file
View File

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

22
applications/solvers/multiphase/compressibleInterDyMFoam/Make/options Normal file
View File

@ -0,0 +1,22 @@
INTERFOAM = $(FOAM_SOLVERS)/multiphase/interFoam
EXE_INC = \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/lnInclude \
-I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
-I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/dynamicMesh/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/dynamicFvMesh/lnInclude
EXE_LIBS = \
-linterfaceProperties \
-lincompressibleTransportModels \
-lincompressibleRASModels \
-lincompressibleLESModels \
-lfiniteVolume \
-ldynamicMesh \
-lmeshTools \
-ldynamicFvMesh

31
applications/solvers/multiphase/compressibleInterDyMFoam/UEqn.H Normal file
View File

@ -0,0 +1,31 @@
surfaceScalarField muf =
twoPhaseProperties.muf()
+ fvc::interpolate(rho*turbulence->nut());
fvVectorMatrix UEqn
(
fvm::ddt(rho, U)
+ fvm::div(rhoPhi, U)
- fvm::laplacian(muf, U)
- (fvc::grad(U) & fvc::grad(muf))
//- fvc::div(muf*(mesh.Sf() & fvc::interpolate(fvc::grad(U)().T())))
);
UEqn.relax();
if (momentumPredictor)
{
solve
(
UEqn
==
fvc::reconstruct
(
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
- ghf*fvc::snGrad(rho)
- fvc::snGrad(pd)
) * mesh.magSf()
)
);
}

76
applications/solvers/multiphase/compressibleInterDyMFoam/alphaEqns.H Normal file
View File

@ -0,0 +1,76 @@
{
word alphaScheme("div(phi,alpha)");
word alpharScheme("div(phirb,alpha)");
surfaceScalarField phir = phic*interface.nHatf();
for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)
{
volScalarField::DimensionedInternalField Sp
(
IOobject
(
"Sp",
runTime.timeName(),
mesh
),
mesh,
dimensionedScalar("Sp", dgdt.dimensions(), 0.0)
);
volScalarField::DimensionedInternalField Su
(
IOobject
(
"Su",
runTime.timeName(),
mesh
),
// Divergence term is handled explicitly to be
// consistent with the explicit transport solution
divU*min(alpha1, scalar(1))
);
forAll(dgdt, celli)
{
if (dgdt[celli] > 0.0 && alpha1[celli] > 0.0)
{
Sp[celli] -= dgdt[celli]*alpha1[celli];
Su[celli] += dgdt[celli]*alpha1[celli];
}
else if (dgdt[celli] < 0.0 && alpha1[celli] < 1.0)
{
Sp[celli] += dgdt[celli]*(1.0 - alpha1[celli]);
}
}
surfaceScalarField phiAlpha1 =
fvc::flux
(
phi,
alpha1,
alphaScheme
)
+ fvc::flux
(
-fvc::flux(-phir, alpha2, alpharScheme),
alpha1,
alpharScheme
);
MULES::explicitSolve(oneField(), alpha1, phi, phiAlpha1, Sp, Su, 1, 0);
surfaceScalarField rho1f = fvc::interpolate(rho1);
surfaceScalarField rho2f = fvc::interpolate(rho2);
rhoPhi = phiAlpha1*(rho1f - rho2f) + phi*rho2f;
alpha2 = scalar(1) - alpha1;
}
Info<< "Liquid phase volume fraction = "
<< alpha1.weightedAverage(mesh.V()).value()
<< " Min(alpha1) = " << min(alpha1).value()
<< " Min(alpha2) = " << min(alpha2).value()
<< endl;
}

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

@ -0,0 +1,43 @@
{
label nAlphaCorr
(
readLabel(piso.lookup("nAlphaCorr"))
);
label nAlphaSubCycles
(
readLabel(piso.lookup("nAlphaSubCycles"))
);
surfaceScalarField phic = mag(phi/mesh.magSf());
phic = min(interface.cAlpha()*phic, max(phic));
volScalarField divU = fvc::div(phi);
if (nAlphaSubCycles > 1)
{
dimensionedScalar totalDeltaT = runTime.deltaT();
surfaceScalarField rhoPhiSum = 0.0*rhoPhi;
for
(
subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
!(++alphaSubCycle).end();
)
{
# include "alphaEqns.H"
rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
}
rhoPhi = rhoPhiSum;
}
else
{
# include "alphaEqns.H"
}
if (oCorr == 0)
{
interface.correct();
}
}

138
applications/solvers/multiphase/compressibleInterDyMFoam/compressibleInterDyMFoam.C Normal file
View File

@ -0,0 +1,138 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 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
compressibleLesInterFoam
Description
Solver for 2 compressible, isothermal immiscible fluids using a VOF
(volume of fluid) phase-fraction based interface capturing approach.
The momentum and other fluid properties are of the "mixture" and a single
momentum equation is solved.
Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "dynamicFvMesh.H"
#include "MULES.H"
#include "subCycle.H"
#include "interfaceProperties.H"
#include "twoPhaseMixture.H"
#include "turbulenceModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createDynamicFvMesh.H"
#include "readEnvironmentalProperties.H"
#include "readControls.H"
#include "initContinuityErrs.H"
#include "createFields.H"
#include "CourantNo.H"
#include "setInitialDeltaT.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.run())
{
#include "readControls.H"
#include "CourantNo.H"
// Make the fluxes absolute
fvc::makeAbsolute(phi, U);
#include "setDeltaT.H"
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime();
// Do any mesh changes
mesh.update();
if (mesh.changing())
{
Info<< "Execution time for mesh.update() = "
<< runTime.elapsedCpuTime() - timeBeforeMeshUpdate
<< " s" << endl;
gh = g & mesh.C();
ghf = g & mesh.Cf();
}
if (mesh.changing() && correctPhi)
{
//***HGW#include "correctPhi.H"
}
// Make the fluxes relative to the mesh motion
fvc::makeRelative(phi, U);
if (mesh.changing() && checkMeshCourantNo)
{
#include "meshCourantNo.H"
}
turbulence->correct();
// --- Outer-corrector loop
for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
{
#include "alphaEqnsSubCycle.H"
solve(fvm::ddt(rho) + fvc::div(rhoPhi));
#include "UEqn.H"
// --- PISO loop
for (int corr=0; corr<nCorr; corr++)
{
#include "pEqn.H"
}
}
rho = alpha1*rho1 + alpha2*rho2;
runTime.write();
Info<< "ExecutionTime = "
<< runTime.elapsedCpuTime()
<< " s\n\n" << endl;
}
Info<< "End\n" << endl;
return(0);
}
// ************************************************************************* //

152
applications/solvers/multiphase/compressibleInterDyMFoam/createFields.H Normal file
View File

@ -0,0 +1,152 @@
Info<< "Reading field pd\n" << endl;
volScalarField pd
(
IOobject
(
"pd",
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<< "Calculating field alpha1\n" << endl;
volScalarField alpha2("alpha2", scalar(1) - alpha1);
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
#include "createPhi.H"
Info<< "Calculating field g.h\n" << endl;
volScalarField gh("gh", g & mesh.C());
surfaceScalarField ghf("ghf", g & mesh.Cf());
Info<< "Reading transportProperties\n" << endl;
twoPhaseMixture twoPhaseProperties(U, phi);
dimensionedScalar rho10
(
twoPhaseProperties.subDict
(
twoPhaseProperties.phase1Name()
).lookup("rho0")
);
dimensionedScalar rho20
(
twoPhaseProperties.subDict
(
twoPhaseProperties.phase2Name()
).lookup("rho0")
);
dimensionedScalar psi1
(
twoPhaseProperties.subDict
(
twoPhaseProperties.phase1Name()
).lookup("psi")
);
dimensionedScalar psi2
(
twoPhaseProperties.subDict
(
twoPhaseProperties.phase2Name()
).lookup("psi")
);
dimensionedScalar pMin(twoPhaseProperties.lookup("pMin"));
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
max
(
(pd + gh*(alpha1*rho10 + alpha2*rho20))
/(1.0 - gh*(alpha1*psi1 + alpha2*psi2)),
pMin
)
);
volScalarField rho1 = rho10 + psi1*p;
volScalarField rho2 = rho20 + psi2*p;
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
alpha1*rho1 + alpha2*rho2
);
// Mass flux
// Initialisation does not matter because rhoPhi is reset after the
// alpha1 solution before it is used in the U equation.
surfaceScalarField rhoPhi
(
IOobject
(
"rho*phi",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
fvc::interpolate(rho)*phi
);
volScalarField dgdt =
pos(alpha2)*fvc::div(phi)/max(alpha2, scalar(0.0001));
// Construct interface from alpha1 distribution
interfaceProperties interface(alpha1, U, twoPhaseProperties);
// Construct incompressible turbulence model
autoPtr<incompressible::turbulenceModel> turbulence
(
incompressible::turbulenceModel::New(U, phi, twoPhaseProperties)
);

77
applications/solvers/multiphase/compressibleInterDyMFoam/pEqn.H Normal file
View File

@ -0,0 +1,77 @@
{
volScalarField rUA = 1.0/UEqn.A();
surfaceScalarField rUAf = fvc::interpolate(rUA);
tmp<fvScalarMatrix> pdEqnComp;
if (transonic)
{
pdEqnComp =
(fvm::ddt(pd) + fvm::div(phi, pd) - fvm::Sp(fvc::div(phi), pd));
}
else
{
pdEqnComp =
(fvm::ddt(pd) + fvc::div(phi, pd) - fvc::Sp(fvc::div(phi), pd));
}
U = rUA*UEqn.H();
surfaceScalarField phiU
(
"phiU",
(fvc::interpolate(U) & mesh.Sf()) + fvc::ddtPhiCorr(rUA, rho, U, phi)
);
phi = phiU +
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
- ghf*fvc::snGrad(rho)
)*rUAf*mesh.magSf();
for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pdEqnIncomp
(
fvc::div(phi)
- fvm::laplacian(rUAf, pd)
);
solve
(
(
max(alpha1, scalar(0))*(psi1/rho1)
+ max(alpha2, scalar(0))*(psi2/rho2)
)
*pdEqnComp()
+ pdEqnIncomp
);
if (nonOrth == nNonOrthCorr)
{
dgdt =
(pos(alpha2)*(psi2/rho2) - pos(alpha1)*(psi1/rho1))
*(pdEqnComp & pd);
phi += pdEqnIncomp.flux();
}
}
U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
U.correctBoundaryConditions();
p = max
(
(pd + gh*(alpha1*rho10 + alpha2*rho20))/(1.0 - gh*(alpha1*psi1 + alpha2*psi2)),
pMin
);
rho1 = rho10 + psi1*p;
rho2 = rho20 + psi2*p;
Info<< "max(U) " << max(mag(U)).value() << endl;
Info<< "min(pd) " << min(pd).value() << endl;
// Make the fluxes relative to the mesh motion
fvc::makeRelative(phi, U);
}

32
applications/solvers/multiphase/compressibleInterDyMFoam/readControls.H Normal file
View File

@ -0,0 +1,32 @@
#include "readPISOControls.H"
#include "readTimeControls.H"
label nAlphaCorr
(
readLabel(piso.lookup("nAlphaCorr"))
);
label nAlphaSubCycles
(
readLabel(piso.lookup("nAlphaSubCycles"))
);
if (nAlphaSubCycles > 1 && nOuterCorr != 1)
{
FatalErrorIn(args.executable())
<< "Sub-cycling alpha is only allowed for PISO, "
"i.e. when the number of outer-correctors = 1"
<< exit(FatalError);
}
bool correctPhi = true;
if (piso.found("correctPhi"))
{
correctPhi = Switch(piso.lookup("correctPhi"));
}
bool checkMeshCourantNo = false;
if (piso.found("checkMeshCourantNo"))
{
checkMeshCourantNo = Switch(piso.lookup("checkMeshCourantNo"));
}

3
applications/solvers/multiphase/compressibleInterFoam/Make/files Normal file
View File

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

15
applications/solvers/multiphase/compressibleInterFoam/Make/options Normal file
View File

@ -0,0 +1,15 @@
INTERFOAM = $(FOAM_SOLVERS)/multiphase/interFoam
EXE_INC = \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/lnInclude \
-I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
-I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
-I$(LIB_SRC)/finiteVolume/lnInclude
EXE_LIBS = \
-linterfaceProperties \
-lincompressibleTransportModels \
-lincompressibleRASModels \
-lincompressibleLESModels \
-lfiniteVolume

34
applications/solvers/multiphase/compressibleInterFoam/UEqn.H Normal file
View File

@ -0,0 +1,34 @@
surfaceScalarField muEff
(
"muEff",
twoPhaseProperties.muf()
+ fvc::interpolate(rho*turbulence->nut())
);
fvVectorMatrix UEqn
(
fvm::ddt(rho, U)
+ fvm::div(rhoPhi, U)
- fvm::laplacian(muEff, U)
- (fvc::grad(U) & fvc::grad(muEff))
//- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
);
UEqn.relax();
if (momentumPredictor)
{
solve
(
UEqn
==
fvc::reconstruct
(
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
- ghf*fvc::snGrad(rho)
- fvc::snGrad(pd)
) * mesh.magSf()
)
);
}

76
applications/solvers/multiphase/compressibleInterFoam/alphaEqns.H Normal file
View File

@ -0,0 +1,76 @@
{
word alphaScheme("div(phi,alpha)");
word alpharScheme("div(phirb,alpha)");
surfaceScalarField phir = phic*interface.nHatf();
for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)
{
volScalarField::DimensionedInternalField Sp
(
IOobject
(
"Sp",
runTime.timeName(),
mesh
),
mesh,
dimensionedScalar("Sp", dgdt.dimensions(), 0.0)
);
volScalarField::DimensionedInternalField Su
(
IOobject
(
"Su",
runTime.timeName(),
mesh
),
// Divergence term is handled explicitly to be
// consistent with the explicit transport solution
divU*min(alpha1, scalar(1))
);
forAll(dgdt, celli)
{
if (dgdt[celli] > 0.0 && alpha1[celli] > 0.0)
{
Sp[celli] -= dgdt[celli]*alpha1[celli];
Su[celli] += dgdt[celli]*alpha1[celli];
}
else if (dgdt[celli] < 0.0 && alpha1[celli] < 1.0)
{
Sp[celli] += dgdt[celli]*(1.0 - alpha1[celli]);
}
}
surfaceScalarField phiAlpha1 =
fvc::flux
(
phi,
alpha1,
alphaScheme
)
+ fvc::flux
(
-fvc::flux(-phir, alpha2, alpharScheme),
alpha1,
alpharScheme
);
MULES::explicitSolve(oneField(), alpha1, phi, phiAlpha1, Sp, Su, 1, 0);
surfaceScalarField rho1f = fvc::interpolate(rho1);
surfaceScalarField rho2f = fvc::interpolate(rho2);
rhoPhi = phiAlpha1*(rho1f - rho2f) + phi*rho2f;
alpha2 = scalar(1) - alpha1;
}
Info<< "Liquid phase volume fraction = "
<< alpha1.weightedAverage(mesh.V()).value()
<< " Min(alpha1) = " << min(alpha1).value()
<< " Min(alpha2) = " << min(alpha2).value()
<< endl;
}

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

@ -0,0 +1,43 @@
{
label nAlphaCorr
(
readLabel(piso.lookup("nAlphaCorr"))
);
label nAlphaSubCycles
(
readLabel(piso.lookup("nAlphaSubCycles"))
);
surfaceScalarField phic = mag(phi/mesh.magSf());
phic = min(interface.cAlpha()*phic, max(phic));
volScalarField divU = fvc::div(phi);
if (nAlphaSubCycles > 1)
{
dimensionedScalar totalDeltaT = runTime.deltaT();
surfaceScalarField rhoPhiSum = 0.0*rhoPhi;
for
(
subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
!(++alphaSubCycle).end();
)
{
#include "alphaEqns.H"
rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
}
rhoPhi = rhoPhiSum;
}
else
{
#include "alphaEqns.H"
}
if (oCorr == 0)
{
interface.correct();
}
}

106
applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C Normal file
View File

@ -0,0 +1,106 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 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
compressibleLesInterFoam
Description
Solver for 2 compressible, isothermal immiscible fluids using a VOF
(volume of fluid) phase-fraction based interface capturing approach.
The momentum and other fluid properties are of the "mixture" and a single
momentum equation is solved.
Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "MULES.H"
#include "subCycle.H"
#include "interfaceProperties.H"
#include "twoPhaseMixture.H"
#include "turbulenceModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "readEnvironmentalProperties.H"
#include "readControls.H"
#include "initContinuityErrs.H"
#include "createFields.H"
#include "CourantNo.H"
#include "setInitialDeltaT.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.run())
{
#include "readControls.H"
#include "CourantNo.H"
#include "setDeltaT.H"
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
// --- Outer-corrector loop
for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
{
#include "alphaEqnsSubCycle.H"
solve(fvm::ddt(rho) + fvc::div(rhoPhi));
#include "UEqn.H"
// --- PISO loop
for (int corr=0; corr<nCorr; corr++)
{
#include "pEqn.H"
}
}
rho = alpha1*rho1 + alpha2*rho2;
turbulence->correct();
runTime.write();
Info<< "ExecutionTime = "
<< runTime.elapsedCpuTime()
<< " s\n\n" << endl;
}
Info<< "End\n" << endl;
return(0);
}
// ************************************************************************* //

152
applications/solvers/multiphase/compressibleInterFoam/createFields.H Normal file
View File

@ -0,0 +1,152 @@
Info<< "Reading field pd\n" << endl;
volScalarField pd
(
IOobject
(
"pd",
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<< "Calculating field alpha1\n" << endl;
volScalarField alpha2("alpha2", scalar(1) - alpha1);
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
#include "createPhi.H"
Info<< "Calculating field g.h\n" << endl;
volScalarField gh("gh", g & mesh.C());
surfaceScalarField ghf("ghf", g & mesh.Cf());
Info<< "Reading transportProperties\n" << endl;
twoPhaseMixture twoPhaseProperties(U, phi);
dimensionedScalar rho10
(
twoPhaseProperties.subDict
(
twoPhaseProperties.phase1Name()
).lookup("rho0")
);
dimensionedScalar rho20
(
twoPhaseProperties.subDict
(
twoPhaseProperties.phase2Name()
).lookup("rho0")
);
dimensionedScalar psi1
(
twoPhaseProperties.subDict
(
twoPhaseProperties.phase1Name()
).lookup("psi")
);
dimensionedScalar psi2
(
twoPhaseProperties.subDict
(
twoPhaseProperties.phase2Name()
).lookup("psi")
);
dimensionedScalar pMin(twoPhaseProperties.lookup("pMin"));
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
max
(
(pd + gh*(alpha1*rho10 + alpha2*rho20))
/(1.0 - gh*(alpha1*psi1 + alpha2*psi2)),
pMin
)
);
volScalarField rho1 = rho10 + psi1*p;
volScalarField rho2 = rho20 + psi2*p;
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
alpha1*rho1 + alpha2*rho2
);
// Mass flux
// Initialisation does not matter because rhoPhi is reset after the
// alpha1 solution before it is used in the U equation.
surfaceScalarField rhoPhi
(
IOobject
(
"rho*phi",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
fvc::interpolate(rho)*phi
);
volScalarField dgdt =
pos(alpha2)*fvc::div(phi)/max(alpha2, scalar(0.0001));
// Construct interface from alpha1 distribution
interfaceProperties interface(alpha1, U, twoPhaseProperties);
// Construct incompressible turbulence model
autoPtr<incompressible::turbulenceModel> turbulence
(
incompressible::turbulenceModel::New(U, phi, twoPhaseProperties)
);

74
applications/solvers/multiphase/compressibleInterFoam/pEqn.H Normal file
View File

@ -0,0 +1,74 @@
{
volScalarField rUA = 1.0/UEqn.A();
surfaceScalarField rUAf = fvc::interpolate(rUA);
tmp<fvScalarMatrix> pdEqnComp;
if (transonic)
{
pdEqnComp =
(fvm::ddt(pd) + fvm::div(phi, pd) - fvm::Sp(fvc::div(phi), pd));
}
else
{
pdEqnComp =
(fvm::ddt(pd) + fvc::div(phi, pd) - fvc::Sp(fvc::div(phi), pd));
}
U = rUA*UEqn.H();
surfaceScalarField phiU
(
"phiU",
(fvc::interpolate(U) & mesh.Sf()) + fvc::ddtPhiCorr(rUA, rho, U, phi)
);
phi = phiU +
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
- ghf*fvc::snGrad(rho)
)*rUAf*mesh.magSf();
for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pdEqnIncomp
(
fvc::div(phi)
- fvm::laplacian(rUAf, pd)
);
solve
(
(
max(alpha1, scalar(0))*(psi1/rho1)
+ max(alpha2, scalar(0))*(psi2/rho2)
)
*pdEqnComp()
+ pdEqnIncomp
);
if (nonOrth == nNonOrthCorr)
{
dgdt =
(pos(alpha2)*(psi2/rho2) - pos(alpha1)*(psi1/rho1))
*(pdEqnComp & pd);
phi += pdEqnIncomp.flux();
}
}
U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
U.correctBoundaryConditions();
p = max
(
(pd + gh*(alpha1*rho10 + alpha2*rho20))/(1.0 - gh*(alpha1*psi1 + alpha2*psi2)),
pMin
);
rho1 = rho10 + psi1*p;
rho2 = rho20 + psi2*p;
Info<< "max(U) " << max(mag(U)).value() << endl;
Info<< "min(pd) " << min(pd).value() << endl;
}

20
applications/solvers/multiphase/compressibleInterFoam/readControls.H Normal file
View File

@ -0,0 +1,20 @@
#include "readPISOControls.H"
#include "readTimeControls.H"
label nAlphaCorr
(
readLabel(piso.lookup("nAlphaCorr"))
);
label nAlphaSubCycles
(
readLabel(piso.lookup("nAlphaSubCycles"))
);
if (nAlphaSubCycles > 1 && nOuterCorr != 1)
{
FatalErrorIn(args.executable())
<< "Sub-cycling alpha is only allowed for PISO, "
"i.e. when the number of outer-correctors = 1"
<< exit(FatalError);
}

35
applications/solvers/multiphase/interFoam/alphaEqn.H Normal file
View File

@ -0,0 +1,35 @@
{
word alphaScheme("div(phi,alpha)");
word alpharScheme("div(phirb,alpha)");
surfaceScalarField phic = mag(phi/mesh.magSf());
phic = min(interface.cAlpha()*phic, max(phic));
surfaceScalarField phir = phic*interface.nHatf();
for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)
{
surfaceScalarField phiAlpha =
fvc::flux
(
phi,
alpha1,
alphaScheme
)
+ fvc::flux
(
-fvc::flux(-phir, scalar(1) - alpha1, alpharScheme),
alpha1,
alpharScheme
);
MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0);
rhoPhi = phiAlpha*(rho1 - rho2) + phi*rho2;
}
Info<< "Liquid phase volume fraction = "
<< alpha1.weightedAverage(mesh.V()).value()
<< " Min(alpha1) = " << min(alpha1).value()
<< " Max(alpha1) = " << max(alpha1).value()
<< endl;
}

35
applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H Normal file
View File

@ -0,0 +1,35 @@
label nAlphaCorr
(
readLabel(piso.lookup("nAlphaCorr"))
);
label nAlphaSubCycles
(
readLabel(piso.lookup("nAlphaSubCycles"))
);
if (nAlphaSubCycles > 1)
{
dimensionedScalar totalDeltaT = runTime.deltaT();
surfaceScalarField rhoPhiSum = 0.0*rhoPhi;
for
(
subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
!(++alphaSubCycle).end();
)
{
# include "alphaEqn.H"
rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
}
rhoPhi = rhoPhiSum;
}
else
{
# include "alphaEqn.H"
}
interface.correct();
rho == alpha1*rho1 + (scalar(1) - alpha1)*rho2;

67
applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H Normal file
View File

@ -0,0 +1,67 @@
{
word alphaScheme("div(phi,alpha)");
word alpharScheme("div(phirb,alpha)");
surfaceScalarField phir("phir", phic*interface.nHatf());
for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)
{
surfaceScalarField phiAlpha =
fvc::flux
(
phi,
alpha1,
alphaScheme
)
+ fvc::flux
(
-fvc::flux(-phir, scalar(1) - alpha1, alpharScheme),
alpha1,
alpharScheme
);
Pair<tmp<volScalarField> > vDotAlphal =
twoPhaseProperties->vDotAlphal();
const volScalarField& vDotcAlphal = vDotAlphal[0]();
const volScalarField& vDotvAlphal = vDotAlphal[1]();
volScalarField Sp
(
IOobject
(
"Sp",
runTime.timeName(),
mesh
),
vDotvAlphal - vDotcAlphal
);
volScalarField Su
(
IOobject
(
"Su",
runTime.timeName(),
mesh
),
// Divergence term is handled explicitly to be
// consistent with the explicit transport solution
divU*alpha1
+ vDotcAlphal
);
//MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0);
//MULES::explicitSolve(oneField(), alpha1, phi, phiAlpha, Sp, Su, 1, 0);
MULES::implicitSolve(oneField(), alpha1, phi, phiAlpha, Sp, Su, 1, 0);
rhoPhi +=
(runTime.deltaT()/totalDeltaT)
*(phiAlpha*(rho1 - rho2) + phi*rho2);
}
Info<< "Liquid phase volume fraction = "
<< alpha1.weightedAverage(mesh.V()).value()
<< " Min(alpha1) = " << min(alpha1).value()
<< " Max(alpha1) = " << max(alpha1).value()
<< endl;
}

53
applications/solvers/multiphase/interPhaseChangeFoam/alphaEqnSubCycle.H Normal file
View File

@ -0,0 +1,53 @@
surfaceScalarField rhoPhi
(
IOobject
(
"rhoPhi",
runTime.timeName(),
mesh
),
mesh,
dimensionedScalar("0", dimensionSet(1, 0, -1, 0, 0), 0)
);
{
label nAlphaCorr
(
readLabel(piso.lookup("nAlphaCorr"))
);
label nAlphaSubCycles
(
readLabel(piso.lookup("nAlphaSubCycles"))
);
surfaceScalarField phic = mag(phi/mesh.magSf());
phic = min(interface.cAlpha()*phic, max(phic));
volScalarField divU = fvc::div(phi);
dimensionedScalar totalDeltaT = runTime.deltaT();
if (nAlphaSubCycles > 1)
{
for
(
subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
!(++alphaSubCycle).end();
)
{
# include "alphaEqn.H"
}
}
else
{
# include "alphaEqn.H"
}
if (nOuterCorr == 1)
{
interface.correct();
}
rho == alpha1*rho1 + (scalar(1) - alpha1)*rho2;
}

19
applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqn.H Normal file
View File

@ -0,0 +1,19 @@
{
fvScalarMatrix alpha1Eqn
(
fvm::ddt(alpha1)
+ fvm::div(phi, alpha1)
- fvm::laplacian(Dab, alpha1)
);
alpha1Eqn.solve();
rhoPhi = alpha1Eqn.flux()*(rho1 - rho2) + phi*rho2;
rho = alpha1*rho1 + (scalar(1) - alpha1)*rho2;
Info<< "Phase 1 volume fraction = "
<< alpha1.weightedAverage(mesh.V()).value()
<< " Min(alpha1) = " << min(alpha1).value()
<< " Max(alpha1) = " << max(alpha1).value()
<< endl;
}