3D navier-stokes equation with divergent velocity

variational formulation
,
1

Good morning everyone ! I hope I posted the question in the right section. I’m a beginner with FEniCS and I’m trying to make this program work which should study the navier stokes equation by applying it to the case of a box with a cylinder inside. When I start the program the speed diverges and I can’t understand why. Below I paste the program with which I created the 3D mesh and the one that should solve the equation.

Thanks to anyone who wants to help me and sorry for my bad English !

MESH-------------------------------

import pygmsh
import gmsh

mesh_size = 1000

geom = pygmsh.occ.Geometry()
model3D = geom.__enter__()

box0 = model3D.add_box([0.0, 0, 0], [5.0, 1.0, 1.0])

center = [2.0, 0.0, 0.5]  
direction = [0.0, 1.0, 0.0]  
radius = 0.2  
height = 1.0  

cylinder = model3D.add_cylinder(center, direction, radius)

union = model3D.boolean_difference([box0], [cylinder])

model3D.synchronize()

model3D.add_physical(union, "Box minus Ball")

gmsh.initialize()
gmsh.option.setNumber('General.Terminal', 2)  # Opzionale: permette la visualizzazione nel terminale
gmsh.model.geo.synchronize()
gmsh.fltk.run()  

geom.generate_mesh(dim=3)

gmsh.write("mesh3D.msh")

import meshio
mesh = meshio.read("mesh3D.msh")
meshio.write("mesh3D.xml", mesh)

model3D.__exit__()

NAVIER STOKES EQ---------------------

from fenics import *
from fenics import plot
import numpy as np
import matplotlib.pyplot as plt
from fenics import set_log_level

T = 5.0  
num_steps = 500  
dt = T / num_steps  
mu = 0.001
rho = 1

mesh_file = "PycharmProjects/pythonProject1/poisson/mesh3D.xml"  
mesh = Mesh(mesh_file)

V = VectorFunctionSpace(mesh, 'Lagrange', 2)
Q = FunctionSpace(mesh, 'Lagrange', 1)

inflow = 'near(x[0], 0) && (x[2] >= 0.0) && (x[2] <= 1.0) && (x[1] >= 0.0) && (x[1] <= 1.0)'

outflow = 'near(x[0], 5.0) && (x[2] >= 0.0) && (x[2] <= 1.0) && (x[1] >= 0.0) && (x[1] <= 1.0)'

walls = 'near(x[1], 0) || near(x[1], 1.0) || near(x[2], 0) || near(x[2], 1.0)'


cylinder = 'on_boundary && exp(-pow(x[0] - 2.0, 2) - pow(x[1] - 0.0, 2)) > 1e-4 && x[2] > 0.0 && x[2] < 1.0'

inflow_profile = ('4.0*1.5*x[1]*(0.41 - x[1]) / pow(0.41, 2)', '0.0', '0.0')

bcu_inflow = DirichletBC(V, Expression(inflow_profile, degree=3), inflow)
bcu_walls = DirichletBC(V, Constant((0, 0, 0)), walls)
bcu_cylinder = DirichletBC(V, Constant((0, 0, 0)), cylinder)
bcp_outflow = DirichletBC(Q, Constant(0), outflow)
bcu = [bcu_inflow, bcu_walls, bcu_cylinder]
bcp = [bcp_outflow]

u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)

u_n = Function(V)
u_  = Function(V)
p_n = Function(Q)
p_  = Function(Q)

U  = 0.5*(u_n + u)
n  = FacetNormal(mesh)
f  = Constant((0, 0, 0))
k  = Constant(dt)
mu = Constant(mu)
rho = Constant(rho)

def epsilon(u):
    return sym(nabla_grad(u))


def sigma(u, p):
    return 2*mu*epsilon(u) - p*Identity(len(u))

F1 = rho*dot((u - u_n) / k, v)*dx \
   + rho*dot(dot(u_n, grad(u_n)), v)*dx \
   + inner(sigma(U, p_n), epsilon(v))*dx \
   + dot(p_n*n, v)*ds - dot(mu*grad(U)*n, v)*ds \
   - dot(f, v)*dx
a1 = lhs(F1)
L1 = rhs(F1)

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

a3 = dot(u, v)*dx
L3 = dot(u_, v) * dx - k * dot(nabla_grad(p_ - p_n), v) * dx

A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)

[bc.apply(A1) for bc in bcu]
[bc.apply(A2) for bc in bcp]

xdmffile_u = XDMFFile('navier_stokes_cylinder/velocity3D.xdmf')
xdmffile_p = XDMFFile('navier_stokes_cylinder/pressure3D.xdmf')

timeseries_u = TimeSeries('navier_stokes_cylinder/velocity_series')
timeseries_p = TimeSeries('navier_stokes_cylinder/pressure_series')

File('navier_stokes_cylinder/cylinder3D.xml.gz') << mesh

progress = Progress('Time-stepping')
set_log_level(20)

