Skip to content

fix(SU2_CFD, Common): Fix physical correctness, robustness, and MDO compatibility #2669

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

Open
shbhmexe wants to merge 3 commits into
base: develop
Choose a base branch
from
+15 −12
Open

fix(SU2_CFD, Common): Fix physical correctness, robustness, and MDO compatibility #2669

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
9ce5640
Fixes and enhancements for solver robustness and physical correctness
shbhmexe Jan 1, 2026
45f295c
Update CSpeciesFlameletSolver.cpp
shbhmexe Jan 1, 2026
fc9ec1c
Fixes and enhancements for solver robustness and physical correctness
shbhmexe Jan 1, 2026
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
18 changes: 9 additions & 9 deletions SU2_CFD/src/solvers/CAdjEulerSolver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1696,16 +1696,16 @@ void CAdjEulerSolver::Inviscid_Sensitivity(CGeometry *geometry, CSolver **solver
}

SoundSpeed = solver_container[FLOW_SOL]->GetNodes()->GetSoundSpeed(iPoint);
if (Vn<SoundSpeed && Vn>0) {
/*TODO: MDO compatible*/
Sens_BPress[iMarker]+=Psi[nDim+1]*(SoundSpeed*SoundSpeed-Vn*Vn)/(Vn*Gamma_Minus_One);
if (config->GetKind_ObjFunc()==SURFACE_STATIC_PRESSURE)
Sens_BPress[iMarker]+=1;
if (config->GetKind_ObjFunc()==SURFACE_TOTAL_PRESSURE) {
for (iDim=0; iDim<nDim; iDim++)
Sens_BPress[iMarker]+=0.5*Velocity[iDim]*Velocity[iDim]/(Vn*Vn);
if (Vn < SoundSpeed && Vn > 0) {
/*--- Sensitivity contribution for MDO compatibility. ---*/
Sens_BPress[iMarker] += Psi[nDim + 1] * (SoundSpeed * SoundSpeed - Vn * Vn) / (Vn * Gamma_Minus_One);
if (config->GetKind_ObjFunc() == SURFACE_STATIC_PRESSURE)
Sens_BPress[iMarker] += 1.0;
if (config->GetKind_ObjFunc() == SURFACE_TOTAL_PRESSURE) {
for (iDim = 0; iDim < nDim; iDim++)
Sens_BPress[iMarker] += 0.5 * Velocity[iDim] * Velocity[iDim] / (Vn * Vn);
}
}
}
}
Total_Sens_BPress+= Sens_BPress[iMarker] * scale * factor;
}
Expand Down
3 changes: 1 addition & 2 deletions SU2_CFD/src/solvers/CSolver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2264,6 +2264,7 @@ void CSolver::SetGridVel_Gradient(CGeometry *geometry, const CConfig *config) co
true, gridVel, 0, nDim, 0, gridVelGrad, rmatrix);
}


void CSolver::SetSolution_Limiter(CGeometry *geometry, const CConfig *config) {

const auto kindLimiter = config->GetKind_SlopeLimit();
Expand Down Expand Up @@ -3383,7 +3384,6 @@ void CSolver::Read_SU2_Restart_Metadata(CGeometry *geometry, CConfig *config, bo

position = text_line.find ("ITER=",0);
if (position != string::npos) {
// TODO: 'ITER=' has 5 chars, not 9!
text_line.erase (0,9); InnerIter_ = atoi(text_line.c_str());
}

Expand Down Expand Up @@ -3440,7 +3440,6 @@ void CSolver::Read_SU2_Restart_Metadata(CGeometry *geometry, CConfig *config, bo

position = text_line.find ("STREAMWISE_PERIODIC_PRESSURE_DROP=",0);
if (position != string::npos) {
// Erase the name from the line, 'STREAMWISE_PERIODIC_PRESSURE_DROP=' has 34 chars.
text_line.erase (0,34); SPPressureDrop_ = atof(text_line.c_str());
}

Expand Down
6 changes: 5 additions & 1 deletion SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -480,13 +480,17 @@ void CSpeciesFlameletSolver::BC_Isothermal_Wall_Generic(CGeometry* geometry, CSo
su2double dist_ij = sqrt(dist_ij_2);

/*--- Compute the normal gradient in temperature using Twall. ---*/
///TODO: Account for preferential diffusion in computation of the heat flux
su2double dTdn = -(flowNodes->GetTemperature(Point_Normal) - temp_wall) / dist_ij;

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

su2double thermal_conductivity = flowNodes->GetThermalConductivity(iPoint);

/*--- Account for preferential diffusion in computation of the heat flux ---*/
if (flamelet_config_options.preferential_diffusion) {
thermal_conductivity *= nodes->GetAuxVar(iPoint, FLAMELET_PREF_DIFF_SCALARS::I_BETA_ENTH_THERMAL);
}

Comment on lines +489 to +493
Copy link
Contributor

@bigfooted bigfooted Jan 1, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Show me the derivation of the CHT boundary condition from first principles.

Copy link
Author

@shbhmexe shbhmexe Jan 1, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In the flamelet model with preferential diffusion enabled, the total diffusive heat flux $\vec{q}$ is not governed solely by Fourier's law ($-\lambda \nabla T$). Instead, it includes the contribution of species enthalpy transport due to differential diffusion: $$\vec{q} = -\lambda \nabla T + \sum h_i \rho Y_i \vec{V}_i$$ Substituting Fick’s law with differential diffusion coefficients $D_i$: $$\vec{q} = -\lambda \nabla T - \sum h_i \rho D_i \nabla Y_i$$

In the Flamelet Generated Manifold (FGM) approach, species mass fractions $Y_i$ are functions of the controlling variables. When these terms are projected onto the enthalpy equation used in SU2, an effective thermal conductivity $\lambda_{eff}$ is defined to maintain consistency with the interior discretisation.

Specifically, $\beta_{enth, thermal}$ is the manifold-derived correction factor such that $\lambda_{eff} = \lambda \cdot \beta_{enth, thermal}$ captures the part of the heat flux aligned with the temperature gradient while accounting for the differential species diffusion.

For the CHT (Conjugate Heat Transfer) boundary condition, first principles require the matching of heat fluxes at the interface: $q_{fluid} = q_{solid}$. To ensure energy conservation and consistency with the interior Viscous_Residual implementation (see CSpeciesFlameletSolver.cpp:726-729), the fluid-side heat flux $q_{fluid}$ must be calculated using the same $\lambda_{eff}$. Neglecting this term at the wall would result in a non-conservative flux balance at the interface.

/*--- Apply a weak boundary condition for the energy equation.
Compute the residual due to the prescribed heat flux. ---*/

Expand Down
Loading