su2code
[WIP] Inconsistencies and improvements to SST model
Proposed Changes
Hi everyone,
I have found some inconsistencies with respect to the literature on the implementation of the Menter's SST model. I would like to use this branch as test bench for any corrections/improvements made to the SST model.
Implementation errors found:
- Use of production limiter constant in the clipping of the cross-diffusion term for the computation of the F1 blending function in CTurbSSTVariable.cpp.
- SST naming conventions (V in V2003m should stay for Vorticity production term.). Not yet assessed.
- Wrong cross-diffusion term included into the residual of w in turb_sources.hpp. It should not be only the positive value as considered in the computation of the F1 blending function. As stated in https://doi.org/10.2514/3.12149. "CDkw is the positive portion of the cross-diffusion term in Eq (A2)" pag. 1604. Moreover, a clipping has been introduced for large negative values of this term, as suggested in Peng, Shia-Hui, Peter Eliasson, and Lars Davidson. "Examination of the shear stress transport assumption with a low-Reynolds number k-omega model for aerodynamic flows." Eq 17.
- Boundary conditions errors at inlets, following the Issue #1851. A fix has been proposed following @pcarruscag suggestions for the supersonic inlet BC.
- Boundary conditions definitions are different than the ones proposed in TMR. Maybe give the user the choice to use the BCs implemented in SU2 or the ones from TMR.
Changes to SST model proposed:
- Inclusion of non modified versions of SST. In this case I think that more work is needed to be sure that the correct terms are taken into account everywhere in the code.
- Give the user the possibility of changing the production limiter for TKE (not essential).
- Change values of default turbulent to laminar viscosity to 1e-2 (NASA TMR reports that it should be in the range of 1e-2 to 1e-5, but in SU2 this is set to 10 as default).
I've seen the proposed changes to the lower limits of k and w in #2323 and I tried implementing it in my branch. However, if the implementation proposed in the respective PR is preferred then I will change mine.
Related Work
#2323 #1851
PR Checklist
Put an X by all that apply. You can fill this out after submitting the PR. If you have any questions, don't hesitate to ask! We want to help. These are a guide for you to know what the reviewers will be looking for in your contribution.
- [X] 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 --allto 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.
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. 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").
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 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.
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 ...
In my experience, by-the-book implementation of SST is very unstable in real-world meshes and applications. For example there is a production limit, but it's proportional to k and w, and thus if those start growing for some reason there is nothing preventing the model from exploding. CD is also tricky due to the 1/w part...
I totally agree. On real meshes, Omega usually drops to extremely low values. In cases where the cross-diffusion term is negative (allowed to be negative), it behaves as a sink term, further amplifying the drop of Omega. A simple addition of this term to the implicit diagonal is insufficient (I tried this). Other more rigorous methods are required (some available in the open literature).
My main argument is that the factor (1-F1) guarantees only the k-epsilon branch outside the boundary layer, in which the cross-diffusion term is already positive. It may not be so because of numerical errors.
Regarding "Fixed Supersonic inlet BC inclusion of TKE," does the energy equation include the contribution of turbulent diffusion?
Hi All,
We are lucky, Florian Menter just replied to me.
He agreed that the factor (1-F1) should activate the CD term only at the "k-epsilon" branch, where the CD term is already positive. Therefore, clipping with zero is in terms of the model and is not incorrect.
Me: A simple boundary layer simulation will reveal that the CD term switch sign is inside the boundary layer and that the F1 function switches from 1 to zero outside the boundary layer. Namely, the CD term should be positive. Florian is unaware of the occasions when it may be negative (even though the model was designed so that this term will be activated when it is positive).
Best wishes, Yair
Regarding "Fixed Supersonic inlet BC inclusion of TKE," does the energy equation include the contribution of turbulent diffusion?
Yeah, I do think so.
I have tested the tripping of the cross-diffusion term also in the residuals on a more complicated test case and indeed it improved convergence. Should I remove this change and use the value of cross-diffusion coming from the computation of F! blending function as was before?
Hi,
May I understand "tripping" as "clipping"?
The cross-diffusion term (that appears in the residual of Omega) is better treated as follows:
2 * (1-F1) * sigma_w2 * rho / Omega * max(grad K * grad Omega, epsilon)
Where "epsilon" could be = 1.e-20 (personally, I'm using epsilon=0)
The term CDkw (that is used to calculate F1) should be implemented in a similar manner:
CDkw = 2 * sigma_w2 * rho / Omega * max(grad K * grad Omega, epsilon)
where epsilon = 1e-20 for the SST1994 and epsilon=1e-10 for the SST2003.
You mentioned that the SU2 implementation of CDkw is not as above?
Yeah, sorry, bad wording from my side. The standard version of SU2 used the value of CDkw as value of the cross-diffusion term in the residual computations. Thus, I think I should revert to that formulation.
Regarding the CDkw, it's implementation in SU2 was made such that the exponent of the epsilon term was equal to minus the production limiter constant of the model. Thus 20 for SST1994 and 10 for SST2003. Is it just a coincidence in the model formulation? In the original paper there is no mention of it being the same variable. This can lead to different results if someone changes the production limiter constant.
Thank you. In the original paper of Menter (Ten Years of Industrial Experience with the SST Turbulence Model) about the SST2003), it is written that the power of epsilon should be minus 10. I think it is a coincidence. You are right; changing the production constant limiter may be problematic. It is better to distinguish between these two constants.
I expect that epsilon to be a simple measure to avoid division by 0, if that lower bound had physical meaning it would have to be multiplied by some reference factors to make its dimensions appropriate, otherwise SST would not give the same results for the same Reynolds obtained with different rho and mu.
You are correct! I have asked Menter about this many, many years back. I don't remember his response, but the formulation given in the TMR is correct (with this fault!)
I have two considerations on the new boundary conditions:
1- should the "domain length" in the Omega farfield BC be adimensionalized by the Reynolds length? Otherwise w gets adimensionalized only by the reference velocity. 2- the farfield value of TKE does not take into account for the freestream turbulence intensity. If we want to match the tke computed via $k_{\infty}=w_{\infty} * (\mu_{\infty} * ViscRatio)/\rho_{\infty}$ and the one computed previously as $k_{\infty}=(3/2)*(U_{\infty}*TI_{\infty})^2$, we can either compute the viscosity ratio combining the two equations or just use the one which explicitly contains the $TI_{\infty}$ disregarding the dependency on $w_{\infty}$. But can we do it? The model as written on the NASA TMR does not explicitly contains the turbulence intensity in the boundary conditions. Maybe it only matters when transition is taken into account.
Hi,
I'm not sure what "new boundary conditions" means. Is someone replacing the present far-field BC of the TKE and Omega? I recognize these "new" BCs as the ones given on the TMR (proposed in the original paper).
I find the present setting in SU2 very comfortable and more "intuitive" than that given on the TMR.
So, the two boundary conditions currently implemented in this branch are:
Standard SU2 BCs $k_{\infty} = (3.0/2.0) * (U_{\infty} * TI_{\infty})^2$ $w_{\infty} = \rho_{\infty} * k_{\infty} / (\mu_{lam,\infty} * ViscRatio)$ $ViscRatio_{default} = 10.0$ $TI_{\infty, default} = 5.0$
NASA TMR BCs for SST1994 and SST 2003 model $U_{\infty}/L \leq w_{\infty} \leq 10 U_{\infty}/L$ $10^{-5}U_{\infty}^2/Re_{L} \leq k_{\infty} \leq 10^{-1}U_{\infty}^2/Re_{L}$
where "L is the approximate length of the computational domain," and the combination of the two farfield values should yield a freestream turbulent viscosity between $10^{-5}$ and $10^{-2}$ times freestream laminar viscosity, as cited on the NASA TMR page.
They give very different values of k and w at the farfield. For example, for a the standard test case of the SWBLI with $ViscRatio = 0.01$ and $TI_{\infty} = 5.0$:
Standard SU2 BCs $k_{\infty} = 2574.55$ $m^2/s^2$ and $w_{\infty} = 3e10$ $s^{-1}$
NASA TMR BCs $k_{\infty} = 0.0013$ $m^2/s^2$ and $w_{\infty} = 1563$ $s^{-1}$
I assume that in the end it doesn't matter since the farfield values drop significantly in the standard version of SST.
I am trying to run again the SWBLI V&V case on the SU2 page. However, as it is provided, it does not work with the SU2 standard boundary conditions. In particular, it always diverges due to NaN residuals in the X-Momentum equation. However, when I use the BCs from the NASA TMR page it works just fine (although only the v2003m model).
Here are the results of the computations performed.
Three configs are used: 1 - With dimensionless limits $k_{lim} = 10^{-20}$ and $w_{lim} = 10^{-10}$. 2 - Without dimensionless limiters. 3 - With dimensionless limits $k_{lim} = 10^{-30}$ and $w_{lim} = 10^{-20}$.
I am reporting the Skin Friction Coefficient on the lower wall against the Experimental results and the trend of the residual for $\rho$ during the simulations. Five different meshes have been considered.
The solutions are pretty much identical for the same mesh level, but the cases with the config 3 achieved convergence across all of the meshes. It is interesting how, on finer meshes, the use of dimensionless limits (config 1) actually had worst convergence than the one with std limits (config 2) (up until the finest mesh, where config 2 diverges).
I am going to try to perform the simulations also with the non m variants of the SST model and the full production term. It may require some work since all of the occurrences of the "ComputeStressTensor" have to be changed. Ideally I think that it should be computed only once and then used that value for all of the other instances. I will focus only on the V2003 model for now as it is the only one that does not diverge.
@YairMO are these results something that you would expect to obtain when changing the lower limits of k and w? With the new updates to the Reynolds Stress computations I obtain mixed results:
On the finest mesh the dimensionless limits diverge regardless of their value (I also tried decreasing them but it didn't work).
Hi,
To begin with, thank you very much for your efforts and the tremendous work you put into improving the SST model!
It is very difficult for me to judge because I'm using a different time integration than the SU2, which may be critical in this regard. Moreover, my code is non-dimensional, so estimating the influence of the "newBC" and the "lower limits" is difficult for me. Can you tell the values of the "lower limits" relative to the freestream values? Please provide many details about the inputs in your config file (Linear solver, linear solver convergence, CFL, inviscid flux, etc')
I found the SWBLI config file. I noticed that it is the NK solver. I would test this case without the NK solver. Moreover, I've noticed that the entropy fix is very low; how about raising it to 0.2 or using a different flux (JST, AUSMPLUSUP2)? It may sound irrelevant, but it is worth the try.
Changing the entropy fix coefficient (EFC) to 0.2 had the largest effects across the board. It actually degraded the convergence rate, but it removed divergence of the solvers for finer meshes.
The converged solutions do not change when changing the EFC, which is desirable.
Great! This is a non-linear example: Can we say that the convergence pattern is not necessarily directly affected by the turbulence model's lower limits? Did you try avoiding the NK solver?
I don't think we can actually say that, or, at least, we can say that it is not the only factor affecting the convergence. The problem is that it is not consistent across the meshes: for sure, using dimensionless limits help, but changing their values has different effects for finer meshes.
The effect on the residual of the turbulent variables is the same. The most problematic one is always the Omega equation where the residuals remain always large.
The problem is that Omega is too large with respect to all the other variables, especially near the wall (where the residuals remain high).
I agree. The NK solver is not well suited to solving the turbulence model. Maybe the linear solver alone will be more robust.?