基于虛網格的格心ADT搜索法

袁 武 閻 超 于 劍 杜若凡

(北京航空航天大學國家計算流體力學實驗室，北京100191)

尋點方法是重疊網格［1－6］的關鍵技術之一，“尋點”是確定網格點在其他網格中的位置、查找網格貢獻單元的統稱.由于在重疊網格的生成過程中需要大量地進行“尋點”操作，因而準確與效率是尋點方法的關鍵.

常見的尋點方法有遍歷法、Stencil Walk［7］、反變換(Inverse Map)［8］、ADT(Alternating Digital Tree)搜索［9－10］等.其中，ADT 搜索方法用笛卡兒坐標下的包圍盒(min，max)對物體進行特征描述，由包圍盒是否相交來判斷物體的相交關系.由于包圍盒相交關系的判斷非常簡單，因此，ADT方法效率很高，且特別適合對幾何不規則目標進行搜索.

本文提出了一種直接在格心網格下操作、基于虛網格概念的ADT搜索方法，能有效提高尋點方法的效率和可靠性，并通過數值算例對新方法進行了驗證.

1 本文方法

1.1 格心ADT搜索

目前，基于格心網格的有限體積法在CFD(Computational Fluid Dynamics)中應用廣泛.在格心網格的尋點問題中，由于網格分區的關系，在分區邊界上，格心網格單元跨越多個子區，單元頂點位于不同子區內，難以給出格心網格單元的唯一位置標識，格心網格單元在邏輯空間上不封閉，這給ADT樹的建立帶來不便.一種解決辦法［11］是，仍然借用格點網格生成ADT樹，然后找到P點的格點貢獻單元，以該格點網格為頂點的格心網格必然比較靠近P點，再結合Stencil Walk方法，在附近找到真正的格心貢獻單元，但Stencil Walk方法存在局部收斂性要求，若初始位置單元不合理，Stencil walk方法有可能收斂不到合理貢獻單元而導致尋點失敗.

另一種思想就是直接在格心網格下進行ADT搜索，搜索結果是可能的格心貢獻單元集合，再利用點面位置關系對結果列表進行逐一判斷，即得到合理貢獻單元.不同于格點網格下搜索結果是解的附近某位置，格心網格下搜索結果是一個包含解的集合，故可以完全摒棄不穩定的Stencil Walk方法，直接對搜索結果列表進行遍歷，有效提高了尋點方法的可靠性.

格心網格下建立ADT樹的關鍵是使搜索空間完整封閉，問題描述如圖1所示，本文圖例均以實線表示格點網格線，虛線表示格心網格線.在分區邊界上，格心單元跨越多個子區，邏輯坐標不連續，無法直接組織ADT樹，導致格心單元搜索區域在計算域內不完整，不能正確搜索落入分區邊界附近的點，如圖中點P.若在原網格系統下封閉格心區域，數據結構復雜，且對各種分區邊界處理十分繁瑣.文獻［12］采用面、線、點順序延拓的方法，重構分區邊界的插值關系，但該方法只解決了完全封閉的計算域內部拓撲，對延拓部分沒有形成理論描述.

圖1 各子區格心單元覆蓋區域(虛線圍成陰影)

1.2 虛網格邊界

本文借鑒CFD計算格式對邊界處理的方法，引入虛網格的思想，較好解決了格心網格下ADT樹封閉問題.虛網格由計算域內格心網格向外擴展得到，圖2是二維的虛網格示意，與圖1情況比較，虛網格在分區邊界的擴展，填補了原格心單元各子區之間的“空白地帶”，確保搜索空間內部完整封閉.除分區邊界外，其余類型邊界也擴展了虛網格，使虛網格單元的實際覆蓋面積大于計算域，在下文中將介紹，非分區邊界虛網格不作為有效貢獻單元，僅在重疊過程中臨時使用.

圖2 二維虛網格示意

虛網格中包含邊界條件信息，如由遠場邊界外推的虛網格記為“遠場”，其余亦然.對邊界搭接的角點，如圖2中點(0，JN)，同時在遠場和對稱面上，則不標記邊界屬性.邊界附近尋點結果包含明確的邊界信息，對重疊過程中正確處理尋點結果有重要意義，如挖洞結束時，洞面點P1可能落在壁面和第1層格心網格之間，傳統方法搜索結果是無貢獻單元，且無具體信息，程序無法區分P1點與遠場外P2點，可能錯誤地解讀為洞邊界已遠離物面，使重疊出錯.由本文方法在計算域中(尤指邊界處格心域小于計算域的情況)均能直接找到頂點含虛網格的貢獻單元，狀態由虛網格邊界屬性確定，為重疊過程中正確處理尋點結果提供有效信息.

使用Stencil Walk方法尋點時，網格點在單一邊界附近有時也能提供邊界信息，但位于多個邊界搭接處時，Stencil Walk方法由初始位置出發，到第1個邊界處，不論是否找到貢獻單元都會停止移動，并記錄停止處單一邊界信息.如圖2中點P4，Stencil Walk方法搜索結果是無貢獻單元，邊界信息可能是“遠場”附近，也可能是“對稱面”附近，由移動過程不可預測地決定.而在重疊過程中，遠場和對稱面附近尋點結果意義是不同的.因此，Stencil Walk在邊界附近尋點存在隱患，可能存在合理貢獻單元但因邊界約束而無法正確查找，或者如前所述無法提供尋點結果的正確邊界信息.使用本文方法時，P4點可以直接找到頂點含虛網格的貢獻單元，虛網格包含“遠場”和“對稱面”信息.

