Align the SST model with the NASA TMR definitions and fix the boundary energy (#1851)

[rois1995]

Proposed Changes

This PR aligns the SST implementation with the definitions of the NASA Turbulence Modeling Resource (TMR, SST page) and Menter's papers, and fixes the turbulent kinetic energy imposed at the flow boundaries (#1851). The default (modified, "m") versions keep their interior formulation; results change mainly because of the boundary fix.

1. SST versions as defined by TMR

• "m" versions (V1994m, V2003m, the default): TMR defines them as ignoring ⅔ρk in the stress tensor and using the production P = μ_t·S². This is what develop already does; unchanged.
• Standard versions (V1994, V2003): they were disabled in develop. They now include ⅔ρk in the stress tensor (momentum and energy, scalar and SIMD viscous fluxes) and use the exact production (including −⅔ρk·div u and −⅔μ_t·(div u)²), as TMR's standard SST.
• k in the energy / equation of state is kept for all versions, as in develop (TMR does not tie it to "m").

2. Model details

• F1 cross-diffusion clip: develop used 10^(−production limiter constant) (CTurbSSTVariable.cpp). The values happened to be right (20 → 1e-20 for SST-1994, 10 → 1e-10 for SST-2003) but were tied to the production limiter. Now the fixed values of each version (Menter 1994 eq. A10; TMR "uses 10⁻¹⁰ rather than 10⁻²⁰" for 2003).
• SST-1994 ω production for the standard version: γ/ν_t·P (Menter 1994 eq. A2), with a guard for μ_t = 0.
• Jacobian of the standard versions: the −⅔·div u term is only added where it increases the diagonal (the residual is unchanged).
• SST-sust ambient values: develop used the free-stream values from the turbulence intensity and viscosity ratio. TMR's SST-sust uses k_amb = 10⁻⁶·U² (Tu = 0.08165 %) and ω_amb = 5·U/L (Spalart & Rumsey, AIAA J. 45(10), 2007). These are now the defaults, with L = REYNOLDS_LENGTH; new options SST_SUST_TKE_AMB and SST_SUST_OMEGA_AMB override them.

3. Boundary energy (#1851)

The flow boundary conditions built the boundary total energy with the free-stream k (Energy += GetTke_Inf()), while the turbulence k at the boundary nodes decays immediately (SST destruction). The flow then received too much energy: the static temperature at inlets was too high and the Mach number too low, especially at high Mach number and turbulence intensity. As proposed in #1851, all flow boundary states that include k (far field, inlet, supersonic inlet, Riemann, outlet, engine, actuator disk) now use the k of the boundary node from the turbulence solver.

Reproduced at the conditions of #1851 (flat plate, M = 6.1, T = 800 K, Re = 4.9·10⁶ 1/m, T_wall = 300 K, Tu = 5 %): inlet temperature 841 K and Mach 5.95 with develop, 793 K and 6.12 with this PR (imposed: 800 K, 6.1).

4. Options and documentation

• SST_OPTIONS: V1994 and V2003 enabled; TMRBC (far-field and inlet ω = 10·U/L_DOMAIN, the original far-field prescription reported by TMR); WALL_OMEGA_LIMIT (clips the ω wall value 60ν/(β₁d²) to the upper limit of ω; off by default).
• config_template.cfg: the SST versions, the new options, and the TMR range for the far-field μ_t/μ (10⁻⁵ to 10⁻², next to FREESTREAM_TURB2LAMVISCRATIO; the default of 10 is unchanged).

Test cases

Reference values updated for the SST tests (serial, parallel, hybrid, AD, tutorials, V&V). Most changes come from the boundary fix; restart cases change most because their solutions were converged with the old boundary k; SST-sust cases change because of the TMR ambient values. A few values could not be reproduced on my machine even with develop, and three tests changed on both sides in the last merge of develop (rae2822_sst_sust, turb_naca0012_sst_restart_mg, bars_SST_2D); these will be taken from the CI logs.

Related Work

• Fixes Turbulent Kinetic Energy(TKE) on energy equation in SST model. #1851.
• The "m" versions were introduced in Some questions about Menter's SST model on SU2 (might be the bugs) #1483.

PR Checklist

• I am submitting my contribution to the develop branch.
• My contribution generates no new compiler warnings (try with --warnlevel=3 when using meson).
• My contribution is commented and consistent with SU2 style (https://su2code.github.io/docs_v7/Style-Guide/).
• I used the pre-commit hook to prevent dirty commits and used pre-commit run --all to format old commits.
• I have added a test case that demonstrates my contribution, if necessary.
• I have updated appropriate documentation (Tutorials, Docs Page, config_template.cpp), if necessary.

🤖 Generated with Claude Code

[bigfooted]

Was this intentionally commented, or...? If it is used, I guess this should go to line 872

Suggested change

	nPrandtl_Lam, /*!< \brief Number of species
	addDoubleOption("FREESTREAM_TURB2LAMVISCRATIO", TurbIntensityAndViscRatioFreeStream[1], 10.0); Prandtl number. */
	nPrandtl_Lam, /*!< \brief Number of species laminar Prandtl number. */
	addDoubleOption("FREESTREAM_TURB2LAMVISCRATIO", TurbIntensityAndViscRatioFreeStream[1], 10.0); /*!<\brief Freestream mu_turb to mu_lam viscosity ratio */

[rois1995]

This is just a wrong copy-paste from me. It should not be there in the first place.

[bigfooted]

Great contribution, Thanks @rois1995 !

[pcarruscag]

If you are looking into robustness aspects too you should get in touch with @emaberman and @YairMO, seems like they have some good ideas and between the free time of 3 people a lot more can get done :)

[YairMO]

Hi,

Regarding the cross-diffusion term (CD) that appears in Omega source term (residual). The SST model (1994/2003) is a high-Reynolds-number model. Namely, It can not predict correctly the sub-layer region (especially the correct profile of the TKE). Therefore, only a positive contribution is required. Moreover, since the SST model was design as a k-w and k-epsilon blending, the CD term "belongs" only to the k-epsilon "branch", that is why the CD term include the factor "1-F1". However, it may happen, that the factor "1-F1" is not a 100% safe guarantee. It may happen that "1-F1" is not zero in region where the CD term is negative (this happen due numerical errors). To avoid such a situation, it is a good idea to clip the CD term with zero. Otherwise, severe numerical robustness issues may rise.
Yes, it is different from Menter publications, but I think that clipping the CD term with zero is completely inline with Menter original idea (that is why, I think, he was including the factor "1-F1". But again the 1-F1 factor is not 100% percent "safe").

[YairMO]

Hi,

The use of an Omega production limiter (about the cross-diffusion term) is correct for low-Reynolds-number (LRN) models (the approach described by Peng et al. is very naive; there are other more rigorous treatments). For high-Reynolds-number (HRN) models, the clipping should be zero, keeping the cross-diffusion term positive; thus, the current implementation is correct.

Indeed, it is not exactly as it appears in Menter's original publication. The factor (1-F1) aimed to promise that the cross-diffusion term will be activated only outside the boundary layer, where it is positive (the cross-diffusion term switches its sign deep inside the boundary layer). This was also recognized by Peng et al. (first paragraph above Eq. 17). However, it may happen that the factor (1-F1)=1 where the cross-diffusion term is negative. Usually, it may happen at the wake, very near the airfoil trailing edge, where the upper and lower boundary layers merge. It is due to the imperfection of the F1 function.

To summarize, the current implementation is correct, and it is perfect for HRN models.

[YairMO]

For the sake of clarity, "current implementation" refers to the current treatment of the production code

[emaberman]

What YairMO is saying, is that allowing negative cross diffusion values is incorrect for high Reynolds models and should not be an option, this is a fix used for low Reynolds models only

[rois1995]

Hi @YairMO, Hi @emaberman ,

thank you very much for your comments. I haven't found any suggestion in literature to clip to only positive values the cross-diffusion term in the w-equation. I understand that it might be more robust, but it is not the standard implementation of the SST model, which is the first thing that we need to achieve. Only then we can build on top of that to improve the robustness of SU2.

Nevertheless, I tried the SWBLI test case and I compared the results across 6 different combinations:

1- develop branch, no changes
2- develop branch, changes to Supersonic_inlet profile as suggested in #1851
3- my branch, with original CDkw implementation (should give exactly the same result as develop+modified BC)
3- my branch, with original CDkw implementation and using boundary conditions from TMR
4- my branch, with original CDkw implementation and using your suggestions for lower limits for k and w.
5- my branch, allowing negative values of CDkw
6- my branch, allowing negative values of CDkw and using boundary conditions from TMR
7- my branch, allowing negative values of CDkw and using your suggestions for lower limits for k and w.
8- my branch, allowing negative values of CDkw, using boundary conditions from TMR and using your suggestions for lower limits for k and w.

When my branch is used, then the changes to the supersonic inlet BC are already in place.

I haven't achieved convergence with 1, 2 and 3. More precisely, 1 diverged right away (after 30 iterations), while 2 and 3 gave "FGMRES - Orthogonalization Failed" after 900ish iterations.

Here you can see the residuals for the different combinations.

Unfortunately I will be busy with the AIAA Conference next week, thus I don't know how much I will be able to work on this. The next test case will be the 2D airfoil near-wake from TMR.

[YairMO]

Hi rois1995,

First of all, enjoy your time in Las Vegas. Any paper that you are presenting?

As for our discussion about the cross-diffusion term, I've emailed the "source" (Menter). I believe he will make it clear.
It may be that he will be able to answer only in a while ...

[rois1995]

Sorry for the long pause on this PR. I have brought it up to date with develop and reduced it to what its title says: aligning the SST implementation with the NASA TMR definitions.

• Following @pcarruscag's request, the default ("m") versions keep their interior formulation (k stays in the energy; TMR's "m" only concerns ⅔ρk in the stress tensor and the production). The standard V1994/V2003 versions are enabled with ⅔ρk and the exact production.
• The boundary energy now uses the k of the boundary node from the turbulence solver, as proposed in Turbulent Kinetic Energy(TKE) on energy equation in SST model. #1851; I reproduced the issue's flat plate case and the inlet temperature and Mach number are now correct.
• SST-sust uses the TMR ambient values by default, with options to change them.
• Everything that was not TMR alignment (a user-defined production limiter, debug outputs, the Suzen-Hoffmann correction) has been removed from this PR and will come separately.

The description is updated with the details.