混合网格重叠方法在多体相对运动中的应用

康忠良 方媛媛

(1.中国建筑科学研究院建研科技股份有限公司，北京 100013；2.北京市燃气集团研究院，北京 100011)

混合网格重叠方法在多体相对运动中的应用

康忠良 方媛媛

(1.中国建筑科学研究院建研科技股份有限公司，北京 100013；2.北京市燃气集团研究院，北京 100011)

本文发展了混合网格的动态重叠方法。基于线性重构方法，给出了一种适用于混合网格间的插值策略。所建立的网格间边界定义方法的性能优秀可靠，网格装配后的插值边界光滑且网格匹配性较好。耦合求解非定常NS方程和刚体6DOF运动方程来模拟多体相对运动问题。外挂物分离问题的数值计算表明，所发展的动网格方法对于处理多体相对运动模拟问题是快速准确可靠的。

动网格；重叠网格；混合网格；多体运动

1 引言

在计算流体力学(CFD)应用中，采用混合网格方法不但对复杂构型具有强大的几何适应能力，网格生成的人工工作量少，而且易于通过采用不同的单元类型及调整网格的疏密来适应不同的流场特征，并容易生成整体网格、整体求解。但是，在实际应用中，混合网格方法也遇到了一些困难，比如对于多体相对运动问题，就需要对网格进行再生，处理非常复杂。

为了解决这一问题，Nakahashi[1]首先提出了非结构网格重叠方法。该方法吸收了结构网格重叠方法[2]的优点，后来又被推广到混合网格重叠及动态重叠应用中[3-5]，得到了较大的发展。混合网格动态重叠方法的优点就是在模拟多体间具有大幅相对运动问题时，在物体运动过程中不需要网格再生，算法实现简单，并且可以获得较高的动网格处理效率。

尽管如此，混合网格重叠方法目前仍然面临诸多难题，其中，网格间边界的定义及其插值方法就是其中一个重要方面，已有研究均基于多层网格节点的流场变量进行插值计算，该思路实现较为繁琐，且容易在高速流动强间断附近产生较大的耗散。另一方面，在实际工程应用中还应该强调动网格处理的自动化程度、最大限度地减少人工干预，否则也不能体现混合网格重叠方法的优势。

本文在优化网格间边界定义的基础上，给出一种适用于任意单元类型的插值策略，并给出一种完全自动化的高效的网格动态装配策略。

2 数值方法

本文算法构造基于任意多面体混合网格单元。算法不考虑具体的网格拓扑，对不同的单元类型统一处理。空间离散采用格心有限体积法，控制体取网格单元，控制面取单元表面，流场变量存储在单元中心。

2.1 控制方程

舍去源项的三维非定常可压缩NS方程组在直角坐标系下的守恒积分形式可表示为

(1)

式中：Ω为控制体；∂Ω为控制面；Q为守恒变矢量；Fc为对流通量；Fv为粘性通量。各项的具体描述详见文献[6]。

2.2 数值离散方法

对流通量离散采用Roe[7]格式近似求解Riemann问题。假设控制面IJ的左右单元分别为I和J,则对流通量可表示为:

(2)

Barth[8]提出的线性重构方法假定解在控制体内呈线性分布，面左右两侧的值可表示为

(3)

式中：U为任意流场变量值；▽UI为单元I体心的梯度；Ψ为限制器函数；r为从体心到面心的矢量。▽UI采用最小二乘法[7]计算，Ψ采用收敛特性非常好的Venkatakrishnan[9]限制器。

另外，粘性通量采用中心格式离散，湍流模型采用Spalart-Allmaras[10]一方程模型，时间离散采用格式。非定常流动控制方程采用双时间步长法[11]求解，每一个物理时间步上的伪时间推进采用LU-SGS隐式格式[6]。

3 混合网格动态重叠方法

与现存的基于格点[3-5]的混合网格重叠方法不同，本文基于格心展开，通过网格间边界定义后，初始网格单元被分为三类：活动单元、插值单元和非活动单元。活动单元在计算区域内部参与流场计算，插值单元分布在网格边界用于子网格间的信息交换，非活动单元作为非计算单元被挖去。另外，本文定义包围插值单元网格中心点的单元为其宿主单元。图1给出了一个网格边界定义示例。其中，子网格G2中单元B(1-2-3)是子网格G1中单元A(a-b-c)的宿主单元，这是因为单元A的格心落在了单元B的内部。

3.1 插值方法

在本文混合网格重叠系统的流场计算中，子网格间的信息交换需要插值单元通过其宿主单元来获取其它子网格的相应边界信息，这一信息交换过程称为插值。与传统的完全依赖于格点流场变量进行插值的方法不同，下面给出一种新型的插值策略。

该方法依赖于插值单元的格心流场变量值、单元梯度值和限制器函数值的计算来完成。比如图1，插值单元A的状态可以通过下式确定

