buoyantPimpleFoam: Added moving and changing mesh capability

applications/solvers/heatTransfer/buoyantPimpleFoam/Make/options

@@ -1,21 +1,27 @@
 EXE_INC = \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I. \
+    -I../../compressible/rhoPimpleFoam \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/radiationModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
-    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/cfdTools \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude
 
 EXE_LIBS = \
-    -lfiniteVolume \
-    -lsampling \
-    -lmeshTools \
-    -lfvOptions \
     -lcompressibleTransportModels \
     -lfluidThermophysicalModels \
     -lradiationModels \
     -lspecie \
     -lturbulenceModels \
-    -lcompressibleTurbulenceModels
+    -lcompressibleTurbulenceModels \
+    -lfiniteVolume \
+    -ldynamicFvMesh \
+    -ltopoChangerFvMesh \
+    -lmeshTools \
+    -lsampling \
+    -lfvOptions

applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -34,11 +34,16 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "rhoThermo.H"
+#include "dynamicFvMesh.H"
+#include "fluidThermo.H"
 #include "turbulentFluidThermoModel.H"
 #include "radiationModel.H"
-#include "fvOptions.H"
 #include "pimpleControl.H"
+#include "pressureControl.H"
+#include "CorrectPhi.H"
+#include "fvOptions.H"
+#include "localEulerDdtScheme.H"
+#include "fvcSmooth.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -48,36 +53,102 @@ int main(int argc, char *argv[])
 
     #include "setRootCaseLists.H"
     #include "createTime.H"
-    #include "createMesh.H"
-    #include "createControl.H"
+    #include "createDynamicFvMesh.H"
+    #include "createDyMControls.H"
+    #include "initContinuityErrs.H"
     #include "createFields.H"
     #include "createFieldRefs.H"
-    #include "initContinuityErrs.H"
-    #include "createTimeControls.H"
-    #include "compressibleCourantNo.H"
-    #include "setInitialDeltaT.H"
+    #include "createRhoUfIfPresent.H"
 
     turbulence->validate();
 
