Skip to content

Update 1T model and catalytic boundary condition handling in NEMO #13

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
raghava-davuluri merged 1 commit into from
Feb 5, 2026
Merged

Update 1T model and catalytic boundary condition handling in NEMO #13

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
36 changes: 19 additions & 17 deletions SU2_CFD/src/fluid/CMutationTCLib.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
* \file CMutationTCLib.cpp
* \brief Source of the Mutation++ 2T nonequilibrium gas model.
* \author C. Garbacz
* \version 8.3.0 "Harrier"
* \version 8.2.0 "Harrier"
*
* SU2 Project Website: https://su2code.github.io
*
Expand Down Expand Up @@ -56,16 +56,13 @@ CMutationTCLib::CMutationTCLib(const CConfig* config, unsigned short val_nDim):
else if (Kind_TransCoeffModel == TRANSCOEFFMODEL::CHAPMANN_ENSKOG)
transport_model = "Chapmann-Enskog_LDLT";

NEWTON_ROBUST = config->Get_Mpp_Temp_Solver_Robust();
NEWTON = config->Get2TNewton();
if (NoneqStateModel == "2T") {
opt.setStateModel("ChemNonEqTTv");
}
else if (NoneqStateModel == "1T"){
opt.setStateModel("ChemNonEq1T");
}
else if (NoneqStateModel == "EQUILIBRIUM"){
opt.setStateModel("Equil");
}

if (frozen) opt.setMechanism("none");

Expand Down Expand Up @@ -141,7 +138,7 @@ void CMutationTCLib::SetTDStateRhosTTv(vector<su2double>& val_rhos, su2double va

Pressure = ComputePressure();

mix->setState(rhos.data(), temperatures.data(), 1, NEWTON_ROBUST);
mix->setState(rhos.data(), temperatures.data(), 1, NEWTON);

}

Expand Down Expand Up @@ -214,24 +211,27 @@ void CMutationTCLib::ChemistryJacobian(unsigned short iReaction, const su2double
const su2double* dTdU, const su2double* dTvedU, su2double **val_jacobian){

unsigned short iVar, jVar, iSpecies;
unsigned short nEve = nSpecies+nDim+1;
unsigned short nVar = nSpecies+nDim+2;

mix->jacobianRho(JacRho.data());
mix->jacobianT(JacT.data());
mix->jacobianTv(JacTv.data());

for(iSpecies = 0; iSpecies < nSpecies; iSpecies++)
for(jSpecies = 0; jSpecies < nSpecies; jSpecies++)
val_jacobian[iSpecies][jSpecies] = JacRho[iSpecies*nSpecies+jSpecies];

for(iSpecies = 0; iSpecies < nSpecies; iSpecies++){
for(jSpecies = 0; jSpecies < nSpecies; jSpecies++){
val_jacobian[iSpecies][jSpecies] = JacRho[iSpecies*nSpecies+jSpecies];}
for(iVar = 0; iVar < nVar; iVar++){
val_jacobian[iSpecies][iVar] += JacT[iSpecies]*dTdU[iVar];
val_jacobian[iSpecies][iVar] += JacTv[iSpecies]*dTvedU[iVar];
}
}
val_jacobian[iSpecies][iVar] += JacT[iSpecies]*dTdU[iVar];}
}


if (NoneqStateModel == "2T") {
mix->jacobianTv(JacTv.data());
for(iSpecies = 0; iSpecies < nSpecies; iSpecies++){
for(iVar = 0; iVar < nVar; iVar++){
val_jacobian[iSpecies][iVar] += JacTv[iSpecies]*dTvedU[iVar];}
}
}

}

