临近空间飞行器多物理场耦合建模的网格映射方法

刘 璟 王 玲，2 胡东飞 铁 鸣，2 吴旭生

1.北京临近空间飞行器系统工程研究所，北京100076

2.空间物理重点实验室，北京100076

临近空间通常是指20～100km的空域［1］。临近空间高超声速飞行器具有快速、机动、灵活的特点，具有极大的发展潜力［2］，对于维护国家安全和利益具有极为重要的战略意义。

临近空间飞行器的研制面临许多复杂的问题和挑战，其功能复杂、机理复杂、工艺复杂，具有多自由度、多变量、非线性、强耦合、参数时变等综合复杂性［3］，而且涉及总体、气动、结构、控制、防热、动力等多个学科和专业领域的相互作用、高度耦合的子系统。传统的各学科独立设计与仿真验证的模式难以体现各学科之间的耦合关系。需要开展多学科、全系统、全流程、多物理场耦合过程的协同仿真研究，验证各学科耦合关系下的更接近于真实情况的飞行性能。学科间存在着大量的耦合与交互关系，其中一些涉及学科间交互的复杂仿真问题需要多个学科领域的仿真模型、软件相互协作共同完成［4］。

本文研究并给出了临近空间高超声速飞行器的各学科之间的一种耦合建模方法。由于计算能力和研制习惯等方面的制约，气动、结构、防热等专业的仿真建模和模型解算工具不同。流体力学领域侧重于使用有限体积法，固体力学领域广泛采用有限元方法，且会按需求选用不同的计算工具，包括ABAQUS，ANSYS，NASTRAN和FORTRAN等，这就带来了不同的网格划分方式。在多物理场耦合仿真计算过程中，学科间复杂的交互关系和异构网格之间的信息传递就成为一个难点。为解决该问题，研究并实现了异构网格耦合界面之间以及耦合域之间的信息传递方法，并应用于临近空间高超声速飞行器的多物理场耦合仿真。

1 临近空间高超声速飞行器多物理场耦合建模

临近空间高超声速飞行器各学科之间的耦合关系复杂，其中包括流固耦合问题［5－6］、热固耦合问题和热流固耦合问题等。

针对计算流体力学(CFD)/计算固体力学(CSD)耦合已经发展了全耦合、紧耦合和松耦合［7］3种耦合方式。飞行器气动与结构之间的耦合关系主要是:在高超声速条件下，气动载荷引起的飞行器结构表面位移不可忽略，而表面位移反过来会影响流场和气动载荷。

飞行器热固领域之间和热流领域之间的耦合关系主要是:在高超声速条件下，气动加热引起的高温使得飞行器结构中产生热应力，进而导致结构应变，温度的变化也使得材料强度和刚度相应改变;飞行器的结构应变以及端头和翼前缘等处的热烧蚀改变了气动外形，进而对飞行器所受的气动力产生影响。热流固的相互耦合关系还会影响飞行器的控制系统精度和可控性。

由于临近空间高超声速飞行器的高度复杂性，目前难以进行全系统的多场紧耦合的仿真建模与验证，因此，采用松耦合方法给出了一种多物理场耦合的协同仿真流程。气动模型基于CFD网格计算飞行器外表面所受载荷，防热模型和结构强度模型分别采用ABAQUS和ANSYS作为网格生成工具和模型解算工具。如图1，在一个仿真步长内，气动、防热和结构强度仿真模型按其耦合关系依次解算，步骤序号说明了解算和数据交互的顺序，双箭头标明了数据传递方向。

图1 耦合仿真流程图

2 异构网格映射方法

在物理场耦合的各学科模型之间的信息交互中，存在以下几类问题:

1)异构网格耦合界面位移映射;

2)异构网格耦合界面载荷映射;

3)异构网格耦合域位移映射;

4)异构网格耦合域载荷映射。

其中，耦合界面上映射方法作为流固耦合等问题的一个重要方面，是一个研究热点。常用的方法有常体积转换法［8］、加权余量法［9］、平板样条法、反距离加权法、径向基函数法等很多种。

2.1 映射方法的原理和思想

异构网格之间信息传递方法应当满足连续条件和守恒条件，即网格信息所体现的外形等物理信息是连续的，网格的划分只是一种离散方法;映射完成后，两套网格的整体或局部的能量都必须是守恒的。连续条件和守恒条件可以作为映射算法的推导依据或验证依据。

