add AUSM scheme and modify KNP scheme in dfHighSpeedFoam

applications/solvers/dfHighSpeedFoam/createFields.H

@@ -103,6 +103,19 @@ volScalarField rho
     thermo.rho()
 );
 
+volScalarField pMax
+(
+    IOobject
+    (
+        "pMax",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+    thermo.p()
+);
+
 volVectorField rhoU
 (
     IOobject
@@ -190,7 +203,7 @@ const word combModelName(mesh.objectRegistry::lookupObject<IOdictionary>("combus
 Info << "Combustion Model Name is confirmed as "<< combModelName << endl;
 
 
-if(combModelName!="ESF" && combModelName!="flareFGM"  && combModelName!="DeePFGM" && combModelName!="FSD")
+if (combModelName != "flareFGM")
 {
     chemistry->correctThermo();
     Info<< "At initial time, min/max(T) = " << min(T).value() << ", " << max(T).value() << endl;

applications/solvers/dfHighSpeedFoam/dfHighSpeedFoam.C

@@ -58,6 +58,9 @@ Description
 #include "PstreamGlobals.H"
 #include "CombustionModel.H"
 
+#include "fluxSchemes/AUSMDVFlux.H"
+#include "fluxSchemes/KNPFlux.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
@@ -104,6 +107,8 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
+        pMax = max(p, pMax);
+
         #include "readTimeControls.H"
 
         //used for AMR

applications/solvers/dfHighSpeedFoam/fluxSchemes/AUSMDVFlux.H

@@ -0,0 +1,106 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2018 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/>.
+Function
+    AUSMDVFlux
+Description
+    A function for solving convective flux based on AUSMDV, which mixes the AUSMD and AUSMV scheme
+Author
+    Daoping Zhang ([email protected])
+\*---------------------------------------------------------------------------*/
+
+void AUSMDVFlux(
+    const surfaceScalarField& meshPhi,
+    const surfaceScalarField& magSf,
+    const surfaceVectorField& normal,
+    const surfaceScalarField& rhoOwn,
+    const surfaceScalarField& rhoNei,
+    const surfaceVectorField& UOwn,
+    const surfaceVectorField& UNei,
+    const surfaceScalarField& HOwn,
+    const surfaceScalarField& HNei,
+    const surfaceScalarField& pOwn,
+    const surfaceScalarField& pNei,
+    const surfaceScalarField& gammaOwn,
+    const surfaceScalarField& gammaNei,
+    const surfaceScalarField& aOwn,
+    const surfaceScalarField& aNei,
+    surfaceVectorField& U12,
+    surfaceScalarField& massFlux,
+    surfaceVectorField& momentumFlux,
+    surfaceScalarField& energyFlux
+)
+{
+
+    dimensionedScalar norP("norP", dimless, 1.0);
+    //dimensionedScalar minU("minU", dimVelocity, SMALL);
+
+    surfaceScalarField UvOwn((UOwn & normal) - meshPhi/magSf);
+
+    surfaceScalarField UvNei((UNei & normal) - meshPhi/magSf);
+
+    // Compute split velocity
+    surfaceScalarField alphaOwn(2*(pOwn/rhoOwn)/(pOwn/rhoOwn + pNei/rhoNei));
+    surfaceScalarField alphaNei(2 - alphaOwn);
+
+    surfaceScalarField cm(max(aOwn , aNei));
+
+    surfaceScalarField uPlus(
+      neg0(mag(UvOwn/cm) - 1)*(alphaOwn*(sqr(UvOwn + cm)/(4*cm) - 0.5*(UvOwn + mag(UvOwn))))
+       + 0.5*(UvOwn + mag(UvOwn))
+    );
+    surfaceScalarField uMinus(
+      neg0(mag(UvNei/cm) - 1)*(alphaNei*(-sqr(UvNei - cm)/(4*cm) - 0.5*(UvNei - mag(UvNei))))
+       + 0.5*(UvNei - mag(UvNei))
+    );
+
+    U12 = (alphaOwn*UOwn + alphaNei*UNei)/2.0;
+
+    surfaceScalarField pPlus(
+       neg0(mag(UvOwn/cm) - 1)*pOwn*sqr(UvOwn/cm + 1.0)*(2.0 - UvOwn/cm)/4.0
+       + Foam::pos(mag(UvOwn/cm) - 1)*pOwn*0.5*(1 + Foam::sign(UvOwn))//pos(mag(UvNei))
+    );
+    surfaceScalarField pMinus(
+       neg0(mag(UvNei/cm) - 1)*pNei*sqr(UvNei/cm - 1.0)*(2.0 + UvNei/cm)/4.0
+       + Foam::pos(mag(UvNei/cm) - 1)*pNei*0.5*(1 - Foam::sign(UvNei))//max(UvNei,minU)
+    );
+
+    surfaceScalarField P12(pPlus + pMinus);
+    surfaceScalarField s(0.5*min(norP , 10.0*mag(pNei - pOwn)/min(pOwn,pNei)));
+
+    surfaceScalarField caseA(Foam::neg(UvOwn - aOwn)*Foam::pos(UvNei - aNei));
+    surfaceScalarField caseB(Foam::neg(UvOwn + aOwn)*Foam::pos(UvNei + aNei));
+
+    massFlux = (uPlus*rhoOwn + uMinus*rhoNei)*magSf
+               -(1-caseA*caseB)*(caseA*0.125*(UvNei - aNei - UvOwn + aOwn)*(rhoNei - rhoOwn)*magSf
+               + (1-caseA)*caseB*0.125*(UvNei + aNei - UvOwn - aOwn)*(rhoNei - rhoOwn)*magSf);
+
+    surfaceVectorField AUSMV((uPlus*rhoOwn*UOwn + uMinus*rhoNei*UNei)*magSf);
+    surfaceVectorField AUSMD(0.5*(massFlux*(UOwn+UNei) - mag(massFlux)*(UNei-UOwn)));
+
+    momentumFlux = (0.5+s)*AUSMV + (0.5-s)*AUSMD + P12*normal*magSf
+                  -(1-caseA*caseB)*(caseA*0.125*(UvNei - aNei - UvOwn + aOwn)*(rhoNei*UNei - rhoOwn*UOwn)*magSf
+                  + (1-caseA)*caseB*0.125*(UvNei + aNei - UvOwn - aOwn)*(rhoNei*UNei - rhoOwn*UOwn)*magSf);
+
+    energyFlux = 0.5*(massFlux*(HOwn+HNei) - mag(massFlux)*(HNei-HOwn)) + meshPhi*P12
+                  -(1-caseA*caseB)*(caseA*0.125*(UvNei - aNei - UvOwn + aOwn)*(rhoNei*HNei - rhoOwn*HOwn)*magSf
+                  + (1-caseA)*caseB*0.125*(UvNei + aNei - UvOwn - aOwn)*(rhoNei*HNei - rhoOwn*HOwn)*magSf);
+
+
+}

applications/solvers/dfHighSpeedFoam/fluxSchemes/KNPFlux.H

@@ -0,0 +1,70 @@
+void KNPFlux(
+    const surfaceScalarField& meshPhi,
+    const surfaceScalarField& magSf,
+    const surfaceVectorField& normal,
+    const surfaceScalarField& rhoOwn,
+    const surfaceScalarField& rhoNei,
+    const surfaceVectorField& UOwn,
+    const surfaceVectorField& UNei,
+    const surfaceScalarField& HOwn,
+    const surfaceScalarField& HNei,
+    const surfaceScalarField& pOwn,
+    const surfaceScalarField& pNei,
+    const surfaceScalarField& gammaOwn,
+    const surfaceScalarField& gammaNei,
+    const surfaceScalarField& aOwn,
+    const surfaceScalarField& aNei,
+    surfaceVectorField& U12,
+    surfaceScalarField& massFlux,
+    surfaceVectorField& momentumFlux,
+    surfaceScalarField& energyFlux
+)
+{
+    surfaceScalarField UvOwn((UOwn & normal) - meshPhi/magSf);
+
+    surfaceScalarField UvNei((UNei & normal) - meshPhi/magSf);
+
+    surfaceScalarField EOwn("EOwn", HOwn - pOwn/rhoOwn);
+    surfaceScalarField ENei("ENei", HNei - pNei/rhoNei);
+
+    surfaceScalarField phiv_pos("phiv_pos", UvOwn*magSf);  //
+    surfaceScalarField phiv_neg("phiv_neg", UvNei*magSf);
+
+
+    surfaceScalarField cSf_pos("cSf_pos",aOwn*magSf);
+    surfaceScalarField cSf_neg("cSf_neg",aNei*magSf);
+    dimensionedScalar v_zero("v_zero", dimVolume/dimTime, 0.0);
+    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);
+
+    surfaceScalarField a_neg("a_neg", 1.0 - a_pos);
+
+    phiv_pos *= a_pos;
+    phiv_neg *= a_neg;
+
+    U12 = a_pos*UOwn + a_neg*UNei;
+
+    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));
+
+    massFlux = aphiv_pos*rhoOwn + aphiv_neg*rhoNei;
+
+    momentumFlux = (aphiv_pos*rhoOwn*UOwn + aphiv_neg*rhoNei*UNei)
+        + (a_pos*pOwn + a_neg*pNei)*magSf*normal;
+
+    energyFlux = aphiv_pos*(rhoOwn*EOwn + pOwn)
+        + aphiv_neg*(rhoNei*ENei + pNei)
+        + aSf*pOwn - aSf*pNei
+        + meshPhi*(a_pos*pOwn + a_neg*pNei);
+
+}

applications/solvers/dfHighSpeedFoam/phiCal.H

@@ -1,4 +1,82 @@
-phi = aphiv_pos*rho_pos + aphiv_neg*rho_neg;
+//phi = aphiv_pos*rho_pos + aphiv_neg*rho_neg;
+
+surfaceVectorField  phiUp
+(
+    IOobject
+    (
+        "phiUp",
+        mesh.time().timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedVector("0", sqr(dimVelocity)*dimDensity*dimArea, Zero)
+);
+
+surfaceVectorField  u12
+(
+    IOobject
+    (
+        "u12",
+        mesh.time().timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedVector("0", dimVelocity, Zero)
+);
+
+surfaceScalarField  phiEp
+(
+    IOobject
+    (
+        "phiEp",
+        mesh.time().timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("0", pow3(dimVelocity)*dimDensity*dimArea, 0)
+);
+
+
+surfaceScalarField  meshPhi
+(
+    IOobject
+    (
+        "meshPhi",
+        mesh.time().timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("0", (dimVelocity)*dimArea, 0)
+);
+if(mesh.moving())
+{
+  meshPhi = meshPhi + mesh.phi();
+}
+
+//-- Computing flux
+fluxMap[fluxScheme](
+    meshPhi,
+    mesh.magSf(),
+    normal,
+    rho_pos,rho_neg,
+    U_pos,U_neg,
+    H_pos,H_neg,
+    p_pos,p_neg,
+    gamma_pos,gamma_neg,
+    a_pos,a_neg,
+    u12,
+    phi,
+    phiUp,
+    phiEp
+);
 
 PtrList<surfaceScalarField> phiYi(nspecies);
 forAll(phiYi,i)
@@ -16,30 +94,12 @@ forAll(phiYi,i)
         IOobject::NO_READ,
         IOobject::NO_WRITE
       ),
-      aphiv_pos*rhoYi_pos[i] + aphiv_neg*rhoYi_neg[i]
+      phi*fvc::interpolate(Y[i], "Yi").ref() 
     )
   );
 }
 
-surfaceVectorField phiUp
-(
-  (aphiv_pos*rhoU_pos + aphiv_neg*rhoU_neg)
-    + (a_pos*p_pos + a_neg*p_neg)*mesh.Sf()
-);
-
-surfaceScalarField phiEp
-(
-    "phiEp",
-    aphiv_pos*(rho_pos*(ea_pos + 0.5*magSqr(U_pos)) + p_pos)
-  + aphiv_neg*(rho_neg*(ea_neg + 0.5*magSqr(U_neg)) + p_neg)
-  + aSf*p_pos - aSf*p_neg
-);
 
-// Make flux for pressure-work absolute
-if (mesh.moving())
-{
-  phiEp += mesh.phi()*(a_pos*p_pos + a_neg*p_neg);
-}
 
 volScalarField muEff("muEff", turbulence->muEff());
 volTensorField tauMC("tauMC", muEff*dev2(Foam::T(fvc::grad(U))));