+    if (!LTS)
+    {
+        #include "compressibleCourantNo.H"
+        #include "setInitialDeltaT.H"
+    }
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
     while (pimple.run(runTime))
     {
-        #include "readTimeControls.H"
-        #include "compressibleCourantNo.H"
-        #include "setDeltaT.H"
+        #include "readDyMControls.H"
+
+        // Store divrhoU from the previous mesh so that it can be mapped
+        // and used in correctPhi to ensure the corrected phi has the
+        // same divergence
+        autoPtr<volScalarField> divrhoU;
+        if (correctPhi)
+        {
+            divrhoU = new volScalarField
+            (
+                "divrhoU",
+                fvc::div(fvc::absolute(phi, rho, U))
+            );
+        }
+
+        if (LTS)
+        {
+            #include "setRDeltaT.H"
+        }
+        else
+        {
+            #include "compressibleCourantNo.H"
+            #include "setDeltaT.H"
+        }
 
         runTime++;
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "rhoEqn.H"
-
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {
+            if (pimple.firstPimpleIter() || moveMeshOuterCorrectors)
+            {
+                // Store momentum to set rhoUf for introduced faces.
+                autoPtr<volVectorField> rhoU;
+                if (rhoUf.valid())
+                {
+                    rhoU = new volVectorField("rhoU", rho*U);
+                }
+
+                // Do any mesh changes
+                mesh.update();
+
+                if (mesh.changing())
+                {
+                    gh = (g & mesh.C()) - ghRef;
+                    ghf = (g & mesh.Cf()) - ghRef;
+
+                    MRF.update();
+
+                    if (correctPhi)
+                    {
+                        // Calculate absolute flux
+                        // from the mapped surface velocity
+                        phi = mesh.Sf() & rhoUf();
+
+                        #include "correctPhi.H"
+
+                        // Make the fluxes relative to the mesh-motion
+                        fvc::makeRelative(phi, rho, U);
+                    }
+
+                    if (checkMeshCourantNo)
+                    {
+                        #include "meshCourantNo.H"
+                    }
+                }
+            }
+
+            if (pimple.firstPimpleIter() && !pimple.simpleRho())
+            {
+                #include "rhoEqn.H"
+            }
+
             #include "UEqn.H"
             #include "EEqn.H"
 

applications/solvers/heatTransfer/buoyantPimpleFoam/correctPhi.H

@@ -0,0 +1,12 @@
+CorrectPhi
+(
+    U,
+    phi,
+    p_rgh,
+    rho,
+    psi,
+    dimensionedScalar("rAUf", dimTime, 1),
+    divrhoU(),
+    pimple,
+    true
+);

applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H

@@ -1,9 +1,16 @@
+#include "createRDeltaT.H"
+
 Info<< "Reading thermophysical properties\n" << endl;
 
-autoPtr<rhoThermo> pThermo(rhoThermo::New(mesh));
-rhoThermo& thermo = pThermo();
+autoPtr<fluidThermo> pThermo
+(
+    fluidThermo::New(mesh)
+);
+fluidThermo& thermo = pThermo();
 thermo.validate(args.executable(), "h", "e");
 
+volScalarField& p = thermo.p();
+
 volScalarField rho
 (
     IOobject
@@ -17,8 +24,6 @@ volScalarField rho
     thermo.rho()
 );
 
-volScalarField& p = thermo.p();
-
 Info<< "Reading field U\n" << endl;
 volVectorField U
 (
@@ -35,7 +40,6 @@ volVectorField U
 
 #include "compressibleCreatePhi.H"
 
-
 Info<< "Creating turbulence model\n" << endl;
 autoPtr<compressible::turbulenceModel> turbulence
 (
@@ -71,20 +75,7 @@ volScalarField p_rgh
 // Force p_rgh to be consistent with p
 p_rgh = p - rho*gh;
 
-label pRefCell = 0;
-scalar pRefValue = 0.0;
-
-if (thermo.incompressible())
-{
-    setRefCell
-    (
-        p,
-        p_rgh,
-        pimple.dict(),
-        pRefCell,
-        pRefValue
-    );
-}
+pressureControl pressureControl(p, rho, pimple.dict(), false);
 
 mesh.setFluxRequired(p_rgh.name());
 

applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H

@@ -1,4 +1,7 @@
-rho = thermo.rho();
+if (!pimple.simpleRho())
+{
+    rho = thermo.rho();
+}
 
 // Thermodynamic density needs to be updated by psi*d(p) after the
 // pressure solution
@@ -8,8 +11,6 @@ volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
 
-surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
-
 surfaceScalarField phiHbyA
 (
     "phiHbyA",
@@ -17,14 +18,19 @@ surfaceScalarField phiHbyA
         fvc::interpolate(rho)*fvc::flux(HbyA)
       + MRF.zeroFilter(rhorAUf*fvc::ddtCorr(rho, U, phi))
     )
-  + phig
 );
 
 MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
+surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
+
+phiHbyA += phig;
+
 // Update the pressure BCs to ensure flux consistency
 constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
 
+fvc::makeRelative(phiHbyA, rho, U);
+
 fvScalarMatrix p_rghDDtEqn
 (
     fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
@@ -35,16 +41,13 @@ fvScalarMatrix p_rghDDtEqn
 
 while (pimple.correctNonOrthogonal())
 {
-    fvScalarMatrix p_rghEqn
-    (
-        p_rghDDtEqn
-      - fvm::laplacian(rhorAUf, p_rgh)
-    );
+    fvScalarMatrix p_rghEqn(p_rghDDtEqn - fvm::laplacian(rhorAUf, p_rgh));
 
-    if (thermo.incompressible())
-    {
-        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
-    }
+    p_rghEqn.setReference
+    (
+        pressureControl.refCell(),
+        pressureControl.refValue()
+    );
 
     p_rghEqn.solve();
 
@@ -67,13 +70,30 @@ while (pimple.correctNonOrthogonal())
 
 p = p_rgh + rho*gh;
 
+bool limitedp = pressureControl.limit(p);
+
+if (limitedp)
+{
+    p_rgh = p - rho*gh;
+}
+
 // Thermodynamic density update
 thermo.correctRho(psi*p - psip0);
 
+if (limitedp)
+{
+    rho = thermo.rho();
+}
+
+#include "rhoEqn.H"
+#include "compressibleContinuityErrs.H"
+
+if (pimple.simpleRho())
+{
+    rho = thermo.rho();
+}
+
 if (thermo.dpdt())
 {
     dpdt = fvc::ddt(p);
 }
-
-#include "rhoEqn.H"
-#include "compressibleContinuityErrs.H"