ENH: overset: Initial release of overset capability.

Adds overset discretisation to selected physics:
- diffusion : overLaplacianDyMFoam
- incompressible steady : overSimpleFoam
- incompressible transient : overPimpleDyMFoam
- compressible transient: overRhoPimpleDyMFoam
- two-phase VOF: overInterDyMFoam

The overset method chosen is a parallel, fully implicit implementation
whereby the interpolation (from donor to acceptor) is inserted as an
adapted discretisation on the donor cells, such that the resulting matrix
can be solved using the standard linear solvers.

Above solvers come with a set of tutorials, showing how to create and set-up
simple simulations from scratch.
This commit is contained in:
mattijs
2017-06-14 09:51:02 +01:00
parent 69deec2e1c
commit fd665b4a3c
374 changed files with 29369 additions and 579 deletions

View File

@ -0,0 +1,53 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
CourantNo
Description
Calculates and outputs the mean and maximum Courant Numbers.
\*---------------------------------------------------------------------------*/
scalar CoNum = 0.0;
scalar meanCoNum = 0.0;
if (mesh.nInternalFaces())
{
surfaceScalarField phiMask(localMin<scalar>(mesh).interpolate(cellMask));
scalarField sumPhi
(
fvc::surfaceSum(mag(phiMask*phi))().internalField()
);
CoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
meanCoNum =
0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();
}
Info<< "Courant Number mean: " << meanCoNum
<< " max: " << CoNum << endl;
// ************************************************************************* //

View File

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

View File

@ -0,0 +1,31 @@
EXE_INC = \
-I. \
-I.. \
-I../../VoF \
-I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/lnInclude \
-I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
-I$(LIB_SRC)/transportModels/immiscibleIncompressibleTwoPhaseMixture/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/dynamicMesh/lnInclude \
-I$(LIB_SRC)/dynamicFvMesh/lnInclude \
-I$(FOAM_SOLVERS)/incompressible/pimpleFoam/overPimpleDyMFoam \
-I$(LIB_SRC)/overset/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude
EXE_LIBS = \
-limmiscibleIncompressibleTwoPhaseMixture \
-lturbulenceModels \
-lincompressibleTurbulenceModels \
-lfiniteVolume \
-ldynamicMesh \
-ldynamicFvMesh \
-ltopoChangerFvMesh \
-loverset \
-lfvOptions \
-lsampling \
-lwaveModels

View File

@ -0,0 +1,33 @@
MRF.correctBoundaryVelocity(U);
fvVectorMatrix UEqn
(
fvm::ddt(rho, U) + fvm::div(rhoPhi, U)
+ MRF.DDt(rho, U)
+ turbulence->divDevRhoReff(rho, U)
==
fvOptions(rho, U)
);
UEqn.relax();
fvOptions.constrain(UEqn);
if (pimple.momentumPredictor())
{
solve
(
UEqn
==
cellMask*fvc::reconstruct
(
(
mixture.surfaceTensionForce()
- ghf*fvc::snGrad(rho)
- fvc::snGrad(p_rgh)
) * mesh.magSf()
)
);
fvOptions.correct(U);
}

View File

@ -0,0 +1,48 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd.
-------------------------------------------------------------------------------
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
continuityErrs
Description
Calculates and prints the continuity errors.
\*---------------------------------------------------------------------------*/
{
volScalarField contErr(interpolatedCells*cellMask*fvc::div(phi));
scalar sumLocalContErr = runTime.deltaTValue()*
mag(contErr)().weightedAverage(mesh.V()).value();
scalar globalContErr = runTime.deltaTValue()*
contErr.weightedAverage(mesh.V()).value();
cumulativeContErr += globalContErr;
Info<< "time step continuity errors : sum local = " << sumLocalContErr
<< ", global = " << globalContErr
<< ", cumulative = " << cumulativeContErr
<< endl;
}
// ************************************************************************* //

View File

