Stokes flow around a sphere¶
Formulation¶
In this example, we numerically reproduce the famous Stokes law (see Landau, Lifschitz: Fluid mechanics, Pergamon Press, Oxford, 1966)
which evaluates the drag on a sphere of radius \(R\) moving at speed \(v\). We load an external mesh which was prepared using Gmsh software.
Instead of moving the sphere in a time-dependant geometry, it is much simpler to fix the sphere and let the fluid move around it. Also, since FEM works on bounded domains, we confine the flow in a finite cylinder of radius \(R_\mathrm{cylinder} = 100R\) and length \(L_\mathrm{cylinder} = 100R\). At the walls and the inflow, we define a Dirichlet boundary condition
homogenous Dirichlet condition on sphere and natural boundary condition at the outflow. For \(R = 1\) and \(\mu = 1\), the drag predicted by theory is \(6\pi\).
This mesh has coarse resolution of the cylinder but is significantly refined near the sphere. Ability to work with non-uniform resolution is a significant advantage of FEM.
Implementation¶
We begin by importing few things. Since we will run this on multiple cores using MPI, we need to access the rank of threads.
[12]:
import dolfin as df
from time import time
from numpy import pi
comm = df.MPI.comm_world
rank = df.MPI.rank(comm)
Load the mesh from a file called mesh_big.h5.
[13]:
mesh = df.Mesh()
hdf = df.HDF5File(mesh.mpi_comm(), "mesh_big.h5", "r")
hdf.read(mesh, "/mesh", False)
Read boundary marker from the file.
[14]:
bdary = df.MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
hdf.read(bdary, "/boundaries")
#1 = INFLOW
#2 = WALLS
#3 = OUTFLOW
#4 = SPHERE
Define function spaces. For extra speed, let us use Crouzeix-Raviart elements.
[16]:
Ev = df.VectorElement("CR", mesh.ufl_cell(), 1)
Ep = df.FiniteElement("DG", mesh.ufl_cell(), 0)
W = df.FunctionSpace(mesh, df.MixedElement([Ev, Ep]))
Then, we define boundary conditions as stated above.
[17]:
V0 = df.Constant((0.0, 0.0, 0.0))
Vin = df.Constant((1.0, 0.0, 0.0))
bc_inflow = df.DirichletBC(W.sub(0), Vin, bdary, 1)
bc_walls = df.DirichletBC(W.sub(0), Vin, bdary, 2)
bc_sphere = df.DirichletBC(W.sub(0), V0, bdary, 4)
bc = [bc_inflow, bc_walls, bc_sphere]
The variational form and solver will be same as in two dimensions. When printing output, we use if rank == 0 condition. (Otherwise, it would be printed by every thread.)
[18]:
nu = df.Constant(1.0)
u, p = df.TrialFunctions(W)
v, q = df.TestFunctions(W)
F = nu*df.inner(df.grad(u), df.grad(v))*df.dx - p*df.div(v)*df.dx - q*df.div(u)*df.dx
w = df.Function(W)
problem = df.LinearVariationalProblem(df.lhs(F), df.rhs(F), w, bc)
solver = df.LinearVariationalSolver(problem)
solver.parameters["linear_solver"] = "mumps"
Solve the linear system:
[19]:
tick = time()
solver.solve()
if rank == 0:
print("ellapsed = ", time() - tick, "s")
ellapsed = Solving linear variational problem.
111.6617522239685 s
Save the result for inspection in Paraview.
[20]:
ufile = df.XDMFFile("results/u.xdmf")
pfile = df.XDMFFile("results/p.xdmf")
u, p = w.split()
u.rename("u", "velocity")
p.rename("p", "pressure")
ufile.write(u)
pfile.write(p)
ufile.close()
pfile.close()
Report the drag. Compare it to the Stokes’ law.
[24]:
n = df.FacetNormal(mesh)
I = df.Identity(mesh.geometry().dim())
ds = df.Measure("ds", subdomain_data=bdary)
T = -p*I + 2.0*nu*df.sym(df.grad(u))
force = df.dot(T, n)
drag = df.assemble(-force[0]*ds(4))
#alt force computation
bc_ = df.DirichletBC(W.sub(0), (1.0, 0.0, 0.0), bdary, 4)
w_ = df.Function(W)
bc_.apply(w_.vector())
drag_ = -df.assemble(df.inner(T, df.grad(w_.sub(0)))*df.dx)
if rank == 0:
print("drag by surface integral = {:.4}".format(drag))
print("drag by volume integral = {:.4}".format(drag_))
print("drag according to Stokes' law = {:.4}".format(6*pi))
print("relative error = {:.4%}".format(drag/(6*pi) - 1.0))
drag by surface integral = 18.92
drag by volume integral = 18.92
drag according to Stokes' law = 18.85
relative error = 0.3954%
Remark¶
In mesh_big.h5, the ratio of cylinder radius to sphere radius is 100:1. This is not an overkill. You can try to recompute this problem with mesh_small.h5, where the ratio is 10:1. You will find that the boundary effect is strong and leads to much higher drag.
Complete Code¶
Run in parallel using the command
$ mpirun -n 4 python3 Stokes3D.py
where 4 is number of threads.
[ ]:
import dolfin as df
from time import time
from numpy import pi
comm = df.MPI.comm_world
rank = df.MPI.rank(comm)
mesh = df.Mesh()
hdf = df.HDF5File(mesh.mpi_comm(), "mesh_big.h5", "r")
hdf.read(mesh, "/mesh", False)
bdary = df.MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
hdf.read(bdary, "/boundaries")
#1 = INFLOW
#2 = WALLS
#3 = OUTFLOW
#4 = SPHERE
Ev = df.VectorElement("CR", mesh.ufl_cell(), 1)
Ep = df.FiniteElement("DG", mesh.ufl_cell(), 0)
W = df.FunctionSpace(mesh, df.MixedElement([Ev, Ep]))
V0 = df.Constant((0.0, 0.0, 0.0))
Vin = df.Constant((1.0, 0.0, 0.0))
bc_inflow = df.DirichletBC(W.sub(0), Vin, bdary, 1)
bc_walls = df.DirichletBC(W.sub(0), Vin, bdary, 2)
bc_sphere = df.DirichletBC(W.sub(0), V0, bdary, 4)
bc = [bc_inflow, bc_walls, bc_sphere]
nu = df.Constant(1.0)
u, p = df.TrialFunctions(W)
v, q = df.TestFunctions(W)
F = nu*df.inner(df.grad(u), df.grad(v))*df.dx - p*df.div(v)*df.dx - q*df.div(u)*df.dx
w = df.Function(W)
problem = df.LinearVariationalProblem(df.lhs(F), df.rhs(F), w, bc)
solver = df.LinearVariationalSolver(problem)
solver.parameters["linear_solver"] = "mumps"
tick = time()
solver.solve()
if rank == 0:
print("ellapsed = ", time() - tick, "s")
ufile = df.XDMFFile("results/u.xdmf")
pfile = df.XDMFFile("results/p.xdmf")
u, p = w.split()
u.rename("u", "velocity")
p.rename("p", "pressure")
ufile.write(u)
pfile.write(p)
n = df.FacetNormal(mesh)
I = df.Identity(mesh.geometry().dim())
ds = df.Measure("ds", subdomain_data=bdary)
T = -p*I + 2.0*nu*df.sym(df.grad(u))
force = df.dot(T, n)
drag = df.assemble(-force[0]*ds(4))
if rank == 0:
print("drag according to simulation = {:.4}".format(drag))
print("drag according to Stokes' law = {:.4}".format(6*pi))
print("relative error = {:.4%}".format(drag/(6*pi) - 1.0))