基于CFD/CSD耦合的火箭跨声速气动阻尼特性分析

李泳德, 郭 力, 季 辰

(中国航天空气动力技术研究院, 北京 100074)

引 言

通常情况下人们认为气动力对火箭的振动起到阻尼作用, 即气动阻尼为正值。然而随着大推力火箭发展, 火箭的长细比逐渐加大, 导致弯曲刚度越来越小, 同时为了满足有效载荷的外形要求, 火箭头部整流罩尺寸不断加大, 后续箱体的直径却保持不变, 形成了典型的锤头体外形。国内外大量的火箭研制经验表明[1-9], 对于此类锤头体外形火箭的气动设计, 必须要进行动态气动载荷与动态气弹稳定性分析, 否则设计的疏忽可能会导致火箭结构出现毁灭性的破坏进而导致发射失败。

目前常用的衡量气弹稳定性的方法是通过风洞试验来获取气动阻尼系数。早在1963年, 美国国家航空航天局Ames研究中心(NASA Ames Research Center)采用半刚性模型开展试验研究[10], 获取火箭头部的气动阻尼来评估其稳定性, 但这只能用来模拟火箭弯曲振型前节点之前部分的结构动力学特性。直到兰利研究中心(NASA Langley Research Center)开发了全弹性模型气动阻尼试验技术, 其可以模拟整体的结构动力学特性以及气动外形, 并应用于多款运载火箭研制[11-15]。国内, 中国航天空气动力技术研究院对气动阻尼问题开展过较多的研究[16-20], 从模型设计方法、 模型制作工艺、 试验机构设计和数据处理等诸多方面, 逐步改进实现了从半刚性模型到全弹性模型的过渡, 并在多个型号上得到验证。

然而通过风洞试验研究气动弹性问题, 技术难度大, 试验成本高, 同时几乎不可能开展全尺寸试验。因此通过数值计算的方法开展相关研究是另一种重要的手段。刘子强等[21]实现了通过数值计算确定气动阻尼系数的技术和方法, 并与试验结果进行对比, 证实了该方法的可靠性。冉景洪等[22]通过模态数据结合准定常理论的方法分析了减阻杆加后体这一弹性结构的气动阻尼, 结果表明减阻杆造成的分离流会对后体的气动阻尼系数产生影响。朱剑等[23]针对新一代捆绑式运载火箭发展了非结构网格下的气动阻尼计算方法, 并分析了攻角、 Mach数等参数对气动阻尼的影响。

本文在之前的计算方法[23]的基础上采用IDDES模型, 考虑脉动压力的影响, 通过强迫振动的方式, 针对捆绑式运载火箭的某一特定模态进行数值计算仿真, 研究前节点位置, 振动振幅, 脉动压力等参数对气动阻尼的影响规律。

1 计算方法

图1为本文所用的捆绑式运载火箭的计算模型, 是典型的锤头体结构。在跨声速阶段, 其头部会产生激波造成激波边界层干扰, 而在锤头体外形的过渡段会出现气流分离。为探究各部分气动阻尼的变化, 将整个箭体分为头部、 过渡段、 弹身3个部分。

图1 表面网格及区域划分Fig. 1 Surface grid and region division

1.1 流场仿真模型

本文分别用Reynolds平均法(Reynolds-averaged Navier-Stokes, RANS)和改进的延迟分离涡模拟(improved delayed detached-eddy simulation, IDDES)[24-25]进行计算, 在RANS方程中, 将变量分为平均值和波动值两部分, 对于速度分量有

(1)

(2)

其中,k和ω分别代表湍流动能和湍流耗散率,Γk和Γω分别代表k和ω的有效扩散系数,Gk和Gω分别代表k和ω的生成率,Yk和Yω分别代表k和ω的耗散率。因此RANS方法只能计算大尺度的平均流动, 本文采用IDDES方法计算脉动压力对气动阻尼的影响。

IDDES方法是由分离涡模拟(detached-eddy simulation, DES)方法改进而来, 其本质思想与DES方法相同, 是想以网格尺度和模型中的特征尺度隐式划分RANS和大涡模拟(large-eddy simulation, LES)区域, 使其既能处理RANS方法无法得到的脉动场, 也能降低LES方法在模拟高Reynolds数流动时所需的计算资源。区别在于当边界层较厚或者分离区域较窄时, DES方法会出现如模型应力损耗(modeled stress depletion, MSD), 网格诱导分离(grid-induced separation, GIS)以及对数层不匹配(logarithmic-layer mismatch, LLM)问题[24], 而IDDES模型通过改良计算区域划分, 结合延迟分离涡模拟(delayed detached-eddy simulation, DDES)和壁面模型大涡模拟(wall-modeled large-eddy simulation, WMLES), 定义新的长度尺度解决了这些问题, 具体公式详见文献[25]。