可将各种异构网格之间信息传递方法的核心思想总结为以下几类:

1)拟合:由网格离散点的信息拟合出连续的形函数或场函数等，并在另一套网格中应用该拟合函数求得各点的值;

2)插值或加权:利用周围各点对目标点影响程度的不同，且相邻网格点的信息是连续变化的特点，求得目标点的信息;

3)守恒:从方法上保证整体和局部各处的信息在映射前后是严格守恒的，例如基于静力等效、能量守恒等原理推导出的映射公式。

2.2 本文研究和应用的映射方法

比较几类映射方法，函数拟合的方法可以在一定程度上还原离散网格点的连续函数，但对于复杂外形或复杂结构的情况，拟合函数的精度会严重下降;守恒类方法推导过程十分严谨，往往存在矩阵求逆过程，难以避免奇异矩阵，在实际工程应用中处理起来比较复杂;插值或加权类方法最为常见，因此采用局部信息的插值或加权方法处理异构网格耦合界面和耦合域的信息传递问题。

信息从一套网格传递到另一套网格，为便于区分，下文把这两套网格分别称作“源网格”和“目标网格”。

2.2.1 表面网格载荷和位移映射

常见的表面网格单元的形状只有三角形和四边形2种，将四边形也作为2个三角形处理。对于耦合界面上的目标网格的每个网格点Ps0都可以唯一确定1个对应的源网格三角单元。将该点Ps0投影在三角单元所在平面上，得到如图2所示的点Ps。

［8］的面积坐标定义，如图2，得到点Ps的面积坐标Ni为:

图2 面积坐标的定义

式中，A为ΔP1P2P3的面积，A1为ΔPsP2P3的面积，A2为ΔPsP1P3的面积，A3为ΔPsP2P1的面积。

处理气动载荷(单位为帕斯卡)映射时，可以将CFD网格三角单元面积坐标视为三角单元各顶点载荷fi对ANSYS网格点Ps的权重，得到Ps点的载荷为:

对于每个ANSYS网格点按上述面积权重法进行处理，即可完成异构网格耦合界面上的载荷映射过程。

处理位移映射时，以常体积转换法为基础，对上述搜索到的三角单元逐个进行计算。按式(3)和式(4)即可得到耦合界面上网格点Ps0进行位移映射后的坐标

式中，加撇的点表示该点位移后的坐标，如A'表示位移后的三角单元面积。

2.2.2 体网格载荷映射

常见的体网格单元形状包括四面体、六面体、三棱柱等多种，需要采用一种与网格形状无关的方法处理异构网格耦合域上的载荷映射问题。

在反距离加权法的基础上，以从ABAQUS网格向ANSYS网格进行热载荷映射为例，得出如下映射方法:

1)将两套网格重叠在一起，在耦合域上搜索每个ANSYS网格点的坐标位于哪个ABAQUS网格单元中。

已知P1，P2，P3，Pn4个点的坐标，其中，点P1，P2，P3构成一个平面。按行列式的计算法则，有

判断点是否在原网格单元内的原理是:

将某个ABAQUS网格单元中心点P0坐标和ANSYS网格点Ps分别代入式(5)，计算结果同号，则可判断点P0与点Ps在平面S123的同侧。

对于一个均质凸n面体(n＞3)，可知其重心(中心)一定在n面体内。因此只要按上述方法判断点Ps与某体网格单元的中心相对于该单元所有的面同侧，即可知该点在该单元内部。

2)体网格单元各顶点对点Ps的权重Ni为:

式中，li表示第i个顶点与Ps之间的距离。

对于在由n个顶点构成的网格单元内的点Ps，其热载荷映射公式为:

需要指出的是反距离加权法使用了除法，且除数有可能为0或接近0，因此在编程实现时需要对两点重合或近似重合的情况进行判断。

2.2.3 体网格位移映射

飞行器结构受力发生变形时，在变形区域，每个微小的网格单元都会随之发生位移、压缩或拉伸。由于网格足够微小，发生扭曲变形的量所造成的误差可以忽略。在此假设基础上，本文提出一种解决体网格位移映射的“中值位移法”，即由每个ABAQUS网格节点所在的由n个节点构成的ANSYS体网格单元各顶点的平均位移求得该ABAQUS网格节点的位移，如式(9):

式中，Ai表示ANSYS体网格单元各节点，Q表示ABAQUS网格节点，加撇表示位移后的节点。

