[WIP] Inconsistencies and improvements to SST model

[rois1995]

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 Turbulent Kinetic Energy(TKE) on energy equation in SST model. #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.

• 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.

[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]

Sure, thank you for checking!

Regarding the SWBLI case, if I understood correctly, you have used your code with NK solver from the start and you changed the number of linear solver iterations, right? Can you also try with SST or is your implementation different only for SA?

[YairMO]

Yes, your understanding is correct.
As for the SST ... checking

[YairMO]

Here are the results obtained using the SST2003 model

[rois1995]

It seems that yoru implementation affects also the SST model, is it correct? Nevertheless, huge improvement wrt both this branch and develop branch results (density residual stagnates around -10)!

[YairMO]

Yes, it also affects the SST model. Thanks

[YairMO]

Hi Rois1995,

I have permission to share the VFE-2 grid; please send me your private email.

[rois1995]

That's perfect, thank you! [email protected]

[rois1995]

I have a few updates. I found out that most of my FGMRES problems in this branch where related to the fact that I was using single precision for the linear system. You can see it in this graphs for the residuals of Rho.

The first two curves are in mixed-precision and they stop long before reaching the minimum residual required due to divergence of the linear solver.

Moreover, it seems that the 2003 model here implemented (which, simply put, considers the full reynolds stress tensor for the computation of the production of k and always considers k in the stress tensor and in the thermodynamic variables) has faster convergence and increases the recirculating zone (maybe due to the reduction of turbulence kinetic energy) which is a feature that has also been seen in [DOI:10.1017/aer.2020.93].

There are some differences between the use of the TMR boundary conditions and the ones used before.

The results improve with respect to the develop in the prediction of the SFC distribution, which has results completely off from the V&V page of SU2

I may keep on checking if the implementation of the v2003 model is correct or not by searching for other test cases (probably coming from DOI:10.1017/aer.2020.93).

[YairMO]

How can one switch between single or double precision?

[rois1995]

when you execute the meson command use the flag -Denable-mixedprec=true

[bigfooted]

I am trying to catch up on the developments in this PR but I've lost the overview. Are there still issues that we need to figure out? It looks like convergence is now recovered and results look OK?
Regarding the symmetry plane, as @pcarruscag mentioned it is improving the accuracy, but not necessarily the convergence. @pcarruscag and I tried to get it to converge as well as before, but no guarantees. Do you still see issues with the new symmetry plane?
Regarding the peak residual that you observed in a single node in the mesh: This is a common observation when using Green-Gauss, which gives O(1) accuracy in mixed mesh nodes where quadrilateral and triangular meshes meet. There are some methods that try to fix this, although none of them are holy grails and the best methods are the most expensive ones.

[rois1995]

The 2003 version of SST seems to be working fine. However, the 1994 version is not. I get Mach 7 values on the upper symmetry boundary even on the coarser mesh. I am trying to understand what is happening.

[rois1995]

Maybe a stupid question, but I noticed that if I print the coordinates of the i-th point in the CTurbSSTSolver class and the problem is 2D, then the z-coordinate is not 0 but has random values. Does this affects anything in the code or it is not important since the problem is 2D?

[pcarruscag]

The z coordinate is not allocated for 3D, so you are either reading out of bounds or reading the x or y value from somewhere else.

[rois1995]

That's what I thought, thanks for confirming it.

[rois1995]

One more thing: should I change the naming scheme for the model version of SST? As it is now, the naming scheme is V+model version but the NASA TMR site says to use the V only for the vorticity form of the production term. Should we remove the V then, or maybe set it as lower case? This will break all of the configs, no problem for me to go through them all and change it though.

[YairMO]

Chris Rumsey made a tremendous effort to standardize popular RANS turbulence models (via the TMR website).
In my opinion, follow his work. Namely, change the present names to exactly what is given on the TMR.

[pcarruscag]

No changes that break backward compatibility.

[rois1995]

That's what I thought.

[bigfooted]

The original description was correct, for species transport we can have a Prandtl number for each of the species

[rois1995]

Yeah I do not recall changing that but for sure it was wrong in my branch. Corrected, thanks!