CN103500144A - 一种基于动态边界条件的助推器再入运动分析方法 - Google Patents

Abstract

本发明是一种预测助推器再入过程气动载荷和运动规律的分析方法。该方法从求解助推器再入过程的三维非定常流场出发，基于多重参考坐标系原理和滑移网格技术，结合动态边界条件，进行助推器再入运动分析。根据牛顿运动定律和动量矩定理编写用于计算助推器再入过程中高度、速度和俯仰角速度等参数的用户程序。该用户程序通过控制旋转域的转动来实现助推器的俯仰运动，通过改变静止域的进口、出口和远场边界条件参数来模拟助推器再入运动过程的物理环境。时间推进求解非定常流场过程中，求解器实时调用该用户程序，以实现边界条件和旋转域转速的实时更新，直至助推器到达地面时停止计算，从而获得助推器的瞬态气动载荷和运动参数。

Description

一种基于动态边界条件的助推器再入运动分析方法

技术领域

本发明是针对预测助推器再入过程气动载荷和运动规律的分析方法。助推器再入至着陆过程中，时间历程和空间跨度相比其运动特征时间和长度都较大，助推器的外流场和气动载荷随时间和姿态变化，非定常流场产生的瞬态载荷直接影响助推器的运动规律。因此，准确获得再入过程的非定常流场和瞬态载荷是进行助推器再入运动分析的首要条件，本发明从求解助推器的三维非定常流场出发，建立运动分析模型，提出助推器再入运动分析方法。

背景技术

可重复使用运载器是降低航天运输成本的有效手段，也是未来航天运载技术发展的必然趋势。目前，国内外大多数的航天运载器都是一次性使用的，现有的技术水平距离航天运载系统完全可重复使用的目标仍有相当大的差距，因此对可重复使用运载器的需求推动了该领域的技术研究。为实现可重复使用的目的，国内外对作为运载器的火箭系统开展了不同层次的运动分析和回收技术研究，其中对助推器的再入运动和回收已成为航天运载领域关注的热点。国外发达国家的航天机构相继开发了再入运动分析工具，以预测助推器再入过程中的载荷和运动规律变化，为回收提供必要参数。如美国NASA开发的DAS和ORSAT软件系统，欧洲ESA开发的SESAM和SCARAB软件系统，日本NASDA在ORSAT基础上开发的ORSAT-J软件系统，这些软件系统能够实现对运载器再入至着陆整个过程的载荷和运动进行预测。相比之下，公开发表的文献表明国内在此领域的系统研究较少。为了获得助推器再入过程中瞬态气动载荷作用下的运动规律和姿态，本发明基于多重参考坐标系原理和滑移网格技术，结合动态边界条件，提出一种助推器再入运动分析方法。

多重参考坐标系原理最早应用于叶轮机械转/静流道的定常/非定常流场计算，其主要思想为：通过交界面将计算域分解为转动域和静止域，对包含旋转壁面的转动域求解转动坐标系下的流动控制方程，对静止计算域求解静止坐标系下的流动控制方程。转动域和静止域之间通过交界面进行数据传递，交界面的处理要保证质量、动量和能量守恒，进而实现各计算域流场同时求解，得到全流场解。

转动/静止坐标系下的流场速度之间的关系是:

$\stackrel{\→}{U_r}=\stackrel{\→}{U}-(\stackrel{\→}{\mathrm{\ω}}\×\stackrel{\→}{r})---\left(1\right)$

其中：

为转动坐标系下的速度矢量；

为静止坐标系下的速度矢量；

为转动角速度矢量；

为相对位置矢量。

转动坐标系下的连续方程为：

$\frac{\∂\mathrm{\ρ}}{\∂t}+\▿\·\left(\mathrm{\ρ}\stackrel{\→}{U_r}\right)=0---\left(2\right)$

在静止计算域中，将（2）式中的速度矢量

替换为

转动坐标系下的动量方程为：

$\frac{\∂\left(\mathrm{\ρ}\stackrel{\→}{U_r}\right)}{\∂t}+\▿\·(\mathrm{\ρ}\stackrel{\→}{U_r}\×\stackrel{\→}{U_r})+\mathrm{\ρ}(2\stackrel{\→}{\mathrm{\ω}}\×\stackrel{\→}{U_r}+\stackrel{\→}{\mathrm{\ω}}\×\stackrel{\→}{\mathrm{\ω}}\×\stackrel{\→}{r}+\stackrel{\→}{a}\×\stackrel{\→}{r})=\▿\·(-\mathrm{p\δ}+\mathrm{\τ})+S_M---\left(3\right)$

