Dear Community
I am trying to solve the Stokes flow problem in a 2D axisymmetric geometry, which uses spherical polar coordinates, which is different from the flow past a cylinder in 2D. The 2D axisymmetric spherical one has a well-known analytical solution, and I thought this would enable validation of my FENICS result. I am about to describe my approach below, and present the FENICS code below that, followed by why I think the numerically computed solution, and would appreciate any input/guidance on this matter.
If you think that the method I have adopted is painful and know of simpler alternatives, I am open to those suggestions as well.
1. Variational formulation
I leveraged my experience in solving the heat-conduction problem in a 2D spherical axisymmetric geometry and set out to similarly derive the variational formulation of the Stokes equation in polar coordinates. The challenge here is that the velocity field is a vector, as compared to the heat conduction case, where we solve for a scalar field. The domain/geometry of the problem is clear from figure (A) below.
(A) Geometry
That is, I am essentially solving
\langle\nabla^2\mathbf{u},\mathbf{v}\rangle_{\Omega} - \langle\nabla p,\mathbf{v}\rangle_{\Omega}=0
\langle\left(\nabla\cdot\mathbf{u}\right),q\rangle_{\Omega}=0
with \mathbf{u} known at the inner and outer boundaries [Dirichlet condition], and
p\mathbf{n}=\left(\nabla\mathbf{u}\right)\cdot\mathbf{n}-\mathbf{g}
as the “do-nothing” boundary condition, and I set \mathbf{g}=0 in the code.
I found it necessary to specify the divergence and gradient operators in polar coordinates, so as to enforce 2D spherical axisymmetry. For instance,
\langle\nabla^2\mathbf{u},\mathbf{v}\rangle_{\Omega}=\int_{r}\int_{\theta}\left[\nabla^2\mathbf{u}\cdot\mathbf{v}\right]rdrd\theta
The variational formulation in polar spherical coordinates is obtained as
and the equivalent Cartesian representation, which is used in the code below, is
Obtaining these equations was a bit of a grind, and the steps are fairly along the lines outlined here. I have checked the math, but it is possible that I may have made a mistake.
2. Code snippet
I assembled the following FEniCS script for solving the problem. I found the tutorial by G Linga and A Bolet to be very useful in this regard.
## Solving Stokes flow in spherical annulus
from dolfin import *
import numpy as np
from mshr import *
Re = 10.
Ri = 1.
num_segments=400
domain = Circle(Point(0., 0.), Re, num_segments) - Circle(Point(0., 0.), Ri, num_segments)
#creating semi-circular domain
domain = (domain - Rectangle(Point(0., 0.), Point(Re, -Re))
- Rectangle(Point(0., 0.), Point(-Re, -Re)))
mesh_resolution = 128
mesh = generate_mesh(domain, mesh_resolution)
x = SpatialCoordinate(mesh)
U_0 = 1.0 # far-field velocity
# Defining conversion to polar coordinates
r = Expression("sqrt(x[0]*x[0]+x[1]*x[1])", degree=2)
theta = Expression("atan2(x[1],x[0])", degree=2)
# Making a mark dictionary
# Note: the values should be UNIQUE identifiers.
mark = {"generic": 0,
"inner": 1,
"outer": 2}
subdomains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
subdomains.set_all(mark["generic"])
class inner_boun(SubDomain):
def inside(self, x, on_boundary):
tol = 1.5*Ri * (1 - np.cos(np.pi / num_segments))
return near(sqrt(x[0]**2 + x[1]**2), Ri, tol) and on_boundary
class outer_boun(SubDomain):
def inside(self, x, on_boundary):
tol = 1.5*Re * (1 - np.cos(np.pi / num_segments))
return near(sqrt(x[0]**2 + x[1]**2), Re, tol) and on_boundary
## Marking boundaries
in_boun = inner_boun()
in_boun.mark(subdomains, mark["inner"])
out_boun = outer_boun()
out_boun.mark(subdomains, mark["outer"])
geom_name="geometry.pvd"
File(geom_name) << subdomains
## Define function spaces
V = VectorElement("CG", triangle, 2)
P = FiniteElement("CG", triangle, 1)
W = FunctionSpace(mesh, V*P)
#
## Define variational problem
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
#
dx = Measure("dx", domain=mesh, subdomain_data=subdomains) # Volume integration
ds = Measure("ds", domain=mesh, subdomain_data=subdomains) # Surface integration
## Applying Dirichlet BCs at the inner and outer boundaries
# No-slip at inner boundary
# Constant velocity at outer boundary
u_outer_boun = Constant((U_0, 0.0))
bc_outer_boun = DirichletBC(W.sub(0), u_outer_boun, subdomains, mark["outer"])
noslip = Constant((0.0, 0.0))
bc_inner_boun = DirichletBC(W.sub(0), noslip, subdomains, mark["inner"])
# fixing the pressure at one node to some reference value
# in order to have a well-posed problem
bc_p = DirichletBC(W.sub(1), 0, "x[0] < -1.0 + DOLFIN_EPS", "pointwise")
bcs = [bc_outer_boun,bc_inner_boun,bc_p]
# Body force:
force = Constant((0.0, 0.0))
##Setting up variational formulation
# Separately extracting the x and y components of u and v
vx, vy = v[0], v[1]
ux, uy = u[0], u[1]
a1 = ((ux.dx(0)*vx.dx(0)) + (uy.dx(0)*vy.dx(0))
+ (ux.dx(1)*vx.dx(1)) + (uy.dx(1)*vy.dx(1)))
a2 = (1./r)*cos(theta)*((ux*vy.dx(1)-uy*vx.dx(1)) + (ux*vx.dx(0)+uy*vy.dx(0))
+ (vx*ux.dx(0)+vy*uy.dx(0)) - (vx*uy.dx(1)-vy*ux.dx(1)))
a3 = (1./r)*sin(theta)*(-(ux*uy.dx(0)-uy*vy.dx(0)) + (ux*vx.dx(1)+uy*vy.dx(1))
+(vx*uy.dx(0)-vy*ux.dx(0)) +(vx*ux.dx(1)+vy*uy.dx(1)))
a4 = (inner(u,v)/(r*r*sin(theta)*sin(theta))
+(x[0]/(x[1]*r*r)) * (x[1]*(vx*ux.dx(0)-vy*uy.dx(0))
-x[0]*(vx*ux.dx(1)-vy*uy.dx(1))))
a_uv = a1+a2+a3+a4
a5 = (vx.dx(0) + vy.dx(1) + cos(theta)*sin(theta)*(vx.dx(0)-vy.dx(0)+(vx/x[1])
+(x[0]/x[1])*(vy/x[1]) + (2*vy/x[0])))*p
a6 = (ux.dx(0) + uy.dx(1) + cos(theta)*sin(theta)*(ux.dx(0)-uy.dx(0)+(ux/x[1])
+(x[0]/x[1])*(uy/x[1]) + (2*uy/x[0])))*q
F = (a_uv-a5-a6)*dx - inner(force, v) * dx
a = lhs(F)
L = rhs(F)
# This method avoids the "too many values to unpack" error
### Solve variational problem
w = Function(W)
solve(a == L, w, bcs)
### Split using deep copy
(u, p) = w.split(True)
## Save solution in VTK format
ufile = File("velocity.pvd")
ufile << u
pfile = File("pressure.pvd")
pfile << p
3. Results and doubts
Given below are the pressure and velocity fields from executing the FENICS script given above. The reason I suspect that something is wrong is because the magnitude of the velocity in certain regions is seen to exceed the far-field value of U_0 (which is set as 1 in this problem). Any help is highly appreciated.
(B) Pressure
(C) Magnitude of velocity
Thank You
Warm Regards