流场网格如图2、 图3所示, 边界层采用棱柱层结构, 并调整第1层网格高度使得y+小于1, 远场部分采用六面体结构网格, 与边界层的过渡层采用非结构网格。整体网格单元数量为4.2×106。

图2 y方向截面网格示意图Fig. 2 Schematic diagram of cross-sectional grid in the y-direction

图3 x方向截面网格示意图Fig. 3 Schematic diagram of cross-sectional grid in the x-direction

物面边界条件为无滑移壁面条件, 远场采用压力远场边界条件, 湍流模型采用SSTk-ω模型, 采用密度基求解, 气体黏性采用Sutherland定律, 空间离散采用2阶迎风格式, 对流通量采用Roe格式。

1.2 结构分析模型

结构与流场耦合分析过程中, 结构部分可以采用模态方法描述。结构模态可以通过有限元方法与结构模态试验方法获得。本文采用有限元分析结果获得的模态, 图4所示为结构的前3阶模态, 本文只分析计算结果中气动阻尼最小的第2阶模态。

(a) f=1.200 Hz

由于火箭结构外形简单, 一般不考虑其扭转影响, 因此可以将其简化为简单的梁模型, 这样就可以给出其模态振动方程

(3)

式中,qi为第i阶模态的广义位移,bi为第i阶模态的结构阻尼系数,ωi为第i阶模态的固有频率,fi为第i阶模态下质量归一化的广义气动力。若将fi按照Taylor展开并略去高阶项, 可以将其转化为气动阻尼项与气动刚度项的形式, 则式(3)可写为

(4)

式中,Bi为气动阻尼系数,Ki为气动刚度系数, 研究表明[26], 气动刚度相对于结构刚度为小量可以忽略不计, 而在计算中结构阻尼往往设置为0, 因此气动阻尼可以直接反映其气弹稳定性。

1.3 气动阻尼分析原理

气动阻尼的分析可以采用强迫振动或者自由振动的方式进行, 这两种方法获得的时域数据不同, 提取气动阻尼的方式也不同。强迫振动方法初始演化过程较短, 因此计算量较小, 同时能够分析某一种振动形式的气动阻尼, 明确该振动形式是收敛还是发散。分析过程中能够获得不同部位与部件的气动阻尼。但是对于多模态相互作用引起的发散(例如颤振)较难预测。自由振动方法需要一定的自由演化时间才能够对时域数据进行分析, 不过自由振动方法能够获得最能够吸收能量的模态及其振动频率。

对于本研究所关注的问题, 气动载荷对结构振动的过程中气动阻尼的影响较大, 而对气动刚度与气动质量影响较小, 即结构的固有振动频率受到来流的影响较小, 其稳定性问题主要由气动阻尼的正、 负引起, 所以采用强迫振动方法分析。

强迫振动下结构做简谐模态振动

qi(t)=Asin(ωit)

式中,A表示振动的振幅, 将其代入计算气动力的公式中[21]并做正交积分可得

(5)

式中,Mi为第i阶模态的模态质量,T为整数倍周期,G为广义气动力。根据式(5)便可以得到局部或分区域的气动阻尼。

1.4 耦合计算流程

首先进行模态分析, 以确定结构的模态频率与振型, 用以设计强迫振动的频率和振幅。非定常流场计算前先进行定常流场计算, 来加快非定常计算的演化速度并增强收敛性, 结构节点位移通过径向基函数(RBF)插值方法[27]映射到气动网格节点上, 来进行网格的变形, 这里径向基函数选用Wendland C2, 如下所示

最后将计算出来的广义力提取出来, 截取演化完毕的整数倍周期, 进行气动阻尼计算。耦合计算流程图如图5所示。

图5 耦合计算流程图Fig. 5 Flow chart of coupled calculation

2 结果分析与讨论

2.1 流场分析结果

计算的来流Mach数范围为0.7～1.2。其中中截面的压力分布如图6所示。可以看出在头部出现了膨胀波以及跨声速激波, 在过渡段存在流动分离, 随着Mach数的增大, 头部低压区域逐渐扩张, 并且能明显看到, 在流动再附的位置产生了再附激波。

