Commit 722363b4 authored by Andrew Heather's avatar Andrew Heather

INT: Updated IH Cantabria code contributions

- updated to streamFunction wave code and tutorial case
- renamed solver to interPorousFoam
- solver rebased on current interFoam solver
- relocated code/cases to appropriate directories
parent 750a7eba
interPorousFoam.C
fvOptions/multiphasePorosity/multiphasePorositySource.C
EXE = $(FOAM_APPBIN)/interPorousFoam
EXE_INC = \
-I. \
-I.. \
-I../../VoF \
-I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/lnInclude \
......
......@@ -2,46 +2,57 @@
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)
tmp<fv::ddtScheme<scalar>> tddtAlpha
(
fv::ddtScheme<scalar>::New
(
mesh,
mesh.ddtScheme("ddt(alpha)")
)
);
const fv::ddtScheme<scalar>& ddtAlpha = tddtAlpha();
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);
}
if
(
alphaRestart
|| mesh.time().timeIndex() > mesh.time().startTimeIndex() + 1
)
{
ocCoeff =
refCast<const fv::CrankNicolsonDdtScheme<scalar>>(ddtAlpha)
.ocCoeff();
}
}
else
{
FatalErrorInFunction
<< "Sub-cycling is not supported "
"with the CrankNicolson ddt scheme"
<< "Only Euler and CrankNicolson ddt schemes are supported"
<< exit(FatalError);
}
ocCoeff =
refCast<const fv::CrankNicolsonDdtScheme<scalar>>(ddtAlpha())
.ocCoeff();
}
else
{
FatalErrorInFunction
<< "Only Euler and CrankNicolson ddt schemes are supported"
<< exit(FatalError);
}
// Set the time blending factor, 1 for Euler
scalar cnCoeff = 1.0/(1.0 + ocCoeff);
// Standard face-flux compression coefficient
......@@ -54,6 +65,14 @@
phic += (mixture.cAlpha()*icAlpha)*fvc::interpolate(mag(U));
}
// Add the optional shear compression contribution
if (scAlpha > 0)
{
phic +=
scAlpha*mag(mesh.delta() & fvc::interpolate(symm(fvc::grad(U))));
}
surfaceScalarField::Boundary& phicBf =
phic.boundaryFieldRef();
......@@ -79,6 +98,8 @@
if (MULESCorr)
{
#include "alphaSuSp.H"
fvScalarMatrix alpha1Eqn
(
(
......@@ -92,6 +113,10 @@
phiCN,
upwind<scalar>(mesh, phiCN)
).fvmDiv(phiCN, alpha1)
// - fvm::Sp(fvc::ddt(dimensionedScalar("1", dimless, 1), mesh)
// + fvc::div(phiCN), alpha1)
==
Su + fvm::Sp(Sp + divU, alpha1)
);
alpha1Eqn.solve();
......@@ -102,28 +127,27 @@
<< " Max(" << alpha1.name() << ") = " << max(alpha1).value()
<< endl;
tmp<surfaceScalarField> talphaPhiUD(alpha1Eqn.flux());
alphaPhi = talphaPhiUD();
tmp<surfaceScalarField> talphaPhi1UD(alpha1Eqn.flux());
alphaPhi10 = talphaPhi1UD();
if (alphaApplyPrevCorr && talphaPhiCorr0.valid())
if (alphaApplyPrevCorr && talphaPhi1Corr0.valid())
{
Info<< "Applying the previous iteration compression flux" << endl;
MULES::correct
(
porosity,
alpha1,
alphaPhi,
talphaPhiCorr0.ref(),
alphaPhi10,
talphaPhi1Corr0.ref(),
zeroField(), zeroField(),
1, 0
);
alphaPhi += talphaPhiCorr0();
alphaPhi10 += talphaPhi1Corr0();
}
// Cache the upwind-flux
talphaPhiCorr0 = talphaPhiUD;
talphaPhi1Corr0 = talphaPhi1UD;
alpha2 = 1.0 - alpha1;
......@@ -133,61 +157,66 @@
for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
{
#include "alphaSuSp.H"
surfaceScalarField phir(phic*mixture.nHatf());
alphaPhiUn =
tmp<surfaceScalarField> talphaPhi1Un
(
fvc::flux
(
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();
}
phiCN(),
cnCoeff*alpha1 + (1.0 - cnCoeff)*alpha1.oldTime(),
alphaScheme
)
+ fvc::flux
(
-fvc::flux(-phir, alpha2, alpharScheme),
alpha1,
alpharScheme
)
);
if (MULESCorr)
{
tmp<surfaceScalarField> talphaPhiCorr(alphaPhiUn - alphaPhi);
tmp<surfaceScalarField> talphaPhi1Corr(talphaPhi1Un() - alphaPhi10);
volScalarField alpha10("alpha10", alpha1);
MULES::correct(alpha1, alphaPhiUn, talphaPhiCorr.ref(), 1, 0);
MULES::correct
(
geometricOneField(),
alpha1,
talphaPhi1Un(),
talphaPhi1Corr.ref(),
Sp,
(-Sp*alpha1)(),
1,
0
);
// Under-relax the correction for all but the 1st corrector
if (aCorr == 0)
{
alphaPhi += talphaPhiCorr();
alphaPhi10 += talphaPhi1Corr();
}
else
{
alpha1 = 0.5*alpha1 + 0.5*alpha10;
alphaPhi += 0.5*talphaPhiCorr();
alphaPhi10 += 0.5*talphaPhi1Corr();
}
}
else
{
alphaPhi = alphaPhiUn;
alphaPhi10 = talphaPhi1Un;
MULES::explicitSolve
(
porosity,
alpha1,
phiCN,
alphaPhi,
zeroField(),
zeroField(),
alphaPhi10,
Sp,
(Su + divU*min(alpha1(), scalar(1)))(),
1,
0
);
......@@ -200,27 +229,37 @@
if (alphaApplyPrevCorr && MULESCorr)
{
talphaPhiCorr0 = alphaPhi - talphaPhiCorr0;
talphaPhi1Corr0 = alphaPhi10 - talphaPhi1Corr0;
talphaPhi1Corr0.ref().rename("alphaPhi1Corr0");
}
else
{
talphaPhi1Corr0.clear();
}
#include "rhofs.H"
if
(
word(mesh.ddtScheme("ddt(rho,U)"))
== fv::EulerDdtScheme<vector>::typeName
|| word(mesh.ddtScheme("ddt(rho,U)"))
== fv::localEulerDdtScheme<vector>::typeName
)
{
rhoPhi = alphaPhi*(rho1 - rho2) + phiCN*rho2;
rhoPhi = alphaPhi10*(rho1f - rho2f) + phiCN*rho2f;
}
else
{
if (ocCoeff > 0)
{
// Calculate the end-of-time-step alpha flux
alphaPhi = (alphaPhi - (1.0 - cnCoeff)*alphaPhi.oldTime())/cnCoeff;
alphaPhi10 =
(alphaPhi10 - (1.0 - cnCoeff)*alphaPhi10.oldTime())/cnCoeff;
}
// Calculate the end-of-time-step mass flux
rhoPhi = alphaPhi*(rho1 - rho2) + phi*rho2;
rhoPhi = alphaPhi10*(rho1f - rho2f) + phi*rho2f;
}
Info<< "Phase-1 volume fraction = "
......
......@@ -121,20 +121,6 @@ if (p_rgh.needReference())
mesh.setFluxRequired(p_rgh.name());
mesh.setFluxRequired(alpha1.name());
// MULES flux from previous time-step
surfaceScalarField alphaPhi
(
IOobject
(
"alphaPhi",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
phi*fvc::interpolate(alpha1)
);
// MULES compressed flux is registered in case scalarTransport FO needs it.
surfaceScalarField alphaPhiUn
(
......@@ -164,7 +150,4 @@ volScalarField porosity
dimensionedScalar("porosity0", dimless, 1)
);
// MULES Correction
tmp<surfaceScalarField> talphaPhiCorr0;
#include "createMRF.H"
......@@ -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
......@@ -22,7 +22,7 @@ License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Application
interFoam
interPorousFoam
Group
grpMultiphaseSolvers
......@@ -51,7 +51,6 @@ Description
#include "pimpleControl.H"
#include "fvOptions.H"
#include "CorrectPhi.H"
#include "localEulerDdtScheme.H"
#include "fvcSmooth.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......@@ -67,6 +66,7 @@ int main(int argc, char *argv[])
#include "createTimeControls.H"
#include "initContinuityErrs.H"
#include "createFields.H"
#include "createAlphaFluxes.H"
#include "createFvOptions.H"
#include "correctPhi.H"
......@@ -110,6 +110,11 @@ int main(int argc, char *argv[])
mixture.correct();
if (pimple.frozenFlow())
{
continue;
}
#include "UEqn.H"
// --- Pressure corrector loop
......
......@@ -30,7 +30,9 @@
{
fvScalarMatrix p_rghEqn
(
fvm::laplacian(rAUf*fvc::interpolate(porosity), p_rgh) == fvc::div(phiHbyA)
fvm::laplacian(rAUf*fvc::interpolate(porosity), p_rgh)
==
fvc::div(phiHbyA)
);
p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
......
mangroves/mangrovesSource.C
multiphasePorosity/multiphasePorositySource.C
LIB = $(FOAM_USER_LIBBIN)/libfvOptions-extra
EXE_INC = \
-DFULLDEBUG -g -O3 \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/solidThermo/lnInclude \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/compressible/lnInclude \
-I$(LIB_SRC)/transportModels/incompressible/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
-I$(LIB_SRC)/fvOptions/lnInclude \
LIB_LIBS = \
-lfiniteVolume \
-lsampling \
-lmeshTools \
-lturbulenceModels \
-lincompressibleTurbulenceModels \
-lcompressibleTurbulenceModels \
-lfvOptions
#!/bin/sh
cd ${0%/*} || exit 1 # run from this directory
wmake interFoam-new
interFoam-new.C
EXE = $(FOAM_USER_APPBIN)/interFoam-new
/*---------------------------------------------------------------------------*\
========= |
\\ / 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;
// ************************************************************************* //
CorrectPhi
(
U,
phi,
p_rgh,
dimensionedScalar("rAUf", dimTime/rho.dimensions(), 1),
geometricZeroField(),
pimple
);
#include "continuityErrs.H"
/*---------------------------------------------------------------------------*\
========= |
\\ / 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;
}
// ************************************************************************* //
{
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)
);