基于动态混合网格的多体分离数值模拟方法

段旭鹏，常兴华，2，张来平，2

(1.空气动力学国家重点实验室，四川 绵阳 621000;2.中国空气动力研究与发展中心计算空气动力学研究所，四川 绵阳 621000)

0 引言

现代飞行武器正向集成化、隐身化、精确化、强突防等方向发展。例如，第四代战机为了提高隐身性能，普遍采用内埋式或半内埋式武器系统;新型战略和战术导弹为了提高突防能力，逐步采用多弹头分导技术。在这些飞行武器实施攻击的过程中，必然会出现多体分离和多体干扰问题。实际上，多体分离和多体干扰是航空航天飞行器研制和飞行过程中普遍存在的问题，例如:战斗机飞行员在紧急情况下抛离座舱、弹射逃生的过程，多级火箭在飞行过程中的级间分离，子母弹的抛壳与子弹分离，大型导弹和火箭的头部整流罩分离等都属于这类问题。

经过CFD工作者的长期努力，多体分离问题的数值模拟方法已经有了长足的进步，从最初的准定常计算逐步发展为非定常模拟。文献调研结果显示，目前，在多体分离过程中采用的动网格方法以重叠动网格技术为主。重叠网格的概念首先是由Steger等［1］提出来的，其最早应用于静态结构网格，以解决复杂外形的定常流数值模拟问题;后来逐步发展应用于静态非结构网格［2－3］;其后一些CFD工作者将其推广应用于多体分离问题的数值模拟［4－6］。这种方法是将计算网格分解为背景网格和运动网格子块，为适应边界的运动，令围绕运动边界的网格子块随分离体一起做刚性运动，运动网格子块与背景网格之间的流动信息通过插值运算进行通讯。重叠网格方法虽然比较容易实现，但是在非定常计算过程中，每一步都需要进行插值运算，一方面影响计算效率，另一方面插值误差可能导致非定常计算过程中的误差累积，并且在物理量间断处的插值精度会降低。

为了克服上述问题，CFD工作者发展了动态变形网格技术［7－8］。在物体运动过程中使用网格变形方法，只改变节点位置，不改变拓扑结构，这样避免了插值，流场的精度也就不会随网格的变形而降低。但是对于大位移的多体分离问题，分离体与母体的初始间距一般很小，大位移之后的变形网格质量难以保证，必须引入局部网格重构的方法以改进网格质量。

为了自动生成大位移或大变形问题的动态网格，作者所在的课题组发展了一种新的动态混合网格生成技术［9］。该方法将基于“弹簧原理”的节点松弛法［8］、基于“Delaunay”背景网格的插值法［7］和局部网格重构法有机结合，在小位移和/或小变形的情况下，利用Delaunay背景网格插值法进行网格的变形，而背景网格本身的变形采用基于“弹簧原理”的节点松弛法。在变形后引入网格质量检测，如果网格质量无法达到计算要求，则在局部进行网格重构，以提高计算网格整体质量。

本文将这一方法应用于多体分离问题的数值模拟，其间发展了六自由度运动方程计算方法以计算分离体运动轨迹，同时与非定常流动方程求解有机集成，形成了一种基于动态混合网格的多体分离问题数值模拟方法。利用该方法数值模拟了典型的多体分离标模问题(机翼/外挂物分离)，得到了与风洞试验一致的计算结果，表明本文方法具有良好的计算精度，具有广阔的应用前景。

1 刚体六自由度运动方程及求解

在多体分离问题中必须利用非定常气动力耦合飞行力学方程(即六自由度运动方程)进行运动轨迹的计算。为了计算方便，这里建立两个坐标系:与机翼固连的惯性坐标系和与外挂物固连的分离体本体坐标系。由于多体分离过程的时间不是很长、外挂物飞行距离不太远、姿态角变化不太大，所以不妨忽略地球形状及转动影响。由经典力学中蔡斯尔定理可知，任何刚体运动，都可以将其分解为刚体质心的平动和绕质心的转动。根据质心动力学、运动学方程以及欧拉动力学、运动学方程，我们可以得到刚体六自由度运动控制方程。

在惯性系下，根据牛顿定律给出刚体平动方程:

其中，m是外挂物质量;V、X分别是外挂物质心在惯性坐标系下的速度矢量和位置矢量;F是外挂物所受合力，包括气动力、重力以及弹射力。

在本体系下，给出刚体转动方程。外挂物的动量矩定义为:

其中I是惯性张量:

根据动量矩定理:

把(3)代入(5)得:

