add two-phase source term in dfHighSpeedFoam and add fluxSchemes

Allwmake

@@ -17,6 +17,7 @@ src/lagrangian/Allwmake $targetType $*
 wmake src/dfCombustionModels
 wmake src/dynamicMesh
 wmake src/dynamicFvMesh
+wmake src/fluxSchemes
 
 wmake applications/solvers/df0DFoam
 wmake applications/solvers/dfLowMachFoam

applications/solvers/dfHighSpeedFoam/CMakeLists.txt

@@ -110,4 +110,4 @@ endif()
 
 # install
 set(CMAKE_INSTALL_PREFIX ${DF_ROOT})
-install(TARGETS ${PROJECT_NAME} DESTINATION bin)
+install(TARGETS ${PROJECT_NAME} DESTINATION bin)

applications/solvers/dfHighSpeedFoam/Make/options

@@ -7,21 +7,33 @@ EXE_INC = -std=c++14 \
     $(PFLAGS) $(PINC) \
     $(if $(LIBTORCH_ROOT),-DUSE_LIBTORCH,) \
     $(if $(PYTHON_INC_DIR),-DUSE_PYTORCH,) \
-    -I$(CANTERA_ROOT)/include \
-    -I$(DF_SRC)/dfCanteraMixture/lnInclude \
-    -I$(DF_SRC)/dfChemistryModel/lnInclude \
-    -I$(DF_SRC)/dfCombustionModels/lnInclude \
     -I$(FOAM_APP)/solvers/compressible/rhoCentralFoam/BCs/lnInclude \
-    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
+    -I$(DF_SRC)/fluxSchemes/lnInclude \
     -I$(LIB_SRC)/finiteVolume/cfdTools \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/Pstream/mpi \
