[MrPGW]

Dear colleagues,
I have just started my CFD journey and am trying to solve a cavity flow problem using SIMPLE-algorithm on a collocated mesh. I have been following the CFD lecture series of Dr. Mazumder on YouTube (https://youtube.com/playlist?list=PL...-3xQDgbB1y9z_T).
Now for my problem:
I have set up all the needed functions (momentum link coefficient and source-term computing, solver for the momentum equations, face velocity computation using Rhie-Chow interpolation, pressure link coefficient computation and pressure-solver as well as pressure-, cell-center-velocity- and face-velocity-corrections). On the first outer iteration, after initializing u, v and p to zero, everything works fine. The solution of the momentum-equation results in a velocity field that is to be expected for a diffusion-only problem (since u and v and therefore the convective terms are zero).
u:
[[0.429717 0.713736 0.770268 0.731336 0.483334]
[0.151675 0.324959 0.389460 0.353109 0.168665]
[0.059850 0.141998 0.183331 0.163704 0.068440]
[0.023683 0.057999 0.078984 0.070698 0.028231]
[0.006708 0.016564 0.023047 0.020666 0.008150]]


The subsequent further steps also yield results that seem plausible.


On the second outer iteration however, the solution for the momentum-equations returns a velocity field that doesn’t seem right and leads to a blow-up of the solution in further iterations. Since the solver works fine for a diffusion-only problem, I’m suspecting that the problem lies in the momentum links and/or source terms. As soon as I introduce non-zero-velocities and the convective terms are activated, the solution becomes unrealistic. I checked the momentum links and source terms again and again, but can’t seem to find the problem. I have also played around with other parameters such as lid-velocity, grid spacing, lengths of the solution domain etc. For demonstration purposes: when I set the initial u-velocity to 1 (except at the boundary faces), the solver returns the following velocity-field:
u:
[[0.001329 0.002657 0.003983 0.006284 1.973763]
[0.000001 0.000002 0.000004 0.000206 0.410125]
[0.000000 0.000000 0.000000 0.000040 0.084236]
[0.000000 0.000000 0.000000 0.000008 0.017135]
[0.000000 0.000000 0.000000 0.000001 0.003217]]


I have included the python-function to calculate the links and sources and would greatly appreciate if someone more experienced could have a look at the calculation or point out where else the problem might lie.