按照定义［10］，刚体的欧拉角(姿态角)是通过绕坐标轴做三次旋转得到的，如图1所示。转轴顺序为y、z、x，依次得到侧滑角β、俯仰角α、滚转角γ。其中各角度的定义为:俯仰角α是弹体轴向(ox)与水平面ox0oy0的夹角;偏航角β是弹体轴向在水平面内的投影与ox0的夹角;滚转角γ是弹体(oy)轴向与包含导弹x轴的铅垂平面的夹角。按右手系定义其转动正方向。从图1可以得到惯性坐标与本体坐标的转换矩阵，利用该矩阵及其逆矩阵可以实现各物理量在两个坐标系下的自由转换:

其中:

图1 欧拉角与三次旋转所得坐标系的关系Fig.1 Relationship of Euler angles and the coordinate systems

另外，由图1还可以得到欧拉角的时间导数与角速度的关系:

将式(6)和式(8)联立就得到含6个未知数的6个一阶常微分方程。至此，刚体转动问题已转化为一阶常微分方程组的初值问题，这里采用四阶Runge-Kutta方法对其求解［7］。

2 数值计算结果与分析

对于多体分离过程中的动态网格生成问题，本文综合利用基于Delaunay背景网格的插值法和局部网格重构法生成多体分离过程中的动态混合网格(具体方法参见文献［9］)。为了提高该方法的适应性和效率，我们对Delaunay背景网格的插值法进行了适当的改进，引入了基于“弹簧原理”的节点松弛法模块，用节点松弛法实现背景网格控制点的运动。通过在适当的位置加入可动的控制点对运动网格进行大范围的调控。相对于以往少量固定控制点的方法，网格质量得到进一步提高。图2显示的是分离过程中3个典型时刻的动态网格图。由于挂架和外挂物间的初始缝隙很小，缝隙间的初始网格量就比较少，导致分离过程不可避免地要对网格进行重构。对于本文计算的外挂物分离问题，在保证网格质量的前提下，整个分离过程的重构次数仅为4次。

图2 分离过程中的动态非结构网格Fig.2 Unstructured dynamic grids during separation

本文采用的流场计算软件是课题组发展的基于混合网格的定常/非定常CFD解算器。该解算器以二阶精度有限体积法为基础;在非定常计算中，采用的是将双时间步方法和基于块LU-SGS(BLU-SGS)的隐式计算方法相结合的非定常计算方法，同时考虑了动网格计算的几何守恒率，具有较好的非定常计算精度和效率，该解算器已经过比较严格的验证与确认，详见参考文献［11］。

在多体分离问题的计算中，我们将前述的动网格生成方法、六自由度运动方程求解和非定常流场计算耦合在一起，用来求解刚体在变化的力和力矩作用下的动力学和运动学问题。具体过程如下:

(1)利用非定常计算方法求解外挂物在某一时刻所受的气动载荷。

(2)以该气动载荷、外加弹射力和自身重力以及外挂物当前运动参数，包括位置、姿态、速度等为输入条件，调用六自由度运动求解器得到其下一时刻的运动参数。

(3)利用基于“弹簧原理”的节点松弛法和基于Delaunay背景网格的插值法的耦合动网格方法对外挂物附近的网格进行变形，得到新时刻的网格。

(4)判断网格质量。当外挂物运动位移很大导致局部网格质量急剧下降而低于某一预定值时，在局部重新生成网格，并通过插值得到新网格上的物理量;否则略过此步。

(5)重复(1)－(4)直至分离结束。

为了验证上述耦合计算方法，我们对前面提到的机翼/外挂物分离标模问题进行了数值模拟。首先利用定常流场解算器，得到分离前的跨声速静态计算结果。然后让外挂物在弹射力、重力和气动力的作用下运动。真实分离历程为0.32s。

计算条件与文献［13］一致(如表1所示)。弹射力是指为使分离干脆、迅速而施加在外挂物前后两端的作用力。在下落过程中，由于受机翼、挂架的干扰，外挂物开始会受到一个强烈的低头力矩，如不加控制就会出现尾部与挂架碰撞的事故。为了克服这个低头力矩，故意使作用于外挂物后端的弹射力比前端弹射力大。弹射力的作用方式也与实验相同:在外挂物本体系下，作用距离为0.33英尺;超过这个距离，弹射力立刻消失。

2.1 静态计算结果

初始静态计算网格的单元数总数约为120万，在流动复杂的局部区域进行了加密。图3给出了初始时刻外挂物两个典型纵向截面上的压力分布计算结果，同时给出了风洞试验结果［12］，可以看到二者吻合较好，说明静态计算结果是可靠的。

表1 计算条件Table 1 Computational condition

