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

MRG: Integrated Foundation code to commit 7d6845d

Browse Source
This commit is contained in:
Andrew Heather
2017-03-23 14:33:33 +00:00
parent 436ec1cf1f 7d6845defa
commit dcb1a95e35
595 changed files with 7175 additions and 11727 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
applications/solvers/multiphase/interFoam/Make/options
View File

@ -1,4 +1,5 @@
EXE_INC = \
-I../VoF \
-I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/lnInclude \

57
applications/solvers/multiphase/interFoam/alphaCourantNo.H
View File

@ -1,57 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ 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 3 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, see <http://www.gnu.org/licenses/>.
Global
alphaCourantNo
Description
Calculates and outputs the mean and maximum Courant Numbers.
\*---------------------------------------------------------------------------*/
scalar maxAlphaCo
(
readScalar(runTime.controlDict().lookup("maxAlphaCo"))
);
scalar alphaCoNum = 0.0;
scalar meanAlphaCoNum = 0.0;
if (mesh.nInternalFaces())
{
scalarField sumPhi
(
mixture.nearInterface()().primitiveField()
*fvc::surfaceSum(mag(phi))().primitiveField()
);
alphaCoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
meanAlphaCoNum =
0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();
}
Info<< "Interface Courant Number mean: " << meanAlphaCoNum
<< " max: " << alphaCoNum << endl;
// ************************************************************************* //

217
applications/solvers/multiphase/interFoam/alphaEqn.H
View File

