CFL condition for incompressible navier-stokes

8 months ago by
Dear all,
I was trying to make some updates on the following tutorial.
The model implements an IPCS for the Navier Stokes equation.
I changed the mesh and the boundary conditions in order to make the calculations of lift and drag of a NACA0012 airfoil.
However, when I try to make a mesh refinement on the boundary of the airfoil, the solver stops to converge. I suspect that due to the CFL condition, the smaller values of the cell h could be a problem.
Thus I wanted to implement an adaptive time-step scheme to solve the problem. However I came into two problems:
  1. Even following the CFL, the solver still do not converge (It converges however when using smaller values for dt)
  2. When I adapt dt considering the residual of the pressure equation to the half of dt's original value, the residual doubles.
For assessing the CFL condition I tried to get the minimum size of a mesh cell with:
thus, I define dt as: 
dt = mesh.hmin()/(2*u_inf)​

where u_inf is the stream velocity which is constant on my inlet profile (I am in doubt if this would be the best procedure as the velocity on the flow could be higher than this, and break the condition).

For the solver, my loop reads:

def residual(A, u, b):
    r = A*u.vector() - b
    return r.norm('l2')

def lin_solve():
    b1 = assemble(L1)
    [bc.apply(b1) for bc in dbcu]
    solve(A1, u_.vector(), b1, 'bicgstab', 'default')
    r1 = residual(A1, u_, b1)
    b2 = assemble(L2)
    [bc.apply(b2) for bc in dbcp]
    solve(A2, p_.vector(), b2, 'bicgstab', 'default')
    r2 = residual(A2, p_, b2)
    b3 = assemble(L3)
    solve(A3, u_.vector(), b3, 'cg', 'sor')
    r3 = residual(A3, u_, b3)
    return r1, r2, r3

for i in range(num_steps):
    dt = mesh.hmin()/(2*u_inf)
    # Update current time
    t += dt
    r1, r2, r3 = lin_solve()
    if r2>1:
        print('Before: '+str(r1)+ ', ' +str(r2)+ ', '+str(r3))
        r1, r2, r3 = lin_solve()
        print('After: '+ str(r1)+ ', ' +str(r2)+ ', '+str(r3))

The form of the equations given in the link are:

# 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:  Tentative velocity step
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: Pressure correction step
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: Velocity correction step
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)

Any help will be hugely appretiated!
All the best, Murilo.

Community: FEniCS Project

2 Answers

7 months ago by
Dear Murilo,

after I made the same mistake as you did on the adaptive time-step of the Navier - Stokes:

after you change your timestep (regarding any condition you like e.g. error or CFL), you must re-assemble your stiffness matrix in the tentative velocity equation. You see, inside the form F1 you have k in the first term. If you don't re-assemble you won't get different error norms, even though the timestep is different. It is like if you were solving the problem with the old time step.

Ofcourse, in the pressure and velocity correction equations, the term with k is only on the Right Hand Side of the equations, so they are correct on accident because you re-assemble b1, b2 and b3 in every time step. But in order for your problem to be correct you must re-assemble your stiffness matrix A1 everytime you adapt your timestep.

So this is what I see as incorrect. After you change this your problems should disappear.

Dear Minas,
Sorry for the long delay in replying.
But tha is it!
It is now working realy nicely!
Thank you very much for your help.
All the best, Murilo
written 7 months ago by Murilo Moreira  
8 months ago by
Hi Murilo,

what are your boundary conditions?

Dear Minas,
I have imposed a no-slip condition on the wall of the airfoil (u = (0,0,0)) and a slip-wall boundary condition on the upper and lower bound of the 2d domain. For reference this is the domain:

Thank you for your interest!
All the best, Murilo.
written 8 months ago by Murilo Moreira  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »