About solving stationary NS eqns

Question

Hi,

I try to solve a stationary Navier-Stokes equation here but some how there is an error that Newton solver did not converge because maximum number of iterations reached.
I changed the viscosity from 1.9E-5 to 0.01 but the results is obviously wrong because the boundary condition of inflow and out flow is supposed to be 0.05m/s.

here is my code:

from dolfin import *

mesh = UnitSquareMesh(40,40);

V = VectorFunctionSpace(mesh,"CG",2)
P = FunctionSpace(mesh, "CG", 1)
W = V * P

v, q = TestFunctions(W)
up = Function(W)
u, p = split(up)
f = Constant((0.0, 0.0)) 
nu = Constant(0.01)


noslip  = DirichletBC(V, (0, 0), "on_boundary")
inlet = DirichletBC(V, Constant((0.0,0.05)), "on_boundary && (x[0] < 0.25 && x[1] == 0)")
outlet = DirichletBC(V, Constant((0.0,0.05)), "on_boundary && (x[0] > 0.75 && x[1] == 1)")

bcu = [noslip,inlet,outlet]

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

solve(F == 0, up, bcu)

plot(u)

Comments

The numerical results look correct given that the regularity requirements of your boundary conditions cannot be satisfied by the CG FE function space. You should also bear in mind that comparisons like x[1] == 0 are bad practise for floating point numbers, look into near(x[1], 0.0) instead. Also, when finding steady state solutions with viscosity as small as 1.9e-5 in the length scales of the unit square mesh, it may be worth investigating stabilised methods, using continuation steps, and so on.

---

thanks for replying. The results seem to be wrong because on the region of outflow and inflow, the velocity is supposed to be (0, 0.05) but the plot of u cannot show that.

---

As I said in my previous comment, the problem runs deeper than that. Your boundary conditions cannot be represented in your CG FE space. They are discontinuous at (x, y) = (0.25, 0.0) and (x, y) = (0.75, 1.0). As a result, you will see unphysical oscillations. Examine also your solution for p. You will probably see singularities. You should consider a smooth parabolic inlet and outlet boundary condition, close to Poiseuille flow. That would provide a much better test of your code.

---

Well , I changed the boundary defination part and it seems the code working. Thanks again!

def right(x, on_boundary):
    return x[0] > (1.0 - DOLFIN_EPS)
def left(x, on_boundary):
    return x[0] < DOLFIN_EPS
def top(x, on_boundary):
    return x[1] > 1.0 - DOLFIN_EPS and x[0] < 0.75 - DOLFIN_EPS
def bottom(x, on_boundary):
    return x[1] < DOLFIN_EPS and x[0] > 0.25 + DOLFIN_EPS
def inlet(x, on_boundary):
    return x[0] < 0.25 and x[0] > DOLFIN_EPS and x[1] < DOLFIN_EPS
def outlet(x, on_boundary):
    return x[0] > 0.75 and x[0] < 1.0 - DOLFIN_EPS and x[1] >1 - DOLFIN_EPS