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. |
…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) <noreply@anthropic.com> 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) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6
Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com> 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) <noreply@anthropic.com> 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) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6
Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com> 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) <noreply@anthropic.com> 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) <noreply@anthropic.com> 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) <noreply@anthropic.com> 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) <noreply@anthropic.com> 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) <noreply@anthropic.com> 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) <noreply@anthropic.com> 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) <noreply@anthropic.com> 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) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6
DENSITY_VELOCITY copies the total energy of the node, which already contains k, and k of the node was then added again to every Riemann case. The new test axi_rans_air_nozzle_density_velocity (air nozzle with a DENSITY_VELOCITY inlet and 20 % turbulence intensity) detects the double count in the inlet static temperature. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6
… selected The compressibility corrections were stored in SSTParsedOptions.production, the field of the production modifiers, so e.g. VORTICITY or KATO-LAUNDER together with COMPRESSIBILITY-WILCOX or -SARKAR silently lost the modified production. The corrections now use only their own flags (compWilcox, compSarkar), and both are applied and reported. Configurations with only one of the two options are unchanged (turb_flatplate_CC_Wilcox/Sarkar give the same values). Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6
…uxes and axisymmetric sources With SST the total energy (and enthalpy) contains k: - ComputeConsistentExtrapolation and LowMachPrimitiveCorrection (non-ideal gas or low-Mach correction) rebuilt the face enthalpy without k. The k of the cells is now added, as the vectorized reconstruction does. Test: an SST flat plate with standard air and with a van der Waals gas of vanishing a, b (same physics, non-ideal code path) gave CD 0.6 % apart, now 0.04 % as with SA. - The speed of sound of Roe (vectorized and scalar), Turkel, general Roe, HLLC (ideal and general) and AUSM was computed from H including k. The (Roe-averaged) k is now removed from the enthalpy where the speed of sound or the thermodynamic state is computed; the energy flux keeps the total enthalpy. The solver passes the k of the cells (interior faces) and of the boundary node (boundary states). - Axisymmetric source: the pressure is recovered without rho k, the -2/3 rho k of the radial normal stress is outside the viscosity factor and only used where k is part of the stress tensor (standard SST versions), and SA no longer passes nu_tilde as k. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6
![github-advanced-security[bot]](https://avatars.githubusercontent.com/in/57789?s=60&v=4){kind=small}
…cobians With k in the total energy, p = (gamma-1)(rho E - |rho u|^2/(2 rho) - rho k), so dp/drho contains -(gamma-1) k (general gas: -Kappa k) and the total enthalpy used in the energy row is gamma E - (gamma-1)(|u|^2/2 + k). GetInviscidProjJac (ideal and general gas) and the vectorized inviscidProjJac take k, which is passed by Roe (vectorized and scalar), Turkel, general Roe, JST/Lax (vectorized), AUSM and HLLC; the HLLC star Jacobians include it in dp/drho. The MUSCL sound-speed check also removes the Roe-averaged k from the enthalpy. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6
| */ | ||
| inline const CVariable* findTurbVars() { return nullptr; } | ||
| template<class T, class... Ts> | ||
| const CVariable* findTurbVars(T& first, Ts&... rest) { |
| */ | ||
| inline const CVariable* findTurbVars() { return nullptr; } | ||
| template<class T, class... Ts> | ||
| const CVariable* findTurbVars(T& first, Ts&... rest) { |
…k each MUSCL face's speed of sound - With k in the total energy, the energy components of the right eigenvectors contain |u|^2/2 + k (contact wave and, through the total enthalpy, acoustic waves), and the left eigenvectors use dp/drho = (gamma-1)(|u|^2/2 - k) (general gas: Chi + Kappa (|u|^2/2 - k)). Without it, a contact at rest with constant pressure and k was partly dissipated as acoustic waves. GetPMatrix / GetPMatrix_inv (ideal and general gas) and the vectorized pMatrix / pMatrixInv take k, passed by Roe (vectorized and scalar, L2Roe, LMRoe), MSW, JST-mat, AUSM, FVS, general Roe and the Riemann boundary conditions. - MUSCL: the speed of sound of each reconstructed side is checked too (HLLC takes the square root of each), not only the one of the Roe average. - Unit test roe_tke_tests: P P^-1 = I with k (2D and 3D), and no mass or energy flux across a stationary contact with constant pressure and k (the previous matrices gave 0.024 and 0.083 in the test). Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6
Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6
…in the Riemann boundary states and Jacobians - Roe (scalar, general gas, SIMD) and JST-mat: the part rho_roe*Delta(k) of the jump of rho*E is not a pressure jump. Remove it before the projection with P^-1 and advect it with the contact eigenvalue (also in L2Roe, LMRoe and the explicit general-gas Roe, which build the wave strengths from the primitive jumps). - BC_TurboRiemann: subtract k from the boundary energy before the fluid model; the imposed total enthalpy is thermodynamic, as in BC_Riemann and BC_Inlet (k is added to the energy only). - BC_Riemann and BC_TurboRiemann: pass k to the flux Jacobian of the boundary state. - BC_Inlet (total conditions), BC_Engine_Exhaust, BC_ActDisk: the total speed of sound used the domain enthalpy with k against an imposed total enthalpy without k. - Unit test: contact with a jump of k (at rest and moving) gives the upwind flux for Roe, L2Roe, LMRoe, general Roe. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6
…tream, fix omega at the SST turbo inlet - SLAU, SLAU2, AUSM: the speed of sound used E - |u|^2/2, which contains k with SST. Unit test ausm_tke_tests.cpp: with the same density, velocity and pressure the mass and momentum fluxes must not depend on k (27/27 checks fail without the fix). - SST: FREESTREAM_TURBULENCEINTENSITY = 0, FREESTREAM_TURB2LAMVISCRATIO = 0 or a zero free-stream velocity gave k = 0 or omega = 0 and mu_t = 0/0 (NaN at the first iteration). Stop with an error that names the options. - SST BC_Inlet_Turbo: omega was computed from the averaged solution k instead of the imposed k, so the imposed eddy-viscosity ratio was not the configured one. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6
…M/SLAU Jacobians, declare k as a preaccumulation input - AUSM+up, AUSM+up2: the critical speed of sound used the total enthalpy with k (and so did its analytic derivative). - Accurate AUSM/SLAU Jacobians (USE_ACCURATE_FLUX_JACOBIANS): dp/drho and dH/drho with k held fixed. - Roe, general Roe, MSW, AUSM/SLAU: turb_ke is used in the flux, so it must be a preaccumulation input for the discrete adjoint (the SIMD schemes register it through gatherVariables). Found by the Codex review of 048b10a. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6
… in the SST wall function - SST has no BC_Engine_Inflow/Exhaust or BC_ActDisk: those faces got no turbulence flux at all (silent no-op). Stop with an error until a validated SST version exists. - SetTurbVars_WF wrote the current solution at the wall node, which is overwritten by the update. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6
…markers after the Euler zones are cleared - MSW: the speed of sound used H - |u|^2/2, which contains k with SST; the P matrices of the weighted states now use the weighted k. MSW added to the k-independence unit test (uniform k, 9 checks fail without the fix). - The SST marker check ran before Euler zones drop the turbulence model, so a multizone Euler zone with MARKER_ACTDISK was rejected. Found by the Codex review of b3ae2c7. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6
…of BC_Giles The interior static energy was E - |u|^2/2 (pressure too high by (gamma-1) rho k), and the boundary enthalpy had no k although the scheme is told that both states contain it. Found in a sweep of the sound speed and pressure recomputations for SST k (all other places already hold k fixed or out). Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com> 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