水下爆炸冲击波载荷沿燃气轮机结构传递特征

2022-10-11姚熊亮熊帮虎王志凯杨娜娜张文启万泽臣

兵工学报 2022年9期

姚熊亮，熊帮虎，王志凯，杨娜娜，张文启，万泽臣

(哈尔滨工程大学船舶工程学院，黑龙江哈尔滨 150001)

0 引言

舰船在服役期间，会受到严酷的海况以及各种水下武器如水雷、鱼雷等的攻击，这些水下武器的攻击会对舰船弯曲外板、龙骨翼板及舱室设备等造成严重损伤。根据水下武器与舰船距离的不同，可将舰船水下爆炸冲击分为两大类：接触爆炸和非接触爆炸。水下接触爆炸会造成舰船局部破损，导致舰船失衡；在水下非接触爆炸的情况下，会对舰船动力系统和设备等造成损伤，导致舰船失去动力。舰用燃气轮机作为舰船的心脏，不仅要拥有优良的动力性能，还需要拥有较强的抗冲击性能。目前对于燃气轮机的抗冲击研究主要集中在结构设计，对燃气轮机的强度设计、爆炸冲击波载荷在燃气轮机内部的传递相关研究较少。为得到水下爆炸冲击载荷作用于舰船或设备的传递特性，谢耀国等基于小波分析，对水下爆炸自由场压力载荷和壁压载荷进行时频域特征分析，得到载荷所包含的频域信号并统计了载荷在每个频段的能量。李丽萍等采用小波包对爆炸原始信号进行分析，并建立起冲击波能量与不同类型炸药、爆距与毁伤目标固有频率之间的关系。

目前得到信号能量的方法多为快速傅里叶变换、小波分析等。快速傅里叶变换多为分析信号为周期性的平稳信号，小波分析在处理非平稳的信号时具有很强的优势而且具有多辨率特性，但选择小波基有异，处理的结果会也不同，而且有限长的小波基也会造成能量的泄露。爆炸冲击载荷具有突变快、持续时间短、含多种噪声干扰的强非线性信号。因此，以上方法不再适合爆炸载荷能量的提取。

本文用经验模态分解(EMD)方法和与之对应的希尔伯特- 黄变换(HHT)，对爆炸冲击波载荷在燃气轮机中的传递特性进行分析。Jin等利用希尔伯特- 黄变换法对节理梁的非线性特性进行研究，并通过对比加速度信号和位移信号验证了数据的正确性。周鹏运用经验模态分解法对比分析了水下爆炸冲击波载荷在单体船和三体船在型宽及型深方向上的分布规律和传递特性，最后得出水下爆炸冲击波信号中高频和低频与舰船冲击响应的关系，表明EMD方法可以分析载荷分布规律和传递路径。相较于舰船研究，现如今对舰船设备的载荷传递特性研究较少，故本文拟采用经验模态分解法对爆炸冲击载荷在燃气轮机中的传递特性分析，并对燃气轮机进行力学模型简化，以验证传递特性的正确性。本文研究内容对燃气轮机抗冲击设计具有参考价值。

1 燃气轮机结构及水下爆炸特征

1.1 燃气轮机结构特征

1.1.1 燃气轮机关键部件类型确定

随着我国海军的迅速发展，高性能动力装备已显得愈发重要，燃气轮机在舰船中应用可提高经济性，在一定的燃料供应下可以提高船舶的航程，其燃气轮机在舰船上的应用已经成为现代化舰船的重要标志。

燃气轮机基座是燃气轮机强有力的支承构件，它的抗冲击特性是国内外科研人员研究的重点。根据韩璐等的试验结果，基座与机匣连接处位置在受到冲击后易发生断裂，是抗冲击的薄弱部位。机匣作为燃气轮机与舰船相连接的主要承力构件，承受由燃气轮机工作时产生的各种气体载荷、冲击载荷、涡动载荷以及惯性载荷等的作用，故机匣不仅仅要满足支承功能，其刚度、强度和可靠性以及稳定性也需要满足需要。转子及轴在外载荷和内载荷双重的冲击下，极易发生径向跳动，产生弯曲现象。这些情况的发生极易使燃气轮机发生故障进而使舰船丧失动力系统。而在仿真计算时，由于传动轴系统零部件多且结构复杂，在建模、设置边界条件及接触的过程中，会面临计算复杂与工作量大的情况。综合上述分析，本文明确燃气轮机关键部件，主要包括：燃气轮机基座、燃气轮机机匣及转子轴系。根据以上关键部件，建立准确可靠的数值模型对于计算燃气轮机在水下爆炸冲击载荷的作用下的动力响应并基于此研究载荷传递路径及系统动力学特性而言是非常关键的。

1.1.2 关键部件力学等效方法

燃气轮机作为一件复杂的船用设备，其内部由大量部件构成，而各个部件间存在复杂的连接方式，同时构件有复杂的曲面，还存在减振器等具有非线性因素的构件，按实际结构建立燃气轮机模型非常困难，即使建立完整的燃机模型，对燃气轮机进行有限元网格划分的质量也不理想。在进行计算时，这样低质量的有限元网格会导致计算难以收敛，且对不需要关注的部件也进行精密计算，造成计算资源的浪费。故需要对燃气轮机实际结构简化，考虑燃气轮机重要部件，对不需要关注的部件可以作为质量点加载在有限元模型上，根据实际结构的力学特点对燃气轮机做适当的简化，并建立易于有限元网格划分及计算的几何模型。

简化几何模型应遵循两条基本原则：1)保留模型的主要尺寸与基本特征；2)能够使得划分的网格质量有明显改善。模型简化后只可以近似反映工程的实际情况，模型正确是有限元计算结果正确的关键，故要充分考虑到结构的实际情况。此外，在建立有限元模型时，在满足一定精度的同时不要划分过多的单元，减少计算时间，即要具有良好的经济性能。

根据燃气轮机转子系统的三维模型可知，转子叶片、轮穀和部分轴段需要简化。运用转动惯量互等的原理，将叶片简化成圆柱，空心圆柱转动惯量参数表示如图1所示，具体计算公式如下：

图1 空心圆柱转动惯量参数表示Fig.1 Parameter representation of moment of inertia of hollow cylinder

(1)

式中:为轮盘叶片数；为轮盘单个叶片质量；为轮盘叶片长度；为叶片重心到转轴中心距离；为设计圆盘外环总质量；为设计圆盘外环外径；为设计圆盘外环内径。

滚动轴承通过鼠笼弹性支承与轴承座弹性连接，滚动轴承- 鼠笼组合支承如图2所示，带有轴的滚动轴承固定在鼠笼内壁上，鼠笼固定在轴承座上，其中与分别为滚动轴承刚度和鼠笼刚度，为滚动轴承和鼠笼的等效质量，为弹性支承相应的刚度，为轴承和鼠笼等效轴承的位移，为轴径的位移。

图2 滚动轴承- 鼠笼组合支承Fig.2 Combined bearing of rolling bearing and squirrel cage

由图2(b)简化为图2(c)的弹性支承模型，其等效质量运动微分方程可以表示为

(2)

设方程的解为

=ei

(3)

式中：为运动幅值；i为复数表示的虚数；为转子涡动频率。则

(4)

建立了轴承座中心坐标与轴径中心之间的关系后，与图2(c)联立，弹性支承相应的刚度可以求得，即

(5)

不是一个常数，与转子的涡动角速度有关。滚动轴承刚度在10N/m级别，鼠笼刚度在10N/m级别，质量较小，在低速转动时，总刚度的变化不大，因此可以把滚动轴承和鼠笼的刚度等效为各部件串联后的刚度，串联后的刚度系数为

(6)

根据文献[19]，可选取滚动轴承刚度为9882 1×10N/m，鼠笼支承刚度为1382 1×10N/m，将刚度值代入刚度公式，可得鼠笼弹性支承刚度为1211 8×10N/m。

1.2 水下爆炸载荷冲击输入

当舰船遭受到水下爆炸载荷冲击作用时，设备基础会受到加速度冲击。在对设备进行瞬态冲击响应分析时，通常采用实测的加速度时历信号或者标准加速度时历曲线作为设备的冲击输入。本文燃气轮机使用的冲击加速度信号是根据德国军用标准 BV 043/85，使用组合三角波的形式对其基础进行加载。根据设备安装位置不同，选择相应的抗冲击标准。由于设备的冲击环境极其复杂，设备的冲击载荷会以响应谱的形式表达。如果把冲击信号以频域表示，则如图3所示设计谱，左右拐点分别为左频率和右频率，低于左频率则为低频响应，此时以激起舰船及设备位移为主，在左右频率中间则以激起舰船及设备速度为主，高于右频率则以激起舰船及设备加速度为主。典型冲击响应谱如图3所示，典型三折线谱如图4所示，低频段的直线代表等位移段，用谱位移来表示，中频段的直线为等速度段，用谱速度来表示，高频段的直线代表等加速度段，用谱加速度来表示，三条直线相交，便形成了三折线设计谱。

图3 典型冲击谱曲线Fig.3 Typical impact spectrum curve

图4 典型设计冲击谱Fig.4 Impact spectrum of typical design

根据规范规定及设备的安装情况，可查得设备的谱位移值、谱速度值和谱加速度值。当舰载设备的质量大于5 t时，需对参数进行折减计算：

(7)

式中：为隔离安装的设备质量(t)；为质量常数，恒等于5 t；表示折减后的谱加速度值(g)；表示折减后的谱速度值(m/s)。

在使用时域分析方法对燃气轮机进行试验时，在试验中测得的时间加速度曲线具有很大的随机性，对结果的分析也会产生很大的影响。故选用德国军用标准BV 04385规范将规定谱值转化为双三角波时间历程曲线，如图5所示。图5中，轴表示加速度(m/s)、轴表示时间(s)，、分别表示与时刻的加速度，、分别表示正三角波积分面积和负三角波积分面积，为正三角波最大值时间、为正负三角波交界时间点、为负三角波取值最大时间点，为负三角波作用时间结束点。

图5 正负三角波历程曲线Fig.5 Positive and negative triangul wave history curves

其转化关系如下：

1)第1个正三角形的加速度峰值，约为最大加速度的06倍；

2)第1个三角形的面积约为最大速度的075倍；

3)第2个负三角形面积应与第1个三角形面积大小相等，从而使基础的最终速度为0 m/s；

4)对加速度历程两次积分、得到位移，此位移比设计谱的最大相对位移要大105倍；

5)再选择和，满足=04和-=06(-)。

根据以上关系，可得到组合三角波各参数的计算公式为

(8)

式中：为冲击谱谱位移。

通过(8)式将冲击谱的谱值转化为组合三角波时间历程曲线，得到的计算输入载荷如表1所示。

表1 组合三角波时历曲线Table 1 Time history curve of combined triangle wave

2 数值计算模型与冲击信号处理方法

2.1 数值计算模型

选取某型燃气轮机为研究对象，燃气轮机的结构支承座分别为前支承、后支承和中间定位支承，其有限元模型如图6所示。

图6 燃气轮机支承有限元模型Fig.6 Finite element model of gas turbine support

模型包括支承系统结构、燃气轮机机匣和转子系统。因为有些原因，只给出支承系统结构。

因为冲击波载荷在箱装体下部近似为均匀分布，故选取燃气轮机前基座处剖面位置进行载荷传递计算，在选择测点位置时避免焊接位置等应力集中点。各测点位置的分布具体如图7所示。使用相同的工况进行计算，载荷加载在箱装体正下方，载荷加载符合德国军用标准BV 085规定。

图7 测点分布位置示意图Fig.7 Distribution of measuring points

2.2 冲击信号处理方法

应力波在舰船及设备中传递引起结构的响应，设备中不同位置处的响应即可反映冲击波载荷的传递情况。因为设备自身的属性以及设备基座等的强非线性特征，设备的冲击响应信号也具有强非线性特征，传统的快速傅里叶变换对此等信号的处理就变得不那么有优势。故本文运用HHT对燃气轮机的加速度响应信号进行时频域分析。该方法主要包含完整集合经验模态叠加法(CEEMD法)以及希尔伯特变换。

EMD是Huang等提出的一种时频域分析方法，其优势在于处理和分析非平稳信号。该方法已被用到设备诊断、生物医学、结构分析等领域的研究中。EMD方法是根据经验把复杂的数据转换成有限个简单数据即本征函数(IMF)的叠加。对于时域信号()首先判别出其极大值、极小值，分别用三次样条曲线组成的包络线将其连接起来，确保样条曲线均经过这些极值点，取包络线均值点，则()与之差记为，即

()-=

(9)

若满足本征值模态函数，则记为1，否则需要进一步分解。重复进行上述步骤，对进行新一轮分解，以此类推，直至分解到满足判别式：

(10)

式中：表示对原函数的第次分解。记第1个本征函数1，将函数分解为若干个本征模态及单调的残差：

()=1+2+…++()

(11)

式中：为满足(10)式时的取值；()为残差。

对仿真数据()做EMD，其分解过程如图8所示，其中表示信号振动幅值。

图8 仿真信号EMD过程Fig.8 Simulating the signal EMD process

限于篇幅，下面以箱装体底板测点做EMD及能量分析：

1)提取时间加速度时域信号，并画出图9所示时间加速度曲线。

图9 时间- 加速度曲线Fig.9 Time-acceleration curve

2)对步骤1中图9所示的时间加速度信号做EMD，分解后的IMF各分量如图10所示，得到的时间加速度信号分解为10个IMF分量，并依次按顺序排列。其中1、2、3、4具有较高的频率，随着分解的慢慢递进，得到的IMF分量的频率逐渐减小，而波长呈现增长趋势。

图10 加速度信号的EMDFig.10 EMD of acceleration signal

分解出的本征函数大都具有清晰的物理意义。为了解载荷沿装备基座到设备内部关键部件的载荷传递情况，对各个本征函数进行能量分析，其计算公式为

(12)

式中：为第个IMF分所对应的能量；()为第个IMF分量；为各个IMF分量各个离散点的幅值；为各个IMF分量离散点的数量。对本征函数能量计算后得到的结果如图11所示，并且能量的分布随着频率的降低能量有变大的趋势。经过快速傅里叶变换得到每个函数的频域信号，可以得到5～10的中心频率分别为329 Hz、179 Hz、139 Hz、59 Hz、29 Hz、19 Hz。这与舰船设备在水下爆炸载荷尤其是冲击波载荷作用下呈现低频响应的理论相符合。10是分解后频率更小的分量，可能是信号固有的，也可能是其他情况下引起的。

图11 本征函数能量分布图Fig.11 Energy distribution of eigenfunction

3 计算结果与分析

3.1 载荷传递分析

在研究载荷沿燃气轮机由基座向内部转子时，分别选取前中基座各测点向涡轮内部传递。选取测点1、2、3、4、5，代表载荷向涡轮轴与转子轴载荷的传递特征。对各测点加速度数值进行经验模态分解，同时分析其各个IMF分量的主频率m及其能量占比，计算结果如表2所示。

表2 燃气轮机加速度测点IMF分量主频率及其能量占比Table 2 Main frequency and energy ratio of IMF component at acceleration measurement points on the gas turbine

由表2中对舰载燃气轮机沿前基座传递路线的结果中可以看出：加速度载荷在箱装体底部所激起的高频响应占总能量较高，约为总能量的0.98；随着载荷沿燃气轮机向内部转子的传递，受阻尼和材料固有属性的影响，载荷所激起高频响应逐渐减低；载荷从燃气轮机底部向内部传递的过程中，虽然总的能量在衰减，但是高频部分的占比逐渐降低，燃气轮机转子轴系相较于基座及箱装体底部的响应却很小，这是因为燃气轮机的前几阶模态都以低阶模态为主，爆炸冲击波载荷所引起设备振动的部分为低频部分，而爆炸冲击波载荷在燃气轮机的传递过程中主要以高频部分存在。

3.2 载荷沿燃气轮机方向传递特征

结合载荷在燃气轮机中的载荷传递特征，可以给出冲击波载荷在燃气轮机中传递路径，其示意图如图12所示。结合云图给出燃气轮机的传递路径。冲击波载荷以加速度的形式加载在箱装体下部的减振元件上，并沿箱装体内底板向机匣下部的每个基座传递，载荷到达燃气轮机机匣与基座连接处位置后还需要以弯曲波的形式沿机匣长度方向传递并沿轴承向转子传递，造成载荷在燃气轮机内部分布的不均匀。为验证上述载荷的传递规律，取燃气轮机前基座附近处的剖面，给出在同一工况下不同时刻的应力云图分布，如图13所示。

图12 燃气轮机载荷传递路径Fig.12 Gas turbine load transfer path

图13 燃气轮机典型剖面应力云图Fig.13 Typical stress cloud picture of gas turbine section

从图12和图13中可以发现：爆炸冲击波载荷作用在燃气轮机上的范围广，作用强度大；在=1 ms 时，冲击波首先作用在箱装体底板，造成箱装体底板的塑性应变，应力最大值达到345 MPa；随着爆炸冲击波载荷向燃气轮机内部的传递，载荷作用在基座支承及基座；在=2.5 ms时基座处出现较大的应力，应力最大处主要出现在前部两边基座的中间转接处位置，最大值为150 MPa，比箱装体底板小，未达到材料的屈服极限，也未出现明显的塑性变形，这是因为基座作为支承机匣和内部转子的强力构件出现了应力集中的现象，在冲击环境中的考核中需要着重注意此类强力构件；在3.75 ms时，机匣及转子开始出现应力，随着载荷的传递机匣和转子的应力开始变大，相较于基座处应力值出现了明显的衰减，但也未出现明显的塑性变形，这是因为基座的载荷传递到燃气轮机内部所引起的。

通过图13中的应力云图变化情况中可以发现，爆炸冲击波载荷在燃气轮机中的传递过程与理论结果一致，在传递过程中能量在衰减，应力响应也在衰减。总之，燃气轮机结构的变化导致冲击波载荷传递路径和分布的不一致，为燃气轮机在水下爆炸载荷冲击作用下的不同响应提供了条件。

3.3 燃气轮机力学模型分析与建立

当使用多自由度系统对燃气轮机进行力学分析时，可以将燃气轮机简化成如下3自由度力学模型，它包含基座、机匣、转子三个质量点和刚度为、、，阻尼为、、的3自由度系统。在多数情况下，阻尼对系统的影响可忽略不计，当冲击的激励力与系统的固有频率相近时，必须考虑阻尼的影响。

图14 燃气轮机力学模型Fig.14 Mechanical model of gas turbine

对如上3自由度系统进行力学分析并建立系统的运动微分方程，运用牛顿第二运动规律建立的方程为

(13)

(14)

在对设备进行抗冲击考核时，只有设备基座会遭受到冲击载荷，因此()=()=0 N。对系统的各参数及爆炸冲击载荷确认后即可对系统的响应、能量传递及载荷传递规律进行分析。其中=5 130 kg，=170 kg，=4 860 kg，=518×10N/m，=17×10N/m，=49×10N/m，系统的阻尼比均为001。在瞬时冲击激励下得到、和的冲击响应。

为验证多自由度系统的有效性，首先假设=0 kg，=0 N/m，则系统为2自由度系统，且结构的初始值

(15)

则此方程为求解初始值已知的常微分方程。通常，用1阶显示的微分方程组来描述：

(16)

(16)式实际为运动方程的状态空间表达形式，即(·)为

[,()]=A+

(17)

式中：

(18)

为单位矩阵。

本文用4阶定步长Runge-Kutta法求解微分方程，先定义4个附加向量：

(19)

式中：为计算步长，一般取常数。

则下一个步长的状态变量为

(20)

在已知初始值的情况下，可用(20)式一次次迭代，即可求解每时刻结构的响应。

在有限元软件中提取机匣处的垂向加速度响应，将信号经过处理后作为2自由度系统的基础激励，其中=5 130 kg，=518×10N/m，=170 kg，=17×10N/m，阻尼比设为001，对设备的响应进行理论分析，用MATLAB软件SimuLink中的仿真模块对2自由度系统进行仿真，求解示意图如图15所示。图15中，、分别为基座、机匣的加速度，、分别为基座、机匣的位移。

图15 SIMULINK 2自由度系统求解框图Fig.15 Block diagram of SIMULINK 2-DoF system solution

在Abaqus软件中建立相同2自由度系统的弹簧质量点模型，冲击载荷力施加在处，得到和处的位移响应。图16为Abaqus 2自由度计算模型。

图16 Abaqus 2自由度计算模型Fig.16 Abaqus 2-DoF calculation model

然后将理论计算结果与有限元仿真结果进行对比，其对比结果如图17所示。

图17 数值仿真结果与理论计算结果对比Fig.17 Comparison between numerical simulation results and theoretical calculation results

从图17中可以看出，理论计算结果与有限元软件仿真结果峰值相同，频率一致，少量误差可能来自有限元软件求解器的原因，从而说明了燃气轮机力学模型理论分析的正确性。

故可以使用MATLAB软件中的SimuLink仿真功能对(13)式进行求解，用变步长4阶、5级Runge-Kutta公式的ODE45算法，求解框图如图18所示，可得到、和的位移、速度及加速度响应。

图18 SIMULINK 3自由度系统求解框图Fig.18 Block diagram of SIMULINK 3-DoF system solution

采用与数值模型相同的信号处理方法对、和处得到的加速度信号进行分析，其分析结果如表3所示。

表3 燃气轮机质量点加速度测点主频率及能量占比Table 3 Main frequency and energy ratio of gas turbine mass points as acceleration measurement points

在研究力学模型载荷传递的规律时，分别选取3个质量点为测点。对测得的数据进行经验模态分解，记录经验模态分解后的主频率以及能量占比情况如表3所示，基本呈现与有限元模型相似的规律。在载荷由基座向内部转子传递的过程中，首先在基座部位激起的高频响应占比较高为0.71，随着载荷的向上传递，机匣的高频响应为0.62，在转子处的高频响应占比为0.54，高频响应呈现衰减的趋势。分析结果与有限元模型分析相一致，表明可以用此力学模型来预先对燃气轮机的响应进行预示。