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
with \mathbf{u} known at the inner and outer boundaries [Dirichlet condition], and

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,


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.
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)

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"])

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))

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