1.3 虛網格構建

明確的邊界定義使計算虛網格物理坐標時可以根據其邊界屬性相應處理.遠場、壁面和對稱面等非分區邊界，虛網格不作為有效貢獻單元，無需計算準確的位置信息，可以直接由內層網格線性外插獲得.一些情況下外插的網格可能出現扭曲或變形，由于ADT方法對包圍盒搜索，只關注目標物理坐標上下限值，與結構的幾何細節無關，故虛網格的扭曲或變形對ADT的搜索結果沒有影響.

分區邊界附近計算跨區單元頂點是建立格心網格ADT搜索體系的難點.在兩個子區對接的簡單情況下，若虛網格有對應的實網格(相鄰分區第1層格心網格)，則使用該實網格位置信息以構成有效貢獻單元，且在重疊過程中虛網格屬性始終指向該實網格.

多個子區連接時，格心單元跨越多個子區，搜索單元頂點十分繁瑣，虛網格的定義較大簡化了處理方法.分兩種情況討論:

1)圖3a所示計算域內部有子區緊密連接時，由P1搜索對角線頂點P4，從P1出發，沿P1—P2—P3—P4或 P1—P5—P6—P4或其他任一條路徑(共有6種最短路徑)通過穿越分區邊界搜索相鄰實網格，均能獲得頂點P4的正確信息;

2)(二維示意可參考圖2中CELL1單元)如圖3b所示，跨區格心單元有非分區邊界頂點P6，此時沿P1—P5—P6—P4搜索，由于P6和P4邏輯坐標沒有關聯，搜索將中斷，只有沿特定路線如P1—P2—P3—P4，才能正確獲得頂點P4信息.

搜索路徑的選擇與算法有關，但對不同的拓撲情況和分區邏輯坐標順序不同，執行時有較大的不確定性.一種做法是對所有可能的路徑都進行分析，判斷跨區單元的虛網格頂點是否有對應的實網格，但這顯然不是一種好的做法.由本文方法對虛網格的定義，注意到頂點含非分區邊界虛網格的格心單元不是有效貢獻單元，該單元只使用邊界信息，不需要準確的位置信息.因此在搜索路徑中，若出現非分區邊界，跨區單元的該頂點可由線性插值近似獲得，無需嚴謹尋找或有的對應實網格，有效簡化了復雜分區邊界上跨區單元頂點搜索的算法.

圖3 邊界上搜索跨區單元頂點

1.4 對稱面處理和擴展的patch面

本文方法的一個重要應用是可以在對稱邊界上進行尋點，擴展了可處理的邊界類型.

先討論傳統方法處理對稱邊界遇到的困難，圖4是半模彈身在背景網格中的截面示意，模型中包括壁面、遠場和公共的對稱面邊界.對于格點系統，彈身網格和背景網格在對稱邊界上重合，格點單元均有效覆蓋半模空間.但在格心系統中，由于對稱面附近網格尺寸不同，兩物體對稱面格心網格線相貫，彈身格心單元的實際覆蓋區域小于背景網格，即意味背景網格在對稱面附近部分網格落在彈身網格域外，傳統方法尋點將找不到貢獻單元，不能作為插值點使用，如圖中格心網格線IS1上端位于彈身遠場外，是正常點，IS1下端被彈身物面挖去，是洞內點，而中間沒有可作重疊邊界的有效插值點，將導致計算中流場信息錯誤地在洞外區域和洞內區域傳遞，使計算結果錯誤甚至發散.

圖4 半模彈身在背景網格中示意

圖4所示情況只是對稱面尋點的一種，對多個物體重疊，或者對稱面附近網格相對關系變化大時，問題更加復雜，但原因均是各物體格心單元覆蓋的半模空間不完整，使對稱面附近ADT搜索區域有缺失，導致部分網格不能相互找點.虛網格的引入解決了上述問題，圖5是按本文方法擴展后的網格體系示意，在彈身和背景網格對稱面上各外推一層虛網格，由虛網格和實網格組成的單元集合完整地覆蓋了半模空間，使對稱面附近任一實網格均能搜索到貢獻單元.貢獻單元頂點中或包含虛網格，程序中作如下處理:

1)網格重疊過程中，只使用貢獻單元中實網格屬性作為洞邊界移動的判據，虛網格提供確保ADT搜索成功的完整空間;

2)重疊結束后建立插值邊界，若貢獻單元頂點含虛網格，則置換為相鄰的實網格單元，以與求解器匹配，此時插值點在實網格單元外部，插值系數在(0，1)區間外，但一般情況下對稱面附近網格匹配較好時，插值系數溢出很少.

圖5 半模彈身在背景網格中的虛網格擴展

