Can anyone lend this noob a hand with his nonlinear Navier-Stokes solver?

1

I’m working on a 2D Navier-Stokes solver that I want to be able to use to model airflow past an airfoil but I’m wrestling with FEniCS a little bit. I’m trying to use a newton solver with a relaxation factor to help it converge but I don’t think this is working the way that I want it to. This code fails to converge even when parameters are chosen that should almost guarantee a stable solution. I’m still a FEniCS noob so I think there must be some stupid error that I’ve made that I’m hoping those with more experience can pick out easily.

#Test and Trial Space
w=Function(W)
w_old=Function(W)

u,p=split(w)
(u_old,p_old)=split(w_old)

# Define functions for solutions at previous and current time steps
(v,q)=TestFunctions(W)

# Define expressions used in variational forms
f = Constant((0, 0))
k = Constant(dt)
R=Constant(1./nu)
tao=Constant(R/dt)

# Define variational problem for step 1
F1= (inner(grad(u),grad(v)))*dx + tao*dot(u,v)*dx - div(v)*p*dx + R*inner(dot(u,grad(u)),v)*dx - tao*dot(u_old,v)*dx
#   a(u_,v)                                 b(v,p_)        c(u_,u_,v)
F2 = div(u)*q*dx
F=F1+F2

# Time-stepping
t = 0
eps=10**-6
maxit=500

j=derivative(F,w)

problem=NonlinearVariationalProblem(F,w,bcs,J=j)
solver=NonlinearVariationalSolver(problem)
prm=solver.parameters
prm['newton_solver']['relaxation_parameter']=0.0001

n=0

for n in range(num_steps): 
    print("Elapsed time: {}".format(n*dt))
    # Update current time
    t += dt
    n+=1
    solver.solve()

    # Update previous solution
    w_old.assign(w)
    (u_old,p_old)=w.split()
    xdmffile_u.write(u_old,t)
    xdmffile_p.write(p_old,t)```
2

One issue I see is that the convection term is incorrect. There are two different gradient operators in UFL. grad places the new index last, i.e., grad(T)_{ij\ldots k}=\partial_k T_{ij\ldots}, while nabla_grad places the new index first, i.e., nabla_grad(T)_{ij\ldots k} = \partial_iT_{j\ldots k}. dot contracts over the last index of its first argument and the first index of its second argument, i.e., dot(T,S)_{ij\ldots kl}=T_{ij\ldots m}S_{m\ldots kl}, so, you either want to replace grad with nabla_grad or switch the order of arguments to dot.