From 77201f0b693f71b99f9ae71de394b386a157ff21 Mon Sep 17 00:00:00 2001
From: hanli94 <han_li@pku.edu.cn>
Date: Wed, 19 Jul 2023 15:07:40 +0800
Subject: [PATCH] Use diffusion correction only when turbulence type is laminar
 for dfHighSpeedFoam

---
 .../solvers/dfHighSpeedFoam/createFields.H    |  1 +
 .../solvers/dfHighSpeedFoam/rhoEEqn.H         | 38 +++++++++++++++----
 .../solvers/dfHighSpeedFoam/rhoYEqn.H         | 17 +++++++--
 3 files changed, 44 insertions(+), 12 deletions(-)

diff --git a/applications/solvers/dfHighSpeedFoam/createFields.H b/applications/solvers/dfHighSpeedFoam/createFields.H
index ee65208b0..852122f45 100644
--- a/applications/solvers/dfHighSpeedFoam/createFields.H
+++ b/applications/solvers/dfHighSpeedFoam/createFields.H
@@ -153,6 +153,7 @@ autoPtr<compressible::turbulenceModel> turbulence
         thermo
     )
 );
+const word turbName(mesh.objectRegistry::lookupObject<IOdictionary>("turbulenceProperties").lookup("simulationType"));
 
 multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 
diff --git a/applications/solvers/dfHighSpeedFoam/rhoEEqn.H b/applications/solvers/dfHighSpeedFoam/rhoEEqn.H
index 2b96e8862..3a76804d4 100644
--- a/applications/solvers/dfHighSpeedFoam/rhoEEqn.H
+++ b/applications/solvers/dfHighSpeedFoam/rhoEEqn.H
@@ -38,9 +38,21 @@ if (!inviscid)
 {
     if((ddtSchemes == "RK2SSP") || (ddtSchemes == "RK3SSP"))
     {
-        rhoE_rhs =   fvc::laplacian(turbulence->alphaEff()*thermo.gamma(), ea)
+        rhoE_rhs =   
+            (
+                turbName == "laminar"
+                ?
+                (
+                    fvc::laplacian(turbulence->alphaEff()*thermo.gamma(), ea)
                   - diffAlphaD 
-                  + fvc::div(hDiffCorrFlux);
+                  + fvc::div(hDiffCorrFlux)
+                )
+                :
+                (
+                    fvc::laplacian(turbulence->alphaEff()*thermo.gamma(), ea)
+                )
+            );
+                  
         
         rhoE +=  rkcoe3[nrk]*rhoE_rhs*runTime.deltaT();
 
@@ -54,12 +66,22 @@ if (!inviscid)
     {
         fvScalarMatrix eEqn
         (
-              fvm::ddt(rho, ea) - fvc::ddt(rho, ea)
-            // alpha in deepflame is considered to calculate h by default (kappa/Cp), so multiply gamma to correct alpha
-            - fvm::laplacian(turbulence->alphaEff()*thermo.gamma(), ea)
-            + diffAlphaD
-          ==
-              fvc::div(hDiffCorrFlux)
+                fvm::ddt(rho, ea) - fvc::ddt(rho, ea)
+            ==
+                (
+                    turbName == "laminar"
+                    ?
+                    (
+                        // alpha in deepflame is considered to calculate h by default (kappa/Cp), so multiply gamma to correct alpha
+                        fvm::laplacian(turbulence->alphaEff()*thermo.gamma(), ea)
+                    -   diffAlphaD
+                    +   fvc::div(hDiffCorrFlux)
+                    )
+                    :
+                    (
+                        fvm::laplacian(turbulence->alphaEff()*thermo.gamma(), ea)
+                    )
+                )
         );
 
         eEqn.solve("ea");
diff --git a/applications/solvers/dfHighSpeedFoam/rhoYEqn.H b/applications/solvers/dfHighSpeedFoam/rhoYEqn.H
index de80a3d54..fac11306c 100644
--- a/applications/solvers/dfHighSpeedFoam/rhoYEqn.H
+++ b/applications/solvers/dfHighSpeedFoam/rhoYEqn.H
@@ -107,7 +107,12 @@ time_monitor_corrDiff += double(end - start) / double(CLOCKS_PER_SEC);
 
                 if((ddtSchemes == "RK2SSP") || (ddtSchemes == "RK3SSP"))
                 {
-                    rhoYi_rhs = fvc::laplacian(DEff(), Yi) - fvc::div(phiUc,Yi,"div(phic,Yi)");
+                    rhoYi_rhs = 
+                        (
+                            turbName == "laminar"
+                            ?  (fvc::laplacian(DEff(), Yi) - fvc::div(phiUc,Yi,"div(phic,Yi)")) 
+                            :  fvc::laplacian(DEff(), Yi)
+                        );
 
                     rhoYi[i] += rkcoe3[nrk]*rhoYi_rhs*runTime.deltaT();
 
@@ -121,9 +126,13 @@ time_monitor_corrDiff += double(end - start) / double(CLOCKS_PER_SEC);
 
                     fvScalarMatrix YiEqn
                     (
-                        fvm::ddt(rho, Yi) - fvc::ddt(rho, Yi)
-                    - fvm::laplacian(DEff(), Yi)
-                    + mvConvection->fvmDiv(phiUc, Yi)
+                            fvm::ddt(rho, Yi) - fvc::ddt(rho, Yi)
+                        -
+                            (
+                                turbName == "laminar"
+                                ?  (fvm::laplacian(DEff(), Yi) + mvConvection->fvmDiv(phiUc, Yi))
+                                :  fvm::laplacian(DEff(), Yi)
+                            )
                     );
 
                     YiEqn.relax();