柔性整流罩地面展开试验仿真分析与飞行状态分离预测

0　引言

运载火箭在飞行过程中处于气动力、气动热、噪声等恶劣环境中，整流罩作为火箭的重要组成部分，可以保护卫星及有效载荷免受外界环境的影响[1]。当运载火箭飞出大气层后，整流罩在分离系统的作用下旋转抛离。随着运载火箭尺寸的不断增大，整流罩分离过程中的呼吸运动也相应增大，从而影响罩内的有效包络空间，对安全分离造成威胁。对此，马忠辉[2]、张大鹏等[3]、雷勇军等[4]基于柔性多体动力学采用商业软件系统地研究了柔性整流罩的分离特性，并分析了主要因素对分离运动与包络空间的影响，但并未考虑空气与结构的耦合效应。

受国内现有真空试验设施的规模限制，难以开展真空条件下的整流罩地面分离试验。现有的试验无法避免空气阻力影响，需要同时考虑刚性与柔性(简称刚柔)耦合和流体与固体(简称流固)耦合问题，通过计算空气阻力作用下的柔性整流罩展开过程，与试验数据进行对比修正，得到准确的整流罩有限元模型，以进行真空环境下的分离预测。张小伟等[5]采用侵入式边界方法求解了非定常气动载荷，并转化为节点压力作用于柔性整流罩上求解柔性体的动态响应。朱学昌等[6]通过计算流体动力学方法计算不同角度下气动网格点上的气动力，并转化到结构网格上进行加载，利用有限元软件进行了刚柔耦合计算。落寿[7]利用动网格法计算了气动阻力，通过逐步耦合的方式获得了分离特性和呼吸变形。顾名坤等[8]提出了多因素计算方法，建立了动力学方程对分离试验进行仿真计算。以上方法均未能实时考虑结构变形对流场的影响，对分离特性的预测也有所偏差。

耦合欧拉- 拉格朗日 (CEL) 算法结合了欧拉网格和拉格朗日网格的优点，能够有效地解决带有较大变形的流固耦合问题。目前，CEL方法已成功应用于水下爆炸[9]、两相流动[10]、材料成型[11-12]、岩土工程[13-14]等方面。李刚等[15]基于CEL算法对整流罩的地面分离试验进行了仿真，通过弹簧效能系数折减法将整流罩简化为刚体，但并未关注整流罩的弹性变形。Liu等[16]利用非线性瞬态动力分析MSC/Dytran软件建立了整流罩的动力学流固耦合模型，采用通用耦合算法分析了气动力对振动频率和质心速度的影响，但计算结果未能与试验结果进行对比。

本文基于CEL算法，利用有限元Abaqus软件对整流罩地面展开试验进行仿真，并通过与现有试验获得的运动特性数据和变形数据进行对比，验证流固耦合模型的正确性；利用验证的模型对实际飞行状态下的分离进行仿真预示，获得了整流罩内的有效包络空间，重点分析了空气阻力和轴向过载对运动特性和呼吸变形的影响。

1　计算方法及模型建立

1.1　流固耦合CEL算法

Noh[17]为解决带有移动流体边界的二维欧拉流体动力学问题首次提出了CEL算法，结合了拉格朗日算法中能够精确追踪物质界面的优点和欧拉算法在处理大变形问题上的优势，适用于计算含有大变形和流体流动的情况。

利用CEL算法时，在第n个时间步tn中，假设已知拉格朗日域的位置和所有材料的状态，并已知二者的边界，为求第n+1个时间步tn+1的状态，通常分为如下3个步骤进行：1)通过当前状态下欧拉材料中的压力求解作用于拉格朗日区域表面的作用力，进而求解拉格朗日网格的运动情况；2)为控制欧拉节点周围的材料体积，应用高斯理论，根据拉格朗日网格的位置构建欧拉- 拉格朗日交界面上不规则形状网格的离散欧拉方程，确定拉格朗日区域外部的欧拉域部分；3)求解离散的欧拉方程，以得到欧拉材料中的压力，进而得到在下一时间步中作用于拉格朗日网格上的作用力。在上述计算过程中，拉格朗日网格边界的速度为欧拉区域提供了动力学约束，同时欧拉材料的压力提供了作用于拉格朗日区域的作用力，并以此实现了耦合。

