potentialFoam: Solve for velocity potential named Phi rather than using the pressure field for this purpose

The Phi field is read if available otherwise created automatically with
boundary conditions obtained automatically from the pressure field if
available (with optional name) otherwise inferred from the velocity
field.  Phi Laplacian scheme and solver specification are required.  See
tutorials for examples.
This commit is contained in:
Henry
2015-02-14 11:03:37 +00:00
parent 1edf292c00
commit 4129560601
16 changed files with 129 additions and 98 deletions

View File

@ -1,21 +1,4 @@
Info<< "Reading field p\n" << endl; Info<< "Reading velocity field U\n" << endl;
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
),
mesh
);
p = dimensionedScalar("zero", p.dimensions(), 0.0);
Info<< "Reading field U\n" << endl;
volVectorField U volVectorField U
( (
IOobject IOobject
@ -51,12 +34,65 @@
} }
label pRefCell = 0; // Default name for the pressure field
scalar pRefValue = 0.0; 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 BCs
wordList pBCTypes
(
U.boundaryField().size(),
fixedValueFvPatchScalarField::typeName
);
forAll(U.boundaryField(), patchi)
{
if (U.boundaryField()[patchi].fixesValue())
{
pBCTypes[patchi] = zeroGradientFvPatchScalarField::typeName;
}
}
Info<< "Constructing pressure field " << pName << nl << endl;
volScalarField p
(
IOobject
(
pName,
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE
),
mesh,
dimensionedScalar(pName, sqr(dimVelocity), 0),
pBCTypes
);
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),
p.boundaryField().types()
);
label PhiRefCell = 0;
scalar PhiRefValue = 0;
setRefCell setRefCell
( (
p, Phi,
potentialFlow, potentialFlow.dict(),
pRefCell, PhiRefCell,
pRefValue PhiRefValue
); );

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -25,29 +25,48 @@ Application
potentialFoam potentialFoam
Description Description
Simple potential flow solver which can be used to generate starting fields Potential flow solver which solves for the velocity potential
for full Navier-Stokes codes. from which the flux-field is obtained and velocity field by reconstructing
the flux.
This application is particularly useful to generate starting fields for
Navier-Stokes codes.
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
#include "fvCFD.H" #include "fvCFD.H"
#include "fvIOoptionList.H" #include "fvIOoptionList.H"
#include "pisoControl.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[]) int main(int argc, char *argv[])
{ {
argList::addBoolOption("writep", "write the final pressure field"); argList::addOption
(
"pName",
"pName",
"Name of the pressure field"
);
argList::addBoolOption argList::addBoolOption
( (
"initialiseUBCs", "initialiseUBCs",
"initialise U boundary conditions" "Initialise U boundary conditions"
);
argList::addBoolOption
(
"writePhi",
"Write the velocity potential field"
); );
#include "setRootCase.H" #include "setRootCase.H"
#include "createTime.H" #include "createTime.H"
#include "createMesh.H" #include "createMesh.H"
#include "readControls.H"
pisoControl potentialFlow(mesh, "potentialFlow");
#include "createFields.H" #include "createFields.H"
#include "createFvOptions.H" #include "createFvOptions.H"
@ -63,55 +82,47 @@ int main(int argc, char *argv[])
adjustPhi(phi, U, p); adjustPhi(phi, U, p);
// Non-orthogonal velocity potential corrector loop
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) while (potentialFlow.correctNonOrthogonal())
{ {
fvScalarMatrix pEqn fvScalarMatrix PhiEqn
( (
fvm::laplacian fvm::laplacian(dimensionedScalar("1", dimless, 1), Phi)
(
dimensionedScalar
(
"1",
dimTime/p.dimensions()*dimensionSet(0, 2, -2, 0, 0),
1
),
p
)
== ==
fvc::div(phi) fvc::div(phi)
); );
pEqn.setReference(pRefCell, pRefValue); PhiEqn.setReference(PhiRefCell, PhiRefValue);
pEqn.solve(); PhiEqn.solve();
if (nonOrth == nNonOrthCorr) if (potentialFlow.finalNonOrthogonalIter())
{ {
phi -= pEqn.flux(); phi -= PhiEqn.flux();
} }
} }
fvOptions.makeAbsolute(phi); fvOptions.makeAbsolute(phi);
Info<< "continuity error = " Info<< "Continuity error = "
<< mag(fvc::div(phi))().weightedAverage(mesh.V()).value() << mag(fvc::div(phi))().weightedAverage(mesh.V()).value()
<< endl; << endl;
U = fvc::reconstruct(phi); U = fvc::reconstruct(phi);
U.correctBoundaryConditions(); U.correctBoundaryConditions();
Info<< "Interpolated U error = " Info<< "Interpolated velocity error = "
<< (sqrt(sum(sqr((fvc::interpolate(U) & mesh.Sf()) - phi))) << (sqrt(sum(sqr((fvc::interpolate(U) & mesh.Sf()) - phi)))
/sum(mesh.magSf())).value() /sum(mesh.magSf())).value()
<< endl; << endl;
// Force the write // Write U and phi
U.write(); U.write();
phi.write(); phi.write();
if (args.optionFound("writep")) // Optionally write Phi
if (args.optionFound("writePhi"))
{ {
p.write(); Phi.write();
} }
runTime.functionObjects().end(); runTime.functionObjects().end();

View File

@ -6,26 +6,6 @@ cd ${0%/*} || exit 1 # Run from this directory
application=`getApplication` application=`getApplication`
# This case uses the #codeStream which is disabled by default. Enable for
# just this case.
MAIN_CONTROL_DICT=`foamEtcFile controlDict`
if [ -f "$MAIN_CONTROL_DICT" ]
then
echo "Modifying ${MAIN_CONTROL_DICT} to enable allowSystemOperations"
# Clean up on termination and on Ctrl-C
trap 'mv ${MAIN_CONTROL_DICT}.$$ ${MAIN_CONTROL_DICT} 2>/dev/null; exit 0' \
EXIT TERM INT
cp ${MAIN_CONTROL_DICT} ${MAIN_CONTROL_DICT}.$$
echo "Enabling allowSystemOperations in ${MAIN_CONTROL_DICT}."
sed \
-e s/"\(allowSystemOperations[ \t]*\)\([0-9]\);"/"\1 1;"/g \
${MAIN_CONTROL_DICT}.$$ > ${MAIN_CONTROL_DICT}
fi
cp -r 0.org 0 > /dev/null 2>&1 cp -r 0.org 0 > /dev/null 2>&1
runApplication blockMesh runApplication blockMesh
runApplication $application runApplication $application

View File

@ -48,7 +48,7 @@ snGradSchemes
fluxRequired fluxRequired
{ {
default no; default no;
p ; Phi ;
} }

View File

@ -17,7 +17,7 @@ FoamFile
solvers solvers
{ {
p Phi
{ {
solver PCG; solver PCG;
preconditioner DIC; preconditioner DIC;

View File

@ -48,7 +48,7 @@ snGradSchemes
fluxRequired fluxRequired
{ {
default no; default no;
p ; Phi ;
} }

View File

@ -17,7 +17,7 @@ FoamFile
solvers solvers
{ {
p Phi
{ {
solver PCG; solver PCG;
preconditioner DIC; preconditioner DIC;

View File

@ -55,7 +55,8 @@ snGradSchemes
fluxRequired fluxRequired
{ {
default no; default no;
p ; p;
Phi;
} }
wallDist wallDist

View File

@ -38,6 +38,11 @@ solvers
relTol 0; relTol 0;
}; };
Phi
{
$p;
}
"(U|nuTilda)" "(U|nuTilda)"
{ {
solver smoothSolver; solver smoothSolver;

View File

@ -53,6 +53,7 @@ fluxRequired
{ {
default no; default no;
p; p;
Phi;
} }
wallDist wallDist

View File

@ -30,6 +30,11 @@ solvers
mergeLevels 1; mergeLevels 1;
} }
Phi
{
$p;
}
U U
{ {
solver smoothSolver; solver smoothSolver;

View File

@ -7,26 +7,6 @@ cd ${0%/*} || exit 1 # Run from this directory
# Get application directory # Get application directory
application=`getApplication` application=`getApplication`
# This case uses the #codeStream which is disabled by default. Enable for
# just this case.
MAIN_CONTROL_DICT=`foamEtcFile controlDict`
if [ -f "$MAIN_CONTROL_DICT" ]
then
echo "Modifying ${MAIN_CONTROL_DICT} to enable allowSystemOperations"
# Clean up on termination and on Ctrl-C
trap 'mv ${MAIN_CONTROL_DICT}.$$ ${MAIN_CONTROL_DICT} 2>/dev/null; exit 0' \
EXIT TERM INT
cp ${MAIN_CONTROL_DICT} ${MAIN_CONTROL_DICT}.$$
echo "Enabling allowSystemOperations in ${MAIN_CONTROL_DICT}."
sed \
-e s/"\(allowSystemOperations[ \t]*\)\([0-9]\);"/"\1 1;"/g \
${MAIN_CONTROL_DICT}.$$ > ${MAIN_CONTROL_DICT}
fi
runApplication blockMesh runApplication blockMesh
runApplication topoSet runApplication topoSet
runApplication refineHexMesh c0 -overwrite runApplication refineHexMesh c0 -overwrite

View File

@ -58,6 +58,7 @@ fluxRequired
{ {
default no; default no;
p; p;
Phi;
} }
wallDist wallDist

View File

@ -68,6 +68,11 @@ solvers
relTol 0; relTol 0;
} }
Phi
{
$p;
}
"(Yi|O2|N2|H2O)" "(Yi|O2|N2|H2O)"
{ {
solver PBiCG; solver PBiCG;
@ -88,7 +93,7 @@ solvers
potentialFlow potentialFlow
{ {
// used for potentialFoam initialisation // Used for potentialFoam initialisation
nNonOrthogonalCorrectors 5; nNonOrthogonalCorrectors 5;
} }

View File

@ -59,6 +59,7 @@ fluxRequired
{ {
default no; default no;
p; p;
Phi;
} }
wallDist wallDist

View File

@ -37,6 +37,11 @@ solvers
maxIter 50; maxIter 50;
}; };
Phi
{
$p;
}
"(U|Yi|h|k|omega)" "(U|Yi|h|k|omega)"
{ {
solver smoothSolver; solver smoothSolver;