洞映射方法的研究和改进

袁 武 阎 超 席 柯

(北京航空航天大学 国家计算流体力学实验室，北京100191)

挖洞是重叠网格［1］的关键技术之一，若某重叠网格单元落入另一网格域的非可透面(如物面或人工指定的挖洞曲面)内，则应被标记“洞内点”，而不参与流场的计算.这一过程被形象地称之为“挖洞”.挖洞的过程，实际是在重叠网格中建立人工插值内边界(即洞边界)的过程，其数学实质等价于解决一个所谓“点与封闭曲面的相对位置关系”问题.

关于挖洞方法的研究主要是如何提高挖洞过程的效率、可靠性和自动化程度，常见的方法有点矢法［2］、射线求交法［3］、洞映射 (hole-map)［4］、Object X-Ray［5］、叉树结构挖洞方法［6］等.其中，洞映射方法通过构建辅助的直角笛卡尔网格，用笛卡尔网格近似物体的挖洞曲面，将点与曲面之间的关系转化为点与笛卡尔网格单元之间的简单关系，因此效率和自动化程度都很高，从而得到广泛的应用，如文献［7］.本文对传统的洞映射方法进行了研究，针对其在存储效率和可靠性方面的不足进行改进，并通过算例对新方法进行了验证.

1 洞映射方法简介

文献［4］认为，对于给定的重叠网格体系，若已知其拓扑结构和流动边界条件，就能够用均匀的笛卡尔网格单元去近似每个网格的挖洞曲面，从而得到该曲面的笛卡尔近似，称之为“洞映射(hole-map)”.

洞映射中的笛卡尔网格是一种与挖洞曲面固联的辅助性网格，根据相对位置的不同，笛卡尔网格单元的属性分为“洞内单元”、“洞外单元”和“边缘单元”.通过构建笛卡尔网格建立了整个挖洞曲面的映射集合后，就可以经简单的下标计算，方便地获得网格点在洞映射单元中的位置及属性.

对于封闭曲面，洞映射方法一般分为4个步骤:离散挖洞曲面、划分洞映射网格、标识边缘单元和区分洞内、外单元.

2 洞映射方法的改进

2.1 扩展的广义封闭和最小洞映射

虽然洞映射的过程涉及到笛卡尔网格，但无须存储这些网格单元的空间坐标，只需要组织一个整型数组保存单元的属性值.尽管如此，由于传统的洞映射方法要求挖洞面必须完全封闭，当物面范围较大或需精细描述时，笛卡尔网格数量很大，此时组织一个三维的整型数组所需存储开销也很大，有时甚至无法实现，这也是洞映射方法为人所诟病的主要原因.

图1是一般外形的导弹弹身+尾舵模型示意，弹身和尾舵分别生成结构网格，不计算底部阻力故弹身只生成前体网格.弹身网格中含对称面和远场边界，由文献［8］提出的“广义封闭”概念，通过适当的处理，是允许挖洞曲面结束于网格边界的.图2是构建的弹身洞映射笛卡尔网格，须对所有物面面元分析，用于计算笛卡尔网格的边界和单元尺寸，其中网格边界是整个弹身的笛卡尔包围盒(min，max).

在对尾舵网格的挖洞中，弹身的洞映射笛卡尔网格中，产生操作的映射单元，全部落在尾舵网格所占据的空间中，其余单元没有参与到挖洞过程.因此，一种十分自然的思想就是在建立弹身洞映射模型时，只在尾舵网格占据的局部区域建立洞映射，弹身物面落在尾舵网格外的应人为地截断，不予分析，这样，新的洞映射边界就由物面、网格边界和截断面构成，显然这是一个可供尾舵网格挖洞的最“紧凑”的区域，本文称为“最小洞映射”.

图1 弹身和尾舵网格示意

图2 弹身和尾舵全模洞映射模型

