Navier-Stokes for cylinder flow does not converge/Newton method

Errors
1

Hi, I wrote a program to solve the flow problem around a cylinder using Newton method. It works but does not converge. Hope to receive help, thanks.

My code:

from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt
from mshr import *
import numpy as np

T = 2.0            # final time
num_steps = 2000   # number of time steps
dt = T / num_steps # time step size
mu = 0.0001         # dynamic viscosity
rho = 1            # density
Re=100              # reynold number

# 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)
elem01 = MixedElement([V.ufl_element(), Q.ufl_element()])
W = FunctionSpace(mesh, elem01)

# 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.15 && x[0]<0.25 && x[1]>0.15 && x[1]<0.25'

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

# Define boundary conditions
bcu_inflow = DirichletBC(W.sub(0), Expression(inflow_profile, degree=2), inflow)
bcu_walls = DirichletBC(W.sub(0), Constant((0, 0)), walls)
bcu_cylinder = DirichletBC(W.sub(0), Constant((0, 0)), cylinder)
bcp_outflow = DirichletBC(W.sub(1), Constant(0), outflow)
bcp= [bcu_inflow, bcu_walls, bcu_cylinder]
bcu = [bcp_outflow]
bcs = [bcu_inflow, bcu_walls, bcu_cylinder, bcp_outflow]

dw = TrialFunction(W)   
v, q = TestFunctions(W)
w = Function(W)
w_n = Function(W)
w_fem = Function(W)

du, dp = split(dw)
u, p = split(w)
u_n, p_n = split(w_n)
u_fem, p_fem = split(w_fem)

class InitialCondition(UserExpression):
    def eval(self, values, x):
 
        if near(x[0], 0):  
            values[0] = 4.0 * 1.5 * x[1] * (0.41 - x[1]) / pow(0.41, 2)  # u
            values[1] = 0  # v
        else:
            values[0] = 0  # u
            values[1] = 0  # v
        values[2] = 0  

    def value_shape(self):
        return (3,)  

initial_condition = InitialCondition(degree=2)
w_n.interpolate(initial_condition)
w.interpolate(initial_condition)

# Netwon solver
class NS(NonlinearProblem):
    def __init__(self, a, L, bcs):
        NonlinearProblem.__init__(self)
        self.a = a
        self.L = L
        self.bcs = bcs

    def F(self, b, x):
        assemble(self.L, tensor=b)
        for bc in self.bcs:  #
            bc.apply(b)
    def J(self, A, x):
        assemble(self.a, tensor=A)
        for bc in self.bcs:  # 
            bc.apply(A)
#weak form


L1 = Re * dot((u - u_n) / dt, v) * dx \
   + Re * dot(nabla_grad(u_n).T*u, v)*dx \
   + inner(nabla_grad(u), nabla_grad(v))*dx \
   - p*div(v)*dx

L2 = div(u) * q * dx

L=L1+L2
a = derivative(L, w, dw)

problem = NS(a, L, bcs)
solver = NewtonSolver()
solver.parameters["linear_solver"] = "lu"
solver.parameters["convergence_criterion"] = "incremental"
solver.parameters["relative_tolerance"] = 1e-6

file = File("NS/output.pvd", "compressed")

# Time-stepping
t=0.0
for n in range(5000):

    t += dt
    w_n.vector()[:] = w.vector()
    solver.solve(problem, w.vector())
    u, p = w.split()
    if n % 1 == 0:
        file << (u, t)

    plot(u)
    plt.draw()
    plt.pause(0.02)
    plt.clf()
2

I believe your boundary conditions are applied wrongly in non-linear problem.

You can start by verifying this claim by using:

    solve(L == 0, w, bcs)

which is inefficient, but shows that the problem should converge.

Why don’t you just use NonlinearVariationalProblem instead of creating your own assembler?
See for instance:
Switching from Fenics to FenicsX for Incompressible Mixed Element Code for an example