在有限元分析软件Abaqus/Explicit中，利用中心差分格式对非线性系统的微分方程进行求解。在每个时间步中，利用基于罚函数法的通用接触算法实现欧拉材料和拉格朗日材料间的耦合作用。罚函数方法本质上是在耦合问题中将两种网格的相互作用关系看作是弹簧- 质量系统，如图1所示，欧拉物质界面上锚定点和拉格朗日单元边界上节点的相互作用力[18]表示为

Fp=kpdp，

(1)

式中：kp为罚刚度系数，与时间步长、材料属性相关；dp为相对位移矢量。仅当拉格朗日节点在欧拉材料中发生贯穿时，才考虑相对作用力，即

d·n>0,

(2)

式中：d为拉格朗日节点的实际位移；n为拉格朗日网格外表面法向量。

图1　CEL算法中的罚函数方法示意图[18]
Fig.1　Illustration of penalty function method for CEL

根据上述罚函数方法得到的相互作用力，可进一步利用欧拉单元基本方程求解第i节点上的欧拉节点力FE,i，以及拉格朗日节点所受的反作用力：

(3)

1.2　柔性整流罩的有限元模型

以直径4.2 m的对称分瓣式整流罩为研究对象，每半罩由端头、锥段、柱段、倒锥段和分离铰链组成，如图2所示，并规定整流罩的轴向为x方向，分离方向为y方向。分离时，半罩在预压缩弹簧的作用下旋转至一定角度，铰链与铰座脱离，完成整流罩的抛罩过程。

图2　整流罩半罩结构示意图
Fig.2　Illustration of half-fairing structure

为了准确地捕捉整流罩各主要部位的变形情况，对整流罩进行了精细建模，包括各部段间的连接框、加强筋、倒锥段的内外桁条等。其中：锥段、柱段由蜂窝铝夹层板构成，其夹层为正交各向异性材料；其他部分均由各项同性材料构成；铰链处的刚度相对较大，尺寸相对较小，且结构复杂，为保证计算效率，将其简化为刚体。夹层板的铝蜂窝部分划分为实体单元，上下面板为壳体；分离铰链划分为实体单元；其他部分均为壳单元。半罩模型中共划分65 190个结构网格单元，半罩的有限元模型如图3所示。

图3　整流罩半罩结构的有限元模型
Fig.3　FE model of half-fairing structure

1.3　分离装置的处理

半罩由预压缩弹簧提供分离动力，使用连接器对分离装置进行建模。分别选取弹簧两端安装位置的对应节点作为参考点，设置两参考点间的作用力方向为连线轴向，作用力大小与两点间的相对位置耦合。预压缩弹簧在压缩时产生两点连线的轴向作用力，为整流罩提供分离推力，恢复原长后无相互作用力。假设弹簧原长为l0，刚度系数为k，弹簧的预压缩量为Δl，则弹簧两端点的初始安装位置距离l为

l=l0-Δl.

(4)

两参考点距离x与弹簧作用力F的关系为

(5)

1.4　流体域的建立

为准确模拟地面大气环境，空气环境采用理想气体状态方程描述[15]如下：

p+pa=ρR(θ-θz)，

(6)

式中：p为气体压强；pa为环境压强，取1.013×105 Pa；ρ为气体密度，取1.17 kg/m3；R为气体常数，取283.0 J/(kg·℃)；θ为气体温度，取室温25 ℃；θz为绝对温度-273 ℃. 此外，设定流体的比热γ为717.0 J/(kg·℃)，黏性系数μ为17.9×10-6 Pa·s.