本文引入了网格截断面的概念，认为挖洞曲面结束于网格的截断面，同样属于封闭情况，这样就将“封闭”的概念扩展到了人为截断的边界上.对于扩展后的广义封闭，仍然可以应用广义封闭用于判断洞映射单元属性的一般方法，在2.2节将详细分析.图3是按最小洞映射重新构建的笛卡尔网格，图中“截断面”由尾舵网格在原弹身洞映射模型中“切除”得到，笛卡尔网格边界由物面、对称面、远场和截断面封闭，组成对尾舵网格挖洞的最小区域.

图3 弹身和尾舵最小洞映射模型

最小洞映射有效缩小了洞映射区域，大量减少了笛卡尔网格单元，节省了存储开销，也为复杂情况下建立更为精细的洞映射提供了必要条件.对图3的洞映射模型，用平均物面面元尺寸作为笛卡尔单元的尺寸，笛卡尔近似曲面显得十分“粗糙”，当弹身和尾舵之间缝隙较小时，图中所示的边缘单元已经错误地将缝隙填补了，从而导致挖洞过程中无法正确地判断缝隙附近的网格属性.图4是在最小洞映射模型上将笛卡尔网格单元的尺寸调整为原来的十分之一，此时缝隙已经能够得到很好地描述.

图4 弹身和尾舵最小洞映射加密模型

表1是对弹身模型分别采用传统方法的全模洞映射与最小洞映射方法构建的笛卡尔网格比较情况，加密指单元尺寸按计算的平均物面面元尺寸缩小至0.1倍，洞映射结果分别见图3和图4，其中全模加密后，因网格数量巨大，在内存为2 GB的微机上已无法实现.由表1可知，采用最小洞映射后，笛卡尔网格的数量下降了一个量级，极大减轻了需精细构建洞映射带来的存储压力.

表1中标记耗时一项是采用2.2节提出的新型洞映射单元标识方法对表中不同笛卡尔网格进行标识，在INTEL CORE2 Q8200的微机上运行，最小加密耗时仅1.26 s，其中全模加密因内存开销过大，无法计算.在2.2节中提到，对加密洞映射，洞映射单元相对网格尺度很小的情况，使用传统的标识方法，在标识尾舵时会发生错误.

表1 不同方法洞映射笛卡尔网格比较

2.2 洞映射单元属性标识

最小洞映射方法允许截断面作为洞映射边界是否合理，关键是能否正确地对洞映射单元属性进行识别.

传统洞映射方法在区分笛卡尔洞内、外单元时，要求物面必须封闭(即洞内单元和洞外单元区域分别单连通)，以此确定笛卡尔网格边缘上的单元为洞外单元(外流问题)，然后依此向网格中心递归推进.否则必须通过复杂的方法将物面封闭，或者由人工来指定区分，因此局限性很大.

广义封闭允许挖洞曲面结束于网格边界，使洞外单元被分割在多个孤立的区域，笛卡尔网格边缘单元的属性也不唯一，传统的标记方法就不再适用.目前广泛使用的是文献［8］提出的解决方案，通过计算物面面元的外法向矢量，经简单的逻辑判断，实现边缘单元邻侧洞映射单元属性的自动识别，原理如图5所示，再根据洞内(或洞外)单元的邻居也一定是洞内(或洞外)单元的原则，遍历洞映射集合内的所有单元，就能方便地判断出其余映射单元的属性.该方法只需对物面信息进行分析，在各个被分割的区域内以物面附近单元作为初始点进行递归推进，而无需考虑边界的情况，因此对本文扩展的广义封闭仍然适用.

图5 物面附近洞映射单元属性标记

但当笛卡尔网格尺寸远小于当地物面面元时，边缘单元两侧的洞映射单元属性判断容易模糊，上述方法就变得极不稳定了.图6是一个简单的二维尖劈外形，图中虚线部分为各物面面元笛卡尔包围盒(min，max)，笛卡尔映射单元P1与下表面面元包围盒唯一相交，被标记为边缘单元，P2是紧邻P1上侧的单元，将由下表面面元外法矢信息判断为洞内单元，显然这是一个错误的判断.在下一步推进中，由于P2须作为初始点向网格内部推进，可能会影响较大范围，甚至使程序陷入死循环.

