Convergence issue of the k-ω SST turbulence model in Nalu-Wind

[BumseokLee]

Hi,

I am reporting a convergence issue with the omega transport equation of the k-ω SST turbulence model in Nalu-Wind. I’ve attached non-linear residual plots for NASA’s canonical case: the subsonic NACA0012 airfoil at Re=6e6. I used the 2-D C meshes and flow conditions from NASA Langley’s Turbulence Modeling Resource (TMR) Website. The freestream Mach number of the case in NASA TMR was 0.15 but I reduced it to 0.1.

In fact, this is not the only case that I’ve observed suspicious or inconsistent behavior in the omega transport equation. I’ve encountered similar issues (very large residuals or inconsistent convergence trend for a different number of Picard iterations) in other cases, including 2-D airfoils and 3-D blade simulations. In my tests, the sustaining terms didn’t resolve the issue.

Despite these behaviors, the residuals of the momentum (mean flow) and k still converge reasonably. Additionally, the predicted aerodynamic coefficients converge, showing good agreement with available experimental data or other numerical results. This might be due to the clipping functions of the model and/or the current implementation in Nalu-Wind.

However, the magnitude of the residuals still looks concerning. The convergence stall of the omega could slow down the overall convergence of the turbulence model and, consequently, the mean flow as well. This issue could be particularly significant in highly unsteady flow simulations, such as IDDES FSI simulations, where sufficient iterations are not available between physical time steps. It could also affect laminar-turbulent transition simulations too because the transition model triggers the transition onset within the boundary layer using the k and omega values, making it highly sensitive to near-wall values. (That’s why I previously fixed the omega boundary condition at the wall, although it doesn’t seem to make a significant difference in fully turbulent simulations.)

I haven’t yet investigated where these large residuals occur, but I believe it is worthwhile to check. I am currently occupied with some other work in my projects, so I cannot attack this issue immediately. However, I can investigate this after completing my urgent work, probably around early or mid-September. If anyone is interested in this, you are welcome. I can share the mesh and input files, or you can directly start from the NASA mesh to reproduce the issue.

[sayerhs]

@BumseokLee I'm curious to learn how you are setting the BCs for inflow/outflow with C-mesh topology when you do non-zero angles of attack.

[BumseokLee]

@sayerhs I simply set the inflow/outflow BCs as below regardless of AoA.

[sayerhs]

@BumseokLee What happens at the horizontal sections on the top and bottom at the back? Depending on angle of attack part of that will be inflow/outflow, correct? Is that being adjusted?

[BumseokLee]

@sayerhs
Do you think the horizontal parts of the top and bottom could result in the convergence stall or oscillations? I am curious.
For some reasons, the convergence issue is only observed for the omega transport equation.

I've also attached the turbulence eddy viscosity at AoA=15deg below. For me, it looks fine.

For the 449 X 129 mesh, AoA=5 deg shows the most significant convergence stall.

[sayerhs]

@BumseokLee It is just an inconsistency to be mindful when executing high angles of attack. This was the reason Ganesh and I had resorted to O-mesh configuration when doing the airfoil runs back in 2018-2019 timeframe.

The omega oscillation is something we have seen with O-mesh as well. From what I remember, once it triggers it never really goes away in Nalu-Wind. We never got to the root cause on that one. My suspicion was the wall BC treatment. But it seems that you've already looked into it.

[sayerhs]

@BumseokLee could you also please share the velocity vectors or streamlines for the same case as the eddy viscosity plot? I am curious to see what happens to the streamlines at those horizontal edges.

[BumseokLee]

@sayerhs

I've attached the velocity vectors near the upper right corner. It looks fine for me. I also checked the velocity contours but didn't see any suspicious behaviors such as reflection waves from the boundaries or sudden changes in the values.

