Hello,
this Fenics sample code for the Navier-Stokes equations with boundary conditions on the pressure runs fine with the default parameters set on the webpage, i.e., L=h=1, \mu = 1, where L and h are the width and the height of the rectangle, and \mu the viscosity.
However, if I let the exact same code run with the same parameters as this other Navier Stokes sampe program, i.e., L= 2.2, h=0.41, \mu = 0.001, it blows up (the velocity has gigantic values).
Please see this minimal working example:
#generate_mesh.py
import numpy
import meshio
import gmsh
import pygmsh
import argparse
parser = argparse.ArgumentParser()
parser.add_argument("resolution")
args = parser.parse_args()
#mesh resolution
resolution = (float)(args.resolution)
# Channel parameters
L = 2.2
h = 0.41
# Initialize empty geometry using the build in kernel in GMSH
geometry = pygmsh.geo.Geometry()
model = geometry.__enter__()
#add the points that define the boundary of the mesh
points = [model.add_point((0, 0, 0), mesh_size=resolution),
model.add_point((L, 0, 0), mesh_size=resolution),
model.add_point((L, h, 0), mesh_size=resolution),
model.add_point((0, h, 0), mesh_size=resolution)
]
channel_lines = [model.add_line(points[0], points[1]),
model.add_line(points[1], points[2]),
model.add_line(points[2], points[3]),
model.add_line(points[3], points[0])]
channel_loop = model.add_curve_loop(channel_lines)
plane_surface = model.add_plane_surface(channel_loop, holes=[])
model.synchronize()
model.add_physical([plane_surface], "Volume")
model.add_physical([channel_lines[3]], "Inflow")
model.add_physical([channel_lines[1]], "Outflow")
geometry.generate_mesh(64)
gmsh.write("membrane_mesh.msh")
gmsh.clear()
geometry.__exit__()
mesh_from_file = meshio.read("membrane_mesh.msh")
def create_mesh(mesh, cell_type, prune_z=False):
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
points = mesh.points[:, :2] if prune_z else mesh.points
out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={
"name_to_read": [cell_data]})
return out_mesh
line_mesh = create_mesh(mesh_from_file, "line", prune_z=True)
meshio.write("line_mesh.xdmf", line_mesh)
triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=True)
meshio.write("triangle_mesh.xdmf", triangle_mesh)
#ft07_navier_stokes_channel.py
"""
FEniCS tutorial demo program: Incompressible Navier-Stokes equations
for channel flow (Poisseuille) on the unit square using the
Incremental Pressure Correction Scheme (IPCS).
u' + u . nabla(u)) - div(sigma(u, p)) = f
div(u) = 0
"""
from __future__ import print_function
from fenics import *
import numpy as np
import argparse
parser = argparse.ArgumentParser()
parser.add_argument("input_directory")
parser.add_argument("output_directory")
parser.add_argument("T")
parser.add_argument("N")
args = parser.parse_args()
T = (float)(args.T)
num_steps = (int)(args.N)
dt = T / num_steps # time step size
mu = 0.001 # kinematic viscosity
rho = 1 # density
# Create XDMF files for visualization output
xdmffile_v = XDMFFile((args.output_directory) + "/v.xdmf")
xdmf_file_p = XDMFFile((args.output_directory) + "/p.xdmf")
# Create mesh and define function spaces
# mesh = UnitSquareMesh(16, 16)
#create mesh
mesh=Mesh()
with XDMFFile((args.input_directory) + "/triangle_mesh.xdmf") as infile:
infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 2)
with XDMFFile((args.input_directory) + "/line_mesh.xdmf") as infile:
infile.read(mvc, "name_to_read")
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)'
# Define boundary conditions
bcu_noslip = DirichletBC(V, Constant((0, 0)), walls)
bcp_inflow = DirichletBC(Q, Constant(8), inflow)
bcp_outflow = DirichletBC(Q, Constant(0), outflow)
bcu = [bcu_noslip]
bcp = [bcp_inflow, 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))
k = Constant(dt)
mu = Constant(mu)
rho = Constant(rho)
# Define strain-rate tensor
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]
# Time-stepping
t = 0
for n in range(num_steps):
# Write the solution to file
xdmffile_v.write(u_n, t)
xdmf_file_p.write(p_n, t)
# 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)
# Step 2: Pressure correction step
b2 = assemble(L2)
[bc.apply(b2) for bc in bcp]
solve(A2, p_.vector(), b2)
# Step 3: Velocity correction step
b3 = assemble(L3)
solve(A3, u_.vector(), b3)
# Update previous solution
u_n.assign(u_)
p_n.assign(p_)
print("\t%.2f %%" % (100.0 * (t / T)), flush=True)
I run with $python3 generate_mesh.py 0.1 to generate the mesh, and run the simulation with $python3 ft07_navier_stokes_channel.py [path where the mesh files were stored by generate_mesh.py] [path where to store the solution] 2 2048. Here is the result for the magnitude of the velocity vector at some intermediate time (notice the scale in the color bar):
Why this blow up? I wonder whether this is related to this comment by @dokken in this post
Thank you 