为避免流场边界对计算结果的影响，首先需要寻求合适的流场尺寸。以整流罩对称轴为基准，固定整流罩在流场中的相对位置，成比例扩大流体域，得到结构- 流场装配效果如图4所示。边界条件为欧拉域外表面法向速度等于0，并在-x轴方向施加重力(取地面重力加速度为9.8 m/s2)。为提高计算效率，选择仅包含整流罩外形尺寸的简化刚体模型进行计算，设定简化刚体壳单元尺寸为150 mm；使用8节点线性欧拉单元对流体域网格进行划分，欧拉单元边长设置为180 mm，计算得到整流罩在1.5 s内头部顶点的x方向位移如图5所示。从图5中可以看出，当流体域尺寸为30 m×30 m×30 m和33 m×33 m×33 m时计算结果差别很小，可忽略流体边界对计算结果的影响。为保证计算效率，选取30 m×30 m×30 m流体域尺寸进行柔性整流罩的流固耦合计算。

图4　流场与整流罩半罩的装配图
Fig.4　Assembly diagram of fluid field and half-fairing

图5　不同流体尺寸下刚性整流罩顶点的x方向位移
Fig.5　Top node displacements of rigid fairing along x-axis for different fluid sizes

为选取合适的欧拉网格尺寸，对6种不同欧拉网格情况进行计算，欧拉单元边长分别为150 mm、180 mm、200 mm、300 mm、400 mm、500 mm；将流场域设定为30 m×30 m×30 m，其他条件不变，得到整流罩在1.5 s内头部顶点的x方向位移如图6所示。由图6可见，欧拉网格尺寸在150 mm和180 mm时，所得x方向位移的差别很小，此时加密网格对于求解结果的精度并没有显著提升，可忽略网格的依赖性[16]。

图6　不同网格下刚性整流罩顶点的x方向位移
Fig.6　Top node displacements of rigid fairing along x-axis in different meshes

2　地面试验与仿真

2.1　地面分离试验

图7　测量系统布局及测点位置示意图
Fig.7　Layout of measuring system and the position of measuring points

基于双目立体视觉原理，采用高速相机、数据处理计算机、光源及其他辅助设备构成的测量系统，对整流罩的整体运动轨迹和局部呼吸运动进行测量,测量系统布局及4个测点A、B、C、D的位置如图7所示。试验中通过追踪整流罩表面具有高反光特性的测量靶标，并进行图像处理和数据解算，成功获取了整流罩分离过程中的质心运动数据、测点呼吸运动数据。

2.2　地面试验仿真计算

综合考虑计算效率和精度，在弹性整流罩的流固耦合计算中，欧拉域尺寸选为30 m×30 m×30 m，划分欧拉单元尺寸为180 mm，欧拉域网格数量为465万。由于计算时间较长，时间步较多，为减小数值耗散导致的误差，采用双精度进行计算，计算时间为3 s.

2.2.1　位移

由于质心不在结构上，可提取整流罩铰链中心点位移数据和头部顶点位移数据，通过坐标转换得到整流罩半罩质心的x、y、z方向的位移与时间曲线，如图8～图10所示。过顶表示整流罩质心达到了最高点，接下来整流罩将依靠自身重力和铰链约束实现分离，而成功过顶说明了分离弹簧推力设计的合理性，因此需要对过顶时刻进行详细分析。通过提取图8中x方向的位移最大值x=55 mm，可得整流罩的过顶时刻为0.81 s，试验所得过顶时刻为0.91 s，相对误差为10.99%. 图10表明，整流罩在z方向也存在着一定程度的侧向晃动，范围在3.5 cm以内。从图8、图9中仿真结果与试验结果的对比可以看出，仿真计算中整流罩的运动稍快，这是由于网格密度不足导致了计算过程中难以捕捉更精密的空气流动细节，造成了误差。

图8　整流罩半罩质心的x方向位移
Fig.8　Displacement of mass center of half-fairing along x-axis

图9　整流罩半罩质心的y方向位移
Fig.9　Displacement of mass center of half-fairing along y-axis

图10　整流罩半罩质心的z方向位移
Fig.10Displacement of mass center of half-fairing along z-axis

2.2.2　转角

