洞映射方法的研究和改進

袁 武 閻 超 席 柯

(北京航空航天大學 國家計算流體力學實驗室，北京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)