Skip to content

Add in Sutherland's law for NEMO problems #1702

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
WallyMaier merged 7 commits into from
Jul 15, 2022
Merged

Add in Sutherland's law for NEMO problems #1702

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
3205992
Add in Sutherland's law for NEMO problems + some clean up to viscous …
Jul 7, 2022
f166ce3
Merge branch 'develop' into feature_NEMO_suth
WallyMaier Jul 10, 2022
6c3eff1
PR comments
Jul 10, 2022
a54a4bb
Merge branch 'feature_NEMO_suth' of https://github.com/su2code/SU2 in…
Jul 10, 2022
cfffec7
Merge branch 'develop' into feature_NEMO_suth
WallyMaier Jul 14, 2022
6cdeb66
Merge branch 'develop' into feature_NEMO_suth
WallyMaier Jul 15, 2022
15ae0bf
Apply suggestions from code review
WallyMaier Jul 15, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions Common/include/option_structure.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -613,11 +613,13 @@ MakePair("ONESPECIES", ONESPECIES)
* \brief types of coefficient transport model
*/
enum class TRANSCOEFFMODEL {
SUTHERLAND,
WILKE,
GUPTAYOS,
CHAPMANN_ENSKOG
};
static const MapType<std::string, TRANSCOEFFMODEL> TransCoeffModel_Map = {
MakePair("SUTHERLAND", TRANSCOEFFMODEL::SUTHERLAND)
MakePair("WILKE", TRANSCOEFFMODEL::WILKE)
MakePair("GUPTA-YOS", TRANSCOEFFMODEL::GUPTAYOS)
MakePair("CHAPMANN-ENSKOG", TRANSCOEFFMODEL::CHAPMANN_ENSKOG)
Expand Down
4 changes: 2 additions & 2 deletions Common/src/CConfig.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3792,8 +3792,8 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
SU2_MPI::Error("Only STANDARD_AIR fluid model can be used with US Measurement System", CURRENT_FUNCTION);
}