-    -I$(LIB_SRC)/sampling/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
-    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/Pstream/mpi \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(DF_SRC)/lagrangian/intermediate/lnInclude \
+    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
+    -I$(DF_SRC)/lagrangian/spray/lnInclude \
+    -I$(LIB_SRC)/lagrangian/spray/lnInclude \
+    -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
+    -I$(DF_SRC)/thermophysicalModels/SLGThermo/lnInclude \
+    -I$(DF_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
+    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
+    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
+    -I$(DF_SRC)/dfCanteraMixture/lnInclude \
+    -I$(DF_SRC)/dfChemistryModel/lnInclude \
+    -I$(DF_SRC)/dfCombustionModels/lnInclude \
+    -I$(CANTERA_ROOT)/include \
     $(if $(LIBTORCH_ROOT),-I$(LIBTORCH_ROOT)/include,) \
     $(if $(LIBTORCH_ROOT),-I$(LIBTORCH_ROOT)/include/torch/csrc/api/include,) \
     $(PYTHON_INC_DIR)
@@ -30,6 +42,13 @@ EXE_LIBS = \
     -lfiniteVolume \
     -lcompressibleTransportModels \
     -lturbulenceModels \
+    -llagrangian \
+    -lregionModels \
+    -ldfSurfaceFilmModels \
+    -ldfThermophysicalProperties \
+    -ldfSLGThermo \
+    -ldfLagrangianTurbulence \
+    -ldfLagrangianSpray \
     -ldynamicFvMesh \
     -ltopoChangerFvMesh \
     -lmeshTools \
@@ -39,6 +58,8 @@ EXE_LIBS = \
     -ldfCanteraMixture \
     -ldfChemistryModel \
     -ldfCombustionModels  \
+    -ldfLagrangianIntermediate \
+    -lfluxSchemes \
     $(CANTERA_ROOT)/lib/libcantera.so \
     $(if $(LIBTORCH_ROOT),$(LIBTORCH_ROOT)/lib/libtorch.so,) \
     $(if $(LIBTORCH_ROOT),$(LIBTORCH_ROOT)/lib/libc10.so,) \

applications/solvers/dfHighSpeedFoam/calculateR.H

@@ -15,7 +15,7 @@ forAll(Y, i)
 {
     volScalarField& Yi = Y[i];
 
-    if (i != inertIndex)
+    // if (i != inertIndex)
     {          
         rhoYi[i].ref() += chemistry->RR(i)*runTime.deltaT();
         Info <<"max reaction rate "<< Yi.name() << " is " << max(chemistry->RR(i)).value() << endl;
@@ -29,9 +29,17 @@ forAll(Y, i)
     }
 }
 
-Y[inertIndex] = scalar(1) - Yt;
-Y[inertIndex].max(0.0);
-rhoYi[inertIndex] = rho*Y[inertIndex];
+// Y[inertIndex] = scalar(1) - Yt;
+// Y[inertIndex].max(0.0);
+// rhoYi[inertIndex] = rho*Y[inertIndex];
+
+forAll(Y, i)
+{
+    Y[i] = Y[i]/Yt;
+    Y[i].max(0.0);
+    rhoYi[i] = rho*Y[i];
+}
+
 
 end = std::clock();
 time_monitor_Y += double(end - start) / double(CLOCKS_PER_SEC);

applications/solvers/dfHighSpeedFoam/centralCourantNo.H

@@ -28,6 +28,80 @@ Description
     Calculates the mean and maximum wave speed based Courant Numbers.
 
 \*---------------------------------------------------------------------------*/
+surfaceScalarField rho_pos(interpolate(rho, pos));
+surfaceScalarField rho_neg(interpolate(rho, neg));
+
+surfaceVectorField rhoU_pos(interpolate(rhoU, pos, U.name()));
+surfaceVectorField rhoU_neg(interpolate(rhoU, neg, U.name()));
+
+volScalarField rPsi("rPsi", 1.0/psi);
+volScalarField c("c", sqrt(thermo.Cp()/thermo.Cv()*rPsi));
+surfaceScalarField rPsi_pos(interpolate(rPsi, pos, T.name()));
+surfaceScalarField rPsi_neg(interpolate(rPsi, neg, T.name()));
+
+// surfaceVectorField U_pos("U_pos", rhoU_pos/rho_pos);
+// surfaceVectorField U_neg("U_neg", rhoU_neg/rho_neg);
+
+surfaceVectorField U_pos(interpolate(U, pos, U.name()));
+surfaceVectorField U_neg(interpolate(U, neg, U.name()));
+
+surfaceScalarField phiv_pos("phiv_pos", U_pos & mesh.Sf());
+surfaceScalarField phiv_neg("phiv_neg", U_neg & mesh.Sf());
+
+// Make fluxes relative to mesh-motion
+if (mesh.moving())
+{
+  phiv_pos -= mesh.phi();
+  phiv_neg -= mesh.phi();
+}
+
+surfaceScalarField cSf_pos
+(
+  "cSf_pos",
+  interpolate(c, pos, T.name())*mesh.magSf()
+);
+surfaceScalarField cSf_neg
+(
+  "cSf_neg",
+  interpolate(c, neg, T.name())*mesh.magSf()
+);
+
+surfaceScalarField ap
+(
+  "ap",
+  max(max(phiv_pos + cSf_pos, phiv_neg + cSf_neg), v_zero)
+);
+surfaceScalarField am
+(
+  "am",
+  min(min(phiv_pos - cSf_pos, phiv_neg - cSf_neg), v_zero)
+);
+
+surfaceScalarField a_pos("a_pos", ap/(ap - am));
+
+surfaceScalarField amaxSf("amaxSf", max(mag(am), mag(ap)));
+
+surfaceScalarField aSf("aSf", am*a_pos);
+
+word fluxScheme("Kurganov");
+mesh.schemesDict().readIfPresent("fluxScheme", fluxScheme);
+if (fluxScheme == "Tadmor")
+{
+  aSf = -0.5*amaxSf;
+  a_pos = 0.5;
+}
+
+surfaceScalarField a_neg("a_neg", 1.0 - a_pos);
+
+phiv_pos *= a_pos;
+phiv_neg *= a_neg;
+
+surfaceScalarField aphiv_pos("aphiv_pos", phiv_pos - aSf);
+surfaceScalarField aphiv_neg("aphiv_neg", phiv_neg + aSf);
+
+// Reuse amaxSf for the maximum positive and negative fluxes
+// estimated by the central scheme
+amaxSf = max(mag(aphiv_pos), mag(aphiv_neg));
 
 if (mesh.nInternalFaces())
 {

applications/solvers/dfHighSpeedFoam/createClouds.H

@@ -0,0 +1,9 @@
+Info<< "\nConstructing reacting cloud" << endl;
+basicSprayCloud parcels
+(
+    "sprayCloud",
+    rho,
+    U,
+    g,
+    slgThermo
+);

applications/solvers/dfHighSpeedFoam/createFields.H

@@ -1,17 +1,24 @@
+#include "readGravitationalAcceleration.H"
 #include "createRDeltaT.H"
 
 word ddtSchemes("Euler");
 if (mesh.schemesDict().readIfPresent("timeScheme", ddtSchemes))
 {
-    if(ddtSchemes == "RK2SSP" || ddtSchemes == "RK3SSP")
+    if(ddtSchemes == "RK2SSP" || ddtSchemes == "RK3SSP" || ddtSchemes == "Euler")
     {
         Info<< "ddtSchemes: " << ddtSchemes << endl;
+        if(ddtSchemes == "RK2SSP" || ddtSchemes == "RK3SSP")
+        {
+            Info<< "!! Note: RK2SSP and RK3SSP are not available for two-phase flow simulaiton. "
+                << "If you want to simulate two-phase flows, please change time scheme to 'Euler' ."
+                << endl;
+        }
     }
     else
     {
         FatalErrorInFunction
             << "This timeScheme is not a valid choice. "
-            << "Please use RK2SSP or RK3SSP scheme."
+            << "Please use Euler, RK2SSP or RK3SSP scheme."
             << abort(FatalError);
     }
 }
@@ -44,14 +51,15 @@ if((ddtSchemes != "RK2SSP") && (ddtSchemes != "RK3SSP"))
     }
 }
 