(a) Ma=0.70

2.2 气动阻尼分布

通过上述流场分析, 可以看出火箭不同部位流动结构并不相同, 在头部与箭身上, 流动主要为附着流动, 而在过渡段会出现较为复杂的波系结构以及流动分离。针对不同的流动结构随流向站位x的变化, 设该位置上广义力与广义位移的相位差为φ(x), 并且简谐振动没有引入其他模态的广义力, 则广义力的表达式为

G(x,t)=Fgen·sin[ωt+φ(x)]+F0

(6)

其中,Fgen为广义力的振动幅度,F0为广义力的常数偏移量。将式(6)代入到式(5)中得到

其中, 广义力的常数偏移量F0的积分为0, 因此省略。通过将等式中的正弦函数部分进行和差化积得到

(7)

式(7)中第1部分在整个周期中的积分为0, 只有第2部分保留, 因此得到

(8)

式(8)中积分部分恒为正值, 决定整个气动阻尼的部分只有相位角φ(x)的正弦值sin[φ(x)], 为了能够更加直观地获得相位角与气动阻尼B之间的关系, 须将符号转化为对应的正弦函数转角, 根据正弦关系, 此转角为π, 因此得到

(9)

图7为气动阻尼变化曲线, 可以看出随着Mach数的增大, 整体气动阻尼先增大后减少, 在Mach数为0.98时达到最大值, 过渡段与箭体的气动阻尼变化趋势与整体基本相同, 而头部区域则不同, 是随着Mach数的增大一直增大, 只是增长速率变缓。

图7 有助推时气动阻尼变化曲线Fig. 7 Aerodynamic damping change curve with boost

根据式(9), 得到相位角与气动阻尼B之间的关系为: 当φ(x)∈(-π, 0)时, 相位角滞后, 气动阻尼B为负值; 当φ(x)∈(0, π), 相位角提前, 气动阻尼B为正值; 为当φ(x)=0时, 无相位角差别, 气动阻尼B为0。

在过渡段上, 复杂的波系结构以及流动分离, 使得气动力与结构位移之间会出现较为明显的迟滞现象, 从而导致相位角φ(x)∈(-π, 0), 由此在过渡段上产生了负的气动阻尼。

计算过程中的广义力与广义位移随时间变化曲线如图8所示, 可以看出所有工况计算结果都表现良好, 需要注意的是在非定常计算初期, 演化的不完全导致广义力存在一些突变异常的结果, 计算气动阻尼时须剔除, 选择后面演化完全的周期。本文计算了9个周期, 剔除了第1个周期出现的错误结果, 采用后8个周期进行气动阻尼分析。强迫运动振幅为芯级直径的0.5%。

(a) Ma=0.70

2.3 气动阻尼影响因素

2.3.1 有无助推对气动阻尼的影响

捆绑式运载火箭相比于传统的运载火箭, 最大的区别就是在尾部四周捆绑了助推器, 使得其流场特性变得复杂, 因此须分析其对气动阻尼的影响。

图7、 图9分别为有无助推时气动阻尼变化曲线, 可以看出随着Mach数的增大整体气动阻尼先增大后减少, 在Mach数为0.98时达到最大值, 过渡段与箭体的气动阻尼变化趋势与整体基本相同, 而头部区域则不同, 是随着Mach数的增大一直增大, 只是增长速率变缓。对比两个图可知, 助推主要起增大气动阻尼的作用。还可以看出有无助推情况下头部的气动阻尼变化很小, 意味着在箭体尾部施加控制很难影响到头部的气动阻尼, 特别是在超声速流场中。

图9 无助推时气动阻尼变化曲线Fig. 9 Aerodynamic damping change curve without boost

2.3.2 前节点位置影响

为了考察前节点位置变化对气动阻尼的影响, 在保持振动频率不变、 头部最大振型位置与振幅不变的条件下移动前节点, 变化后的振型如图10所示。

(a) Front node after the transition region

根据对计算结果的分析分别获得了不同前节点位置的整体气动阻尼对比与过渡段气动阻尼对比, 如图11、 图12所示, 可以看出前节点位置的改变并没有影响整体气动阻尼随Mach数增大而增大的趋势, 且前节点在过渡段上与过渡段前的整体气动阻尼相差不大, 而前节点在过渡段后的整体气动阻尼要高于另两种情况, 因此过渡段与头部放在同一侧有助于提高气动阻尼。过渡段的气动阻尼会随着前节点的变化发生剧烈改变, 前节点在过渡段前后随Mach数增大的变化规律相反, 节点前后的振动相位变化导致不同节点位置过渡段的振动相位不同, 进而导致气动阻尼发生变化。

