0 引　言

船舶轴系不但需要承担自身的重量，还需要承受来自船舶柴油机的负荷、船舶螺旋桨的推力以及船体本身的振动变形等动态以及静态的载荷^[1]。除此之外，随着船舶技术快速发展，船体不断地大型化，例如船舶轴承直径和长度这些尺寸数据不断的增大，这进一步提升了船舶轴承承受的载荷的复杂性，所以需要监测船舶轴承负荷情况，以避免船舶故障的发生^[2]。船舶轴系工作状态通过轴承动载荷来体现，作为船舶轴系工作状态的重要指标，轴承动载荷不但能够对船舶轴系的工作状态进行检测，还能够对船舶的轴系校中质量进行评价，因此对船舶的轴系的性能评价有着重要意义^[3]。由于缺乏测量工具，当前对于船舶轴承动载荷的测量工作进行得比较少^[4]。本文借助有限元仿真技术，对船舶的轴承承载进行仿真分析，这有助于我国船舶轴承技术的快速发展。

1 有限元仿真技术 1.1 有限元计算理论

有限元技术以Newmark法为基础，对刚度矩阵进行转置操作，以此来达到增量转换的目的，并且模型的线性解算通过逼近原理来实现。在隐式算法中需要对刚度矩阵进行逆运算，因此要求刚度矩阵为非奇异矩阵，否则，在大变形情况下，计算结果无法收敛^[5]。对于高动态模型，通常基于微积分对增量步进行解算，以便能够得到模型的动态状况。当模型处于运动状态的时候，其全局系统下的离散方程如下式：

$ {\boldsymbol{M}}{\ddot u_n} + K\left( {{u_n}} \right){u_n} = {R_n}。$

(1)

式中：u_n为位移矢量；M为质量矩阵。基于中心差分变换，能够得到下式：

$ {u_{n + 1}} = {u_n} + \Delta {t_{n + 1}}{\dot u_n} + \frac{1}{2}\Delta t_{n + 1}^2{\ddot u_n}\text{，} $

(2)

$ {\dot u_{n + 1}} = {\dot u_n} + \frac{1}{2}\Delta {t_{n + 1}}\left( {{{\ddot u}_{n + 1}} + {{\ddot u}_n}} \right)\text{，} $

(3)

$ K_nu_n=\sum\limits_i\int_{V^{\left(l\right)}}B_n^{\mathrm{T}}S_n\mathrm{d}V。$

(4)

根据上述公式可得：

$ \left( {\frac{1}{{\Delta {t^2}}}M} \right){u_{n + 1}} = {R_n} - \sum\limits_i {F_n^i} - \frac{1}{{\Delta {t^2}}}\left( {{u_{n - 1}} - 2{u_n}} \right)\text{。} $

(5)

式中：$F^i_n $为节点的反作用力。

基于式（5），则可以获得形变位移矩阵在t时刻的数学表达式：

$ _0^tB_L^{\left(k\right)}=_0^tB_{L_00}^{\left(k\right)t}X^{\mathrm{T}}。$

(6)

通过对式（6）中的变形梯度${}^t_0X $进行高斯积分解算，则可以获得节点的反作用力，其计算方法如式（7）所示。通常采用显示求解的方法对非线性方程进行解算，这样则无需对方程的收敛问题进行考虑，但是由于该方法的仿真步长很小，因此对于很大的模型而言，该方法需要的时间比较长，所以一般采用显示求解的方法对模型的瞬态问题进行解算。

$ {}^t{\hat F^{\left( m \right)}} = \int_{{0_V}} {{}_0^t} B_L^{{\mathrm{T}}_0^t}\tilde S{{\mathrm{d}}^0}V。$

(7)

在进行有限元分析的时候，一般采用节点接触法计算离散化的接触问题，该方法是把单边接触转变为面接触，利用离散化的方法进行解算，如下式：

$ {\hat X^m}\left( \zeta \right) = X_1^m + \left( {X_2^m - X_1^m} \right)\zeta 。$

