用于推进轴系振动分析的改进数值组装法

2022-08-16 09:50谢溪凌张志谊

振动与冲击 2022年15期

巫頔，谢溪凌，张志谊, 2

(1. 上海交通大学机械系统与振动国家重点实验室，上海 200240； 2. 上海交通大学振动、冲击、噪声研究所，上海 200240)

船舶推进轴系振动校核是一项重要的设计内容，高精度的动力学模型和计算方法是获得较精确的振动特性和振动响应的基础。有限元法是目前广泛使用的数值方法，但对于诸如梁的一维结构，可以通过振动微分方程在梁的自然分段内精确描述动力响应，得到精度更高的计算结果。目前，已有多种变截面梁横向振动精确求解方法，如传递矩阵法、谱元法、数值组装法等，其中传递矩阵法在船舶轴系振动分析中应用最为广泛。传递矩阵法中一般取梁单元端点挠度、转角、弯矩与剪力作为梁单元的状态变量，通过梁单元内部的传递关系和节点连续性条件，将各单元传递矩阵按序累乘[1]。传递矩阵法缩减了系统方程阶数，计算速度快，但是系统矩阵的病态易导致特征方程的数值不稳定，且梁单元内部振动信息无法获得。Wu等[2]从界面连续性条件出发，提出连续质量传递矩阵法(continuous-mass transfer matrix method, CTMM)，并用于求解多跨变截面Timoshenko梁横向振动特征值问题。CTMM方法采用梁单元振动微分方程的4个系数作为梁单元的状态变量，便于连续表示梁单元内部振动。Zhang等[3]在CTMM中考虑了梁的横向与扭转弹簧约束、附加弹簧-质量系统、弹簧-质量-弹簧约束以及联轴器连接等各种常见边界条件。徐颖蕾[4]推导了CTMM的模态正交性，并使用模态叠加法求解轴系的强迫振动响应。Xie等[5]在界面连续性条件中加入了外载荷向量，并在频域直接求解变截面Timoshenko梁的强迫振动响应，避免了模态叠加过程中的截断误差。CTMM方法融合了传统传递矩阵法的优势，但依然存在高频数值发散问题，影响特征频率与强迫振动响应的求解。

为避免固有频率求解中的数值发散问题，Riccati传递矩阵法[6]将经典传递矩阵法中两端状态变量的对应关系转换为在每一截面上位移与载荷的对应关系，减小了矩阵累乘误差，提高了系统行列式的数值稳定性，但奇点干扰和多次求逆带来的数值误差仍会影响固有频率的求解精度。王正[7]对Riccati传递矩阵法在求特征值时遇到的奇点问题进行了讨论，并提出了消除奇点干扰后的频率方程式；黄娟等[8]在计算汽轮发电机轴系时使用消去奇点干扰的改进Riccati方法。Bestle等[9]为提高Euler梁横向振动的计算频率范围，提出了拆分梁段与变换局部坐标系方向等方法，但对高频数值发散的抑制效果有限，无法显著提高求解频率范围。

Wu等[10]提出的数值组装法(numerical assembly method，NAM)以梁单元微分方程的4个系数作为状态变量，将n个梁单元对应的4n个方程一次性求解。由于矩阵方程只求解一次，避免了传递矩阵法多次求解矩阵方程带来的误差放大问题，因此计算频率范围相较之有一定提高。Wu等解决了梁上分布附加弹簧-质量块的均匀Euler-Bernoulli梁特征值问题。Lin等[11]利用NAM计算了附加集中质量、集中转动惯量、弹簧与扭簧等多种集中元素的多跨Timoshenko梁的固有频率和模态振型。Yesilce等[12]使用NAM分析附加多个弹簧与质量块的多跨变截面梁，考虑了轴向作用力对系统横向振动固有频率的影响。Xu等[13]使用NAM计算Euler-Bernoulli阶梯梁，通过局部坐标系拓宽计算频带。NAM避免了误差的累乘放大，因而能够提高计算频率范围。但对于梁段单元参数差异过大的系统，通过NAM法建立的系统矩阵仍存在病态严重的问题，制约了轴系高频段振动的求解。

