Dear Community
I am trying to solve the Stokes flow past a sphere in axisymmetric spherical coordinates using the streamfunction (\psi)-vorticity (w) framework. The geometry is an annular domain as shown below.
As a way of testing my problem setup, I am adopting the method of manufactured solutions, by using the analytically known values of the streamfunction and vorticity as Dirichlet boundary conditions at the surface of the inner and outer spheres.
1. Doubt
I considered two values of the outer sphere radius, 5 and 40, with the radius of the inner sphere fixed at unity. At both the values, the maximum difference between the computed and exact vorticity is of the order of 10^{-5} to 10^{-4} [error higher for the larger radius case], which I deem to acceptable. However, the difference between the computed and exact stream function are orders of magnitude greater than in the case of vorticity, and is also higher for the case of larger radius, as shown in the plots below.
I am not sure how to interpret this behavior, and how to troubleshoot the issue. I would like to think that the coupled system of partial differential equations have been written in the variational form properly and are being solved correctly, but the trend shown in the pictures worries me. The mathematical setup and the relevant code snippet is given below.
I would really appreciate any help/guidance on the matter.
Thank You
Warm Regards
2. Mathematical setup
The coupled governing equations are
w r\sin\theta+E^2\psi=0\\
\nabla^2 w-\dfrac{w}{r^2\sin^2\theta}=0
where
\nabla^2\equiv\dfrac{1}{r^2}\dfrac{\partial}{\partial r}\left(r^2\dfrac{\partial}{\partial r}\right)+\dfrac{1}{r^2\sin\theta}\dfrac{\partial}{\partial \theta}\left(\sin\theta\dfrac{\partial}{\partial \theta}\right)\\[10pt]
E^2\equiv\dfrac{\partial^2}{\partial r^2} + \dfrac{\sin\theta}{r^2}\dfrac{\partial}{\partial \theta}\left(\dfrac{1}{\sin\theta}\dfrac{\partial}{\partial \theta}\right)
and the following exact solutions are applied as boundary conditions on r=1 and r=R_o
\psi_{\text{exact}}(r,\theta)=U\sin^2\theta\left[\dfrac{r^2}{2}-\dfrac{3r}{4}+\dfrac{1}{4r}\right]
w_{\text{exact}}(r,\theta)=-\dfrac{3U\sin\theta}{2r^2}
The variational form is as follows
\int\int\Biggl[\left(\dfrac{\partial w}{\partial x}\dfrac{\partial v}{\partial x}+\dfrac{\partial w}{\partial y}\dfrac{\partial v}{\partial y}\right)-\dfrac{v}{r\sin\theta}\dfrac{\partial w}{\partial y}+\dfrac{wv}{r^2\sin^2\theta}\Biggr]dxdy=0\\[5pt]
\int\int\Biggl[\left(\dfrac{\partial \psi}{\partial x}\dfrac{\partial u}{\partial x}+\dfrac{\partial \psi}{\partial y}\dfrac{\partial u}{\partial y}\right)+\dfrac{u}{r\sin\theta}\dfrac{\partial \psi}{\partial y}-wr\sin\theta u\Biggr]dxdy=0
where w and \psi are the trial functions and v and u are the respective test functions defined on the mixed function space W=\Omega\times\Psi such that (w,v)\in\Omega and (\psi,u)\in\Psi.
The code segment for obtaining the results discussed in this post is given below
3. FEniCS script
from dolfin import *
from mshr import *
import math
import numpy as np
Re = 40. #parameter critical to discussion above
Ri = 1.
rect = Rectangle(Point(0., 0.), Point(Re, Re))
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 = 512
mesh = generate_mesh(domain, mesh_resolution)
x = SpatialCoordinate(mesh)
strm = FiniteElement("CG",triangle,2)
vort = FiniteElement("CG",triangle,2)
W = FunctionSpace(mesh,strm*vort)
r = Expression("sqrt(x[0]*x[0]+x[1]*x[1])", degree=2)
theta = Expression("atan2(x[1],x[0])", degree=2)
U = Constant(1.0)
## Constructing known solution for \psi and vorticity and applying on boundaries
psi_mfg = Expression("U*sin(theta)*sin(theta)*(0.5*r*r-0.75*r+(0.25/r))",r=r,
theta=theta,U=U,degree=1)
omeg_mfg = Expression("-1.5*U*sin(theta)/(r*r)",r=r,
theta=theta,U=U,degree=1)
### Defining boundary conditions
### and marking subdomains
class Inside(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 Outside(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
sub_domains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
sub_domains.set_all(0)
Inside().mark(sub_domains, 1)
Outside().mark(sub_domains, 2)
## Assembling the BCs
inner_circ_psi = DirichletBC(W.sub(0), psi_mfg, sub_domains, 1)
outer_circ_psi = DirichletBC(W.sub(0), psi_mfg, sub_domains, 2)
inner_circ_omeg = DirichletBC(W.sub(1), omeg_mfg, sub_domains, 1)
outer_circ_omeg = DirichletBC(W.sub(1), omeg_mfg, sub_domains, 2)
bcs = [inner_circ_psi, outer_circ_psi, inner_circ_omeg, outer_circ_omeg]
### Defining test and trial functions
(psi, omeg) = TrialFunctions(W)
(u, v) = TestFunctions(W)
## Assembling bilinear equation
a1 = (inner(grad(omeg),grad(v))-(Constant(1.0)/(r*sin(theta)))*v*(omeg.dx(1))
+(omeg*v/(r*r*sin(theta)*sin(theta))))
a2 = (inner(grad(psi),grad(u))+(Constant(1.0)/(r*sin(theta)))*u*(psi.dx(1))
-(omeg*r*sin(theta)*u))
a = (a1+a2)
f = Constant(0.0)
g = Constant(0.0)
a = (a1+a2)*dx
L = f*v*dx+g*u*dx
## Compute solution
w = Function(W)
solve(a == L, w, bcs)
# split using deep copy
(psi, omeg) = w.split(True)
# Error calculation
error_L2_psi = errornorm(psi_mfg,psi, 'L2')
print ('L2 Error in psi : ', error_L2_psi)
error_L2_omeg = errornorm(omeg_mfg,omeg, 'L2')
print ('L2 Error in omega : ', error_L2_omeg)
# compute maximum error at vertices
vertex_values_psi_e = psi_mfg.compute_vertex_values(mesh)
vertex_values_psi = psi.compute_vertex_values(mesh)
error_max_psi = np.max(np.abs(vertex_values_psi_e - vertex_values_psi))
print('error_max in psi : ',error_max_psi)
vertex_values_omeg_e = omeg_mfg.compute_vertex_values(mesh)
vertex_values_omeg = omeg.compute_vertex_values(mesh)
error_max_omeg = np.max(np.abs(vertex_values_omeg_e - vertex_values_omeg))
print('error_max in omega : ',error_max_omeg)
#writing output to file
psi_file = XDMFFile("psi.xdmf")
psi_file.write(psi)
omeg_file = XDMFFile("omega.xdmf")
omeg_file.write(omeg)
exact_psi = project(psi_mfg,W.sub(0).collapse())
exact_omeg = project(omeg_mfg,W.sub(1).collapse())
diff_psi = project((psi_mfg-psi),W.sub(0).collapse())
diff_omeg = project((omeg_mfg-omeg),W.sub(1).collapse())
ex_psi_file = XDMFFile("ex_psi.xdmf")
ex_psi_file.write(exact_psi)
ex_omeg_file = XDMFFile("ex_omeg.xdmf")
ex_omeg_file.write(exact_omeg)
diff_psi_file = XDMFFile("diff_psi.xdmf")
diff_psi_file.write(diff_psi)
diff_omeg_file = XDMFFile("diff_omeg.xdmf")
diff_omeg_file.write(diff_omeg)