Navier Stokes equiation 3D

[Sebastian]

Hi, I’m new to fenics and have been working with the fenics tutorial to solve the Navier Stokes equations in 2D.
I want to solve the problem in 3D for the flow around a cylinder. I have modified the tutorial code but when I run it the speed values do not evolve. Any idea how to solve it? Some error in my code?

from __future__ import print_function
from fenics import *
from mshr import *
import numpy as np

T = 5.0            # final time
num_steps = 5000   # number of time steps
dt = T / num_steps # time step size
mu = 0.001         # dynamic viscosity
rho = 1            # density

# Create mesh
channel = Box(Point(0, 0, 0), Point(2.2, 0.41, 0.41))
cylinder = Cylinder(Point(0.2, 0.2, 0.4), Point(0.2, 0.2, 0.01), 0.05, 0.05, 32)
domain = channel - cylinder
mesh = generate_mesh(domain, 64)

# Define function spaces
V = VectorFunctionSpace(mesh, 'P', 2)
Q = FunctionSpace(mesh, 'P', 1)

# Define boundaries
inflow   = 'near(x[0], 0)'
outflow  = 'near(x[0], 2.2)'
walls    = 'near(x[1], 0) || near(x[1], 0.41) || near(x[2], 0) || near(x[2], 0.41)'
cylinder = 'on_boundary && x[0]>0.1 && x[0]<0.3 && x[1]>0.1 && x[1]<0.3'

# Define inflow profile
inflow_profile = ('16.0*1.5*x[1]*x[2]*(0.41 - x[1])*(0.41 - x[2]) / pow(0.41, 4)', '0', '0')


# Define boundary conditions
bcu_inflow = DirichletBC(V, Expression(inflow_profile, degree=2), 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]

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

# Define functions for solutions at previous and current time steps
u_n = Function(V)
u_  = Function(V)
p_n = Function(Q)
p_  = Function(Q)

# Define expressions used in variational forms
U  = 0.5*(u_n + u)
n  = FacetNormal(mesh)
f  = Constant((0, 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
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
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
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)

# Apply boundary conditions to matrices
[bc.apply(A1) for bc in bcu]
[bc.apply(A2) for bc in bcp]

# Create XDMF files for visualization output
xdmffile_u = XDMFFile('navier_stokes_cylinder/velocity.xdmf')
xdmffile_p = XDMFFile('navier_stokes_cylinder/pressure.xdmf')

# Create time series (for use in reaction_system.py)
timeseries_u = TimeSeries('navier_stokes_cylinder/velocity_series')
timeseries_p = TimeSeries('navier_stokes_cylinder/pressure_series')

# Save mesh to file (for use in reaction_system.py)
File('navier_stokes_cylinder/cylinder.xml.gz') << mesh

# Create progress bar
# progress = Progress('Time-stepping')
progress = Progress('Time-stepping', num_steps)
# set_log_level(PROGRESS)

# Time-stepping
t = 0
for n in range(num_steps):

    # Update current time
    t += dt

    # Step 1: Tentative velocity step
    b1 = assemble(L1)
    [bc.apply(b1) for bc in bcu]
    solve(A1, u_.vector(), b1, 'bicgstab', 'hypre_amg')

    # Step 2: Pressure correction step
    b2 = assemble(L2)
    [bc.apply(b2) for bc in bcp]
    solve(A2, p_.vector(), b2, 'bicgstab', 'hypre_amg')

    # Step 3: Velocity correction step
    b3 = assemble(L3)
    solve(A3, u_.vector(), b3, 'cg', 'sor')

    # Plot solution
    # plot(u_, title='Velocity')
    # plot(p_, title='Pressure')

    # Save solution to file (XDMF/HDF5)
    xdmffile_u.write(u_, t)
    xdmffile_p.write(p_, t)

    # Save nodal values to file
    timeseries_u.store(u_.vector(), t)
    timeseries_p.store(p_.vector(), t)

    # Update previous solution
    u_n.assign(u_)
    p_n.assign(p_)

    # Update progress bar
    # progress.update(t / T)
    set_log_level(LogLevel.PROGRESS)
    progress += 1
    set_log_level(LogLevel.ERROR)
    print('u max:', u_.vector().get_local().max())

[dokken]

Start by inspecting your mesh.
As you can see by visualizing the problem in Paraview, your mesh is not as you would expect:

For instance have a look at how to create your mesh with gmsh (You can find a version of this mesh created with gmsh in my github repo).
To convert this mesh to XDMF, and read it into dolfin, consider: Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio or Converter from GMSH to XDMF (with physical groups)

[Sebastian]

OK thank you very much
I was able to create the new mesh as seen in the image but now I have a problem reading the “facets.xdmf” file and defining the boundary conditions, i use this code for create the mesh https://github.com/jorgensd/dolfinx_ipcs/blob/a356280a0aa49d439187a52ce13b1a4bdef0eda2/create_and_convert_3D_mesh.py. when i run the code i get this error:

Traceback (most recent call last):
File “naviercilindro3d-2-1.py”, line 69, in
bcu_walls = DirichletBC(V, Constant((0, 0, 0)), boundaries)
File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/dirichletbc.py”, line 131

This is my code:

from __future__ import print_function
from fenics import *
from mshr import *
import numpy as np
    
    
T = 5.0            # final time
num_steps = 5000   # number of time steps
dt = T / num_steps # time step size
mu = 0.001         # dynamic viscosity
rho = 1            # density


# Read mesh 
mesh = Mesh()

with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)

with XDMFFile("facets.xdmf") as infile:
    infile.read(mesh)

tdim = mesh.topology().dim()

# Create facet function over mesh, unitialized values
boundaries = MeshFunction('size_t', mesh, tdim-1)

# Set all values to 0
boundaries.set_all(0)

