incompressible/pimpleFoam, solver + tutorials

applications/solvers/incompressible/pimpleFoam/UEqn.H

@@ -1,43 +1,40 @@
-// Solve the Momentum equation
-
-tmp<fvVectorMatrix> UEqn
+// Convection-diffusion matrx
+fvVectorMatrix HUEqn
 (
-    fvm::ddt(U)
-  + fvm::div(phi, U)
+    fvm::div(phi, U)
   + turbulence->divDevReff(U)
 );
 
+// Time derivative matrix
+fvVectorMatrix ddtUEqn(fvm::ddt(U));
+
+// Get under-relaxation factor
+scalar UUrf = mesh.solutionDict().relaxationFactor(U.name());
+
 if (oCorr == nOuterCorr - 1)
 {
     if (mesh.solutionDict().relax("UFinal"))
     {
-        UEqn().relax(mesh.solutionDict().relaxationFactor("UFinal"));
+        UUrf = mesh.solutionDict().relaxationFactor("UFinal");
     }
     else
     {
-        UEqn().relax(1);
+        UUrf = 1;
     }
 }
-else
-{
-    UEqn().relax();
-}
-
-volScalarField rUA = 1.0/UEqn().A();
 
 if (momentumPredictor)
 {
-    if (oCorr == nOuterCorr-1)
-    {
-        solve(UEqn() == -fvc::grad(p), mesh.solutionDict().solver("UFinal"));
-    }
-    else
-    {
-        solve(UEqn() == -fvc::grad(p));
-    }
+    solve
+    (
+        ddtUEqn
+      + relax(HUEqn, UUrf)
+     ==
+      - fvc::grad(p)
+    );
 }
 else
 {
-    U = rUA*(UEqn().H() - fvc::grad(p));
+    U = (ddtUEqn.H() + HUEqn.H() - fvc::grad(p))/(HUEqn.A() + ddtUEqn.A());
     U.correctBoundaryConditions();
 }

applications/solvers/incompressible/pimpleFoam/createFields.H

@@ -40,3 +40,19 @@ 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/incompressible/pimpleFoam/pEqn.H

@@ -1,53 +1,67 @@
-U = rUA*UEqn().H();
-
-if (nCorr <= 1)
 {
-    UEqn.clear();
-}
+    p.boundaryField().updateCoeffs();
 
-phi = (fvc::interpolate(U) & mesh.Sf())
-    + fvc::ddtPhiCorr(rUA, U, phi);
+    // Prepare clean Ap without time derivative contribution and
+    // without contribution from under-relaxation
+    // HJ, 26/Oct/2015
+    aU = HUEqn.A();
 
-adjustPhi(phi, U, p);
+    // Store velocity under-relaxation point before using U for the flux
+    // precursor
+    U.storePrevIter();
 
-// Non-orthogonal pressure corrector loop
-for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
-{
-    // Pressure corrector
-    fvScalarMatrix pEqn
-    (
-        fvm::laplacian(rUA, p) == fvc::div(phi)
-    );
+    U = HUEqn.H()/aU;
+    phi = (fvc::interpolate(U) & mesh.Sf());
 
-    pEqn.setReference(pRefCell, pRefValue);
+    // Global flux balance
+    adjustPhi(phi, U, p);
 
-    if
-    (
-        oCorr == nOuterCorr - 1
-     && corr == nCorr - 1
-     && nonOrth == nNonOrthCorr
-    )
+    for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++)
     {
-        pEqn.solve(mesh.solutionDict().solver("pFinal"));
-    }
-    else
-    {
-        pEqn.solve();
+        fvScalarMatrix pEqn
+        (
+            fvm::laplacian(1/aU, p) == fvc::div(phi)
+        );
+
+        pEqn.setReference(pRefCell, pRefValue);
+
+        if
+        (
+            corr == nCorr - 1
+         && nonOrth == nNonOrthCorr
+        )
+        {
+            pEqn.solve
+            (
+                mesh.solutionDict().solver(p.name() + "Final")
+            );
+        }
+        else
+        {
+            pEqn.solve(mesh.solutionDict().solver(p.name()));
+        }
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            phi -= pEqn.flux();
+        }
     }
 
-    if (nonOrth == nNonOrthCorr)
+    // Explicitly relax pressure for momentum corrector
+    if (oCorr != nOuterCorr - 1)
     {
-        phi -= pEqn.flux();
+        p.relax();
     }
+
+#   include "movingMeshContinuityErrs.H"
+
+    U = UUrf*
+        (
+            1.0/(aU + ddtUEqn.A())*
+            (
+                U*aU - fvc::grad(p) + ddtUEqn.H()
+            )
+        )
+      + (1 - UUrf)*U.prevIter();
+    U.correctBoundaryConditions();
 }
-
-#include "continuityErrs.H"
-
-// Explicitly relax pressure for momentum corrector except for last corrector
-if (oCorr != nOuterCorr-1)
-{
-    p.relax();
-}
-
-U -= rUA*fvc::grad(p);
-U.correctBoundaryConditions();

applications/solvers/incompressible/pimpleFoam/pimpleFoam.C

@@ -30,6 +30,11 @@ Description
 
     Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
 
+    Consistent formulation without time-step and relaxation dependence by Jasak
+
+Author
+    Hrvoje Jasak, Wikki Ltd.  All rights reserved
+
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"

tutorials/incompressible/pimpleFoam/t-junction/system/fvSchemes

@@ -43,7 +43,7 @@ laplacianSchemes
 {
     default         none;
     laplacian(nuEff,U) Gauss linear corrected;
-    laplacian((1|A(U)),p) Gauss linear corrected;
+    laplacian((1|aU),p) Gauss linear corrected;
     laplacian(DkEff,k) Gauss linear corrected;
     laplacian(DepsilonEff,epsilon) Gauss linear corrected;
     laplacian(DREff,R) Gauss linear corrected;