[panagiat]

Context

I have made an incompressible FVM 2D Navier Stokes solver using the SIMPLE algorithm on a colocated grid. I am using Gauss seidel as the solver. I am trying to implement under-relaxation to it. My sources are Fluid Mechanics 101's video on implicit under-relaxation and Versteeg and Malalasekera's book on FVM. For my benchamrk problem I have a plain lid-driven cavity.

Under-relaxation

Regarding the under-relaxation, I am applying it:

Implicitly in the momentum equations in the following way:

• In both of my sources this equation pops

  a_P / v_relaxation_factor = Sum(a_nb * v_nb) + S + (1 - v_relaxation_factor) * (a_P / v_relaxation_factor) * v_old

  I am implementing it by
• Scaling the main diagonal momentum coefficients (a_P) by dividing them by the velocity relaxation factor right after calculating the momentum coefficients
• Adding the term (1 - v_relaxation_factor) * (a_P / v_relaxation_factor) * v_old to the source when solving the equation itself.
• Since the a_P coefficients are modified, the pressure correction equation coefficients are indirectly going to be affected by this. The same can be said for the velocity nodal/face correction values as well since they make use of the modified coefficients. Again, the only actual changes I made here are the scaling of the a_P momentum coefficients during the coeffient generation and the rest gets propagated.

Explicitly in the pressure field right after solving the pressure correction equation in the following way:

p_new = p_old + p_relaxation_factor * p_correction

Question 1

Is this the right approach to implement implicit under-relaxation to the SIMPLE algorithm? Am I missing some other stuff? I honestly don't understand why there is so much confusion regarding this topic with so much wrong information spread throughout the internet.

Confusion 1

My confusion lies in the term v_old. From my understanding of under-relaxation as a technique v_old refers to the old value of velocity, in another way, the previous value of velocity. In the context of the SIMPLE algorithm, where we have both outer and inner iterations, v_old would refer to the inner iterations since this is where the Gauss seidel takes place.

Though, by doing it this way I am diverging really fast. What I noticed is that if instead of adding the term (1 - v_relaxation_factor) * a_P * v_old to the source term when solving the equation, I add this term right after the momentum coefficient generation and thus keeping this "extra source" constant throughout the inner iterations, my code is suddenly converging. Thus, in this case, v_old refers to the starting velocity value, or in other words the value at the outer iteration right before even starting to solve the momentum equations. On a deeper search, I noticed that this is exactly what a lot of people on the internet do in their codes.

Is this indeed the correct approach and if so why? My knowledge is that under-relaxing the Gauss seidel algorithm includes using old values from the previous iterations, not just using the starting value.

Confusion 2

Even in the case I am managing to converge, by doing what I explained in the Confusion 1 part, I am getting something weird regarding mass imbalance. While monitoring the global mass imbalance (RMS of the divergence of each node before and after correcting the velocities) , it seems that it is getting lower and lower with each iteration as expected. My confusion lies on the fact that in each outer iteration when I am calculating the global mass imbalance there are always some individual cells that end up getting their mass imbalance increased. Is this normal behaviour? Are they increasing their mass imbalance for the greater good of the global mass imbalance or is this something that should not be happening at all and mass imbalance should be getting lower and lower within each outer iteration in every cell?

Thanks for your attention!