Conversation
- Added production terms to SST
- Added user-defined production limiter constant for TKE
rois1995
changed the title
- Updated boundary conditions as in TMR page
| nPrandtl_Lam, /*!< \brief Number of species | ||
| addDoubleOption("FREESTREAM_TURB2LAMVISCRATIO", TurbIntensityAndViscRatioFreeStream[1], 10.0); Prandtl number. */ |
There was a problem hiding this comment.
Was this intentionally commented, or...? If it is used, I guess this should go to line 872
| 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 */ |
There was a problem hiding this comment.
This is just a wrong copy-paste from me. It should not be there in the first place.
|
Great contribution, Thanks @rois1995 ! |
|
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 :) |
|
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. |
- Given option for cross diffusion limiting in W residual
|
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. |
|
For the sake of clarity, "current implementation" refers to the current treatment of the production code |
|
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 |
|
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 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. |
|
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. |
- Cross diffusion: develop computes the full term in the viscous fluxes, so it is removed from the source term; the Suzen-Hoffmann omega correction is kept as a separate source term, only when that option is used. - Wall BC: develop's rough-wall models, with the upper limit of omega of this branch on smooth walls. Inlet BC: the TMR-based option (newBC) moved to the code shared by compressible and incompressible inlets. - k in the flow equations: this branch includes it only for the non-modified SST versions. The 2/3 rho k term of the stress tensor is kept in develop's SIMD viscous flux; the tke in the SIMD convective fluxes is not ported to develop's refactored upwind schemes. With UQ the 2/3 rho k term is not added twice (it is already in the perturbed Reynolds stress). - Free-stream values: develop's protection against division by zero. Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]> Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6
- Remove the GetStrainMag overloads added to CVariable: one duplicated the one of develop, and both had no return statement (undefined behavior if called); CFlowVariable provides them for the flow variables. - Remove the using declaration of cp, no longer a member of the SIMD viscous flux base class. Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]> Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6
ProdDistr was only set when the wall distance is larger than 1e-10, so the points on the wall kept the values of the previous point. Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]> Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6
…osity The pressure-dilatation term of the omega equation, -(1-F1) p'd'/nu_t, was divided by the dynamic eddy viscosity, which has the wrong units (a factor of density is missing). Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]> Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6
The modified (m) SST versions ignore k in the total energy, but the free stream energy of the incompressible, NEMO and FEM-DG solvers and the energy of the engine inflow/exhaust and actuator disk boundaries still added it. Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]> Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6
LoadSurfaceData, shared by the compressible, incompressible and NEMO outputs, sets NEIGHBORDIST, but it was only registered by the compressible one, so incompressible and NEMO runs with viscous walls stopped with "Cannot find output field with name NEIGHBORDIST". Also fix its description. Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]> Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6
SST_SUST_TKE_AMB and SST_SUST_OMEGA_AMB set the ambient k and omega of the sustaining terms, which are also the free-stream values. By default (<= 0) they are those of Spalart and Rumsey (2007) used by the NASA TMR SST-sust, k = 1e-6 U^2 and omega = 5 U / L with L = REYNOLDS_LENGTH, as before. The values were hard-coded in the four flow solvers and in the SST solver; they are now computed by CConfig, non-dimensionalized with the reference velocity and omega, and the SST solver takes the free-stream values set by the flow solver. The default results are unchanged. Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]> Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6
The production of omega, gamma/nu_t * P, was computed as P / mu_t, which is 0/0 where the eddy viscosity vanishes (e.g. k = 0). Expand P / mu_t term by term: the strain (or vorticity) squared, the full-production divergence term and -2/3 rho k div(u) / mu_t, which is bounded because mu_t is proportional to k. The modified (m) versions are again identical to develop. Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]> Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6
The derivative of -2/3 rho k div(u) was added to the diagonal of the k equation, which reduces it in compression (div(u) < 0) and can make the implicit system unstable. Add it only where it increases the diagonal; the residual is unchanged, so the converged solution is the same. Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]> Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6
This PR only aligns SST with the NASA TMR. The COMPRESSIBILITY-SH option and its terms are moved to the local branch feature_SST_SuzenHoffmann (based on develop). Also set the codi and medi submodules back to those of develop. Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]> Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6
…t boundaries The NASA TMR "m" versions only ignore 2/3 rho k in tau_ij and approximate P = mu_t S^2; they say nothing about k in the energy, which develop includes for all SST versions. Restore that (free streams, boundary conditions and primitive variables of all flow solvers), keeping the 2/3 rho k term in tau_ij only for the standard versions. The k added to the energy of a boundary state is now the one the turbulence boundary condition imposes: - far field: free-stream k at inflow, interior k at outflow; - inlet and supersonic inlet: the k of the turbulence inlet (profile, or intensity and viscosity ratio of the marker) for the same inlet state, computed by CTurbSSTSolver::GetInletTke, shared with its BC_Inlet; - Riemann: the inlet k of the exterior state at inflows, the interior k at STATIC_PRESSURE; k is also removed from the static energy of the interior and boundary states (it was not); - outlet: interior k (the turbulence outlet is Neumann). Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]> Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6
…LPROD NASA TMR: the standard SST uses the exact production P = tau_ij du_i/dx_j = mu_t (S^2 - 2/3 div(u)^2) - 2/3 rho k div(u); the "m" versions use P = mu_t S^2. With the vorticity and Kato-Launder forms only the -2/3 rho k div(u) term is added. The FULLPROD option, which selected the -2/3 mu_t div(u)^2 term, is no longer needed. Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]> Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6
Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]> Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6
SST_OPTIONS= WALL_OMEGA_LIMIT clips Menter's wall value 60 nu / (beta_1 d^2) to the upper limit of omega used in the rest of the domain. By default the wall value is not clipped, as in develop. Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]> Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6
… branch This PR only aligns SST with the NASA TMR. Moved to the local branch feature_SST_extras: the user-defined production limiter (PRODLIM, PROD_LIM_CONST), the debug volume outputs (production and destruction of k and omega, production limiter flag, CDkw, F1, F2, wall distance, strain magnitude, velocity gradient, NEIGHBORDIST) and their storage, and unused code in CFlowTractionInterface. Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]> Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6
Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]> Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6
Write the 2/3 rho k term as a separate addend, so that without it (all models but the standard SST versions) the result is bit-identical to develop. Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]> Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6
The standard (non-m) SST versions add 2/3 rho k to the normal stresses after the stress tensor of develop, which is otherwise unchanged; the overload of stressTensor with density and k is no longer needed. Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]> Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6
…tes (#1851) As diagnosed in #1851, adding the free-stream k to the total energy of a boundary state gives the flow too much energy, because k decays right after the boundary (SST destruction): at inflows the static temperature is too high and the Mach number too low, most visibly at high Mach number and high turbulence intensity. Following the fix proposed there, read k from the turbulence solver at the boundary node in every flow boundary condition that builds a total energy: far field, inlet (total conditions and mass flow), supersonic inlet, Riemann (also removed from the static energy of the boundary state), outlet, engine inflow and exhaust, actuator disks. This replaces the previous choice of the free-stream or inlet k at inflows, which had the same problem; the shared GetInletTke is no longer needed. Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]> Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6
Serial, parallel, hybrid, AD, tutorials and V&V SST tests. The values change because of the k of the boundary node in the boundary energy (#1851), the NASA TMR SST-sust free stream, the exact production of the standard versions, and UQ no longer counting 2/3 rho k twice. Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]> Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6
rois1995
changed the title
|
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.
The description is updated with the details. |
The reference values of rae2822_sst_sust, turb_naca0012_sst_restart_mg and bars_SST_2D, changed on both sides, are taken from develop for now and will be updated from the CI. Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]> Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6
It sets the far-field values of the NASA TMR, so the name says where they come from. Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]> Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6
…ean up whitespace The 2/3 rho k term is selected inside SetStressTensor, so the callers are as in develop. Remove stray blank lines and extra spaces, and use beta_star in the omega Jacobian. Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]> Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6
SST cases change with the TMR alignment and the boundary k (#1851). The laminar, FSI and SA cases (poiseuille_profile, fsi_cht, pywrapper_unsteadyFSI, turb_naca0012_sa, turb_flatplate) move only in the last digits, from the rounding of the optimized build. Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]> Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6
This branch has not been deployed
6 participants
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
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.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.2. Model details
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).REYNOLDS_LENGTH; new optionsSST_SUST_TKE_AMBandSST_SUST_OMEGA_AMBoverride 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:V1994andV2003enabled;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 toFREESTREAM_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
PR Checklist
pre-commit run --allto format old commits.🤖 Generated with Claude Code