Lid driven cavity not converging for Reynolds > 600

1

Hi all,

I’m trying to solve the lid driven cavity problem in steady state for high Reynolds numbers. I didn’t manage to get the solver converging for Re > 600 even with more refined meshes.

Does any one have an idea on how to get it converging ?
Do I need to impose the pressure somehow ?

Thanks in advance,

Rem

from fenics import *

# ##### Mesh parameters ##### #

mesh = UnitSquareMesh(40, 40, diagonal="crossed")

# ##### Define function spaces ##### #

V_ele = VectorElement("CG", mesh.ufl_cell(), 2)
Q_ele = FiniteElement("CG", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, MixedElement([V_ele, Q_ele]))

# ##### Boundary conditions ##### #
lid_location = "near(x[1],  1.)"
fixed_wall_locations = "near(x[0], 0.) | near(x[0], 1.) | near(x[1], 0.)"

noslip = DirichletBC(W.sub(0), (0, 0), fixed_wall_locations)

u_x = 1
velocity_x = Constant((u_x, 0.0))
top_velocity = DirichletBC(W.sub(0), velocity_x,
                           lid_location)

bcu = [top_velocity, noslip]

# ##### Define variational parameters ##### #

v, q = TestFunctions(W)
up = Function(W)
u, p = split(up)

f = Constant((0, 0))

rho_0 = 1
mu = 1
u_x = 1
Reynolds = rho_0*u_x/mu

velocity_x.assign(Constant((u_x, 0)))

# Solver

F = (
     rho_0*inner(grad(u)*u, v)*dx + mu*inner(grad(u), grad(v))*dx
     - inner(p, div(v))*dx + inner(q, div(u))*dx
     + inner(f, v)*dx
)

F += inner(q, div(u))*dx

solve(F == 0, up, bcu)

u, p = up.split()
XDMFFile("Re={:.1f}_u.xdmf".format(Reynolds)).write(u)
XDMFFile("Re={:.1f}_p.xdmf".format(Reynolds)).write(p)
2

There are several points to consider:

  • First, the Bubnov–Galerkin formulation is inherently unstable at a high element Reynolds number (i.e., the Reynolds number taking the element diameter as a length scale), so it will become impractical to use at a high global Reynolds number (which would force a very small element size for stability). This affects the accuracy of the solution even if you can coax the nonlinear solver into converging, but it can also make the equation system harder to solve. The remedy for this is some form of stabilization, e.g., SUPG/PSPG, also considered a type of variational multiscale (VMS) method. (See Chapter 6 of my course notes.)

  • Second, the algebraic residual is difficult to converge, even for stabilized formulations! The usual checklist of things to try is:

    • Use a line search scheme, which can be accessed by using a PETSc SNES through a NonlinearVariationalSolver (as demonstrated in FEniCS here).
    • Use a continuation method, e.g., solve once with a low Reynolds number, use that solution as the initial guess for a slightly higher Reynolds number, etc.
    • If nothing else works, solve an unsteady problem and try to drive it to a steady state. (However, this may not always work well for Navier–Stokes, because if turbulence develops it may never settle into a steady solution.)
3 Likes
3

Thank you for the very useful feedback!

  • For some reason using a line search scheme with the provided method didn’t work.
    from fenics import *

    # ##### Mesh parameters ##### #

    mesh = UnitSquareMesh(40, 40, diagonal="crossed")

    # ##### Define function spaces ##### #

    V_ele = VectorElement("CG", mesh.ufl_cell(), 2)
    Q_ele = FiniteElement("CG", mesh.ufl_cell(), 1)
    W = FunctionSpace(mesh, MixedElement([V_ele, Q_ele]))

    # ##### Boundary conditions ##### #
    lid_location = "near(x[1],  1.)"
    fixed_wall_locations = "near(x[0], 0.) | near(x[0], 1.) | near(x[1], 0.)"

    noslip = DirichletBC(W.sub(0), (0, 0), fixed_wall_locations)

    u_x = 1
    velocity_x = Constant((u_x, 0.0))
    top_velocity = DirichletBC(W.sub(0), velocity_x,
                               lid_location)

    bcu = [top_velocity, noslip]

    # ##### Define variational parameters ##### #

    v, q = TestFunctions(W)
    up = Function(W)
    u, p = split(up)

    f = Constant((0, 0))

    rho_0 = 1
    mu = 1
    u_x = 1
    Reynolds = rho_0*u_x/mu

    velocity_x.assign(Constant((u_x, 0)))

    # Solver

    F = (
         rho_0*inner(grad(u)*u, v)*dx + mu*inner(grad(u), grad(v))*dx
         - inner(p, div(v))*dx + inner(q, div(u))*dx
         + inner(f, v)*dx
    )

    F += inner(q, div(u))*dx

    JF = derivative(F, up, TrialFunction(W))
    problem = NonlinearVariationalProblem(F, up, bcu, JF)
    solver = NonlinearVariationalSolver(problem)
    print(solver.parameters.keys())
    solver.parameters['nonlinear_solver'] = 'snes'
    solver.parameters['snes_solver']['line_search'] = 'bt'

    solver.solve()

    u, p = up.split()
    XDMFFile("Re={:.1f}_u.xdmf".format(Reynolds)).write(u)
    XDMFFile("Re={:.1f}_p.xdmf".format(Reynolds)).write(p)
  • However, we’ve managed to implement a continuation method solving low Re numbers first and using them as a first guess. Worked like a charm! This is the code we’ve come up with:
from fenics import *

# ##### Mesh parameters ##### #
N = 128
mesh = UnitSquareMesh(N, N, diagonal="crossed")

# ##### Define function spaces ##### #

V_ele = VectorElement("CG", mesh.ufl_cell(), 2)
Q_ele = FiniteElement("CG", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, MixedElement([V_ele, Q_ele]))

# ##### Boundary conditions ##### #
lid_location = "near(x[1],  1.)"
fixed_wall_locations = "near(x[0], 0.) | near(x[0], 1.) | near(x[1], 0.)"

noslip = DirichletBC(W.sub(0), (0, 0), fixed_wall_locations)

u_x = 1
velocity_x = Constant((u_x, 0.0))
top_velocity = DirichletBC(W.sub(0), velocity_x,
                           lid_location)

bcu = [top_velocity, noslip]

# ##### Define variational parameters ##### #

v, q = TestFunctions(W)
up = Function(W)
u_export = Function(W)
up_old = Function(W)
u, p = split(up)

f = Constant((0, 0))


for top_velocity in [600, 1000, 2000, 4000, 4500, 5000, 7000, 9000, 10000]:

    up.assign(up_old)
    rho_0 = 1
    mu = 1
    # u_x = 1
    Reynolds = rho_0*top_velocity/mu
    print("Doing it for Re={}".format(Reynolds))
    velocity_x.assign(Constant((top_velocity, 0)))

    # Solver

    F = (
        rho_0*inner(grad(u)*u, v)*dx + mu*inner(grad(u), grad(v))*dx
        - inner(p, div(v))*dx + inner(q, div(u))*dx
        + inner(f, v)*dx
    )

    F += inner(q, div(u))*dx

    solve(F == 0, up, bcu)

    up_old.assign(up)
    u_export.assign(up)
    u_export, p_export = u_export.split()
    XDMFFile("Results/Re={:.1f}_u.xdmf".format(Reynolds)).write(u_export)
    XDMFFile("Results/Re={:.1f}_p.xdmf".format(Reynolds)).write(p_export)
1 Like