(8)

归一化处理之后的切线矢量的计算方法如下式：

$ {t_m} = \frac{1}{e}\hat X (\xi ) 。$

(9)

式（9）中的$\hat X$(ζ)和e的计算方法如下式：

$ \mathop X\limits^ \wedge (\xi ) = X_2^m - X_1^m \text{，} $

(10)

$ e = \left\| {X_2^m - X_1^m} \right\|。$

(11)

从面节点和主面之间的最小间隔：

$ {g_N} = \left[ {{X^s} - \left( {1 - \bar \zeta } \right)X_1^m - \bar \zeta X_2^m} \right]{n_m}。$

(12)

式中：n_m为单元的法向矢量。并且从面节点的投影坐标计算方法如下：

$ \bar \zeta = \frac{{\left( {{X^s} - X_1^m} \right){t_m}}}{e}。$

(13)

将点和面之间的接触转变成平衡方程的模型，如式（14）所示，式中a_N为罚函数系数。同时法向力矢量可以通过式（15）来计算。

$ {t_N} = {a_N}{g_N}\text{，} $

(14)

$ t = {t_N}{N_m}\text{。} $

(15)

1.2 材料模型理论

船舶的轴承是由高强度以及高模型的刚性材料构成的，并且不容易变形^[6]。因此为了能够得到精确的船舶轴承模型特征，则需要对船舶轴承的有限元模型匹配最精确的模型参数，这是船舶轴承有限元解算的关键步骤，对后续有限元仿真的结果的精确有着决定性的作用^[7]。船舶轴承的应变数学模型可以用下式表示：

$ S = \frac{{\partial U\left( F \right)}}{{\partial F}}。$

(16)

式中：S为形变量。船舶轴承的本构方程有多种形式，其广义的表现形式如下式：

$ U = \sum\limits_{i + j}^N {{C_{ij}}{{\left( {{{\bar I}_1} - 3} \right)}^i}{{\left( {{{\bar I}_2} - 3} \right)}^j} + \sum\limits_{i = 1}^N {\frac{1}{{{D_i}}}{{\left( {{J_{el}} - 1} \right)}^{2i}}} } 。$

(17)

式中：C_ij为剪切性能；D_i为可压缩性系数。

以熵统计理论为基础，可以构建出船舶轴承的本构模型，其剪切模型的表达式如下式：

$ U' = \mu \left[ {\lambda _m^2\ln \left( \eta \right)\alpha \left( {\frac{{\bar I}}{2}} \right)} \right] + \frac{1}{D}\left[ {\frac{{J_{el}^2 - 1}}{2} - \ln {J_{el}}} \right]。$

(18)

式中：$\bar I$和η的计算方法分别如式（19）和式（20）所示；μ的数值大小和应力方向相关；λ_m的数值大小和应变方向相关；α的数值大小和模型的形状相关；β的数值大小和形状的变化相关。

$ \bar I = \left( {1 - \beta } \right){\bar I_1} + \beta {\bar I_2}\text{，} $

(19)

$ \eta = \sqrt {\frac{{\bar I - 3}}{{\lambda _m^2 - 3}}}。$

(20)

船舶轴承长时间工作之后会出现Mullins效应。在Mullins效应方程中，应变能密度的计算方法如下式：

$ W\left( F \right) = \left( {1 - d} \right){W_0}F。$

(21)

式中：F为船舶轴承的形变梯度；d为船舶轴承的损伤系数，其数值大小如下式：

