Report heat fluxes as imposed by wall boundary conditions by pcarruscag · Pull Request #2394 · su2code/SU2 · GitHub
Skip to content

Report heat fluxes as imposed by wall boundary conditions #2394

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
pcarruscag merged 8 commits into from
Dec 22, 2024
Merged

Report heat fluxes as imposed by wall boundary conditions #2394

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
bef2a1e
measure the heat flux used in BCs instead of computing from the norma…
pcarruscag Dec 8, 2024
016df92
Merge remote-tracking branch 'upstream/develop' into improve_heat_flux
pcarruscag Dec 9, 2024
397c228
recover old sign convention and use helper function
pcarruscag Dec 15, 2024
82ac791
update tests
pcarruscag Dec 15, 2024
d476c8b
update file regression
pcarruscag Dec 16, 2024
2212d57
Merge remote-tracking branch 'upstream/develop' into improve_heat_flux
pcarruscag Dec 21, 2024
0fc99e9
re-update
pcarruscag Dec 22, 2024
91fdb3f
Merge branch 'develop' into improve_heat_flux
pcarruscag Dec 22, 2024
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
8 changes: 8 additions & 0 deletions Common/include/toolboxes/geometry_toolbox.hpp
View file
Open in desktop
47 changes: 34 additions & 13 deletions SU2_CFD/include/solvers/CFVMFlowSolverBase.inl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2357,12 +2357,11 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
const su2double *Coord = nullptr, *Coord_Normal = nullptr, *Normal = nullptr;
const su2double minYPlus = config->GetwallModel_MinYPlus();

string Marker_Tag, Monitoring_Tag;

const su2double Alpha = config->GetAoA() * PI_NUMBER / 180.0;
const su2double Beta = config->GetAoS() * PI_NUMBER / 180.0;
const su2double RefLength = config->GetRefLength();
const su2double RefHeatFlux = config->GetHeat_Flux_Ref();
const su2double RefTemperature = config->GetTemperature_Ref();
const su2double Gas_Constant = config->GetGas_ConstantND();
auto Origin = config->GetRefOriginMoment(0);

Expand Down Expand Up @@ -2393,15 +2392,17 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr

for (iMarker = 0; iMarker < nMarker; iMarker++) {

Marker_Tag = config->GetMarker_All_TagBound(iMarker);
if (!config->GetViscous_Wall(iMarker)) continue;
const auto Marker_Tag = config->GetMarker_All_TagBound(iMarker);

const bool py_custom = config->GetMarker_All_PyCustom(iMarker);

/*--- Obtain the origin for the moment computation for a particular marker ---*/

const auto Monitoring = config->GetMarker_All_Monitoring(iMarker);
if (Monitoring == YES) {
for (iMarker_Monitoring = 0; iMarker_Monitoring < config->GetnMarker_Monitoring(); iMarker_Monitoring++) {
Monitoring_Tag = config->GetMarker_Monitoring_TagBound(iMarker_Monitoring);
const auto Monitoring_Tag = config->GetMarker_Monitoring_TagBound(iMarker_Monitoring);
if (Marker_Tag == Monitoring_Tag) Origin = config->GetRefOriginMoment(iMarker_Monitoring);
}
}
Expand Down Expand Up @@ -2506,26 +2507,47 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr

/*--- Compute total and maximum heat flux on the wall ---*/

su2double dTdn = -GeometryToolbox::DotProduct(nDim, Grad_Temp, UnitNormal);

if (!nemo){

if (!nemo) {
if (FlowRegime == ENUM_REGIME::COMPRESSIBLE) {

Cp = (Gamma / Gamma_Minus_One) * Gas_Constant;
thermal_conductivity = Cp * Viscosity / Prandtl_Lam;
}
if (FlowRegime == ENUM_REGIME::INCOMPRESSIBLE) {
if (!energy) dTdn = 0.0;
thermal_conductivity = nodes->GetThermalConductivity(iPoint);
}
HeatFlux[iMarker][iVertex] = -thermal_conductivity * dTdn * RefHeatFlux;

if (config->GetMarker_All_KindBC(iMarker) == BC_TYPE::HEAT_FLUX) {
if (py_custom) {
HeatFlux[iMarker][iVertex] = -geometry->GetCustomBoundaryHeatFlux(iMarker, iVertex);
} else {
HeatFlux[iMarker][iVertex] = -config->GetWall_HeatFlux(Marker_Tag);
if (config->GetIntegrated_HeatFlux()) {
HeatFlux[iMarker][iVertex] /= geometry->GetSurfaceArea(config, iMarker);
}
}
} else if (config->GetMarker_All_KindBC(iMarker) == BC_TYPE::ISOTHERMAL) {
su2double Twall = 0.0;
if (py_custom) {
Twall = geometry->GetCustomBoundaryTemperature(iMarker, iVertex) / RefTemperature;
} else {
Twall = config->GetIsothermal_Temperature(Marker_Tag) / RefTemperature;
}
iPointNormal = geometry->vertex[iMarker][iVertex]->GetNormal_Neighbor();
Coord_Normal = geometry->nodes->GetCoord(iPointNormal);
const su2double dist_ij = GeometryToolbox::NormalDistance(nDim, UnitNormal, Coord, Coord_Normal);
const su2double There = nodes->GetTemperature(iPointNormal);
HeatFlux[iMarker][iVertex] = thermal_conductivity * (There - Twall) / dist_ij * RefHeatFlux;
} else {
su2double dTdn = GeometryToolbox::DotProduct(nDim, Grad_Temp, UnitNormal);
if (FlowRegime == ENUM_REGIME::INCOMPRESSIBLE && !energy) dTdn = 0.0;
HeatFlux[iMarker][iVertex] = thermal_conductivity * dTdn * RefHeatFlux;
}
} else {

const auto& thermal_conductivity_tr = nodes->GetThermalConductivity(iPoint);
const auto& thermal_conductivity_ve = nodes->GetThermalConductivity_ve(iPoint);

const su2double dTdn = -GeometryToolbox::DotProduct(nDim, Grad_Temp, UnitNormal);
const su2double dTvedn = -GeometryToolbox::DotProduct(nDim, Grad_Temp_ve, UnitNormal);

/*--- Surface energy balance: trans-rot heat flux, vib-el heat flux ---*/
Expand Down Expand Up @@ -2653,8 +2675,7 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
/*--- Compute the coefficients per surface ---*/

for (iMarker_Monitoring = 0; iMarker_Monitoring < config->GetnMarker_Monitoring(); iMarker_Monitoring++) {
Monitoring_Tag = config->GetMarker_Monitoring_TagBound(iMarker_Monitoring);
Marker_Tag = config->GetMarker_All_TagBound(iMarker);
const auto Monitoring_Tag = config->GetMarker_Monitoring_TagBound(iMarker_Monitoring);
if (Marker_Tag == Monitoring_Tag) {
SurfaceViscCoeff.CL[iMarker_Monitoring] += ViscCoeff.CL[iMarker];
SurfaceViscCoeff.CD[iMarker_Monitoring] += ViscCoeff.CD[iMarker];
Expand Down
27 changes: 12 additions & 15 deletions SU2_CFD/src/solvers/CIncNSSolver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -341,7 +341,6 @@ void CIncNSSolver::BC_Wall_Generic(const CGeometry *geometry, const CConfig *con
/*--- Variables for streamwise periodicity ---*/
const bool streamwise_periodic = (config->GetKind_Streamwise_Periodic() != ENUM_STREAMWISE_PERIODIC::NONE);
const bool streamwise_periodic_temperature = config->GetStreamwise_Periodic_Temperature();
su2double Cp, thermal_conductivity, dot_product, scalar_factor;

/*--- Identify the boundary by string name ---*/

Expand Down Expand Up @@ -434,15 +433,17 @@ void CIncNSSolver::BC_Wall_Generic(const CGeometry *geometry, const CConfig *con
/*--- With streamwise periodic flow and heatflux walls an additional term is introduced in the boundary formulation ---*/
if (streamwise_periodic && streamwise_periodic_temperature) {

Cp = nodes->GetSpecificHeatCp(iPoint);
thermal_conductivity = nodes->GetThermalConductivity(iPoint);
const su2double Cp = nodes->GetSpecificHeatCp(iPoint);
const su2double thermal_conductivity = nodes->GetThermalConductivity(iPoint);

/*--- Scalar factor of the residual contribution ---*/
const su2double norm2_translation = GeometryToolbox::SquaredNorm(nDim, config->GetPeriodic_Translation(0));
scalar_factor = SPvals.Streamwise_Periodic_IntegratedHeatFlow*thermal_conductivity / (SPvals.Streamwise_Periodic_MassFlow * Cp * norm2_translation);
const su2double scalar_factor =
SPvals.Streamwise_Periodic_IntegratedHeatFlow*thermal_conductivity /
(SPvals.Streamwise_Periodic_MassFlow * Cp * norm2_translation);

/*--- Dot product ---*/
dot_product = GeometryToolbox::DotProduct(nDim, config->GetPeriodic_Translation(0), Normal);
const su2double dot_product = GeometryToolbox::DotProduct(nDim, config->GetPeriodic_Translation(0), Normal);

LinSysRes(iPoint, nDim+1) += scalar_factor*dot_product;
} // if streamwise_periodic
Expand Down Expand Up @@ -472,18 +473,17 @@ void CIncNSSolver::BC_Wall_Generic(const CGeometry *geometry, const CConfig *con

const auto Coord_i = geometry->nodes->GetCoord(iPoint);
const auto Coord_j = geometry->nodes->GetCoord(Point_Normal);
su2double Edge_Vector[MAXNDIM];
GeometryToolbox::Distance(nDim, Coord_j, Coord_i, Edge_Vector);
su2double dist_ij_2 = GeometryToolbox::SquaredNorm(nDim, Edge_Vector);
su2double dist_ij = sqrt(dist_ij_2);
su2double UnitNormal[MAXNDIM] = {0.0};
for (auto iDim = 0u; iDim < nDim; ++iDim) UnitNormal[iDim] = Normal[iDim] / Area;
const su2double dist_ij = GeometryToolbox::NormalDistance(nDim, UnitNormal, Coord_i, Coord_j);

/*--- Compute the normal gradient in temperature using Twall ---*/

su2double dTdn = -(nodes->GetTemperature(Point_Normal) - Twall)/dist_ij;
const su2double dTdn = -(nodes->GetTemperature(Point_Normal) - Twall)/dist_ij;

/*--- Get thermal conductivity ---*/

su2double thermal_conductivity = nodes->GetThermalConductivity(iPoint);
const su2double thermal_conductivity = nodes->GetThermalConductivity(iPoint);
Comment thread Fixed

/*--- Apply a weak boundary condition for the energy equation.
Compute the residual due to the prescribed heat flux. ---*/
Expand All @@ -493,10 +493,7 @@ void CIncNSSolver::BC_Wall_Generic(const CGeometry *geometry, const CConfig *con
/*--- Jacobian contribution for temperature equation. ---*/

if (implicit) {
su2double proj_vector_ij = 0.0;
if (dist_ij_2 > 0.0)
proj_vector_ij = GeometryToolbox::DotProduct(nDim, Edge_Vector, Normal) / dist_ij_2;
Jacobian.AddVal2Diag(iPoint, nDim+1, thermal_conductivity*proj_vector_ij);
Jacobian.AddVal2Diag(iPoint, nDim+1, thermal_conductivity * Area / dist_ij);
}
break;
} // switch
Expand Down
12 changes: 4 additions & 8 deletions SU2_CFD/src/solvers/CNSSolver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -672,22 +672,19 @@ void CNSSolver::BC_Isothermal_Wall_Generic(CGeometry *geometry, CSolver **solver

const auto Coord_i = geometry->nodes->GetCoord(iPoint);
const auto Coord_j = geometry->nodes->GetCoord(Point_Normal);

su2double dist_ij = GeometryToolbox::Distance(nDim, Coord_i, Coord_j);
const su2double dist_ij = GeometryToolbox::NormalDistance(nDim, UnitNormal, Coord_i, Coord_j);

/*--- Store the corrected velocity at the wall which will
be zero (v = 0), unless there is grid motion (v = u_wall)---*/

if (dynamic_grid) {
nodes->SetVelocity_Old(iPoint, geometry->nodes->GetGridVel(iPoint));
}
else {
} else {
su2double zero[MAXNDIM] = {0.0};
nodes->SetVelocity_Old(iPoint, zero);
}

for (auto iDim = 0u; iDim < nDim; iDim++)
LinSysRes(iPoint, iDim+1) = 0.0;
for (auto iDim = 0u; iDim < nDim; iDim++) LinSysRes(iPoint, iDim+1) = 0.0;
nodes->SetVel_ResTruncError_Zero(iPoint);

/*--- Get transport coefficients ---*/
Expand All @@ -708,8 +705,7 @@ void CNSSolver::BC_Isothermal_Wall_Generic(CGeometry *geometry, CSolver **solver
if (cht_mode) {
Twall = GetCHTWallTemperature(config, val_marker, iVertex, dist_ij,
thermal_conductivity, There, Temperature_Ref);
}
else if (config->GetMarker_All_PyCustom(val_marker)) {
} else if (config->GetMarker_All_PyCustom(val_marker)) {
Twall = geometry->GetCustomBoundaryTemperature(val_marker, iVertex) / Temperature_Ref;
}

Expand Down
34 changes: 17 additions & 17 deletions TestCases/hybrid_regression.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -103,32 +103,32 @@ def main():
flatplate.cfg_dir = "navierstokes/flatplate"
flatplate.cfg_file = "lam_flatplate.cfg"
flatplate.test_iter = 100
flatplate.test_vals = [-7.679131, -2.206953, 0.001084, 0.036233, 2.361500, -2.325300, -1.984700, -1.984700]
flatplate.test_vals = [-7.679131, -2.206953, 0.001084, 0.036233, 2.361500, -2.325300, 0, 0]
test_list.append(flatplate)

# Laminar cylinder (steady)
cylinder = TestCase('cylinder')
cylinder.cfg_dir = "navierstokes/cylinder"
cylinder.cfg_file = "lam_cylinder.cfg"
cylinder.test_iter = 25
cylinder.test_vals = [-8.266513, -2.783904, -0.019899, 1.615668, -0.010207]
cylinder.test_vals = [-8.266513, -2.783904, -0.019899, 1.615668, 0]
test_list.append(cylinder)

# Laminar cylinder (low Mach correction)
cylinder_lowmach = TestCase('cylinder_lowmach')
cylinder_lowmach.cfg_dir = "navierstokes/cylinder"
cylinder_lowmach.cfg_file = "cylinder_lowmach.cfg"
cylinder_lowmach.test_iter = 25
cylinder_lowmach.test_vals = [-6.830996, -1.368850, -0.143956, 73.963354, 0.002457]
cylinder_lowmach.test_vals_aarch64 = [-6.830996, -1.368850, -0.143956, 73.963354, 0.002457]
cylinder_lowmach.test_vals = [-6.830996, -1.368850, -0.143956, 73.963354, 0]
cylinder_lowmach.test_vals_aarch64 = [-6.830996, -1.368850, -0.143956, 73.963354, 0]
test_list.append(cylinder_lowmach)

# 2D Poiseuille flow (body force driven with periodic inlet / outlet)
poiseuille = TestCase('poiseuille')
poiseuille.cfg_dir = "navierstokes/poiseuille"
poiseuille.cfg_file = "lam_poiseuille.cfg"
poiseuille.test_iter = 10
poiseuille.test_vals = [-5.046131, 0.652984, 0.008355, 13.735818, -2.142500]
poiseuille.test_vals = [-5.046131, 0.652984, 0.008355, 13.735818, 0]
test_list.append(poiseuille)

# 2D Poiseuille flow (inlet profile file)
Expand Down Expand Up @@ -157,15 +157,15 @@ def main():
rae2822_sa.cfg_dir = "rans/rae2822"
rae2822_sa.cfg_file = "turb_SA_RAE2822.cfg"
rae2822_sa.test_iter = 20
rae2822_sa.test_vals = [-2.020123, -5.269339, 0.807147, 0.060499, -80603.000000]
rae2822_sa.test_vals = [-2.020123, -5.269339, 0.807147, 0.060499, 0]
test_list.append(rae2822_sa)

# RAE2822 SST
rae2822_sst = TestCase('rae2822_sst')
rae2822_sst.cfg_dir = "rans/rae2822"
rae2822_sst.cfg_file = "turb_SST_RAE2822.cfg"
rae2822_sst.test_iter = 20
rae2822_sst.test_vals = [-0.510363, 4.872736, 0.815617, 0.060920, -73391.000000]
rae2822_sst.test_vals = [-0.510363, 4.872736, 0.815617, 0.060920, 0]
test_list.append(rae2822_sst)

# RAE2822 SST_SUST
Expand All @@ -189,24 +189,24 @@ def main():
turb_oneram6.cfg_dir = "rans/oneram6"
turb_oneram6.cfg_file = "turb_ONERAM6.cfg"
turb_oneram6.test_iter = 10
turb_oneram6.test_vals = [-2.408523, -6.662833, 0.238333, 0.158910, -52718]
turb_oneram6.test_vals = [-2.408523, -6.662833, 0.238333, 0.158910, 0]
test_list.append(turb_oneram6)

# NACA0012 (SA, FUN3D finest grid results: CL=1.0983, CD=0.01242)
turb_naca0012_sa = TestCase('turb_naca0012_sa')
turb_naca0012_sa.cfg_dir = "rans/naca0012"
turb_naca0012_sa.cfg_file = "turb_NACA0012_sa.cfg"
turb_naca0012_sa.test_iter = 5
turb_naca0012_sa.test_vals = [-12.098325, -14.149988, 1.057665, 0.022971, 20.000000, -2.292707, 0.000000, -12.068169, -44.871000]
turb_naca0012_sa.test_vals_aarch64 = [-12.098325, -14.149988, 1.057665, 0.022971, 20.000000, -2.292707, 0.000000, -12.068169, -44.871000]
turb_naca0012_sa.test_vals = [-12.098325, -14.149988, 1.057665, 0.022971, 20.000000, -2.292707, 0.000000, -12.068169, 0]
turb_naca0012_sa.test_vals_aarch64 = [-12.098325, -14.149988, 1.057665, 0.022971, 20.000000, -2.292707, 0.000000, -12.068169, 0]
test_list.append(turb_naca0012_sa)

# NACA0012 (SST, FUN3D finest grid results: CL=1.0840, CD=0.01253)
turb_naca0012_sst = TestCase('turb_naca0012_sst')
turb_naca0012_sst.cfg_dir = "rans/naca0012"
turb_naca0012_sst.cfg_file = "turb_NACA0012_sst.cfg"
turb_naca0012_sst.test_iter = 10
turb_naca0012_sst.test_vals = [-12.105781, -15.277738, -6.210248, 1.049757, 0.019249, -2.807857, -38.976000]
turb_naca0012_sst.test_vals = [-12.105781, -15.277738, -6.210248, 1.049757, 0.019249, -2.807857, 0]
test_list.append(turb_naca0012_sst)

# NACA0012 (SST_SUST, FUN3D finest grid results: CL=1.0840, CD=0.01253)
Expand Down Expand Up @@ -250,7 +250,7 @@ def main():
axi_rans_air_nozzle_restart.cfg_dir = "axisymmetric_rans/air_nozzle"
axi_rans_air_nozzle_restart.cfg_file = "air_nozzle_restart.cfg"
axi_rans_air_nozzle_restart.test_iter = 10
axi_rans_air_nozzle_restart.test_vals = [-12.070954, -7.407644, -8.698118, -4.008751, -3572.100000]
axi_rans_air_nozzle_restart.test_vals = [-12.070954, -7.407644, -8.698118, -4.008751, 0]
test_list.append(axi_rans_air_nozzle_restart)

#################################
Expand Down Expand Up @@ -381,8 +381,8 @@ def main():
inc_poly_cylinder.cfg_dir = "incomp_navierstokes/cylinder"
inc_poly_cylinder.cfg_file = "poly_cylinder.cfg"
inc_poly_cylinder.test_iter = 20
inc_poly_cylinder.test_vals = [-7.851512, -2.093420, 0.029974, 1.921595, -175.300000]
inc_poly_cylinder.test_vals_aarch64 = [-7.851510, -2.093419, 0.029974, 1.921595, -175.300000]
inc_poly_cylinder.test_vals = [-7.827942, -2.061513, 0.029525, 1.953498, -174.780000]
inc_poly_cylinder.test_vals_aarch64 = [-7.827942, -2.061513, 0.029525, 1.953498, -174.780000]
test_list.append(inc_poly_cylinder)

# X-coarse laminar bend as a mixed element CGNS test
Expand Down Expand Up @@ -452,8 +452,8 @@ def main():
square_cylinder.cfg_dir = "unsteady/square_cylinder"
square_cylinder.cfg_file = "turb_square.cfg"
square_cylinder.test_iter = 3
square_cylinder.test_vals = [-2.560839, -1.173497, 0.061188, 1.399403, 2.220575, 1.399351, 2.218781, -0.584690]
square_cylinder.test_vals_aarch64 = [-2.557902, -1.173574, 0.058050, 1.399794, 2.220402, 1.399748, 2.218604, -0.453270]
square_cylinder.test_vals = [-2.560839, -1.173497, 0.061188, 1.399403, 2.220575, 1.399351, 2.218781, 0]
square_cylinder.test_vals_aarch64 = [-2.557902, -1.173574, 0.058050, 1.399794, 2.220402, 1.399748, 2.218604, 0]
square_cylinder.unsteady = True
test_list.append(square_cylinder)

Expand Down Expand Up @@ -503,7 +503,7 @@ def main():
ddes_flatplate.cfg_dir = "ddes/flatplate"
ddes_flatplate.cfg_file = "ddes_flatplate.cfg"
ddes_flatplate.test_iter = 10
ddes_flatplate.test_vals = [-2.714786, -5.882652, -0.215041, 0.023758, -617.470000]
ddes_flatplate.test_vals = [-2.714786, -5.882652, -0.215041, 0.023758, 0]
ddes_flatplate.unsteady = True
test_list.append(ddes_flatplate)

Expand Down
4 changes: 2 additions & 2 deletions TestCases/hybrid_regression_AD.py
View file
Open in desktop
Loading