國家數值風洞(NNW)通用軟件同構混合求解器設計

陳堅強, 馬燕凱, 閔耀兵, 趙 鐘, 何先耀, 何 琨

(中國空氣動力研究與發展中心 計算空氣動力研究所, 綿陽 621000)

0 引 言

計算流體力學(Computational Fluid Dynamics，CFD)，是當前飛行器設計過程中，用以獲得氣動數據的重要手段，已從早期的僅用于校驗輔助設計，發展到能為設計包線中的部分狀態直接提供氣動數據。然而，盡管近二三十年來，CFD發揮的作用在持續增加，但計算結果可信度僅局限于無大分離流動的巡航狀態(圖1)[1]，而在飛行包線的邊緣區，仍然存在較大的不確定度。

圖1 CFD在飛行包線中發揮的作用[1]

造成此種窘境的關鍵因素之一，在于當前使用“單一”的CFD方法模擬“復雜”的非線性流動：(1)單一的網格類型。不管是廣泛使用的In-House代碼，還是通用商業CFD軟件，都采用單一的網格類型，要么是基于結構網格、要么是基于非結構網格，相應的CFD求解器也只能是二選一。(2)單一的精度。當前工程中廣泛應用的CFD軟件，基本都是二階精度，受限于高階精度方法較低的健壯性，高階精度方法還難以在工程中成熟應用。(3)單一的學科。大多是求解僅考慮空氣動力學完全氣體效應的流動模擬。(4)單一的物理模型。鑒于直接數值模擬方法(DNS)和大渦模擬方法(LES)的巨量計算資源需求，工程中一般采用全湍流的雷諾平均NS方法(RANS)。

然而，現實中的空氣動力學卻是“多元”的：(1)多元的網格。相對結構網格，雖然因非結構網格生成難度的降低而應用越來越多，但在高超聲速、高精度等對網格質量要求較高的領域，由于結構網格具有計算效率高、對邊界層的模擬更精確、易于處理各向異性網格等優勢，依然有其存在的價值，二者仍將共存，且在各自領域發揮自身的優點。(2)多元的精度。當前，越來越多的高精度方法研究者意識到，即使在將來，高精度方法仍無法完全替代二階精度方法，而是互補的關系。高精度方法在大分離流動、強間斷等關鍵流動區域扮演關鍵角色，而二階方法在巡航等常規狀態模擬中具有較高的性價比。(3)多學科模擬。真實飛行環境中往往涉及“聲、光、電、磁、熱”等多學科耦合，單一的流體力學學科難以準確地模擬真實環境，毋庸置疑，未來的CFD必定是多學科耦合模擬。(4)多元的物理模型。由于復雜流動中伴隨層流、轉捩、湍流、大范圍分離等多種流動，采用單一的湍流模擬方法要么性價比低，要么難以準確捕捉流場結構，而在不同區域采用不同的湍流模擬方法，是一種現實可行的求解策略。

采用“混合求解器”，是解決當前的“單一”化與“多元”化需求矛盾的有效途徑之一。混合求解器，即在不同的區域，采用不同的求解器，發揮各自的優勢，以實現有工程價值的復雜CFD模擬。例如，不同網格類型的混合求解，不同精度的混合求解，不同學科的混合求解。

早在20世紀90年代，就有研究者實現了結構/非結構網格的Euler混合計算[2]，應用于飛機外形，通過在兩種類型網格的交界面上采用拼接網格，實現信息交換。2010年，研究者用混合計算方法模擬氣動聲學的多學科問題，在笛卡爾直角網格上采用有限差分法，在非結構網格上采用高階間斷有限元(DG)方法，以實現不大幅增加計算量的前提下，完成對復雜外形的高精度模擬[3-4]。知名In-House CFD代碼elsA-Hybrid，原本是一款采用多塊結構網格、基于有限體積方法的求解器elsA，為了充分利用結構網格高效高精度和非結構網格靈活性的優點，采用面向對象編程，在原有的結構求解代碼基礎上增加了非結構求解模塊，開發了新一代混合求解器elsA-Hybrid[5-6]，實現了三維高升力外形、翼身組合體等外形的二階混合模擬，兩種求解算法盡可能做到代碼的統一，例如在通量計算層次基本實現了代碼統一，基本實現了在細粒度層面的混合求解。在空間計算域融合方面，在作為CREATE-AV子項之一的Helios軟件系統中[7]，對旋翼類外形(圖2)，Helios在空間中采用不同的網格類型離散，分別調用不同類型的求解器，如在外場空間笛卡爾網格上調用SAMARC求解器，在旋翼附近調用結構求解器OverFLOW，在機身附近調用非結構求解器FUN3D，不同的求解器之間通過Python實現數據共享。作者所在團隊也在此前開展了結構/非結構混合計算相關研究[8-9]。