图6 二维尖劈物面附近洞映射单元属性标记

上述问题在三维问题如尾舵等物面有拐折且网格没有连续过渡要求的结构上常见，本文2.1节研究的尾舵模型，在洞映射单元加密时就会发生类似的错误.即使不采用加密策略，由于传统方法计算的是全模的平均物面面元尺寸，在局部地方仍有可能出现上述问题.

如何能既避免可靠性低的使用物面信息判断的方法，又能实现方法与边界无关，以满足广义封闭，同时希望方法自动化程度高.针对这一问题，本文提出了一种新的特别适合于广义封闭的标识方法.

对广义封闭，洞映射单元属性识别的困难在于洞外单元被物面、网格边界和截断面分割在多个孤立的区域，如何保证在这些区域中都能找到合适的属性已明确的洞映射单元作为初始点进行推进是解决问题的关键，使用物面外法矢信息进行判断是其中一种策略.本文注意到这样一个事实:

准则1 洞映射笛卡尔网格的洞外单元必在计算网格中.

准则1的成立是十分明显的，洞映射的原理就是建立一个由主网格到挖洞网格的中间映射，其洞内单元对应主网格的洞内区域，洞外单元对应主网格的计算域.

对准则1的考虑是，广义封闭时，各个被分割的洞外区域里，洞外单元也必然是在主网格计算网格中，反之，计算网格也必然在洞映射的非洞内单元区域(洞外单元或边缘单元).如图7所示弹身最小洞映射模型的剖面网格，洞外点单元被分割在互不连通的两个区域内，按传统方法须在所有互不连通的区域人为地指定推进的初始点，这对复杂外形来说是比较繁琐的，依据准则1则十分简单，通过对计算网格进行分析和简单的下标计算，可以得到计算网格结点所在的笛卡尔洞映射单元，这些洞映射单元必然是洞外单元或边缘单元，且分布在各个分割区域.这一过程与用洞映射单元标记挖洞网格的洞内、外单元属性思路正好相反，故本文称新方法为“Inverse mark”.

图7 弹身最小洞映射模型剖面网格

本文建议在“Inverse mark”中使用计算网格的结点对笛卡尔洞映射单元进行标识，而不是使用更为常用的网格包围盒(min，max)，后者虽然可以完整地标识整个洞外区域，但在物面附近网格扭曲较大时，网格包围盒可能越过物面面元的包围盒，错误地对洞内单元进行标识.

“Inverse mark”使广义封闭的各个分割的洞外区域均有洞映射单元被标识为洞外单元，因此可以作为推进的初始点.本文认为以下准则成立:

准则2 洞映射笛卡尔网格的边缘单元连续，沿网格方向，边缘单元只能自封闭或终止于网格边界.

准则3 I、J、K方向上相邻的8个映射单元，不可能同时出现洞内单元和洞外单元.

准则2是对物面连续的要求，边缘单元一般由物面面元包围盒(min，max)标注，必然连续.依据准则2，笛卡尔网格的洞内单元与洞外单元分别多连通，洞外(洞内)单元的推进过程，只终止于边缘单元或网格边界，不会相互干扰，因此可以只使用一种如洞外单元进行推进.准则3是对相邻点原则的加强，认为洞外(洞内)单元的邻侧和对角单元也是洞外(洞内)单元.准则3为本文洞外单元的快速推进方法提供了理论依据.同时，准则3也能作为检查洞映射单元属性标记是否正确的标准.

基于上述3条准则，“Inverse mark”方法实现步骤如下:

1)按最小洞映射计算笛卡尔网格边界及网格单元尺寸，初始属性标记为临时单元;

2)计算物面面元包围盒(min，max)，标记笛卡尔网格边缘单元;

3)对计算网格进行遍历，若网格结点落在最小洞映射边界内，则分析该结点所在的笛卡尔网格单元，若该映射单元为临时单元，则标记为过渡单元;

