CFL condition for incompressible navier-stokes
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:
- Even following the CFL, the solver still do not converge (It converges however when using smaller values for dt)
- When I adapt dt considering the residual of the pressure equation to the half of dt's original value, the residual doubles.
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) k.assign(dt) # Update current time t += dt r1, r2, r3 = lin_solve() refi_r1.append(r1) refi_r2.append(r2) refi_r3.append(r3) if r2>1: print('Before: '+str(r1)+ ', ' +str(r2)+ ', '+str(r3)) k.assign(dt/2) 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.
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.
what are your boundary conditions?