圖2 Create項目中的多求解器混合計算[7]

對于混合計算技術，目前通常有兩種不同的技術途徑。一種是“異構”混合計算，以多求解器間松耦合為特征，如上面提及的CREATE中的Helios分系統；另一種是“同構”混合計算，以多算法在同一個平臺上緊耦合模擬為特征，如上文提及的混合求解器elsA-Hybrid。“異構”混合計算能很方便地集成已有軟件，開發成本低，風險小。但因各求解器之間數據結構不同，并行框架不相容，重復模塊多，使其可擴展性、可維護性低。反之，“同構”混合計算是在同一個框架上開發，在數據結構、并行框架等方面都可以大量復用，因而效率高、可擴展性、可維護性強。但因需統一多種求解器，因而風險大、工作量大、難度大。“異構”混合計算模式適用于集成存量，而“同構”混合計算是一種面向下一代的CFD軟件的挑戰性模式。

“國家數值風洞(NNW)”是我國正在自主研發的CFD軟件系統，目標是建設國際領先的CFD高性能計算機系統，建成擁有自主知識產權、國內開放共享、達到世界一流水平的空氣動力數值模擬平臺。“風雷(PHengLEI)”軟件[10-11]，是NNW工程的通用CFD軟件，設計了面向結構/非結構的、適用于大規模并行的軟件架構，其典型特色是可實現任意求解器的“同構”混合計算，以解決未來多學科模擬中對“多元”化CFD模擬的需求。

WCNS高階精度有限差分方法由鄧小剛等人提出[12]并經過二十多年的持續發展，其對于流動結構的分辨能力較高[13]。在滿足幾何守恒律的基礎上[14-15]，WCNS方法在復雜網格計算中的適應能力得以大幅提升。近年來，在全局守恒條件約束下，鄧小剛等人又發展了與內點相匹配的WCNS邊界計算方法[16]。當前制約WCNS高階方法應用于工程復雜外形流動計算的主要難點在于高質量的多塊對接結構網格的生成，以及復雜邊界處理，其直接影響到計算魯棒性，是制約工程應用推廣的瓶頸。

本文就PHengLEI軟件中的混合求解器，介紹了混合求解策略及軟件設計方法，并通過典型算例，驗證結構/非結構網格耦合求解、二階有限體積/WCNS高階耦合求解這兩類混合求解方法的可行性。有關PHengLEI軟件的整體架構、并行方法，可參考文獻[10-11]。

1 域、交界面與數據交換

在正式介紹混合計算策略前，先介紹與混合計算緊密聯系的一些概念(如域、交界面)，以及相關的數據結構。

1.1 域和交界面的抽象

CFD模擬一般基于計算網格進行，整個計算域由多個計算區域組成，每個計算域被離散為若干網格單元。采用結構網格時，首先用線、面將整體計算域劃分為多個子域，再在每個子域生成網格單元。采用非結構網格時，雖然初始整體計算域作為一個整體存在，但是在并行計算時，一個計算域被分解為多個子域。我們將每個子域及其所屬的網格單元，定義為域(Zone)，不同域之間的交界面，定義為交界面(Interface)。圖3的示例中，整個計算域被劃分為2個域，Zone1和Zone2，為示普遍性，分別為非結構和結構網格，二者間是互不重疊且對接(點一一對應)的交界面。

圖3 域與交界面

域，是混合計算中不同求解器的載體。在混合計算時，將在每個域上存儲各自的求解器。在每一迭代步中，遍歷所有的子域求解。

交界面，是混合求解器的核心。不同域(其上可能加載不同的求解器)間的數據交換，通過交界面實施。實際上，這里的交界面是廣義的，包含三種類型：

1)多塊結構網格中，塊與塊之間的交界面；

2)非結構網格中，并行計算分區后，分區之間的交界面；

3)在計算域中同時含有結構、非結構網格時，兩種網格類型間的交界面。

混合求解，即是每個交界面兩側的兩個單元(分別是結構、非結構網格單元)交換流場數據，進而實現不同域之間耦合的過程。

1.2 數據結構與信息交換

PHengLEI軟件中定義的“交界面”是以上三種類型的合集，即將三種交界面抽象為統一的“Interface”。無論是何種類型，都采用相同的數據結構存儲、交換。因之，交界面的數據結構是核心。

為此，風雷軟件中設計了“三合一”的數據底層DataContainer，用于交換網格塊間的各類數據。任意數據經DataContainer標準化為數據底層后，在完全相同的機制下實現以下“三合一”數據交換：

