Dear Mgu Ant

These links might be useful, have a look:

1. https://www.mathworks.com/matlabcentral/fileexchange/68348-2d-lid-driven-cavity-flow-using-simple-algorithm

2. https://www.youtube.com/watch?v=z3OWibc-630

3. https://github.com/mathworks/2D-Lid-Driven-Cavity-Flow-Incompressible-Navier-Stokes-Solver/blob/master/docs_part1/vanilaCavityFlow_EN.md

4. https://www.youtube.com/watch?v=iyqAib3QXY8

Kind Regards

Qamar Ul Islam

I recommend that you test and compare each part of your code with manufactured solution, ie. a known solution from tha simplest case to the more complex. For example, you can begin with only the laplacian (viscous term) with dirichlet boundary conditions and compare with u1=sin(pix)sin(piy) in the domain (0,)x(0,1). Then, incorporate the time derivative and compare with u2=exp(-at)*u1, with a being a real number controlling the decay of solution. Then, you can consider steady convection and diffusion, later transient convection-diffusion to, finally, consider the full Navier-Stokes system.

Kind regards.

thank you very much for your reply Qamar Ul Islam and Roberto Carlos Cabrales .

A steady convection-diffusion equation is already working. The purpose is to have a steady SIMPLE algorithm so the transient part wont be an issue.

I have the full SIMPLE algorithm already implemented, the problem is to know where the problems are since I already used the information available on the internet to create the algorithm.

The only suspicion I have is on the Continuity and Correction equations, since the mass inbalance is almost the same after and before the correction.

The continuity equation is a laplacian equation without source from boundaries (source term = ue-uw+un-us) where sum(p') is forced to be =0

File indexed explains how these equations are written in Matlab.

Pleasure Mgu Ant Stay Happy Stay Healthy

Hi Mgu Ant, I recommend you to review your calculations with those in section 6.6 from the classical book "Numerical heat transfer and fluid flows" of Suhas Patankar. Kind regards. Roberto C.

PD. The equation to solve the pressure correction is not a Laplacian equation (for the projection method of Chorin and Temam it is). From the matricial point of view, this equation is the following product:

D*diag(Mv)^{-1}G,

where D and G are the discrete representations of the diverence and the gradient, respectively and diag(Mv) is the diagonal of the coefficient matrix of the velocities including the discretization of the viscous and convective terms and the (if any) time derivative.

I did that in cpp. It works like charm. Please make sure to create the flux values before coefficient creation in momentum equation and then correct those flux values based on the pressure correction.

In comparison to the staggered arrangement, collocated form gives faster results and improved rate of convergence. However, the first iterations are more instable.