式(9)为中值位移法的数学表述，该方法的前处理过程，即确定节点与单元的位置关系的方法与式(5)所述相同。该方法的优点是步骤简单、易于实现，计算效率高。

2.3 异构网格信息传递算法库

将异构网格之间信息传递的各种方法整合到一起，形成异构网格信息传递算法库，如图3所示。算法库由搜索算法库、映射算法库、前处理模块和后处理模块组成。前处理模块按用户选择的网格类型动态生成信息传递所需的所有数据结构。后处理模块将映射结果按相应的网格文件格式输出，并释放内存。搜索算法库内包含:节点与表面网格单元位置关系搜索算法、节点与体网格单元位置关系搜索算法。映射算法库可以处理各种异构网格之间的信息传递问题，有包括上文所介绍的常体积转换法、面积权重法、中值位移法、反距离加权法在内的多种映射算法，并且可以不断丰富。

2.4 算法库精度分析

为了验证算法库的算法精度，选择了半径为0.1m的半球体和长0.2m、半径0.05m的圆柱体模型，各建立2套网格，如图4和图5所示。

给半球模型建立温度场分布函数，如式(10)所示:

其耦合界面与耦合域上的载荷映射平均误差见表1。

给圆柱模型表面加载位移分布函数，如式(11)所示:

δx沿径向加载，表面以下各层网格等比例压缩。其耦合界面与耦合域上的载荷映射平均误差见表1。

可见以上算例的精度符合工程计算要求，可应用于多物理场耦合建模与仿真试验。

3 结论

研究了临近空间高飞行器多物理场耦合建模方法，为解决领域间因不同设计习惯、使用不同的工具软件等造成复杂数据交互问题，研究并实现了适用于复杂外形和结构的各种异构网格耦合界面及耦合域之间信息传递的方法，封装成算法库，并可应用于临近空间飞行器多物理场耦合仿真。

图3 异构网格信息传递算法库

表1 算法库验证结果

参考文献

［1］ Stephens H.Near-space［J］.Air Force Magazine，2005，88(7):31.

［2］ 李怡勇，李智，沈怀荣.临近空间飞行器发展与应用分析［J］.装备指挥技术学院学报，2008，19(2):61-65.(LI Yi-yong，LI Zhi，SHEN Huai-rong.Analysis on Development and Application of Near Space Vehicle［J］.Journal of the Academy of Equipment Command＆Technology，2008，19(2):61-65.)

［3］ 铁鸣，王玲，孙学功.临近空间飞行器一体化仿真分析平台体系研究［C］.临近空间飞行器技术论坛，2009.(TIE Ming，WANG Ling，SUN Xue-gong.Study on Simulation Analysis Platform System for Near Space Vehicle［C］.Nearspace Vehicle’s Technology Forums，2009.)

［4］ 铁鸣，胡东飞，王玲，吴旭生.基于HLA的临近空间飞行器虚拟飞行试验方法［J］.信息与电脑，2010，7:51-52.(TIE Ming，HU Dong-fei，WANG Ling，WU Xu-Sheng.Virtual Flight Testing Method of Near Space Vehicle Based on HLA［J］.China Computer＆Communication，2010，7:51-52.)

［5］ 钱若军，董石麟，袁行飞.流固耦合理论研究进展［J］.空间结构，2008，14(1):3-15.(QIAN Ruo-jun，DONG Shi-lin，YUAN Xing-fei.Advances in Research on Fluidstructure Interaction Theory［J］.Spatial Structures，2008，14(1):3-15.)

［6］ 安效民，徐敏，陈士橹.多场耦合求解非线性气动弹性的研究综述［J］.力学进展，2009，39(3):284-298.(AN Xiao-min，XU Min，CHEN Shi-lu.An Overview of CFD/ CSD Coupled Solution for Nonlinear Aeroelasticity［J］.Advances in Mechanics，2009，39(3):284-298.)

［7］ Hurka J，Ballmann J.Elastic Panels in Transonic Flow［J］.AIAA Paper 2001-2722，2001.

［8］ Goura G S，Badcock K J，et al.A Data Exchange Method for Fluid-structure Interaction Problems［J］.The Aeronautical Journal，2001，105:215-221.

［9］ De Boer C，Van Zuijlen，Van A H，Bijl H.Review of Coupling Methods for Non-matching Meshes［J］.Computer Methods in Applied Mechanics and Engineering，2007，196:1515-1525.