Stokes equations in 2D spherical axisymmetric coordinates

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