@ -0,0 +1,125 @@
{
if (mesh.changing())
{
volVectorField::Boundary& bfld = U.boundaryFieldRef();
forAll(bfld, patchi)
{
if (bfld[patchi].fixesValue())
{
bfld[patchi].initEvaluate();
}
}
surfaceScalarField::Boundary& phiBfld = phi.boundaryFieldRef();
forAll(bfld, patchi)
{
if (bfld[patchi].fixesValue())
{
bfld[patchi].evaluate();
phiBfld[patchi] =
bfld[patchi]
& mesh.Sf().boundaryField()[patchi];
}
}
}
wordList pcorrTypes
(
p_rgh.boundaryField().size(),
zeroGradientFvPatchScalarField::typeName
);
for (label i=0; i<p_rgh.boundaryField().size(); i++)
{
if (p_rgh.boundaryField()[i].fixesValue())
{
pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
}
}
volScalarField pcorr
(
IOobject
(
"pcorr",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedScalar("pcorr", p_rgh.dimensions(), 0.0),
pcorrTypes
);
if (pcorr.needReference())
{
fvc::makeRelative(phi, U);
adjustPhi(phi, U, pcorr);
fvc::makeAbsolute(phi, U);
}
mesh.setFluxRequired(pcorr.name());
dimensionedScalar rAUf("rAUf", dimTime/rho.dimensions(), 1.0);
const cellCellStencilObject& overlap = Stencil::New(mesh);
const labelList& cellTypes = overlap.cellTypes();
const labelIOList& zoneIDs = overlap.zoneID();
while (pimple.correctNonOrthogonal())
{
label nZones = gMax(zoneIDs)+1;
//label refCellI2 = -1;
labelList refCells(nZones, -1);
labelList refZones(nZones, -1);
forAll(zoneIDs, cellI)
{
label zoneId = zoneIDs[cellI];
if
(
refCells[zoneId] == -1
&& cellTypes[cellI] == cellCellStencil::CALCULATED
&& refZones[zoneId] == -1
)
{
refCells[zoneId] = cellI;
refZones[zoneId] = zoneId;
}
}
fvScalarMatrix pcorrEqn
(
fvm::laplacian(rAUf, pcorr) == fvc::div(phi)
);
//pcorrEqn.setReference(refCellI2, 0, true);
scalarList values(nZones, 0.0);
pcorrEqn.setReferences(refCells, values, true);
const dictionary& d = mesh.solver
(
pcorr.select
(
pimple.finalInnerIter()
)
);
//Bypass virtual layer
mesh.fvMesh::solve(pcorrEqn, d);
if (pimple.finalNonOrthogonalIter())
{
phi -= pcorrEqn.flux();
}
}
if (runTime.outputTime())
{
volScalarField("contPhiPcorr", fvc::div(phi)).write();
pcorr.write();
}
}

View File

@ -0,0 +1,25 @@
bool correctPhi
(
pimple.dict().lookupOrDefault<Switch>("correctPhi", true)
);
bool checkMeshCourantNo
(
pimple.dict().lookupOrDefault<Switch>("checkMeshCourantNo", false)
);
bool moveMeshOuterCorrectors
(
pimple.dict().lookupOrDefault<Switch>("moveMeshOuterCorrectors", false)
);
bool massFluxInterpolation
(
pimple.dict().lookupOrDefault("massFluxInterpolation", false)
);
bool ddtCorr
(
pimple.dict().lookupOrDefault("ddtCorr", true)
);

View File

@ -0,0 +1,171 @@
#include "createRDeltaT.H"
Info<< "Reading field p_rgh\n" << endl;
volScalarField p_rgh
(
IOobject
(
"p_rgh",
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"
//- Overset specific
// Add solver-specific interpolations
{
dictionary oversetDict;
oversetDict.add("U", true);
oversetDict.add("p", true);
oversetDict.add("HbyA", true);
oversetDict.add("p_rgh", true);
oversetDict.add("alpha1", true);
oversetDict.add("minGradP", true);
const_cast<dictionary&>
(
mesh.schemesDict()
).add
(
"oversetInterpolationRequired",
oversetDict,
true
);
}
// Mask field for zeroing out contributions on hole cells
#include "createCellMask.H"
// Create bool field with interpolated cells
#include "createInterpolatedCells.H"
Info<< "Reading transportProperties\n" << endl;
immiscibleIncompressibleTwoPhaseMixture mixture(U, phi);
volScalarField& alpha1(mixture.alpha1());
volScalarField& alpha2(mixture.alpha2());
const dimensionedScalar& rho1 = mixture.rho1();
const dimensionedScalar& rho2 = mixture.rho2();
// Need to store rho for ddt(rho, U)
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT
),
alpha1*rho1 + alpha2*rho2
);
rho.oldTime();
// Mass flux
surfaceScalarField rhoPhi
(
IOobject
(
"rhoPhi",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
fvc::interpolate(rho)*phi
);
// Construct incompressible turbulence model
autoPtr<incompressible::turbulenceModel> turbulence
(
incompressible::turbulenceModel::New(U, phi, mixture)
);
#include "readGravitationalAcceleration.H"
#include "readhRef.H"
#include "gh.H"
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
p_rgh + rho*gh
);
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell
(
p,
p_rgh,
pimple.dict(),
pRefCell,
pRefValue
);
if (p_rgh.needReference())
{
p += dimensionedScalar
(
"p",
p.dimensions(),
pRefValue - getRefCellValue(p, pRefCell)
);
p_rgh = p - rho*gh;
}
mesh.setFluxRequired(p_rgh.name());
mesh.setFluxRequired(alpha1.name());
// MULES compressed flux is registered in case scalarTransport FO needs it.
surfaceScalarField alphaPhiUn
(
IOobject
(
"alphaPhiUn",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedScalar("zero", phi.dimensions(), 0.0)
);
#include "createMRF.H"

View File

@ -0,0 +1,209 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2016-2017 OpenCFD Ltd.
-------------------------------------------------------------------------------
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/>.
Application
overInterDyMFoam
Group
grpMultiphaseSolvers grpMovingMeshSolvers
Description
Solver for 2 incompressible, isothermal immiscible fluids using a VOF
(volume of fluid) phase-fraction based interface capturing approach,
with optional mesh motion and mesh topology changes including adaptive
re-meshing.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "dynamicFvMesh.H"
#include "CMULES.H"
#include "EulerDdtScheme.H"
#include "localEulerDdtScheme.H"
#include "CrankNicolsonDdtScheme.H"
#include "subCycle.H"
#include "immiscibleIncompressibleTwoPhaseMixture.H"
#include "turbulentTransportModel.H"
#include "pimpleControl.H"
#include "fvOptions.H"
#include "CorrectPhi.H"
#include "fvcSmooth.H"
#include "cellCellStencilObject.H"
#include "localMin.H"
#include "interpolationCellPoint.H"
#include "transform.H"
#include "fvMeshSubset.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "postProcess.H"
#include "setRootCase.H"
#include "createTime.H"
#include "createDynamicFvMesh.H"
#include "initContinuityErrs.H"
#include "createControl.H"
#include "createTimeControls.H"
#include "createDyMControls.H"
#include "createFields.H"
#include "createAlphaFluxes.H"
#include "createFvOptions.H"
volScalarField rAU
(
IOobject
(
"rAU",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("rAUf", dimTime/rho.dimensions(), 1.0)
);
#include "correctPhi.H"
#include "createUf.H"
turbulence->validate();
if (!LTS)
{
#include "CourantNo.H"
#include "setInitialDeltaT.H"
}
#include "setCellMask.H"
#include "setInterpolatedCells.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.run())
{
#include "readControls.H"
if (LTS)
{
#include "setRDeltaT.H"
}
else
{
#include "CourantNo.H"
#include "alphaCourantNo.H"
#include "setDeltaT.H"
}
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())
{
if (pimple.firstIter() || moveMeshOuterCorrectors)
{
scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime();
mesh.update();
if (mesh.changing())
{
Info<< "Execution time for mesh.update() = "
<< runTime.elapsedCpuTime() - timeBeforeMeshUpdate
<< " s" << endl;
// Do not apply previous time-step mesh compression flux
// if the mesh topology changed
if (mesh.topoChanging())
{
talphaPhiCorr0.clear();
}
gh = (g & mesh.C()) - ghRef;
ghf = (g & mesh.Cf()) - ghRef;
// Update cellMask field for blocking out hole cells
#include "setCellMask.H"
#include "setInterpolatedCells.H"
}
if ((mesh.changing() && correctPhi) || mesh.topoChanging())
{
// Calculate absolute flux from the mapped surface velocity
// Note: temporary fix until mapped Uf is assessed
Uf = fvc::interpolate(U);
// Calculate absolute flux from the mapped surface velocity
phi = mesh.Sf() & Uf;
#include "correctPhi.H"
// Make the flux relative to the mesh motion
fvc::makeRelative(phi, U);
mixture.correct();
}
if (mesh.changing() && checkMeshCourantNo)
{
#include "meshCourantNo.H"
}
}
#include "alphaControls.H"
#include "alphaEqnSubCycle.H"
mixture.correct();
#include "UEqn.H"
// --- Pressure corrector loop
while (pimple.correct())
{
#include "pEqn.H"
}
if (pimple.turbCorr())
{
turbulence->correct();
}
}
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //

View File

@ -0,0 +1,96 @@
{
rAU = 1.0/UEqn.A();
surfaceScalarField faceMask(localMin<scalar>(mesh).interpolate(cellMask));
surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
volVectorField HbyA("HbyA", U);
//HbyA = rAU*UEqn.H();
HbyA = constrainHbyA(rAU*UEqn.H(), U, p_rgh);
if (massFluxInterpolation)
{
#include "interpolatedFaces.H"
}
surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA));
if (ddtCorr)
{
phiHbyA += fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, Uf);
}
MRF.makeRelative(phiHbyA);
if (p_rgh.needReference())
{
fvc::makeRelative(phiHbyA, U);
adjustPhi(phiHbyA, U, p_rgh);
fvc::makeAbsolute(phiHbyA, U);
}
surfaceScalarField phig
(
(
mixture.surfaceTensionForce()
- ghf*fvc::snGrad(rho)
)*faceMask*rAUf*mesh.magSf()
);
phiHbyA += phig;
// Update the pressure BCs to ensure flux consistency
constrainPressure(p_rgh, U, phiHbyA, rAUf, MRF);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix p_rghEqn
(
fvm::laplacian(faceMask*rAUf, p_rgh) == fvc::div(phiHbyA)
);
p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
if (pimple.finalNonOrthogonalIter())
{
phi = phiHbyA - p_rghEqn.flux();
p_rgh.relax();
// Reconstruct body forces (-grad(p) and gh etc)
volVectorField minGradP
(
"minGradP",
fvc::reconstruct((phig - p_rghEqn.flux())/rAUf)
);
U = HbyA + rAU*cellMask*minGradP;
U.correctBoundaryConditions();
fvOptions.correct(U);
}
}
#include "continuityErrs.H"
{
Uf = fvc::interpolate(U);
surfaceVectorField n(mesh.Sf()/mesh.magSf());
Uf += n*(phi/mesh.magSf() - (n & Uf));
}
// Make the fluxes relative to the mesh motion
fvc::makeRelative(phi, U);
p == p_rgh + rho*gh;
if (p_rgh.needReference())
{
p += dimensionedScalar
(
"p",
p.dimensions(),
pRefValue - getRefCellValue(p, pRefCell)
);
p_rgh = p - rho*gh;
}
}

View File

@ -0,0 +1,15 @@
#include "readTimeControls.H"
correctPhi = pimple.dict().lookupOrDefault<Switch>("correctPhi", true);
checkMeshCourantNo =
pimple.dict().lookupOrDefault<Switch>("checkMeshCourantNo", false);
moveMeshOuterCorrectors =
pimple.dict().lookupOrDefault<Switch>("moveMeshOuterCorrectors", false);
massFluxInterpolation =
pimple.dict().lookupOrDefault("massFluxInterpolation", false);
ddtCorr =
pimple.dict().lookupOrDefault("ddtCorr", true);