图11 不同节点位置的整体气动阻尼Fig. 11 Overall aerodynamic damping at different node positions

图12 不同节点位置的过渡段气动阻尼Fig. 12 Aerodynamic damping of the transition region at different node positions

2.3.3 强迫振动振幅大小对气动阻尼的影响

为了考察强迫振动振幅大小对气动阻尼的影响, 在保证流场结构不发生改变的前提下, 振动振幅分别为原来的一半和两倍, 根据工程经验, 如果振幅超过芯级直径的5%, 则须考虑流场结构改变所造成的影响。图13、 图14分别为不同振幅下的整体与头部气动阻尼。

图13 不同振幅下整体气动阻尼Fig. 13 Overall aerodynamic damping at different amplitudes

图14 不同振幅下头部气动阻尼Fig. 14 Aerodynamic damping of the head region at different amplitudes

可以发现改变振幅无论是对整体气动阻尼还是头部气动阻尼来说变化都很小, 这意味着气动阻尼的大小主要取决于气动力与结构振动的相位差, 不依赖于振动幅度的大小。

2.3.4 脉动压力对气动阻尼的影响

为了模拟出脉动压力的影响, 采用IDDES方法对火箭气动阻尼进行计算, 计算来流Mach数为0.92, 计算过程中的广义力与广义位移如图15所示, 相较于图8可以看出广义力随时间变化曲线并不光滑, 脉动压力的存在导致广义力由多个频率叠加而成。

图15 基于IDDES的广义力与广义位移变化曲线Fig. 15 Variation cures of generalized force and generalized displacement based on IDDES

由于第2阶模态的频率为2.46 Hz, 而由分离流、 激波振荡等引起的脉动压力频率往往远大于此频率, 因此这里选择3.5 Hz为分界, 将高于3.5 Hz的部分视为由抖振脉动压力引起的广义力, 低于3.5 Hz的部分视为强迫振动引起的广义力, 通过低通滤波把高于3.5 Hz的广义力滤掉, 可以获得由强迫振动引起的广义力与广义位移变化曲线, 如图16所示, 通过此广义力计算的气动阻尼为2.08‰。同样地, 进行高通滤波将低于3.5 Hz的广义力滤掉, 可以获得由抖振脉动压力引起的气动阻尼为(2.94×10-3)‰, 由此得到脉动压力引起的气动阻尼变化为0.14%, 可以忽略不计。同时使用RANS方法计算的气动阻尼为2.07‰, 与IDDES的计算结果相比误差约为(2.94×10-3+2.08-2.07)/2.07≈0.48%, 这说明针对气动阻尼的模拟, 抖振引起的脉动压力对气动阻尼的计算结果影响很小, 起主要作用的还是广义力的变化, 该变化由强迫振动引起的结构边界变化所导致。

图16 滤波后的广义力与广义位移变化曲线Fig. 16 Variation cures of generalized force and generalized displacement variation curve after filtering

3 结论

本文通过数值计算方法研究了火箭的气动阻尼特性。根据流动特征分析与理论推导, 发现火箭过渡段几何外形的收缩导致该区域出现复杂的分离与激波结构, 从而造成了气动力相对于结构振动相位的滞后, 导致了该区域为气动负阻尼, 即气动不稳定性的主要来源。

在此机理的基础上, 分析了前节点位置、 振动振幅、 脉动压力等因素对气动阻尼的影响规律。可以得出以下结论:

1) 助推增加了正阻尼区域的面积, 从而相对于没有助推的构型起到了增加气动阻尼的作用。

2) 前节点位置的改变对过渡段气动阻尼影响很大, 节点前后的振动方向相反, 导致节点在过渡段前后的气动阻尼变化规律也截然相反, 将过渡段与头部区域放在节点的同一侧有助于增加气动阻尼。

3) 在不改变流场结构的前提下, 改变振动的振幅, 气动力也会产生相应幅度的变化, 因此结构振幅对气动阻尼的影响可忽略不计。

4) 高频部分的广义力对气动阻尼的贡献很小, 即结构振动引起的广义力变化对气动阻尼起主要作用, 而脉动压力对计算气动阻尼影响不大, 可忽略不计。