▽UA=▽UB

ΨA=ΨB

(4)

通过上式计算，插值单元A可以充分利用宿主单元B的所有相邻单元的信息。这是因为式(4)与式(3)相类似，▽UB和ΨB都是通过单元B的所有相邻单元构造得来。显然，这里要求宿主单元B在流场重构时能够构造出足够多的模板。

当网格间边界重叠区附近流动梯度变化较小时，该插值方法与线性重构一样，可使空间离散精度达到二阶。当边界附近流动变化剧烈时，也可以通过调整重叠区网格的分辨率来提高插值精度，并且这时满足插值单元与其宿主单元的尺寸匹配尤为重要。

因此，上述插值方法应用简洁方便并且能够保证计算精度。它只需要一层插值单元，就可以充分利用其宿主及所有相邻单元的信息。并且，因为借助了上一时间步已经确定的梯度值，所以同梯度计算方法一样，该插值方法适用于任意网格单元类型间的插值计算。

图1 两子网格间边界定义示例

另外，应用本文插值方法，插值边界可以同其它边界条件一样，在复杂求解器中可以做到透明处理。即迭代计算中，插值单元与其它边界虚网格单元功能类似。

3.2 网格间边界定义方法

本文中网格间边界定义过程可归纳为两步：网格分类和边界优化。首先通过以物面距离作为网格分类参数，将网格单元分为活动单元和非活动单元，形成初始分类边界。然后在初始分类的基础上，对边界进行优化，将网格单元分为活动单元、插值单元和非活动单元，最终形成插值边界。边界优化过程主要考虑如下几点因素：

(1)保证插值单元必须存在宿主单元。因为插值单元的流场信息是从它在其它网格的宿主单元获得的，所以只有保证存在宿主单元，流场计算才能正确进行。

(2)保证插值单元只存在一个宿主单元。因为实际应用中可能是任意多个子网格的任意单元类型重叠在一起，所以对一给定插值单元可能存在多个宿主单元。

(3)清除孤点。所谓孤点就是与周围单元的类别不相同的局部小部分网格单元，它可能是活动单元，也可能是非活动单元。在流场计算中，孤点的存在不但会损害计算精度而且还会降低收敛效率。

(4)适当把插值边界外移，即保证适当的重叠区域，以保证流场计算中宿主单元重构时足够的模板需求，从而有利于提高流场重构精度。

(5)光滑插值边界，且尽量保持插值单元和宿主单元的网格尺寸一致，从而有利于提高插值精度。

图2给出了六圆柱和五球体不同单元类型的边界定义结果。可见，本文方法定义的网格间边界光滑，网格尺寸大小匹配，网格重叠区域大小适当，说明本文给出的网格间边界定义方法是优秀可靠的。

3.3 混合网格动态重叠方法

本文混合网格动态重叠方法具体算法如下：

(1)获取所有子区域初始网格；

(2)采用2.2节的混合网格重叠方法装配网格；

(3)对所有活动区域，在当前物理时间步上迭代求解非定常NS方程，并积分出运动物体受到的气动力和气动力矩；

(4)基于运动物体受到的气动力及气动力矩、其它外力及力矩、重力和其它相关已知量等，在当前物理时间步求解刚体6DOF运动方程，得到下一时刻物体的位移和姿态角；

(5)根据步骤(4)的计算结果，移动运动物体所在子网格；

(6)判断非定常计算是否达到要求，如果达到，则结束计算；否则返回步骤(2)。

在上述计算流程中，步骤(2)到(6)不需要任何人工干预，完全是自动完成的。并且实践中发现，网格装配消耗时间与网格量以及子网格分布结构有关，但是一般情况下，网格装配一次所用时间与伪时间步上迭代一次消耗的时间基本是同一量级。相对于整个非定常计算过程而言，网格装配过程所消耗时间是一个小量。

图2 重叠网格边界定义示例

图3 三维机翼与外挂物分离网格

图4 三维外挂物分离不同时刻重叠网格及压力分布云图(x=5.3截面)

图5 三维外挂物分离不同时刻壁面压力分布云图

因此，本文采用的基于混合网格动态重叠方法的动网格技术具有如下优点：算法简单可靠，网格装配质量好效率高，处理复杂外形及复杂运动方式能力强，且可以完全实现自动化。

4 算例验证

本节通过模拟一个三维机翼与外挂物的分离过程对本文数值方法进行验证。该算例是一个考核动网格方法的标准算例，具有详实可靠的风洞实验数据，其模型由机翼、挂架和外挂物组成，详细几何参数及实验数据见参考文献[12]。

如图3所示，采用的重叠网格包括三个子网格，分别是背景网格(共65万)、绕机翼的子网格(共193万)和绕外挂物的子网格(共287万)。背景网格采用均匀的六面体单元覆盖整个分离区，其它外围区域采用四面体填充，绕机翼和绕外挂物的子网格在壁面附近均采用各向异性网格单元。