a)進程內的結構多塊網格交界面數據交換；

b)進程間的相鄰網格塊(既可以是結構，也可以是非結構)交界面數據交換，主要用于MPI通信；

c)節點內任意網格塊交界面數據交換，主要用于混合求解器計算。

數據標準化的原理很簡單：對于結構、非結構兩種不同的數據存儲結構，將二者統一壓縮為字符型數據，以實現兩類數據的標準化[10]。

2 混合求解策略

下文以結構/非結構混合計算、二階/高階混合計算，說明混合求解器的耦合策略。

2.1 結構/非結構耦合策略

PHengLEI軟件框架既適應于結構網格、非結構網格，還能實現結構/非結構網格下的求解器耦合計算。結構/非結構兩種求解器可獨立運行(即在流場中同時含有結構網格和非結構網格的情況下，在結構網格上調用結構求解器，在非結構網格上調用非結構求解器)，并能進行結構求解器和非結構求解器的同步耦合計算。

結構/非結構耦合計算的關鍵技術是兩種類型網格交界面的數據通信。在兩種網格內部分別調用結構/非結構求解器計算，兩種網格間的交界面形成各自求解器的邊界條件，將這種邊界條件抽象為InterFace類。在結構網格和非結構網格之間，通過交界面建立一一對應的面映射關系(圖4)，計算模板相應擴展：對于結構求解器，在內場單元i處將另一側的非結構網格體心值視為其虛擬單元(i+1)；對于非結構求解器，在內場單元le處將另一側的結構網格體心值視為其虛擬單元(re)。對于結構求解器，往往需要在邊界處拓寬虛擬單元模板，除了i+1點外，還需要i+2，……，i+n點信息，為簡便起見，混合計算時一般采用一層虛擬單元。

(a)結構求解器模板 (b)非結構求解器模板

結構/非結構網格交界面InterFace的處理方式和并行分區交界面、多塊結構網格間交界面完全一樣，采用“三合一”數據底層進行數據交換。在計算過程中，通過交界面交換流場變量q和梯度q。

2.2 二階/高階精度耦合策略

PHengLEI軟件框架中，二階/高階精度耦合策略是指，二階精度求解器采用結構或非結構的有限體積方法，高階精度求解器采用基于結構網格的WCNS高階精度有限差分方法。研究表明[17]，WCNS求解器分辨率較二階有限體積求解器高，對關鍵流動結構的模擬能力明顯優于傳統的二階求解器。但WCNS高階精度計算方法對網格質量要求較高，其適應復雜外形計算的能力不如非結構網格的二階求解器。二階非結構解算器適用于復雜外形計算，但不具備WCNS高階精度求解器對流場結構的分辨率。因此，從WCNS求解器工程實用化觀點看，在關鍵流動特征區域采用高質量的結構網格基于WCNS求解器計算，而在其他復雜區域采用非結構網格基于二階精度非結構求解器計算，使高階精度結構求解器和二階精度非結構求解器優勢互補，不失為一種解決工程實際復雜問題的可行方案。

二階/高階精度求解器耦合策略與結構/非結構耦合策略相同，都是交換流場變量、梯度等(圖5)。

圖5 混合求解器信息交換

單一的WCNS求解器應用于多塊網格時，需根據精度的需要，交換緊鄰交界面的多層單元數據(三階精度需交換2層)。很顯然，非結構網格中由于沒有網格層，只向WCNS求解器提供1層信息，相當于在交界面處退化為二階精度。實際上，對于混合求解器來說，在非結構網格上已經全部是二階精度，即使在結構網格上有1層網格精度退化，也無大礙。

3 混合計算案例

為考察本文混合求解耦合策略的有效性，以二維圓柱和三維雙橢球的高超聲速黏性繞流為例，對二階非結構/二階結構混合計算、二階非結構/高階WCNS混合求解方法可行性與正確性進行驗證。

3.1 二維圓柱繞流

計算條件[18]：層流，來流馬赫數8.03，來流單位雷諾數1.835×105，來流溫度124.94 K，壁面溫度294.44 K。采用三套網格(圖6)進行計算：第一套為結構網格，網格數目為121×87(流線×法向)，采用結構WCNS高階精度求解器；第二套是混合網格，將第一套網格劃分為內外兩個塊得到，靠近壁面的是結構網格，采用WCNS高階精度求解器，遠離壁面的結構網格被非結構化，采用二階精度非結構有限體積求解器(Unstr-FVM)；第三套是混合網格，將第二套網格中外場塊四邊形網格三角化得到，內場結構網格上依然采用WCNS高階求解器，外場三角形網格采用Unstr-FVM求解器。

