Implementing Penalty FEM for Lid-Driven Cavity Problem

I am using FEATool Multiphysics for a lid-driven cavity problem and struggling to implement a penalty FEM formulation that responds correctly to changes in the penalty constant, 𝛾. My goal is to enforce incompressibility, but increasing 𝛾 does not improve the results or lead to convergence.

Here’s a brief on my approach:

I modified the momentum equations to include a penalty term linked to the velocity divergence.
Despite these changes, the results do not significantly improve with higher 𝛾 values, and the model does not converge.
Here are the equations I used:

rho_ns*u' - (miu_ns + miuT_ns)*(2*ux_x + uy_y + vx_y) + rho_ns*(u*ux + v*uy) - gamaa*(ux_x + vy_x) = Fx_ns

rho_ns*v' - (miu_ns + miuT_ns)*(vx_x + uy_x + 2*vy_y) + rho_ns*(u*vx + v*vy) - gamaa*(ux_y + vy_y) = Fy_ns

Any help or pointers would be greatly appreciated. Thank you in advance for your support!

Hi, I believe the sign for the penalty term is inverted, and from my understanding of the Navier-Stokes penalty method is that you use it to avoid solving for the pressure and continuity equation, so it is better to use a custom equation so as not to have the extra pressure dependent variable. See the attached model script file (note that it is very sensitive to mesh configuration and element type, so probably would need more stabilization).

navier-stokes-penalty-method.fes

Thank you for your reply and for providing the model script. I have tried different formulations and observed that the penalty formulations are highly sensitive to the element type. Specifically, linear elements do not converge to the normal Navier-Stokes results, but quadratic elements yield acceptable results.

Additionally, I found that formulations with either a negative or positive sign for the penalty term produce similar results at high values of gamma. However, at low values of gamma, the formulation with a negative sign provides more rational results.

Thank you again for your support and insights.