su2double CMutationTCLib::ComputeEveSourceTerm(){
Expand All @@ -246,6 +246,8 @@ su2double CMutationTCLib::ComputeEveSourceTerm(){
void CMutationTCLib::GetEveSourceTermJacobian(const su2double *V, const su2double *eve, const su2double *cvve, const su2double *dTdU,
const su2double* dTvedU, su2double **val_jacobian){

if (NoneqStateModel == "1T") return;

unsigned short iVar, jVar, iSpecies;
unsigned short nEve = nSpecies+nDim+1;
unsigned short nVar = nSpecies+nDim+2;
Expand Down Expand Up @@ -299,7 +301,7 @@ vector<su2double>& CMutationTCLib::ComputeTemperatures(vector<su2double>& val_rh
energies[0] = rhoE - rhoEvel;
energies[1] = rhoEve;

mix->setState(rhos.data(), energies.data(), 0, NEWTON_ROBUST);
mix->setState(rhos.data(), energies.data(), 0, NEWTON);

mix->getTemperatures(temperatures.data());

Expand Down
18 changes: 12 additions & 6 deletions SU2_CFD/src/fluid/CNEMOGas.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
* \file CNEMOGas.cpp
* \brief Source of the nonequilibrium gas model.
* \author C. Garbacz, W. Maier, S. R. Copeland
* \version 8.3.0 "Harrier"
* \version 8.2.0 "Harrier"
*
* SU2 Project Website: https://su2code.github.io
*
Expand Down Expand Up @@ -232,8 +232,8 @@ void CNEMOGas::ComputedPdU(const su2double *V, const vector<su2double>& val_eves
val_dPdU[nSpecies+nDim] = conc*Ru / rhoCvtr;

/*--- Vib.-el energy derivative ---*/
val_dPdU[nSpecies+nDim+1] = -val_dPdU[nSpecies+nDim] +
rho_el*Ru/MolarMass[0]*1.0/rhoCvve;
val_dPdU[nSpecies+nDim+1] = -val_dPdU[nSpecies+nDim];
if (rhoCvve != 0.) val_dPdU[nSpecies+nDim+1] += rho_el*Ru/MolarMass[0]*1.0/rhoCvve;

}

Expand Down Expand Up @@ -280,15 +280,21 @@ void CNEMOGas::ComputedTvedU(const su2double *V, const vector<su2double>& val_ev
su2double rhoCvve = V[RHOCVVE_INDEX];

/*--- Species density derivatives ---*/
for (iSpecies = 0; iSpecies < nSpecies; iSpecies++) {
val_dTvedU[iSpecies] = -val_eves[iSpecies]/rhoCvve;
if(rhoCvve == 0.) {
for (iSpecies = 0; iSpecies < nSpecies; iSpecies++)
val_dTvedU[iSpecies] = 0.0;
} else {
for (iSpecies = 0; iSpecies < nSpecies; iSpecies++)
val_dTvedU[iSpecies] = -val_eves[iSpecies]/rhoCvve;
}

/*--- Momentum derivatives ---*/
for (iDim = 0; iDim < nDim; iDim++)
val_dTvedU[nSpecies+iDim] = 0.0;

/*--- Energy derivatives ---*/
val_dTvedU[nSpecies+nDim] = 0.0;
val_dTvedU[nSpecies+nDim+1] = 1.0 / rhoCvve;
if (rhoCvve == 0.) val_dTvedU[nSpecies+nDim+1] = 0.0;
else val_dTvedU[nSpecies+nDim+1] = 1.0 / rhoCvve;

}
18 changes: 12 additions & 6 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 @@ -4,7 +4,7 @@
* Contains methods for common tasks, e.g. compute flux
* Jacobians.
* \author S.R. Copeland, W. Maier, C. Garbacz
* \version 8.3.0 "Harrier"
* \version 8.2.0 "Harrier"
*
* SU2 Project Website: https://su2code.github.io
*
Expand Down Expand Up @@ -257,7 +257,7 @@ void CNEMONumerics::GetViscousProjFlux(const su2double *val_primvar,
/*--- Pre-compute mixture quantities ---*/ //TODO
su2double Vector[MAXNDIM] = {0.0};
for (auto iDim = 0; iDim < nDim; iDim++) {
for (auto iSpecies = nEl; iSpecies < nSpecies; iSpecies++) {
for (auto iSpecies = nEl; iSpecies < nHeavy; iSpecies++) {
Vector[iDim] += rho*Ds[iSpecies]*GV[RHOS_INDEX+iSpecies][iDim];
}
}
Expand All @@ -270,11 +270,11 @@ void CNEMONumerics::GetViscousProjFlux(const su2double *val_primvar,

/*--- Species diffusion velocity ---*/
if (nEl == 1) Flux_Tensor[0][iDim] = 0.0; // Ambipolar diffusion for electron is handled differently than heavy species below
for (auto iSpecies = nEl; iSpecies < nSpecies; iSpecies++) {
for (auto iSpecies = nEl; iSpecies < nHeavy; iSpecies++) {
Flux_Tensor[iSpecies][iDim] = rho*Ds[iSpecies]*GV[RHOS_INDEX+iSpecies][iDim]
- V[RHOS_INDEX+iSpecies]*Vector[iDim];
if (nEl == 1){
Flux_Tensor[0][iDim] += -1.0 * Ms[0] * Flux_Tensor[iSpecies][iDim] * Cs[iSpecies]/ Ms[iSpecies];
Flux_Tensor[0][iDim] += -1.0 * Ms[0] * Flux_Tensor[iSpecies][iDim] / Ms[iSpecies];
}
}

Expand All @@ -286,7 +286,7 @@ void CNEMONumerics::GetViscousProjFlux(const su2double *val_primvar,
}

/*--- Diffusion terms ---*/
for (auto iSpecies = 0ul; iSpecies < nSpecies; iSpecies++) {
for (auto iSpecies = 0ul; 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];
}
Expand Down Expand Up @@ -343,8 +343,14 @@ void CNEMONumerics::GetViscousProjJacs(const su2double *val_Mean_PrimVar,
return COMPUTE_VISCOUS_JACS(9, 5);

case 10:
return COMPUTE_VISCOUS_JACS(10, 5);
switch (nSpecies) {
case 5: return COMPUTE_VISCOUS_JACS(10, 5);

case 6: return COMPUTE_VISCOUS_JACS(10, 6);

default: SU2_MPI::Error("nVar and nSpecies mismatch.", CURRENT_FUNCTION);
}
break;
case 11:
return COMPUTE_VISCOUS_JACS(11, 7);

Expand Down
84 changes: 72 additions & 12 deletions SU2_CFD/src/solvers/CNEMONSSolver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
* \file CNEMONSSolver.cpp
* \brief Headers of the CNEMONSSolver class
* \author S. R. Copeland, F. Palacios, W. Maier.
* \version 8.3.0 "Harrier"
* \version 8.2.0 "Harrier"
*
* SU2 Project Website: https://su2code.github.io
*
Expand Down Expand Up @@ -532,9 +532,13 @@ void CNEMONSSolver::BC_IsothermalNonCatalytic_Wall(CGeometry *geometry,
const bool ionization = config->GetIonization();
su2double UnitNormal[MAXNDIM] = {0.0};

//if (ionization) {
// SU2_MPI::Error("NEED TO TAKE A CLOSER LOOK AT THE JACOBIAN W/ IONIZATION",CURRENT_FUNCTION);
//}
if (ionization) {
SU2_MPI::Error("NEED TO TAKE A CLOSER LOOK AT THE JACOBIAN W/ IONIZATION",CURRENT_FUNCTION);
}

/*--- Extract required indices ---*/
const unsigned short RHOCVTR_INDEX = nodes->GetRhoCvtrIndex();
const unsigned short RHO_INDEX = nodes->GetRhoIndex();

/*--- Define 'proportional control' constant ---*/
const su2double C = 5;
Expand Down Expand Up @@ -769,7 +773,7 @@ void CNEMONSSolver::BC_IsothermalCatalytic_Wall(CGeometry *geometry,
// Temperature
for (auto iSpecies = 0ul; iSpecies < nSpecies; iSpecies++) {
for (auto jSpecies = 0ul; jSpecies < nSpecies; jSpecies++) {
Jacobian_j[nSpecies+nDim][iSpecies] += Jacobian_j[jSpecies][iSpecies]*hs[iSpecies];
Jacobian_j[nSpecies+nDim][iSpecies] += Jacobian_j[jSpecies][iSpecies]*hs[jSpecies];
}
Jacobian_j[nSpecies+nDim][nSpecies+nDim] += Res_Visc[iSpecies]/Area*(Ru/Ms[iSpecies] +
Cvtrs[iSpecies] );
Expand All @@ -779,7 +783,7 @@ void CNEMONSSolver::BC_IsothermalCatalytic_Wall(CGeometry *geometry,
// Vib.-El. Temperature
for (auto iSpecies = 0ul; iSpecies < nSpecies; iSpecies++) {
for (auto jSpecies = 0ul; jSpecies < nSpecies; jSpecies++)
Jacobian_j[nSpecies+nDim+1][iSpecies] += Jacobian_j[jSpecies][iSpecies]*eves[iSpecies];
Jacobian_j[nSpecies+nDim+1][iSpecies] += Jacobian_j[jSpecies][iSpecies]*eves[jSpecies];
Jacobian_j[nSpecies+nDim+1][nSpecies+nDim+1] += Res_Visc[iSpecies]/Area*Cvve[iSpecies];
}

Expand All @@ -795,10 +799,6 @@ void CNEMONSSolver::BC_IsothermalCatalytic_Wall(CGeometry *geometry,

} else {

if (implicit) {
SU2_MPI::Error("Implicit not currently available for partially catalytic wall.",CURRENT_FUNCTION);
}

/*--- Identify the boundary ---*/
string Marker_Tag = config->GetMarker_All_TagBound(val_marker);

Expand All @@ -821,11 +821,71 @@ void CNEMONSSolver::BC_IsothermalCatalytic_Wall(CGeometry *geometry,
int Index = SU2_TYPE::Int(RxnTable(iSpecies,1));
Res_Visc[iSpecies] = RxnTable(iSpecies,0)*factor*Vi[Index]/Vi[RHO_INDEX]*sqrt(1/Ms[Index]);
}

if (implicit) {
/*--- Initialize the transformation matrix ---*/
for (auto iVar = 0ul; iVar < nVar; iVar++) {
for (auto jVar = 0ul; jVar < nVar; jVar++) {
dVdU[iVar][jVar] = 0.0;
Jacobian_j[iVar][jVar] = 0.0;
Jacobian_i[iVar][jVar] = 0.0;
}
}

for (auto iSpecies = 0ul; iSpecies < nSpecies; iSpecies++) {
for (auto jSpecies = 0ul; jSpecies < nSpecies; jSpecies++) {
dVdU[iSpecies][iSpecies] = 1.0;
}
}
for (auto iVar = 0ul; iVar < nVar; iVar++) {
dVdU[nSpecies+nDim][iVar] = dTdU[iVar];
dVdU[nSpecies+nDim+1][iVar] = dTvedU[iVar];
}

const auto& Cvtrs = FluidModel->GetSpeciesCvTraRot();
const auto& Cvve = nodes->GetCvve(iPoint);

// Species density
for (auto iSpecies = 0ul; iSpecies < nSpecies; iSpecies++) {
int Index = SU2_TYPE::Int(RxnTable(iSpecies,1));
for (auto jSpecies = 0ul; jSpecies < nSpecies; jSpecies++) {
if (jSpecies == Index)
Jacobian_j[iSpecies][jSpecies] = RxnTable(iSpecies,0)*gam*sqrt(RuSI*Tw/2/PI_NUMBER/Ms[Index]);
}
}

// Temperature
for (auto iSpecies = 0ul; iSpecies < nSpecies; iSpecies++) {
for (auto jSpecies = 0ul; jSpecies < nSpecies; jSpecies++) {
Jacobian_j[nSpecies+nDim][iSpecies] += Jacobian_j[jSpecies][iSpecies]*hs[jSpecies];
}
Jacobian_j[nSpecies+nDim][nSpecies+nDim] += Res_Visc[iSpecies]/Area*(Ru/Ms[iSpecies] +
Cvtrs[iSpecies] );
Jacobian_j[nSpecies+nDim][nSpecies+nDim+1] += Res_Visc[iSpecies]/Area*Cvve[iSpecies];
}

// Vib.-El. Temperature
for (auto iSpecies = 0ul; iSpecies < nSpecies; iSpecies++) {
for (auto jSpecies = 0ul; jSpecies < nSpecies; jSpecies++)
Jacobian_j[nSpecies+nDim+1][iSpecies] += Jacobian_j[jSpecies][iSpecies]*eves[jSpecies];
Jacobian_j[nSpecies+nDim+1][nSpecies+nDim+1] += Res_Visc[iSpecies]/Area*Cvve[iSpecies];
}

/*--- Multiply by the transformation matrix and store in Jac. ii ---*/
for (auto iVar = 0ul; iVar < nVar; iVar++)
for (auto jVar = 0ul; jVar < nVar; jVar++)
for (auto kVar = 0ul; kVar < nVar; kVar++)
Jacobian_i[iVar][jVar] += Jacobian_j[iVar][kVar]*dVdU[kVar][jVar]*Area;

/*--- Apply to the linear system ---*/
Jacobian.SubtractBlock(iPoint, iPoint, Jacobian_i);
}

}

for (auto iSpecies = 0ul; iSpecies < nSpecies; iSpecies++) {
Res_Visc[nSpecies+nDim] += (Res_Visc[iSpecies]*hs[iSpecies])*Area;
Res_Visc[nSpecies+nDim+1] += (Res_Visc[iSpecies]*eves[iSpecies])*Area;
Res_Visc[nSpecies+nDim] += (Res_Visc[iSpecies]*hs[iSpecies]);
Res_Visc[nSpecies+nDim+1] += (Res_Visc[iSpecies]*eves[iSpecies]);
}

/*--- Viscous contribution to the residual at the wall ---*/
Expand Down
Loading