Skip to content

Commit 629902a

Browse files
jtneedelsjtneedels
authored
Merge pull request #1162 from su2code/feature_nemo_axi_viscous
Addition of Source Terms for NEMO Axisymmetric Flows
2 parents a785353 + 925bd40 commit 629902a

3 files changed

Lines changed: 166 additions & 23 deletions

File tree

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

Lines changed: 57 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -296,7 +296,25 @@ 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 *Ds, **GV, ktr, kve;
301+
302+
/*--- Rename for convenience ---*/
303+
Ds = Diffusion_Coeff_i;
304+
ktr = Thermal_Conductivity_i;
305+
kve = Thermal_Conductivity_ve_i;
306+
rho = V_i[RHO_INDEX];
307+
T = V_i[T_INDEX];
308+
Tve = V_i[TVE_INDEX];
309+
GV = PrimVar_Grad_i;
310+
RuSI= UNIVERSAL_GAS_CONSTANT;
311+
Ru = 1000.0*RuSI;
312+
313+
auto& Ms = fluidmodel->GetSpeciesMolarMass();
314+
315+
bool viscous = config->GetViscous();
316+
bool rans = (config->GetKind_Turb_Model() != NONE);
317+
hs = fluidmodel->ComputeSpeciesEnthalpy(T, Tve, eve_i);
300318

