Code Programming: Problems with skewed cells at transient calculations

Corby · March 27, 2014, 12:43

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

agd · March 27, 2014, 22:05

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.

Corby · March 28, 2014, 05:01

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.

agd · March 29, 2014, 20:47

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

Corby · March 31, 2014, 05:18

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.

March 27, 2014, 12:43	Code Programming: Problems with skewed cells at transient calculations	#1
Corby New Member Join Date: Feb 2013 Posts: 5 Rep Power: 13	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 rhoVol/deltaTime and the source term is modified by rhoVol/deltaTimeV_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

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
snappyhexmesh remove blockmesh geometry	philipp1	OpenFOAM Running, Solving & CFD	2	December 12, 2014 11:58
Highly Skewed Cells	chrisoturner	FLUENT	7	July 22, 2010 07:43
ghost cells in axisymmetric problems	ma	Main CFD Forum	2	May 17, 2004 14:52
Difference between steady and transient problems	Jianping Zhang	Main CFD Forum	9	April 14, 2003 07:17
memory problems with a code on unix	vince	Main CFD Forum	2	November 5, 2002 04:48

March 27, 2014, 22:05		#2
agd Senior Member Join Date: Jul 2009 Posts: 358 Rep Power: 19	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.

March 28, 2014, 05:01		#3
Corby New Member Join Date: Feb 2013 Posts: 5 Rep Power: 13	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.

March 29, 2014, 20:47		#4
agd Senior Member Join Date: Jul 2009 Posts: 358 Rep Power: 19	If you are using a staggered grid, then what velocities are you interpolating?

March 31, 2014, 05:18		#5
Corby New Member Join Date: Feb 2013 Posts: 5 Rep Power: 13	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.