與對稱面尋點類似，若有兩個邊界面貼合在一起，傳統方法的處理也存在尋點失敗的問題.圖6a是超音速來流條件下，彈身前體在背景網格中示意，模型中有公共的后場面，與上例類似的是后場上格點網格共面，格心網格位置交錯，使后場附近尋點搜索空間不完整.本文方法處理時，給后場網格面建立新的邊界類型“patch”，表示有網格面在此處貼合.圖6b是擴展后的網格體系，虛網格由相鄰實網格線性外推，用于封閉搜索空間，使尋點在后場附近能有效執行.程序中對頂點含“patch”虛網格的貢獻單元處理與對稱面情況一致.

圖6 前體在背景網格中的網格示意

2 算例研究

本文發展的尋點方法已在 MI-GRID［11，13-14］中得到應用.MI-GRID是北航閻超課題組研制的重疊網格軟件，從2000年開始研制，目前仍在不斷改進和發展.該重疊網格軟件計算效率高、可靠性好，先后參與了國內多個航空航天型號研制工作，包括復雜外形飛行器、子母彈拋撒、助推級分離、折疊翼打開等項目，得到了有效考核.

為驗證本文發展的方法，對美國大力神四號(TitanⅣ)大型捆綁式運載火箭的超音速繞流問題［15］進行數值模擬研究，計算條件為:M∞=1.6，ReL=1.1 ×107，α =0°.分別生成芯級和助推級的半模計算網格，網格數目為190萬和40萬.算例在配置INTEL CORE2 Q8200，2 GB內存的微機上運行，經測試，格點搜索法和格心搜索法生成ADT樹所需存儲開銷相當，但格點法在查找中需執行大量Stencil Walk操作，故執行時間要稍多于格心法.

圖7是TitanⅣ運載火箭半模模型示意，模型中包括公共的對稱面和后場patch面.圖8是重疊結果，在對稱面和patch面上重疊邊界整齊、重疊形式合理，說明本文方法在復雜邊界條件下能正確應用.

圖7 半模模型的邊界類型

圖9是火箭芯級中心線上的壓力分布與實驗的比較，計算結果與實驗值吻合較好.因為在風洞實驗中，芯級與助推級間存在連接機構，導致實驗的峰值點略高.

圖8 半模模型重疊網格示意

圖9 芯級中心線壓力分布

3 結論

本文提出了一種直接在格心網格下操作、基于虛網格邊界的ADT搜索方法，研究結果表明:

1)通過擴展虛網格，使格心網格覆蓋區域完整封閉，成功實現格心網格下的ADT數據結構.與格點搜索法相比，格心搜索法直接給出貢獻單元可能的解集，再進行局部遍歷，完全摒棄了不穩定的Stencil Walk方法，有效提高了尋點方法的可靠性.

2)虛網格邊界的定義使尋點在邊界附近的處理更為靈活，可以準確給出邊界附近貢獻單元的有效信息，同時，有效簡化了虛網格系統的構建，尤其是簡化了各種復雜拓撲下分區邊界附近跨區單元頂點的搜索算法.

3)虛網格的重要應用是擴展了重疊邊界類型，傳統方法處理對稱面時存在理論缺陷，新方法創建的搜索集合能完整覆蓋網格空間，解決了對稱面重疊問題.

References)

［1］Steger J L，Dougherty F C，Benek J A.A Chimera grid scheme［C］//ASME Mini-Symposium on Advances in Grid Generation.Houston:［s.n.］，1982:59 －69

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

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

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

［5］Prewitt N C，Belk D M，Shyy W.Parallel computing of overset grids for aerodynamic problems with moving grids［M］.Gainesville:［s.n.］，1999:117 －172

［6］Morgan P E，Visbal M R.Chimera-based parallelization of an implicit Navier-Stokes solver with applications［R］.AIAA-2001-1088，2001

［7］Belk D M，Maple R C.Automated assembly of structured grids for moving body problems［R］.AIAA-95-1680-CP，1995

［8］Meaking R L.A new method for establishing intergrid communication among systems of overset grids［R］.AIAA-91-1586，1991

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

［10］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

［11］李亭鶴.重疊網格自動生成方法研究［D］.北京:北京航空航天大學航空科學與工程學院，2004 Li Tinghe.Investigation of chimera grid au-tomatic generation algorithm［D］.Beijing:School of Aeronautic Science and Engineering，Beihang University，2004(in Chinese)

［12］王建濤.直升機旋翼前緣下垂控制動態失速數值模擬研究［D］.綿陽:中國空氣動力研究與發展中心研究生部，2008 Wang Jiantao.Numerical simulation of dynamic stall control on helicopter rotor wings with variable droop leading edge method［D］.Mianyang:Graduate School，China Aerodynamics Research and Development Center，2008(in Chinese)

［13］閻超.計算流體力學方法及應用［M］.北京:北京航空航天大學出版社，2006:200－220 Yan Chao.The methodology and application of computational fluid dynamics［M］.Beijing:Beihang University Press，2006:200－220(in Chinese)

［14］范晶晶.復雜重疊網格方法研究及多體運動的非定常流動模擬［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，Beihang University，2010(in Chinese)

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