Block-coupled p-U solver improvements

applications/solvers/coupled/MRFPorousFoam/MRFPorousFoam.C

@@ -90,7 +90,7 @@ int main(int argc, char *argv[])
         U.correctBoundaryConditions();
         p.correctBoundaryConditions();
 
-        phi = (fvc::interpolate(U) & mesh.Sf()) + pEqn.flux() + presSource;
+        phi = (fvc::interpolate(U) & mesh.Sf()) + tpEqn().flux() + tpresSource;
 
         // Make flux relative in rotating zones
         mrfZones.relativeFlux(phi);

applications/solvers/coupled/MRFPorousFoam/UEqn.H

@@ -1,15 +1,94 @@
-    // Momentum equation
-    tmp<fvVectorMatrix> UEqn
+{
+    volScalarField divPhi
     (
-        fvm::ddt(U)
-      + fvm::div(phi, U)
+        "divPhi",
+        fvc::div(phi)
+    );
+
+    // Momentum equation
+    fvVectorMatrix UEqn
+    (
+        fvm::div(phi, U)
       + turbulence->divDevReff(U)
     );
 
-    // Add MRF and porous sources
-    mrfZones.addCoriolis(UEqn());
-    pZones.addResistance(UEqn());
+    // Add MRF sources
+    mrfZones.addCoriolis(UEqn);
 
-    UEqn().relax();
+    // Add porous sources
+    tmp<volTensorField> tTU;
 
-    UpEqn.insertEquation(0, UEqn());
+    if (addPorosity)
+    {
+        tTU = tmp<volTensorField>
+        (
+            new volTensorField
+            (
+                IOobject
+                (
+                    "TU",
+                    runTime.timeName(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                mesh,
+                dimensionedTensor("zero", dimless/dimTime, tensor::zero)
+            )
+        );
+        volTensorField& TU = tTU();
+
+        pZones.addResistance(UEqn, TU);
+
+        trTU = inv(TU + tensor(I)*UEqn.A());
+        trTU().rename("rAU");
+    }
+    else
+    {
+        trAU = 1.0/UEqn.A();
+        trAU().rename("rAU");
+    }
+
+    // Insert the additional components. Note this will destroy the H and A
+
+    UEqn += fvm::SuSp(-divPhi, U) + divPhi*U;
+    UEqn.relax();
+
+    // Insert momentum equation
+    UpEqn.insertEquation(0, UEqn);
+
+    if (addPorosity)
+    {
+        // Manually over-ride the 3x3 block to handle the off-diagonal
+        // part of the Ap coefficient
+        const tensorField& TUIn = tTU().internalField();
+
+        CoeffField<vector4>::squareTypeField& DD = UpEqn.diag().asSquare();
+
+        const scalarField& V = mesh.V().field();
+
+        // Note: this insertion should only happen in porous cell zones
+        // A rewrite of the porousZones class interface is needed.
+        // HJ, 14/Mar/2016
+        forAll (TUIn, cellI)
+        {
+            const scalar& cellV = V[cellI];
+
+            const tensor& cellTU = TUIn[cellI];
+
+            CoeffField<vector4>::squareType& cellDD = DD[cellI];
+
+            cellDD(0, 0) += cellV*cellTU.xx();
+            cellDD(0, 1) += cellV*cellTU.xy();
+            cellDD(0, 2) += cellV*cellTU.xz();
+
+            cellDD(1, 0) += cellV*cellTU.yx();
+            cellDD(1, 1) += cellV*cellTU.yy();
+            cellDD(2, 2) += cellV*cellTU.yz();
+
+            cellDD(2, 0) += cellV*cellTU.zx();
+            cellDD(2, 1) += cellV*cellTU.zy();
+            cellDD(2, 2) += cellV*cellTU.zz();
+        }
+    }
+}

applications/solvers/coupled/MRFPorousFoam/createFields.H

@@ -50,3 +50,9 @@ volVector4Field Up
     mesh,
     dimensionedVector4("zero", dimless, vector4::zero)
 );
+
+Info<< "Creating field rAU\n" << endl;
+// Note: depending on the porosity treatment, only one of the rAU fields
+// will be used.  HJ, 1/Mar/2016
+tmp<volScalarField> trAU;
+tmp<volTensorField> trTU;

applications/solvers/coupled/MRFPorousFoam/createPorous.H

@@ -1 +1,8 @@
     porousZones pZones(mesh);
+
+    bool addPorosity = false;
+
+    if (!pZones.empty())
+    {
+        addPorosity = true;
+    }

applications/solvers/coupled/MRFPorousFoam/pEqn.H

@@ -1,25 +1,48 @@
 // Pressure parts of the continuity equation
-surfaceScalarField rUAf
-(
-    "rUAf",
-    fvc::interpolate(1.0/UEqn().A())
-);
+tmp<fvScalarMatrix> tpEqn;
+tmp<surfaceScalarField> tpresSource;
 
-UEqn.clear();
+if (addPorosity)
+{
+    // Collect pressure source with tensorial 1/Ap
+    const volTensorField& rTU = trTU();
 
-surfaceScalarField presSource
-(
-    "presSource",
-    rUAf*(fvc::interpolate(fvc::grad(p, "grad(pSource)")) & mesh.Sf())
-);
+    tpresSource =
+    (
+        (mesh.Sf() & fvc::interpolate(rTU))
+        & fvc::interpolate(fvc::grad(p, "grad(pSource)"))
+    );
+    const surfaceScalarField& presSource = tpresSource();
 
-fvScalarMatrix pEqn
-(
-  - fvm::laplacian(rUAf, p)
- ==
-  - fvc::div(presSource)
-);
+    // Assemble pressure matrix with tensorial 1/Ap
+    tpEqn =
+    (
+      - fvm::laplacian(rTU, p)
+     ==
+      - fvc::div(presSource)
+    );
+}
+else
+{
+    // Collect pressure source with scalar 1/Ap
+    const volScalarField& rAU = trAU();
 
-pEqn.setReference(pRefCell, pRefValue);
+    tpresSource =
+    (
+        fvc::interpolate(rAU)*
+        (mesh.Sf() & fvc::interpolate(fvc::grad(p, "grad(pSource)")))
+    );
+    const surfaceScalarField& presSource = tpresSource();
 
-UpEqn.insertEquation(3, pEqn);
+    // Assemble pressure matrix with tensorial 1/Ap
+    tpEqn =
+    (
+        - fvm::laplacian(rAU, p)
+       ==
+        - fvc::div(presSource)
+    );
+}
+
+tpEqn().setReference(pRefCell, pRefValue);
+
+UpEqn.insertEquation(3, tpEqn());

applications/solvers/coupled/pUCoupledFoam/UEqn.H

@@ -1,11 +1,23 @@
-// Momentum equation
-tmp<fvVectorMatrix> UEqn
-(
-    fvm::ddt(U)
-  + fvm::div(phi, U)
-  + turbulence->divDevReff(U)
-);
+{
+    volScalarField divPhi
+    (
+        "divPhi",
+        fvc::div(phi)
+    );
 
-UEqn().relax();
+    // Momentum equation
+    fvVectorMatrix UEqn
+    (
+        fvm::div(phi, U)
+      + turbulence->divDevReff(U)
+    );
 
-UpEqn.insertEquation(0, UEqn());
+    rAU = 1.0/UEqn.A();
+
+    // Insert the additional components. Note this will destroy the H and A
+
+    UEqn += fvm::SuSp(-divPhi, U) + divPhi*U;
+    UEqn.relax();
+
+    UpEqn.insertEquation(0, UEqn);
+}

applications/solvers/coupled/pUCoupledFoam/createFields.H

@@ -50,3 +50,18 @@ volVector4Field Up
     mesh,
     dimensionedVector4("zero", dimless, vector4::zero)
 );
+
+Info<< "Creating field rAU\n" << endl;
+volScalarField rAU
+(
+    IOobject
+    (
+        "rAU",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    runTime.deltaT()
+);

applications/solvers/coupled/pUCoupledFoam/pEqn.H

@@ -1,21 +1,14 @@
 // Pressure parts of the continuity equation
-surfaceScalarField rUAf
-(
-    "rUAf",
-    fvc::interpolate(1.0/UEqn().A())
-);
-
-UEqn.clear();
-
 surfaceScalarField presSource
 (
     "presSource",
-    rUAf*(fvc::interpolate(fvc::grad(p, "grad(pSource)")) & mesh.Sf())
+    fvc::interpolate(rAU)*
+    (fvc::interpolate(fvc::grad(p, "grad(pSource)")) & mesh.Sf())
 );
 
 fvScalarMatrix pEqn
 (
-  - fvm::laplacian(rUAf, p)
+  - fvm::laplacian(rAU, p)
  ==
   - fvc::div(presSource)
 );