301319
/*--- Initialize residual and Jacobian arrays ---*/
302320
for (iVar = 0; iVar < nVar; iVar++) {
@@ -312,6 +330,7 @@ CNumerics::ResidualType<> CSource_NEMO::ComputeAxisymmetric(const CConfig *confi
312330
rhov = U_i[nSpecies+1];
313331
H = V_i[H_INDEX];
314332
vel2 = 0.0;
333+
315334
for (iDim = 0; iDim < nDim; iDim++)
316335
vel2 += V_i[VEL_INDEX+iDim]*V_i[VEL_INDEX+iDim];
317336
for (iSpecies = 0; iSpecies < nSpecies; iSpecies++)
@@ -324,6 +343,43 @@ CNumerics::ResidualType<> CSource_NEMO::ComputeAxisymmetric(const CConfig *confi
324343
residual[nSpecies+2] = yinv*rhov*H*Volume;
325344
residual[nSpecies+3] = yinv*rhov*U_i[nSpecies+nDim+1]/rho*Volume;
326345

346+
if (viscous) {
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+
su2double qy_ve = kve*GV[TVE_INDEX][1];
362+
363+
/*--- Enthalpy and vib-el energy transport due to y-direction diffusion---*/
364+
for (iSpecies = 0; iSpecies < nHeavy; iSpecies++) {
365+
sumJhs_y += (rho*Ds[iSpecies]*GV[RHOS_INDEX+iSpecies][1] - V_i[RHOS_INDEX+iSpecies]*Vector[1]) * hs[iSpecies];
366+
sumJeve_y += (rho*Ds[iSpecies]*GV[RHOS_INDEX+iSpecies][1] - V_i[RHOS_INDEX+iSpecies]*Vector[1]) * eve_i[iSpecies];
367+
}
368+
369+
for (iSpecies = 0; iSpecies < nSpecies; iSpecies++)
370+
residual[iSpecies] -= 0.0;
371+
residual[nSpecies] -= Volume*(yinv*total_viscosity_i*(PrimVar_Grad_i[nSpecies+2][1]+PrimVar_Grad_i[nSpecies+3][0])
372+
-TWO3*AuxVar_Grad_i[0][0]);
373+
residual[nSpecies+1] -= Volume*(yinv*total_viscosity_i*2*(PrimVar_Grad_i[nSpecies+3][1]-v*yinv)
374+
-TWO3*AuxVar_Grad_i[0][1]);
375+
residual[nSpecies+2] -= Volume*(yinv*(- sumJhs_y + total_viscosity_i*(u*(PrimVar_Grad_i[nSpecies+3][0]+PrimVar_Grad_i[nSpecies+2][1])
376+
+v*TWO3*(2*PrimVar_Grad_i[nSpecies+2][1]-PrimVar_Grad_i[nSpecies+2][0]
377+
-v*yinv+rho*turb_ke_i))
378+
-total_conductivity_i*PrimVar_Grad_i[nSpecies][1])
379+
-TWO3*(AuxVar_Grad_i[1][1]+AuxVar_Grad_i[2][1]));
380+
residual[nSpecies+3] -= Volume*(yinv*(-sumJeve_y -qy_ve));
381+
}
382+
327383
// if (implicit) {
328384
//
329385
// /*--- Initialize ---*/

‎SU2_CFD/src/solvers/CNEMOEulerSolver.cpp‎

Lines changed: 102 additions & 22 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,25 +1153,102 @@ 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+
}
1181+
}
1182+
1183+
/*--- loop over points ---*/
1184+
SU2_OMP_FOR_DYN(omp_chunk_size)
1185+
for (iPoint = 0; iPoint < nPointDomain; iPoint++) {
1186+
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 necessary, set variables needed for viscous computation ---*/
1200+
if (viscous) {
1201+
1202+
/*--- Set gradient of primitive variables ---*/
1203+
numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), nodes->GetGradient_Primitive(iPoint));
1204+
1205+
/*--- Set gradient of auxillary variables ---*/
1206+
numerics->SetAuxVarGrad(nodes->GetAuxVarGradient(iPoint), nullptr);
1207+
1208+
/*--- Set diffusion coefficient ---*/
1209+
numerics->SetDiffusionCoeff(nodes->GetDiffusionCoeff(iPoint), nodes->GetDiffusionCoeff(iPoint));
1210+
1211+
/*--- Laminar viscosity ---*/
1212+
numerics->SetLaminarViscosity(nodes->GetLaminarViscosity(iPoint), nodes->GetLaminarViscosity(iPoint));
1213+
1214+
/*--- Eddy viscosity ---*/
1215+
numerics->SetEddyViscosity(nodes->GetEddyViscosity(iPoint), nodes->GetEddyViscosity(iPoint));
1216+
1217+
/*--- Thermal conductivity ---*/
1218+
numerics->SetThermalConductivity(nodes->GetThermalConductivity(iPoint), nodes->GetThermalConductivity(iPoint));
1219+
1220+
/*--- Vib-el. thermal conductivity ---*/
1221+
numerics->SetThermalConductivity_ve(nodes->GetThermalConductivity_ve(iPoint), nodes->GetThermalConductivity_ve(iPoint));
1222+
1223+
/*--- Vib-el energy ---*/
1224+
numerics->SetEve(nodes->GetEve(iPoint), nodes->GetEve(iPoint));
1225+
1226+
/*--- Set turbulence kinetic energy ---*/
1227+
if (rans){
1228+
CVariable* turbNodes = solver_container[TURB_SOL]->GetNodes();
1229+
numerics->SetTurbKineticEnergy(turbNodes->GetSolution(iPoint,0), turbNodes->GetSolution(iPoint,0));
1230+
}
1231+
}
1232+
1233+
auto residual = numerics->ComputeAxisymmetric(config);
1234+
1235+
/*--- Check for errors before applying source to the linear system ---*/
1236+
err = false;
1237+
for (iVar = 0; iVar < nVar; iVar++)
1238+
if (residual[iVar] != residual[iVar]) err = true;
1239+
if (implicit)
1240+
for (iVar = 0; iVar < nVar; iVar++)
1241+
for (jVar = 0; jVar < nVar; jVar++)
1242+
if (Jacobian_i[iVar][jVar] != Jacobian_i[iVar][jVar]) err = true;
1243+
1244+
/*--- Apply the update to the linear system ---*/
1245+
if (!err) {
1246+
LinSysRes.AddBlock(iPoint, residual);
1247+
if (implicit)
1248+
Jacobian.AddBlock(iPoint, iPoint, Jacobian_i);
1249+
}else
1250+
eAxi_local++;
11691251
}
1170-
else
1171-
eAxi_local++;
11721252
}
11731253

11741254
if(!monoatomic){

‎SU2_CFD/src/variables/CNEMOEulerVariable.cpp‎

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

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

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

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

0 commit comments

Comments
 (0)