Unrealistic results 3D Navier Stokes

Question

I want to solve a 3D (stationary) Navier Stokes problem. Since I am not very skilled in fenics I modified the 2D (instationary) example demo_navier_stokes.py. I generated the geometry and mesh with Gmsh and converted to xml format with dolfin-convert.py. Since I could not tranfer boundaries from msh to xml I had to redefine the boundaries in fenics.

For example a cylindrical part:

# class GasInletBoundary(SubDomain):
        def inside(self,x,on_boundary):
        return on_boundary and \
        x[2]<eps and math.hypot(x[0],x[1])<=d1/2
gamma_1 = GasInletBoundary()
gamma_1.mark(boundary_parts,1)
#

Etc...for a total of 8 boundaries.
There are 2 boundaries with prescribed velocities and 1 with prescribed pressure. Also many surfaces with noslip boundary conditions. I don't think that the boundary conditions are conflicting.

# boundary conditions VectorFunctionSpace V (velocity)
bc_1    = DirichletBC(V,(0,0,v_gas),boundary_parts,1)
bc_2    = DirichletBC(V,(0,0,v_air),boundary_parts,2)

# boundary conditions for (scalar) FunctionSpace Q (pressure)
bc_3    = DirichletBC(Q, p_out,boundary_parts,3)

# noslip boundary conditions VectorFunctionSpace V (velocity)
bc_4    = DirichletBC(V,(0,0,0),boundary_parts,4)
bc_5    = DirichletBC(V,(0,0,0),boundary_parts,5)

All boundaries are put in a list bcu and bcp like the example.
Solution the same as in the example.

When I choose the timing parameters such that the cycle is executed oly once like dt = 0.2 and total time T = 0.3 than the result is a reasonable solution that may or may not be right. Probably not.

If I increase the total time T and the loop is executed several times then the (some) velocities increase more or less exponentially. The final solution has a number of velocities order of magnitude 1.0E+100 or so. It looks as if in certain areas the boundary conditions get detached since when I print the velocity at the boundary I also get other values returned than the initial condition. The plot of the velocity reveals that only in certain areas the velocity has large magnitude.

There are no errors exept "Requested preconditioner is not available for uBAS Krylov solver" Using ILU.

I can send code and mesh if you contact me and/or a plot of the velocity.

Hope that somebody can help me.

Answer 1

I am not sure if have clearly understand your Problem. But you want to use an instationary solver (2D) to solve the stationary Problem (3D). Thats probably not a good idea.
I would try to use the Nonlinear variational Problem

    # build derivative
    dw = TrialFunction(self.W)
    dF = derivative(F, w, dw)

    # solve the problem
    nsproblem = NonlinearVariationalProblem(F, w, bc, dF)
    solver = NonlinearVariationalSolver(nsproblem)  
    solver.solve()

where F is your weak formulation in Mixed FunctionSpace W=V*Q. If you get no convergence you can compute a solution for a low Reynoldsnumber and take this as a initial w for the next iteration. Be sure that your mesh is fine enough.

Comments

after that you can split your solution w in velocity and pressure components

(u, p) = w.split(deepcopy=True)

---

Thank you Max for your quick reply. I understand your reservations but I could not find a stationary 3D Navier Stokes example where I only needed to change the geometry, mesh and boundary conditions. I will study your solution and modify my code (if I can). Maybe for clarification I copy some more code of my 3D NS problem. I indeed refined the mesh a couple of times but without result.

Function Spaces:

# Define function spaces (P2-P1)
V = VectorFunctionSpace(mesh, "CG", 2)
Q = FunctionSpace(mesh, "CG", 1)

# Define trial and test functions
u = TrialFunction(V)
p = TrialFunction(Q)
v = TestFunction(V)
q = TestFunction(Q)

Boundary conditions etc... I have explained in the question.

Create Functions:

# Create functions
 u0 = Function(V)
 u1 = Function(V)
 p1 = Function(Q)

Then the solution is more or less like in the example:

# Tentative velocity step
F1 = (1/k)*inner(u - u0, v)*dx + inner(grad(u0)*u0, v)*dx + \
nu*inner(grad(u), grad(v))*dx - inner(f, v)*dx
a1 = lhs(F1)
L1 = rhs(F1)

# Pressure update
a2 = inner(grad(p), grad(q))*dx
L2 = -(1/k)*div(u1)*q*dx

 # Velocity update
 a3 = inner(u, v)*dx
 L3 = inner(u1, v)*dx - k*inner(grad(p1), v)*dx

# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)

# Time-stepping
t = dt
while t < T + DOLFIN_EPS:

# Compute tentative velocity step
begin("Computing tentative velocity")
b1 = assemble(L1)
[bc.apply(A1, b1) for bc in bcu]
solve(A1, u1.vector(), b1, "gmres", "default")
end()

# Pressure correction
begin("Computing pressure correction")
b2 = assemble(L2)
[bc.apply(A2, b2) for bc in bcp]
solve(A2, p1.vector(), b2, "gmres", "default")
end()

# Velocity correction
begin("Computing velocity correction")
b3 = assemble(L3)
[bc.apply(A3, b3) for bc in bcu]
solve(A3, u1.vector(), b3, "gmres", "default")
end()