Navier Stokes Problem (convergence error)

Hi I’m new to fenics and working through the fenics tutorial to solve the Navier-Stokes equations.
I want to solve a problem for flow around a cylinder. As is known, the flow around a cylinder depends on the Reynolds number:
I want to increase the speed, for example, instead of 1.5 to 4, to get a Karman vortex street, turbulent mode, but I get a convergence error. How to fix it?

CODE:

from __future__ import print_function

from fenics import *
from mshr import *
import numpy as np

T = 5.0            # final time
num_steps = 5000   # number of time steps
dt = T / num_steps # time step size
mu = 0.001         # dynamic viscosity
rho = 1            # density

# Create mesh
channel = Rectangle(Point(0, 0), Point(2.2, 0.41))
cylinder = Circle(Point(0.2, 0.2), 0.05)
domain = channel - cylinder
mesh = generate_mesh(domain, 64)

# Define function spaces
V = VectorFunctionSpace(mesh, 'P', 2)
Q = FunctionSpace(mesh, 'P', 1)

# Define boundaries
inflow   = 'near(x[0], 0)'
outflow  = 'near(x[0], 2.2)'
walls    = 'near(x[1], 0) || near(x[1], 0.41)'
cylinder = 'on_boundary && x[0]>0.1 && x[0]<0.3 && x[1]>0.1 && x[1]<0.3'

# Define inflow profile
inflow_profile = ('4.0*1.0*x[1]*(0.41 - x[1]) / pow(0.41, 2)', '0')

# Define boundary conditions
bcu_inflow = DirichletBC(V, Expression(inflow_profile, degree=2), inflow)
bcu_walls = DirichletBC(V, Constant((0, 0)), walls)
bcu_cylinder = DirichletBC(V, Constant((0, 0)), cylinder)
bcp_outflow = DirichletBC(Q, Constant(0), outflow)
bcu = [bcu_inflow, bcu_walls, bcu_cylinder]
bcp = [bcp_outflow]

# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)

# Define functions for solutions at previous and current time steps
u_n = Function(V)
u_  = Function(V)
p_n = Function(Q)
p_  = Function(Q)

# Define expressions used in variational forms
U  = 0.5*(u_n + u)
n  = FacetNormal(mesh)
f  = Constant((0, 0))
k  = Constant(dt)
mu = Constant(mu)
rho = Constant(rho)

# Define symmetric gradient
def epsilon(u):
    return sym(nabla_grad(u))

# Define stress tensor
def sigma(u, p):
    return 2*mu*epsilon(u) - p*Identity(len(u))

# Define variational problem for step 1
F1 = rho*dot((u - u_n) / k, v)*dx \
   + rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx \
   + inner(sigma(U, p_n), epsilon(v))*dx \
   + dot(p_n*n, v)*ds - dot(mu*nabla_grad(U)*n, v)*ds \
   - dot(f, v)*dx
a1 = lhs(F1)
L1 = rhs(F1)

# Define variational problem for step 2
a2 = dot(nabla_grad(p), nabla_grad(q))*dx
L2 = dot(nabla_grad(p_n), nabla_grad(q))*dx - (1/k)*div(u_)*q*dx

# Define variational problem for step 3
a3 = dot(u, v)*dx
L3 = dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_n), v)*dx

# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)

# Apply boundary conditions to matrices
[bc.apply(A1) for bc in bcu]
[bc.apply(A2) for bc in bcp]

# Create XDMF files for visualization output
xdmffile_u = XDMFFile('navier_stokes_cylinder_4/velocity.xdmf')
xdmffile_p = XDMFFile('navier_stokes_cylinder_4/pressure.xdmf')

# Create time series (for use in reaction_system.py)
timeseries_u = TimeSeries('navier_stokes_cylinder_4/velocity_series')
timeseries_p = TimeSeries('navier_stokes_cylinder_4/pressure_series')

# Save mesh to file (for use in reaction_system.py)
File('navier_stokes_cylinder_4/cylinder.xml.gz') << mesh

# Create progress bar
progress = Progress('Time-stepping', num_steps)
# set_log_level(LogLevel.PROGRESS)

# Time-stepping
t = 0
for n in range(num_steps):

    # Update current time
    t += dt

    # Step 1: Tentative velocity step
    b1 = assemble(L1)
    [bc.apply(b1) for bc in bcu]
    solve(A1, u_.vector(), b1, 'bicgstab', 'hypre_amg')

    # Step 2: Pressure correction step
    b2 = assemble(L2)
    [bc.apply(b2) for bc in bcp]
    solve(A2, p_.vector(), b2, 'bicgstab', 'hypre_amg')

    # Step 3: Velocity correction step
    b3 = assemble(L3)
    solve(A3, u_.vector(), b3, 'cg', 'sor')

    # Plot solution
    plot(u_, title='Velocity')
    plot(p_, title='Pressure')

    # Save solution to file (XDMF/HDF5)
    xdmffile_u.write(u_, t)
    xdmffile_p.write(p_, t)

    # Save nodal values to file
    timeseries_u.store(u_.vector(), t)
    timeseries_p.store(p_.vector(), t)

    # Update previous solution
    u_n.assign(u_)
    p_n.assign(p_)

    # Update progress bar
    # progress.update(t / T)
    set_log_level(LogLevel.PROGRESS)
    progress += 1
    set_log_level(LogLevel.ERROR)
    vertex_values_u_ = u_.compute_vertex_values(mesh)
    # print('u max:', np.max(vertex_values_u_))
# Hold plot
interactive()

Please note that your scheme has to satisfy the CFL condition, as mentioned in: Test problem 2: Flow past a cylinder (DFG 2D-3 benchmark) — FEniCSx tutorial
referencing: Courant–Friedrichs–Lewy condition - Wikipedia
which gives a relationship between the fluid velocity, the spatial and temporal discretization.

If I counted correctly, then it seems like everything satisfies

How can I increase the speed without violating the condition CFL?

Decreasing the time step.
Have you looked at the solution (or does it fail at the first time step)?
If you look at the solution up to the point of non-convergence, you can usually see the unstability forming.

Okay, I’ll try to reduce the time step

increased the number of steps for speed 6 to 35000 and still there is a convergence error

As I asked in the previous post, is the instability at the first time step, or does it happen after a while?

If it happens after some time steps, how does the velocity and pressure look prior to the non-convergence?

using paraview to show you the speed and pressure before the discrepancy?

Yes, that is what Im asking for. Does the solution look sensible after the first time step? If not, can you see an instability?

You should also check if your problem is solved by using direct solvers, as opposed to iterative solvers.


honestly I don’t understand what is not stable and how to fix it

Here the speed is set to 6

and what is a direct solver?

See for instance: How to choose the optimal solver for a PDE problem? - #2 by nate
or for even more details: KSP: Linear System Solvers — PETSc v3.17.1-435-gd0fabe2265 documentation

I would also suggest to compute the Reynolds number of your problem, as described in 4.2 of: Advanced partial differential equation models | SpringerLink
and consider:
A FEniCS-based programming framework for modeling turbulent flow by the Reynolds-averaged Navier–Stokes equations - ScienceDirect
and the tutorial I posted in the previous post, that uses a slightly more advanced splitting scheme.

1 Like