ENH: Up to date icoReactingMultiphaseInterFoam solver and libs

Adding tutorials and other minor changes

applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/Allwclean

 #!/bin/sh
 cd ${0%/*} || exit 1    # Run from this directory
-set -x
 
 wclean libso phasesSystem
-wclean libso meltingEvaporationModel
+wclean libso massTransferModels
 wclean libso CompressibleMultiPhaseTurbulenceModels
 wclean libso laserDTRM
 wclean

applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/Allwmake

 #!/bin/sh
 cd ${0%/*} || exit 1    # Run from this directory
-set -x
 
-wmakeLnInclude meltingEvaporationModel
-wmake libso phasesSystem
-wmake libso meltingEvaporationModel
-wmake libso CompressibleMultiPhaseTurbulenceModels
-wmake libso laserDTRM
+targetType=libso
+. $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments
+
+wmakeLnInclude massTransferModels
+wmake phasesSystem
+wmake massTransferModels
+wmake CompressibleMultiPhaseTurbulenceModels
+wmake laserDTRM
 wmake
 
 

applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/Make/options

@@ -28,7 +28,8 @@ EXE_LIBS = \
     -lfluidThermophysicalModels \
     -lIncompressibleMultiphaseSystems \
     -lCompressibleMultiPhaseTurbulenceModels \
-    -lmeltingEvaporationModels \
+    -lmassTransferModels \
     -lsolidThermo \
     -lsolidSpecie \
-    -ltwoPhaseProperties
+    -ltwoPhaseProperties \
+    -llaserDTRM

applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/TEqn.H

 {
     radiation->correct();
+    rhoCp =  rho*fluid.Cp();
 
-    tmp<volScalarField> texpSource
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "texpSource",
-                mesh.time().timeName(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh,
-            dimensionedScalar("zero", dimTemperature/dimTime, 0),
-            zeroGradientFvPatchScalarField::typeName
-        )
-    );
-    volScalarField& expSource = texpSource.ref();
+    const surfaceScalarField rhoCpPhi(fvc::interpolate(fluid.Cp())*rhoPhi);
 
-    tmp<volScalarField> tkappaEff
+    const volScalarField kappaEff
     (
-        new volScalarField
-        (
-            IOobject
-            (
-                "kappaEff",
-                mesh.time().timeName(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh,
-            dimensionedScalar("zero", sqr(dimLength)/dimTime, 0),
-            zeroGradientFvPatchScalarField::typeName
-        )
+        "kappaEff",
+        fluid.kappa() + fluid.Cp()*turbulence->mut()/fluid.Prt()
     );
 
-    volScalarField& kappaEff = tkappaEff.ref();
-
-    //const surfaceScalarField rhoTempPhi("phi", rhoPhi/fvc::interpolate(rho));
-
-    const surfaceScalarField rhoTempPhi("phi", fluid.phi());
-
-    const volScalarField divU(fvc::div(phi));
-
-    forAllIter(UPtrList<phaseModel>, fluid.phases(), iter)
-    {
-        phaseModel& phase = iter();
-        const volScalarField& alpha = phase;
-        const volScalarField DDtAlpha(fvc::DDt(phi, alpha));
-
-        const volScalarField invCpRho(1.0/phase.rho()/phase.Cp());
-
-        if (fluid.dpdt())
-        {
-            const volScalarField ddtp(fvc::ddt(p));
-            expSource += (DDtAlpha*p + alpha*(p*divU + ddtp))*invCpRho;
-        }
-        else
-        {
-            expSource += (DDtAlpha*p + alpha*(p*divU))*invCpRho;
-        }
-
-        kappaEff += alpha*phase.kappa()*invCpRho;
-
-        DebugVar(max(alpha*phase.kappa()));
-        DebugVar(max(alpha*phase.Cp()));
-    }
-
-    kappaEff += turbulence->nut()/fluid.Prt();
-
-    if (mesh.time().outputTime())
-    {
-        expSource.write();
-        kappaEff.write();
-    }
-
-    //dimensionedScalar S("S", dimEnergy/dimVolume/dimTime, 1.225e8);
-
     fvScalarMatrix TEqn
     (
-        fvm::ddt(T)
-      + fvm::div(rhoTempPhi, T)
+        fvm::ddt(rhoCp, T)
+      + fvm::div(rhoCpPhi, T, "div(phi,T)")
+      - fvm::Sp(fvc::ddt(rhoCp) + fvc::div(rhoCpPhi), T)
       - fvm::laplacian(kappaEff, T,  "laplacian(kappa,T)")
       ==
         fluid.heatTransfer(T)
-      + expSource
-      + radiation->ST(fluid.Cp()*rho, T)
-//      + S/Cp/rho
+      + radiation->ST(T)
+      + fvOptions(rhoCp, T)
     );
 
     TEqn.relax();
+
+    fvOptions.constrain(TEqn);
+
     TEqn.solve();
 
+    fvOptions.correct(T);
+
     fluid.correct();
 
     Info<< "min/max(T) = "

applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/UEqn.H

-    fvVectorMatrix UEqn
-    (
-        fvm::ddt(rho, U)
-      + fvm::div(rhoPhi, U)
-      + turbulence->divDevRhoReff(U)
-      ==
-        fvOptions(rho, U)
-    );
-
-    UEqn.relax();
+fvVectorMatrix UEqn
+(
+    fvm::ddt(rho, U)
+  + fvm::div(rhoPhi, U)
+    //- fvm::Sp(fvc::ddt(rho) + fvc::div(rhoPhi), U)
+  + turbulence->divDevRhoReff(U)
+  ==
+    fvOptions(rho, U)
+);
 
-    fluid.addInterfacePorosity(UEqn);
+UEqn.relax();
 
-    if (pimple.momentumPredictor())
-    {
+fluid.addInterfacePorosity(UEqn);
 
-
-        solve
+if (pimple.momentumPredictor())
+{
+    solve
+    (
+        UEqn
+      ==
+        fvc::reconstruct
         (
-            UEqn
-         ==
-            fvc::reconstruct
             (
-                (
-                    fluid.surfaceTensionForce()
-                  - ghf*fvc::snGrad(rho)
-                  - fvc::snGrad(p_rgh)
-                ) * mesh.magSf()
-            )
-        );
+                fluid.surfaceTensionForce()
+                - ghf*fvc::snGrad(rho)
+                - fvc::snGrad(p_rgh)
+            ) * mesh.magSf()
+        )
+    );
 
-        fvOptions.correct(U);
-        K = 0.5*magSqr(U);
-    }
+    fvOptions.correct(U);
+    K = 0.5*magSqr(U);
+}

applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/YEqns.H

 {
-    // Semi-implicit mass transfer for species
-    // Initilize dmdt for alpha Eq's for mass transfers driven by species
-    autoPtr<phaseSystem::massTransferTable>
-        massTransferPtr(fluid.massTransfer(T));
-
-    //phaseSystem::massTransferTable& massTransfer(massTransferPtr());
-
     forAllIter(UPtrList<phaseModel>, fluid.phases(), iter)
     {
         phaseModel& phase = iter();
         PtrList<volScalarField>& Y = phase.Y();
-        //const surfaceScalarField& alphaPhi = phase.alphaPhi();
 
         if (!Y.empty())
         {
-            const volScalarField& alpha = phase;
-
-            label inertIndex = -1;
-            volScalarField Yt(0.0*Y[0]);
+            //- Su phase source terms
+            PtrList<volScalarField::Internal> Sus(Y.size());
+            //- Sp phase source terms
+            PtrList<volScalarField::Internal> Sps(Y.size());
 
-            forAll(Y, i)
+            forAll(Sus, i)
             {
-                tmp<fvScalarMatrix> YiEqn(phase.YiEqn(Y[i]));
-
-                if (YiEqn.valid())
-                {
-                    YiEqn.ref() =
+                Sus.set
+                (
+                    i,
+                    new volScalarField::Internal
                     (
-                        YiEqn()
-                      - fvm::laplacian
+                        IOobject
                         (
-                            alpha*turbulence->nuEff(),
-                            Y[i]
-                        )
-                    ==
-                    //  (*massTransfer[Y[i].name()])(/phase.rho()
-                        fvc::ddt(alpha)
-                       *pos
+                            "Su" + phase.name(),
+                            mesh.time().timeName(),
+                            mesh
+                        ),
+                        mesh,
+                        dimensionedScalar("zero", dimless/dimTime, 0.0)
+                    )
+                );
+                Sps.set
+                (
+                    i,
+                    new volScalarField::Internal
+                    (
+                        IOobject
                         (
-                           fluid.dmdtYi(Y[i].name())
-                        - dimensionedScalar("zero", dimDensity/dimTime, 1e-3)
-                        )
-                        //fluid.dmdtYi(Y[i].name())/phase.rho()
-                       //explicit mass sources (P or T)
-                    );
-
-                    YiEqn->relax();
-                    YiEqn->solve(mesh.solver("Yi"));
-                    Y[i].max(0.0);
-                    Y[i].min(1.0);
-                    Yt += Y[i];
-                }
-                else
-                {
-                    inertIndex = i;
-                }
-
-//                 if (mesh.time().outputTime())
-//                 {
-//                     volScalarField dmdtYi("dmdtYi", pos(fluid.dmdtYi(Y[i].name())));
-//                     dmdtYi.write();
-//                 }
+                            "Sp" + phase.name(),
+                            mesh.time().timeName(),
+                            mesh
+                        ),
+                        mesh,
+                        dimensionedScalar("zero", dimless/dimTime, 0.0)
+                    )
+                );
+            }
 
-                Info << "Min/Max : " << min(Y[i]) << " " << max(Y[i]) << endl;
-                Info<< "Max dmdtYi : "
-                    << max(fluid.dmdtYi(Y[i].name())().internalField()) << endl;
-                Info<< "Min dmdtYi : "
-                    << min(fluid.dmdtYi(Y[i].name())().internalField()) << endl;
+            forAll(Y, i)
+            {
+                // Calculate mass exchange for species consistent with
+                // alpha's source terms.
+                fluid.massSpeciesTransfer(phase, Sus[i], Sps[i], Y[i].name());
             }
-            Y[inertIndex] = scalar(1) - Yt;
-            Y[inertIndex].max(0.0);
+            phase.solveYi(Sus, Sps);
         }
     }
 }

applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/alphaCourantNo.H

@@ -34,8 +34,20 @@ scalar maxAlphaCo
     readScalar(runTime.controlDict().lookup("maxAlphaCo"))
 );
 
+scalar maxAlphaDdt
+(
+    runTime.controlDict().lookupOrDefault("maxAlphaDdt", GREAT)
+);
+
+scalar maxDi
+(
+    runTime.controlDict().lookupOrDefault<scalar>("maxDi", GREAT)
+);
+
 scalar alphaCoNum = 0.0;
 scalar meanAlphaCoNum = 0.0;
+scalar ddtAlphaNum = 0.0;
+scalar DiNum = 0.0;
 
 if (mesh.nInternalFaces())
 {
@@ -49,9 +61,17 @@ if (mesh.nInternalFaces())
 
     meanAlphaCoNum =
         0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();
+
+    ddtAlphaNum = fluid.ddtAlphaMax().value()*runTime.deltaTValue();
+
+    DiNum = fluid.maxDiffNo();
 }
 
 Info<< "Interface Courant Number mean: " << meanAlphaCoNum
     << " max: " << alphaCoNum << endl;
 
+Info<< "Maximum ddtAlpha : " << ddtAlphaNum << endl;
+
+Info<< "Maximum DiffNum : " << DiNum << endl;
+
 // ************************************************************************* //

applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/createFields.H

@@ -28,8 +28,9 @@
 
 
     Info<< "Calculating field g.h\n" << endl;
-    volScalarField gh("gh", g & mesh.C());
-    surfaceScalarField ghf("ghf", g & mesh.Cf());
+    #include "readGravitationalAcceleration.H"
+    #include "readhRef.H"
+    #include "gh.H"
 
     volScalarField  p
     (
@@ -41,7 +42,7 @@
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        p_rgh //+ rho*gh
+        p_rgh
     );
 
     Info<< "Creating multiphaseSystem\n" << endl;
@@ -49,7 +50,13 @@
 
     multiphaseSystem& fluid = fluidPtr();
 
-    //volScalarField& e = fluid.he();
+    if (!fluid.incompressible())
+    {
+        FatalError << "One or more phases are not incompressible. " << nl
+            << "This is a incompressible solver." << abort(FatalError);
+    }
+
+
     volScalarField& T = fluid.T();
 
     // Need to store rho for ddt(rho, U)
@@ -115,23 +122,22 @@
     (
         radiation::radiationModel::New(T)
     );
-/*
-    Info<< "Calculating field kappaEff\n" << endl;
-    volScalarField kappaEff
+
+    Info<< "Calculating field rhoCp\n" << endl;
+    volScalarField rhoCp
     (
         IOobject
         (
-            "kappaEff",
+            "rhoCp",
             runTime.timeName(),
             mesh,
             IOobject::NO_READ,
-            IOobject::AUTO_WRITE
+            IOobject::NO_WRITE
         ),
-        fluid.kappa()
+        fluid.rho()*fluid.Cp()
     );
+    rhoCp.oldTime();
 
-    kappaEff.correctBoundaryConditions();
-*/
 
     Info<< "Creating field kinetic energy K\n" << endl;
     volScalarField K("K", 0.5*magSqr(U));

applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/icoReactingMultiphaseInterFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2015 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -29,8 +29,7 @@ Group
 
 Description
     Solver for n incompressible, non-isothermal immiscible fluids with
-    phase-change (evaporation-condensation).  Uses a VOF (volume of fluid)
-    phase-fraction based interface capturing approach.
+    phase-change.  Uses a VOF (volume of fluid) phase-fraction based interface capturing approach.
 
     The momentum, energy and other fluid properties are of the "mixture" and a
     single momentum equation is solved.
@@ -40,7 +39,6 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "CMULES.H"
 #include "subCycle.H"
 #include "multiphaseSystem.H"
 #include "turbulentFluidThermoModel.H"
@@ -50,6 +48,7 @@ Description
 #include "radiationModel.H"
 #include "HashPtrTable.H"
 #include "fvcDDt.H"
+#include "zeroField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -61,7 +60,6 @@ int main(int argc, char *argv[])
 
     pimpleControl pimple(mesh);
 
-    #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "createFvOptions.H"
     #include "createTimeControls.H"

applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/laserDTRM/DTRMParticle/DTRMParticle.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2015 OpenCFD Ltd
+    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -28,6 +28,7 @@ License
 #include "physicoChemicalConstants.H"
 
 using namespace Foam::constant;
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::DTRMParticle::DTRMParticle
@@ -38,19 +39,40 @@ Foam::DTRMParticle::DTRMParticle
     const scalar I,
     const label cellI,
     const scalar dA,
-    const label reflectedId,
-    const scalar Imin,
-    bool doCellFacePt
+    const label transmissiveId
+)
+:
+    particle(mesh, position, cellI),
+    p0_(position),
+    p1_(targetPosition),
+    I0_(I),
+    I_(I),
+    dA_(dA),
+    transmissiveId_(transmissiveId)
+{}
+
+
+Foam::DTRMParticle::DTRMParticle
+(
+    const polyMesh& mesh,
+    const barycentric& coordinates,
+    const label celli,
+    const label tetFacei,
+    const label tetPti,
+    const vector& position,
+    const vector& targetPosition,
+    const scalar I,
+    const scalar dA,
+    const label transmissiveId
 )
 :