@ -1,217 +0,0 @@
{
word alphaScheme("div(phi,alpha)");
word alpharScheme("div(phirb,alpha)");
tmp<fv::ddtScheme<scalar>> ddtAlpha
(
fv::ddtScheme<scalar>::New
(
mesh,
mesh.ddtScheme("ddt(alpha)")
)
);
// Set the off-centering coefficient according to ddt scheme
scalar ocCoeff = 0;
if
(
isType<fv::EulerDdtScheme<scalar>>(ddtAlpha())
|| isType<fv::localEulerDdtScheme<scalar>>(ddtAlpha())
)
{
ocCoeff = 0;
}
else if (isType<fv::CrankNicolsonDdtScheme<scalar>>(ddtAlpha()))
{
if (nAlphaSubCycles > 1)
{
FatalErrorInFunction
<< "Sub-cycling is not supported "
"with the CrankNicolson ddt scheme"
<< exit(FatalError);
}
ocCoeff =
refCast<const fv::CrankNicolsonDdtScheme<scalar>>(ddtAlpha())
.ocCoeff();
}
else
{
FatalErrorInFunction
<< "Only Euler and CrankNicolson ddt schemes are supported"
<< exit(FatalError);
}
scalar cnCoeff = 1.0/(1.0 + ocCoeff);
// Standard face-flux compression coefficient
surfaceScalarField phic(mixture.cAlpha()*mag(phi/mesh.magSf()));
// Add the optional isotropic compression contribution
if (icAlpha > 0)
{
phic *= (1.0 - icAlpha);
phic += (mixture.cAlpha()*icAlpha)*fvc::interpolate(mag(U));
}
surfaceScalarField::Boundary& phicBf =
phic.boundaryFieldRef();
// Do not compress interface at non-coupled boundary faces
// (inlets, outlets etc.)
forAll(phic.boundaryField(), patchi)
{
fvsPatchScalarField& phicp = phicBf[patchi];
if (!phicp.coupled())
{
phicp == 0;
}
}
tmp<surfaceScalarField> phiCN(phi);
// Calculate the Crank-Nicolson off-centred volumetric flux
if (ocCoeff > 0)
{
phiCN = cnCoeff*phi + (1.0 - cnCoeff)*phi.oldTime();
}
if (MULESCorr)
{
fvScalarMatrix alpha1Eqn
(
(
LTS
? fv::localEulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
: fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
)
+ fv::gaussConvectionScheme<scalar>
(
mesh,
phiCN,
upwind<scalar>(mesh, phiCN)
).fvmDiv(phiCN, alpha1)
);
alpha1Eqn.solve();
Info<< "Phase-1 volume fraction = "
<< alpha1.weightedAverage(mesh.Vsc()).value()
<< " Min(" << alpha1.name() << ") = " << min(alpha1).value()
<< " Max(" << alpha1.name() << ") = " << max(alpha1).value()
<< endl;
tmp<surfaceScalarField> talphaPhiUD(alpha1Eqn.flux());
alphaPhi = talphaPhiUD();
if (alphaApplyPrevCorr && talphaPhiCorr0.valid())
{
Info<< "Applying the previous iteration compression flux" << endl;
MULES::correct(alpha1, alphaPhi, talphaPhiCorr0.ref(), 1, 0);
alphaPhi += talphaPhiCorr0();
}
// Cache the upwind-flux
talphaPhiCorr0 = talphaPhiUD;
alpha2 = 1.0 - alpha1;
mixture.correct();
}
for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
{
surfaceScalarField phir(phic*mixture.nHatf());
alphaPhiUn =
(
fvc::flux
(
phi,
alpha1,
alphaScheme
)
+ fvc::flux
(
-fvc::flux(-phir, alpha2, alpharScheme),
alpha1,
alpharScheme
)
);
// Calculate the Crank-Nicolson off-centred alpha flux
if (ocCoeff > 0)
{
alphaPhiUn =
cnCoeff*alphaPhiUn + (1.0 - cnCoeff)*alphaPhi.oldTime();
}
if (MULESCorr)
{
tmp<surfaceScalarField> talphaPhiCorr(alphaPhiUn - alphaPhi);
volScalarField alpha10("alpha10", alpha1);
MULES::correct(alpha1, alphaPhiUn, talphaPhiCorr.ref(), 1, 0);
// Under-relax the correction for all but the 1st corrector
if (aCorr == 0)
{
alphaPhi += talphaPhiCorr();
}
else
{
alpha1 = 0.5*alpha1 + 0.5*alpha10;
alphaPhi += 0.5*talphaPhiCorr();
}
}
else
{
alphaPhi = alphaPhiUn;
MULES::explicitSolve(alpha1, phiCN, alphaPhi, 1, 0);
}
alpha2 = 1.0 - alpha1;
mixture.correct();
}
if (alphaApplyPrevCorr && MULESCorr)
{
talphaPhiCorr0 = alphaPhi - talphaPhiCorr0;
talphaPhiCorr0.ref().rename("alphaPhiCorr0");
}
else
{
talphaPhiCorr0.clear();
}
if
(
word(mesh.ddtScheme("ddt(rho,U)"))
== fv::EulerDdtScheme<vector>::typeName
)
{
rhoPhi = alphaPhi*(rho1 - rho2) + phiCN*rho2;
}
else
{
if (ocCoeff > 0)
{
// Calculate the end-of-time-step alpha flux
alphaPhi = (alphaPhi - (1.0 - cnCoeff)*alphaPhi.oldTime())/cnCoeff;
}
// Calculate the end-of-time-step mass flux
rhoPhi = alphaPhi*(rho1 - rho2) + phi*rho2;
}
Info<< "Phase-1 volume fraction = "
<< alpha1.weightedAverage(mesh.Vsc()).value()
<< " Min(" << alpha1.name() << ") = " << min(alpha1).value()
<< " Max(" << alpha1.name() << ") = " << max(alpha1).value()
<< endl;
}

41
applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H
View File

@ -1,41 +0,0 @@
if (nAlphaSubCycles > 1)
{
dimensionedScalar totalDeltaT = runTime.deltaT();
surfaceScalarField rhoPhiSum
(
IOobject
(
"rhoPhiSum",
runTime.timeName(),
mesh
),
mesh,
dimensionedScalar("0", rhoPhi.dimensions(), 0)
);
tmp<volScalarField> trSubDeltaT;
if (LTS)
{
trSubDeltaT =
fv::localEulerDdt::localRSubDeltaT(mesh, nAlphaSubCycles);
}
for
(
subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
!(++alphaSubCycle).end();
)
{
#include "alphaEqn.H"
rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
}
rhoPhi = rhoPhiSum;
}
else
{
#include "alphaEqn.H"
}
rho == alpha1*rho1 + alpha2*rho2;

