From b89450fedd49def26047dc64324bb14760f9b3f9 Mon Sep 17 00:00:00 2001 From: Raghava Davuluri Date: Wed, 4 Feb 2026 22:18:34 -0500 Subject: [PATCH] Update 1T and Catalytic boundary condition --- SU2_CFD/src/fluid/CMutationTCLib.cpp | 36 ++++----- SU2_CFD/src/fluid/CNEMOGas.cpp | 18 +++-- SU2_CFD/src/numerics/NEMO/CNEMONumerics.cpp | 18 +++-- SU2_CFD/src/solvers/CNEMONSSolver.cpp | 84 ++++++++++++++++++--- 4 files changed, 115 insertions(+), 41 deletions(-) diff --git a/SU2_CFD/src/fluid/CMutationTCLib.cpp b/SU2_CFD/src/fluid/CMutationTCLib.cpp index a38607b965fb..e86c67293413 100644 --- a/SU2_CFD/src/fluid/CMutationTCLib.cpp +++ b/SU2_CFD/src/fluid/CMutationTCLib.cpp @@ -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 * @@ -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"); @@ -141,7 +138,7 @@ void CMutationTCLib::SetTDStateRhosTTv(vector& val_rhos, su2double va Pressure = ComputePressure(); - mix->setState(rhos.data(), temperatures.data(), 1, NEWTON_ROBUST); + mix->setState(rhos.data(), temperatures.data(), 1, NEWTON); } @@ -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(){ @@ -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; @@ -299,7 +301,7 @@ vector& CMutationTCLib::ComputeTemperatures(vector& 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()); diff --git a/SU2_CFD/src/fluid/CNEMOGas.cpp b/SU2_CFD/src/fluid/CNEMOGas.cpp index cfce16b052e5..c457776f79e4 100644 --- a/SU2_CFD/src/fluid/CNEMOGas.cpp +++ b/SU2_CFD/src/fluid/CNEMOGas.cpp @@ -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 * @@ -232,8 +232,8 @@ void CNEMOGas::ComputedPdU(const su2double *V, const vector& 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; } @@ -280,15 +280,21 @@ void CNEMOGas::ComputedTvedU(const su2double *V, const vector& 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; } diff --git a/SU2_CFD/src/numerics/NEMO/CNEMONumerics.cpp b/SU2_CFD/src/numerics/NEMO/CNEMONumerics.cpp index 19a68a64b5fd..ff2906447620 100644 --- a/SU2_CFD/src/numerics/NEMO/CNEMONumerics.cpp +++ b/SU2_CFD/src/numerics/NEMO/CNEMONumerics.cpp @@ -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 * @@ -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]; } } @@ -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]; } } @@ -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]; } @@ -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); diff --git a/SU2_CFD/src/solvers/CNEMONSSolver.cpp b/SU2_CFD/src/solvers/CNEMONSSolver.cpp index 7220351028cd..8933c8287dd3 100644 --- a/SU2_CFD/src/solvers/CNEMONSSolver.cpp +++ b/SU2_CFD/src/solvers/CNEMONSSolver.cpp @@ -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 * @@ -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; @@ -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] ); @@ -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]; } @@ -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); @@ -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 ---*/