4)对笛卡尔网格单元进行遍历，若有过渡单元，则标记为洞外单元，并检查其邻侧和对角单元，若有临时单元，则标记该临时单元为过渡单元，直至没有过渡单元;

5)遍历笛卡尔网格单元，剩余的临时单元全标记为洞内单元.

图8 “Inverse mark”方法标识洞映射单元

图8是“Inverse mark”方法实现的一个简单示意图，图8分别对应上述步骤2)～步骤5).由图8可知，初始推进点由计算网格产生，与边界或物面外法矢信息无关，因此完全满足广义封闭情况，同时可靠性较高.

3 计算方法

3.1 重叠网格方法

本文发展的最小洞映射和“Inverse mark”方法已在MI-GRID中得到应用.MI-GRID是北航阎超课题组研制的重叠网格软件，核心模块采用洞映射和割补法技术，包含了孤点清除、体积优化、壁面重叠、动态重叠等方法［8－9］.该重叠网格软件计算效率高、可靠性好、使用方便，先后参与了国内多个航空航天型号研制工作，包括复杂外形飞行器、子母弹抛撒、助推级分离、折叠翼打开等项目，得到了有效考核.

3.2 数值方法

本文算例均求解雷诺平均N-S控制方程，空间离散采用Roe的FDS(Flux Difference Splitting)格式，MUSCL(Monotone Upstream-centred Schemes for Conservation Laws)插值方法和Van Albada限制器用于获得二阶空间离散精度;湍流模型采用剪切应力输运(SST，Shea-Stress Transport)模型;时间离散采用稳定性高的LU-SGS(Lower-Upper Symmeffic Gauss-Seider)隐式计算方法.关于数值方法具体可参考文献［10］.

4 算例研究

4.1 TitanⅣ运载火箭

本文采用上述方法，对美国大力神四号(TitanⅣ)大型捆绑式运载火箭的超音速绕流问题［11］进行数值模拟研究.计算条件为:M∞=1.6，ReL=1.1×107，α=0°，由于流动条件对称，故使用半模计算.采用重叠网格方法分别生成芯级和助推级的计算网格，网格数目为190万和40万.

图9是建立的芯级和助推级最小洞映射模型，网格边界包括远场、对称面和截断面，笛卡尔单元数目共计925 580，在INTEL CORE2 Q8200的微机上运行，从建立洞映射到单元标识完成耗时0.45 s.洞映射结果表明，本文方法实现了对缝隙的精细刻画，标识方法在处理各种边界时可靠性较高.

图10是重叠结果和流场计算结果的示意，由于洞映射结果能准确描述芯级和助推级之间的缝隙，重叠结果正确，流场结构清晰合理.

图11是火箭芯级中心线上的压力分布与实验值的比较，计算结果与实验值吻合很好，因为在风洞实验中，芯级与助推级间存在连接机构，导致图11中实验的峰值点略高.

图9 TitanⅣ运载火箭最小洞映射模型

图10 对称面重叠网格和等马赫线流场图

图11 芯级中心线压力分布

4.2 子母弹抛撒模拟

子母弹抛撒的空气动力学问题是一类典型的超声速或高超声速多体干扰、非定常复杂流动问题.采用重叠网格方法求解子母弹抛撒问题，已成为国内外最常用的方法［12－13］.本文对典型外形的子母弹抛撒问题进行了研究，计算条件为:抛撒马赫数1.6，抛撒高度10 km，母弹攻角0°.母弹和8枚子弹分别生成计算网格，其中母弹网格数目为330万，子弹网格数目为35万.

图12是建立的母弹和子弹的最小洞映射模型，笛卡尔单元数目共计2945560，从建立洞映射到单元标识完成耗时1.25 s.母弹洞映射模型在弹仓端部被截断，“Inverse mark”方法能正确地进行标识.

图13、图14分别是对称面和某截面上重叠结果和流场图，图中重叠边界整齐、重叠形式合理，流场等值线在重叠区衔接较为光滑，流场结构清晰，说明本文使用的重叠网格方法和CFD求解方法对子母弹抛撒流场具有较强的解算能力.