其中应力张量为：

$\mathrm{\τ}=\mathrm{\μ}(\▿\stackrel{\→}{U_r}+{(\▿\stackrel{\→}{U_r})}^T-\frac{2}{3}\mathrm{\δ}\▿\·\stackrel{\→}{U_r})---\left(4\right)$

在静止计算域中，将（3）式中的速度

替换为

并去掉因旋转所引起的离心力项科氏力项

和欧拉力项

转动坐标系下的能量方程为：

$\frac{\∂(\mathrm{\ρ}{H}_r-p)}{\∂t}+\▿\·\left(\mathrm{\ρ}\stackrel{\→}{U_r}{H}_r\right)=\▿\·(k\▿T+\mathrm{\τ}\·\stackrel{\→}{U_r})+S_E---\left(5\right)$

在静止计算域中，将（5）式中的速度

替换为转焓

替换为静止坐标系下的总焓 $h^{*}=h+\frac{1}{2}{\stackrel{\→}{U}}^2.$

滑移网格技术多用于解决旋转域相对静止域转动的非定常流场计算问题，如叶轮机械动/静流道非定常流场数值模拟。在非定常流场计算时，转动域网格随着离散时间步的变化相对静止域沿旋转轴转动，转动计算域网格不变形。每个时间步内，转动域和静止域交界面两侧网格节点的流场参数通过滑移网格技术进行传递和处理。对于每个物理时间步，都要确定新的网格空间位置和交界面属性，网格空间位置、交界面属性更新之后，在多重参考坐标系下进行流场计算，来获得各物理时间步的流场解。

发明内容

考虑到助推器的运动特征，本发明针对助推器的三维流场建立旋转和静止计算域，将内流叶轮机械转/静流道非定常流场的计算方法应用到助推器的外流场和姿态求解中，即多重参考坐标系系原理和滑移网格技术，结合动态边界条件，通过时间推进求解非定常流场来获得助推器的瞬态气动载荷和再入运动参数。

本发明根据牛顿运动定律和动量矩定理编写用于计算助推器再入过程中海拔高度、下落速度和俯仰角速度等参数的用户程序，采用的主要公式如（6）和（7）所示。公式（6）通过控制旋转域的转动来实现助推器的俯仰运动；公式（7）通过改变静止域的进口、出口和远场边界条件参数来模拟助推器下落运动过程的物理环境。非定常流场计算过程中，CFD求解器实时调用该用户程序，以实现边界条件和旋转域转速的实时更新。