Yes, I'm also aware of that this simple BCs with C-mesh can be problematic for very high angles of attack such as AoA=90deg. That's why Shreyas Bidadi generated O-meshes at each AoA in his hybrid RANS/LES simulations, converging AoAs from 0 to 360 degrees. However, all of the NASA’s and AIAA’s meshes are C-meshes. Additionally, the far-field boundary is located at 500 chord lengths from the airfoil surface, and the AoA is just 15deg, so I just assumed it should be fine.

I have a quick question. Is the far-field Riemann boundary condition (or similar type boundary condition) not applicable in incompressible flow solvers? My background is still mainly in a structured compressible flow solver (Maryland solver, OVERTURNS), so this might be a naive question, but I'm curious if there's a possibility of Riemann type boundary conditions for the far-field, so we don’t need to worry about the correct far-field boundary conditions for airfoil cases.

I already checked the wall boundary condition of the omega transport equation and fixed an issue with the default value, but it didn't resolve the issue. Given the magnitude of the values, I suspect 1) there's a division by very small numbers (close to zero) that's being clipped by a function within the model or the current implementation or 2) there might be some numerical issues leading to negative omega values. You may remember that the original S-A turbulence model had problems with unphysical negative turbulent eddy viscosity. To fix this, Spalart and Allmaras proposed the positivity in the implicit operator in the SA-92 model and then later proposed SA-neg models in 2012. A similar issue was also observed in the Menter’s gamma transition model, which exhibited negative intermittency, so I applied a similar approach to the implicit operator to fix that.

If the current issue is caused by simply inaccurate applications of the boundary conditions, then, it would be great because nothing to fix in the codes. However, if not the case, I think I need to investigate what happens in Nalu-Wind.

[sayerhs]

Thanks for sharing the velocity vector plot. You will also want to check the pressure field along with the velocity field as the pressure will adjust if there is any change in the streamline direction.

Is the far-field Riemann boundary condition (or similar type boundary condition) not applicable in incompressible flow solvers?

OpenFOAM has an inflowOutflow BC which is probably the closest in behavior to what you are looking for. This is something we had talked about implementing in Nalu-Wind, but never got around to doing it. This would switch the BC type (inflow or outflow) depending on the local flow direction.

[BumseokLee]

@sayerhs

I've attached the pressure and velocity contours. I intentionally adjusted the range of the contour levels to see any suspicious behaviors. For example, the freestream velocity is 34.1m/s, and the contour level for the velocity contour is from 34.0 to 34.2. Some effects of the BCs are observed but they don't look suspicious for me.

While I was checking the contours, I realized that the highest omega values occur at the wall. This is due to the boundary condition, which includes the division by the wall distance square.

Below figure shows the omega near the trailing edge. It is seen that very high magnitudes of the omega are observed at the first few cells in the wall normal direction and at the wake average mesh. I may also need to compare this against converged cases and see how they look different.

[sayerhs]

@BumseokLee Thanks for plotting the pressure and velocity and showing them side by side. I think that creating similar plots for all the angle of attack cases would be useful to understand the issue.

Here is what I observe in your results. While benign, there is a wavy variation in pressure field, and a corresponding change in the velocity field as it nears the top horizontal section (you can see that as a faint bluish blob in the top right corner). Also it looks like pressure jumps to zero at outflow across one cell. You'll want to interrogate the points along the top boundary and the interior to be sure. I believe this is because we are enforcing a dirichlet velocity (inlet BC) at the top horizontal section and the flow is adjusting the pressure field to match this requirement as the exit velocity is not the same as the inlet velocity.

$\omega$ being maximum at the airfoil boundary is expected as $\omega \to \infty$ near walls. However, the $\omega$ at the wake cut is spurious and needs further examination. Are these separate sidesets in Nalu-Wind? And, if so, are you explicitly setting BCs there? Or are they coincident nodes?

I would also recommend that you restart from a converged solution and run a few timesteps and output solution at each step. That would help you visualize where in the flowfield $\omega$ is fluctuating, and if it is having any effect on the pressure/velocity field.

It might also be useful to compare the solution against some of the $O-$mesh simulations you mentioned that you have available.