Incompressible dynamic mesh solver with Immersed Boundary support

applications/solvers/immersedBoundary/pimpleDyMIbFoam/Make/files

@@ -0,0 +1,3 @@
+pimpleDyMIbFoam.C
+
+EXE = $(FOAM_APPBIN)/pimpleDyMIbFoam

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

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())))
+        );
+    }

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"
+}
+

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)
+);

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<incompressible::turbulenceModel> 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
+    );

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();
+}

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 <http://www.gnu.org/licenses/>.
+
+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;
+}
+
+
+// ************************************************************************* //

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);