Update of interIsoFoam and isoAdvection extending to work with porous media

Extension of isoAdvection and interIsoFoam to deal with a porous medium occupying part of the cell volumes and described by a porosity volScalarField taking values larger than 0 (cell full of solid) and less than or equal to one (cell full of fluid).

The current implementation deals with both the advection part of the problem and the modification of the UEqn (p_rghEqn needs no changes).

The aim of the implementation is to make it invisible to the user who does not work with porosity both in terms of efficiency, memory usage and code complexity. It would be preferable if the required interIsoFaom changes could be done without touching the solver files, e.g. implementing the UEqn modifications as an fvOption. This seems to be impossible because we need to modify the existing terms in the UEqn and the fvMatrix is build from right to left meaning the fvOption cannot touch e.g. the fvm::div(rhoPhi,U)-term. The chosen solution is therefore to add an #include "UEqnAddPorosity.H" line to the UEqn.H file AFTER the UEqn is constructed.

The activation of porosity treatment is controlled by a bool parameter called porosityEnabled in a constant/porosityProperties dictionary. This is disabled by default. If it is enabled the user must provide a porosity volScalarField in the 0 time dir. The isoAdvection class has two new data members: porosityPtr_ to the porosity field, is nullptr by default, and a bool called porosityEnabled_ (default false) read from the constant/porosityProperties by the isoAdvection constructor.

It has been necessary to introduce "if (porosityEnabled_)" lines various places in the code but most places this is not inside loops over cells or faces so it is cheap. An exception is in the isoAdvection bounding functions but these are only called for very few cells so this should not be very expensive.

At this point I am not sure that the interplay with the source terms Su and Sp. These should probably also be divided by porosity when alpha is updated.

The interIsoFoam solver now reads the porosityEnabled bool in constant/ porosityProperties,reads the porosity field if needed and the porosity parameters for the Darcy-Forchheimer force and added mass force from the constant/ porosityProperties dict.

I have introduced a porousAlphaCourantNo.H and porousCourantNo.H file as replacements for alphaCourantNo.H and CourantNo.H to correct the adaptive time stepping so it accounts for the increased flow velocity in porous cells if porosity is enabled.

Also the "Phase-1 volume fraction" write out has been modified to take porosity properly into account.

There are two test cases included:

1. discInConstantPorousFlow which is the same as discInConstantFlow but with a porous zone in the middle of the domain where the disc elongates and moves faster.
2. A porousDamBreak testing porosity implementation and comparing with exp. data.

This code contribution is joined work between:

• Konstantinos Missios
• Niels Gjøl Jacobsen
• Henning Scheufler
• Johan Roenby

Johan Roenby, December 2021

applications/solvers/multiphase/interIsoFoam/UEqn.H

+    MRF.correctBoundaryVelocity(U);
+
+    fvVectorMatrix UEqn
+    (
+        fvm::ddt(rho, U) + fvm::div(rhoPhi, U)
+      + MRF.DDt(rho, U)
+      + turbulence->divDevRhoReff(rho, U)
+     ==
+        fvOptions(rho, U)
+    );
+
+    #include "UEqnAddPorosity.H"
+
+    UEqn.relax();
+
+    fvOptions.constrain(UEqn);
+
+    if (pimple.momentumPredictor())
+    {
+        solve
+        (
+            UEqn
+         ==
+            fvc::reconstruct
+            (
+                (
+                    mixture.surfaceTensionForce()
+                  - ghf*fvc::snGrad(rho)
+                  - fvc::snGrad(p_rgh)
+                ) * mesh.magSf()
+            )
+        );
+
+        fvOptions.correct(U);
+    }