0 引　言

船舶航行过程中，推进轴系作为船尾主要旋转部件，轴系旋转及其驱动机构（如螺旋桨、轴传动装置等）会受到各种激励力，从而产生复杂振动，研究表明，90%以上船舶尾部振动来自动力装置、螺旋桨以及轴系自身耦合激励^[1−2]。

目前，针对船舶推进轴系的研究主要集中在确定方法下的轴系振动特性和减振方面。尹红升等^[3]基于有限元方法建立了弹性支撑的桨－轴－筏架－基座模型，研究分析了弹性支撑结构对轴系振动的能量传递和衰减；何江洋等^[4]基于传递矩阵法建立推力轴承边界条件改变后的轴系纵振模型，研究集成减振系统设计参数对传递特性的影响规律；Lai等^[5]以轴承垂向位移为优化目标，建立基于遗传算法的轴系多学科多目标优化方法。但随着船舶尺寸的增大和载荷的增加，船舶推进轴系振动的不确定性显著影响船舶航行过程中轴系的关键性能和可靠性，因此，考虑工程结构的不确定因素以及轴系激励的随机性对船舶推进轴系的精确动态评估显得尤为重要。

在船舶推进轴系运行过程中，各种来源都会产生不确定性。从传统的物理角度来看，这些不确定性主要存在于模型参数不确定性、边界条件不确定性、外部载荷不确定性、故障不确定性和模型不确定性^[6]。Kang等^[7]研究了由入流湍流引起的弹性螺旋桨轴系系统的随机振动和力传递特性；Zhang等^[8]开发了一个随机模型，用于受螺旋桨力和轴承刚度不确定性影响的船舶轴系的不确定性量化；Shang等^[6]基于小样本条件下的时变不确定激励数学模型和相应的不确定振动响应分析方法，分析螺旋桨的不确定性振动响应；Gaidai等^[9]基于路径积分法研究了6自由度轴系横向随机振动问题；Pietkiewicz等^[10]针对滑动轴承支撑的转子系统，研究了随机载荷作用横向响应问题；Ma等^[11]利用具有未知确定性系数的非侵入式广义多项式混沌展开用于表示转子动力学中不确定性的传递，分析了考虑随机分布特性的参数不确定性对转子系统随机振动特性的影响；Lu等^[12]基于随机矩阵理论，将多源不确定性纳入轴系动力学方程，研究结果表明，考虑不确定性时，时域波形在规定范围内波动，不平衡和错位加剧，导致稳定性差。

尽管以上文献已对船舶轴系的随机轴向振动问题开展研究，但更多关注的是螺旋桨推力或柴油机曲轴某个单方面不确定因素对轴系轴向振动影响，而且忽略了曲轴刚度和轴承刚度与阻尼对振动的影响。

1 理论基础 1.1 轴系离散模型

在计算船舶轴系扭振和轴向振动时，实际将轴系简化成一个由一系列质量-弹簧元件组成的离散模型，该模型每个元件的质量体现轴系振动的惯性特性、轴向刚度体现轴系振动弹性特性，船舶推进轴系和离散模型如图1所示。

在这种情况下，作为激励力和轴向共振振动计算对象的轴系动态特性是由其质量相对于其在激励力下的静态平衡位置随时间运动决定。为了描述船舶推进轴系运动，使用第二类拉格朗日方程得到以下二阶恒定系数微分方程组的矩阵形式：

$ \left[ m \right]\left\{ {\ddot x} \right\} + \left[ B \right]\left\{ {\dot x} \right\} + \left[ C \right]\left\{ x \right\} = \left\{ T \right\} 。$

(1)

其中：[m]和[B]为质量和耗散对角线矩阵；[C]为二次刚度矩阵；$ \left\{ x \right\} $、$ \left\{ {\dot x} \right\} $、$ \left\{ {\ddot x} \right\} $为广义坐标、速度和加速度列矩阵；$ \left\{ T \right\} $为激励力矩阵。

由于以上二阶微分方程没有标准解，因此，目前普遍求解方式是将方程(1)转化为一阶微分方程，并逐步实现四阶Runge-Kutt的标准求解，求解过程如下：

首先，用x=y₁和$ \dot x $=y₂代换，得到由2个一阶方程组成的系统：