通过坐标转换和几何关系，可得整流罩半罩绕铰链转动的转角与时间曲线(见图11)。由图11可见，在0.5 s内整流罩转动较为平缓，在重力作用下转动逐渐加快。

图11　绕铰链转动转角与时间曲线
Fig.11　Angle of rotation along hinge axis

2.2.3　角速度

根据转角数据，可确定整流罩半罩绕铰链旋转的角速度如下：

(7)

图12　绕铰链转动角速度与时间曲线
Fig.12　Angular velocity of rotation along hinge axis

式中：ω为角速度；θn为第n步的转角；Δt为时间步长。由转角与时间曲线和(7)式得到半罩旋转角速度与时间曲线，如图12所示。由图12可见：初始阶段在预压缩弹簧作用下，半罩的角速度增大；行至整流罩过顶时，角速度趋于平缓；随后在重力作用下角速度持续增大；过顶时刻0.81 s对应的整流罩半罩的过顶角速度为16.78 °/s，试验中过顶时刻0.91 s所对应的过顶角速度为16.52 °/s，相对误差为1.57%.

2.3　呼吸运动

提取与试验中4个呼吸运动测点相对应的有限元模型节点的z方向位移数据，排除整流罩在z方向的刚性位移部分，得到整流罩的呼吸变形如图13所示。由图13可以看出：在运动初始时刻，由于弹簧预应力的存在和气动阻力影响，倒锥段的变形较大；随着整流罩的转动，呼吸变形幅度逐渐减小；半罩倒锥段下端振幅较大。

图13　4个测点的呼吸变形曲线
Fig.13　Breathing deformations of 4 measuring points

将流固耦合仿真得到的测点呼吸变形最大值与试验数据进行对比，得到相对误差如表1所示。由表1可知，试验与仿真的误差在15%以内，说明该有限元模型能够较好地刻画整流罩半罩的实际变形情况和变形特点，在真空分离预测中具有一定的通用性。

3　真空状态下分离预测

在实际飞行的分离过程中，整流罩近似处于真

表1　测点的呼吸运动最大值对比
Tab.1　Comparison on maximum breathing deformations of 4 measuring points

空状态，无空气阻力影响，地面试验很难预测其分离特性。采用经CEL算法验证过的有限元模型对真空状态下的分离情况进行预测，可获取可靠的分离特性，并对整流罩分离的天地差异进行分析。

3.1　模态分析

为分析整流罩半罩的整体振动特性，对自由状态下半罩模型的固有频率和振型进行分析(见表2)。其中，1阶振型为扭转振动，2阶振型为整体呼吸振动，3阶振型为上半部分、下半部分交错的呼吸振动，4阶振型为扭转、呼吸振动的叠加。

表2　半罩前4阶模态

Tab.2　First 4 rank vibration modes of half-fairing

3.2　分离特性对比

计算轴向过载为1 g时真空状态下的分离转角与时间曲线、角速度与时间曲线(见图14、图15)，与流固耦合仿真结果对比可知：在空气阻力影响下，过顶时刻由0.60 s减缓至0.81 s，推迟了0.21 s；过顶角速度由27.20 °/s减缓至16.78 °/s，降低了38.31%；在2.13 s时，半罩的角速度曲线出现拐点，角加速度发生突变，说明该时刻为整流罩半罩的脱钩时刻。

图14　有无空气情况下半罩的转角与时间曲线
Fig.14　Angles of rotation of half-fairing along hinge axis with and without air resistance

图15　有无空气情况下半罩的角速度与时间曲线
Fig.15　Angular velocities of rotation of half-fairing along hinge axis with and without air resistance

3.3　呼吸运动对比

由于脱钩后整流罩倒锥段的呼吸变形对包络空间影响不大，故对2.0 s内有无空气情况下的呼吸运动进行对比分析，如图16所示。由图16可以看出，4个测点的运动特征一致，都在刚体运动的平衡位置进行内外往复的呼吸运动。真空状态下，整流罩倒锥段运动周期(图中记为Tv)为0.30 s，频率为3.33 Hz，与1阶固有频率较为接近，说明整流罩的低阶振型对呼吸运动起主导作用；在流固耦合状态下，倒锥段运动周期(图中记为TCEL)为0.42 s，频率为2.38 Hz，呼吸运动频率降低了28.50%. 这是由于在计算流固耦合问题时，空气作为附加质量计入了系统振动方程中的惯性项，对系统的运动频率造成了影响。这也间接地体现了流固耦合系统的湿模态特性。

