How to implement Non linear Navier Stokes with SUPG, PSPG and LSIC stabilizers?

variational formulation
1

Dear Fenics community,

I am new to fenics and PDEs overall. I am using dolphin version 2019.2.0.dev0. and I am trying to solve a simple time dependant lid-driven cavity flow problem using a non linear Navier Stokes with SUPG, PSPG and LSIC stabilizers to handle big Reynolds numbers, but I feel that the code works better without the stabilizers. I want to know if I am doing something wrong.

Thank you in advance!

from dolfin import *

# Set up mesh and finite element spaces
N = 32
mesh = UnitSquareMesh(N, N)
V = VectorElement("CG", triangle, 2)
P = FiniteElement("CG", triangle, 1)
W = FunctionSpace(mesh, V * P)

# Time-stepping parameters
dt = 0.01  # Time step size
T = 1.0    # Total simulation time

# Define boundary conditions for velocity and pressure
noslip = Constant((0, 0))
bc_wall = DirichletBC(W.sub(0), noslip, 'on_boundary && x[1] > DOLFIN_EPS')

lid_velocity = Constant((1, 0))
bc_lid = DirichletBC(W.sub(0), lid_velocity, 'near(x[1], 1)')

bc_p = DirichletBC(W.sub(1), 0, "x[0] < DOLFIN_EPS && x[1] < DOLFIN_EPS", "pointwise")

bcs = [bc_wall, bc_lid, bc_p]

# Define functions for the current and previous time step
w = Function(W)  # Current solution
w_n = Function(W)  # Previous solution

# Split mixed functions
u, p = split(w)
u_n, p_n = split(w_n)
v, q = TestFunctions(W)

# Define Reynolds number and viscosity
Reynolds = 1_000.0
nu = 1.0 / Reynolds

# Define mesh size and velocity magnitude for stabilization parameters
h = CellDiameter(mesh)
u_magnitude = sqrt(dot(u_n, u_n))

# Stabilization parameters (tau)
tau_SUPG = 1.0 / sqrt((2.0 / dt)**2 + (u_magnitude / h)**2 + (4.0 * nu / h**2)**2)
tau_PSPG = tau_SUPG  # Often, PSPG uses the same tau as SUPG
tau_LSIC = h**2 / nu  # Typically related to the mesh size and viscosity

# Define stabilization terms
# SUPG stabilization of convective term
SUPG_term = tau_SUPG * inner(dot(grad(u), u), dot(grad(v), u)) * dx

# PSPG stabilization of pressure term
PSPG_term = tau_PSPG * inner(grad(q), grad(p)) * dx

# LSIC stabilization to enforce incompressibility
LSIC_term = tau_LSIC * inner(div(u), div(v)) * dx

# Define variational problem with stabilization terms
F = (inner(u - u_n, v) / dt * dx  # Time derivative
     + nu * inner(grad(u), grad(v)) * dx  # Diffusion
     + inner(dot(grad(u), u), v) * dx  # Non-linear convection
     - div(v) * p * dx  # Pressure term
     + q * div(u) * dx  # Continuity equation
     # Add stabilization terms
     + SUPG_term  # SUPG for convective stabilization
     + PSPG_term  # PSPG for pressure stabilization
     + LSIC_term)  # LSIC for incompressibility

# Compute Jacobian
J = derivative(F, w)

# Time-stepping loop
t = 0.0
cont = 0
while t < T:
    t += dt

    # Solve for the current time step
    solve(F == 0, w, bcs, J=J)

    # Update the previous solution
    w_n.assign(w)

    cont +=1

    if ( cont % 10 == 0 ):
        print(f"Time: {t}")
        # Optionally plot the results
        plot(u, title="Velocity")
        plt.show()
        plot(p, title="Pressure")
        plt.show()
2

I see two issues:

First, PSPG is meant to enable the use of equal-order elements for velocity and pressure, but you are still using Taylor Hood elements (2nd order velocity, 1st order pressure). Either use PSPG with both elements as equal order (i.e., 1), or use Taylor Hood elements.

2nd, I don’t think your stabilization terms are correct. I believe the tau terms are ok, but I believe that for PSPG, the stabilizing term should be just tau_PSPG * grad(q), and for SUPG, it should just be tau_SUPG * dot(grad(v), u). Could be wrong.

But a bigger thing is that SUPG and PSPG are residual-based terms, and should be taken as inner products with the strong form residual, for example, with SUPG:

SUPG_term = inner( tau_SUPG * res_strong, dot(grad(v), u)) * dx

and for PSPG:

PSPG_term = inner( tau_PSPG * res_strong, grad(q)) * dx

Where res_strong is a the full residual for the Navier Stokes equations.

Since LSIC is already using the incompressibility constraint div(u) in its formulation, this should go to zero as the solution converges – so I think that is ok in your current formulation.

This is a more complicated example, but you can look at Kamensky’s repo using legacy FEniCS for a fluid-structure interaction problem that uses SUPG/PSPG/LSIC for time-dependent Navier-Stokes:

3

Hi fesolutions,

