Skip to content
Snippets Groups Projects
Commit 985a18df authored by sergio's avatar sergio
Browse files

Adding overPotentialFoam and overRhoSimpleFoam and tutorials

parent a6d88f2f
No related merge requests found
Showing
with 978 additions and 0 deletions
potentialFoam.C
EXE = $(FOAM_APPBIN)/overPotentialFoam
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/dynamicFvMesh/lnInclude \
-I$(LIB_SRC)/overset/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude
EXE_LIBS = \
-lfiniteVolume \
-lmeshTools \
-lsampling \
-loverset
const dictionary& potentialFlow
(
mesh.solutionDict().subDict("potentialFlow")
);
const int nNonOrthCorr
(
potentialFlow.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0)
);
Info<< "Reading velocity field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
// Initialise the velocity internal field to zero
U = dimensionedVector("0", U.dimensions(), Zero);
surfaceScalarField phi
(
IOobject
(
"phi",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
fvc::flux(U)
);
if (args.optionFound("initialiseUBCs"))
{
U.correctBoundaryConditions();
phi = fvc::flux(U);
}
// Construct a pressure field
// If it is available read it otherwise construct from the velocity BCs
// converting fixed-value BCs to zero-gradient and vice versa.
word pName("p");
// Update name of the pressure field from the command-line option
args.optionReadIfPresent("pName", pName);
// Infer the pressure BCs from the velocity
wordList pBCTypes
(
U.boundaryField().size(),
fixedValueFvPatchScalarField::typeName
);
forAll(U.boundaryField(), patchi)
{
if (U.boundaryField()[patchi].fixesValue())
{
pBCTypes[patchi] = zeroGradientFvPatchScalarField::typeName;
}
}
// Note that registerObject is false for the pressure field. The pressure
// field in this solver doesn't have a physical value during the solution.
// It shouldn't be looked up and used by sub models or boundary conditions.
Info<< "Constructing pressure field " << pName << nl << endl;
volScalarField p
(
IOobject
(
pName,
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE,
false
),
mesh,
dimensionedScalar(pName, sqr(dimVelocity), 0),
pBCTypes
);
// Infer the velocity potential BCs from the pressure
wordList PhiBCTypes
(
p.boundaryField().size(),
zeroGradientFvPatchScalarField::typeName
);
forAll(p.boundaryField(), patchi)
{
if (p.boundaryField()[patchi].fixesValue())
{
PhiBCTypes[patchi] = fixedValueFvPatchScalarField::typeName;
}
}
Info<< "Constructing velocity potential field Phi\n" << endl;
volScalarField Phi
(
IOobject
(
"Phi",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE
),
mesh,
dimensionedScalar("Phi", dimLength*dimVelocity, 0),
PhiBCTypes
);
label PhiRefCell = 0;
scalar PhiRefValue = 0;
setRefCell
(
Phi,
potentialFlow.dict(),
PhiRefCell,
PhiRefValue
);
mesh.setFluxRequired(Phi.name());
#include "createMRF.H"
// Add overset specific interpolations
{
dictionary oversetDict;
oversetDict.add("Phi", true);
oversetDict.add("U", 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"
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2017 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 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
potentialFoam
Group
grpBasicSolvers
Description
Potential flow solver which solves for the velocity potential, to
calculate the flux-field, from which the velocity field is obtained by
reconstructing the flux.
\heading Solver details
The potential flow solution is typically employed to generate initial fields
for full Navier-Stokes codes. The flow is evolved using the equation:
\f[
\laplacian \Phi = \div(\vec{U})
\f]
Where:
\vartable
\Phi | Velocity potential [m2/s]
\vec{U} | Velocity [m/s]
\endvartable
The corresponding pressure field could be calculated from the divergence
of the Euler equation:
\f[
\laplacian p + \div(\div(\vec{U}\otimes\vec{U})) = 0
\f]
but this generates excessive pressure variation in regions of large
velocity gradient normal to the flow direction. A better option is to
calculate the pressure field corresponding to velocity variation along the
stream-lines:
\f[
\laplacian p + \div(\vec{F}\cdot\div(\vec{U}\otimes\vec{U})) = 0
\f]
where the flow direction tensor \f$\vec{F}\f$ is obtained from
\f[
\vec{F} = \hat{\vec{U}}\otimes\hat{\vec{U}}
\f]
\heading Required fields
\plaintable
U | Velocity [m/s]
\endplaintable
\heading Optional fields
\plaintable
p | Kinematic pressure [m2/s2]
Phi | Velocity potential [m2/s]
| Generated from p (if present) or U if not present
\endplaintable
\heading Options
\plaintable
-writep | write the Euler pressure
-writePhi | Write the final velocity potential
-initialiseUBCs | Update the velocity boundaries before solving for Phi
\endplaintable
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "pisoControl.H"
#include "dynamicFvMesh.H"
#include "cellCellStencilObject.H"
#include "localMin.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
argList::addOption
(
"pName",
"pName",
"Name of the pressure field"
);
argList::addBoolOption
(
"initialiseUBCs",
"Initialise U boundary conditions"
);
argList::addBoolOption
(
"writePhi",
"Write the final velocity potential field"
);
argList::addBoolOption
(
"writep",
"Calculate and write the Euler pressure field"
);
argList::addBoolOption
(
"withFunctionObjects",
"execute functionObjects"
);
#include "setRootCase.H"
#include "createTime.H"
#include "createNamedDynamicFvMesh.H"
pisoControl potentialFlow(mesh, "potentialFlow");
#include "createFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< nl << "Calculating potential flow" << endl;
mesh.update();
surfaceScalarField faceMask(localMin<scalar>(mesh).interpolate(cellMask));
// Since solver contains no time loop it would never execute
// function objects so do it ourselves
runTime.functionObjects().start();
MRF.makeRelative(phi);
adjustPhi(phi, U, p);
// Non-orthogonal velocity potential corrector loop
while (potentialFlow.correct())
{
phi = fvc::flux(U);
while (potentialFlow.correctNonOrthogonal())
{
fvScalarMatrix PhiEqn
(
fvm::laplacian(faceMask, Phi)
==
fvc::div(phi)
);
PhiEqn.setReference(PhiRefCell, PhiRefValue);
PhiEqn.solve();
if (potentialFlow.finalNonOrthogonalIter())
{
phi -= PhiEqn.flux();
}
}
MRF.makeAbsolute(phi);
Info<< "Continuity error = "
<< mag(fvc::div(phi))().weightedAverage(mesh.V()).value()
<< endl;
U = fvc::reconstruct(phi);
U.correctBoundaryConditions();
Info<< "Interpolated velocity error = "
<< (sqrt(sum(sqr(fvc::flux(U) - phi)))/sum(mesh.magSf())).value()
<< endl;
}
// Write U and phi
U.write();
phi.write();
// Optionally write Phi
if (args.optionFound("writePhi"))
{
Phi.write();
}
// Calculate the pressure field from the Euler equation
if (args.optionFound("writep"))
{
Info<< nl << "Calculating approximate pressure field" << endl;
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell
(
p,
potentialFlow.dict(),
pRefCell,
pRefValue
);
// Calculate the flow-direction filter tensor
volScalarField magSqrU(magSqr(U));
volSymmTensorField F(sqr(U)/(magSqrU + SMALL*average(magSqrU)));
// Calculate the divergence of the flow-direction filtered div(U*U)
// Filtering with the flow-direction generates a more reasonable
// pressure distribution in regions of high velocity gradient in the
// direction of the flow
volScalarField divDivUU
(
fvc::div
(
F & fvc::div(phi, U),
"div(div(phi,U))"
)
);
// Solve a Poisson equation for the approximate pressure
while (potentialFlow.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
fvm::laplacian(p) + divDivUU
);
pEqn.setReference(pRefCell, pRefValue);
pEqn.solve();
}
p.write();
}
runTime.functionObjects().end();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //
rhoSimpleFoam.C
EXE = $(FOAM_APPBIN)/overRhoSimpleFoam
EXE_INC = \
-I.. \
-I$(LIB_SRC)/transportModels/compressible/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
-I$(LIB_SRC)/finiteVolume/cfdTools \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude \
-I$(LIB_SRC)/dynamicFvMesh/lnInclude \
-I$(LIB_SRC)/dynamicMesh/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/overset/lnInclude
EXE_LIBS = \
-lcompressibleTransportModels \
-lfluidThermophysicalModels \
-lspecie \
-lturbulenceModels \
-lcompressibleTurbulenceModels \
-lfiniteVolume \
-lsampling \
-lmeshTools \
-lfvOptions \
-loverset \
-ldynamicFvMesh \
-ltopoChangerFvMesh
// Solve the Momentum equation
MRF.correctBoundaryVelocity(U);
tmp<fvVectorMatrix> tUEqn
(
fvm::div(phi, U)
+ MRF.DDt(rho, U)
+ turbulence->divDevRhoReff(U)
==
fvOptions(rho, U)
);
fvVectorMatrix& UEqn = tUEqn.ref();
UEqn.relax();
fvOptions.constrain(UEqn);
if (simple.momentumPredictor())
{
solve(UEqn == -cellMask*fvc::grad(p));
}
fvOptions.correct(U);
const volScalarField& psi = thermo.psi();
Info<< "Reading thermophysical properties\n" << endl;
autoPtr<fluidThermo> pThermo
(
fluidThermo::New(mesh)
);
fluidThermo& thermo = pThermo();
thermo.validate(args.executable(), "h", "e");
volScalarField& p = thermo.p();
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
thermo.rho()
);
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
#include "compressibleCreatePhi.H"
pressureControl pressureControl(p, rho, simple.dict());
mesh.setFluxRequired(p.name());
Info<< "Creating turbulence model\n" << endl;
autoPtr<compressible::turbulenceModel> turbulence
(
compressible::turbulenceModel::New
(
rho,
U,
phi,
thermo
)
);
dimensionedScalar initialMass = fvc::domainIntegrate(rho);
#include "createMRF.H"
//- Overset specific
// Add solver-specific interpolations
{
dictionary oversetDict;
oversetDict.add("U", true);
oversetDict.add("p", true);
oversetDict.add("HbyA", true);
oversetDict.add("grad(p)", true);
oversetDict.add("rho", true);
const_cast<dictionary&>
(
mesh.schemesDict()
).add
(
"oversetInterpolationRequired",
oversetDict,
true
);
}
// Mask field for zeroing out contributions on hole cells
#include "createCellMask.H"
#include "createInterpolatedCells.H"
bool adjustFringe
(
simple.dict().lookupOrDefault("oversetAdjustPhi", false)
);
Info<< "Create dynamic mesh for time = "
<< runTime.timeName() << nl << endl;
autoPtr<dynamicFvMesh> meshPtr
(
dynamicFvMesh::New
(
IOobject
(
dynamicFvMesh::defaultRegion,
runTime.timeName(),
runTime,
IOobject::MUST_READ
)
)
);
dynamicFvMesh& mesh = meshPtr();
// Calculate initial mesh-to-mesh mapping. Note that this should be
// done under the hood, e.g. as a MeshObject
mesh.update();
{
surfaceScalarField faceMask(localMin<scalar>(mesh).interpolate(cellMask));
volScalarField rAU(1.0/UEqn.A());
surfaceScalarField rhorAUf("rhorAUf", faceMask*fvc::interpolate(rho*rAU));
volVectorField HbyA("HbyA", U);
HbyA = constrainHbyA(cellMask*rAU*UEqn.H(), U, p);
tUEqn.clear();
bool closedVolume = false;
surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(rho)*fvc::flux(HbyA));
MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
// Update the pressure BCs to ensure flux consistency
constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
if (simple.transonic())
{
surfaceScalarField phid
(
"phid",
(fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
);
phiHbyA -= fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
while (simple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
fvc::div(phiHbyA)
+ fvm::div(phid, p)
- fvm::laplacian(rhorAUf, p)
==
fvOptions(psi, p, rho.name())
);
// Relax the pressure equation to ensure diagonal-dominance
pEqn.relax();
pEqn.setReference
(
pressureControl.refCell(),
pressureControl.refValue()
);
pEqn.solve();
if (simple.finalNonOrthogonalIter())
{
phi = phiHbyA + pEqn.flux();
}
}
}
else
{
closedVolume = adjustPhi(phiHbyA, U, p);
if (adjustFringe)
{
oversetAdjustPhi(phiHbyA, U);
}
while (simple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
fvc::div(phiHbyA)
- fvm::laplacian(rhorAUf, p)
==
fvOptions(psi, p, rho.name())
);
pEqn.setReference
(
pressureControl.refCell(),
pressureControl.refValue()
);
pEqn.solve();
if (simple.finalNonOrthogonalIter())
{
phi = phiHbyA + pEqn.flux();
}
}
}
#include "incompressible/continuityErrs.H"
// Explicitly relax pressure for momentum corrector
p.relax();
volVectorField gradP(fvc::grad(p));
U = HbyA - rAU*cellMask*gradP;
U.correctBoundaryConditions();
fvOptions.correct(U);
bool pLimited = pressureControl.limit(p);
// For closed-volume cases adjust the pressure and density levels
// to obey overall mass continuity
if (closedVolume)
{
p += (initialMass - fvc::domainIntegrate(psi*p))
/fvc::domainIntegrate(psi);
}
if (pLimited || closedVolume)
{
p.correctBoundaryConditions();
}
rho = thermo.rho();
if (!simple.transonic())
{
rho.relax();
}
}
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2017 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 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
overRhoSimpleFoam
Group
grpCompressibleSolvers
Description
Overset steady-state solver for turbulent flow of compressible fluids.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "dynamicFvMesh.H"
#include "fluidThermo.H"
#include "turbulentFluidThermoModel.H"
#include "simpleControl.H"
#include "pressureControl.H"
#include "fvOptions.H"
#include "cellCellStencilObject.H"
#include "localMin.H"
#include "oversetAdjustPhi.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#define CREATE_MESH createUpdatedDynamicFvMesh.H
#include "postProcess.H"
#include "setRootCase.H"
#include "createTime.H"
#include "createUpdatedDynamicFvMesh.H"
#include "createControl.H"
#include "createFields.H"
#include "createFieldRefs.H"
#include "createFvOptions.H"
#include "initContinuityErrs.H"
turbulence->validate();
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (simple.loop())
{
Info<< "Time = " << runTime.timeName() << nl << endl;
// Pressure-velocity SIMPLE corrector
#include "UEqn.H"
#include "EEqn.H"
#include "pEqn.H"
turbulence->correct();
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //
#!/bin/sh
cd ${0%/*} || exit 1 # run from this directory
. $WM_PROJECT_DIR/bin/tools/CleanFunctions
# Extrude mesh around cylinder
(cd cylinderAndBackground && ./Allclean)
# Add background mesh
(cd cylinderMesh && foamCleanTutorials)
# ----------------------------------------------------------------- end-of-file
#!/bin/sh
. $WM_PROJECT_DIR/bin/tools/RunFunctions
# Extrude mesh around cylinder
(cd cylinderMesh && ./Allrun.pre)
# Add background mesh
(cd cylinderAndBackground && ./Allrun)
#!/bin/sh
. $WM_PROJECT_DIR/bin/tools/RunFunctions
# Extrude mesh around cylinder
(cd cylinderMesh && ./Allrun.pre)
# Add background mesh
(cd cylinderAndBackground && ./Allrun.pre)
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: plus |
| \\ / A nd | Web: www.OpenFOAM.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "include/initialConditions"
dimensions [0 1 -1 0 0 0 0];
internalField uniform $flowVelocity;
boundaryField
{
#includeEtc "caseDicts/setConstraintTypes"
overset
{
type overset;
}
walls
{
type movingWallVelocity;
value uniform (0 0 0);
}
inlet
{
type fixedValue;
value $internalField;
}
outlet
{
type inletOutlet;
inletValue uniform (0 0 0);
value $internalField;
}
topAndBottom
{
type zeroGradient;
}
}
// ************************************************************************* //
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: plus |
| \\ / A nd | Web: www.OpenFOAM.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
inlet
{
type fixedValue;
value $internalField;
}
// ************************************************************************* //
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: plus |
| \\ / A nd | Web: www.OpenFOAM.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
flowVelocity (1 0 0);
pressure 0;
#inputMode merge
// ************************************************************************* //
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: plus |
| \\ / A nd | Web: www.OpenFOAM.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "include/initialConditions"
dimensions [0 2 -2 0 0 0 0];
internalField uniform $pressure;
boundaryField
{
#includeEtc "caseDicts/setConstraintTypes"
overset
{
type overset;
}
wall
{
type zeroGradient;
}
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue; //calculated;
value $internalField;
}
".*"
{
type zeroGradient;
}
}
// ************************************************************************* //
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment