火箭多台发动机并联推力传递结构动静耦合计算方法

2023-12-14曾耀祥

宇航总体技术 2023年6期

李斌，刘晖，曾耀祥，崔苗，吕军

(1. 大连理工大学工业装备结构分析国家重点实验室，大连 116024；2. 北京宇航系统工程研究所，北京 100076)

0 引言

大推力发动机多机并联技术带来了复杂结构的静力和动态力耦合作用及在联合载荷作用下的非线性传递问题[1-2]，成为了大型运载火箭研究的关键技术。目前对于大推力多台并联火箭发动机的力学研究较少。在已有的研究中，一般从火箭发动机传力结构的静力分析和动力学响应两方面开展。国外一些工程人员对Falcon 9重型火箭的八爪式机架、土星-5重型火箭的十字横梁结构分别从静力承载和动力学响应进行了设计与分析[3-4]。文献[5]研究了750 kN推力火箭发动机的静态和动态结构响应，分析了静态载荷对应力分布和结构自由振动的影响。我国对于大推力火箭的研究多数都针对于四机及以下并联火箭。文献[6]建立了某双台液体火箭发动机的整机动力学模型，通过有限元计算模态频率并与试验对比，结果吻合良好。文献[7]针对大推力液体火箭发动机结构对启动冲击环境的适应性问题，基于冲击动力学理论，提出了研究发动机启动冲击力学环境效应的非线性瞬态动力学分析方法。文献[8]采用子结构试验建模综合技术，对四机并联液体火箭发动机的动力学响应进行了研究。另外还有一些学者[9-11]从动力学优化角度，对大推力液体火箭发动机结构进行了研究。上述这些研究中只以静或动单方面的特性作为参考指标。文献[12-13]则从温度、应力及流体附加质量等影响因素，详细介绍了大推力火箭结构动力学分析所要考虑的问题。多台发动机并联火箭在工作过程中，由于机架之间存在静态力和动态力的复杂传递，巨大的静态推力势必会对结构的动力学特性产生影响。因此，若需要考虑机架和发动机在多台发动机工作下载荷及力学环境适应性问题，就要通过动静耦合的仿真计算方法考量结构的动静强度，进一步降低后续点火试验设计风险。

本文开展火箭7台并联发动机推力传递结构动静耦合计算方法的研究。通过研究在大推力作用下的几何非线性表达式，分析由发动机推力产生的预应力对机架结构频率的影响。在静力分析的基础上，开展振动模态分析及试验工况下随机振动响应分析，并对结果进行了分析，对火箭多台发动机并联推力传递结构动静联合试验有一定的参考意义。

1 动静耦合分析方法

1.1 几何非线性刚度阵

多台发动机并联火箭的推力传递结构是一种基于燃料贮箱箱底与火箭壳体结构联合传力结构，每台发动机又是通过由高强度材料焊接而成的杆系结构与上方贮箱和壳段相连。在发动机工作时，巨大的推力载荷容易使杆系结构发生非线性变形，进而影响其结构刚度，并使动力学特性发生变化。

非线性问题中主要包括材料非线性、几何非线性和状态非线性3类。其中几何非线性一般是由于变形后的状态与变形之前相差较大，必须在变形后的状态下重新建立物体的平衡方程，来保证问题的准确性，其本质就是结构的刚度阵不断变化。本文仅考虑几何非线性条件下机架梁结构的非线性动力学行为。几何非线性刚度阵的组成分为3个部分，即小位移线性刚度矩阵、初位移刚度矩阵和几何刚度矩阵的叠加构成。这里以两节点组成的杆单元为例，来模拟机架中杆系结构的非线性行为。

假设机架单根杆系结构长度为l，弹性模量为E，且两端节点位移矩阵可表示为q={ui，vi，uj，vj}T。受到轴向力产生的位移可分为轴向位移u和横向位移v，考虑杆的非线性，因此所产生的轴向应变εx表示为[14]

(1)

d={u,v}=[NiI,NjI]q

(2)

式中，I为与节点数量相关的单位矩阵，这里为二阶单位矩阵。

将式(2)代入式(1)中，可得

(3)

(4)

根据弹性力学基本方程，将各参数代入式(3)中，可以得到小位移线性刚度阵K0、几何刚度矩阵Kσ和初位移刚度矩阵KL分别为

(5)

(6)

(7)

式中，F为机架杆在火箭发动机工作时所受的轴向力。

根据式(5)～式(7)，可以推导出预应力和非线性条件下的单元刚度矩阵KT为

KT=K0+KL+Kσ=

(8)

综上，通过考虑机架杆系结构的大位移、初始应力两种非线性行为的影响，在理论上推导了非线性刚度阵的具体形式。

1.2 动力学分析理论

动静耦合计算方法的一般流程是先需要开展静力载荷的非线性分析，提取非线性刚度矩阵，以此刚度阵计算结构模态，然后分析随机振动的功率谱，并处理模态解，最后对计算结果进行后处理。因此，在有限元计算方法中，几何非线性作为分析的初始条件会被继承于模态分析中，模态分析中特征值的计算是一个线性过程。下面对考虑静载的随机振动计算过程进行推导。

根据上述非线性刚度阵代入到结构动力学自由振动方程中，可得

(9)

式中，M为结构的质量矩阵，C为结构的阻尼矩阵，KT为结构的刚度矩阵，x为结构的位移列阵。

由式(9)可得预应力条件下的结构特征值方程为

(KT-ω2M)Φ=0

(10)

式中，Φ为特征向量，ω为结构固有频率。

将位移列阵x分为自由和约束两部分，分别为xf和xr，结构动力学运动方程可分解为[15-17]

(11)

式中，q为节点载荷向量。

结构的自由位移向量xf可分解为拟静位移和动位移两项，即xf=xs+xd，xs可由下式求得

(12)

(13)

利用模态叠加法解耦可得

(14)

根据上述理论，首先求解出等效模态载荷，进而代入到解耦后的单自由度方程求解广义位移。

2 动静耦合仿真计算

2.1 多台发动机并联火箭推力传递结构

火箭7台发动机并联安装，排布方式如图1(a)所示。其传力方式为：发动机—多个并联小机架—主梁—箱底和过渡段—燃料贮箱，如图1(b)所示，每个发动机通过1个小机架共同连接在圆形主梁上，而每个小机架又由1个传力环和6根高强度钢管构成。主梁内圈与上端的锥形箱底连接，外圈与后过渡段连接，共同将推力传递到贮箱壁上。

(a)7台发动机排布方式

(b)推力传递结构各部件及传力方式图1 7台并联火箭布局和传力方式示意图Fig.1 The seven parallel rockets layout and force transmission method

多台发动机并联火箭有限元模型的建立中，考虑结构形式和网格适应性问题，箱底、过渡段等薄壁结构采用梁壳单元组合加筋形式进行建模；非薄壁结构比如小机架的传力环采用实体单元建模。此外，该模型中小机架和箱底壳段的材料分别采用高强度钢和铝合金，两种材料物性参数如表1所示。

表1 材料物性参数

该结构在发射条件下所受静载与约束可进行简化，见图2。将每个小机架传力环的底端面固定，在燃料贮箱上方端口加载向下的压力等效为发动机总推力，并在箱底处施加内压等效为贮箱内油压。

图2 火箭推力传递结构静态边界条件Fig.2 The static boundary conditions of rocket thrust transmission structure

2.2 多台发动机并联火箭推力传递结构预应力模态分析

考虑静载对于火箭机架结构固有频率的影响，首先开展无静载下七机并联火箭模态分析。只将7个传力环的底端面固定，计算仅在约束条件下的各阶典型振型下的频率；接着进行模态分析，底端的约束不变，计算在考虑静载时各阶典型振型下的频率，两种工况下的频率结果如表2所示。对比两种结构在各阶典型振型下的频率结果可以看出，结构在整体变形的振型下频率降低，其本质在于机架杆在静载作用下部分杆件产生了非线性大变形，几何刚度矩阵发生变化。箱底呼吸模态作为其中的局部模态，在内压的作用下，频率较无内压时升高，且升高幅度较大，说明内压对于局部呼吸模态的影响较为明显。另外，考虑静载后的各阶典型振型较无静载时无明显变化，考虑静载后的各阶振型如图3所示。

表2 有无静载时推力传递结构频率

(a)贮箱一阶弯曲

(b)贮箱扭转

(c)箱底呼吸

(d)机架一阶弯曲

(e)整体Z向拉伸图3 推力传递结构各阶振型图Fig.3 The modal of thrust transmission structure

2.3 多台发动机并联火箭推力传递结构随机响应分析

开展考虑静载的随机振动响应分析，轴向随机振动功率谱密度(Power Spectral Density，PSD)如表3所示，设定阻尼系数为0.02。基于上述模态分析，分别计算有无静载时随机振动响应的轴向均方根应力与加速度，计算结果如表4所示。

表 3 轴向随机功率谱密度

表4 各部件轴向最大均方根应力与加速度

以上计算结果可以看出，在静载下，机架整体的轴向最大均方根应力和加速度相比于无静载时略有提升。其中，小机架杆的提升幅度最为明显，原因是小机架杆作为发动机和上方贮箱的连接结构，其刚性对于整体模态的影响较大，当发生压力静载时，小机架杆的变形最为明显。因此，固有频率的变化导致小机架杆的响应变化最为明显。各部件的均方根应力云图如图4所示。主梁均方根应力最大位置均在与上方壳段连接的位置，均方根应力值约为105 MPa。小机架中最大均方根应力为根部位置，约为136 MPa。传力环在每个立筋处应力相对最大，其中，最大均方根应力在靠外侧的立筋根部，应力约为333 MPa。因此，在随机载荷下应注意结构中应力较大的地方，加载静载之后，尤其是一些承力主体的根部连接部位，易产生应力集中，需进一步加强设计。