近年来, 诸多国家开展了一子级垂直起降的重复使用火箭相关研究[1], 美国太空探索公司的猎鹰 9 号(Falcon-9)运载火箭和蓝色起源公司的新谢泼德(New Shepard)火箭已经成功地实现火箭一子级的垂直回收和重复使用。火箭的一子级在上升段完成分离后, 通过气动减速和发动机反推减速, 将速度降至较低水平[2]。但是, 由于偏差的存在, 难以控制到零速, 使得火箭必须采用着陆缓冲机构来降低冲击载荷, 提高着陆稳定性[3]。猎鹰 9 号火箭和新谢泼德火箭均采用腿式着陆缓冲机构(简称着陆腿), 其中猎鹰 9 号火箭采用倒三角式着陆腿。垂直起降运载火箭的着陆问题与星球着陆探测器的着陆问题在动力学上具有一定的相似性, 倒三角式着陆腿在 Luna 系列探测器[4]、Surveyor 1月球探测器[5]和 Viking-1 火星探测器[6]等星球着陆探测器中也被广泛采用, 如图1 所示。
图1 倒三角着陆腿的应用
Fig.1 Application of inverted triangle landing leg
着陆过程中, 着陆腿产生横向振动, 并导致较大的弯矩, 该弯矩的计算是着陆腿设计的关键。工程上常采用有限元法对着陆过程中的动力学响应进行数值计算[7–8]。有限元法可实现较高的计算精度,但无法得到着陆腿载荷与着陆腿设计参数的解析表达式。在着陆腿构型参数设计阶段, 需要着陆腿参数设计与有限元分析反复迭代, 反复建立不同参数下的有限元模型, 导致设计效率较低, 且参数的影响也不够直观。为提高参数设计效率, 近年来众多学者使用代理模型方法优化着陆腿参数。Yue 等[9]和雷波等[10]利用径向基(radial basis function)代理模型, 对着陆腿可采用的多种缓冲器的参数进行优化。Li 等[11]、王家俊等[12]、吴宏宇等[13]以及刘学翱[14]利用响应面(response surface)代理模型, 对着陆腿的缓冲器参数和构型参数等多个设计参数进行优化。代理模型的精度依赖于样本数量的提高, 样本通常采用数值仿真得到, 故随着优化变量的增多,需要的样本数量和计算量迅速增大。并且, 代理模型通常基于统计得到, 物理意义较不明确, 结果的可解释性较差。
着陆腿在着陆冲击作用下的横向冲击响应问题, 可等效为简支梁在惯性力激励下的横向振动问题。火箭的着陆腿为多节薄壁圆筒结构, 直径沿轴向阶梯变化, 且在节间连接处有端部加强和锁定机构等, 可将其简化为分段等截面和含集中质量的简支梁。针对非均匀简支梁的横向振动问题, 众多学者推导了多种截面变化形式的简支梁横向振动的解析解, 包括含集中质量的等截面简支梁[15–16]、半径随轴向坐标多种幂次变化形式的简支梁[17–18]以及截面沿轴向阶梯变化的简支梁[19]等。上述解析方法均未给出分段等截面和含集中质量的简支梁横向振动的解析解。
面对复杂梁截面形式下梁的横向振动问题时,一般采用假设模态法, 近似地求解振型函数和频率。选取基函数时, 由于非均匀梁横向振动的解析解形式较为复杂, 难以进一步推导半解析解, 故通常采用相同边界条件下 Euler 梁自由振动的振型作为基函数。Cheng 等[20]采用假设模态法建立梁的横向和纵向耦合振动模型, 研究旋转悬臂梁问题中角速度对悬臂梁固有频率的影响。潘旦光等[21]以等截面 Euler 梁的自由振动模态为基函数, 采用假设模态法分析变截面箱型梁振动中的剪力滞效应, 结果表明取前 12 阶基函数时已经具有较高的精度。Kong 等[22]采用假设模态法建立含有两个感应电机的简支梁的横向振动模型, 分析电机参数对系统振动的稳定性的影响。Luo 等[23]采用假设模态法建立弯曲梁的面内和面外振动模型, 并利用该模型研究列车–轨道系统中轨道曲率对系统动特性的影响。
基于上述假设模态法对非均匀梁振动的研究,本文提出一种计算效率高、可解释性强的半解析方法。针对倒三角式着陆腿, 分析着陆过程中着陆腿上的载荷。在轴向, 考虑含缓冲的着陆过程中的静力平衡给出轴向应力, 并给出不同摩擦力方向下的最大轴向应力。在横向, 将着陆过程的激励等效为刚体运动引起的分布惯性力, 通过假设模态法, 近似地求解主腿的振动模态, 利用模态叠加法求解动力学响应, 并考虑横向弯曲变形引起的纵向缩短导致的耦合, 给出横向振动弯矩的半解析计算方法。随后, 通过实验来验证本文提出的计算方法。最后,以满足承载力下着陆腿系统的质量最小为目标, 利用该方法, 分析着陆腿的最优构型, 并进一步分析材料强度对着陆腿最优构型的影响。
1.1 倒三角式着陆腿设计
重复使用火箭采用 4 套倒三角式着陆腿, 每套着陆腿由一条主腿、两条副腿、一个足垫及箭体连接构件等组成。主腿为多节套筒式伸缩设计, 在上升段, 着陆腿处于收拢状态, 紧靠箭体, 在着陆前,着陆腿展开锁定, 实现着陆缓冲的功能, 其构型如图2 所示。左侧着陆腿为收拢状态, 中间和右侧着陆腿为展开状态。主腿、副腿与箭体、足垫均通过球铰相连, 在静力状态下, 主腿、副腿均为轴向受力杆件。主腿共有 4 节, 每节为薄壁圆管, 在收拢状态下, 下方着陆腿套于上方着陆腿内, 展开时套筒逐节展开。主腿的各节间有挡环和锁定机构, 在展开到位后实现两节间的锁定。主腿与足垫相连的节内装有金属蜂窝缓冲器, 在碰撞时, 通过压溃的塑性变形吸收能量, 降低冲击加速度。副腿为等截面薄壁圆管。
图2 倒三角式着陆腿示意图
Fig.2 Figure of inverted tripod landing gear
1.2 几何约束
火箭着陆前, 着陆腿已经完成展开锁定, 展开后箭体、主腿、副腿和足垫组成静定结构, 等效力学模型如图3 所示。图3(a)为主腿与箭体轴线构成的平面, 图3(b)为俯视图。设着陆体质量为 MB, 箭体上主腿、副腿安装点与箭体轴线的距离分别为rL1 和 rL2, 着陆腿跨距为 xB。设展开后主腿的长度为 l1, 主腿与箭体轴线的夹角为 β。设副腿的长度为 l2, 副腿与箭体的连接点靠近箭体底部, 展开后的支撑高度为 h, 主腿、副腿与箭体连接点在 y 方向的距离为 h12, 副腿与箭体轴线的夹角为 α, 副腿与图3(a)中 x-o-y 平面的夹角为 γ。为便于分析, 令
代表在 x-o-y 平面内副腿与箭体轴线的夹角。着陆腿上述参数之间的几何约束为
设主腿第 i 节的长度为 l1i, 半径为 r1i, 壁厚均为 t, 则各节的面积和惯性矩分别为 A 1i =2 π r1it,Ili = 
。将主腿各节间的挡环和锁定机构等效为附加集中质量, 如图3 中 m11~m13 所示, 将缓冲器也等效为附加集中质量, 如图3 中 m14 所示。设副腿的半径为 r2, 壁厚为 t, 与主腿相同, 则其截面积为A2 = 2πr2t, 截面惯性矩为 I 2 =
。
图3 着陆腿几何参数
Fig.3 Geometry parameters of landing leg
1.3 着陆过程的刚体运动
火箭着陆时, 足垫首先与地面发生碰撞, 碰撞后足垫与地面保持接触, 主腿和副腿分别绕其与箭体的连接点发生转动, 主腿内的缓冲器压溃, 导致主腿长度缩短, 火箭向下运动, 且火箭质心加速度向上, 直至减速至箭体质心静止。
由于足垫与地面的碰撞无缓冲, 碰撞时间很短,足垫与地面碰撞前后, 着陆腿构型参数变化可以忽略[24–25]。由于缓冲器的作用, 箭体质心加速度为有限值, 故足垫与地面碰撞前后箭体的速度变化同样可忽略。箭体减速下降的过程中, 令 v 表示箭体的下落速度, v0 表示触地时刻的速度, 由于着陆过程箭中体向下运动, 定义 v 向下为正。分析机构的刚体运动, 并忽略弹性变形, 着陆过程中着陆腿和箭体连接点速度与箭体相同, 着陆腿与足垫连接点的速度与足垫相同, 则 l2, rL1, rL2 和 γ 保持不变, l1, xB,
, β 和 h 由于刚体运动而发生变化。则由式(1)得副腿角速度
、主腿角速度β˙和主腿长度变化率
:
令a=-
˙代表箭体 y 方向的加速度, 向上为正,由式(1)和(2)可得足垫与地面保持接触。箭体减速下降的过程中, 副腿的角加速度α˙˙和主腿的角加速度
为
着陆过程中, 主腿缓冲器的压溃距离Δl1相对于主腿长度 l1 较小, 通常有 Δ l1 /l1＜ 1 0%[10,12–14], 即着陆腿系统的构型变化较小, 故采用式(2)和(3)计算角速度和角加速度时, 忽略角度
和β的变化,且近似认为Δl1<<l1。由式(2)可得, 火箭减速过程的下降高度
由于金属蜂窝缓冲器压溃应力稳定, 近似认为火箭着陆过程中加速度为 a 常值, 则有 a =v02/2Δh, 进一步得 a>>v02/l1。由于着陆过程中火箭做减速运动, v≤v0, 故认为 a>>v2/l1, a>>v2/l2, 则角加速度
˙和
˙可近似为
着陆时, 主腿通过压缩金属蜂窝实现对主腿轴向力的缓冲, 降低箭体的加速度。本文分别分析碰撞过程中的轴向力和横向振动弯矩, 计算着陆过程中载荷着陆腿的最大应力。
2.1 轴向力模型
设 T1, T21 和 T22 分别为主腿和两条副腿的轴向力, 为拉力时取正值。N 为地面对足垫的支撑力, f为摩擦力, 令 μ = f / N 为摩擦系数, 其上界为 μmax, 即满足 0≤μ≤μmax。在图3 所示的坐标系下, 支撑力N 沿 y 方向, 摩擦力在 z-o-x 平面内, 设 φ ∈[-1 80°,180°]为摩擦力与 x 轴正方向的夹角, 如图4 所示。着陆腿展开后, 箭体、主腿、两条副腿和足垫组成静定结构, 在足垫节点处的静力平衡关系为
图4 静力平衡示意图
Fig.4 Diagram of the axle force equilibrium
主腿内装有金属蜂窝缓冲器, 通过金属蜂窝的异面压缩, 降低着陆过程的轴向力。金属蜂窝具有压溃应力稳定的特点[26–27], 则着陆过程中, 主腿轴向力 T1<0 可近似为常值, |T1|为金属蜂窝的压溃力,为待设计参数。由式(5)可得支撑力 N 和副腿 1 的拉力T21 分别为
由于 T22 与 T21 对称, 故仅需分析 T21。令 φ0 代表使T21 取极值的角度 φ, 则 φ0 满足
将式(8)的解代入式(7), 可得 T21 的极大值T2,max 和极小值 T2,min, 并且 T2,max 和 T2,min 均关于 μ 单调, 即 ∂ T 2,max /∂μ>0, ∂T 2,min/∂μ<0, 故计算副腿的最大拉力 T2,max 和最大压力-T2,min 时, 取 μ = μmax。
由式(6)可知, 着陆腿需满足
＞μm ax , 否则φ=180°方向主腿的轴向力为拉力, 难以通过缓冲器实现缓冲。在满足该条件的情况下, φ = 180°, μ =μmax时N取极大值。令amax代表着陆冲过程中纵向的容许加速度, 则支撑力 N=MB(amax+g)/4, 主腿蜂窝压溃力设计值为-T1, 计算公式为
由于触地前着陆腿的轴向力仅来自着陆腿系统的重力, 故似为 0。着陆过程中, 由于缓冲器的作用, 轴向力近似为定值 T, 则着陆腿的轴向加载过程近似为阶跃激励。阶跃激励下, 最大响应不大于 2 倍激励, 即 2T。由于着陆腿在轴向为多节结构,其刚度和振动模态受节间锁定机构的影响较大, 因此参数设计阶段难以对锁定机构准确建模, 导致难以精确地计算其在冲击下的响应, 故在强度设计中考虑最大轴向力为 2T, 该假设将导致应力的计算结果偏大。
2.2 横向振动模型
将着陆过程横向振动的激励考虑为两部分: 1)足垫与地面碰撞前, 着陆腿相对箭体静止, 碰撞后着陆腿在很短时间内发生转动, 引起着陆腿横向振动, 该部分称为“冲击激励”; 2)足垫与地面碰撞后,在着陆腿缓冲器的作用下, 箭体做减速运动, 同时着陆腿的转动过程具有角加速度, 上述刚体运动加速度产生的惯性力引起着陆腿横向振动, 该部分称为“稳态激励”。
本文将着陆腿的横向振动问题等效为简支梁的横向振动问题进行分析。建立单条着陆腿的局部坐标系, 以着陆腿与箭体的连接点为原点, 梁的轴线方向为 x 轴, 垂直于着陆腿斜向上的方向为 y 轴, 设着陆腿长度为 l, 则着陆腿与足垫的连接点在 x=l处。由于着陆腿的截面尺寸沿轴向变化, 设等效简支梁的截面积为 A(x), 截面惯性矩为 I(x), 材料的密度为 ρ, 模量为 E。以主腿为例, 考虑附加集中质量后, 上述等效简支梁模型如图5 所示。主腿每段的长度为 l1i, 则节间等效质量块的位置坐标为 x1j =l11+…+l1j, j = 1~3, 缓冲器等效质量 m14 的位置为x14 = x13+l14,1。由于位于梁端点的集中力不产生弯矩, 故 1.2 节建立的几何模型中忽略了靠近端点处构件的质量对弯矩的影响。
图5 着陆腿的等效简支梁动力学模型
Fig.5 Equivalent beam dynamic model of landing leg
由于副腿为等截面结构且无附加集中质量, 故本文以主腿为例建立横向振动的等效简支梁模型。分析副腿时, 采用主腿的等效简支梁模型, 令附加集中质量为 0, A(x)和 I(x)为副腿的截面参数 A2和 I2,并将着陆腿长度、角速度和角加速度替换为副腿对应的参数。
首先, 建立“冲击激励”的方程。我们认为足垫与地面的碰撞时间很短, 碰撞过程着陆腿横向振动的位移为零, 采用冲量形式描述碰撞过程, 且碰撞过程常规力的冲量可忽略[24–25]。假设足垫与地面的碰撞时刻为 t0, 则碰撞过程主腿的角加速度为
其中δ()为 Dirac-delta 函数, 
为足垫与地面碰撞后主腿的角速度, 设触地时箭体下降速度为 v0, 由式(2)可得
将碰撞过程中着陆腿视为非惯性系, 则碰撞引起的横向惯性力激励见图6。记碰撞引起的惯性力为qc(x):
图6 冲击激励的等效模型
Fig.6 Equivalent model of the impact excitation caused by the touchdown
然后, 建立“稳态激励”的方程。着陆过程中,火箭质心的加速度为 a, 足垫与地面碰撞后到箭体静止的过程中, 主腿的角加速度
如式(4)所示, 将着陆过程中着陆腿视为非惯性系, 则质心加速度、角加速度和重力的合力在着陆腿横向的分量记为q1(x):
其中第一项由角加速度引起, 第二项由箭体质心平动和重力引起, 由于缓冲器的作用减速过程中 a 近似为定值, 由式(4)可得
˙正比于 a, 故忽略减速过程中sinβ的变化后, q1(x)对于时间可近似为定值。
着陆过程中, 着陆腿所受的横向激励为 q(x) =qc(x)+q1(x)。对图5 所示的振动模型, 忽略系统的阻尼, 其横向振动方程为
其中, E 为着陆腿材料的弹性模量, I(x)为着陆腿截面的惯性矩。
本文通过模态叠加法求解横向的冲击响应。由于主腿的等效简支梁模型截面非均匀， 且包含集中质量, 其自由振动模态难以解析表达, 故采用假设模态法近似求解横向振动模态, 基函数采用等截面简支梁的自由振动的模态, 第 ξ 阶为 Y ξ( x ) = sin(ξπxl1)。令 Y(x)表示基函数列向量, η(t)表示广义坐标列向量, 则主腿的横向运动表示为
由于系统具有无穷维自由度, 应取n→∞, 但在梁振动问题中, 取少量低阶基函数即可获得精度较高的解[21–22]。本文取 n = 20。
考虑横向弯曲变形引起的纵向缩短导致的轴向力与横向振动的耦合[28], 令 ul 代表缩短的长度, ul>0 代表轴向缩短, 当采用式(13)表示 y(x, t)时, ul 可表示为
将轴向力 T 考虑为恒定值, 以(t)为广义坐标,则系统动能 Ek 和势能 Ep 分别为
其中, Ep 考虑了 ul 导致的耦合项, M 和 K 分别代表广义质量和广义刚度。式中对积分项进行了无量纲化,
= x/l1, A0 和 I0 分别代表 x = 0 处的截面积 A(0)和截面惯性矩 I(0)。由式(15)和(16)可知, mς ξ = m ξς ,kς ξ = kξς, 即矩阵 M 与 K 对称。
由式(16)可知, 当轴向力为压力(即 T<0)时, 轴向力将降低 K 的特征值, 即降低系统的刚度。K 的最小特征值小于 0 代表存在一个特征向量, 使得随着该向量的模增大, 系统的势能 Ep 降低, 即着陆腿发生轴压失稳。因此, 若轴向力 T < 0, 则假设着陆腿横向振动过程中 T 始终存在且恒定; 若 T > 0, 则假设着陆腿横向振动过程中 T = 0。
分布力 q(x)在虚位移 δ 上的虚功为
将 Ek 和 Ep 和 δW 代入 Lagrange 方程, 可得 η 表示的n 自由度振动方程:
式中, pη为广义力, 表达式为
可见广义力 pη 与 q(x)线性相关, 令 pη,c 表示冲击激励 qc(x)对应的广义力, pη,1 表示稳态激励 qc(x)对应的广义力, 则有 pη=pη,c+pη,1。pη,c 和 pη,1的表达式为
其中, 
和
为不显含
, 
a 和 g 的向量, 各阶的表达式为
当式(18)的右端项为 0 时, 其通解为广义特征值问题(K - ω2M)η0=0。K 满足轴压稳定时, 该问题有 n个正实数特征值 ω12, L, ωn2 及其对应的单位特征向量 z1, L, zn。令 Dω=diag{ω12, L, ωn2}, Z={z1, L,zn}, 则 M -1KZ=ZDω。进一步, 令 Z(t)=Z-1η(t), 并代入式(18), 可得 v˙˙+Dωv=(MZ)-1pη, υ 与 η 同为 n 维广义坐标。令向量 pv(0)=(MZ)-1pη(0), pv(1)=(MZ)-1pη(1),则式(18)可解耦为
着陆前着陆腿相对箭体处于静止状态, 故初始状态 v=0, 
˙=0。则由式(22)可得 v 在碰撞后满足
, 其中 Ac,ξ 和 AI,ξ分别表示第ξ阶模态由冲击激励和稳态激励产生的振幅:
则着陆腿横向振动的位移y(x, t)为
对于薄壁圆管, 其弯矩 M(x,t)和弯曲应力的大小
分别为
式中, E 为着陆腿材料的弹性模量; r(x)和 I(x)分别为着陆腿截面的半径和惯性矩, 考虑着陆腿的变截面特征后为截面轴向位置x的函数。
2.3 着陆腿的最大应力
我们考虑轴向力和横向弯曲, 给出着陆腿的最大应力计算方法。
主腿需满足轴压稳定, 即广义刚度矩阵 K 正定。主腿横向振动弯曲应力
由式(25)计算, 式中轴向力 T=T1, T1 由蜂窝压溃强度确定。主腿最大应力的绝对值 σ 1, m ax为
副腿的最大应力应考虑轴向力为正和为负两种情况。首先由静力平衡条件(式(6)和(8)), 得到副腿的最大轴向力 T2max 和最小轴向力 T2min。针对轴向力为正的情况, 由于此时轴向力可降低振动弯矩,故计算横向振动时轴向力取 0, 由式(25)得横向振动的最大应力 σ 2,M,max|T=0。若 T 2min＜ 0, 应再考虑轴向力为负的情况, 再将轴向力取 T 2m in, 由式(25)得
则副腿的最大应力绝对值 σ 2 ,m ax 为
3.1 单腿着陆冲击实验
采用单套着陆腿的着陆冲击实验来验证本文提出的计算模型, 实验装置的示意图见图7。采用重量 M = 5 吨的模拟箭体, 模拟箭体沿刚性壁面沿高度方向滑动。着陆腿的主腿采用 4 节式伸缩结构,在第 4 节腿内装有铝蜂窝缓冲器, 其直径为 318 mm, 压溃强度为 5 MPa, 则主腿轴力 T1= -397 kN。主腿和副腿的材料均为钢材, 材料的模量 E=206 GPa, 密度 ρ=7800 kg/m3。着陆腿上装有整流罩, 与副腿和足垫连接。由于摩擦系数μ与碰撞过程界面的相对运动有关[29–30], 故针对本文的研究对象难以精确地建模。为降低摩擦建模精度对实验结果的影响, 在足垫与“地面”的接触点下设置滑轨, 近似地取摩擦系数μ=0。
图7 着陆冲击实验装置示意图
Fig.7 Device of landing leg impact experiment
着陆腿几何参数: 跨距 xB=8.25 m, 支撑高度 h=2.8 m, 副腿、主腿安装点与箭体轴线的距离 rL1=rL2=1.9 m, 主腿与箭体的夹角 β=47.5°, 着陆腿壁厚为 t=5 mm。副腿为等截面圆管, 直径为 r2=90 mm。将上述几何参数代入几何约束式(1), 得副腿的几何参数α= 67.9°, γ = 10.2°, l2=7.573 m。主腿采用如图5所示的等效模型, 总长度 l1=8.614 m, 各段长度l1i、半径 r1i、质量块位置 x1i 和质量块质量 m1i 如表1 所示, 其中 x14 和 m14 为缓冲器。
表1 主腿等效模型的参数
Table 1 Parameters of equivalent model of the mean leg
编号 l1i / m r1i / m x1i / m m1i / kg 1 2.206 207.5 2.206 107.8 2 2.020 195 4.226 99.4 3 2.020 182.5 6.246 91.6 4 2.368 170 8.008 110.4
采用高速摄像对触地速度进行测量, 采样频率为 1 kHz。采用动态应变仪测量着陆过程中主腿的应变, 采样频率为 10 kHz。在主腿的第 1 节、第 2节、第 3 节和第 4 节各布置一个测点, 测量轴向应变, 如图7 所示。每个测点中对称地布置 4 片轴线方向的应变片, 在图5 所示的主腿局部坐标系中,4个测点的坐标分别为 1.381, 3.401, 5.421 和 6.753 m, 其中测点 4 位于缓冲器与主腿连接处的下方。进行两次着陆冲击实验, 通过调整下落高度来调整触地 速度。
两次实验中, 由高速摄像测得的触地速度如图8 所示, 图中时间零点为启动摄像机时刻。由图8可见, 实验 1 中触地速度为 1.490 m/s, 实验 2 中触地速度为 2.658 m/s。箭体减速下降过程中, 由于刚性壁面与模拟箭体间存在摩擦力, 导致难以通过支撑力准确地计算箭体加速度, 由式(6)可得支撑力为N =154.9 kN。由图8 可见实验 1 中加速度约为 83.1 m/s2, 实验 2 中为 167.6 m/s2, 说明摩擦力的影响较大。因此, 计算着陆腿的动力学响应时, 实验 1 和实验 2 中加速度 a 分别取 83.1 和 167.6 m/s2, 并在速度降为 0 后取 a = 0, 针对当前各阶模态的位置和速度, 重新计算模态的振幅和相位。
图8 触地速度测量结果
Fig.8 Measurement of the touchdown velocity
3.2 理论模型与实验结果的对比分析
首先检验计算角加速度时由式(3)到式(4)的近似是否合理。式(3)中取速度为触地速度 v0 = 2.658 m/s, α˙˙中 省 略 的 项 cosα/sin3α·v02/l22=0.057 rad/s2,该值显著小于 α/l2sinα=11.8 rad/s2; β˙˙中省略项的值为 2.3×10-4 rad/s2, 未省略项的值为 9.76 rad/s2。被省略的项均远小于未被省略的项, 说明这种近似是合理的。
由广义速度表示的冲击响应式(23)和角速度的几何约束式(2)可得, 主腿轴压力确定时, 各阶模态的响应 Ac,ξ 与触地速度成正比, AI,ξ 与箭体加速度线性相关。取 20 阶模态坐标, 采用 2.2 节的方法计算主腿前 6 阶模态的频率及在触地速度 1 m/s 下各阶模态的响应, 结果见表2, 其中频率考虑了轴向压力的耦合。由理论值可见 1 阶模态的响应显著高于其他高阶模态, 即振动弯矩主要由 1 阶模态引起。由于 6 阶频率已达到 1 阶频率的 42.7 倍, 在阻尼作用下衰减速度更快, 且响应 Ac,6 和 AI,6 均远小于一阶模态的响应 Ac,1 和 AI,1, 故由式(25)计算弯矩和应力时仅考虑前 6 阶模态。
表2 主腿各阶模态频率和单位响应
Table 2 Frequency and response for each mode
阶数 fξ/Hz Ac,ξ/mm AI,ξ/mm 1 9.07 –11.7 –2.45 2 38.26 –1.44 0.053 3 86.33 0.448 9.50×10–3 4 157.5 –0.295 8.05×10–4 5 272.4 –0.0358 –14.5×10–4 6 387.5 –0.0290 5.02×10–4
将测点处的应变数据转化为测点处的弯矩和轴向力, 分别对比实验结果与理论模型计算结果, 对本文提出的计算模型进行验证。由于轴向振动受节间锁定机构的影响较大, 2.1 节中未对轴向冲击的动力学响应进行精确的建模, 故本文在轴向仅对比最大轴向压力, 结果见表3。可见最大轴向力绝对值均未超过|2T1| = 794 kN, 表明 2.1 节的假设是合理的。由于测点 4 位于缓冲器与主腿连接处的下方,故测点 4 处的轴向力较小。
表3 主腿各测点的最大轴向力
Table 3 Maximum value of the axle force of each measured section of the main leg
实验 最大轴向力/ kN测点1 测点2 测点3 测点4实验1 –389.8 –406.2 –393.4 –104.7实验2 –437.1 –477.2 –532.1 –158.7
由于考虑高阶模态时阻尼的影响不可忽略, 本文将阻尼考虑为模态阻尼比。在金属材料伸缩式机构的振动问题中, 模态阻尼比取值通常为 0.01~0.09[31–33], 本文取各阶模态阻尼比均为 0.09。通过对实验测量结果进行傅里叶变换, 计算动力学响应的频率, 得到实验 1 中前 3 阶的频率分别为 4.78,19.91 和 44.59 Hz, 实验 2 中前 3 阶频率分别为 4.51,
20.36 和 44.25 Hz。实验测得的频率低于表2 中频率的计算值(约为计算值的 0.523 倍), 原因是实际着陆腿为多节结构, 节间连接降低了着陆腿的横向振动的刚度, 进而降低了频率。4个测点处弯矩的对比结果如图9 所示。为修正频率偏差的影响, 做弯矩–时间曲线时对“计算值”的时间轴进行修正, 修正系数为 1/0.523。图9 中时间零点为碰撞发生时刻, 时间范围为碰撞前 60 ms 至碰撞后 80 ms。最大弯矩的测量值与计算值见表4。
表4 主腿各测点的最大弯矩(kN·m)
Table 4 Maximum moment of the measured section (kN·m)
测点 实验1 实验2测量值 计算值 测量值 计算值1 45.2 44.3 71.9 79.4 2 65.5 68.6 104.7 118.0 3 64.8 68.7 105.2 116.9 4 59.4 53.7 89.8 92.7
图9 实验测量结果与等效模型计算结果
Fig.9 Results of experiments and equivalent model analysis
由图9 可见, 在 4个测量截面, 实验结果与本文方法得到的弯矩变化曲线的趋势较为接近, 但在碰撞开始时刻及弯矩变化较为剧烈的时刻存在一定的差异。表4 显示, 最大弯矩计算值的偏差为 13%,说明本文提出的方法可较精确地计算低阶振动, 但对高阶振动的计算以及对阻尼比的估计存在一定的误差。综合上述分析可知, 计算结果得到的最大弯矩与实验结果相近, 该半解析的理论模型可较准确地计算着陆冲击引起的着陆腿横向振动。
4 着陆缓冲机构参数影响分析
火箭着陆腿的跨距 xB 由着陆稳定决定[25], 展开后支撑高度 h 由发动机羽流的流场设计决定[34], 主腿长度和着陆腿的半径为可设计的参数。利用本文提出的理论模型, 通过强度要求计算着陆腿半径,进而确定着陆腿质量, 分析最优主腿长度。由几何约束式(1)可知, 主腿长度 l1 可由 β 唯一确定, 故本文分析最优角 β。
火箭一级着陆体的质量 MB=20 t, 箭体直径为3.35 m。采用 4 套着陆腿, 取 xB=8.25 m, h=2.8 m, rL1=rL2 = 1.9 m, 则α= 67.9°, 上述参数与 3.1 节实验中采用的着陆腿相同。着陆腿材料采用 LS-3k 0°/90°碳纤维复合材料[35], 模量 E = 51.4 GPa, 密度 ρ=1600 kg/m3, 抗拉强度为 393.3 MPa, 许可应力取抗拉强度的 70%, 即[σ]= 275.3 MPa, 着陆腿壁厚 t = 6 mm, 足垫与地面摩擦系数上限 μmax= 0.35, 满足μmax
取触地速度 v = 3 m/s, 容许加速度 amax = 5 g, 分析着陆腿质量 ML, m1, m2 与 β 的关系, 结果如图10所示。可见, 着陆腿质量 ML关于角 β存在极小值,该构型如下: 主腿与箭体夹角 βopt = 41.6°, 第一节主腿半径 r11opt = 210.9 mm, 副腿半径 r2opt=100.4 mm,主腿设计轴力为 T1opt = -527.1 kN, 对应 ML 的极小值为 MLopt=458.6 kg。分析最优构型中主腿的轴压稳定性, 令 λk,1(T1)表示主腿刚度矩阵 K 的最小特征值, 最优构型下 λk,1(T1) = 0.272λk,1(0), 临界轴向力T1,cri = -724 kN, 可见此时着陆腿主腿半径由轴压稳定性和横向振动共同决定。
图10 着陆腿质量与夹角β的关系
Fig.10 Leg weight under different angle β
提升材料的许可应力[σ]可以降低着陆腿质量,且主腿与箭体的最优夹角 βopt 也发生改变。下面分析最优构型的着陆腿质量 MLopt, βopt 与[σ]的关系。假设仅改变[σ], 材料的模量 E 和密度 ρ 保持不变如图11 所示, 随[σ]提高, MLopt 下降, βopt 增大。当[σ]<490 MPa 时, 随[σ]提高, MLopt 迅速下降(图11 中虚线左侧); 当[σ]>490 MPa 时, 随[σ]提高, MLopt和 βopt 的变化均较为缓慢(图11 中虚线右侧)。当[σ]= 490 MPa 和[σ]=800 MPa 时, 最优构型的主腿夹角 βopt、主腿质量 m1opt、单条副腿质量 m2opt、主腿轴向压力 T1 和主腿轴压失稳的临界轴向压力 T1,cri如表5 所示, 可见许可应力由 490 MPa 提升至800 MPa 后, 主腿质量和主腿与箭体夹角变化均较小,副腿质量显著降低。原因在于, 当强度较小时, 主腿截面半径主要由强度决定, 提升[σ]可降低主腿半径, 增大角 βopt(即减小主腿长度); 随主腿强度提高,主腿截面半径转而由轴压稳定决定, 此时降低主腿长度将提高主腿轴向压力 T1, 进而增大主腿的截面半径, 主腿质量变化不明显。由于副腿的截面半径主要由轴向拉力决定, 故[σ]>490 MPa 时, 提高材料强度仅可降低副腿截面半径, 进而降低副腿质量。由于着陆腿系统的质量主要来自主腿, 故此时提高材料强度总质量变化较小。
图11 着陆腿质量、主腿最优夹角与材料强度曲线
Fig.11 Leg weight and optimal angle under different strength
表5 典型材料强度下的最优构型参数
Table 5 Configuration parameters of typical strength
[σ]/MPa βopt/(°) m1opt/kg M2opt/kg T1/kN T1,cri /kN 490 48.9 321.6 28.5 –717 –758 800 50.4 319.5 17.8 –776 –817
5 结论
针对垂直起降重复使用火箭采用的倒三角式着陆腿, 为在着陆腿构型设计阶段准确、快速地计算着陆冲击作用下的振动载荷, 本文提出一种半解析计算方法。在分析着陆冲击导致的横向振动时, 将横向激励等效为足垫与地面碰撞产生的冲击激励和缓冲后刚体运动加速度产生的稳态激励。进一步地, 采用假设模态法建立着陆腿横向振动模态的半解析解, 采用模态叠加法完成着陆冲击作用下系统响应的求解, 并在模型中考虑了轴向力与横向振动的耦合及着陆腿的截面特征。单套着陆腿着陆冲击实验结果表明, 本文提出的方法可较精确地计算着陆过程中着陆腿的弯矩。最后, 基于提出的冲击载荷近似计算方法, 研究了着陆腿构型参数对着陆腿设计的影响, 以着陆腿系统的总质量最小为目标,分析了最优构型, 并进一步讨论了最优构型随着陆腿材料强度的变化规律。
[1]王辰, 王小军, 张宏剑, 等.可重复使用运载火箭发展研究.飞航导弹, 2018(9): 18–26
[2]牟宇, 孙冀伟, 秦旭东.猎鹰 9 火箭 block 5 构型首次飞行任务解析.宇航总体技术, 2018, 2(5): 1–7
[3]宋征宇, 吴义田, 徐珊姝, 等.CZ-8: 长征火箭系列商业化与智慧化的先行者.深空探测学报.2021,8(1): 1–14
[4]刘志全, 黄传平.月球探测器软着陆机构发展综述.中国空间科学技术, 2006(1): 33–39
[5]王闯, 邓宗全, 高海波, 等.国内外月球着陆器研究状况.导弹与航天运载技术, 2006(4): 31–36
[6]Seiff A.Post-Viking models for the structure of the summer atmosphere of Mars.Advances in Space Research, 1982, 2(2): 3–17
[7]梁绍敏, 王永滨, 王立武, 等.月球着陆器着陆过程的 DEM-FEM 耦合分析.固体力学学报, 2019,40(1): 39–50
[8]万峻麟, 聂宏, 李立春, 等.月球着陆器着陆性能及多因素影响分析.南京航空航天大学学报, 2010,42(3): 288–293
[9]Yue Shuai, Nie Hong, Zhang Ming, et al.Design and landing dynamic analysis of reusable landing leg for a near-space manned capsule.Acta Astronautica, 2018,147: 9–26
[10]雷波, 张明, 岳帅.可重复使用运载器的耐坠毁缓冲装置的设计优化.宇航学报, 2019, 40(9): 996–1005
[11]Li Meng, Deng Zongquan, Guo Hongwei, et al.Optimizing crashworthiness design of square honeycomb structure.Journal of Central South University, 2014,21(3): 912–919
[12]王家俊, 王春洁, 宋顺广.基于响应面法的月球着陆器软着陆性能优化.北京航空航天大学学报,2014, 40(5): 707–711
[13]吴宏宇, 王春洁, 丁宗茂, 等.着陆姿态不确定下的着陆器缓冲机构优化设计.宇航学报, 2018, 39(12): 1323–1331
[14]刘学翱, 吴宏宇, 王春洁, 等.着陆器变阻尼缓冲器性能分析及参数优化.北京航空航天大学学报,2018, 44(10): 2149–2155
[15]Low K H.On the methods to derive frequency equations of beams carrying multiple masses.International Journal of Mechanical Sciences, 2001, 43: 871–881
[16]Low K H.Frequencies of beams carrying multiple masses: Rayleigh estimation versus eigenanalysis solutions.Journal of Sound and Vibration, 2003, 268:843–853
[17]周叮.一类变截面梁横向自由振动的精确解析解.振动与冲击, 1996, 15(3): 12–15
[18]Chen D W, Wu J S.The exact solution for the natural frequencies and mode shapes of non-uniform beams with multiple spring-mass system.Journal of Sound and Vibration, 2002, 255(2): 299–322
[19]Jang S K, Bert C W.Free vibration of stepped beams:higher mode frequencies and effects of steps on frequency.Journal of Sound and Vibration, 1989,132(1): 164–168
[20]Cheng Jianlian, Xu Hui, Yan Anzhi.Frequency analysis of a rotating cantilever beam using assumed mode method with coupling effect.Mechanics Based Design of Structures and Machines, 2006, 34(1): 25–47
[21]潘旦光, 付相球, 韦杉杉, 等.变高度悬臂箱梁剪力滞效应的半解析解.工程力学, 2018, 35(9): 207–213
[22]Kong Xiangxi, Jiao Jiang, Zhou Chong, et al.Sommerfeld effect and synchronization analysis in a simply supported beam system excited by two nonideal induction motors.Nonlinear Dynamics, 2020,100(3): 2047–2070
[23]Luo Jun, Zhu Shengyang, Zhai Wanming.Formulation of curved beam vibrations and its extended application to train-track spatial interactions.Mechanical Systems and Signal Processing, 2022, 165: 108393
[24]Liu Caishan, Zhao Zhen, Brogliato B.Frictionless multiple impacts in multibody systems.I.Theoretical framework.Proc R Soc A, 2008, 464: 3193–3211
[25]袁晗, 王小军, 张宏剑, 等.重复使用火箭着陆结构稳定性分析.力学学报, 2020, 52(4): 1007–1023
[26]陈金宝, 钱佳程, 贾山, 等.新型多级铝蜂窝缓冲器的缓冲性能研究与分析.西北工业大学学报,2020, 38(S1): 1–6
[27]Yang Xianfeng, Sun Yuxin, Yang Jialing, et al.Outof-plane crashworthiness analysis of bio-inspired aluminum honeycomb patterned with horseshoe mesostructure.Thin-Walled Structures, 2018, 125: 1–11
[28]范纪华, 陈立威, 王明强, 等.旋转中心刚体-FGM梁刚柔热耦合动力学特性研究.力学学报, 2019,51(6): 1905–1917
[29]鲁建东, 赵振.振动引起的界面切向相对运动对摩擦的影响.北京大学学报(自然科学版), 2020, 56(5):777–784
[30]王辰, 张宏剑, 王小军, 等.自旋摩阻动力学模型与实验研究.北京大学学报(自然科学版), 2018,54(5): 915–920
[31]李青, 任德鹏, 王闯, 等.月球探测器力学环境研究及试验条件制定.航天器工程, 2018, 27(1): 137–142
[32]安洋, 林宝军, 刘佳伟.框架面板式构型卫星冲击响应特性分析方法.中国科学: G 辑, 2021, 51(1):15–25
[33]朱冬梅, 范占贝, 刘海平, 等.一种周期结构隔振器力学性能研究.湖南大学学报(自然科学版),2020, 47(4): 32–39
[34]严立, 王平阳, 欧阳华.月面环境发动机羽流冲击力效应模拟计算.上海交通大学学报(自然版),2012, 46(8): 1310–1314
[35]Xiong Jian, Feng Lina, Ghosh R, et al.Fabrication and mechanical behavior of carbon fiber composite sandwich cylindrical shells with corrugated cores.Composite Structures, 2015, 156: 307–319