Pressure Singularity - Boundary layer problem

mathematics
#1

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

p(r).PNG
p(r).PNG3968×3704 1.12 MB

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

Hello @Niall_hanevy I’m getting a similar problem.

For axisymmetric flows I think the boundary conditions should be u=0 and \partial_r w=0. These make sense physically (u is non physical at r=0 and w must be an even function with respect to r, requiring \partial_rw=0 to obtain a smooth function).

I derived my boundary conditions specifically by looking at the terms of highest order in 1/r and forcing them to cancel at r=0 (mass and r momentum yield u=0, \partial_rw=0 comes from Z momentum). From your code, I think you impose the same.

This was enough for me in the axisymmetric case (recall FENICs imposes free flow implicitly because of the IBP, so \nabla U\cdot n=pn\Leftrightarrow \partial_r U=pe_r\Rightarrow \partial_rw=0 naturally).

By the way there’s a typo for your pressure term in the Z direction, it has a minus sign too.

I think the problem flows from the fact that the coordinate transformation is not enough to cancel out all boundary terms at r=0. The diffusion term still reads something like (U=(u\;w)^T velocity vector, v vector of test functions) :
\frac{1}{Re}(\nabla U,\nabla v)=\frac{1}{Re}\int\frac{u}{r}\frac{v_r}{r}rdrd\theta dz+...=\frac{1}{Re}\int\frac{uv_r}{r}drd\theta dz+...

It still contains a 1/r term, so it’s still singular. I tried using a transformation v\rightarrow rv to fix this issue, but it seems UFL does not properly do the symbolic cancellation of the various r terms. Personally, I’m still getting errors in a slightly more complicated azimuthal decomposition case because of this.

#4

Hi,

So my weak formulation was the source of my error.
The term -\sigma_{\theta \theta}/r in the r- momentum equation above is

x[0] v[0]*(-p + 2*u[0]/x[0]) *dx
whereas is should be
-v[0]*(-p + 2*u[0]/x[0]) *dx
So I didn’t cancel the 1/r with the coordinate system transformation r dr dz.

This might not solve your problem but I hope it helps.
I see the typo for p in my z equations. Apologies but the code is fine.

1 Like