$ \left\{ \begin{aligned} &\left\{ {{{\dot y}_1}} \right\} = \left\{ {{y_2}} \right\}，\\ &[m]\left\{ {{{\dot y}_2}} \right\} + [B]\left\{ {{y_2}} \right\} + [C]\left\{ {{y_1}} \right\} = T(t) 。\\ \end{aligned} \right. $

(2)

通过正则化求解，轴系某一位置的位移和速度分别是：

$ \left\{ \begin{gathered} {{\dot y}_1} = {F_1}\left( {t,{y_1},{y_2}} \right)，\\ {{\dot y}_2} = {F_2}\left( {t,{y_1},{y_2}} \right)。\\ \end{gathered} \right. $

(3)

$ \left\{ \begin{gathered} {F_1}\left( {t,{y_1},{y_2}} \right) = {y_2}，\\ {F_2}\left( {t,{y_1},{y_2}} \right) = \frac{{T(t)}}{m} - \frac{B}{m}{y_2} - \frac{C}{m}{y_1}。\\ \end{gathered} \right. $

(4)

使用上述数值方法的标准解，可以离散获得y₁和y₂。此外，使用上述任何数值方法的规范公式，可以在离散时刻获得y₁和y₂的值。当然，在计算过程开始之前，必须指定积分步骤和初始条件，例如t=t₀=0，y₁(t₀)=y₁₀=x₀，y₂(t₀)=y₂₀=$ {\dot x_0} $。

将降序过程扩展到原始方程组(2)，得到一阶微分方程组，该方程组的解可以用标准Runge-Kutt算法实现。为此，定义：

$\begin{split} &y_{1}=x_{1}\text{，}y_{2}= {\dot x_1} \text{，}y_{3}=x_{2}\text{，}y_{4}= {\dot x_2} \text{，} \cdots \text{，}\\ &y_{m-1}=x_{n+1}\text{，} y_{m}= {\dot x_{n + 1}}。\\[-1pt] \end{split}$

(5)

这些解可以在给定的初始条件下通过Runge-Kutt递归公式得到，则计算过程的实现将基于以下递归依赖关系：

$ {y_{j,i + 1}} = {y_{j,i}} + \frac{1}{6}\left( {{K_{1,j}} + 2{K_{2,j}} + 3{K_{3,j}} + {K_{4,j}}} \right)。$

(6)

式中：j为等式的数字，j=1，2，3，…，m；m=2(n+1)；i为积分步数，i=0，1，2，3，…；K_1,j、K_2,j、K_3,j、K_4,j系数根据公式计算：

$ {\left\{\begin{split}&{K}_{1,j}=h{F}_{j}({t}_{i},{y}_{1,i},{y}_{2,i},\cdots ,{y}_{m,i}), \\ &{K}_{2,j}=h{F}_{j}({t}_{i}+{h}/2,{y}_{1,i}+{K}_{1,1}/2, {y}_{2,i}+{K}_{1,2}/2,\cdots ,{y}_{m,i}+{K}_{1,m}/2),\\ &{K}_{3,j}=h{F}_{j}({t}_{i}+{h}/2,{y}_{1,i}+{K}_{2,1}/2, {y}_{2,i}+{K}_{2,2}/2,\cdots ,{y}_{m,i}+{K}_{2,m}/2),\\ &{K}_{4,j}=h{F}_{j}({t}_{i}+{h},{y}_{1,i}+{K}_{3,1}, {y}_{2,i}+{K}_{3,2},\cdots ,{y}_{m,i}+{K}_{3,m})。\end{split} \right.}$

(7)

1.2 曲轴刚度和传递函数系数

有限元法的存在使得可以通过更精确的简化模型确定曲轴刚度和传递函数系数。利用Workbench建立曲轴中段有限元模型，将曲轴中段右侧框架颈部端面（A节点）固定约束，在左框架颈端面截面中心（C点），施加轴向作用力P，在给定的边界条件和轴/径向力作用下，曲轴C点轴向变形为Δl_p；根据刚度定义，曲轴轴向刚度e₀=P/Δl_p。

船舶推进轴系轴向振动很大程度上是由柴油机曲柄在径向载荷作用下的变形引起的。通过对曲柄连杆机构的结构动力学进行分析计算，在曲柄中间段的颈部中间施加径向力Z，在这种力的作用下，曲轴弯曲，弯曲变形伴随着曲轴的纵向缩短或延长，等效轴力T=Δl_x/e₀，传递函数$ {k_e} $=T/Z。

1.3 推力轴承刚度和阻尼

柴油机内部的推力轴承是在液体摩擦模式下运行的滑动轴承。众所周知，在这种模式下，摩擦表面被一层油隔开，它们之间没有金属接触。在这种情况下，螺旋桨振动通过油膜传递给船体，因此轴向振动的发展模式在很大程度上取决于油层流体动力学定律。

雷诺兹方程是确定流体动态压力及其间隙变化的理论基础。如果油膜厚度较低，且油粘度恒定与柴油机进口温度相对应，则此圆柱形坐标方程为^[13]:

$ \begin{split} & \frac{\partial }{{\partial r}}\left( {r{h^3}\frac{{\partial p}}{{\partial r}}} \right) + \frac{\partial }{{\partial \phi }}\left( {{h^3}\frac{{\partial p}}{{r\partial r}}} \right) = 12\mu {\vartheta _z}r + \\ &6\mu \frac{\partial }{{\partial r}}\left( {rh{\vartheta _r}} \right) + 6\mu \frac{\partial }{{\partial \phi }}\left( {h{\vartheta _\phi }} \right) 。\end{split} $

(8)

式中：p为轴瓦表面压力分布；$ \mu $为油的动态粘度；h为油膜厚度；$ {\vartheta _r} $、$ {\vartheta _\phi } $、$ {\vartheta _z} $分别为轴瓦圆周径向、切向和轴向点速度的分量（见图2）。

对于倾斜工作面的流体动力自定位止推轴承，轴向间隙取决于半径r和角度$ \varphi $（见图3）。在这个基础上，厚度h的依赖关系来自纯粹的几何推理：

$ h\left(r,\varphi\right)=h_0+r\mathrm{sin}\varphi\mathrm{tan}\gamma=h_0+r\gamma\mathrm{sin}\varphi。$

(9)

式中：h₀为安装间隙；$ \gamma $为轴瓦的倾斜角度，$ \tan\gamma = \gamma $。

式(8)解决了确定一个推力轴承垫块油层流体动力压力曲线的问题。通过集成压力表达式，可以根据公式计算所需的推力轴承刚度和阻尼系数：

$ \begin{gathered} K = z\int_0^\theta {\int_{{R_1}}^{{R_2}} {\left( {\frac{{\partial p}}{{\partial h}}} \right)} } r{\rm{d}}r{\rm{d}}\phi，\\ B = z\int_0^\theta {\int_{{R_1}}^{{R_2}} {\left( {\frac{{\partial p}}{{\partial {\vartheta _z}}}} \right)} } r{\rm{d}}r{\rm{d}}\phi 。\\ \end{gathered} $

(10)

式中：z为推力轴承轴瓦的数量。

为求解式(10)，将Reynolds方程转换为有限差分形式，就得到了一个线性代数方程组，用最小二乘法近似得到的积分解。

2 轴系振动随机因素模拟 2.1 柴油机随机激励模拟

柴油机是船舶轴向振动的主要来源。从轴向振动激发的角度来看，柴油机是一个由热力学、气体动力学和机械运动组成的复杂系统。柴油机径向力的表达式^[14]为：

$ Z=zF=pF\frac{\mathrm{cos}\left(\alpha+\beta\right)}{\mathrm{cos}\beta}=pFw 。$

(11)

式中：p为驱动力；w为径向力函数，$ w = \cos\alpha \left( {1 + \lambda \cos\alpha } \right) - \lambda $；F为活塞面积；$ \lambda $为常数，曲柄半径R与连杆长度L之比。驱动力p是作用于活塞的气体压力、运动质量的惯性和这些质量的重力等力的总和，表达式为：

$ \left\{\begin{aligned} & p\left(\alpha\right)=p_e-p_x+q_s+g_s，\\ & q_s=-m_sJ\approx-m_sR\omega^2\left(\mathrm{cos}\alpha+\lambda\mathrm{cos}2\alpha\right) ，\\ &g_s=m_sg。\end{aligned} \right.$

(12)

式中：$ {p_e} $为燃烧气体对活塞的压力；$ {p_x} $为活塞腔中的压力；$ {q_s} $为运动质量的惯性力；$ {g_s} $为重力；$ {q_s} $为与活塞面积相关的运动质量；$ \omega $为曲轴旋转角速度；g为重力加速度。

可以看出，柴油机径向力的随机性质是由气体压力、体积和质量决定的，以上因素均随机，其中气体压力的随机性是由制造和磨损引起的燃油喷嘴不同、活塞气腔存在不同程度的间隙以及燃油泄漏共同造成燃料循环的不均匀性，从而导致气体压力均有随机性；另外，体积和质量以及前述的零件间间隙均与零件尺寸直接相关，由于几何结构与设计尺寸存在随机偏差，因此体积和质量即使在公差范围内仍具有一定的随机性。

因此，为研究柴油机径向力的随机性，必须对每个参数进行统计处理，确定每个参数的分布规律。经文献调研分析：

1）文献[15]对曲柄连杆机构中气缸、曲柄、连杆、活塞等主要零件尺寸和质量进行了概率统计分析，以上参数均满足正态分布，且其变异系数为0.002；

2）文献[16 - 18]，压缩过程结束时压力的变化服从正态分布，且其变异系数为0.155；

3）文献[16 - 18]，船用柴油机额定运行模式下燃油供应服从正态分布，且变异系数为0.06～0.08；

4）文献[16 - 18]，平均压力符合正态分布，变异系数为0.148。

从本质上讲，以上文献数据已经证明了船舶轴向振动具有偶然性，其偶然性是由径向力的随机性决定的。径向力与算术平均值的随机偏差由变异系数决定，变异系数不仅反映了样本的相对散射量，而且还反映了样本的均匀性。如果变异系数V<0.33则变量序列被认为均匀，结合概率统计理论，该统计数据服从于正态分布定律，即对于均匀统计数据，随机因素的作用是独立的和等价的，并不存在一个因素占主导地位的情况。实际径向力与平均值的偏差遵循正态分布定律，其数值随曲轴转角变化而不同，实际径向力的变异系数为0.309～0.026。

利用蒙特卡洛方法模拟柴油发动机在最大持续功率模式下运行的模拟结果，经过150次的模拟，得出如图4所示随机数据源下的相对径向力，随机性的最大影响发生在燃料燃烧时刻，其他时间的随机性比例微不足道。

同时，为充分说明径向力的随机性符合正态分布，对径向力分布经验密度和理论密度进行对比，径向力正态分布的假设得到了充分的证实，如图5所示。

2.2 螺旋桨随机激励模拟

由于螺旋桨空间构造复杂，实际制造过程中难免存在尺寸误差，因此，螺旋桨在如桨叶厚度、倾角等方面存在一定的随机误差。根据文献[19]，螺旋桨轴直径、盘面比、螺距、螺旋桨直径以及桨叶倾角均服从正态分布，各参数的变异系数分别为0.0012、0.0023、0.0024、0.0018、0.4。一般来说，由于螺旋桨在不均匀的水流中工作而产生的可变力和力矩的大小取决于许多因素，其中包括船尾末端的形状、螺旋桨在船尾视图中的位置、轴线的倾角、叶片形状、吃水、船舷和垂直摆动。另外，螺旋桨进气口的斜面或螺旋桨轴与进气口速度的倾斜也会导致螺旋桨叶片的非稳定周期性流动条件，即使进气口的速度场均匀，因此，在斜流中工作的单个螺旋桨叶片上的推力和力矩在一次旋转中周期性地变化。

根据螺旋桨推力和扭矩影响因素的分析，并结合推力系数和扭矩系数近似公式，推导得出螺旋桨随机推力和扭矩可表示为：

$ P_{gw}=\frac{2P_{gwm}\left(1-\displaystyle\frac{V_{cyo}}{\omega R_0}\mathrm{tan}\phi\cos z\alpha\right)^2}{2+\left(\displaystyle\frac{V_{cyo}\mathrm{tan}\phi}{\omega R_0}\right)^2}，$

(13)

$ M_{gw}=\frac{2M_{gwm}\left(1-\displaystyle\frac{V_{cyo}}{\omega R_0}\mathrm{tan}\phi\cos z\alpha\right)^2}{2+\left(\displaystyle\frac{V_{cyo}\mathrm{tan}\phi}{\omega R_0}\right)^2} 。$

(14)

式中：V_cyo为船舶航速，kn；$\phi $为桨叶倾角，＜20°；α为螺旋桨转角；R₀为桨叶半径；P_gw、M_gw分别为推力和扭矩的随机值，P_gwm、M_gwm分别为推力和扭矩的平均值。

式(13)、式(14)可根据螺旋桨几何参数及其流体动力脉动的随机偏差，确定每艘船的推力和扭矩的周期性。图6为螺旋桨推力变化图的计算结果，其中包含蒙特卡洛方法生成的参数。获得了螺旋桨推力的统计数据。图7显示了时间为0.1 s的单个截面对给定图形依赖关系的统计处理结果。研究表明，螺旋桨作为触发船轴随机摆动的源，遵循正态分布定律，变异系数在0.067～0.01。

3 考虑随机因素的轴系振动计算 3.1 计算模型及参数

本文研究对象船用动力装置是一个双轴驱动，扭矩直接传递到螺旋桨上。每台发动机为6NVD48A2U发动机，额定功率485 kW，320 r/min，配有硅减震器和内置推力轴承、轴线和四叶片固定节距螺旋桨。柴油机推进离散轴系示意图如图8所示。其中，1是硅阻尼器A-710，惯性矩为15.73 km²；2～7为轴颈直径为215 mm的曲轴；8是推力轴承；9是惯性矩302.57 kgm²，质量1889.5 kg的飞轮；10是前轴、中轴和尾轴，直径D分别为155、170以及180 mm，其中推进器轴锥区域直径为200 mm；11是螺旋桨，质量786.02 kg，直径1550 mm，圆盘比0.55。所有质量都由相应的轴向柔性连接。该模型的参数见表1。

根据表2提供的柴油机推力轴承参数，利用式(10)计算得出柴油机推力轴承的轴向刚度为1.632×10⁹ N/m和阻尼系数3.817×10⁶ kg/s。同时，利用Ansys有限元仿真，计算得出柴油机曲轴轴向刚度为1.91×10⁹ N/m，曲轴轴向传递系数为0.124，有限元曲轴轴向和径向变形计算结果如图9所示。

3.2 计算结果分析

自由振荡的计算表明，在研究范围内0～200 Hz有3种自由振荡频率：1）单节点振荡形式f₁ =52.84 Hz；2）双节点振荡形式f₂ = 74.10 Hz；3）三节振荡形式f₃ = 156.14 Hz。共振计算表明，在运行转速范围内存在11阶～14阶轴向共振，其转速分别为288、276、264、254、244、235、226 r/min。

以转速为264 r/min的共振振荡为例，同时为评估考虑随机因子对船舶轴向振动的影响，在此分别采用能量法和考虑随机激励的动力学方程求解（即1.1中描述的Runge-Kutt求解方法），2种计算方法在以上共振转速下的轴向振幅如表3所示，2种方法求解得出的振幅差异非常明显。例如，对于共振转速264 r/min，这种差异为30.07%，这种差异在图10所示的曲轴前端振动曲线得到了证实。可以看出，能量法只考虑一个强迫力谐波，而本文提出的考虑随机因子的Runge-Kutt求解方法不仅考虑了螺旋桨的随机激励还考虑柴油机激励的随机性。

此外，利用随机数据生成柴油机工作过程参数的能力，模拟了所考虑的振荡过程，其持续时间等于相关间隔，该过程曲轴前端随机轴向振动如图11所示。

通过201次模拟统计试验，得出了所有共振模式的平均振幅和变异系数，如表4所示。

分析可知：

1）轴振动模拟中考虑随机因素的效果是存在的，随着轴线转速的增加而增加；

2）考虑随机性可以减少振幅，例如，转速264 r/min，共振振幅从0.094 mm降低到0.093 mm，降低1.0%，相对平均值的变异系数为0.143；

3）在船舶运行期间，轴向振幅可达到0.108 mm，小于船级社允许的0.15 mm。

4 结　语

1）建立了一种基于概率方法计算轴系轴向随机激励、并利用Runge-Kutt数值方法求解对二阶微分方程的轴系轴向振动计算方法；

2）实现了利用有限元法确定曲轴刚度特性，并利用油层流体动力学计算推力轴承刚度和阻尼参数的方法；

3）基于统计数据结合蒙特卡罗方法实现了对螺旋桨推力和柴油机径向力概率特性的模拟。