t = 0
for n in range(num_steps):

    t += dt

    b1 = assemble(L1)
    [bc.apply(b1) for bc in bcu]
    solve(A1, u_.vector(), b1, 'cg', 'amg')
    b2 = assemble(L2)
    [bc.apply(b2) for bc in bcp]
    solve(A2, p_.vector(), b2, 'cg', 'ilu')
 
    b3 = assemble(L3)
    solve(A3, u_.vector(), b3, 'cg', 'ilu')
  
    xdmffile_u.write(u_, t)
    xdmffile_p.write(p_, t)

    timeseries_u.store(u_.vector(), t)
    timeseries_p.store(p_.vector(), t)

    u_n.assign(u_)
    p_n.assign(p_)

    print('u max:', u_.vector().get_local().max())

MESSAGE ERROR--------------

u max: 2.7359762055116166
u max: 4.654361694080788
u max: 5.447495587843657
u max: 6.6365560406239155
u max: 9.366411572418862
u max: 11.650689977047024
u max: 13.358632734799954
u max: 14.56366184477287
u max: 15.466757431052363
u max: 16.298928921488013
u max: 17.244402998964883
u max: 18.398908770510193
u max: 19.76655044863949
u max: 21.298394084285274
u max: 57.56885650832666
u max: 648.3555059972953
u max: 72366.29021211487
u max: 818675553.0731914
u max: 9.920378582775226e+16
u max: 1.433316426121824e+33
u max: 2.99352001478459e+65
u max: 1.3244332715306184e+130
Traceback (most recent call last):
  File "PycharmProjects/pythonProject1/poisson/EQplot.py", line 306, in <module>
    solve(A1, u_.vector(), b1, 'cg', 'amg')
  File "/Users/muz/anaconda3/envs/fenics-env/lib/python3.8/site-packages/dolfin/fem/solving.py", line 227, in solve
    return dolfin.la.solver.solve(*args)
  File "/Users/muz/anaconda3/envs/fenics-env/lib/python3.8/site-packages/dolfin/la/solver.py", line 72, in solve
    return cpp.la.solve(A, x, b, method, preconditioner)
RuntimeError:
2

Looking at F1, the boundary integrals seem off to me. Usually, you substitute something in place of those integrals to enforce a condition on the traction (or v = 0 when you prescribe a Dirichlet condition on u). That could be a major source of instability. What happens when you remove those terms?

( Along the same line, a Dirichlet condition on the pressure is not a well-posed condition for the Navier-Stokes equation with velocity and pressure solution variables. Usually you can get away with this though, so I don’t expect this to be the perpetrator. )

Where exactly did you get this formulation from? Some background on the equations that you are trying to implement would be important to be able to help you effectively.

3

I may be wrong, but I think this is Chorin’s projection method. There’s a section dedicated to it in @dokken’s tutorial. This legacy DOLFIN code likely derives from an old tutorial, e.g., here.

4

And a better splitting scheme is described here:
https://jsdokken.com/dolfinx-tutorial/chapter2/ns_code2.html
with proper mathematical notes at:
https://computationalphysiology.github.io/oasisx/splitting_schemes.html

1 Like
6

I’m actually trying to apply Chorin’s projection method. I had found a 2D old example and i tried to modify it in 3D but with no success. Unfortunately I can’t find that exact tutorial. The things I changed are the mesh (originally created with mshr) and the part where I define inflow, outflow and walls trying to adapt them to the 3D case.

7

Yes exactly ! I noticed that more up-to-date tutorials, like the one you recommended, work with fenicsx. Do you recommend working with fenicsx instead of fenics? Sorry if I ask stupid questions but I’m really a beginner.

8

Yes, I do recommend fenicsx over fenics.
See:

for context

9

Thank you very much, I will try fenicsx. Thanks again everyone for your replies!

10

Did you try removing the terms that I mentioned?

All the sources given by @dokken and @nate on Chorins method do not have the boundary terms dot(p_n*n, v)*ds - dot(mu*grad(U)*n, v)*ds in the predictor step (your F1).

For 99% of cases, those integrals should not be there. Unless you have a reason to add them that I’m not seeing, they make your problem ill posed. Hence my earlier question on where you got your variational formulation from.

11

Yes I tried to do as you say and that piece of F1 code was effectively useless. Even if you eliminate it, however, the error remains the same. Another thing I did was plot this navier_stokes_cylinder/velocity3D.xdmf with pyvista. But this comes out. Could it be useful to understand what’s wrong?

Screenshot 2024-11-22 alle 17.54.22|690x474

12

A screenshot without corresponding code doesn’t help

13

I’m sorry! pyvista for my convenience I preferred to start it directly from the terminal using this

python -c "import pyvista as pv; pv.start_xvfb(); mesh = pv.read('velocity3D.xdmf'); mesh.plot()"
14

I would then start with checking that the data is visualized correctly in Paraview.

This also seems like a pure Paraview question, as you are asking how to visualize velocity data from an XDMFFile with Paraview. Such questions are best asked at: pyvista/pyvista · Discussions · GitHub

15

Ok I will try to check where you suggest, thanks again for your help