Commit 3aadb63c authored by andy's avatar andy
Browse files

ENH: Added new potentialFreeSurfaceFoam solver - single phase with wave height inclusion

parent 7db0461d
potentialFreeSurfaceFoam.C
EXE = $(FOAM_APPBIN)/potentialFreeSurfaceFoam
EXE_INC = \
-I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude
EXE_LIBS = \
-lincompressibleTransportModels \
-lincompressibleTurbulenceModel \
-lincompressibleRASModels \
-lincompressibleLESModels \
-lfiniteVolume \
-lmeshTools
tmp<fvVectorMatrix> UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
+ turbulence->divDevReff(U)
==
sources(U)
);
UEqn().relax();
sources.constrain(UEqn());
if (pimple.momentumPredictor())
{
solve(UEqn() == -fvc::grad(p_gh));
}
Info<< "Reading field p (kinematic)\n" << endl;
volScalarField p
(
IOobject
(
"p",
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"
singlePhaseTransportModel laminarTransport(U, phi);
autoPtr<incompressible::turbulenceModel> turbulence
(
incompressible::turbulenceModel::New(U, phi, laminarTransport)
);
#include "readGravitationalAcceleration.H"
Info<< "\nReading freeSurfaceProperties\n" << endl;
IOdictionary freeSurfaceProperties
(
IOobject
(
"freeSurfaceProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
)
);
word freeSurfacePatch(freeSurfaceProperties.lookup("freeSurfacePatch"));
label freeSurfacePatchI = mesh.boundaryMesh().findPatchID(freeSurfacePatch);
if (freeSurfacePatchI < 0)
{
FatalErrorIn(args.executable())
<< "Patch " << freeSurfacePatch << " not found. "
<< "Available patches are:" << mesh.boundaryMesh().names()
<< exit(FatalError);
}
Info<< "Creating field refLevel\n" << endl;
volVectorField refLevel
(
IOobject
(
"refLevel",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedVector("zero", dimLength, vector::zero)
);
refLevel.boundaryField()[freeSurfacePatchI]
== mesh.C().boundaryField()[freeSurfacePatchI];
Info<< "Creating field zeta\n" << endl;
volVectorField zeta
(
IOobject
(
"zeta",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mesh,
dimensionedVector("zero", dimLength, vector::zero)
);
Info<< "Creating field p_gh\n" << endl;
volScalarField p_gh
(
IOobject
(
"p_gh",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
// Force p_gh to be consistent with p
// Height is made relative to field 'refLevel'
p_gh = p - (g & (mesh.C() + zeta - refLevel));
label p_ghRefCell = 0;
scalar p_ghRefValue = 0.0;
setRefCell(p_gh, pimple.dict(), p_ghRefCell, p_ghRefValue);
IObasicSourceList sources(mesh);
volScalarField rAU(1.0/UEqn().A());
surfaceScalarField rAUf(rAU.name() + 'f', fvc::interpolate(rAU));
U = rAU*UEqn().H();
if (pimple.nCorrPISO() <= 1)
{
UEqn.clear();
}
phi = (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rAU, U, phi);
adjustPhi(phi, U, p_gh);
// Non-orthogonal pressure corrector loop
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix p_ghEqn
(
fvm::laplacian(rAUf, p_gh) == fvc::div(phi)
);
p_ghEqn.setReference(p_ghRefCell, p_ghRefValue);
p_ghEqn.solve(mesh.solver(p_gh.select(pimple.finalInnerIter())));
if (pimple.finalNonOrthogonalIter())
{
phi -= p_ghEqn.flux();
}
}
#include "continuityErrs.H"
// Explicitly relax pressure for momentum corrector
p_gh.relax();
p = p_gh + (g & (mesh.C() + zeta - refLevel));
U -= rAU*fvc::grad(p_gh);
U.correctBoundaryConditions();
/*---------------------------------------------------------------------------*\
========= |
\\ / 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/>.
Application
potentialFreeSurfaceFoam
Description
Incompressible Navier-Stokes solver with inclusion of a wave height field
to enable single-phase free-surface approximations
Wave height field, zeta, used by pressure boundary conditions
Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "singlePhaseTransportModel.H"
#include "turbulenceModel.H"
#include "pimpleControl.H"
#include "IObasicSourceList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
pimpleControl pimple(mesh);
#include "createFields.H"
#include "initContinuityErrs.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.run())
{
#include "readTimeControls.H"
#include "CourantNo.H"
#include "setDeltaT.H"
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())
{
#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;
}
// ************************************************************************* //
Supports Markdown
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