$ d = \left\{ {\begin{array}{*{20}{l}} {d\left( x \right)}，{x > 0} ，\\ {d\left( {\max \left[ {x\left( s \right)} \right]} \right)}，{{\mathrm{others}}} 。\end{array}} \right. $

(22)

船舶轴承的Mullins效应可以通过引入损伤参数κ来模拟，其计算方法如下式：

$ \kappa = 1 - \frac{1}{r}erf\left[ {\frac{1}{m}\left( {{W_m} - \bar W\left( {{\lambda _1},{\lambda _2}} \right)} \right)} \right]\text{。} $

(23)

船舶轴承应力应变拟合曲线如图1所示。可知，随着应变的变大，船舶轴承受到的应力也变大。

2 船舶轴承动力学仿真分析

在仿真过程中，船舶轴承的转动利用电机来驱动，并且电机的驱动扭矩通过估算的方法确定。本文利用吸工装置来吸收轴承传递给螺旋桨的驱动力，以此来替代螺旋桨。激励频率一般通过轴频和电机级数相乘得到，并且电机在工作过程中的扭振激励可以表示为：

$ {M_m} = \beta {M_e}\sin \left( {p\omega t + \varepsilon } \right)。$

(24)

船舶螺旋桨的激振扭矩可以以函数级数的形式来描述，并且函数的基频可以表示为：

$ {M_x} = {M_0} + \sum\limits_{k = 1}^\infty {{M_{k{Z_p}}}\sin \left( {k{Z_p}\omega t + {\varepsilon _{k{Z_p}}}} \right)}。$

(25)

考虑到式中的相位角取决于船舶尾部的伴流系数，但是船舶尾部的伴流系数很难确定，因此螺旋桨的激振扭矩的经验计算式为：

$ {M_{{Z_p}}} = \beta {M_0}。$

(26)

船舶螺旋桨的激励扭矩计算式为：

$ {M_{v{z_p}}} = 9549.3\beta \frac{p}{{v{n_p}}}{\left( {\frac{{vn}}{{v{n_p}}}} \right)^2}\text{。} $

(27)

利用有限元仿真软件Ansys对船舶轴承动力学模型进行构建。在50 r/min工况下，船舶轴承的仿真结果如图2～图4所示。由图2可知，船舶的角加速度在0.0125～0.01565°/s²之间波动。由图3可知，船舶的角速度在300°/s上下波动。可以看出，船舶的轴承角速度并没有无限制的往上涨，而是维持在一定的数值上下波动，这是因为船舶轴承在转动过程中需要克服阻力做功。由图4可知，受力最大值可达37 N。

3 船舶轴承承载仿真分析

模型中钢的参数主要涉及轴承密度、泊松比等，设置的划分大小为100 mm，并且对船舶轴承模型添加相应的约束以及载荷，包括轴承的约束、螺旋桨的驱动力以及重力。添加完船舶轴承系统的外部激励之后，船舶轴承振幅随频率变化的曲线如图5所示，可以看出最大振幅可达0.35 mm。

任何机械结构均存在固有频率，通常情况下机械结构的固有频率和其自身的质量以及刚度相关。船舶轴承的固有频率的计算方法如式（28）所示，可知，影响船舶轴承固有频率的因素主要有螺旋桨质量和轴承的刚度。仿真得到的船舶纵振频率随轴承刚度的变化曲线如图6所示。可知，纵振频率随着轴承钢度的增大而变大。

$ \omega = \frac{1}{{2{\text{π}}}}\sqrt {\frac{{{K_{th}}}}{m}}。$

(28)

图7为船舶轴承垂向承受的最大压力随转速的变化曲线。可知，在船舶轴承转速为150 r/min时，船舶轴承的垂向受力的最大值出现了最小的情况。

图8为船舶轴承受力随轴承弹性模量的变化曲线。可知，随着弹性模量的增大，船舶轴承受到的应力也随之变大。

4 结　语

船舶的轴系负荷太大或者振动太强的时候，船舶的轴承会出现部分磨损甚至损坏，严重的情况下可能会出现船舶轴承断裂。各类船舶朝着大型化的方向发展，这使得船舶轴承的各项尺寸均会变大，并且随着船舶尾部的刚度降低，船舶轴承受到的应力负荷将会变得更加复杂，因此船舶轴系的振动会产生更加严重的后果。本文基于有限元仿真技术，对船舶轴承承载进行了仿真，这有助于促进我国船舶轴承技术的快速进步。