-
 Info<< "Reading thermophysical properties\n" << endl;
 
 // psiThermo* pThermo = new hePsiThermo<psiThermo, CanteraMixture>(mesh, word::null);
 // psiThermo& thermo = *pThermo;
 rhoThermo* pThermo = new heRhoThermo<rhoThermo, CanteraMixture>(mesh, word::null);
 rhoThermo& thermo = *pThermo;
 
+SLGThermo slgThermo(mesh, thermo);
+
 //move from creatFieldRefs.H to createFields.H
 //p needed to be created before e
 volScalarField& p = thermo.p();
@@ -116,6 +124,51 @@ volVectorField rhoU
     rho*U
 );
 
+surfaceScalarField phi("phi", fvc::flux(rhoU));
+
+Info<< "Creating flux properties\n" << endl;
+surfaceScalarField rhoPhi
+(
+    IOobject
+    (
+        "rhoPhi",
+        mesh.time().timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("0", dimDensity*dimVelocity*dimArea, 0.0)
+);
+
+surfaceVectorField rhoUPhi
+(
+    IOobject
+    (
+        "rhoUPhi",
+        mesh.time().timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedVector("0", dimDensity*sqr(dimVelocity)*dimArea, Zero)
+);
+
+surfaceScalarField rhoEPhi
+(
+    IOobject
+    (
+        "rhoEPhi",
+        mesh.time().timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("0", dimDensity*pow3(dimVelocity)*dimArea, 0.0)
+);
+
 surfaceScalarField pos
 (
     IOobject
@@ -140,7 +193,7 @@ surfaceScalarField neg
     dimensionedScalar(dimless, -1.0)
 );
 
-surfaceScalarField phi("phi", fvc::flux(rhoU));
+autoPtr<Foam::fluxScheme> fluxSchemeFields(Foam::fluxScheme::New(mesh));
 
 Info<< "Creating turbulence model\n" << endl;
 autoPtr<compressible::turbulenceModel> turbulence
@@ -153,6 +206,7 @@ autoPtr<compressible::turbulenceModel> turbulence
         thermo
     )
 );
+
 const word turbName(mesh.objectRegistry::lookupObject<IOdictionary>("turbulenceProperties").lookup("simulationType"));
 
 multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
@@ -186,15 +240,8 @@ volScalarField rhoE
     rho*(ea + 0.5*magSqr(U))
 );
 
-const word combModelName(mesh.objectRegistry::lookupObject<IOdictionary>("combustionProperties").lookup("combustionModel"));
-Info << "Combustion Model Name is confirmed as "<< combModelName << endl;
-
-
-if(combModelName!="ESF" && combModelName!="flareFGM"  && combModelName!="DeePFGM" && combModelName!="FSD")
-{
-    chemistry->correctThermo();
-    Info<< "At initial time, min/max(T) = " << min(T).value() << ", " << max(T).value() << endl;
-}
+chemistry->correctThermo();
+Info<< "At initial time, min/max(T) = " << min(T).value() << ", " << max(T).value() << endl;
 
 forAll(Y, i)
 {
@@ -224,6 +271,9 @@ forAll(rhoYi,i)
     );
 }
 
+#include "createMRF.H"
+#include "createClouds.H"
+
 const scalar Sct = chemistry->lookupOrDefault("Sct", 1.);
 volScalarField diffAlphaD
 (
@@ -264,3 +314,29 @@ volVectorField sumYDiffError
     mesh,
     dimensionedVector("sumYDiffError", dimDynamicViscosity/dimLength, Zero)
 );
+
+volScalarField gradP
+(
+    IOobject
+    (
+        "gradP",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mag(fvc::grad(thermo.p()))
+);
+
+volScalarField pMax
+(
+    IOobject
+    (
+        "pMax",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    thermo.p()
+);