图3 Φ=90°，275°处的压力系数分布Fig.3 Pressure coefficient;Φ =90°，275°

2.2 非定常计算结果

图4给出了分离过程中典型时刻的压力等值云图。纵观分离全过程可以发现:分离之初，外挂物与机翼的相互影响很明显，外挂物尾部存在明显的低压区，同时在弹射力的作用下外挂物起初有抬头趋势;随着时间的推移，外挂物随之下落，二者的影响逐渐减小，低压区慢慢消失;随着弹射力的撤销，外挂物恢复低头趋势，其姿态渐趋水平。

图5给出了分离过程中的运动姿态、运动速度和气动力的变化曲线，其中的点标为风洞试验结果，Cal.是本文的计算结果，Ref.是文献［12］的计算结果。从图中可以看出，本文的计算结果与风洞试验结果吻合较好，尤其是在分离初期。在分离后期，由于计算中仅采用无粘模型，与风洞试验结果偏差较大。图5(a)显示的是质心三个方向的运动轨迹，可以看到外挂物的纵向位移显著，其原因是弹射力与重力之和远大于下落过程所受升力，因此向下的合力起主导作用。从图5(b-c)中可以看到，俯仰速度、俯仰角速度在0.05s前后发生突变，这是因为此刻外挂物的运动距离刚刚超过0.33英尺，弹射力突然消失导致了质心速度和角速度的突变。一般来说，力矩系数是最难计算的气动参数，由于本文采用无粘计算，忽略了粘性影响，从图5(f)中可以看到，初始时刻的俯仰力矩比实验值略低，进而导致最终结果与真实情况有一定的差异，这与文献［12］中的无粘计算结果是一致的。

图4 分离过程压力云图Fig.4 Surface pressure contours

图5 分离过程中物理量随时间的变化曲线Fig.5 Kinetic and aerodynamic parameters during separation

3 结论

本文在课题组以往发展的动态混合网格技术和非定常计算方法的基础上，发展了刚体六自由度运动计算模块，实现了运动轨迹、动网格生成和非定常流场的耦合计算，用来预测复杂多体分离问题。机翼/外挂物分离问题的数值计算结果表明本方法具有较好的精度，并对复杂外形具有良好的适应性，具有广阔的应用前景。下一步我们将开展非定常粘性流动的耦合数值模拟，进一步提高本文方法的计算精度和效率。

［1］STEGER J L，DOUGHERTY F C and BENEK J A.A chimera grid scheme［A］.ASME Mini-Symposium on Advances in Grid Generation［C］.1982.

［2］NAKAHASHI K，TOGASHI F.An intergrid-boundary definition method for overset unstructured grid approach［R］.AIAA －99－3304.

［3］LOHNER R，SHAROV D，LUO H，RAMAMURTI R.Overlapping unstructured grids［R］.AIAA－01－0439.

［4］WANG Z J and PARTHASARATHY V.A fully automated chimera methodology for multiple moving body problems［J］.Int.J.Numer.Meth.Fluids，2000，33:919 －938.

［5］李亭鹤.重叠网格自动生成方法研究［D］.［博士学位论文］.北京航空航天大学，2004.

［6］PPEWITT N C，BELK D M，SHYY W.Parallel computing of overset grids for aerodynamic problems with moving objects［J］.Progress in Aerospace Sciences，2000，36:117 －172.

［7］LIU X Q，QIN N，XIA H.Fast dynamic grid deformation based on Delaunay graph mapping［J］.J.Comput.Phys.，2006，211:405－42.

［8］BLOM F J.Consideration on the spring analogy［J］.Int.J.Numer.Meth.Fluids，2000，32:647 －668.

［9］张来平，段旭鹏，常兴华，张涵信.基于Delaunay背景网格插值方法和局部网格重构的动态混合网格生成技术［J］.空气动力学学报，2009，27(1):26－32.

［10］钱杏芳.导弹飞行力学［M］.北京:北京理工大学出版社，2006，7.

［11］ZHANG L P，WANG Z J.A block LU-SGS implicit dual time-stepping algorithm for hybrid dynamic meshes［J］.Computers＆ Fluids.2004，33:891－916.

［12］HEIM E R.CFD wing/pylon/finned store mutual interference wind tunnel experiment［R］.AFATIJFXA，Eglin Air Force Base，FL 32542 －5000，February，1991.

［13］HALL L H，PARTHASARATHY V.Validation of an automated chimera/6-DOF methodology for multiple moving body problems［A］.The 6th AIAA Aerospace Sciences Meeting and Exhibit［C］.January 12 －15，1998.