Turbulent Kinetic Energy(TKE) on energy equation in SST model.
sun5k opened this issue · 21 comments
Describe the bug
I found a bug while testing the SST model.
I'm trying to simulate the high Mach number flat plate case. The flow condition is below.
| Temperature (K) | Unit Re (1/m) | Mach | Pressure | Wall Temperature(K) |
|---|---|---|---|---|
| 800 | 4.9×10^6 | 6.1 | 12100 | 300 |
After calculation, I check the flow field, including the boundary. But, the inlet boundary is not the same flow conditions I imposed.
Here is the comparison of other models using the same flow conditions.
So, I tried to find the cause.
When deriving the RANS equation from the energy equation, the turbulent kinetic energy(TKE) is derived from the kinetic energy of Total Energy(E) as shown below.
In the boundary condition and initialization part, k is considered depending on whether or not the SST model is used for Total Energy.
In CEulerSolver.cpp
bool tkeNeeded = (turbulent && config->GetKind_Turb_Model() == TURB_MODEL::SST);
if (tkeNeeded) { Energy_FreeStream += Tke_FreeStream; }; config->SetEnergy_FreeStream(Energy_FreeStream);
if (tkeNeeded) Energy_FreeStreamND += Tke_FreeStreamND;
But In the flux calculation part, when calculating the speed of sound(a), using the total enthalpy(H). But, in this part consider only the kinetic energy.
Expressed mathematically
but in the SU2 code using below.
The speed of sound always has a slightly high value because TKE is not subtracted.
See "roe.cpp" line 119 ~ 149.
RoeSoundSpeed2 = (Gamma-1)*(RoeEnthalpy-0.5*sq_vel);
This line might be changed roughly like :
RoeSoundSpeed2 = (Gamma-1)*(RoeEnthalpy-0.5*sq_vel-k);
There seem to be two possible ways to this problem.
- Subtract TKE when the speed of sound calculation
- Store e(static energy) or k in primitive variables
I'm not sure which is clear and easy.
Bug report checklist
Please make sure that you have followed the checklist below, many common problems can be solved by:
- Consulting the build instructions (https://su2code.github.io/docs_v7/Build-SU2-Linux-MacOS/).
- Looking for similar problems on GitHub or CFD-Online (https://www.cfd-online.com/Forums/su2/).
- Updating to the newest version of SU2 (either develop, master, or the pre-built executables https://su2code.github.io/download.html).
Desktop (please complete the following information):
- OS: [e.g., Linux (Redhat 6.6)]
- C++ compiler and version: [e.g., g++ (GCC) 4.8.5]
- MPI implementation and version: [e.g., OpenMPI 3.0.0]
- SU2 Version: [e.g., v7.4.0]
This is a good point, as you know we've tried to move to the "m" versions of SST (#1483) (dropping k from viscous terms) but I don't think we have looked too deeply into the convective terms.
However, the Roe speed of sound is only used by dissipation terms, the consistent part of the flux is based on primitive variables, and those account for tke, so I suspect there is some other inconsistency at play when defining farfield conditions or the farfield boundary itself.
Also note that you are probably using very high values of free stream turbulence intensity and viscosity ratio, the defaults are (5% and 10x). I did not see this type of problem in https://su2code.github.io/vandv/swbli/
This is a good point, as you know we've tried to move to the "m" versions of SST (#1483) (dropping k from viscous terms) but I don't think we have looked too deeply into the convective terms. However, the Roe speed of sound is only used by dissipation terms, the consistent part of the flux is based on primitive variables, and those account for tke, so I suspect there is some other inconsistency at play when defining farfield conditions or the farfield boundary itself. Also note that you are probably using very high values of free stream turbulence intensity and viscosity ratio, the defaults are (5% and 10x). I did not see this type of problem in https://su2code.github.io/vandv/swbli/
Thank you for reply.
I just ran the V&V SWBLI case. Both k-w SST and SA model results were no problem.
In the k-w SST case, When I arbitrarily changed free stream turbulent intensity 5e-3 to 5e-2, the flow condition at the inlet was changed. You can see the below figure(The S-A and kw(5e-3) is on the same line).
The problem seems to arise when you have a large TKE.
@sun5k @pcarruscag I noticed this issue when rearranging NEMO to handle SST cases. I believe several users on the CFD-online have noticed these issues for high-mach flows.
Ok, the issue is that the inlet values in SST tend to decay very quickly (as you can find in the literature https://turbmodels.larc.nasa.gov/sst.html) so far as I can tell we are adding the right amount of TKE to the flow and turbulence solver inlets.
However, in the turbulence solver this is destroyed immediately (note that destruction is proportional to TKE) and so TKE at the inlet nodes is not TKE_Inf.
Effectively this means we put too much energy into the flow, which is not matched by the turbulence solver due to the immediate decay, and therefore static temperature appears high, and Mach number low.
E.g. we put 10 + 0.1 into the flow solver inlet, and 0.1 into the turbulence solver, but this 0.1 decays to 0.01 right away. Then, when we compute the temperature for this inlet node we have 10 + 0.1 total energy (which does not decay in the flow solver since it does not have destruction terms) from which we subtract the decayed 0.01 from the turbulence solver and end up with 10.09 instead of 10.
I tested this by turning on the SST sustaining terms, which prevent decay below freestream values, and I got the right Mach number.
So what we could do is use the TKE from the turbulence solver in the flow solver inlets, instead of the freestream TKE which the turbulence solver is not guaranteed to sustain.
Since this is not really a bug but a result of turbulence boundary conditions that are far from their equilibrium values, and since we can now apply boundary conditions for turbulence quantities directly using MARKER_INLET_TURBULENT, I think we can close this issue.
I think we can do better.
HI ~
Recently, I read a some paper [DOI:10.1017/aer.2020.93]
The main topic of the above paper is considering the TKE term in the RANS equation for supersonic and hypersonic.
In the paper, the original RANS equation is:
The paper said CFL3D and VULCAN codes ignored every TKE term, and
“Since different authors treated the k terms in different ways, sometimes it is difficult to obtain a uniform evaluation from different papers for the same turbulence models on the performance of supersonic and hypersonic turbulent flow predictions.”
So, the paper classified 3 types of RANS equations (C1, C2, and C3) like below:
The paper’s result about k-w SST
Mach 5 flat plate :
Mach 2.25 impinging shock :
Mach 8.18 impinging shock :
Mach 2.9 compression corner :
Mach 7.05 compression corner :
It looks like the attached boundary layer region(like a flat plate) is almost the same, while the separation region(impinging shock and compression corner) shows slightly different results. The results indicate that C1 is quite reasonable, but according to the RANS approach, C3 is right. So I'm not sure which one is correct.
Return to our issue, I quickly try to C1 (ignore all “if(tkeneed) ~~”) in SU2.
Here are results
Mach 6 flat plate :
Mach 5 impinging shock (v&v case) :
C1 and SA results are the same, and "Current Tu 5e-4" is almost the same but slightly different.
We could make C1, C2, and C3 options in SU2 later. But now, I recommend temporarily applying the C1 combination because our boundary condition can’t persist imposed flow conditions when turbulence intensity is high.
Do you have some pictures of what the pressure looks like inside the boundary layer with the approaches you tried?
You are proposing the exact opposite of the conclusion of the paper:
"From the above findings, it is recommended that all three of these terms be included when
running hypersonic, or even supersonic, turbulent flow simulations, especially for flows with
shock wave-induced separations."
And they clearly say this:
"While the full inclusion of these terms does not always result in predictions that agree better
with DNS/experimental data, this is likely caused by the fact that their exclusion cancels out
effects of other flaws in the RANS models employed."
If your strategy is to get a better match with experiments by neglecting physics terms, then you should rethink your strategy.
Do you have some pictures of what the pressure looks like inside the boundary layer with the approaches you tried?
Mach 6 flat plate case, the flow conditions are (same as the above, but rewritten below):
| Temperature (K) | Unit Re (1/m) | Mach | Pressure | Wall Temperature(K) |
|---|---|---|---|---|
| 800 | 4.9×10^6 | 6.1 | 12100 | 300 |
FLUENT results are just for comparison of reference profile shape.
The other result is almost the same except for Inlet and pressure profiles.
If you want SWBLI problem profiles, I can show you.
You are proposing the exact opposite of the conclusion of the paper:
"From the above findings, it is recommended that all three of these terms be included when running hypersonic, or even supersonic, turbulent flow simulations, especially for flows with shock wave-induced separations."
And they clearly say this: "While the full inclusion of these terms does not always result in predictions that agree better with DNS/experimental data, this is likely caused by the fact that their exclusion cancels out effects of other flaws in the RANS models employed."
If your strategy is to get a better match with experiments by neglecting physics terms, then you should rethink your strategy.
Thank you for your comment @bigfooted
The above paper is not presented to improve the current k-w SST model.
As you can see in the first post, there is a problem with the high Mach number and freestream turbulence intensity case. If high turbulence kinetic energy(TKE) and Mach number condition, the boundary condition cannot be maintained the imposed value. (I think there seems to be a bug in the temperature calculation using total energy when including the TKE)
The introduction of C1 was intended to provide a 'temporary' solution at the level of first aid (simply commenting out conditional statements in code).
@sun5k maybe you can try implementing the alternative I suggested of obtaining turbulence kinetic energy at inlets from the turbulence solver instead of assuming that the free stream value is imposed exactly?
In the SWBLI case the SST-m versions (without divergence terms in the stress tensor) seem to underpredict the separation region https://su2code.github.io/vandv/swbli/
Thank you for replying while busy preparing the High Lift Prediction Workshop 5.
Sorry, I didn't clearly understand. From what I understand, I can suggest another energy equation calculating method instead of the current SU2 method(reading calculated TKE from inlet boundary condition and using it as an energy equation). Is that right? If not, could you please explain in more detail?
@sun5k maybe you can try implementing the alternative I suggested of obtaining turbulence kinetic energy at inlets from the turbulence solver instead of assuming that the free stream value is imposed exactly? In the SWBLI case the SST-m versions (without divergence terms in the stress tensor) seem to underpredict the separation region https://su2code.github.io/vandv/swbli/
In places like if (tkeNeeded) Energy += GetTke_Inf(); we would instead read tke from the turbulence solver.
Hi everyone,
I am interested in this problem. Are there any news from your side @sun5k ?
I may try to help in running some test cases or changing the code if needed.
Hi @rois1995
For now, I'm ignoring all TKE in Total Energy in personal research.
I don't remember the details clearly.
The problem was that the enthalpy added TKE was stored in the primitive variables.
When I tried to fix it, the problem was when the enthalpy added TKE was stored in the primitive variables. For the Roe scheme in convective flux calculations (not sure about other flux scheme), the Roe speed of sound is calculated using enthalpy. But as I mentioned above, the stored enthalpy is higher than other simulation because of TKE. I thought it make a problem.
Thank you for your answer. If you could share the configs and meshes that you are using I can try following the approach suggested by @pcarruscag and use the TKE from the turbulence solver instead of the freestream one.
Plus, I have some doubts on the default value of the turbulent to laminar viscosity ratio which is equal to 10 in SU2, although on the NASA page suggests to be in the range of 10^-2 to 10^-5. However, I think that this is another issue.
Here is SWBLI config and mesh file I used. The original freestream turbulence intensity is a 5.0e-4, I was changed arbitrary to 5.0e-2 for test. Also you can easily observe any high Mach number simulation with high TKE.
I tried changing the supersonic inlet function and take the value of TKE coming from the turbulent nodes instead of using the freestream one and the Mach number at the inlet is much closer to the specified one.
Here I am plotting for the SWBLI case the Mach number at the inlet for SA, SST with Turbulence Intensity (TI) 5e-3 with and without Fix, and SST with TI = 5e-2 with fix (without fix the Mach number was way off). There is still some discrepancy wrt the SA results, but it is an improvement.
Thank you for your answer. If you could share the configs and meshes that you are using I can try following the approach suggested by @pcarruscag and use the TKE from the turbulence solver instead of the freestream one.
Plus, I have some doubts on the default value of the turbulent to laminar viscosity ratio which is equal to 10 in SU2, although on the NASA page suggests to be in the range of 10^-2 to 10^-5. However, I think that this is another issue.
@rois1995 Hi ~
Can you give me a link to the NASA TMR guide for setting the viscosity ratio? I'm not sure where to find it.
Sorry, I found it!
@sun5k it is on the SST page (https://turbmodels.larc.nasa.gov/sst.html) at the end of the explanation of the standard model.