From cc774a93bd509546c4ce50544f9e0824d4bff6b3 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Fri, 29 Dec 2017 13:57:35 +0000 Subject: [PATCH] Incompressible dynamic mesh solver with Immersed Boundary support --- .../pimpleDyMIbFoam/Make/files | 3 + .../pimpleDyMIbFoam/Make/options | 23 +++ .../immersedBoundary/pimpleDyMIbFoam/UEqn.H | 26 ++++ .../pimpleDyMIbFoam/correctPhi.H | 31 ++++ .../pimpleDyMIbFoam/createControls.H | 11 ++ .../pimpleDyMIbFoam/createFields.H | 59 ++++++++ .../immersedBoundary/pimpleDyMIbFoam/pEqn.H | 58 ++++++++ .../pimpleDyMIbFoam/pimpleDyMIbFoam.C | 138 ++++++++++++++++++ .../pimpleDyMIbFoam/readControls.H | 5 + 9 files changed, 354 insertions(+) create mode 100644 applications/solvers/immersedBoundary/pimpleDyMIbFoam/Make/files create mode 100644 applications/solvers/immersedBoundary/pimpleDyMIbFoam/Make/options create mode 100644 applications/solvers/immersedBoundary/pimpleDyMIbFoam/UEqn.H create mode 100644 applications/solvers/immersedBoundary/pimpleDyMIbFoam/correctPhi.H create mode 100644 applications/solvers/immersedBoundary/pimpleDyMIbFoam/createControls.H create mode 100644 applications/solvers/immersedBoundary/pimpleDyMIbFoam/createFields.H create mode 100644 applications/solvers/immersedBoundary/pimpleDyMIbFoam/pEqn.H create mode 100644 applications/solvers/immersedBoundary/pimpleDyMIbFoam/pimpleDyMIbFoam.C create mode 100644 applications/solvers/immersedBoundary/pimpleDyMIbFoam/readControls.H diff --git a/applications/solvers/immersedBoundary/pimpleDyMIbFoam/Make/files b/applications/solvers/immersedBoundary/pimpleDyMIbFoam/Make/files new file mode 100644 index 000000000..46c7c1f2a --- /dev/null +++ b/applications/solvers/immersedBoundary/pimpleDyMIbFoam/Make/files @@ -0,0 +1,3 @@ +pimpleDyMIbFoam.C + +EXE = $(FOAM_APPBIN)/pimpleDyMIbFoam diff --git a/applications/solvers/immersedBoundary/pimpleDyMIbFoam/Make/options b/applications/solvers/immersedBoundary/pimpleDyMIbFoam/Make/options new file mode 100644 index 000000000..a011e27ec --- /dev/null +++ b/applications/solvers/immersedBoundary/pimpleDyMIbFoam/Make/options @@ -0,0 +1,23 @@ +EXE_INC = \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/immersedBoundary/immersedBoundary/lnInclude \ + -I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \ + -I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \ + -I$(LIB_SRC)/transportModels \ + -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel + +EXE_LIBS = \ + -lfiniteVolume \ + -limmersedBoundary \ + -ldynamicFvMesh \ + -ltopoChangerFvMesh \ + -ldynamicMesh \ + -lmeshTools \ + -lincompressibleTransportModels \ + -lincompressibleTurbulenceModel \ + -lincompressibleRASModels \ + -lincompressibleLESModels \ + -llduSolvers \ + -L$(MESQUITE_LIB_DIR) -lmesquite diff --git a/applications/solvers/immersedBoundary/pimpleDyMIbFoam/UEqn.H b/applications/solvers/immersedBoundary/pimpleDyMIbFoam/UEqn.H new file mode 100644 index 000000000..44acc8af0 --- /dev/null +++ b/applications/solvers/immersedBoundary/pimpleDyMIbFoam/UEqn.H @@ -0,0 +1,26 @@ + // Convection-diffusion matrix + fvVectorMatrix HUEqn + ( + fvm::div(phi, U) + + turbulence->divDevReff() + ); + + // Time derivative matrix + fvVectorMatrix ddtUEqn(fvm::ddt(U)); + + // Get under-relaxation factor + scalar UUrf = + mesh.solutionDict().equationRelaxationFactor(U.select(pimple.finalIter())); + + if (pimple.momentumPredictor()) + { + // Solve momentum predictor + solve + ( + ddtUEqn + + relax(HUEqn, UUrf) + == + - fvc::grad(p), + mesh.solutionDict().solver((U.select(pimple.finalIter()))) + ); + } diff --git a/applications/solvers/immersedBoundary/pimpleDyMIbFoam/correctPhi.H b/applications/solvers/immersedBoundary/pimpleDyMIbFoam/correctPhi.H new file mode 100644 index 000000000..eda2e9eda --- /dev/null +++ b/applications/solvers/immersedBoundary/pimpleDyMIbFoam/correctPhi.H @@ -0,0 +1,31 @@ +{ + volScalarField pcorr("pcorr", p); + + // Initialise flux with interpolated velocity + phi = fvc::interpolate(U) & mesh.Sf(); + + adjustPhi(phi, U, pcorr); + + mesh.schemesDict().setFluxRequired(pcorr.name()); + + while (pimple.correctNonOrthogonal()) + { + fvScalarMatrix pcorrEqn + ( + fvm::laplacian(1/aU, pcorr) == fvc::div(phi) + ); + + pcorrEqn.setReference(pRefCell, pRefValue); + pcorrEqn.solve(); + + if (pimple.finalNonOrthogonalIter()) + { + phi -= pcorrEqn.flux(); + } + + // Fluxes are corrected to absolute velocity and further corrected + // later. HJ, 6/Feb/2009 + } +# include "continuityErrs.H" +} + diff --git a/applications/solvers/immersedBoundary/pimpleDyMIbFoam/createControls.H b/applications/solvers/immersedBoundary/pimpleDyMIbFoam/createControls.H new file mode 100644 index 000000000..ca5e25906 --- /dev/null +++ b/applications/solvers/immersedBoundary/pimpleDyMIbFoam/createControls.H @@ -0,0 +1,11 @@ +#include "createTimeControls.H" + +bool correctPhi +( + pimple.dict().lookupOrDefault("correctPhi", false) +); + +bool checkMeshCourantNo +( + pimple.dict().lookupOrDefault("checkMeshCourantNo", false) +); diff --git a/applications/solvers/immersedBoundary/pimpleDyMIbFoam/createFields.H b/applications/solvers/immersedBoundary/pimpleDyMIbFoam/createFields.H new file mode 100644 index 000000000..74a461355 --- /dev/null +++ b/applications/solvers/immersedBoundary/pimpleDyMIbFoam/createFields.H @@ -0,0 +1,59 @@ + Info<< "Reading field p\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" + + + label pRefCell = 0; + scalar pRefValue = 0.0; + setRefCell(p, pimple.dict(), pRefCell, pRefValue); + mesh.schemesDict().setFluxRequired(p.name()); + + singlePhaseTransportModel laminarTransport(U, phi); + + autoPtr turbulence + ( + incompressible::turbulenceModel::New(U, phi, laminarTransport) + ); + + Info<< "Reading field aU if present\n" << endl; + volScalarField aU + ( + IOobject + ( + "aU", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + mesh, + 1/runTime.deltaT(), + zeroGradientFvPatchScalarField::typeName + ); diff --git a/applications/solvers/immersedBoundary/pimpleDyMIbFoam/pEqn.H b/applications/solvers/immersedBoundary/pimpleDyMIbFoam/pEqn.H new file mode 100644 index 000000000..1c02221ba --- /dev/null +++ b/applications/solvers/immersedBoundary/pimpleDyMIbFoam/pEqn.H @@ -0,0 +1,58 @@ +{ + p.boundaryField().updateCoeffs(); + + // Prepare clean Ap without time derivative contribution and + // without contribution from under-relaxation + // HJ, 26/Oct/2015 + aU = HUEqn.A(); + + // Store velocity under-relaxation point before using U for the flux + // precursor + U.storePrevIter(); + + U = HUEqn.H()/aU; + phi = (fvc::interpolate(U) & mesh.Sf()); + + // Global flux balance + adjustPhi(phi, U, p); + + while (pimple.correctNonOrthogonal()) + { + fvScalarMatrix pEqn + ( + fvm::laplacian(1/aU, p) == fvc::div(phi) + ); + + pEqn.setReference(pRefCell, pRefValue); + pEqn.solve + ( + mesh.solutionDict().solver(p.select(pimple.finalInnerIter())) + ); + + if (pimple.finalNonOrthogonalIter()) + { + phi -= pEqn.flux(); + } + } + + // Explicitly relax pressure for momentum corrector except for last corrector + if (!pimple.finalIter()) + { + p.relax(); + } + + // Make the fluxes relative to the mesh motion + fvc::makeRelative(phi, U); + +# include "movingMeshContinuityErrs.H" + + U = UUrf* + ( + 1.0/(aU + ddtUEqn.A())* + ( + U*aU - fvc::grad(p) + ddtUEqn.H() + ) + ) + + (1 - UUrf)*U.prevIter(); + U.correctBoundaryConditions(); +} diff --git a/applications/solvers/immersedBoundary/pimpleDyMIbFoam/pimpleDyMIbFoam.C b/applications/solvers/immersedBoundary/pimpleDyMIbFoam/pimpleDyMIbFoam.C new file mode 100644 index 000000000..f50b32a1c --- /dev/null +++ b/applications/solvers/immersedBoundary/pimpleDyMIbFoam/pimpleDyMIbFoam.C @@ -0,0 +1,138 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 4.0 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend 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. + + foam-extend 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 foam-extend. If not, see . + +Application + pimpleDyMFoam.C + +Description + Transient solver for incompressible, flow of Newtonian fluids + with dynamic mesh using the PIMPLE (merged PISO-SIMPLE) algorithm. + + Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected. + + Consistent formulation without time-step and relaxation dependence by Jasak + + Support for immersed boundary + +Author + Hrvoje Jasak, Wikki Ltd. All rights reserved + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "singlePhaseTransportModel.H" +#include "turbulenceModel.H" +#include "dynamicFvMesh.H" +#include "pimpleControl.H" + +#include "immersedBoundaryPolyPatch.H" +#include "immersedBoundaryFvPatch.H" +#include "emptyFvPatch.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ +# include "setRootCase.H" + +# include "createTime.H" +# include "createDynamicFvMesh.H" + + pimpleControl pimple(mesh); + +# include "initContinuityErrs.H" +# include "createIbMasks.H" +# include "createFields.H" +# include "createControls.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< "\nStarting time loop\n" << endl; + + while (runTime.run()) + { +# include "readControls.H" +# include "CourantNo.H" +# include "setDeltaT.H" + + runTime++; + + Info<< "Time = " << runTime.timeName() << nl << endl; + + bool meshChanged = mesh.update(); + + U.correctBoundaryConditions(); + p.correctBoundaryConditions(); + aU.correctBoundaryConditions(); + phi.correctBoundaryConditions(); + turbulence->correct(); + +# include "updateIbMasks.H" +# include "volContinuity.H" + + if (runTime.outputTime()) + { + volScalarField divMeshPhi("divMeshPhi", mag(fvc::surfaceIntegrate(mesh.phi()))); + divMeshPhi.write(); + } + + if (checkMeshCourantNo) + { +# include "meshCourantNo.H" + } + + // Fluxes will be corrected to absolute velocity + // HJ, 6/Feb/2009 +# include "correctPhi.H" + + // Make the fluxes relative to the mesh motion + fvc::makeRelative(phi, U); + + // --- PIMPLE loop + while (pimple.loop()) + { +# include "UEqn.H" + + // --- PISO loop + while (pimple.correct()) + { +# include "pEqn.H" + } + + turbulence->correct(); + } + + runTime.write(); + + Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" + << " ClockTime = " << runTime.elapsedClockTime() << " s" + << nl << endl; + } + + Info<< "End\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/applications/solvers/immersedBoundary/pimpleDyMIbFoam/readControls.H b/applications/solvers/immersedBoundary/pimpleDyMIbFoam/readControls.H new file mode 100644 index 000000000..9f982e260 --- /dev/null +++ b/applications/solvers/immersedBoundary/pimpleDyMIbFoam/readControls.H @@ -0,0 +1,5 @@ +#include "readTimeControls.H" + +correctPhi = pimple.dict().lookupOrDefault("correctPhi", false); + +checkMeshCourantNo = pimple.dict().lookupOrDefault("checkMeshCourantNo", false);