# Actually use facet labels from gmsh
boundaries = MeshFunction('size_t', mesh, mesh.domains())


# Define function spaces
V = VectorFunctionSpace(mesh, 'P', 2)
Q = FunctionSpace(mesh, 'P', 1)

# Define boundaries
inflow   = 'near(x[0], 0)'
outflow  = 'near(x[0], 2.2)'
"""walls    = 'near(x[1], 0) || near(x[1], 0.41) || near(x[2], 0) || near(x[2], 0.41)'
cylinder = 'on_boundary && x[0]>0.1 && x[0]<0.3 && x[1]>0.1 && x[1]<0.3'"""

# Define inflow profile
inflow_profile = ('16.0*1.5*x[1]*x[2]*(0.41 - x[1])*(0.41 - x[2]) / pow(0.41, 4)', '0', '0')


# Define boundary conditions
bcu_inflow = DirichletBC(V, Expression(inflow_profile, degree=3), inflow)
bcu_walls = DirichletBC(V, Constant((0, 0, 0)), boundaries)
bcu_cylinder = DirichletBC(V, Constant((0, 0, 0)), boundaries)
bcp_outflow = DirichletBC(Q, Constant(0), outflow)
bcu = [bcu_inflow, bcu_walls, bcu_cylinder]
bcp = [bcp_outflow]

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

# Define functions for solutions at previous and current time steps
u_n = Function(V)
u_  = Function(V)
p_n = Function(Q)
p_  = Function(Q)

# Define expressions used in variational forms
U  = 0.5*(u_n + u)
n  = FacetNormal(mesh)
f  = Constant((0, 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
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
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
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)

# Apply boundary conditions to matrices
[bc.apply(A1) for bc in bcu]
[bc.apply(A2) for bc in bcp]

# Create XDMF files for visualization output
xdmffile_u = XDMFFile('navier_stokes_cylinder/velocity.xdmf')
xdmffile_p = XDMFFile('navier_stokes_cylinder/pressure.xdmf')

# Create time series (for use in reaction_system.py)
timeseries_u = TimeSeries('navier_stokes_cylinder/velocity_series')
timeseries_p = TimeSeries('navier_stokes_cylinder/pressure_series')

# Save mesh to file (for use in reaction_system.py)
File('navier_stokes_cylinder/cylinder.xml.gz') << mesh

# Create progress bar
# progress = Progress('Time-stepping')
progress = Progress('Time-stepping', num_steps)
# set_log_level(PROGRESS)

# Time-stepping
t = 0
for n in range(num_steps):

    # Update current time
    t += dt

    # Step 1: Tentative velocity step
    b1 = assemble(L1)
    [bc.apply(b1) for bc in bcu]
    solve(A1, u_.vector(), b1, 'bicgstab', 'hypre_amg')

    # Step 2: Pressure correction step
    b2 = assemble(L2)
    [bc.apply(b2) for bc in bcp]
    solve(A2, p_.vector(), b2, 'bicgstab', 'hypre_amg')

    # Step 3: Velocity correction step
    b3 = assemble(L3)
    solve(A3, u_.vector(), b3, 'cg', 'sor')

    # Plot solution
    # plot(u_, title='Velocity')
    # plot(p_, title='Pressure')

    # Save solution to file (XDMF/HDF5)
    xdmffile_u.write(u_, t)
    xdmffile_p.write(p_, t)

    # Save nodal values to file
    timeseries_u.store(u_.vector(), t)
    timeseries_p.store(p_.vector(), t)

    # Update previous solution
    u_n.assign(u_)
    p_n.assign(p_)

    # Update progress bar
    # progress.update(t / T)
    set_log_level(LogLevel.PROGRESS)
    progress += 1
    set_log_level(LogLevel.ERROR)
    print('u max:', u_.vector().get_local().max())

[dokken]

Please add the full error message.

[Sebastian]

raceback (most recent call last):
File “naviercilindro3d-2-1.py”, line 69, in
bcu_walls = DirichletBC(V, Constant((0, 0, 0)), boundaries)
File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/dirichletbc.py”, line 131, in init
super().init(*args)
TypeError: init(): incompatible constructor arguments. The following argument types are supported:
1. dolfin.cpp.fem.DirichletBC(arg0: dolfin.cpp.fem.DirichletBC)
2. dolfin.cpp.fem.DirichletBC(V: dolfin.cpp.function.FunctionSpace, g: dolfin.cpp.function.GenericFunction, sub_domain: dolfin.cpp.mesh.SubDomain, method: str = ‘topological’, check_midpoint: bool = True)
3. dolfin.cpp.fem.DirichletBC(V: dolfin.cpp.function.FunctionSpace, g: dolfin.cpp.function.GenericFunction, sub_domains: dolfin.cpp.mesh.MeshFunctionSizet, sub_domain: int, method: str = ‘topological’)

Invoked with: <dolfin.cpp.function.FunctionSpace object at 0x7f78fb5449a0>, <dolfin.cpp.function.Constant object at 0x7f78fbb30330>, <dolfin.cpp.mesh.MeshFunctionSizet object at 0x7f78fbab8170>, ‘topological’

[dokken]

You Need to supply Which tag on your boundary you are using from the meshfunction

[volkerk]

An alternative to the mesh generation with gmsh might be:

from fenics import *
from mshr import *
from vedo.dolfin import plot

# Create mesh
channel = Rectangle(Point(0, 0), Point(2.2, 0.41))
cylinder = Circle(Point(0.2, 0.2), 0.05, 32)
domain = Extrude2D(channel - cylinder, 0.41)

mesh = generate_mesh(domain, 64, 'cgal')
plot(mesh)

mesh = generate_mesh(domain, 64, 'tetgen')
plot(mesh)