Added solver for pseudo-steady conditions.

applications/solvers/cfdemSolverRhoSimple/Make/files

@@ -1,3 +0,0 @@
-cfdemSolverRhoSimple.C
-
-EXE=$(CFDEM_APP_DIR)/cfdemSolverRhoSimple

applications/solvers/cfdemSolverRhoSteadyPimple/Make/files

@@ -0,0 +1,3 @@
+cfdemSolverRhoSteadyPimple.C
+
+EXE=$(CFDEM_APP_DIR)/cfdemSolverRhoSteadyPimple

applications/solvers/cfdemSolverRhoSteadyPimple/cfdemSolverRhoSteadyPimple.C

@@ -31,11 +31,12 @@ Description
 #include "psiThermo.H"
 #include "turbulentFluidThermoModel.H"
 #include "bound.H"
-#include "simpleControl.H"
+#include "pimpleControl.H"
 #include "fvOptions.H"
 #include "localEulerDdtScheme.H"
 #include "fvcSmooth.H"
 
+
 #include "cfdemCloudEnergy.H"
 #include "implicitCouple.H"
 #include "clockModel.H"
@@ -68,13 +69,21 @@ int main(int argc, char *argv[])
     #include "checkModelType.H"
 
     turbulence->validate();
+  //        #include "compressibleCourantNo.H"
+  //  #include "setInitialDeltaT.H"
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
-    while (simple.loop())
+    while (runTime.run())
     {
+        #include "readTimeControls.H"
+        #include "compressibleCourantNo.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+
         particleCloud.clockM().start(1,"Global");
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
@@ -103,23 +112,38 @@ int main(int argc, char *argv[])
         particleCloud.clockM().stop("Coupling");
 
         particleCloud.clockM().start(26,"Flow");
-
         volScalarField rhoeps("rhoeps",rho*voidfraction);
-        // Pressure-velocity SIMPLE corrector
+        while (pimple.loop())
+        {
+            // if needed, perform drag update here
+            if (pimple.nCorrPIMPLE() <= 1)
+            {
+                #include "rhoEqn.H"
+            }
+
+
+            // --- Pressure-velocity PIMPLE corrector loop
 
             #include "UEqn.H"
-
-
-        // besides this pEqn, OF offers a "simple consistent"-option
-        #include "pEqn.H"
-        rhoeps=rho*voidfraction;
-
             #include "EEqn.H"
 
-        turbulence->correct();
+            // --- Pressure corrector loop
+            while (pimple.correct())
+            {
+                // besides this pEqn, OF offers a "pimple consistent"-option
+                #include "pEqn.H"
+                rhoeps=rho*voidfraction;
+            }
 
-        particleCloud.clockM().start(32,"postFlow");
+            if (pimple.turbCorr())
-        if(hasEvolved) particleCloud.postFlow();
+            {
+                turbulence->correct();
+            }
+        }
+        particleCloud.clockM().stop("Flow");
+
+        particleCloud.clockM().start(31,"postFlow");
+        particleCloud.postFlow();
         particleCloud.clockM().stop("postFlow");
 
         runTime.write();
@@ -129,7 +153,7 @@ int main(int argc, char *argv[])
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
             << nl << endl;
 
-        particleCloud.clockM().stop("Flow");
+
         particleCloud.clockM().stop("Global");
     }
 

applications/solvers/cfdemSolverRhoSteadyPimple/createFields.H

@@ -158,7 +158,7 @@ Info<< "Reading thermophysical properties\n" << endl;
         dimensionedScalar::lookupOrDefault
         (
             "rhoMax",
-            simple.dict(),
+            pimple.dict(),
             dimDensity,
             GREAT
         )
@@ -169,7 +169,7 @@ Info<< "Reading thermophysical properties\n" << endl;
         dimensionedScalar::lookupOrDefault
         (
             "rhoMin",
-            simple.dict(),
+            pimple.dict(),
             dimDensity,
             0
         )
@@ -189,7 +189,7 @@ Info<< "Reading thermophysical properties\n" << endl;
 
     label pRefCell = 0;
     scalar pRefValue = 0.0;
-    setRefCell(p, simple.dict(), pRefCell, pRefValue);
+    setRefCell(p, pimple.dict(), pRefCell, pRefValue);
 
     mesh.setFluxRequired(p.name());
 

applications/solvers/cfdemSolverRhoSteadyPimple/pEqn.H

@@ -14,7 +14,7 @@ volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
 surfaceScalarField phiUs("phiUs", fvc::interpolate(rhoeps*rAU*Ksl*Us)& mesh.Sf());
 
 
-if (simple.transonic())
+if (pimple.transonic())
 {
 //  transonic version not implemented yet
 }
@@ -34,7 +34,7 @@ else
     // Update the pressure BCs to ensure flux consistency
     constrainPressure(p, rhoeps, U, phi, rhorAUf);
 
-    while (simple.correctNonOrthogonal())
+    while (pimple.correctNonOrthogonal())
     {
         // Pressure corrector
         fvScalarMatrix pEqn
@@ -49,7 +49,7 @@ else
 
         pEqn.solve();
 
-        if (simple.finalNonOrthogonalIter())
+        if (pimple.finalNonOrthogonalIter())
         {
             phi += pEqn.flux();
         }

etc/solver-list.txt

@@ -7,6 +7,6 @@ cfdemSolverRhoPimple/dir
 cfdemSolverPisoMS/dir
 cfdemSolverPiso/dir
 cfdemSolverRhoPimple/dir
-cfdemSolverRhoSimple/dir
+cfdemSolverRhoSteadyPimple/dir
 cfdemSolverIB/dir
 cfdemSolverPisoScalar/dir