-    particle(mesh, position, cellI, doCellFacePt),
+    particle(mesh, coordinates, celli, tetFacei, tetPti),
     p0_(position),
     p1_(targetPosition),
     I0_(I),
     I_(I),
     dA_(dA),
-    reflectedId_(reflectedId),
-    Imin_(Imin)
+    transmissiveId_(transmissiveId)
 {}
 
 
@@ -62,14 +84,14 @@ Foam::DTRMParticle::DTRMParticle(const DTRMParticle& p)
     I0_(p.I0_),
     I_(p.I_),
     dA_(p.dA_),
-    reflectedId_(p.reflectedId_),
-    Imin_(p.Imin_)
+    transmissiveId_(p.transmissiveId_)
 {}
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 bool Foam::DTRMParticle::move
 (
+    Cloud<DTRMParticle>& spc,
     trackingData& td,
     const scalar trackTime
 )
@@ -77,106 +99,120 @@ bool Foam::DTRMParticle::move
     td.switchProcessor = false;
     td.keepParticle = true;
 
-    const polyBoundaryMesh& pbMesh = mesh_.boundaryMesh();
-
-    while (td.keepParticle && !td.switchProcessor)
+    while (td.keepParticle && !td.switchProcessor && stepFraction() < 1)
     {
-        point p0 = position();
-        label cell0 = cell();
-
-        scalar dt = trackToFace(p1_, td);
-
-        // Consider the cell between f0(start of tracking) and f1
-        label celli = cell();
-
-        const vector dsv = position() - p0;
-        const scalar ds = mag(dsv);