From 6ff341b5f5c33986d8fcb70b32c31d11ebb6c0e7 Mon Sep 17 00:00:00 2001
From: sergio <sergio>
Date: Mon, 19 Jun 2017 11:10:19 -0700
Subject: [PATCH] Adding alphaEqn.H with interpolation method. Adding special
 alphaCourantNo for overlaping Adding bounded term to UEq.H for
 overInterDyMFoam Changing to NO_WRITE for the cellMask field Changing
 twoSimpleRotors tutorial to open domain

---
 .../interFoam/overInterDyMFoam/UEqn.H         |   1 +
 .../overInterDyMFoam/alphaCourantNo.H         |  59 ++++++++
 .../interFoam/overInterDyMFoam/alphaEqn.H     | 131 ++++++++++++++++++
 src/overset/include/createCellMask.H          |   2 +-
 .../twoSimpleRotors/0.orig/U                  |  22 +--
 .../twoSimpleRotors/0.orig/p                  |  29 ++--
 .../twoSimpleRotors/system/blockMeshDict      |  16 +++
 .../twoSimpleRotors/system/controlDict        |   2 +-
 .../background/system/controlDict             |   4 +-
 .../floatingBody/background/system/fvSolution |   4 +-
 10 files changed, 233 insertions(+), 37 deletions(-)
 create mode 100644 applications/solvers/multiphase/interFoam/overInterDyMFoam/alphaCourantNo.H
 create mode 100644 applications/solvers/multiphase/interFoam/overInterDyMFoam/alphaEqn.H

