I want to model the flow through a cylindrical air bearing such as the New Way Flat round air bearing, so that I can play around with modifying parameters (porous and fluid medium type and geometry) and computing physical parameters (stiffness, pressure distribution, flow).
I've chosen a 2D axisymmetric approach, modeling half of the cross-section through the axis of rotation. The cross-section extends from the axis of rotation to the outer shell of the cylinder. On this section a rectangular grid of elements is defined. Each element has a permeability with its direct neighbours (up, down, left, right) defined, and each element has a pressure.
The basic idea is that, for each element, flow in and out is determined by its permeability and pressure difference with its neighbours, under the constraint of flow continuity (no accumulation). On two points, the inlet and outlet, the pressures are defined.
The 'bottom' row of elements differ from the other elements. They represent the flow in the gap between the sliding surface of the air bearing and the contra-surface.
The permeabilities in the porous medium are modeled as Darcy flow, the flow in the gap as plane Poiseuille flow.
The resulting equations are collected and solved as a matrix-vector equation. Code below.
A few observations: The result give reasonable (to me) answers on bearing load and flow.* Example: 250 N for a 50 mm OD bearing and 7078 N for a 250 mm OD bearing. Ideal load as specified: 330 N and 12233 N respectively.
*when tuning permeability so that the nominal flow works out. I don't have a reference value for porous carbon and depending on number of elements.
My questions:
- There is zero sensitivity of the load bearing capacity to changing the gap, whereas the contact in reality is stiff. Example: a 50 mm OD bearing with nominal load bearing capacity of 333 N and a gap height of 5 µm has a stiffness of 58 N/µm. My code computes a 282 N load bearing capacity at 5 µm gap, but also at 50 µm and 500 µm gap. It needs a gap of 5000 µm for the load bearing capacity to change to 216 N. This is obviously non-sensical. Any suggestions?
Look at the permeability in horizontal direction, including in the bottom row the flow in the gap, one sees in the bottom row values that are orders of magnitude lower than in the bulk material. I interpret that as a 'layer of insulation', but I struggle to then visualized where the flow is supposed to go. I would figure through the gap, but is the permeability of the bulk is 10000x higher, it would go through there. Not what I expected. See pic below.
- There is a strong dependency of the outcome to the number of nodes. At lower numbers of elements (<10 per direction), the load bearing capacity and flow can be twice of that at higher number (>30 per direction. Can you help me understand?
import numpy as np
import math
import matplotlib.pyplot as plt
r = 0.025
p_inlet = 5e5
p_outlet = 1e5
kmat = 6e-10 # permeability
t = 5e-6
dynvisc = 1.8e-5
z = 0.01 # length of bearing
# number of nodes
nr = 20
nz = 20
rs = np.linspace(0, r, nr, endpoint=False)+0.5*r/nr
# problem formulation
# {flow} = [c] * {p}
# right hand side. BC and flow continuity
RHS = np.zeros((nz*nr,1))
# coefficient matrix
c = np.zeros((nr*nz, nr*nz))
klm = np.zeros((nz, nr))
krm = np.zeros((nz, nr))
kum = np.zeros((nz, nr))
kdm = np.zeros((nz, nr))
areaslm = np.zeros((nz, nr))
areasrm = np.zeros((nz, nr))
areasum = np.zeros((nz, nr))
areasdm = np.zeros((nz, nr))
# permeabilities
for i in range(nz):
for j in range(nr):
kl = 0
kr = 0
ku = 0
kd = 0
arear = 0
areal = 0
areau = 0
aread = 0
# flow from/to right
if j==nr-1:
# outer cylindrical surface of porous medium, blocked
kr = 0
else:
arear = 2*math.pi*(rs[j]+0.5*r/nr)*z/nz
if i==nz-1:
# gap flow
kr = -t**3/(12*dynvisc*r/nr)
else:
kr = -kmat/(dynvisc*r/nr)
krm[i,j] = kr
areasrm[i,j] = arear
# flow from/to left
if j==0:
# axis of rotation, blocked
kl = 0
else:
areal = 2*math.pi*(rs[j]-0.5*r/nr)*z/nz
if i==nz-1:
# gap flow
kl = -t**3/(12*dynvisc*r/nr)
else:
kl = -kmat/(dynvisc*r/nr)
klm[i,j] = kl
areaslm[i,j] = areal
# flow from/to up
if i==0:
# top row is blocked
pass
else:
areau = math.pi*((rs[j]+0.5*r/nr)*(rs[j]+0.5*r/nr)-(rs[j]-0.5*r/nr)*(rs[j]-0.5*r/nr))
ku = -kmat/(dynvisc*z/nz)
kum[i,j] = ku
areasum[i,j] = areau
# flow from/to down
if i==nz-1:
# bottom row is bearing surface, blocked
pass
else:
aread = math.pi*((rs[j]+0.5*r/nr)*(rs[j]+0.5*r/nr)-(rs[j]-0.5*r/nr)*(rs[j]-0.5*r/nr))
kd = -kmat/(dynvisc*z/nz)
kdm[i,j] = kd
areasdm[i,j] = aread
#fig, axs = plt.subplots(2,2)
#axs[0,0].imshow(krm, cmap='hot', interpolation='nearest')
#axs[1,0].imshow(klm, cmap='hot', interpolation='nearest')
#axs[0,1].imshow(kum, cmap='hot', interpolation='nearest')
#axs[1,1].imshow(kdm, cmap='hot', interpolation='nearest')
#
#fig1, axs1 = plt.subplots(2,2)
#axs1[0,0].imshow(areasrm, cmap='hot', interpolation='nearest')
#axs1[1,0].imshow(areaslm, cmap='hot', interpolation='nearest')
#axs1[0,1].imshow(areasum, cmap='hot', interpolation='nearest')
#axs1[1,1].imshow(areasdm, cmap='hot', interpolation='nearest')
for i in range(nz):
for j in range(nr):
# central node of stencil
c[i+j*nz, i+j*nz] += krm[i,j]*areasrm[i,j]
c[i+j*nz, i+j*nz] += klm[i,j]*areaslm[i,j]
c[i+j*nz, i+j*nz] += kum[i,j]*areasum[i,j]
c[i+j*nz, i+j*nz] += kdm[i,j]*areasdm[i,j]
# right
if j==nr-1:
pass
else:
c[i+j*nz, i+(j+1)*nz] -= krm[i,j]*areasrm[i,j]
# left
if j==0:
pass
else:
c[i+j*nz, i+(j-1)*nz] -= klm[i,j]*areaslm[i,j]
# up
if i==0:
pass
else:
c[i+j*nz, i-1+j*nz] -= kum[i,j]*areasum[i,j]
# down
if i==nz-1:
pass
else:
c[i+j*nz, i+1+j*nz] -= kdm[i,j]*areasdm[i,j]
K = c.copy()
c[0,:] = np.zeros((1,nz*nr))
c[0,0] = 1
RHS[0,0] = p_inlet
c[-1,:] = np.zeros((1,nz*nr))
c[-1,-1] = 1
RHS[-1,0] = p_outlet
p = np.linalg.solve(c, RHS)
p_mtrx = p.reshape([nr,nz]).transpose()
Q = np.dot(K,p)
Q_mtrx = Q.reshape([nr,nz]).transpose()
F = np.dot(areasum[-1,:], p_mtrx[-1,:].T)
print('load: {:0.0f} N'.format(F))
print('flow: {:0.1E} NL/min'.format(Q[0,0]*1000/60))
#plt.imshow(p_mtrx, cmap='hot', interpolation='nearest')
#plt.imshow(Q_mtrx, cmap='hot', interpolation='nearest')
#plt.plot(rs, p_mtrx[-1,:])
#plt.imshow(c[0,:].reshape([nr,nz]).transpose(), cmap='hot', interpolation='nearest')