Thank you for your answer!
I apply the 2 changes you mentioned and I also found this other paper Deep reinforcement learning for the control of conjugate heat transfer, formula (13)-(14), and I think I now have more accurate results. However, it still fail for Reynolds values higher than 1_000. I would like to know if this has something to do with the implementation or if it is just that the stabilizer cannot handle bigger numbers?

Here is the new code

from dolfin import *

# Set up mesh and finite element spaces
N = 50
mesh = UnitSquareMesh(N, N)
V = VectorElement("CG", triangle, 1)
P = FiniteElement("CG", triangle, 1)
W = FunctionSpace(mesh, V * P)

# Time-stepping parameters
dt = 0.1  # Time step size
T = 1.0    # Total simulation time

# Define boundary conditions for velocity and pressure
noslip = Constant((0, 0))
bc_wall = DirichletBC(W.sub(0), noslip, 'on_boundary && x[1] > DOLFIN_EPS')

lid_velocity = Constant((1, 0))
bc_lid = DirichletBC(W.sub(0), lid_velocity, 'near(x[1], 1)')

bc_p = DirichletBC(W.sub(1), 0, "x[0] < DOLFIN_EPS && x[1] < DOLFIN_EPS", "pointwise")

bcs = [bc_wall, bc_lid, bc_p]

# Define functions for the current and previous time step
w = Function(W)  # current solution
u, p = split(w)
v, q = TestFunctions(W)

h = CellDiameter(mesh)
delta = .2*h**(3./2.) # stabilization coefficient

w_n = Function(W)  # Previous solution
u_n, p_n = split(w_n)

Reynolds = 100.0
nu = 1.0 / Reynolds

# Define the stress tensor
def sigma(u, p):
    return 2 * nu * sym(grad(u)) - p * Identity(len(u))

# Variational problem for momentum conservation
F_momentum = (inner((u - u_n), v)/dt * dx # Time derivative
              + inner(dot(u, grad(u)), v) * dx
              + div(v) * p * dx + nu * inner(grad(u), grad(v)) * dx)

# Variational form for the incompressibility
F_continuity = div(u) * q * dx

# Parameters for stabilazer
h = CellDiameter(mesh)
vnorm = sqrt(inner(u_n,u_n))
#
# Residuo of momentum
R = (u - u_n)/dt + dot(u_n,grad(u)) + grad(p)#
#PSPG
tau_pspg=((2.0/dt)**2 + (2.0*vnorm/h)**2 + 9.0*(4.0*nu/(h**2))**2)**(-0.5)
F_pspg = inner(tau_pspg*R, grad(q))*dx
#SUPG
tau_supg = ((2.0/dt)**2+(2.0*vnorm/h)**2 + 9.0*(4.0*nu/(h**2))**2)**(-0.5)
F_supg = inner(tau_supg*R, dot(u_n, grad(v)))*dx
#LSIC
tau_lsic = (vnorm*h)/2.0
F_lsic = tau_lsic*inner(div(u), div(v))*dx # div instead of nabla_div ?
# Add all the stabilyzer terms
F_momentum += F_continuity - F_pspg - F_supg - F_lsic

# Jacobian
#F1= action(F_momentum, w)
J=derivative(F_momentum, w)
# Time-stepping loop
t = 0.0
cont = 0
while t < T:
    t += dt
    print(f"Time: {t}")

    # Solve for the current time step
    solve(F_momentum==0, w, bcs, J=J)
    # Update the previous solution
    w_n.assign(w)

    cont +=1

    if ( cont % 1 == 0 ):

      print(f"Time: {t}")

      # Optionally plot the results
      plot(u, title="Velocity")
      plt.show()
      plot(p, title="Pressure")
      plt.show()
4

This is where “continuation methods” come in. Because you are solving a nonlinear system, a good initial solution (guess) will speed up convergence. If the problem is complex/highly nonlinear, you may need an initial guess to converge. Continuation methods use similar results as initial guesses for more complex problems.

So here, you could first compute a lower Reynolds number solution that does converge, and feed that solution into your solver as the starting point. You could then continue this chain of solution → next initial guess to compute higher Re flows.

Also, I’m not sure that your variational formulation is correct yet. For example, your strong form residual seems to be missing the viscous stress from the stress tensor.

5

You are using the first-order method in time. So, you need to use smaller time steps to get accurate solutions, assuming your mesh is fine enough.
For Re=1000, try with a finer mesh and either with smaller time steps or implement a second-order scheme.

I have implemented the generalised alpha scheme, which is second-order, here: GitHub - chennachaos/CFDwithFEniCSx: Computational Fluid Dynamics with FEniCSx library. This is without stabilisation, but you can add it.

Another thing is that you don’t need the Taylor-Hood element when using PSPG stabilisation. You can use equal-order spaces with PSPG stabilisation; that’s the whole idea of it.

6

Hello,

Thank you for your answers.

At the end, the simplest solution I implemented was to use a penalizer in the variational form instead of the stabilizers.

+ Constant(DOLFIN_EPS)*p*q*dx

I dont think this is the best possible solution, but at least it solves my problem.