diff --git a/applications/solvers/multiphase/interFoam/overInterDyMFoam/UEqn.H b/applications/solvers/multiphase/interFoam/overInterDyMFoam/UEqn.H
index 600c94829cc..fdfe1347397 100644
--- a/applications/solvers/multiphase/interFoam/overInterDyMFoam/UEqn.H
+++ b/applications/solvers/multiphase/interFoam/overInterDyMFoam/UEqn.H
@@ -3,6 +3,7 @@
     fvVectorMatrix UEqn
     (
         fvm::ddt(rho, U) + fvm::div(rhoPhi, U)
+      - fvm::Sp(fvc::ddt(rho) + fvc::div(rhoPhi), U)
       + MRF.DDt(rho, U)
       + turbulence->divDevRhoReff(rho, U)
      ==
diff --git a/applications/solvers/multiphase/interFoam/overInterDyMFoam/alphaCourantNo.H b/applications/solvers/multiphase/interFoam/overInterDyMFoam/alphaCourantNo.H
new file mode 100644
index 00000000000..8666b664ac1
--- /dev/null
+++ b/applications/solvers/multiphase/interFoam/overInterDyMFoam/alphaCourantNo.H
@@ -0,0 +1,59 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Global
+    alphaCourantNo
+
+Description
+    Calculates and outputs the mean and maximum Courant Numbers.
+
+\*---------------------------------------------------------------------------*/
+
+scalar maxAlphaCo
+(
+    readScalar(runTime.controlDict().lookup("maxAlphaCo"))
+);
+
+scalar alphaCoNum = 0.0;
+scalar meanAlphaCoNum = 0.0;
+
+if (mesh.nInternalFaces())
+{
+    surfaceScalarField phiMask(localMin<scalar>(mesh).interpolate(cellMask));
+
+    scalarField sumPhi
+    (
+        mixture.nearInterface()().internalField()
+       *fvc::surfaceSum(mag(phiMask*phi))().internalField()
+    );
+
+    alphaCoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
+
+    meanAlphaCoNum =
+        0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();
+}
+
+Info<< "Interface Courant Number mean: " << meanAlphaCoNum
+    << " max: " << alphaCoNum << endl;
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/interFoam/overInterDyMFoam/alphaEqn.H b/applications/solvers/multiphase/interFoam/overInterDyMFoam/alphaEqn.H
new file mode 100644
index 00000000000..ac40176f08e
--- /dev/null
+++ b/applications/solvers/multiphase/interFoam/overInterDyMFoam/alphaEqn.H
@@ -0,0 +1,131 @@
+{
+    word alphaScheme("div(phi,alpha)");
+    word alpharScheme("div(phirb,alpha)");
+
+    // Standard face-flux compression coefficient
+    surfaceScalarField phic(mixture.cAlpha()*mag(phi/mesh.magSf()));
+
+    // Add the optional isotropic compression contribution
+    if (icAlpha > 0)
+    {
+        phic *= (1.0 - icAlpha);
+        phic += (mixture.cAlpha()*icAlpha)*fvc::interpolate(mag(U));
+    }
+
+    surfaceScalarField::Boundary& phicBf =
+        phic.boundaryFieldRef();
+
+    // Do not compress interface at non-coupled boundary faces
+    // (inlets, outlets etc.)
+    forAll(phic.boundaryField(), patchi)
+    {
+        fvsPatchScalarField& phicp = phicBf[patchi];
+
+        if (!phicp.coupled())
+        {
+            phicp == 0;
+        }
+    }
+
+    tmp<surfaceScalarField> tphiAlpha;
+
+    if (MULESCorr)
+    {
+        mesh.interpolate(alpha1);
+
+        fvScalarMatrix alpha1Eqn
+        (
+            fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
+          + fv::gaussConvectionScheme<scalar>
+            (
+                mesh,
+                phi,
+                upwind<scalar>(mesh, phi)
+            ).fvmDiv(phi, alpha1)
+        );
+
+        alpha1Eqn.solve();
+
+        tmp<surfaceScalarField> tphiAlphaUD(alpha1Eqn.flux());
+        alphaPhi = tphiAlphaUD();
+
+        if (alphaApplyPrevCorr && tphiAlphaCorr0.valid())
+        {
+            MULES::correct(alpha1, alphaPhi, tphiAlphaCorr0.ref(), 1, 0);
+
+            alphaPhi += tphiAlphaCorr0();
+        }
+
+        // Cache the upwind-flux
+        tphiAlphaCorr0 = tphiAlphaUD;
+
+        alpha2 = 1.0 - alpha1;
+
+        mixture.correct();
+    }
+
+    for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
+    {
+        surfaceScalarField phir(phic*mixture.nHatf());
+
+        tmp<surfaceScalarField> tphiAlphaUn
+        (
+            fvc::flux
+            (
+                phi,
+                alpha1,
+                alphaScheme
+            )
+          + fvc::flux
+            (
+               -fvc::flux(-phir, alpha2, alpharScheme),
+                alpha1,
+                alpharScheme
+            )
+        );
+
+        if (MULESCorr)
+        {
+            tmp<surfaceScalarField> tphiAlphaCorr(tphiAlphaUn() - alphaPhi);
+            volScalarField alpha10("alpha10", alpha1);
+
+            MULES::correct(alpha1, tphiAlphaUn(), tphiAlphaCorr.ref(), 1, 0);
+            mesh.interpolate(alpha1);
+
+            // Under-relax the correction for all but the 1st corrector
+            if (aCorr == 0)
+            {
+                alphaPhi += tphiAlphaCorr();
+            }
+            else
+            {
+                alpha1 = 0.5*alpha1 + 0.5*alpha10;
+                alphaPhi += 0.5*tphiAlphaCorr();
+            }
+        }
+        else
+        {
+            alphaPhi = tphiAlphaUn;
+            MULES::explicitSolve(alpha1, phi, alphaPhi, 1, 0);
+            mesh.interpolate(alpha1);
+        }
+
+        alpha2 = 1.0 - alpha1;
+
+        mixture.correct();
+    }
+
+
+    rhoPhi = alphaPhi*(rho1 - rho2) + phi*rho2;
+
+    if (alphaApplyPrevCorr && MULESCorr)
+    {
+        tphiAlphaCorr0 = alphaPhi - tphiAlphaCorr0;
+    }
+
+    Info<< "Phase-1 volume fraction = "
+        << alpha1.weightedAverage(mesh.Vsc()).value()
+        << "  Min(alpha1) = " << min(alpha1).value()
+        << "  Max(alpha1) = " << max(alpha1).value()
+        << endl;
+}
diff --git a/src/overset/include/createCellMask.H b/src/overset/include/createCellMask.H
index 61de1068858..1351e1738fa 100644
--- a/src/overset/include/createCellMask.H
+++ b/src/overset/include/createCellMask.H
@@ -41,7 +41,7 @@ volScalarField cellMask
         runTime.timeName(),
         mesh,
         IOobject::NO_READ,
-        IOobject::AUTO_WRITE
+        IOobject::NO_WRITE
     ),
     mesh,
     dimensionedScalar("cellMask", dimless, 1.0),
