You have to check the skewness and make sure that the maximum value is less than 0.85

it seems you are using low cell number I suggest to use more than 4 millions

You have to check the skewness and make sure that the maximum value is less than 0.85

Make sure about the number of the mesh, you have to make mesh dependence on knowing the best number of nodes.

Thank you M.f Abd and Ghasan Fadhil.

May I ask Mr. M.F. Abdo do you mean to increase the number of nodes in the local mesh? I am now using a mesh of 1 L and 0.5 W with about 10000 elements. U know I am using spectral element method that makes the mesh finer based on number of polynomial order. Do u mean the original mesh should be increased up to 4 million meshes?

Mr.Ghassan Fadhil, I have checked the mesh dependence for 10000 nodes it was fine, could you tell me what do you mean by best number of nodes and is if this skewness lower than 0.8 should be checked with my original mesh? you know I use spectral element method and it makes the mesh finer itself, so this skewness would be changed after employing it, should this limit be  checked with the original mesh or the final mesh that spectral method would apply? 

Thanks. 

firstly, grid independent test must be conducted to select the perfect grid size then go to apply the boundary conditions for solution. 

Thanks for both responses.

Dear Mehdi !

If you are interesting in alternative solution, pls, read my paper :

“T-Flow Lab” Modification of Elmer FEM Package for Vortex and Complex Flows Modeling .

I think the anisotropic model of the turbulence will provide non trivial, may be :),  solution for you.

Best regards,

Serge.

You may try to under-relax the energy and density terms. This will help to dampen any sharp increase in the bouyancy flow gradients that may be causing divergence as the flow sets up after some iterations.

I assume you want to simulate the problem in a DNS formulation.

For natural convection, the Ra number is the forcing term and it should be seen as same as the Reynolds number considered in forced convection. Therefore, you have to evaluate the cell Reynolds number based on the characteristic buoyancy velocity. If such number is O(1) everywhere, then you really have a grid resolution to perform a DNS. Otherwise, you are performing an unresolved DNS (or also denoted as LES no-model) and spectral element methods generally encounters numerical instability (at least if you do not filter out high wavenumbers components).

Furthermore, try using an implicit time integration for the buoyancy term in the momentum, that should increase the stability.

If you are not able to use a correct fine grid to perform DNS I suggest consider a LES formulation.

Could you please explain precisely "DNS simulation diverges"? Does the temperature grow unbounded at t increases? What time discretization are you using?   Are the velocity & temp. equations solved coupled & implicitly or ...? To have a chance the (spacial) mesh must be fine enough near the walls with imposed temperatures to have at least 1 mesh point within the relevant boundary layer AND the time discretization & timestep must be well chosen for the problem. [Note: these are minimal conditions for just having an approximate temperature that doesnt grow exponentially-accuracy requires more than these- it requires solve the eqns right and solve the right equations (see Prof. Denaro's astute comments)].

Thank you very much for all the responses specially Fillipo Maria Denaro  and William J Layton. 

Mr. Fillipo Maria with spectral element method how can I calculate cell Reynilds number when this method itself makes the original mesh finer, should this calculation be done for the original mesh or the final mesh results form spectral method?

Mr. William:Temperature suddenly becomes infinite at the final time step. for time discretization third order operator-splitting multi step backward differentiation is used for time integration. Equation velocity and temperature are coupled and are solved implicitly.

Finally, I add that for 3 dimensional modeling, spectral element Fourier method is used which means 2d mesh is used for the spectral plane and it is extended on the z direction by using Fourier expansion. (we didn't create any 3d mesh)

The lenght you use in the cell Reynolds number is the computational size actually used by the method. I don't think you can have a difference of one (or more) order of magnitude within the spectral element

Maybe I am missing something, but if the temp. goes from O(1) to overflow in 1  timestep,  you must be dividing by zero somewhere. Even using an unstable method would only give exponential growth.  I would suggest putting in checks for small  divisors. A related possibility  is that the problem is split so a subproblem solved at every timestep is singular due to BCs for it or a sign error. Good luck!

I agree, such a behavior in a single time-step is likely due to a bug in the code rather than to a numerical instability.

Thanks for all the responses.