Code Programming: Problems with skewed cells at transient calculations
 

Hello everybody,

I implemented an own RANS code in Matlab with the following main specifications:

• 2D, incompressible RANS
• structured, staggered grid on curvilinear coordinates
• pressure-based coupled solver (implicit)
• turbulence model: Menter SST

I already achieved good results for steady state flow problems both for orthogonal and non-orthogonal meshes. As I want to solve separation driven problems I have to introduce the transient terms now to increase the probability to get accurate results. I followed the approach presented by Yang et al.: "General Strong Conservation Formulation of Navier-Stokes-Equations in Nonorthogonal Curvilinear Coordinates", AIAA Journal, Vol. 32, No. 5, May 1994. The dependent variables in the momentum equations are the contravariant velocity components projected on the axis of the resprective curvilinear direction. As it seems to be standard the transient term is approximated by first order. After integration a_p is enlarged by rho*Vol/deltaTime and the source term is modified by rho*Vol/deltaTime*V_curvilinear_n=0.



Everything appears to work fine on orthogonal grids. But when I want to calculate an unsteady solution (CFL between 1 and 1000) on a nonorthogonal grid I get non-physical results preferably in skewed cells (max 10° deviation from 90°). Since one week I am trying to figure out the problem but now I have no more idea where to search the error.


So here is my question: Did anybody else have a similar problem and could give me hint? As all the necessary parameters (control volume, contravariant velocity components from previous/initial time step) seem to be ok, there must be a fundamental fault in my understanding of the equations.

Thank you very much in advance!

Corby

Have you done a uniform flow check (turn off boundary conditions, initialize the flowfield to uniform values, run the solver) to see if the solver maintains a uniform flow? Discretization of the grid metrics can show up as a mass source if done incorrectly.

Thanks for reply!
I checked the uniform flow case according to your proposal. The flow stays uniform inside the numerical accuracy limits, i.e. maximum deviation from 1 is 1e-29.
Other orthogonal test cases showed reasonable results. I think I would have noticed if there is a bug in the metrics implementation that can be seen also on orthogonal cases.

I have the idea, that the bug must be hided in the way I perform the velocity interpolation on the curvilinear grid which is not obvious for me at all. Until now I did not find literature on this topic that is not written by a mathematician, i.e. comprehensible by an engineer and CFD dummy like me.
However, my implemented interpolation method that is also necessary for the steady calculation (fluxes on staggered grid cell borders) seems to work well.

If you are using a staggered grid, then what velocities are you interpolating?

In general you have to interpolate the velocities from the previous iteration to calculate the fluxes at the borders of the staggered grid cells for the momentum equations. For non-orthogonal meshes an "exact" interpolation is quite tricky. Here, I simply calculated the weighted linear interpolation based on the distance to the next resprective node. With curvilinear coordinates the direction of the coordinate axis and thus the contravariant velocity components change. As the direction of the mid point must be somewhere in between the simple linear interpolation of the contravariant velocity values does not seem to be catastrophically wrong. It seems to work well for the steady case.
But for transient calculations you have to introduce the interpolation implicitly into the matrix to be solved. Here, it is necessary to interpolate the second velocity component (in 2D) on the node where you evaluate the momentum of the first velocity component. For example, for the momentum along xi (horizontal) you consider the respective cell where the vertical components along eta are (unfortunately) stored on the corners. Since you actually need them on the cell centre one has to inpolate. This is at least my view of things. So the chosen interpolation method (if there are really several) implicitly affects the solution only through the transient term. Even if I revisited my understanding of the equations and my implemented discretisation several times, I cannot understand, why the transient solution is locally wrong at skewed cells where the flow is unphysically deviated.