diff --git a/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/0.orig/U b/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/0.orig/U
index b1fc09c1253..45c9f5a085e 100644
--- a/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/0.orig/U
+++ b/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/0.orig/U
@@ -16,7 +16,7 @@ FoamFile
 
 dimensions      [0 1 -1 0 0 0 0];
 
-internalField   uniform (0 0 0);
+internalField   uniform (0.1 0 0);
 
 boundaryField
 {
@@ -36,16 +36,16 @@ boundaryField
         value           uniform (0 0 0);
     }
 
-//    left1
-//    {
-//        type            pressureInletOutletVelocity;
-//        value           uniform (0 0 0);
-//    }
-//    left1
-//    {
-//        type            fixedValue;
-//        value           $internalField;
-//    }
+   outlet
+   {
+       type            pressureInletOutletVelocity;
+       value           uniform (0 0 0);
+   }
+   inlet
+   {
+       type            fixedValue;
+       value           $internalField;
+   }
 //
 //    right1
 //    {
diff --git a/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/0.orig/p b/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/0.orig/p
index 55e4c8b692e..29afa2eaee6 100644
--- a/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/0.orig/p
+++ b/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/0.orig/p
@@ -28,26 +28,15 @@ boundaryField
         type            zeroGradient;
     }
 
-//    left1
-//    {
-//        type            uniformTotalPressure;
-//        pressure        table
-//        (
-//            (0 10)
-//            (1 40)
-//        );
-//        p0              40; // only used for restarts
-//        U               U;
-//        phi             phi;
-//        rho             none;
-//        psi             none;
-//        gamma           1;
-//        value           uniform 40;
-//    }
-//    left1
-//    {
-//        type            zeroGradient;
-//    }
+   outlet
+   {
+       type            fixedValue;
+       value           uniform 0;
+   }
+   inlet
+   {
+       type            zeroGradient;
+   }
 //
 //    right1
 //    {
diff --git a/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/system/blockMeshDict b/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/system/blockMeshDict
index 5456fdbf9c0..e2fe1595c63 100644
--- a/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/system/blockMeshDict
+++ b/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/system/blockMeshDict
@@ -94,7 +94,23 @@ boundary
         (
             (3 7 6 2)
             (1 5 4 0)
+        );
+    }
+
+    inlet
+    {
+        type patch;
+        faces
+        (
             (0 4 7 3)
+        );
+    }
+
+    outlet
+    {
+        type patch;
+        faces
+        (
             (2 6 5 1)
         );
     }
diff --git a/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/system/controlDict b/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/system/controlDict
index 591998cc3f6..7e93ade3783 100644
--- a/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/system/controlDict
+++ b/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/system/controlDict
@@ -24,7 +24,7 @@ startTime       0;
 
 stopAt          endTime;
 
-endTime         0.06;
+endTime         0.6;
 
 deltaT          0.00025;
 
diff --git a/tutorials/multiphase/overInterDyMFoam/floatingBody/background/system/controlDict b/tutorials/multiphase/overInterDyMFoam/floatingBody/background/system/controlDict
index a21e0935f31..86e963f08af 100644
--- a/tutorials/multiphase/overInterDyMFoam/floatingBody/background/system/controlDict
+++ b/tutorials/multiphase/overInterDyMFoam/floatingBody/background/system/controlDict
@@ -20,7 +20,7 @@ libs            ("liboverset.so");
 
 application     overInterDyMFoam ;
 
-startFrom       startTime;
+startFrom       latestTime;//startTime;
 
 startTime       0.0;
 
@@ -32,7 +32,7 @@ deltaT          0.001;
 
 writeControl    adjustableRunTime;
 
-writeInterval   0.05;
+writeInterval   0.1;
 
 purgeWrite      0;
 
diff --git a/tutorials/multiphase/overInterDyMFoam/floatingBody/background/system/fvSolution b/tutorials/multiphase/overInterDyMFoam/floatingBody/background/system/fvSolution
index d786f239265..6aff6062e0c 100644
--- a/tutorials/multiphase/overInterDyMFoam/floatingBody/background/system/fvSolution
+++ b/tutorials/multiphase/overInterDyMFoam/floatingBody/background/system/fvSolution
@@ -80,8 +80,8 @@ solvers
 PIMPLE
 {
     momentumPredictor   no;
-    nOuterCorrectors    3;
-    nCorrectors         1;
+    nOuterCorrectors    2;
+    nCorrectors         2;
     nNonOrthogonalCorrectors 0;
 
     ddtCorr         yes;
-- 
GitLab