本算例模拟物理时间步长取为0.001s。图4给出了不同时刻重叠后的网格及压力分布云图，图5给出了不同时刻壁面压力分布云图。可见，随着时间的推进，三个子网格在不同时刻自动装配，并且由于均匀背景子网格的存在，重叠区附近保证了网格单元尺寸大小匹配。计算的压力在重叠区附近过度光滑，分布合理。

图6给出沿外挂物轴向压力取值的几个不同位置，图7给出了初始时刻不同位置的壁面压力分布对比，可见数值模拟结果与实验结果比较一致。

图6 外挂物压力取值的几个不同位置

图7 外挂物初始时刻壁面压力分布

图8 外挂物运动参数随时间变化

图8给出了外挂物运动参数随时间的变化曲线。从图中可以清楚地看到弹射力的消失对外挂物质心速度和角速度的影响位置。计算的外挂物质心位移、速度、角速度、滚转角、偏航角及俯仰角均与实验结果吻合较好。

5 结论

(1)本文插值方法依赖于插值单元的格心流场变量值、单元梯度值和限制器函数值的计算来完成，只需要一层插值单元，就可以充分利用其宿主及所有相邻单元的信息。

(2)本文给出的网格间边界定义方法是优秀可靠的，所定义的网格间边界光滑，网格尺寸大小匹配，网格重叠区域大小适当。

(3)本文采用的基于混合网格动态重叠方法的动网格技术具有如下优点：算法简单可靠，网格装配质量好效率高，处理复杂外形及复杂运动方式能力强，且可以完全实现自动化。

(4)三维外挂物分离问题的数值计算表明，本文所发展的动网格方法对于处理多体相对运动模拟问题是快速准确可靠的。

[1]Kazuhiro Nakahash, Fumiya Togash, Dmitri Sharov. An intergrid-boundary definition method for overset unstructured grid approach. AIAA-1999-3304, 1999.

[2]Benek J A, Steger J L, Dougherty F.A flexible grid embedding technique with applications to the Euler equations.AIAA-1983-1944, 1983.

[3]Hong Luo, Dmitri Sharov,Joseph D Baum.Anoverlappingunstructured gridmethod for viscous flows.AIAA-2001-2603, 2001.

[4]Fumiya Togashi, Yasushi Ito,Kazuhiro Nakahashi. Extensions of overset unstructured grids to multiple bodies in contact.Journalof Aircraft 2006; 43(1): 52-57.

[5]Shun Takahashi, IchieMonjugawa, and Kazuhiro Nakahashi. Unsteady flow computation around moving multiple bodies using overset unstructured grids. AIAA-2006-2839, 2006.

[6]J Blazek.Computational fluid dynamics principles and application. Oxford: Elsevier Sci Ltd, 2001.

[7]Roe P L. Approximate Riemann solvers, parameter vectors, and difference schemes.Journal of Computational Physics 1981; 43: 357-372.

[8]Barth T J, Jespersen D C. The design and application of upwind schemes on unstructured meshes.AIAA-1989-0366, 1989.

[9]Venkatakrishnan V. On the accuracy of limiters and convergence to steady state solutions.AIAA-1993-0880, 1993.

[10]P R Spalart, S R Allmaras. A one-equation turbulence model for aerodynamic flows.AlAA-1992-0439, 1992.

[11]Jameson A..Time-dependent calculations using multigrid with application to unsteady flows past airfoils and wings[R]. AIAA-1991-1596.

[12]HallL. H.,ParthasarathV.Y..Validation of an automated chimera/6-dof methodology for multiple moving body problems[R]. AIAA-1998-0767.

A Dynamic Overset Mixed Grids Approach for Multiple Relative Motion

Kang ZhongLiang1, Fang Yuanyuan2

(1.CABRTECHNOLOGYCO.,LTD.，ChinaAcademyofBuildingResearch,Beijing100013,China；2.AcademyofBeijingGasGroupCo.,Ltd.,Beijing100011,China)

A dynamic overset mixed grids algorithm is proposed. An interpolation method based on the linear reconstruction is developed for the overset mixedgrids approach. The optimized intergrid-boundary definition method is shown to be excellent and reliable. The intergrid boundaries are very smooth and have similar cell sizes.The flows around multiple bodies in relative motion can be simulated by solving the unsteady Navier-Stokes equations and the 6DOF equations. The numerical results of store separation prove that the present dynamic overset grids method is fast, accurate and reliable.

Dynamic Grids; Overset Grids; Mixed Grids; Relative Motion

“十二五”国家科技支撑计划(2012BAJ09B04)

康忠良(1981-)，工学博士。主要从事CFD技术和BIM技术研究。

V 211.3

A

1674-7461(2015)01-0069-06