1 引言

近年来，随着航空航天、高速铁路以及高速精密机床等领域呈现出高速、精密、环保的发展趋势，磁流体轴承也随之得到广泛应用。磁流体轴承是一种以磁流体作为润滑介质的新型轴承技术，相比于传统轴承，在外加磁场的作用下，磁流体轴承具有良好的自密封性、可控性以及高回转精度、振动小、噪声低、寿命长等性能优点^{[1, 2]}。

磁流体是一种将磁性纳米颗粒经过表面活性处理，并与基液均匀混合以改变流体性质并增强润滑性能的磁性胶体溶液，因此其同时具备流体特性和磁性材料的类似性质^[3]。由于润滑剂中磁性纳米颗粒的作用以及磁场效应的影响，使其表现出非牛顿特性，这主要是由应力偶效应所导致^[4]。

国内外学者对磁流体轴承的研究大多集中在单一影响效应的研究。Wang等^[5]对应力偶流体润滑的动载轴承的性能进行数值研究；Lin^[6]对应力偶流体润滑的有限长轴承挤压膜行为进行理论研究；Osman等^[7]推导出广义的磁流体润滑雷诺方程，并研究不同磁场模型对磁流体轴承运行特性的影响；Shimpi等^[8]通过考虑轴承变形对基于铁磁流体的旋转弯曲粗糙多孔圆板的挤压膜特性进行研究；Hsu等^[9]综合考虑表面粗糙度和磁场的影响，研究磁流体润滑短轴承的性能；Chawla等^[10]通过构建理论模型来研究表面粗糙度效应对应力偶流体润滑的多孔枢轴滑块轴承性能的影响。上述这些研究并没有综合考虑磁场效应、应力偶效应以及表面粗糙形貌效应的组合影响。因此，笔者在前人研究的基础上，通过应力偶效应和粗糙表面形貌效应的双重影响，对非牛顿磁流体轴承的润滑特性进行可靠性预测。

本文首先采用蒙特卡罗方法对高斯随机粗糙表面进行表征，并与实际加工表面形貌图进行对比，同时基于Christensen随机模型和Stokes微连续理论，推导出随机广义非牛顿Reynolds方程，通过有限差分法对微尺度条件下磁流体轴承进行数值模拟，对比研究横向和纵向粗糙模型下不同应力偶参数、磁力系数、动力参数等对其润滑性能的影响。

2 Gauss随机粗糙表面的模拟

基于蒙特卡罗（Monte Carlo）方法对Gauss型随机粗糙面进行模拟。由于实际的工程粗糙表面是由许多不同频率的谐波叠加而成，并且各谐波幅度均可由Gauss分布随机变量获得。因此，首先在频域用功率谱对其进行滤波，再做逆快速傅里叶变换得到粗糙面的起伏高度^{[11, 12]}。

2.1 Gauss型粗糙表面数学模型

以二维粗糙表面建模为例，假设表征的二维粗糙面在x和y方向上的长度分别为L_x和L_y，相邻两点间的间隔为Δx和Δy，离散点数分别为M和N，因此有L_x=M·Δx，L_y=N·Δy，则粗糙面离散点高度可以表示为

$\begin{array}{l} f\left( {{x_m}, {y_n}} \right) = \frac{1}{{{L_x}{L_y}}}\mathop {\mathop {\mathop \sum \limits^{M/2} }\limits_{{m_k} = - M/2 + 1} {\mkern 1mu} }\limits^{} {\mkern 1mu} \mathop {\mathop {\mathop \sum \limits^{N/2} }\limits_{{n_k} = - N/2 + 1} {\mkern 1mu} }\limits^{} {\mkern 1mu} F\left( {{k_{{m_k}}}, {k_{{n_k}}}} \right) \cdot \\ {\rm{exp}}\left[{{\rm{i}}\left( {{k_{{m_k}}}{x_m} + {k_{{n_k}}}{y_n}} \right)} \right] \end{array} $

(1)

式中

