SA Model testcase

[zhan1457]

Hi,

I try the SA model on a flat plate but found the nu_tilde is pretty small. Do you have any testcase for RANS SA model to share? I noticed that there is a NACA0012 case on your paper, but could not find it in the source code.

Thank you.

[jewatkins]

We are still working on developing a tutorial for running a case using the SA model. This feature is still under development but you can reproduce the results from the paper by constructing an appropriate mesh and using the code provided.

[zhan1457]

Hi Jerry,

Thanks for reply. Actually, I found some inconsistent of the code and your paper with the reference (Moro, 2011) you used. Besides, the code will face a NaN (0/0) problem when psi and S_tilde are both zero. Have you noticed them when you developing the tutorial for SA model?

[jewatkins]

The reference to (Moro, 2011) is only used to cite the modification used in the code to prevent large negative values of the modified eddy viscosity from causing the solution to become unbounded or creating anti-dissipation. To be more specific, it prevents psi from becoming negative. There are other modifications taken from (Oliver, 2008) so the code does not mimic the equations in (Moro, 2011) exactly.

If the modified eddy viscosity becomes negative infinity, it is possible for psi to become zero and the ratio of psi over S_tilde to become NaN if there is no vorticity. In this case, the reference length, r, is limited to a max value of 10.0 (see turbmodels.larc.nasa.gov/spalart.html, under the "Standard SA model") since there is no vorticity. This is the only location where the ratio psi over S_tilde is used so this should prevent the NaN from polluting the solution.

This feature is still under development and so encountering a NaN is highly possible. If you do encounter a NaN, please let us know.

[zhan1457]

Hi,

I checked the Oliver's paper you mentioned and find that he used a different way to prevent the negative turbulent viscosity by forcing the nu_t to be no less than zero. Maybe that's why you did not use psi to compute f_v1 as Moro did. However, the f_v2 is also dependent on f_v1, so I am not sure if it is OK to use both of the modifications together for neglect viscosity. The modification of r you mentioned does not prevent the NaN spreads in the code, since you use a min(NaN,10) in the code and it will return NaN in my testcase.

In addition, the compressible turbulent viscosity equation in Oliver's paper is different from the equation from Allamaras' (Modifications and Clarifications for the Implementation of the Spalart-Allamaras turbulence model). It misses one term generated from the diffusion term. I have not checked how important this term is, but these two equations should be different. Allamaras et al. consider the S-A equation is applicable for both incompressible and compressible flows. The conservation equation is constructed by combining S-A equation with mass conservation eqn, but not simply add the density into the S-A equation.

Besides, the formula of S_bar(line 82 in source.cpp) uses a square of nu_tilde which must be a typo. And not all the boundary conditions are modified if RANS turbulence models are used.

[jewatkins]

Hi,

The exact equations I used were taken from "Robust Computation of Turbulent flows using a Discontinuous Galerkin Method" by Burgess/Mavriplis. When I wrote source.cpp, my main goal was to try to reproduce some of their results from this paper and another paper called "High-order Discontinuous Galerkin Methods for Turbulent High-lift Flows". I can see the discrepancies you mentioned including the formulas for f_v1 and S_bar. It might be interesting to run a few cases to see how this affects the final solution.

The min(NaN,10) does indeed return NaN. In that case, I don't think I have run into the case where the ratio of psi over S_tilde is NaN. Nevertheless, this should be taken to account. I'll have to look into this more thoroughly. If you're experiencing this NaN, try adding a check in source.cpp to ensure r doesn't become NaN. Something like: if(isnan(r)) {r = 10.0}; Let me know if that fixes the problem.

I have yet to add the modification you mentioned from Allmaras' paper. I wanted to first verify that the "Standard" SA model was implemented correctly.

I have yet to modify all the boundary conditions to work with RANS turbulence models. I was experiencing a few problems with boundary conditions so I may create new boundary conditions specific to RANS.

Thank you for all your input. Please let me now if you find anything else or if you have anymore questions.

-Jerry

[zhan1457]

Hi,

I have modified the S_tilde to be no little than 1e-10 when psi is zero as SU2 does and it works. I think the equation in Allmaras' paper should be the standard form, because only some math work is involved.

Recently, my time step is pretty small to avoid divergence. It seems not natural. Could you share the input configure file for SA model? What time step and how many elements did you use?

Thank you.

Wanjia

[zhan1457]

Hi,

I have tried only Moro's modification and it works now. In addition, the standard SA model you mentioned in the NASA website has applied trip terms. Does it matter not having it in your code? However, still waiting for the result. Is the local time stepping ready to use? It is kind of not user-friendly to stop and restart the simulation each time, since the history.plt will be overwrite and we need calculate the total time by ourselves. Thank you.

Wanjia

[jewatkins]

Hi Wanjia,

I'm glad you got something working. The time step is small because it's extremely limited by the mesh size for the explicit time stepping schemes used in the code. I can send you a configuration file if you still need it. The time step I used for the 6,539 element NACA0012 airfoil mesh at Re = 6mil was as small as 1.0e-7. You can use much larger time steps in other codes because they use implicit time stepping to converge to steady state.

The trip terms are not in the code. The code is only setup to run fully turbulent flow. Local time stepping is working but it is still incomplete and hasn't been fully tested. We are still working on methods to improve convergence for steady state cases so you can look for that in the future.

If you set "n_steps" in the configuration file to a large number, you can just leave the code running in the background. I've also ran into the frustration of history.plt being overwritten. I'll look into changing that. By total time, do you mean physical time or wall clock time?

-Jerry

[zhan1457]

Hi Jerry,

Thank you for providing the information. Mine is the same order of magnitude. I also found the polynomial order will influence the converge time step a lot as well.

I appreciate your effort to improve the code. I am looking forward to more good news from yours.

I mean physical time, the output file only shows the number of iterations.

Wanjia

[jewatkins]

Wanjia,

I used a polynomial order of 3. From my experience, a lower polynomial order will allow you to use a greater time step. A higher polynomial order will require a smaller time step.

Thanks for taking time to look at the code. I'll add physical time to the list of things that need to be fixed in the history file output. Let me know if you have any other questions.

-Jerry