From 3aadb63c4b4cd5867f94e933aec2cf472592377a Mon Sep 17 00:00:00 2001 From: andy Date: Fri, 18 Nov 2011 10:07:46 +0000 Subject: [PATCH] ENH: Added new potentialFreeSurfaceFoam solver - single phase with wave height inclusion --- .../potentialFreeSurfaceFoam/Make/files | 3 + .../potentialFreeSurfaceFoam/Make/options | 14 ++ .../potentialFreeSurfaceFoam/UEqn.H | 18 +++ .../potentialFreeSurfaceFoam/createFields.H | 121 ++++++++++++++++++ .../potentialFreeSurfaceFoam/pEqn.H | 43 +++++++ .../potentialFreeSurfaceFoam.C | 101 +++++++++++++++ 6 files changed, 300 insertions(+) create mode 100644 applications/solvers/incompressible/potentialFreeSurfaceFoam/Make/files create mode 100644 applications/solvers/incompressible/potentialFreeSurfaceFoam/Make/options create mode 100644 applications/solvers/incompressible/potentialFreeSurfaceFoam/UEqn.H create mode 100644 applications/solvers/incompressible/potentialFreeSurfaceFoam/createFields.H create mode 100644 applications/solvers/incompressible/potentialFreeSurfaceFoam/pEqn.H create mode 100644 applications/solvers/incompressible/potentialFreeSurfaceFoam/potentialFreeSurfaceFoam.C diff --git a/applications/solvers/incompressible/potentialFreeSurfaceFoam/Make/files b/applications/solvers/incompressible/potentialFreeSurfaceFoam/Make/files new file mode 100644 index 0000000000..a027148d26 --- /dev/null +++ b/applications/solvers/incompressible/potentialFreeSurfaceFoam/Make/files @@ -0,0 +1,3 @@ +potentialFreeSurfaceFoam.C + +EXE = $(FOAM_APPBIN)/potentialFreeSurfaceFoam diff --git a/applications/solvers/incompressible/potentialFreeSurfaceFoam/Make/options b/applications/solvers/incompressible/potentialFreeSurfaceFoam/Make/options new file mode 100644 index 0000000000..ecb86aecbc --- /dev/null +++ b/applications/solvers/incompressible/potentialFreeSurfaceFoam/Make/options @@ -0,0 +1,14 @@ +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 diff --git a/applications/solvers/incompressible/potentialFreeSurfaceFoam/UEqn.H b/applications/solvers/incompressible/potentialFreeSurfaceFoam/UEqn.H new file mode 100644 index 0000000000..c869df39a2 --- /dev/null +++ b/applications/solvers/incompressible/potentialFreeSurfaceFoam/UEqn.H @@ -0,0 +1,18 @@ +tmp 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)); +} diff --git a/applications/solvers/incompressible/potentialFreeSurfaceFoam/createFields.H b/applications/solvers/incompressible/potentialFreeSurfaceFoam/createFields.H new file mode 100644 index 0000000000..4437f4983d --- /dev/null +++ b/applications/solvers/incompressible/potentialFreeSurfaceFoam/createFields.H @@ -0,0 +1,121 @@ + 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 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); diff --git a/applications/solvers/incompressible/potentialFreeSurfaceFoam/pEqn.H b/applications/solvers/incompressible/potentialFreeSurfaceFoam/pEqn.H new file mode 100644 index 0000000000..a5adeb9e91 --- /dev/null +++ b/applications/solvers/incompressible/potentialFreeSurfaceFoam/pEqn.H @@ -0,0 +1,43 @@ +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(); diff --git a/applications/solvers/incompressible/potentialFreeSurfaceFoam/potentialFreeSurfaceFoam.C b/applications/solvers/incompressible/potentialFreeSurfaceFoam/potentialFreeSurfaceFoam.C new file mode 100644 index 0000000000..35a787e97b --- /dev/null +++ b/applications/solvers/incompressible/potentialFreeSurfaceFoam/potentialFreeSurfaceFoam.C @@ -0,0 +1,101 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 . + +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; +} + + +// ************************************************************************* //