3
applications/solvers/multiphase/interFoam/alphaSuSp.H Normal file
View File

@ -0,0 +1,3 @@
zeroField Su;
zeroField Sp;
zeroField divU;

1
applications/solvers/multiphase/interFoam/interDyMFoam/Make/options
View File

@ -1,6 +1,7 @@
EXE_INC = \
-I. \
-I.. \
-I../../VoF \
-I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/lnInclude \

7
applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
View File

@ -47,7 +47,6 @@ Description
#include "pimpleControl.H"
#include "fvOptions.H"
#include "CorrectPhi.H"
#include "localEulerDdtScheme.H"
#include "fvcSmooth.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -129,7 +128,11 @@ int main(int argc, char *argv[])
<< " s" << endl;
// Do not apply previous time-step mesh compression flux
talphaPhiCorr0.clear();
// if the mesh topology changed
if (mesh.topoChanging())
{
talphaPhiCorr0.clear();
}
gh = (g & mesh.C()) - ghRef;
ghf = (g & mesh.Cf()) - ghRef;

2
applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H
View File

@ -28,7 +28,7 @@
phiHbyA += phig;
// Update the pressure BCs to ensure flux consistency
constrainPressure(p_rgh, U, phiHbyA, rAUf);
constrainPressure(p_rgh, U, phiHbyA, rAUf, MRF);
while (pimple.correctNonOrthogonal())
{

3
applications/solvers/multiphase/interFoam/interFoam.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -51,7 +51,6 @@ Description
#include "pimpleControl.H"
#include "fvOptions.H"
#include "CorrectPhi.H"
#include "localEulerDdtScheme.H"
#include "fvcSmooth.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

1
applications/solvers/multiphase/interFoam/interMixingFoam/Make/options
View File

@ -1,6 +1,7 @@
EXE_INC = \
-I. \
-I.. \
-I../../VoF \
-I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
-IimmiscibleIncompressibleThreePhaseMixture \
-IincompressibleThreePhaseMixture \

0
applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqns.H → applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqn.H
View File

4
applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqnsSubCycle.H → applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqnSubCycle.H
View File

@ -19,7 +19,7 @@ if (nAlphaSubCycles > 1)
!(++alphaSubCycle).end();
)
{
#include "alphaEqns.H"
#include "alphaEqn.H"
rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
}
@ -27,7 +27,7 @@ if (nAlphaSubCycles > 1)
}
else
{
#include "alphaEqns.H"
#include "alphaEqn.H"
}
{

4
applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -96,7 +96,7 @@ int main(int argc, char *argv[])
while (pimple.loop())
{
#include "alphaControls.H"
#include "alphaEqnsSubCycle.H"
#include "alphaEqnSubCycle.H"
mixture.correct();

2
applications/solvers/multiphase/interFoam/rhofs.H Normal file
View File

@ -0,0 +1,2 @@
const dimensionedScalar& rho1f(rho1);
const dimensionedScalar& rho2f(rho2);

53
applications/solvers/multiphase/interFoam/setDeltaT.H
View File

@ -1,53 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ 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 3 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, see <http://www.gnu.org/licenses/>.
Global
setDeltaT
Description
Reset the timestep to maintain a constant maximum courant Number.
Reduction of time-step is immediate, but increase is damped to avoid
unstable oscillations.
\*---------------------------------------------------------------------------*/
if (adjustTimeStep)
{
scalar maxDeltaTFact =
min(maxCo/(CoNum + SMALL), maxAlphaCo/(alphaCoNum + SMALL));
scalar deltaTFact = min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
runTime.setDeltaT
(
min
(
deltaTFact*runTime.deltaTValue(),
maxDeltaT
)
);
Info<< "deltaT = " << runTime.deltaTValue() << endl;
}
// ************************************************************************* //

136
applications/solvers/multiphase/interFoam/setRDeltaT.H
View File

@ -1,136 +0,0 @@
{
volScalarField& rDeltaT = trDeltaT.ref();
const dictionary& pimpleDict = pimple.dict();
scalar maxCo
(
pimpleDict.lookupOrDefault<scalar>("maxCo", 0.9)
);
scalar maxAlphaCo
(
pimpleDict.lookupOrDefault<scalar>("maxAlphaCo", 0.2)
);
scalar rDeltaTSmoothingCoeff
(
pimpleDict.lookupOrDefault<scalar>("rDeltaTSmoothingCoeff", 0.1)
);
label nAlphaSpreadIter
(
pimpleDict.lookupOrDefault<label>("nAlphaSpreadIter", 1)
);
scalar alphaSpreadDiff
(
pimpleDict.lookupOrDefault<scalar>("alphaSpreadDiff", 0.2)
);
scalar alphaSpreadMax
(
pimpleDict.lookupOrDefault<scalar>("alphaSpreadMax", 0.99)
);
scalar alphaSpreadMin
(
pimpleDict.lookupOrDefault<scalar>("alphaSpreadMin", 0.01)
);
label nAlphaSweepIter
(
pimpleDict.lookupOrDefault<label>("nAlphaSweepIter", 5)
);
scalar rDeltaTDampingCoeff
(
pimpleDict.lookupOrDefault<scalar>("rDeltaTDampingCoeff", 1.0)
);
scalar maxDeltaT
(
pimpleDict.lookupOrDefault<scalar>("maxDeltaT", GREAT)
);
volScalarField rDeltaT0("rDeltaT0", rDeltaT);
// Set the reciprocal time-step from the local Courant number
rDeltaT.ref() = max
(
1/dimensionedScalar("maxDeltaT", dimTime, maxDeltaT),
fvc::surfaceSum(mag(rhoPhi))()()
/((2*maxCo)*mesh.V()*rho())
);
if (maxAlphaCo < maxCo)
{
// Further limit the reciprocal time-step
// in the vicinity of the interface
volScalarField alpha1Bar(fvc::average(alpha1));
rDeltaT.ref() = max
(
rDeltaT(),
pos(alpha1Bar() - alphaSpreadMin)
*pos(alphaSpreadMax - alpha1Bar())
*fvc::surfaceSum(mag(phi))()()
/((2*maxAlphaCo)*mesh.V())
);
}
// Update tho boundary values of the reciprocal time-step
rDeltaT.correctBoundaryConditions();
Info<< "Flow time scale min/max = "
<< gMin(1/rDeltaT.primitiveField())
<< ", " << gMax(1/rDeltaT.primitiveField()) << endl;
if (rDeltaTSmoothingCoeff < 1.0)
{
fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff);
}
if (nAlphaSpreadIter > 0)
{
fvc::spread
(
rDeltaT,
alpha1,
nAlphaSpreadIter,
alphaSpreadDiff,
alphaSpreadMax,
alphaSpreadMin
);
}
if (nAlphaSweepIter > 0)
{
fvc::sweep(rDeltaT, alpha1, nAlphaSweepIter, alphaSpreadDiff);
}
Info<< "Smoothed flow time scale min/max = "
<< gMin(1/rDeltaT.primitiveField())
<< ", " << gMax(1/rDeltaT.primitiveField()) << endl;
// Limit rate of change of time scale
// - reduce as much as required
// - only increase at a fraction of old time scale
if
(
rDeltaTDampingCoeff < 1.0
&& runTime.timeIndex() > runTime.startTimeIndex() + 1
)
{
rDeltaT = max
(
rDeltaT,
(scalar(1.0) - rDeltaTDampingCoeff)*rDeltaT0
);
Info<< "Damped flow time scale min/max = "
<< gMin(1/rDeltaT.primitiveField())
<< ", " << gMax(1/rDeltaT.primitiveField()) << endl;
}
}