if (Kind_FluidModel == SU2_NONEQ && Kind_TransCoeffModel != TRANSCOEFFMODEL::WILKE ) {
SU2_MPI::Error("Only WILKE transport model is stable for the NEMO solver using SU2TClib. Use Mutation++ instead.", CURRENT_FUNCTION);
if (Kind_FluidModel == SU2_NONEQ && (Kind_TransCoeffModel != TRANSCOEFFMODEL::WILKE && Kind_TransCoeffModel != TRANSCOEFFMODEL::SUTHERLAND) ) {
SU2_MPI::Error("Only WILKE and SUTHERLAND transport models are stable for the NEMO solver using SU2TClib. Use Mutation++ instead.", CURRENT_FUNCTION);

@jtneedels jtneedels Jul 11, 2022

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How hard would it be to add this routine so it is supported by MPP?

Also, just to confirm, Sutherland is just used for viscosity then WBE is used for diffusion coeffs and thermal conductivity?

@WallyMaier WallyMaier Jul 11, 2022

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Great question. I don't know if MPP has sutherland's law, but if it does....it should be fairly simple to implement.
I guess with that thought, the Sutherlands law would be the identical implementation, but it would need additional errors if someone tried to used not AR, N2, or Air5/7.

Yes, Sutherland is used for viscosity and thermal conductivity. Blottner part of WBE is still used for diffusion coefficients.

}

if (Kind_FluidModel == MUTATIONPP && (Kind_TransCoeffModel != TRANSCOEFFMODEL::WILKE && Kind_TransCoeffModel != TRANSCOEFFMODEL::CHAPMANN_ENSKOG)) {
Expand Down
29 changes: 23 additions & 6 deletions SU2_CFD/include/fluid/CSU2TCLib.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -64,6 +64,13 @@ class CSU2TCLib : public CNEMOGas {
phis, mus, /*!< \brief Auxiliary vectors to be used in Wilke/Blottner/Eucken model */
A; /*!< \brief Auxiliary vector to be used in net production rate computation */

std::array<su2double,1> mu_ref; /*!< \brief Vector containing reference viscosity for Sutherland's law */
std::array<su2double,1> k_ref; /*!< \brief Vector containing reference thermal conducivities for Sutherland's law */
std::array<su2double,1> Sm_ref; /*!< \brief Vector containing Sutherland's constant for viscosity */
std::array<su2double,1> Sk_ref; /*!< \brief Vector containing Sutherland's constant for thermal conductivities */

const su2double T_ref_suth = 273.15; /*!<\brief Reference temperature for Sutherland's model [K] */

su2activematrix CharElTemp, /*!< \brief Characteristic temperature of electron states. */
ElDegeneracy, /*!< \brief Degeneracy of electron states. */
RxnConstantTable, /*!< \brief Table of chemical equiibrium reaction constants */
Expand Down Expand Up @@ -203,35 +210,45 @@ class CSU2TCLib : public CNEMOGas {
void ComputeKeqConstants(unsigned short val_Reaction);

/*!
* \brief Get species diffusion coefficients with Wilke/Blottner/Eucken transport model.
* \brief Calculate species diffusion coefficients with Wilke/Blottner/Eucken transport model.
*/
void DiffusionCoeffWBE();

/*!
* \brief Get viscosity with Wilke/Blottner/Eucken transport model.
* \brief Calculate viscosity with Wilke/Blottner/Eucken transport model.
*/
void ViscosityWBE();

/*!
* \brief Get T-R and V-E thermal conductivities vector with Wilke/Blottner/Eucken transport model.
* \brief Calculate T-R and V-E thermal conductivities vector with Wilke/Blottner/Eucken transport model.
*/
void ThermalConductivitiesWBE();

/*!
* \brief Get species diffusion coefficients with Gupta-Yos transport model.
* \brief Calculate species diffusion coefficients with Gupta-Yos transport model.
*/
void DiffusionCoeffGY();

/*!
* \brief Get viscosity with Gupta-Yos transport model.
* \brief Calculate viscosity with Gupta-Yos transport model.
*/
void ViscosityGY();

/*!
* \brief Get T-R and V-E thermal conductivities vector with Gupta-Yos transport model.
* \brief Calculate T-R and V-E thermal conductivities vector with Gupta-Yos transport model.
*/
void ThermalConductivitiesGY();

/*!
* \brief Calculate viscosity with Sutherland's transport model.
*/
void ViscositySuth();

/*!
* \brief Calculate T-R and V-E thermal conductivities vector with Sutherland's transport model.
*/
void ThermalConductivitiesSuth();

/*!
* \brief Get reference temperature.
*/
Expand Down
80 changes: 74 additions & 6 deletions SU2_CFD/src/fluid/CSU2TCLib.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -52,7 +52,7 @@ CSU2TCLib::CSU2TCLib(const CConfig* config, unsigned short val_nDim, bool viscou
eve_eq.resize(nSpecies,0.0);
eve.resize(nSpecies,0.0);

if(viscous){
if (viscous) {
MolarFracWBE.resize(nSpecies,0.0);
phis.resize(nSpecies,0.0);
mus.resize(nSpecies,0.0);
Expand Down Expand Up @@ -110,6 +110,14 @@ CSU2TCLib::CSU2TCLib(const CConfig* config, unsigned short val_nDim, bool viscou
ElDegeneracy(0,5) = 5;
ElDegeneracy(0,6) = 15;

if (viscous) {
//F.M. White, Viscous Fluid Flow, 3rd ed., McGraw-Hill, 2006.
mu_ref[0] = 2.125E-5;

@jtneedels jtneedels Jul 11, 2022

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Worth adding a comment with a reference for these values?

@WallyMaier WallyMaier Jul 11, 2022

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good point. I can add a quick reference.

k_ref[0] = 0.0163;
Sm_ref[0] = 114.0;
Sk_ref[0] = 170;
}

} else if (gas_model == "N2"){
/*--- Check for errors in the initialization ---*/
if (nSpecies != 2) {
Expand Down Expand Up @@ -252,6 +260,14 @@ CSU2TCLib::CSU2TCLib(const CConfig* config, unsigned short val_nDim, bool viscou
Omega11(1,0,0) = -8.3493693E-03; Omega11(1,0,1) = 1.7808911E-01; Omega11(1,0,2) = -1.4466155E+00; Omega11(1,0,3) = 1.9324210E+03;
Omega11(1,1,0) = -7.7439615E-03; Omega11(1,1,1) = 1.7129007E-01; Omega11(1,1,2) = -1.4809088E+00; Omega11(1,1,3) = 2.1284951E+03;

if (viscous) {
//F.M. White, Viscous Fluid Flow, 3rd ed., McGraw-Hill, 2006.
k_ref[0] = 0.0242;
mu_ref[0] = 1.663E-5;
Sm_ref[0] = 107.0;
Sk_ref[0] = 150.0;
}

} else if (gas_model == "AIR-5"){

/*--- Check for errors in the initialization ---*/
Expand Down Expand Up @@ -598,6 +614,14 @@ CSU2TCLib::CSU2TCLib(const CConfig* config, unsigned short val_nDim, bool viscou
Omega11(4,3,0) = -5.0478143E-03; Omega11(4,3,1) = 1.0236186E-01; Omega11(4,3,2) = -9.0058935E-01; Omega11(4,3,3) = 4.4472565E+02;
Omega11(4,4,0) = -4.2451096E-03; Omega11(4,4,1) = 9.6820337E-02; Omega11(4,4,2) = -9.9770795E-01; Omega11(4,4,3) = 8.3320644E+02;

if (viscous) {
//F.M. White, Viscous Fluid Flow, 3rd ed., McGraw-Hill, 2006.
k_ref[0] = 0.0241;
mu_ref[0] = 1.716E-5;
Sm_ref[0] = 111.0;
Sk_ref[0] = 194.0;
}

} else if (gas_model == "AIR-7"){

/*--- Check for errors in the initialization ---*/
Expand Down Expand Up @@ -1027,6 +1051,14 @@ CSU2TCLib::CSU2TCLib(const CConfig* config, unsigned short val_nDim, bool viscou
Omega11(4,2,0) = -1.0066279E-03; Omega11(4,2,1) = 1.1029264E-02; Omega11(4,2,2) = -2.0671266E-01; Omega11(4,2,3) = 8.2644384E+01;
Omega11(4,3,0) = -5.0478143E-03; Omega11(4,3,1) = 1.0236186E-01; Omega11(4,3,2) = -9.0058935E-01; Omega11(4,3,3) = 4.4472565E+02;
Omega11(4,4,0) = -4.2451096E-03; Omega11(4,4,1) = 9.6820337E-02; Omega11(4,4,2) = -9.9770795E-01; Omega11(4,4,3) = 8.3320644E+02;

if (viscous) {
//F.M. White, Viscous Fluid Flow, 3rd ed., McGraw-Hill, 2006.
k_ref[0] = 0.0241;
mu_ref[0] = 1.716E-5;
Sm_ref[0] = 111.0;
Sk_ref[0] = 194.0;
}
}

if (ionization) { nHeavy = nSpecies-1; nEl = 1; }
Expand Down Expand Up @@ -1598,6 +1630,8 @@ vector<su2double>& CSU2TCLib::GetDiffusionCoeff(){
DiffusionCoeffWBE();
if(Kind_TransCoeffModel == TRANSCOEFFMODEL::GUPTAYOS)
DiffusionCoeffGY();
if(Kind_TransCoeffModel == TRANSCOEFFMODEL::SUTHERLAND)
DiffusionCoeffWBE();

@WallyMaier WallyMaier Jul 7, 2022

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The diffusion coefficient used in the Sutherland model is still computed using the WBE routine


return DiffusionCoeff;

Expand All @@ -1609,6 +1643,8 @@ su2double CSU2TCLib::GetViscosity(){
ViscosityWBE();
if(Kind_TransCoeffModel == TRANSCOEFFMODEL::GUPTAYOS)
ViscosityGY();
if(Kind_TransCoeffModel == TRANSCOEFFMODEL::SUTHERLAND)
ViscositySuth();

return Mu;

Expand All @@ -1620,6 +1656,8 @@ vector<su2double>& CSU2TCLib::GetThermalConductivities(){
ThermalConductivitiesWBE();
if(Kind_TransCoeffModel == TRANSCOEFFMODEL::GUPTAYOS)
ThermalConductivitiesGY();
if(Kind_TransCoeffModel == TRANSCOEFFMODEL::SUTHERLAND)
ThermalConductivitiesSuth();

return ThermalConductivities;

Expand Down Expand Up @@ -1802,7 +1840,7 @@ void CSU2TCLib::DiffusionCoeffGY(){
//}

/*--- Assign species diffusion coefficient ---*/
DiffusionCoeff[iSpecies] = gam_t*gam_t*Mi*(1-Mi*gam_i) / denom;
DiffusionCoeff[iSpecies] = (denom > EPS) ? (gam_t*gam_t*Mi*(1-Mi*gam_i) / denom) : su2double(0.0);
}
// if (ionization) {
//TODO: Update correct iElectron....
Expand Down Expand Up @@ -1971,18 +2009,48 @@ void CSU2TCLib::ThermalConductivitiesGY(){
}

/*--- Translational contribution to thermal conductivity ---*/
ThermalCond_tr += (15.0/4.0)*kb*gam_i/denom_t;
ThermalCond_tr += (denom_t > EPS) ? ((15.0/4.0)*kb*gam_i/denom_t) : su2double(0.0);

/*--- Translational contribution to thermal conductivity ---*/
if (RotationModes[iSpecies] != 0.0)
ThermalCond_tr += kb*gam_i/denom_r;
if (RotationModes[iSpecies] != 0.0) ThermalCond_tr += (denom_r > EPS) ? (kb*gam_i/denom_r) : su2double(0.0);

/*--- Vibrational-electronic contribution to thermal conductivity ---*/
ThermalCond_ve += kb*Cvve/R*gam_i / denom_r;
ThermalCond_ve += (denom_r > EPS) ? (kb*Cvve/R*gam_i / denom_r) : su2double(0.0);
}

ThermalConductivities[0] = ThermalCond_tr;
ThermalConductivities[1] = ThermalCond_ve;
}

void CSU2TCLib::ViscositySuth(){

su2double T_nd = T / T_ref_suth;

/*--- Calculate mixture laminar viscosity ---*/
Mu = mu_ref[0] * T_nd * sqrt(T_nd) * ((T_ref_suth + Sm_ref[0]) / (T + Sm_ref[0]));

}

void CSU2TCLib::ThermalConductivitiesSuth(){

/*--- Compute mixture quantities ---*/
su2double mass = 0.0, rho = 0.0;
for (unsigned short ii=0; ii<nSpecies; ii++) rho += rhos[ii];
for (unsigned short ii=0; ii<nSpecies; ii++) mass += rhos[ii]/rho*MolarMass[ii];

su2double Cvtr = ComputerhoCvtr()/rho;
su2double Cvve = ComputerhoCvve()/rho;

/*--- Compute simple Kve scaling factor ---*/
su2double scl = Cvve/Cvtr;

/*--- Compute k's using Sutherland's law ---*/
su2double T_nd = T / T_ref_suth;
su2double k = k_ref[0] * T_nd * sqrt(T_nd) * ((T_ref_suth + Sk_ref[0]) / (T + Sk_ref[0]));
su2double kve = scl*k;

ThermalConductivities[0] = k;
ThermalConductivities[1] = kve;

}

Expand Down
70 changes: 32 additions & 38 deletions SU2_CFD/src/numerics/NEMO/CNEMONumerics.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -228,61 +228,58 @@ void CNEMONumerics::GetViscousProjFlux(const su2double *val_primvar,
su2double val_therm_conductivity_ve,
const CConfig *config) {

if (ionization) {
SU2_MPI::Error("NEED TO IMPLEMENT IONIZED FUNCTIONALITY!!!",CURRENT_FUNCTION);
}

// Requires a slightly non-standard primitive vector:
// Assumes - V = [Y1, ... , Yn, T, Tve, ... ]
// and gradient GV = [GY1, ... , GYn, GT, GTve, ... ]
// rather than the standard V = [r1, ... , rn, T, Tve, ... ]

unsigned short iSpecies, iVar, iDim, jDim;
su2double mu, ktr, kve, rho, T, Tve, RuSI, Ru;
auto& Ms = fluidmodel->GetSpeciesMolarMass();

su2activematrix Flux_Tensor(nVar,nDim);

/*--- Initialize ---*/
for (iVar = 0; iVar < nVar; iVar++) {
for (auto iVar = 0; iVar < nVar; iVar++) {
Proj_Flux_Tensor[iVar] = 0.0;
for (iDim = 0; iDim < nDim; iDim++)
for (auto iDim = 0; iDim < nDim; iDim++)
Flux_Tensor[iVar][iDim] = 0.0;
}

/*--- Rename for convenience ---*/
/*--- Rename variables for convenience ---*/
const auto& Ms = fluidmodel->GetSpeciesMolarMass();
const auto& Ds = val_diffusioncoeff;
mu = val_lam_viscosity+val_eddy_viscosity;
ktr = val_therm_conductivity;
kve = val_therm_conductivity_ve;
rho = val_primvar[RHO_INDEX];
T = val_primvar[T_INDEX];
Tve = val_primvar[TVE_INDEX];
const su2double mu = val_lam_viscosity+val_eddy_viscosity;
su2double ktr = val_therm_conductivity;
su2double kve = val_therm_conductivity_ve;
const su2double rho = val_primvar[RHO_INDEX];
const su2double T = val_primvar[T_INDEX];
Comment on lines +255 to +256

@pcarruscag pcarruscag Jul 7, 2022

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

👍 nice

const su2double Tve = val_primvar[TVE_INDEX];
const auto& V = val_primvar;
const auto& GV = val_gradprimvar;
RuSI= UNIVERSAL_GAS_CONSTANT;
Ru = 1000.0*RuSI;

const su2double Ru = 1000.0*UNIVERSAL_GAS_CONSTANT;
const auto& hs = fluidmodel->ComputeSpeciesEnthalpy(T, Tve, val_eve);

/*--- Scale thermal conductivity with turb visc ---*/
// TODO: Need to determine proper way to incorporate eddy viscosity
// This is only scaling Kve by same factor as ktr
// NOTE: V[iSpecies] is == Ys.
su2double Mass = 0.0;
su2double tmp1, scl, Cptr;
for (iSpecies=0;iSpecies<nSpecies;iSpecies++)
for (auto iSpecies = 0;iSpecies<nSpecies;iSpecies++)
Mass += V[iSpecies]*Ms[iSpecies];
Cptr = V[RHOCVTR_INDEX]/V[RHO_INDEX]+Ru/Mass;
tmp1 = Cptr*(val_eddy_viscosity/Prandtl_Turb);
scl = tmp1/ktr;

su2double Cptr = V[RHOCVTR_INDEX]/V[RHO_INDEX]+Ru/Mass;
su2double tmp1 = Cptr*(val_eddy_viscosity/Prandtl_Turb);
su2double scl = tmp1/ktr;
ktr += Cptr*(val_eddy_viscosity/Prandtl_Turb);
kve = kve*(1.0+scl);
//Cpve = V[RHOCVVE_INDEX]+Ru/Mass;
//kve += Cpve*(val_eddy_viscosity/Prandtl_Turb);

/*--- Pre-compute mixture quantities ---*/

/*--- Pre-compute mixture quantities ---*/ //TODO
su2double Vector[MAXNDIM] = {0.0};

for (iDim = 0; iDim < nDim; iDim++) {
for (iSpecies = 0; iSpecies < nHeavy; iSpecies++) {
for (auto iDim = 0; iDim < nDim; iDim++) {
for (auto iSpecies = 0; iSpecies < nHeavy; iSpecies++) {
Vector[iDim] += rho*Ds[iSpecies]*GV[RHOS_INDEX+iSpecies][iDim];
}
}
Expand All @@ -291,37 +288,34 @@ void CNEMONumerics::GetViscousProjFlux(const su2double *val_primvar,
ComputeStressTensor(nDim,tau,val_gradprimvar+VEL_INDEX, mu);

/*--- Populate entries in the viscous flux vector ---*/
for (iDim = 0; iDim < nDim; iDim++) {
for (auto iDim = 0; iDim < nDim; iDim++) {

/*--- Species diffusion velocity ---*/
for (iSpecies = 0; iSpecies < nHeavy; iSpecies++) {
for (auto iSpecies = 0; iSpecies < nHeavy; iSpecies++) {
Flux_Tensor[iSpecies][iDim] = rho*Ds[iSpecies]*GV[RHOS_INDEX+iSpecies][iDim]
- V[RHOS_INDEX+iSpecies]*Vector[iDim];
}
if (ionization) {
SU2_MPI::Error("NEED TO IMPLEMENT IONIZED FUNCTIONALITY!!!",CURRENT_FUNCTION);
}

/*--- Shear stress related terms ---*/
/*--- Shear-stress/momentum related terms ---*/
Flux_Tensor[nSpecies+nDim][iDim] = 0.0;
for (jDim = 0; jDim < nDim; jDim++) {
for (auto jDim = 0; jDim < nDim; jDim++) {
Flux_Tensor[nSpecies+jDim][iDim] = tau[iDim][jDim];
Flux_Tensor[nSpecies+nDim][iDim] += tau[iDim][jDim]*val_primvar[VEL_INDEX+jDim];
}

/*--- Diffusion terms ---*/
for (iSpecies = 0; iSpecies < nHeavy; iSpecies++) {
for (auto iSpecies = 0; iSpecies < nHeavy; iSpecies++) {
Flux_Tensor[nSpecies+nDim][iDim] += Flux_Tensor[iSpecies][iDim] * hs[iSpecies];
Flux_Tensor[nSpecies+nDim+1][iDim] += Flux_Tensor[iSpecies][iDim] * val_eve[iSpecies];
}

/*--- Heat transfer terms ---*/
Flux_Tensor[nSpecies+nDim][iDim] += ktr*GV[T_INDEX][iDim] +
kve*GV[TVE_INDEX][iDim];
Flux_Tensor[nSpecies+nDim][iDim] += ktr*GV[T_INDEX][iDim] + kve*GV[TVE_INDEX][iDim];
Flux_Tensor[nSpecies+nDim+1][iDim] += kve*GV[TVE_INDEX][iDim];
}

for (iVar = 0; iVar < nVar; iVar++) {
for (iDim = 0; iDim < nDim; iDim++) {
for (auto iVar = 0; iVar < nVar; iVar++) {
for (auto iDim = 0; iDim < nDim; iDim++) {
Proj_Flux_Tensor[iVar] += Flux_Tensor[iVar][iDim]*val_normal[iDim];
}
}
Expand Down
12 changes: 4 additions & 8 deletions SU2_CFD/src/numerics/NEMO/NEMO_diffusion.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -128,14 +128,10 @@ CNumerics::ResidualType<> CAvgGrad_NEMO::ComputeResidual(const CConfig *config)
PrimVar_j[iSpecies] = V_j[iSpecies]/V_j[RHO_INDEX];
Mean_PrimVar[iSpecies] = 0.5*(PrimVar_i[iSpecies] + PrimVar_j[iSpecies]);
for (auto iDim = 0; iDim < nDim; iDim++) {
Mean_GradPrimVar[iSpecies][iDim] = 0.5*(1.0/V_i[RHO_INDEX] *
(PrimVar_Grad_i[iSpecies][iDim] -
PrimVar_i[iSpecies] *
PrimVar_Grad_i[RHO_INDEX][iDim]) +
1.0/V_j[RHO_INDEX] *
(PrimVar_Grad_j[iSpecies][iDim] -
PrimVar_j[iSpecies] *
PrimVar_Grad_j[RHO_INDEX][iDim]));
Mean_GradPrimVar[iSpecies][iDim] = 0.5*(1.0/V_i[RHO_INDEX] * (PrimVar_Grad_i[iSpecies][iDim] -
PrimVar_i[iSpecies] * PrimVar_Grad_i[RHO_INDEX][iDim]) +
1.0/V_j[RHO_INDEX] * (PrimVar_Grad_j[iSpecies][iDim] -
PrimVar_j[iSpecies] * PrimVar_Grad_j[RHO_INDEX][iDim]));
}
}

Expand Down
Loading