圖6 二維層流圓柱算例的網格分布及各區域所加載的求解器示意圖

這三套網格(三種計算策略)的目的是：第一套網格是全高精度計算，用于對比參考；第二套網格與第一套相同，僅僅是將四邊形非結構化，用于比較混合求解器帶來的差異；第三套是混合網格、混合求解器，是測試的目的。

從計算收斂性看(圖7)，三套網格殘差均收斂至機器“零”，且收斂后的氣動力系數保持一致。從物面壓力和熱流分布看(圖8)，三套網格上的結果幾乎重合，混合計算與全結構網格WCSN高階精度模擬結果相差無幾。該算例表明，混合計算既能保持高精度格式的方法，又適應非結構網格。

(a)殘差

(a)壁面壓力分布

3.2 雙橢球繞流

雙橢球繞流是典型的高超聲速模擬算例[19-21]。計算條件為：層流，0°攻角，來流馬赫數8.15，來流單位雷諾數1.67×107，來流溫度56 K，壁面溫度288 K，不考慮底部流動。采用兩套網格計算(圖9)：第一套為純結構網格，根據流場特點對強激波、激波與邊界層干擾、激波與激波干擾的區域進行了局部網格加密，網格總數為468萬，運行WCNS高階精度求解器；第二套是結構/非結構混合網格，壁面附近為結構網格，運行WCNS求解器，遠離壁面區域為非結構網格，運行二階精度非結構有限體積求解器(Unstr-FVM)。

(a)純結構網格

從殘差收斂看(圖10)，WCNS/FVM混合求解器能保持純WCNS高階求解器的收斂性。圖11是兩種求解器計算所得的表面極限流線，混合求解器保持了WCNS高階求解器對二次分離流的清晰捕捉能力。

圖10 混合求解器與純高階求解器殘差收斂對比

圖11 混合求解器與純高階求解器極限流線對比

從物面的壓強和熱流分布看(圖12)，兩種求解器計算得到的壓強分布與試驗數據有很好一致性，流經第二橢圓之后，壓力與熱流均顯著上升，達到峰值。同樣地，該算例表明了，在三維情況下，混合求解器仍然能保持高精度方法的特性。

(a)對稱面壓力分布

3.3 混合求解器效率比較

之所以發展“同構”混合求解，原因之一是其能吸收結構求解器的高效率優點。對上述兩個算例，對比分析了結構、非結構、混合求解三種策略的計算效率。對每個算例，三種策略都采用同一套結構網格，非結構求解器將其全局非結構化計算，混合求解器將外場部分非結構化計算。

在保證計算穩定收斂的前提下，三種策略均采用最大CFL數。在殘差收斂曲線圖中(圖13、圖14)，迭代步代表算法效率，墻上CPU時間則反應實際計算效率。從算法層面看(迭代步vs.殘差)，結構求解器效率最高，混合計算次之，非結構最慢；從實際墻上CPU時間看，規律一樣。即，兼顧了結構網格的計算精度和非結構網格的網格生成便利的混合求解策略，在計算效率上要優于完全的非結構網格計算。

(a)殘差收斂曲線(迭代步數)

(a)殘差收斂曲線(迭代步數)

4 結 論

基于風雷軟件，設計了“同構”混合求解策略，以及相應的交界面抽象方法、數據結構。作為“同構”混合計算的應用案例，以解決高階精度方法的工程實用化為切入點，通過在結構網格上運行WCNS高階求解器，在非結構網格上運行二階有限體積(FVM)求解器，以達到在流動關鍵區域高階模擬、在復雜區域非結構方便過渡的目的。

對二維圓柱和三維雙橢球高超聲速流動問題進行了模擬，結果表明，混合求解器能保持高階求解器的特性，驗證了風雷軟件中區域混合計算策略的可行性，同時說明了混合計算的效率優勢。

二階/高階混合求解的目標，是在對流動大范圍分離、激波等典型特征處采用高階方法，當前的工作僅是針對高超聲速流動，在高超聲速標模外形上進行的串行計算模擬，要實現在復雜外形上，高效、高精度模擬關鍵流動區域的最終目標，必須要實現混合求解器的并行化。風雷軟件中純結構、非結構求解器均有較強的并行可擴展性[10]，但要實現混合求解的并行計算，必須解決結構/非結構網格分區負載平衡問題。下一步，將著力解決求解器的并行化問題，以實現工程實用能力，為關鍵氣動問題模擬提供技術支撐。

致謝:感謝國家數值風洞工程軟件開發團隊所有人員的付出和努力，感謝風雷(PHengLEI)軟件前期開發者們的貢獻，也感謝工程兩委會的指導。