Hi,
I am trying to solve the steady state Navier Stokes equations in cylindrical co-ordinates
for the following boundary layer problem where the flow is assumed to be axisymmetric.
\begin{align}
\frac{1}{ r} u + \frac{\partial u }{\partial r} + \frac{\partial w }{\partial Z} = 0 , \\
-\frac{\partial p }{\partial r} + \frac{1}{ r} \frac{\partial u }{\partial r} + \frac{\partial^2 u }{\partial r^2} + \frac{\partial^2 u }{\partial Z^2}- \frac{ u }{ r^2} = 0 , \\
\frac{\partial p }{\partial Z} + \frac{1}{ r} \frac{\partial w }{\partial r} + \frac{\partial^2 w }{\partial r^2} + \frac{\partial^2 w }{\partial Z^2} = 0, \\
w = 0 ,\quad u = 0, \quad R = 0 \mbox{ at } Z = 0, \\
p =- 1, \quad u = -\frac{ r}{2} \mbox{ at } Z = \infty , \\
| u |, \; | w |, \; | p | \mbox{ finite and symmetric at } r = 0, \\
u = \frac{\mathrm{d} R}{\mathrm{d} Z } \mbox{ at } r = 1 , \\
- p + 2\frac{\partial u }{\partial r} = 0 \mbox{ at } r = 1 ,\\
\frac{\partial u }{\partial Z} + \frac{\partial w }{\partial r} = 0 \mbox{ at } r = 1 .
\end{align}
However I get a pressure singularity at (r = 0, z =0).
The finite and symmetric condition does not get directly imposed but all
boundary terms at r = 0 in the integration by parts dissapear from the co-ordinate system transformation ie. dx dy \to r dr dz.
I additionally imposed that the radial velocity is zero at r = 0, which removed singularities for the radial velocity and reduced the magnitude of the spikiness at the origin for the pressure.
I also tried moving my domain from r = 0 \to r = 0.001 but it didn’t resolve the issue.
Any ideas how I might avoid this non physical pressure?
I’m using Taylor-Hood elements
Relevant Code for reference
from dolfin import *
import numpy as np
L = 10
Z_m = 80*L
R_m = 80
tol = 0.00
# Define mesh and geometry -
mesh = RectangleMesh(Point(tol, 0), Point(1, L), R_m, Z_m )
n = FacetNormal(mesh)
# Define Taylor--Hood function space W
V = VectorElement("Lagrange", triangle, 2)
Q = FiniteElement("Lagrange", triangle, 1)
S = FiniteElement("Lagrange", triangle, 2)
W = FunctionSpace(mesh, MixedElement([V, Q, S]))
# Define Function and TestFunction(s)
w = Function(W)
# speed convergence
w = interpolate(Expression(('x[1]*0','x[1]','-1*exp(0*x[1])','-0.5*x[1]'), degree=2) ,W)
(u, p, r) = split(w)
(v, q, s) = split(TestFunction(W))
u_out = Expression(('-x[0]*0.5'), degree=2)
p_out = Expression(('-1'), degree=1)
u_in = Expression(('0', '(1-x[0]*x[0])'), degree=2)
out = 'near(x[1],' + str(L) + ')'
FB = 'near(x[0], 1.0)'
cent = 'near(x[0],' + str(tol) + ')'
inl = 'near(x[1], 0.0)'
bcu_inflow = DirichletBC(W.sub(0), (0,0), inl)
bcu_cent = DirichletBC(W.sub(0).sub(0), (0.0), cent)
bcu_outflow_u = DirichletBC(W.sub(0).sub(0), u_out, out)
bcu_outflow_p = DirichletBC(W.sub(1), p_out, out)
bcr_inflow = DirichletBC(W.sub(2), 0,inl)
#bcu_inflow_p = DirichletBC(W.sub(1), -1, inl)
bcs = [bcu_inflow, bcr_inflow ,bcu_outflow_u,bcu_outflow_p,bcu_cent]
def div_cyl(v):
return (1/x[0])*((x[0]*v[0])).dx(0) + v[1].dx(1)
colors = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
colors.set_all(0)
CompiledSubDomain(out).mark(colors, 1)
CompiledSubDomain("near(x[0], 1.0)").mark(colors, 3) # wall want ds(3) for free boundary
# Create the measure
ds = Measure("ds", subdomain_data=colors)
rm = (-v[0].dx(0)*x[0]*(-p+2*u[0].dx(0)) - x[0]*v[0].dx(1)*(u[1].dx(0) + u[0].dx(1) ) -x[0]*v[0]*(-p + 2*u[0]/x[0]) )*dx \
zm = (-v[1].dx(0)*x[0]*(u[1].dx(0) + u[0].dx(1)) - v[1].dx(1)*x[0]*(-p+2*u[1].dx(1)) )*dx \
+ (v[1]*x[0]*(-p + 2*u[1].dx(1)))*ds(1)
ce = div_cyl(u)*q*x[0]*dx
fb = (s*(r.dx(0))*x[0] )*dx + ( s*(u[0] - r.dx(1))*x[0] )*ds(3)
F = rm + zm + ce + fb
# Solve problem
solve(F == 0, w, bcs)
# Plot solutions
(u, p,r) = w.split()
import matplotlib.pyplot as plt
tol = 0.001 # avoid hitting points outside the domain
y = np.linspace(0,2, 100)
pe = np.array([( 1,y_) for y_ in y]) # 2D points
xe = pe[:,0]
ye = pe[:,1]
p2 = np.array([p(point) for point in pe])
veloc = np.array([u(point) for point in pe])
r = np.linspace(tol, 1, 100)
c = 0
w_match = np.empty( (100, 2*L) )
w_match[:] = np.nan
u_match = np.copy(w_match)
p_match = np.copy(w_match)
x = np.empty( (100, 10*L) )
x[:] = np.nan
y = np.empty( (100, 10*L) )
y[:] = np.nan
while c < 2*L:
pts = np.array([(r_, c/(L)) for r_ in r])
x[:,c] = pts[:,0]
y[:,c] = pts[:,1]
vel = np.array([u(pt) for pt in pts])
w_match[:,c] = vel[:,1]
u_match[:,c] = vel[:,0]
p_match[:,c] = np.array([p(pt) for pt in pts])
c = c+1
from mpl_toolkits import mplot3d
fig1 = plt.figure()
ax1 = plt.axes(projection='3d')
fig2 = plt.figure()
ax2 = plt.axes(projection='3d')
fig3 = plt.figure()
ax3 = plt.axes(projection='3d')
#Data for a three-dimensional line
c = 0
while c < 2*L:
#w
ax1.plot3D(x[:,c], y[:,c], u_match[:,c],'black')
ax1.plot3D(-x[:,c], y[:,c], u_match[:,c],'black')
#U
ax2.plot3D(x[:,c], y[:,c], w_match[:,c],'black')
ax2.plot3D(-x[:,c], y[:,c], w_match[:,c],'black')
#P
ax3.plot3D(x[:,c], y[:,c], p_match[:,c],'black')
ax3.plot3D(-x[:,c], y[:,c], p_match[:,c],'black')
c = c + 1
# y[ye > 1] = np.nan
ax1.plot3D(-xe, ye, veloc[:,0], 'red' )
ax1.plot3D(xe, ye, veloc[:,0], 'red')
ax2.plot3D(-xe, ye, veloc[:,1], 'red' )
ax2.plot3D(xe, ye, veloc[:,1], 'red')
ax3.plot3D(-xe, ye, p2, 'red' )
ax3.plot3D(xe, ye, p2, 'red')
ax1.set_xlabel('r')
ax1.set_ylabel('z')
# ax.set_zlabel('$\bar u $')
ax3.xaxis.pane.fill = False
ax3.yaxis.pane.fill = False
ax3.zaxis.pane.fill = False
ax3.zaxis.set_rotate_label(False) # disable automatic rotation
ax3.set_zlabel('$\\bar p$',labelpad=25, fontsize = 30)
ax3.set_xlabel('$\\bar r$',labelpad=25, fontsize = 30)
ax3.set_ylabel('$\\bar Z$',labelpad=25, fontsize = 30)
ax3.zaxis.set_tick_params(labelsize=20)
ax3.xaxis.set_tick_params(labelsize=20)
ax3.yaxis.set_tick_params(labelsize=20)
fig3.savefig('p(r)', dpi=400, bbox_inches='tight')
ax2.xaxis.pane.fill = False
ax2.yaxis.pane.fill = False
ax2.zaxis.pane.fill = False
ax2.zaxis.set_rotate_label(False) # disable automatic rotation
ax2.set_zlabel('$\\bar w$',labelpad=25, fontsize = 30)
ax2.set_xlabel('$\\bar r$',labelpad=25, fontsize = 30)
ax2.set_ylabel('$\\bar Z$',labelpad=25, fontsize = 30)
ax2.zaxis.set_tick_params(labelsize=20)
ax2.xaxis.set_tick_params(labelsize=20)
ax2.yaxis.set_tick_params(labelsize=20)
fig2.savefig('w(r)bound', dpi=400, bbox_inches='tight')
ax1.xaxis.pane.fill = False
ax1.yaxis.pane.fill = False
ax1.zaxis.pane.fill = False
ax1.zaxis.set_rotate_label(False) # disable automatic rotation
ax1.set_zlabel('$\\bar u$',labelpad=25, fontsize = 30)
ax1.set_xlabel('$\\bar r$',labelpad=25, fontsize = 30)
ax1.set_ylabel('$\\bar Z$',labelpad=25, fontsize = 30)
ax1.zaxis.set_tick_params(labelsize=20)
ax1.xaxis.set_tick_params(labelsize=20)
ax1.yaxis.set_tick_params(labelsize=20)