-
Notifications
You must be signed in to change notification settings - Fork 842
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.
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
[WIP] Isothermal streamwise periodicity #1799
base: develop
Are you sure you want to change the base?
Changes from 28 commits
514fd95
449c16a
97820a2
fa1f38a
15a5eb2
fd9afe2
988b5a9
978a87f
26f3339
9799d82
c6e1277
0425778
3eaa05e
e383836
5651b3d
f0a817c
71c3d95
2f4888f
a777bd6
19e3f41
757c6b2
8ef2141
02f61e5
c931b98
03fe2ad
78b0ab0
0f4bb42
50360de
d67666b
99e39e8
bec2016
5e4b0a8
c0ae9fd
8c5f3f4
b1f7acb
2f95ece
fc0caaf
5ef89a5
abfa537
688d8ea
6f80a62
0f2e642
93532b2
52fc485
804d89d
ce08528
2b6a2ea
33df3db
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -120,7 +120,7 @@ void CIncNSSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container | |
void CIncNSSolver::GetStreamwise_Periodic_Properties(const CGeometry *geometry, | ||
CConfig *config, | ||
const unsigned short iMesh) { | ||
|
||
const bool turbulent = (config->GetKind_Turb_Model() != TURB_MODEL::NONE); | ||
/*---------------------------------------------------------------------------------------------*/ | ||
// 1. Evaluate massflow, area avg density & Temperature and Area at streamwise periodic outlet. | ||
// 2. Update delta_p is target massflow is chosen. | ||
|
@@ -191,15 +191,16 @@ void CIncNSSolver::GetStreamwise_Periodic_Properties(const CGeometry *geometry, | |
SPvals.Streamwise_Periodic_BoundaryArea = Area_Global; | ||
SPvals.Streamwise_Periodic_AvgDensity = Average_Density_Global; | ||
|
||
su2double HeatFlow_Local = 0.0, HeatFlow_Global = 0.0; | ||
su2double dTdn_Local = 0.0, dTdn_Global = 0.0; | ||
|
||
if (config->GetEnergy_Equation()) { | ||
/*---------------------------------------------------------------------------------------------*/ | ||
/*--- 3. Compute the integrated Heatflow [W] for the energy equation source term, heatflux ---*/ | ||
/*--- boundary term and recovered Temperature. The computation is not completely clear. ---*/ | ||
/*--- Here the Heatflux from all Boundary markers in the config-file is used. ---*/ | ||
/*---------------------------------------------------------------------------------------------*/ | ||
|
||
su2double HeatFlow_Local = 0.0, HeatFlow_Global = 0.0; | ||
|
||
/*--- Loop over all heatflux Markers ---*/ | ||
for (auto iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { | ||
|
||
|
@@ -222,13 +223,91 @@ void CIncNSSolver::GetStreamwise_Periodic_Properties(const CGeometry *geometry, | |
HeatFlow_Local += FaceArea * (-1.0) * Wall_HeatFlux/config->GetHeat_Flux_Ref();; | ||
} // loop Vertices | ||
} // loop Heatflux marker | ||
|
||
if (config->GetMarker_All_KindBC(iMarker) == ISOTHERMAL) { | ||
/*--- Identify the boundary by string name and retrive ISOTHERMAL from config ---*/ | ||
|
||
for (auto iVertex = 0ul; iVertex < geometry->nVertex[iMarker]; iVertex++) { | ||
|
||
const auto iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); | ||
|
||
if (!geometry->nodes->GetDomain(iPoint)) continue; | ||
|
||
/*--- Compute wall heat flux (normal to the wall) based on computed temperature gradient ---*/ | ||
const auto AreaNormal = geometry->vertex[iMarker][iVertex]->GetNormal(); | ||
|
||
su2double GradT[MAXNDIM] = {0.0,0.0,0.0}; | ||
for (auto iDim = 0u; iDim < nDim; iDim++) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. So I have to revert this as GetGradient_Primitive gives a Matrix su2double which does not comply with DotProduct. I saw at other parts of the code for loop being used to construct a variable and then do DotProduct. Do you know a better way to handle this? There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Sorry, slightly wrong suggestion, do this instead please:
|
||
GradT[iDim] = nodes->GetGradient_Primitive(iPoint, prim_idx.Temperature(), iDim); | ||
|
||
dTdn_Local += nodes->GetThermalConductivity(iPoint) * GeometryToolbox::DotProduct(nDim, GradT, AreaNormal); | ||
} // loop Vertices | ||
} // loop Isothermal Marker | ||
} // loop AllMarker | ||
|
||
/*--- MPI Communication sum up integrated Heatflux from all processes ---*/ | ||
SU2_MPI::Allreduce(&HeatFlow_Local, &HeatFlow_Global, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); | ||
SU2_MPI::Allreduce(&dTdn_Local, &dTdn_Global, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); | ||
|
||
su2double Volume_Temp_Local = 0.0, Volume_Temp_Global = 0.0; | ||
su2double Volume_TempS_Local = 0.0, Volume_TempS_Global = 0.0; | ||
su2double Volume_Local = 0.0, Volume_Global = 0.0; | ||
su2double Volume_VTemp_Local = 0.0, Volume_VTemp_Global = 0.0; | ||
su2double turb_b1_coeff_Local = 0.0, turb_b1_coeff_Global = 0.0; | ||
|
||
/*--- Loop over all heatflux Markers ---*/ | ||
for (unsigned long iPoint = 0; iPoint < geometry->GetnPointDomain(); iPoint++) { | ||
|
||
const su2double volume = geometry->nodes->GetVolume(iPoint); | ||
|
||
const su2double Temp = nodes->GetTemperature(iPoint); | ||
|
||
Volume_Local += volume; | ||
|
||
Volume_TempS_Local += volume * Temp; | ||
|
||
Volume_Temp_Local += volume * Temp * nodes->GetThermalConductivity(iPoint); | ||
|
||
// coeff_b1 for turbulence | ||
if (turbulent && (config->GetnMarker_Isothermal() != 0)) | ||
turb_b1_coeff_Local += Temp * nodes->GetAuxVarGradient(iPoint, 0, 0) * config->GetSpecific_Heat_Cp() * volume / config->GetPrandtl_Turb(); | ||
|
||
Volume_VTemp_Local += volume * Temp * nodes->GetVelocity(iPoint, 0) * nodes->GetDensity(iPoint); | ||
|
||
} // points | ||
|
||
/*--- MPI Communication sum up integrated Heatflux from all processes ---*/ | ||
SU2_MPI::Allreduce(&Volume_Temp_Local, &Volume_Temp_Global, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); | ||
SU2_MPI::Allreduce(&Volume_VTemp_Local, &Volume_VTemp_Global, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); | ||
SU2_MPI::Allreduce(&Volume_Local, &Volume_Global, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); | ||
SU2_MPI::Allreduce(&Volume_TempS_Local, &Volume_TempS_Global, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); | ||
if (turbulent && (config->GetnMarker_Isothermal() != 0)) | ||
SU2_MPI::Allreduce(&turb_b1_coeff_Local, &turb_b1_coeff_Global, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); | ||
|
||
/*--- Set the solver variable Integrated Heatflux ---*/ | ||
SPvals.Streamwise_Periodic_IntegratedHeatFlow = HeatFlow_Global; | ||
if (config->GetnMarker_HeatFlux() > 0) | ||
SPvals.Streamwise_Periodic_IntegratedHeatFlow = HeatFlow_Global; | ||
|
||
if (config->GetnMarker_Isothermal() > 0) { | ||
SPvals.Streamwise_Periodic_ThetaScaling = Volume_TempS_Global/Volume_Global; | ||
|
||
for (unsigned long iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) | ||
nodes->SetStreamwise_Periodic_RecoveredTemperature(iPoint, nodes->GetTemperature(iPoint)/SPvals.Streamwise_Periodic_ThetaScaling); | ||
|
||
/*--- Set the solver variable Lambda_L for iso-thermal BCs ---*/ | ||
const su2double b0_coeff = Volume_Temp_Global; | ||
const su2double b1_coeff = Volume_VTemp_Global * config->GetSpecific_Heat_Cp() + turb_b1_coeff_Global; | ||
const su2double b2_coeff = -dTdn_Global; | ||
|
||
/*--- Find the value of Lambda L by solving the quadratic equation ---*/ | ||
const su2double pred_lambda = (- b1_coeff + sqrt(b1_coeff * b1_coeff - 4 * b0_coeff * b2_coeff))/(2 * b0_coeff); | ||
if (!config->GetDiscrete_Adjoint()) | ||
SPvals.Streamwise_Periodic_LambdaL -= 0.01 * (SPvals.Streamwise_Periodic_LambdaL - pred_lambda); | ||
else | ||
SPvals.Streamwise_Periodic_LambdaL = pred_lambda; | ||
|
||
config->SetStreamwise_Periodic_LamdaL(SPvals.Streamwise_Periodic_LambdaL); | ||
} // if isothermal | ||
} // if energy | ||
} | ||
|
||
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Do you see any edge case that might pass this?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Looks ok but @TobiKattmann should know better