图13 子母弹重叠网格示意图

5 结论

本文对传统的洞映射方法进行了研究，并针对传统方法存储大可靠性差的缺点提出改进，主要有:

1)对“广义封闭”的概念进行扩展，允许截断面作为笛卡尔网格边界，在此基础上，提出了最小洞映射方法.

2)最小洞映射方法有效缩小了洞映射区域，节省了存储的开销，有利于细化洞映射单元，避免了一些常见的由于洞映射过于粗糙引起的问题.

3)分析了不同洞映射单元标识方法的优劣，指出使用物面信息判断洞映射单元属性在笛卡尔网格较密时可靠性较差.

4)提出了一种新的特别适合广义封闭的标识方法:“Inverse mark”，利用计算网格对部分洞外单元进行自动标识，再作为初始点向网格内部递归推进.经算例验证，该方法可靠性好，自动化程度高.

5)提出广义封闭的洞映射笛卡尔网格满足的三条准则，其中准则3可以作为检查洞映射单元是否被正确识别的依据.

References)

［1］Ralph W N，Jeffrey P S.A summary of the 2004 overset symposium on composite grids and solution technology［R］.AIAA-2005-921，2005

［2］Benek J A，Buning P G，Steger J L.A 3-D Chimera grid embedding technique［R］.AIAA-85-1523，1985

［3］Bonet J，Peraire J.An alternating digital tree(ADT)algorithm for 3D geometric searching and intersection problems［J］.International Journal for Numerical Methods in Engineering，1991，31(1):1－17

［4］Chiu I T，Meakin R L.On automating domain connectivity for overset grids［R］.AIAA-95-0854，1995

［5］Meakin R.Object X-ray for cutting holes in composite overset structured grids［R］.AIAA-2001-2537，2001

［6］Sudharsun Jagannathan.A methodology for assembling overset generalized grids［D］.Mississippi:Computational Engineering in the College of Engineering，Mississippi State University，2004

［7］Sogers S E，Suhs N E，Dietz W E.PEGSUS 5:an automated pre-processor for overset-grid computational fluid dynamics［J］.AIAA Journal，2003，41(6):1037－1045

［8］李亭鹤.重叠网格自动生成方法研究［D］.北京:北京航空航天大学航空科学与工程学院，2004

Li Tinghe.Investigation of chimera grid automatic generation algorithm［D］.Beijing:School of Aeronautic Science and Engineering，Beijing University of Aeronautics and Astronautics，2004(in Chinese)

［9］范晶晶.复杂重叠网格方法研究及多体运动的非定常流动模拟［D］.北京:北京航空航天大学航空科学与工程学院，2010

Fan Jingjing.Enhancement of complex overset grid assembly and numerical simulation of unsteady multi-body movement［D］.Beijing:School of Aeronautic Science and Engineering，Beijing University of Aeronautics and Astronautics，2010(in Chinese)

［10］阎超.计算流体力学方法及应用［M］.北京:北京航空航天大学出版社，2006:245-248

Yan Chao.The methodology and application of computational fluid dynamics［M］.Beijing:Beihang University Press，2006:245－248(in Chinese)

［11］Taylor S，Johnson C T.Launch-vehicle simulations using a concurrent，implicit Navier-Stokes solver［R］.AIAA-95-0223，1995

［12］Spinetti R L，Jolly B A.Time-accurate numerical simulation of GBU-38s separating from the B-1B aircraft with various ejector forces，store properties，and load-out configurations［R］.AIAA-2008-187，2008

［13］张辉，重叠网格技术在飞行器多体分离问题中的应用［D］.北京:北京航空航天大学航空科学与工程学院，2010

Zhang Hui.The overlapping grid technology in the application of multi-body separation problem in the aerocraft［D］.Beijing:School of Aeronautic Science and Engineering，Beijing University of Aeronautics and Astronautics，2010(in Chinese)