$\begin{array}{l} F\left( {{k_{{m_k}}}, {k_{{n_k}}}} \right) = 2{\rm{ \mathit{ π} }}{\left[{{L_x}{L_y}S\left( {{k_{{m_k}}}, {k_{{n_k}}}} \right)} \right]^{1/2}} \times \\ \left\{ \begin{array}{l} \frac{1}{{\sqrt 2 }}\left[{N\left( {0, 1} \right) + iN\left( {0, 1} \right)} \right]{\rm{\;\;\;\;\;\;\;\;\;\;}}{m_k} \ne 0, M/2;{n_k} \ne 0, N/2\\ N\left( {0, 1} \right)\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{m_k} = 0, M/2;{n_k} = 0, N/2 \end{array} \right. \end{array} $

式中，k_mk=2πm_k/L_x, k_nk=2πn_k/L_y, S(k_x, k_y)为二维随机粗糙面功率谱密度，同时为了使f(x_m, y_n)为实数，则需满足条件：F(k_mk, k_nk)=F^*(-k_mk, k_nk)，F(k_mk, -k_nk)=F^*(-k_mk, k_nk)。因此，二维Gauss粗糙面的功率谱密度为

$S\left( {{k_x}, {k_y}} \right) = {\mathit{\zeta }^2}\frac{{{l_x}{l_y}}}{{4{\rm{ \mathit{ π} }}}}{\rm{exp}}\left( { - \frac{{k_x^2l_x^2 + k_y^2l_y^2}}{4}} \right)$

(2)

式中ζ为均方根高度；l_x和l_y分别为x，y方向上的相关长度，其中l_x=l_y时，表面呈各向同性；l_x≠l_y时，表面呈各向异性。

2.2 模型仿真与分析

对二维Gauss随机粗糙面进行模拟时, 分别在X, Y坐标轴上设置采样点个数为M=N=160, 取相关长度l_x=l_y=0.2μm, 均方根高度ζ=0.1μm, 生成的各向同性Gauss粗糙表面如图 1所示。

然而, 由于加工工艺的原因, 实际的工程粗糙表面往往会呈现出各向异性, 使其加工表面具有规则的窄脊和谷的形式, 在对各向异性二维粗糙面进行模拟时, 取均方根高度ζ=0.2μm, l_x和l_y取为0.2μm、0.8μm或0.8μm、0.2μm, 模拟结果如图 2所示, 同时给出了扫描电子显微镜(Scanning electronic microscope, SEM)下加工表面实际粗糙形貌图, 如图 3所示。

经过对比图 2, 图 3可知, 模拟表面和实际加工表面存在相似的纹理结构, 因此, 本文主要研究各向异性粗糙表面纹理结构对磁流体轴承润滑性能的影响, 其模拟结果更能接近于实际工况。

3 数学模型 3.1 物理模型

图 4示出了粗糙形貌下的磁流体轴承物理模型。其中O_b和O_j分别为轴承中心和轴颈中心, 轴颈以角速度ω做逆时针旋转, e为偏心距, φ和θ分别为偏位角和极角, h(θ)为轴承内表面与轴颈外表面之间的间隙, h为润滑膜厚度, h_s为光滑部分膜厚, δ为由随机粗糙引起的膜厚, W为轴颈所受的外载荷, W_ε, W_φ分别为偏心方向和垂直于偏心方向的静载荷分量。

3.2 非牛顿磁流体润滑动载Reynolds方程

根据Stokes的微连续理论，存在应力偶效应的不可压缩流体的连续性方程和动量方程分别为^[13]

$\nabla \cdot \mathit{\boldsymbol{V}} = 0 $

(3)

$\rho \frac{{{\rm{d}}{\boldsymbol{V}}}}{{{\rm{d}}t}} = - \nabla p + {\mathit{\boldsymbol{F}}_{\rm{m}}} + \frac{1}{2}\nabla \mathit{\boldsymbol{B}} + \left( {\mu - \eta {\nabla ^2}} \right){\nabla ^2}\mathit{\boldsymbol{V}} $

(4)

式中矢量V, F_m和B分别为流速、单位体积的体力和体力偶; ρ为流体密度(kg/m³); t为时间(s); p为压力(Pa); μ为动力粘度(Pa·s); η为与应力偶流体物性有关的材料常数。

基于流体动力润滑理论，由于轴承液膜厚度较小，并且忽略惯性力、体力及体力偶，同时把磁场力视作外部力，因此，笛卡尔坐标系下的磁流体控制方程为

$\left\{ {\begin{array}{*{20}{l}} {\frac{{\partial p}}{{\partial x}} = {F_{{\rm{m}}x}} + \mu \frac{{{\partial ^2}u}}{{\partial {y^2}}} - \eta \frac{{{\partial ^4}u}}{{\partial {y^4}}}}\\ {\frac{{\partial p}}{{\partial y}} = 0}\\ {\frac{{\partial p}}{{\partial z}} = {F_{{\rm{m}}z}} + \mu \frac{{{\partial ^2}w}}{{\partial {y^2}}} - \eta \frac{{{\partial ^4}w}}{{\partial {y^4}}}}\\ {\frac{{\partial u}}{{\partial x}} + \frac{{\partial v}}{{\partial y}} + \frac{{\partial w}}{{\partial z}} = 0} \end{array}} \right.$

(5)

式中F_mx和F_mz分别为磁场力的周向和轴向分量, 其表达式为

$\left\{ {\begin{array}{*{20}{l}} {{F_{{\rm{m}}x}} = {\mu _0}{X_{\rm{m}}}{h_{\rm{m}}}\frac{{\partial {h_{\rm{m}}}}}{{\partial x}}}\\ {{F_{{\rm{m}}z}} = {\mu _0}{X_{\rm{m}}}{h_{\rm{m}}}\frac{{\partial {h_{\rm{m}}}}}{{\partial z}}} \end{array}} \right.$

(6)

式中u, v, w分别为x, y, z方向上的速度分量; μ₀为真空磁导率, 取4π×10^-7AT/m; X_m为磁流体的磁化率; h_m为磁场强度。

为了简化计算，假设磁场分布沿轴向呈抛物线状，其磁场控制方程为

${h_{\rm{m}}}\left( z \right) = {h_{{\rm{mo}}}} - \left( {{h_{{\rm{mo}}}} - {h_{{\rm{me}}}}} \right){(z/{L_{\rm{b}}})^2}$

(7)

通过对方程(5)进行积分并引入速度边界条件u(h)=ωR, v(h)=w(h)=0, u(0)=v(0)=w(0)=0, 最终获得修正的Reynolds方程为

$\begin{array}{*{20}{l}} {\frac{\partial }{{\partial x}}\left( {g\left( {h, l} \right)\frac{{\partial p}}{{\partial x}}} \right) + \frac{\partial }{{\partial z}}\left( {g\left( {h, l} \right)\frac{{\partial p}}{{\partial z}}} \right) = 6U\mu \frac{{\partial h}}{{\partial x}} + }\\ {\frac{\partial }{{\partial z}}\left( {g\left( {h, l} \right){F_{{\rm{m}}z}}} \right) + 12\mu \frac{{\partial h}}{{\partial t}}} \end{array}$

(8)

其中, g(h, l)=h³-12l²h+24l³tanh$ \left( \frac{h}{2l} \right)$

式中h_mo和h_me分别为轴承中部和端部磁场强度; L_b为轴承宽度(m); R为轴颈半径(m); U为轴颈表面切向速度分量(m/s), U=ωR; l为应力偶参数, 其中$l=\sqrt{\left( \mathit{\eta }/\mathit{\mu } \right)} $。

3.3 广义的随机非牛顿Reynolds方程

基于Christensen随机粗糙模型，对方程（8）两边取数学期望值，则考虑粗糙模型的随机Reynolds方程为^[14]

$\begin{array}{*{20}{l}} {\frac{\partial }{{\partial x}}E\left( {g\left( {h, l} \right)\frac{{\partial p}}{{\partial x}}} \right) + \frac{\partial }{{\partial z}}E\left( {g\left( {h, l} \right)\frac{{\partial p}}{{\partial z}}} \right) = }\\ {6U\mu E\left( {\frac{{\partial h}}{{\partial x}}} \right) + \frac{\partial }{{\partial z}}E\left( {g\left( {h, l} \right){F_{{\rm{m}}z}}} \right) + 12\mu E\left( {\frac{{\partial h}}{{\partial t}}} \right)} \end{array}$

(9)

式中期望值运算符E(·)被定义为

$E\left( \cdot \right) = \int_{ - \infty }^\infty {\left( \cdot \right)} f\left( \delta \right){\rm{d}}\delta $

(10)

其中f(δ)为随机变量δ的概率密度函数, 由于大多数工程粗糙表面近似服从Gauss分布, 因此可以通过以下多项式近似代替Gauss分布^[15]

$f\left( \delta \right) = \left\{ {\begin{array}{*{20}{l}} {\frac{{35}}{{32{c^7}}}{{({c^2} - {\delta ^2})}^3}, {\rm{\;\;\;\;}} - c \le \delta \le c}\\ {0, {\rm{\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;others}}} \end{array}} \right. $

(11)

式中c是随机油膜厚度变化范围的一半；σ为标准差；其中函数的边界点为c = ±3σ。

根据Christensen随机粗糙模型原理，随机非牛顿Reynolds方程（9）主要取决于粗糙表面纹理结构类型，以下主要针对横向和纵向两种粗糙表面纹理结构进行研究^[16]：

（1）对于横向粗糙模型（Transverse roughness），假设粗糙纹理在轴向存在延伸的窄脊和谷的形式。因此，润滑膜厚度方程和Reynolds方程分别为

$h = {h_{\rm{s}}}\left( {x, z} \right) + \delta \left( {x, \xi } \right)$

(12)

$\begin{array}{*{20}{l}} {\frac{\partial }{{\partial x}}\left( {\frac{1}{{E\left( {1/g\left( {h, l} \right)} \right)}}\frac{{\partial \bar p}}{{\partial x}}} \right) + \frac{\partial }{{\partial z}}\left( {E\left( {g\left( {h, l} \right)} \right)\frac{{\partial \bar p}}{{\partial z}}} \right) = }\\ {6U\mu \frac{\partial }{{\partial x}}\frac{{E\left( {1/{h^2}} \right)}}{{E\left( {1/{h^3}} \right)}} + \frac{\partial }{{\partial z}}\left( {E\left( {g\left( {h, l} \right)} \right){F_{{\rm{m}}z}}} \right) + 12\mu \frac{{\partial E\left( h \right)}}{{\partial t}}} \end{array}$

(13)

（2）对于纵向粗糙模型（Longitudinal roughness），假设粗糙纹理在周向上存在延伸的窄脊和谷的形式。因此，润滑膜厚度方程和Reynolds方程分别为

$h = {h_{\rm{s}}}\left( {x, z} \right) + \delta \left( {z, \xi } \right)$

(14)

$\begin{array}{*{20}{l}} {\frac{\partial }{{\partial x}}\left( {E\left( {g\left( {h, l} \right)} \right)\frac{{\partial \bar p}}{{\partial x}}} \right) + \frac{\partial }{{\partial z}}\left( {\frac{1}{{E\left( {1/g\left( {h, l} \right)} \right)}}\frac{{\partial \bar p}}{{\partial z}}} \right) = }\\ {6U\mu \frac{{\partial E\left( h \right)}}{{\partial x}} + \frac{\partial }{{\partial z}}\left( {E\left( {g\left( {h, l} \right)} \right){F_{{\rm{m}}z}}} \right) + 12\mu \frac{{\partial E\left( h \right)}}{{\partial t}}} \end{array}$

(15)

为了便于计算，令有效角速度^{[17, 18]}${{\omega }^{*}}=\mathit{\omega -2}\frac{\text{d}\varphi }{\text{d}t} $，并引入的无量纲变量及参数为

$\left\{ \begin{array}{l} \theta = x/R, Z = z/{L_{\rm{b}}}, \lambda = {L_{\rm{b}}}/\left( {2R} \right)\\ H = {h_{\rm{s}}}/C = 1 + \varepsilon {\rm{cos}}\theta, \Lambda = c/C\\ \bar P = \bar p{\left( {C/R} \right)^2}/\left( {\mu \omega } \right), L = l/C\\ q = \left( {2/{\omega ^{\rm{*}}}} \right)\left( {{\rm{d}}\varepsilon /{\rm{d}}t} \right), {\delta ^{\rm{*}}} = \delta /C\\ {H_{\rm{m}}} = {h_{\rm{m}}}/{h_{{\rm{mo}}}} = 1 - \left( {1 - \beta } \right){Z^2} \end{array} \right.$

(16)

通过对方程（13），（15）进行量纲一化处理，最终得到横向和纵向粗糙模型的无量化Reynolds方程分别为

$\begin{array}{l} \frac{\partial }{{\partial \theta }}\left( {{G_1}\left( {H, L, \Lambda } \right)\frac{{\partial \bar P}}{{\partial \theta }}} \right) + \frac{1}{{4{\lambda ^2}}}\frac{\partial }{{\partial Z}}\left( {{G_2}\left( {H, L, \Lambda } \right)\frac{{\partial \bar P}}{{\partial Z}}} \right) = \\ 6\frac{\partial }{{\partial \theta }}{G_3}\left( {H, L, \Lambda } \right) + \alpha \frac{\partial }{{\partial Z}}\left[{{G_2}\left( {H, L, \Lambda } \right){H_{\rm{m}}}\frac{{\partial {H_{\rm{m}}}}}{{\partial Z}}} \right] + q{\rm{cos}}\theta \end{array}$

(17)

$\begin{array}{*{20}{l}} {\frac{\partial }{{\partial \theta }}\left( {{G_2}\left( {H, L, \Lambda } \right)\frac{{\partial \bar P}}{{\partial \theta }}} \right) + \frac{1}{{4{\lambda ^2}}}\frac{\partial }{{\partial Z}}\left( {{G_1}\left( {H, L, \Lambda } \right)\frac{{\partial \bar P}}{{\partial Z}}} \right) = }\\ {6\frac{{\partial H}}{{\partial \theta }} + \alpha \frac{\partial }{{\partial Z}}\left[{{G_1}\left( {H, L, \Lambda } \right){H_{\rm{m}}}\frac{{\partial {H_{\rm{m}}}}}{{\partial Z}}} \right] + q{\rm{cos}}\theta } \end{array}$

(18)

其中，

$\begin{array}{l} \left\{ \begin{array}{l} {G_1}\left( {H, L, \Lambda } \right) = 1/T\left( {H, L, \Lambda } \right)\\ {G_2}\left( {H, L, \Lambda } \right) = {H^3} + \frac{1}{3}H{{\rm{\Lambda }}^2} - \\ 12{L^2}H + 24{L^3}K\left( {H, L, \Lambda } \right)\\ {G_3}\left( {H, L, \Lambda } \right) = 1 + \frac{1}{3}\frac{{{\Lambda ^2}}}{{{H^2}}}\\ \alpha = {({h_{{\rm{mo}}}})^2}{\mu _0}{X_{\rm{m}}}{c^2}/\left( {\mu \omega L_{\rm{b}}^2} \right) \end{array} \right.\\ \end{array} $

(19)

$\begin{array}{*{20}{l}} {T\left( {H, L, \mathit{\Lambda }} \right) = \frac{{35}}{{32{\mathit{\Lambda }^7}}} \times }\\ {\int_{ - \mathit{\Lambda }}^\mathit{\Lambda } {\frac{{{{({\mathit{\Lambda }^2} - {\delta ^{\rm{*}}}^2)}^3}{\rm{d}}{\delta ^{\rm{*}}}}}{{{{(H + {\delta ^{\rm{*}}})}^3} - 12{L^2}\left( {H + {\delta ^{\rm{*}}}} \right) + 24{L^3}{\rm{tanh}}\left( {\frac{{H + {\delta ^{\rm{*}}}}}{{2L}}} \right)}}} } \end{array} $

(20)

$K\left( {H, L, \mathit{\Lambda }} \right) = \frac{{35}}{{32{{\rm{\Lambda }}^7}}}\int_{ - \mathit{\Lambda }}^\mathit{\Lambda } {{\rm{tanh}}} (\frac{{H + {\delta ^{\rm{*}}}}}{{2L}}){({\Lambda ^2} - {\delta ^{\rm{*}}}^2)^3}{\rm{d}}{\delta ^{\rm{*}}}$

(21)

式中Z为无量纲轴向坐标; λ为宽径比; C为半径间隙(m), 其中偏心率ε=e/C; H和p分别为无量纲液膜厚度和压力; Λ为粗糙度参数; L为无量纲应力偶参数, 表征润滑剂中添加剂分子链长度; q为动力参数; δ^*为无量纲随机膜厚; H_m为无量纲磁场强度; β为h_mo和h_me两者比值, 取为0.5;α为磁力系数。

为了求解随机广义Reynolds方程，假设沿周向润滑膜破裂初始角度为θ_c，因此，引入无量纲Reynolds边界条件为

$\left\{ {\begin{array}{*{20}{l}} {\bar p\left| {_{Z = \pm 1/2}} \right. = 0}\\ {\bar p\left| {_{\theta = 0} = 0} \right., P\left| {_{\theta = {\theta _{\rm{c}}}} = 0} \right., \frac{{{\rm{d}}\bar p}}{{{\rm{d}}\theta }}\left| {_{\theta = {\theta _c}}} \right. = 0} \end{array}} \right.$

(22)

3.4 数值求解方法

采用有限差分法对随机非牛顿Reynolds方程（17），（18）进行离散求解，如图 5所示，分别沿着周向θ方向和轴向Z方向进行网格划分，网格节点用i, j编号，其压力值用p_{i, j}表示，并对各节点（i，j）利用中心差分近似求导^[19]，则

$\begin{array}{*{20}{l}} {{{\bar p}_{i, j}} = ({A_{i, j}}{{\bar p}_{i + 1, j}} + {B_{i, j}}{{\bar p}_{i - 1, j}} + {C_{i, j}}{{\bar p}_{i, j + 1}} + }\\ {{D_{i, j}}{{\bar p}_{i, j - 1}} - {F_{i, j}} - {G_{i, j}})/{E_{i, j}}} \end{array}$

(23)

其中${A_{i, j}} = {G_{i + 1/2, j}}$，${B_{i, j}} = {G_{i - 1/2, j}} $，${C_{i, j}} = \frac{1}{{4{\lambda ^2}}}{\left( {\frac{{{\rm{\Delta }}\theta }}{{{\rm{\Delta }}Z}}} \right)^2}\cdot{G_{i, j + 1/2}}$，${D_{i, j}} = \frac{1}{{4{\lambda ^2}}}{\left( {\frac{{{\rm{\Delta }}\theta }}{{{\rm{\Delta }}Z}}} \right)^2}{G_{i, j - 1/2}} $，${E_{i, j}} = {A_{i, j}} + {B_{i, j}} + {C_{i, j}} + {D_{i, j}}, {F_{i, j}} = 6{\rm{\Delta }}\theta \left( {{H_{i + 1/2, j}} - {H_{i - 1/2, j}}} \right) + q{\left( {{\rm{\Delta }}\theta } \right)^2}{\rm{cos}}\theta $，${G_{i, j}} = \alpha {\left( {{\rm{\Delta }}\theta } \right)^2}\left[{{G_{i, j}}\left( {{H_{\rm{m}}}\frac{{{\partial ^2}{H_{\rm{m}}}}}{{\partial {Z^2}}} + {{\left( {\frac{{\partial {H_{\rm{m}}}}}{{\partial Z}}} \right)}^2}} \right)} \right] $。

采用超松弛迭代法进行数值求解，并在迭代过程中，当网格节点迭代误差满足以下收敛条件时，则迭代过程终止^[20]。

$\frac{{\mathop {\mathop {\mathop \sum \limits^{m - 1} }\limits_{i = 2} {\mkern 1mu} }\limits^{} {\mkern 1mu} \mathop {\mathop {\mathop \sum \limits^{n - 1} }\limits_{j = 2} {\mkern 1mu} }\limits^{} {\mkern 1mu} \left| {\bar p_{i, j}^k - \bar p_{i, j}^{k - 1}} \right|}}{{\mathop {\mathop {\mathop \sum \limits^{} }\limits_{i = 2} {\mkern 1mu} }\limits^{m - 1} {\mkern 1mu} \mathop {\mathop {\mathop \sum \limits^{} }\limits_{j = 2} {\mkern 1mu} }\limits^{n - 1} {\mkern 1mu} \left| {\bar p_{i, j}^k} \right|}} \le \gamma $

(24)

$\bar p_{i, j}^k = \tau p_{i, j}^k + \left( {1 - \tau } \right)p_{i, j}^{k - 1}$

(25)

式中p_{i, j}^k为新压力；p_{i, j}^k-1为当前压力；p_{i, j}^k为本轮迭代得到的压力；γ为收敛精度，本文取10^-4；τ为加速系数，一般取1<τ<2。

最后输出液膜内部节点无量纲压力分布，从而通过式（26）得到液膜无量纲承载力和偏位角

$\left\{ \begin{array}{l} W_{\rm{ \mathit{ ε} }}^{\rm{*}} = 2\int_0^{{L_b}/2} {\int_0^{{\theta _c}} {\bar p{\rm{cos}}\theta {\rm{d}}\theta {\rm{d}}Z} } \\ W_{\rm{ \mathit{ φ} }}^{\rm{*}} = 2\int_0^{{L_b}/2} {\int_0^{{\theta _c}} {\bar p{\rm{sin}}\theta {\rm{d}}\theta {\rm{d}}Z} } \\ {W^{\rm{*}}} = \sqrt {W{{_{\rm{ \mathit{ ε} }}^{\rm{*}}}^2} + W{{_{\rm{ \mathit{ φ} }}^{\rm{*}}}^2}}, \varphi = {\rm{ta}}{{\rm{n}}^{ - 1}}\left( { - \frac{{W_{\rm{ \mathit{ φ} }}^{\rm{*}}}}{{W_{\rm{ \mathit{ ε} }}^{\rm{*}}}}} \right) \end{array} \right.$

(26)

式中W_ε^*，W_φ^*分别为偏心方向和垂直于偏心方向的无量纲静载荷分量；W^*为无量纲承载力；φ为偏位角。

3.5 模型验证及精度选择

为了验证本文数值模拟的准确性及可靠性，并且考虑到磁流体润滑理论尚不完善以及自身实验条件的限制，暂选取文献[2]中模拟结果进行对比验证。如图 6所示，为了验证算法的有效性，以横向粗糙模型为例，计算过程中关键参数的设置与文献中保持一致，并通过不断调整收敛精度γ，从而保证模拟结果的准确性以及提高运算效率。

从上图可以发现，当收敛精度γ = 10^-4时，计算结果与文献结果具有较好的一致性，且最大误差为6.13%，因此本文算法有效，数值模拟较为准确。

4 计算结果与分析 4.1 粗糙纹理下不同参数对无量纲压力分布的影响

选取偏心率ε=0.5，应力偶参数L=0.1，动力参数q=0.2，磁力系数α=0.2时，计算得到粗糙度参数Λ分别为0.0、0.3和宽径比λ分别为0.5，1.0，2.0时轴承中截面周向无量纲压力及最大压力平面处轴向无量纲压力分布曲线，如图 7所示。

从图 7可知，粗糙纹理结构类型并不会改变液膜压力分布规律，其中最大无量纲压力随宽径比的增加而显著增大，同时横向和纵向粗糙模型对压力分布的影响趋势将发生变化。当宽径比较小时，纵向粗糙模型能够增强液膜整体压力水平，而横向粗糙模型将削弱压力水平；当宽径比较大时，纵向和横向粗糙模型的影响趋势会发生反转。由图 7（a）知道，当λ=0.5时，相比于光滑表面的情况，纵向粗糙模型产生的最大压力值将提升3.9%；当λ=2.0时，横向粗糙模型产生的最大压力值将提升5.7%。

选取宽径比λ=2.0，偏心率ε=0.5，计算所采用的磁力系数α分别为0.2，0.5，动力参数q分别为0.5、1.0，应力偶参数L分别为0.1，0.2，0.3，粗糙度参数Λ分别为0.0，0.3。求解得到不同粗糙模型下各参数对周向及轴向无量纲压力分布的影响如图 8，9所示。

从图 8，9可以发现，磁力系数、动力参数以及应力偶参数的增加均能不同程度的增强液膜压力水平，同时并不会改变轴向压力的分布规律，使得轴向压力呈抛物线状。从图 8（a）可知，随着磁力系数α的增大，液膜最大压力显著增加，同时横向和纵向粗糙模型的影响效应也显著增强；以α分别为0.2和0.5为例，相比于光滑表面，横向粗糙模型产生的最大压力分别提升33.2%，44.9%，纵向粗糙模型产生的最大压力分别降低24.9%，31.2%。从图 8（b）可知，动力参数q增加时，横向和纵向粗糙模型的影响趋势相反，影响程度几乎相同。当θ<1.46rad，其压力值几乎无变化；当θ>1.46rad时，随着动力参数的增加，其压力峰值点向右移动。从图 8（c）发现，应力偶参数L对液膜压力水平的影响更加显著，随着应力偶参数的增加，不同粗糙模型对压力影响效应增强；相比于粗糙模型，应力偶效应将起到主导作用。

4.2 粗糙纹理下不同参数对无量纲承载力和偏位角的影响

当偏心率ε=0.5，磁力系数α=0.2，动力参数q=0.2时，计算时所采用应力偶参数L分别为0.2，0.5，粗糙度参数Λ分别为0.0，0.5，求解得到无量纲承载力随宽径比的变化规律如图 10所示。

从图上可以看出，随着宽径比的增加，液膜无量纲承载力有所增大，同时应力偶效应将显著增强轴承的承载性能；当宽径比较小或较大时，横向和纵向粗糙模型对承载力的影响效应相反，以L=0.2为例，当λ<1.70时，纵向粗糙模型将产生较大的承载力；当λ > 1.70时，横向粗糙模型则产生较大的承载力。然而，当应力偶参数变化时，承载力随宽径比变化的分界点也会变化，例如L=0.5时，分界点右移至λ = 1.95。

当宽径比λ=2.0，偏心率ε=0.5时，求解得到不同参数对液膜承载力及偏位角的影响，如图 11，12所示。

从图 11（a）知道，无量纲承载力随偏心率的增加而明显增大，尤其在高偏心率、高应力偶参数情况下；相比于光滑表面情况，两种粗糙模型对承载力的影响程度不可忽略，例如当ε=0.5，L=0.4时，相比于光滑情况，横向粗糙模型使承载力提升9.9%，纵向粗糙模型使其降低9.3%。从图 11（b）可以发现，磁力系数和应力偶参数越大，承载力增加幅度越显著，粗糙模型对承载力的影响效应也就越强。

从图 12（a）可知，轴承偏位角随着偏心率的增加先增大后减小，同时粗糙模型也会影响偏位角的变化；图 12（b）显示，应力偶参数和磁力系数越大，轴承偏位角也就越小，同时有效承载力也会显著提高。

5 结论

本文对粗糙表面形貌下的非牛顿磁流体轴承进行数值研究，得出以下结论：

（1）大多数实际工程粗糙表面服从高斯分布，由于加工工艺的原因，其表面会呈现出各向异性；基于Christensen随机模型来研究横向和纵向粗糙纹理对磁流体轴承性能的影响，其能够获得较为准确的可靠性预测结果。

（2）不同粗糙模型对承载性能的影响趋势与宽径比有关；当粗糙度参数为0.3，轴承宽径比为0.5时，纵向粗糙模型将使最大无量纲压力提升3.9%，而横向粗糙模型将会削弱压力水平；当轴承宽径比增至2.0时，横向粗糙模型将使其提升5.7%从而改善承载性能。

（3）动力参数、磁力系数以及应力偶参数的增加均能不同程度地增强液膜压力水平，其中后两者的增加将使粗糙模型的影响效应显著增强，同时承载力增大且偏位角减小；当应力偶参数从0.2增至0.5时，粗糙模型对承载力的影响趋势随宽径比变化的分界点从1.70增至1.95。

（4）综合考虑磁流体应力偶效应（非牛顿性）及轴承轴颈内外表面粗糙形貌效应，对磁流体轴承进行研究更加接近于实际工况，相比于光滑表面牛顿流体的润滑情况，在合理的控制下，粗糙纹理结构下的非牛顿磁流体轴承将显著改善轴承的润滑性能。