$\left\{\begin{array}{c}{M}_x^i=J_x\·{\mathrm{\α}}_x^i\\ \mathrm{\Δ}{\mathrm{\θ}}^i={\mathrm{\ω}}_x^i\·\mathrm{\Δt}\\ \mathrm{\θ}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\mathrm{\Δ}{\mathrm{\θ}}^i\\ {\mathrm{\ω}}_x^{i+1}={\mathrm{\ω}}_x^i+{\mathrm{\α}}_x^i\·\mathrm{\Δt}\end{array}\right.---\left(6\right)$

$\left\{\begin{array}{c}\mathrm{mg}-f_z^i=m\·a_z^i\\ \mathrm{\Δ}{H}^i=V_z^i\·\mathrm{\Δt}+0.5\·a_z^i\·{\left(\mathrm{\Δt}\right)}^2\\ H=H_0-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\mathrm{\Δ}{H}^i\\ V_z^{i+1}=V_z^i+a_z^i\·\mathrm{\Δt}\\ -f_y^i=m\·a_y^i\\ V_y^{i+1}=V_y^i+a_y^i\·\mathrm{\Δt}\end{array}\right.---\left(7\right)$

其中：M_x为俯仰力矩，J_x为俯仰转动惯量，ω_x和α_x分别代表俯仰角速度与俯仰角加速度，θ为累积俯仰转角，m为助推器质量，f_y和f_z为水平与竖直方向气动力，V_y和V_z为水平运动与竖直下落速度，a_y和a_z为水平运动与竖直下落的加速度，H为再入高度，初始高度H₀≤20km，i和Δt代表物理时间步与时间步长，时间历程t_total≤5min。

附图说明

图1.计算域和滑移交界面2D示意图；

图中Stationary domain和Rotating domain代表静止和旋转计算域；Sliding interface代表滑移转/静交界面；ω代表旋转域转速，即助推器的俯仰角速度；Inlet、Outlet和Far-field代表计算域的边界。

图2.再入运动分析模型3D示意图；

图中圆柱形旋转计算域的中心即为助推器的质心。

图3.某助推器下落过程高度时程曲线；

图4.某助推器下落过程俯仰力矩时程曲线；

图5.某助推器下落过程俯仰角速度时程曲线；

图6.某助推器下落过程累积俯仰转角时程曲线；

图7.某助推器下落过程气动力时程曲线；

图8.某助推器下落过程速度时程曲线。

具体实施方式

为更清楚地描述本发明，具体给出助推器再入运动分析的主要步骤：

（1）采用建模软件建立助推器的几何模型，并确定质心的空间位置，根据助推器的质量分布计算俯仰转动惯量等相关参数；

（2）确定静止和转动计算域的尺寸，采用分网软件划分助推器外流场计算网格；

（3）确定再入初始条件，如初始高度、初始速度、初始俯仰角速度等，编写用于确定助推器再入过程中海拔高度、速度和俯仰角速度变化的用户程序；

（4）在CFD求解器中设置求解参数，实时调用（3）中编写的用户程序，以实现非定常流场求解过程中边界条件和旋转域转速的实时更新，时间推进迭代求解非定常流场，直至助推器到达地面时停止计算；

（5）对计算结果进行分析处理，获得助推器再入过程中瞬态气动载荷作用下的运动规律。

Claims (2)

1.一种助推器再入运动分析方法，其特征是，从时间推进求解助推器再入过程的非定常流场出发，基于多重参考坐标系原理和滑移网格技术，结合动态边界条件，建立运动分析模型，模拟和预测助推器再入过程的气动载荷和运动规律。

2.用于计算助推器再入过程中海拔高度、下落速度和俯仰角速度等参数的用户程序，该用户程序根据牛顿运动定律和动量矩定理编写，采用的主要公式如（1）和（2）所示。公式（1）通过控制旋转域的转动来实现助推器的俯仰运动；公式（2）通过改变静止域的进口、出口和远场边界条件参数来模拟助推器下落运动过程的物理环境。非定常流场计算过程中，CFD求解器实时调用该用户程序，以实现边界条件和旋转域转速的实时更新。

$\left\{\begin{array}{c}{M}_x^i=J_x\·{\mathrm{\α}}_x^i\\ \mathrm{\Δ}{\mathrm{\θ}}^i={\mathrm{\ω}}_x^i\·\mathrm{\Δt}\\ \mathrm{\θ}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\mathrm{\Δ}{\mathrm{\θ}}^i\\ {\mathrm{\ω}}_x^{i+1}={\mathrm{\ω}}_x^i+{\mathrm{\α}}_x^i\·\mathrm{\Δt}\end{array}\right.---\left(1\right)$

$\left\{\begin{array}{c}\mathrm{mg}-f_z^i=m\·a_z^i\\ \mathrm{\Δ}{H}^i=V_z^i\·\mathrm{\Δt}+0.5\·a_z^i\·{\left(\mathrm{\Δt}\right)}^2\\ H=H_0-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\mathrm{\Δ}{H}^i\\ V_z^{i+1}=V_z^i+a_z^i\·\mathrm{\Δt}\\ -f_y^i=m\·a_y^i\\ V_y^{i+1}=V_y^i+a_y^i\·\mathrm{\Δt}\end{array}\right.---\left(2\right)$

其中：M_x为俯仰力矩，J_x为俯仰转动惯量，ω_x和α_x分别代表俯仰角速度与俯仰角加速度，θ为累积俯仰转角，m为助推器质量，f_y和f_z为水平与竖直方向气动力，V_y和V_z为水平运动与竖直下落速度，a_y和a_z为水平运动与竖直下落的加速度，H为再入高度，初始高度H₀≤20km，i和Δt代表物理时间步与时间步长，时间历程t_total≤5min。

Images