Skip to content

Commit 2f91c39

Browse files
jtneedelsjtneedels
committed
Implemented additional source residual terms for NEMO axisymmetric simulation
Signed-off-by: jtneedels <jneedels@stanford.edu>
1 parent 1f1e6e3 commit 2f91c39

3 files changed

Lines changed: 166 additions & 44 deletions

File tree

‎SU2_CFD/src/numerics/NEMO/NEMO_sources.cpp‎

Lines changed: 58 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -296,7 +296,24 @@ CNumerics::ResidualType<> CSource_NEMO::ComputeVibRelaxation(const CConfig *conf
296296
CNumerics::ResidualType<> CSource_NEMO::ComputeAxisymmetric(const CConfig *config) {
297297

298298
unsigned short iDim, iSpecies, iVar;
299-
su2double rho, rhov, vel2, H, yinv;
299+
su2double rho, rhov, vel2, H, yinv, T, Tve, qy_ve, Ru, RuSI;
300+
su2double *V, *Ds, **GV, mu, ktr, kve, div_vel;
301+
302+
auto& Ms = fluidmodel->GetSpeciesMolarMass();
303+
304+
bool viscous = config->GetViscous();
305+
bool rans = (config->GetKind_Turb_Model() != NONE);
306+
hs = fluidmodel->ComputeSpeciesEnthalpy(T, Tve, eve_i);
307+
308+
Ds = Diffusion_Coeff_i;
309+
ktr = Thermal_Conductivity_i;
310+
kve = Thermal_Conductivity_ve_i;
311+
rho = V_i[RHO_INDEX];
312+
T = V_i[T_INDEX];
313+
Tve = V_i[TVE_INDEX];
314+
GV = PrimVar_Grad_i;
315+
RuSI= UNIVERSAL_GAS_CONSTANT;
316+
Ru = 1000.0*RuSI;
300317

301318
/*--- Initialize residual and Jacobian arrays ---*/
302319
for (iVar = 0; iVar < nVar; iVar++) {
@@ -312,6 +329,7 @@ CNumerics::ResidualType<> CSource_NEMO::ComputeAxisymmetric(const CConfig *confi
312329
rhov = U_i[nSpecies+1];
313330
H = V_i[H_INDEX];
314331
vel2 = 0.0;
332+
315333
for (iDim = 0; iDim < nDim; iDim++)
316334
vel2 += V_i[VEL_INDEX+iDim]*V_i[VEL_INDEX+iDim];
317335
for (iSpecies = 0; iSpecies < nSpecies; iSpecies++)
@@ -324,6 +342,45 @@ CNumerics::ResidualType<> CSource_NEMO::ComputeAxisymmetric(const CConfig *confi
324342
residual[nSpecies+2] = yinv*rhov*H*Volume;
325343
residual[nSpecies+3] = yinv*rhov*U_i[nSpecies+nDim+1]/rho*Volume;
326344

345+
if (viscous) {
346+
347+
if (!rans){ turb_ke_i = 0.0; }
348+
349+
su2double sumJhs_y = 0.0;
350+
su2double sumJeve_y = 0.0;
351+
su2double Mass = 0.0;
352+
353+
for (iSpecies=0; iSpecies<nSpecies; iSpecies++)
354+
Mass += V_i[iSpecies]*Ms[iSpecies];
355+
356+
su2double heat_capacity_cp_i = V_i[RHOCVTR_INDEX]+Ru/Mass;
357+
su2double total_viscosity_i = Laminar_Viscosity_i + Eddy_Viscosity_i;
358+
su2double total_conductivity_i = ktr + kve + heat_capacity_cp_i*Eddy_Viscosity_i/Prandtl_Turb;
359+
su2double u = V_i[VEL_INDEX];
360+
su2double v = V_i[VEL_INDEX+1];
361+
362+
qy_ve = kve*GV[TVE_INDEX][1];
363+
364+
/*--- Enthalpy and vib-el energy transport due to y-direction diffusion---*/
365+
for (iSpecies = 0; iSpecies < nHeavy; iSpecies++) {
366+
sumJhs_y += (rho*Ds[iSpecies]*GV[RHOS_INDEX+iSpecies][1] - V_i[RHOS_INDEX+iSpecies]*Vector[1]) * hs[iSpecies];
367+
sumJeve_y += (rho*Ds[iSpecies]*GV[RHOS_INDEX+iSpecies][1] - V_i[RHOS_INDEX+iSpecies]*Vector[1]) * eve_i[iSpecies];
368+
}
369+
370+
for (iSpecies = 0; iSpecies < nSpecies; iSpecies++)
371+
residual[iSpecies] -= 0.0;
372+
residual[nSpecies] -= Volume*(yinv*total_viscosity_i*(PrimVar_Grad_i[nSpecies+2][1]+PrimVar_Grad_i[nSpecies+3][0])
373+
-TWO3*AuxVar_Grad_i[0][0]);
374+
residual[nSpecies+1] -= Volume*(yinv*total_viscosity_i*2*(PrimVar_Grad_i[nSpecies+3][1]-v*yinv)
375+
-TWO3*AuxVar_Grad_i[0][1]);
376+
residual[nSpecies+2] -= Volume*(yinv*(- sumJhs_y + total_viscosity_i*(u*(PrimVar_Grad_i[nSpecies+3][0]+PrimVar_Grad_i[nSpecies+2][1])
377+
+v*TWO3*(2*PrimVar_Grad_i[nSpecies+2][1]-PrimVar_Grad_i[nSpecies+2][0]
378+
-v*yinv+rho*turb_ke_i))
379+
-total_conductivity_i*PrimVar_Grad_i[nSpecies][1])
380+
-TWO3*(AuxVar_Grad_i[1][1]+AuxVar_Grad_i[2][1]));
381+
residual[nSpecies+3] -= Volume*(yinv*(-sumJeve_y -qy_ve));
382+
}
383+
327384
// if (implicit) {
328385
//
329386
// /*--- Initialize ---*/

‎SU2_CFD/src/solvers/CNEMOEulerSolver.cpp‎

Lines changed: 100 additions & 43 deletions
Original file line numberDiff line numberDiff line change
@@ -1108,16 +1108,19 @@ bool CNEMOEulerSolver::CheckNonPhys(const su2double *V) const {
11081108

11091109
void CNEMOEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_container, CNumerics **numerics_container, CConfig *config, unsigned short iMesh) {
11101110

1111-
unsigned short iVar;
1111+
unsigned short iVar, jVar;
11121112
unsigned long iPoint;
11131113
unsigned long eAxi_local, eChm_local, eVib_local;
11141114
unsigned long eAxi_global, eChm_global, eVib_global;
11151115

11161116
/*--- Assign booleans ---*/
11171117
bool err = false;
1118-
//bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT);
1118+
bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT);
11191119
bool frozen = config->GetFrozen();
11201120
bool monoatomic = config->GetMonoatomic();
1121+
bool viscous = config->GetViscous();
1122+
bool ideal_gas = (config->GetKind_FluidModel() == STANDARD_AIR) || (config->GetKind_FluidModel() == IDEAL_GAS);
1123+
bool rans = (config->GetKind_Turb_Model() != NONE);
11211124

11221125
CNumerics* numerics = numerics_container[SOURCE_FIRST_TERM];
11231126

@@ -1137,11 +1140,11 @@ void CNEMOEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_con
11371140
numerics->SetPrimitive (nodes->GetPrimitive(iPoint), nodes->GetPrimitive(iPoint) );
11381141

11391142
/*--- Pass supplementary information to CNumerics ---*/
1140-
numerics->SetdPdU (nodes->GetdPdU(iPoint), nodes->GetdPdU(iPoint));
1141-
numerics->SetdTdU (nodes->GetdTdU(iPoint), nodes->GetdTdU(iPoint));
1143+
numerics->SetdPdU(nodes->GetdPdU(iPoint), nodes->GetdPdU(iPoint));
1144+
numerics->SetdTdU(nodes->GetdTdU(iPoint), nodes->GetdTdU(iPoint));
11421145
numerics->SetdTvedU(nodes->GetdTvedU(iPoint), nodes->GetdTvedU(iPoint));
1143-
numerics->SetEve (nodes->GetEve(iPoint), nodes->GetEve(iPoint));
1144-
numerics->SetCvve (nodes->GetCvve(iPoint), nodes->GetCvve(iPoint));
1146+
numerics->SetEve(nodes->GetEve(iPoint), nodes->GetEve(iPoint));
1147+
numerics->SetCvve(nodes->GetCvve(iPoint), nodes->GetCvve(iPoint));
11451148

11461149
/*--- Set volume of the dual grid cell ---*/
11471150
numerics->SetVolume(geometry->nodes->GetVolume(iPoint));
@@ -1150,49 +1153,103 @@ void CNEMOEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_con
11501153

11511154
/*--- Compute axisymmetric source terms (if needed) ---*/
11521155
if (config->GetAxisymmetric()) {
1153-
auto residual = numerics->ComputeAxisymmetric(config);
11541156

1155-
/*--- Check for errors before applying source to the linear system ---*/
1156-
err = false;
1157-
for (iVar = 0; iVar < nVar; iVar++)
1158-
if (residual[iVar] != residual[iVar]) err = true;
1159-
//if (implicit)
1160-
// for (iVar = 0; iVar < nVar; iVar++)
1161-
// for (jVar = 0; jVar < nVar; jVar++)
1162-
// if (Jacobian_i[iVar][jVar] != Jacobian_i[iVar][jVar]) err = true;
1157+
if (viscous) {
11631158

1164-
/*--- Apply the update to the linear system ---*/
1165-
if (!err) {
1166-
LinSysRes.AddBlock(iPoint, residual);
1167-
//if (implicit)
1168-
// Jacobian.AddBlock(iPoint, iPoint, Jacobian_i);
1159+
for (iPoint = 0; iPoint < nPoint; iPoint++) {
1160+
1161+
su2double yCoord = geometry->nodes->GetCoord(iPoint, 1);
1162+
su2double yVelocity = nodes->GetVelocity(iPoint,1);
1163+
su2double xVelocity = nodes->GetVelocity(iPoint,0);
1164+
su2double Total_Viscosity = nodes->GetLaminarViscosity(iPoint) + nodes->GetEddyViscosity(iPoint);
1165+
1166+
if (yCoord > EPS){
1167+
su2double nu_v_on_y = Total_Viscosity*yVelocity/yCoord;
1168+
nodes->SetAuxVar(iPoint, 0, nu_v_on_y);
1169+
nodes->SetAuxVar(iPoint, 1, nu_v_on_y*yVelocity);
1170+
nodes->SetAuxVar(iPoint, 2, nu_v_on_y*xVelocity);
1171+
}
1172+
}
1173+
1174+
/*--- Compute the auxiliary variable gradient with GG or WLS. ---*/
1175+
if (config->GetKind_Gradient_Method() == GREEN_GAUSS) {
1176+
SetAuxVar_Gradient_GG(geometry, config);
1177+
}
1178+
if (config->GetKind_Gradient_Method() == WEIGHTED_LEAST_SQUARES) {
1179+
SetAuxVar_Gradient_LS(geometry, config);
1180+
}
11691181
}
1170-
else
1171-
eAxi_local++;
1172-
}
11731182

1174-
if(!monoatomic){
1175-
if(!frozen){
1176-
/*--- Compute the non-equilibrium chemistry ---*/
1177-
auto residual = numerics->ComputeChemistry(config);
1183+
/*--- loop over points ---*/
1184+
SU2_OMP_FOR_DYN(omp_chunk_size)
1185+
for (iPoint = 0; iPoint < nPointDomain; iPoint++) {
11781186

1179-
/*--- Check for errors before applying source to the linear system ---*/
1180-
err = false;
1181-
for (iVar = 0; iVar < nVar; iVar++)
1182-
if (residual[iVar] != residual[iVar]) err = true;
1183-
//if (implicit)
1184-
// for (iVar = 0; iVar < nVar; iVar++)
1185-
// for (jVar = 0; jVar < nVar; jVar++)
1186-
// if (Jacobian_i[iVar][jVar] != Jacobian_i[iVar][jVar]) err = true;
1187-
1188-
/*--- Apply the chemical sources to the linear system ---*/
1189-
if (!err) {
1190-
LinSysRes.SubtractBlock(iPoint, residual);
1191-
//if (implicit)
1192-
// Jacobian.SubtractBlock(iPoint, iPoint, Jacobian_i);
1193-
} else
1194-
eChm_local++;
1187+
/*--- Set solution ---*/
1188+
numerics->SetConservative(nodes->GetSolution(iPoint), nodes->GetSolution(iPoint));
1189+
1190+
/*--- Set control volume ---*/
1191+
numerics->SetVolume(geometry->nodes->GetVolume(iPoint));
1192+
1193+
/*--- Set y coordinate ---*/
1194+
numerics->SetCoord(geometry->nodes->GetCoord(iPoint), geometry->nodes->GetCoord(iPoint));
1195+
1196+
/*--- Set primitive variables for viscous terms and/or generalised source ---*/
1197+
numerics->SetPrimitive(nodes->GetPrimitive(iPoint), nodes->GetPrimitive(iPoint));
1198+
1199+
if (viscous) {
1200+
1201+
/*--- Set gradient of primitive variables ---*/
1202+
numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), nodes->GetGradient_Primitive(iPoint));
1203+
1204+
/*--- Set gradient of auxillary variables ---*/
1205+
numerics->SetAuxVarGrad(nodes->GetAuxVarGradient(iPoint), nullptr);
1206+
1207+
/*--- Set diffusion coefficient ---*/
1208+
numerics->SetDiffusionCoeff(nodes->GetDiffusionCoeff(iPoint), nodes->GetDiffusionCoeff(iPoint));
1209+
1210+
/*--- Laminar viscosity ---*/
1211+
numerics->SetLaminarViscosity(nodes->GetLaminarViscosity(iPoint), nodes->GetLaminarViscosity(iPoint));
1212+
1213+
/*--- Eddy viscosity ---*/
1214+
numerics->SetEddyViscosity(nodes->GetEddyViscosity(iPoint), nodes->GetEddyViscosity(iPoint));
1215+
1216+
/*--- Thermal conductivity ---*/
1217+
numerics->SetThermalConductivity(nodes->GetThermalConductivity(iPoint), nodes->GetThermalConductivity(iPoint));
1218+
1219+
/*--- Vib-el. thermal conductivity ---*/
1220+
numerics->SetThermalConductivity_ve(nodes->GetThermalConductivity_ve(iPoint), nodes->GetThermalConductivity_ve(iPoint));
1221+
1222+
/*--- Vib-el energy ---*/
1223+
numerics->SetEve(nodes->GetEve(iPoint), nodes->GetEve(iPoint));
1224+
1225+
/*--- Set turbulence kinetic energy ---*/
1226+
if (rans){
1227+
CVariable* turbNodes = solver_container[TURB_SOL]->GetNodes();
1228+
numerics->SetTurbKineticEnergy(turbNodes->GetSolution(iPoint,0), turbNodes->GetSolution(iPoint,0));
1229+
}
1230+
}
1231+
1232+
auto residual = numerics->ComputeAxisymmetric(config);
1233+
1234+
/*--- Check for errors before applying source to the linear system ---*/
1235+
err = false;
1236+
for (iVar = 0; iVar < nVar; iVar++)
1237+
if (residual[iVar] != residual[iVar]) err = true;
1238+
if (implicit)
1239+
for (iVar = 0; iVar < nVar; iVar++)
1240+
for (jVar = 0; jVar < nVar; jVar++)
1241+
if (Jacobian_i[iVar][jVar] != Jacobian_i[iVar][jVar]) err = true;
1242+
1243+
/*--- Apply the update to the linear system ---*/
1244+
if (!err) {
1245+
LinSysRes.AddBlock(iPoint, residual);
1246+
if (implicit)
1247+
Jacobian.AddBlock(iPoint, iPoint, Jacobian_i);
1248+
}
1249+
else
1250+
eAxi_local++;
11951251
}
1252+
11961253
}
11971254

11981255
/*--- Compute vibrational energy relaxation ---*/

‎SU2_CFD/src/variables/CNEMOEulerVariable.cpp‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -78,6 +78,14 @@ CNEMOEulerVariable::CNEMOEulerVariable(su2double val_pressure,
7878

7979
Res_TruncError.resize(nPoint,nVar) = su2double(0.0);
8080

81+
/*--- Size Grad_AuxVar for axiysmmetric ---*/
82+
83+
if (config->GetAxisymmetric()){
84+
nAuxVar = 3;
85+
Grad_AuxVar.resize(nPoint,nAuxVar,nDim,0.0);
86+
AuxVar.resize(nPoint,nAuxVar) = su2double(0.0);
87+
}
88+
8189
/*--- Only for residual smoothing (multigrid) ---*/
8290

8391
for (unsigned long iMesh = 0; iMesh <= config->GetnMGLevels(); iMesh++) {

0 commit comments

Comments
 (0)