由图16可知，空气阻力的存在增加了呼吸运动的振幅。提取测点的呼吸变形最大值进行对比，并计算最大值的差值百分比，结果如表3所示。由表3可知，4个测点对应的有限元节点呼吸运动幅度均有所增大，最低增幅为46.44%. 因此，在计算大型结构的振动时，需要考虑空气耦合对振动特性的影响。

表3　有无空气情况下测点的呼吸运动最大值对比
Tab.3　Comparison of maximum breathing deformations of 4 measuring points with and without air resistance

4　轴向过载的影响

整流罩抛罩时的过载取决于运载火箭的总体设计。下面通过计算不同轴向过载下的整流罩分离特性和呼吸运动情况，对影响规律进行研究。

4.1　轴向过载对分离特性的影响

采用相同的整流罩半罩模型，分别计算轴向过载为1.0 g、1.2 g、1.4 g、1.6 g、1.8 g时的分离特性，得到转角与时间曲线(见图17)和角速度与时间曲线(见图18)。由图17和图18可见：轴向过载越大，过顶角速度越小，过顶时刻逐渐延迟；过顶之后，轴向过载越大，转动速度越快，脱钩时间逐渐提前。

4.2　轴向过载对呼吸运动的影响

图19为上述5个工况整流罩半罩测点的呼吸变形情况。由图19可见，4个测点所体现的运动频率均为3.3 Hz，并未发生改变。4个测点呼吸变形最大值与轴向过载的关系如图20所示。从图20可见，随着轴向过载的增大，呼吸变形幅度逐渐增大，且最大值均发生在分离后的第1个振动周期内。在1.0 g～1.8 g的较小范围内，呼吸变形最大值近似呈线性增加。其中，靠近底端的C、D两点与A、B两点相比，幅值变化更为明显，这与整流罩半罩的1阶振型较为相似。

图16　有无空气情况下半罩的呼吸运动对比
Fig.16　Comparison of breathing deformations of half-fairing with and without air resistance

图17　不同轴向过载下的转角与时间曲线
Fig.17　Angle of rotation along hinge axis under different axial overloasds

图18　不同轴向过载下的角速度与时间曲线
Fig.18　Angular velocities of rotation along hinge axis under different axial overloads

5　结论

本文采用Abaqus软件对大型柔性整流罩进行了分离动力学建模，并基于CEL算法对地面展开试验进行了流固耦合仿真计算，通过与试验数据的对比验证了模型的正确性。利用验证的模型对飞行(真空)状态下的展开情况进行预测，并分析了空气阻力和轴向过载对整流罩分离的影响。本文的技术方案综合考虑了运动、变形和流固耦合作用，得到的模型更具通用性。所得结论如下：

1)采用CEL算法能求解空气和整流罩的耦合运动问题，可用于验证整流罩分离方案的正确性。

2)在空气阻力的作用下，整流罩的过顶角速度降低38.31%，过顶时刻和脱钩时刻均有所延迟；由于耦合时空气的附加质量效应，整流罩呼吸运动频率降低、幅值增大；在真空状态下，整流罩的呼吸运动频率与1阶固有频率近似。

图19　不同轴向过载下的呼吸运动对比
Fig.19　Comparison of breathing deformations of half-fairing under different axial overloads

图20　呼吸变形最大值与轴向过载的关系
Fig.20　Relation between maximum value of breathing deformation and axial overload

3)随着轴向过载的增加，脱钩时刻提前，呼吸运动幅值明显增大。因此，在飞行器的总体设计阶段，为避免整流罩与罩内的有效载荷产生碰撞，需要综合考虑振型和过载对可用包络空间的影响。