本文首先分析了现有计算方法的误差来源，然后使用局部坐标系推导多跨Timoshenko阶梯梁的横向振动方程，并考虑推进轴系上的质量、惯量、弹簧及联轴器等集中元素，进而给出行归一化措施，降低矩阵方程的条件数，提高推进轴系高频振动计算精度。

1 基于Timoshenko梁理论的改进数值组装法

1.1 Timoshenko梁横向振动方程

考虑转动惯量和剪切变形对横向振动的影响，设第i梁段横向振动位移为yi(x,t)，转角为φi(x,t)，对应的运动方程可表示为

(1)

(2)

式中：ρi为材料密度；Ii为截面转动惯量；Ei为弹性模量；Gi为剪切模量；Ai为截面面积；κi为Timoshenko梁剪切系数。设yi(x,t)=Yi(x)ejωt，φi(x,t)=Φi(x)ejωt，将转角φi与挠度yi解耦可得

(3)

(4)

式(4)为4阶齐次微分方程，挠度Yi(x)的解为振型函数，如式(5)所示

Yi(x)=C1icoshλi(x-Xi-1)+C2isinhλi(x-Xi-1)+

(5)

类似的，转角的振型函数Φi(x)为

Φi(x)=C1iqisinhλi(x-Xi-1)+C2iqicoshλi(x-

(6)

其中

1.2 界面连续性条件

在轴系横向振动时，梁段节点满足挠度、转角、弯矩与剪力的协调条件，如式(7)

(7)

式中：KRi为梁段交界处转角弹簧刚度；KTi为线弹簧刚度。当第i节点有集中质量mi和集中转动惯量Ji时，KRi=-ω2Ji，KTi=-ω2mi。在局部坐标系下，第i梁段的坐标原点定义为Xi-1,R，Xi,L为第i梁段的终点。将Xi,R=0，Xi,L=li代入式(6)，并将挠度、转角、弯矩与剪力的表达式展开，可得矩阵方程

(8)

其中

方程式(8)可简写为AcriC(i)=Acli+1C(i+1)，Acri即为第i梁段右端参数矩阵，Acli+1为第i+1梁段左端参数矩阵，C(i)为第i梁段微分方程解的系数向量。若第i节点有联轴器，则节点位移与转角不连续，式(7)中的转角与位移协调条件变为

Yi(Xi,L)-Yi+1(Xi,R)=

(9)

(10)

式中：KITi为联轴器横向刚度；KIRi为转角刚度。相应的，矩阵方程式(8)的形式变为

(11)

1.3 边界条件

对于n段变截面梁，共有n+1个节点，其中有n-1个内节点，每个节点存在4个连续性条件构成的界面连续性条件方程组；两个边界节点，各有两个边界条件方程。不同边界条件对应的方程组形式不同，在无位移约束，两端有集中质量、集中转动惯量或弹簧支承的条件下，梁段边界的剪力与弯矩平衡条件表示为

左端

(12)

右端

(13)

式中：KTL，KTR分别为左右两端的横向弹簧；KRL，KRR为左右两端的扭转弹簧；ML，MR，JL，JR为各端的集中质量与集中转动惯量。利用式(12)、式(13)的边界条件方程组可以得到边界节点的内力表达式

BLC1=0,BRCn=0

(14)

式中：BL与BR为2×4矩阵，为自由-自由条件下边界节点处的剪力和弯矩为零。矩阵BL与BR见附录A所示。在无外力作用的条件下，将各梁段方程依次组装并求解，即为NAM的矩阵方程

(15)

记方程式(15)为T(ω)C(ω)=0，通过求解T(ω)的特征方程T(ω)=0可求得系统固有振动特性；求解强迫振动响应问题时，外力向量在方程右侧，矩阵方程可简写为T(ω)C(ω)=F。

1.4 行归一化加权阵

随着频率升高，NAM的特征方程(T(ω))曲线会在高频出现数值发散现象。观察方程式(14)可以发现，矩阵T(ω)的每一行与系数C的乘积都有明显的物理意义，各子块的四行与系数向量的乘积分别对应梁段的挠度、转角、弯矩与剪力。随着频率的升高，挠度与转角所对应的行元素会与剪力、弯矩所对应的行元素存在数量级的差距，导致矩阵条件数升高，影响方程对舍入误差的敏感度。当矩阵病态严重，系统特征方程曲线就会发散。因此，采用适当的对角加权阵W(ω)左乘T(ω)，以改善矩阵的病态问题。加权矩阵W(ω)可取为

(16)

式中，m为矩阵T的阶数。对T(ω)左乘加权矩阵，可平衡各行的数量级差距，有效降低矩阵的条件数，减弱方程式(8)对数值误差的敏感度，在保证求解数值稳定性的同时，提高了计算频率范围。因此，本方法实际上为改进的数值组装法(modified NAM, m-NAM)。左乘W(ω)后，特征方程式变为

Δ=W(ω)T(ω)=0

(17)

为避免数值过大对频率方程式求根带来的不便，参考Wu等的研究中对CTMM系统特征方程的变换，对Δ做如下变换

(18)

对变换后的特征方程δ(ω)求根，即可求得系统固有频率。

2 数值算例验证

2.1 算例1(等截面解析算例)

对于均质等截面Timoshenko梁，在自由-自由边界条件下，系统矩阵行列式的展开式为

(19)

式(19)即为系统特征方程，对其求根可得系统固有频率的解析解。以一段长度10 m、直径0.06 m的圆截面钢梁为例，分别使用CTMM、NAM、m-NAM与解析表达式计算系统的特征方程曲线。将梁分为3段，长度分别为1 m，8 m，1 m，计算结果如图1所示。

图1 系统特征方程曲线

CTMM也是针对梁段微分方程的系数进行求解，将式(8)改写为

(20)

通过递推式，得到第一段系数C1与第n段系数Cn的关系，并利用边界条件构建系统矩阵

(21)

方程式(21)即为CTMM法的系统特征方程，其中包含各梁段子块累乘的部分(第三、第四行)。在数值仿真中发现，分段数与分段长度对CTMM系统矩阵数值大小没有影响，而各段的舍入误差不断累积并放大，且系统矩阵Z病态严重，方程的解易在微小数值扰动下发散。

NAM中梁段的误差没有经过累乘放大，且由于矩阵病态主要来源于矩阵元素的数量级差距。在梁段子块中，随着频率升高，双曲余弦、双曲正弦元素与正弦、余弦元素的数值(一、二列与三、四列)量级差距非常大，在运算过程中舍入误差会对矩阵方程的求解造成扰动。从图1中可以看出，CTMM法的特征方程曲线在200 Hz以上完全失去信息，无法判断固有频率，NAM的计算频率范围也止于350 Hz，而m-NAM能够很好地保持数值稳定性。从图2中的对比可以看出，m-NAM的前150阶固有频率与对式(19)求根得到的固有频率解析解吻合良好，固有频率的相对误差在0.5/10 000以内。

图2 模态频率计算结果及相对误差

2.2 算例2(推进轴系算例)

为进一步验证m-NAM在轴系动力学特性计算上的准确性，使用有限元方法建立算例轴系模型。轴系总长15 m，采用实心圆截面，材料为钢，杨氏模量E=2.1×1011Pa，密度ρ=7 850 kg/m3，泊松比ν=0.3。

表1 算例系统轴段长度与截面参数

轴系共有4个支承，分别轴承k1，轴承k2，轴承k3以及轴承k4。轴系由支承位置与截面变化共分为11段，共12个节点，4个支承布置于节点2，5，8，11处，如图3所示。k1=5×108N/m，k2=k3=k4=1×109N/m，各支承置于刚性基础上。各轴段对应长度与截面半径如表1所示。考虑到较高的计算频率范围，设置单元尺寸10 mm，共计15 000个单元。