NNW-HyFLOW高超聲速流動(dòng)模擬軟件框架設(shè)計(jì)

李鵬，陳堅(jiān)強(qiáng)，丁明松，何先耀，趙鐘，董維中,*

1. 中國空氣動(dòng)力研究與發(fā)展中心 空氣動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室，綿陽 621000

2. 中國空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所，綿陽 621000

在高超聲速(馬赫數(shù)Ma≥5)飛行條件下，飛行器繞流流場中將出現(xiàn)熱化學(xué)非平衡效應(yīng)(高溫氣體內(nèi)能模態(tài)激發(fā)與松弛，以及組分間的復(fù)雜化學(xué)反應(yīng)過程等[1])，對高超聲速飛行器的氣動(dòng)力特性、氣動(dòng)熱環(huán)境、流場結(jié)構(gòu)、光輻射特性、等離子體環(huán)境、電磁散射特性和電磁通信等均有重要影響[2]。受運(yùn)行成本和技術(shù)條件的制約，地面高焓測試設(shè)備可獲取的熱化學(xué)非平衡流場信息較為有限，難以實(shí)現(xiàn)真實(shí)飛行條件下對飛行器高溫氣體效應(yīng)的系統(tǒng)性研究。因此，常采用數(shù)值模擬方法對高超聲速飛行器的高溫氣體效應(yīng)進(jìn)行預(yù)測和評估。

自20世紀(jì)八九十年代以來，美國和歐盟等在高溫氣體非平衡效應(yīng)研究方面做了大量的理論和應(yīng)用研究，相關(guān)軟件的發(fā)展從未中斷，具備對高超聲速飛行器高溫非平衡流動(dòng)及其相關(guān)的氣動(dòng)力、氣動(dòng)熱環(huán)境以及等離子體分布等進(jìn)行數(shù)值模擬和分析的強(qiáng)大功能，例如美國NASA蘭利研究中心的LAURA[3]、美國Aerosoft公司的GASP[4]和GUST[5]、美國NASA蘭利研究中心HAPB小組開發(fā)的 VULCAN[6]、美國NPARC聯(lián)盟的WIND[7]以及法國ONERA的MSD[8]等。這些軟件功能十分完備，應(yīng)用面廣，但在國內(nèi)應(yīng)用受限，難以獲取技術(shù)支持。例如，GASP和GUST軟件可以供北美的公司和大學(xué)免費(fèi)使用，但是國內(nèi)用戶卻很難得到它的PDF形式的使用手冊。LAURA和VULCAN軟件同樣很難獲取，國內(nèi)多是通過公開發(fā)表的學(xué)術(shù)文獻(xiàn)了解它們的相關(guān)動(dòng)態(tài)。因此，國內(nèi)許多大學(xué)和研究單位也在發(fā)展自己的高溫非平衡流數(shù)值模擬軟件，如北京航空航天大學(xué)的ACANS[9-10]、山東大學(xué)的PHAROS[1-2]、中國空氣動(dòng)力研究與發(fā)展中心的AEROPH_Flow[11-12]和HyperFLOW[13-14]等。這些In-house軟件，在物理化學(xué)模型完整度、適應(yīng)復(fù)雜流動(dòng)狀態(tài)、多種物理機(jī)制耦合以及多種網(wǎng)格適應(yīng)性等方面，或多或少存在困難，很難得到有效推廣；其代碼低層次重復(fù)較多，可擴(kuò)展性差，不利于軟件的進(jìn)一步開發(fā)。

為了建設(shè)功能完善、適用面廣、可擴(kuò)展性強(qiáng)、自由共享的國產(chǎn)高超聲速流動(dòng)數(shù)值模擬軟件，在國家數(shù)值風(fēng)洞(National Numerical Windtunnel，NNW)工程[15]的支持下，中國空氣動(dòng)力研究與發(fā)展中心計(jì)算空氣動(dòng)力研究所聯(lián)合國內(nèi)相關(guān)單位，基于風(fēng)雷(PHengLEI)開源軟件[16-17]基礎(chǔ)框架開發(fā)以高溫氣體熱化學(xué)非平衡效應(yīng)及其氣動(dòng)力、氣動(dòng)熱和氣動(dòng)物理特性模擬分析為主要功能，具有自主知識產(chǎn)權(quán)的新型高性能工業(yè)應(yīng)用CFD軟件，即國家數(shù)值風(fēng)洞高超聲速流動(dòng)模擬軟件(NNW-HyFLOW)。

區(qū)別于國內(nèi)傳統(tǒng)In-house軟件，NNW-HyFLOW定位于國內(nèi)共享和自主創(chuàng)新，其研制將面向更廣的用戶群體，包括但不限于高校、研究機(jī)構(gòu)、相關(guān)工業(yè)部門和公司企業(yè)，滿足高超聲速流動(dòng)和高溫氣體熱化學(xué)非平衡效應(yīng)及其相關(guān)的氣動(dòng)力、氣動(dòng)熱、氣體輻射以及等離子體環(huán)境等特性數(shù)值模擬和綜合分析的需求。該軟件將具備如下特點(diǎn)：① 集成國內(nèi)外在化學(xué)反應(yīng)模型、熱力學(xué)模型、壁面催化、燒蝕、稀薄滑移、非平衡湍流模型、高效數(shù)值算法等方面的理論創(chuàng)新和先進(jìn)方法；② 實(shí)現(xiàn)網(wǎng)格架構(gòu)的多樣性、功能的完整性、代碼的可擴(kuò)展性和執(zhí)行高效性的協(xié)調(diào)統(tǒng)一。為了兼顧商業(yè)計(jì)算和科學(xué)研究，建立適用于結(jié)構(gòu)、非結(jié)構(gòu)、結(jié)構(gòu)/非結(jié)構(gòu)混合、重疊網(wǎng)格等多種網(wǎng)格類型的計(jì)算框架，滿足用戶對不同網(wǎng)格類型的偏好和選擇；為了擴(kuò)展復(fù)雜物理現(xiàn)象綜合評估的功能需求，設(shè)計(jì)了多種物理效應(yīng)清晰完備的接口定義，滿足非平衡解算器與多種物理化學(xué)模型和氣動(dòng)物理等多學(xué)科軟件耦合的計(jì)算需求；為了適應(yīng)不同物理化學(xué)模型和氣體屬性等對象的多樣性和高效計(jì)算需求，大量采用了面向?qū)ο箢惡透咝ㄓ玫亩嗑S數(shù)組等數(shù)據(jù)結(jié)構(gòu)設(shè)計(jì)，以減少代碼的低效重復(fù)，增加程序模塊的可擴(kuò)展性，提升軟件執(zhí)行的計(jì)算效率。軟件研發(fā)進(jìn)度安排為：2017—2023年完成軟件主體框架構(gòu)建和常規(guī)高溫氣體非平衡流動(dòng)模擬模塊開發(fā)，形成軟件1.0版本；2023—2027年將集成國內(nèi)相關(guān)優(yōu)勢單位在高溫非平衡效應(yīng)研究方面的理論創(chuàng)新和先進(jìn)方法，形成功能完善、適用面廣、可擴(kuò)展性強(qiáng)、高效的高超聲速流動(dòng)及其相關(guān)物理效應(yīng)的綜合模擬應(yīng)用軟件。目前已完成軟件的主體框架和數(shù)據(jù)結(jié)構(gòu)設(shè)計(jì)，形成軟件測試版本。

本文將主要介紹NNW-HyFLOW的總體設(shè)計(jì)框架、特點(diǎn)和軟件研發(fā)的相關(guān)進(jìn)展。首先從框架設(shè)計(jì)、數(shù)據(jù)結(jié)構(gòu)、耦合算法、并行計(jì)算以及接口設(shè)計(jì)等方面分別陳述軟件的主要設(shè)計(jì)框架及特點(diǎn)。其次，對求解器解算模型及其理論方法進(jìn)行簡要介紹，包括模型方程、物理化學(xué)模型以及核心數(shù)值算法等。最后，選取典型的熱化學(xué)非平衡流算例開展模擬測試，檢驗(yàn)當(dāng)前軟件測試版本的正確性，以及對復(fù)雜工程外形熱化學(xué)非平衡流場數(shù)值模擬的適用性和可靠性。

1 軟件設(shè)計(jì)

NNW-HyFLOW軟件基于風(fēng)雷開源軟件(風(fēng)雷軟件或PHengLEI)的基礎(chǔ)框架進(jìn)行設(shè)計(jì)開發(fā)，結(jié)合“黑匣子”原理，實(shí)現(xiàn)軟件各模塊的封裝，以適應(yīng)多種網(wǎng)格類型、多種物理機(jī)制、多種物理化學(xué)模型以及多種流動(dòng)狀態(tài)的模擬需求，因而具備適用性廣、功能完整、擴(kuò)展性好等優(yōu)點(diǎn)。以下從框架設(shè)計(jì)、數(shù)據(jù)結(jié)構(gòu)、耦合算法、并行計(jì)算以及接口設(shè)計(jì)等方面，詳細(xì)介紹軟件設(shè)計(jì)框架及特點(diǎn)。

1.1 框架設(shè)計(jì)

作為計(jì)算流體力學(xué)(CFD)軟件，計(jì)算框架是NNW-HyFLOW軟件設(shè)計(jì)的核心。如圖1所示，NNW-HyFLOW的計(jì)算框架由前處理、求解器和后處理等3個(gè)系統(tǒng)模塊構(gòu)成。其中，前處理負(fù)責(zé)前臺數(shù)據(jù)接收和轉(zhuǎn)換，包括不同網(wǎng)格類型的格式轉(zhuǎn)換、并行計(jì)算要求的網(wǎng)格分區(qū)轉(zhuǎn)換以及化學(xué)模型參數(shù)讀取與轉(zhuǎn)換等。基于前置網(wǎng)格處理模塊可將外部商業(yè)軟件生成的網(wǎng)格格式，如Gridgen的grd、Pointwise的cgns、Fluent的cas以及GMSH的msh等，轉(zhuǎn)換為風(fēng)雷開源軟件定義的混合格式，可滿足不同類型網(wǎng)格的耦合計(jì)算需求。另外，前處理模塊基于文獻(xiàn)[17]所述的網(wǎng)格分區(qū)方法實(shí)現(xiàn)均衡的網(wǎng)格分區(qū)劃分，能滿足大規(guī)模并行計(jì)算的需要。與風(fēng)雷開源軟件的最大不同，NNW-HyFLOW前處理增加了化學(xué)模型參數(shù)讀取與轉(zhuǎn)換功能模塊，可以滿足對不同類型化學(xué)反應(yīng)模型(包括自定義化學(xué)模型)的計(jì)算需求。

NNW-HyFLOW求解器包含黏性模型、化學(xué)模型、氣體模型、壁面條件、數(shù)值格式以及控制條件等功能模塊。為避免底層代碼低效的重復(fù)利用，方便程序的進(jìn)一步拓展與維護(hù)，其求解器繼承了風(fēng)雷開源軟件面向?qū)ο蟮木幊趟枷搿⒍鄬哟蔚臄?shù)據(jù)結(jié)構(gòu)及其在并行底層、邊界處理和結(jié)果顯示方面的優(yōu)勢，采用分層模塊設(shè)計(jì)思想進(jìn)行開發(fā)。如圖2所示，各功能模塊按照內(nèi)核層、算法層、功能層和應(yīng)用層進(jìn)行分層構(gòu)建，其中非平衡求解器基于已有的Navier-Stokes(N-S)求解器進(jìn)行擴(kuò)充和完善，在數(shù)值算法層面增加和封裝多種物理化學(xué)模型及函數(shù)，在邊界條件中增加多種非平衡流動(dòng)特定的壁面條件，如輻射平衡溫度、壁面催化、稀薄滑移、壁面材料燒蝕等條件類型，可以滿足對多種物理機(jī)制、多種物理效應(yīng)計(jì)算模擬的需要。為了兼容原有的完全氣體CFD求解器，在數(shù)值算法層面、功能層面以及應(yīng)用層面擴(kuò)展相應(yīng)模塊，通過類的重載、繼承和派生及其熱化學(xué)流動(dòng)空間、模板的定義，實(shí)現(xiàn)數(shù)據(jù)結(jié)構(gòu)拓展和功能性擴(kuò)充。

圖2 功能模塊分層設(shè)計(jì)示意圖Fig.2 Hierarchical design of major functional modules

NNW-HyFLOW后處理主要用于氣動(dòng)力、熱環(huán)境以及等離子體環(huán)境等非平衡流場特性的輸出和顯示?；陲L(fēng)雷開源軟件的后置處理模塊，NNW-HyFLOW增加了非平衡流動(dòng)相關(guān)的結(jié)果處理函數(shù)，包括但不限于非平衡流場特性、氣動(dòng)力特性、氣動(dòng)熱環(huán)境、等離子體分布特性等參數(shù)的提取和輸出功能。同時(shí)增加了與氣體輻射場、電磁場、固體溫度場等模塊耦合計(jì)算的接口模塊，能夠滿足多種物理機(jī)制耦合效應(yīng)模擬的需要。

NNW-HyFLOW對標(biāo)LAURA軟件進(jìn)行開發(fā)，其在框架設(shè)計(jì)上幾乎覆蓋了LAURA軟件的全部功能點(diǎn)，以LAURA[3]版本為例，兩者具備的主要功能模塊如表1所示。

表1 NNW-HyFLOW與LAURA的功能對比

1.2 數(shù)據(jù)結(jié)構(gòu)

從類結(jié)構(gòu)設(shè)計(jì)和數(shù)組兩個(gè)方面分別介紹NNW-HyFLOW軟件的底層數(shù)據(jù)結(jié)構(gòu)。

NNW-HyFLOW采用面向?qū)ο缶幊趟枷朐O(shè)計(jì)和封裝各功能模塊，利用C++類的重載、繼承和派生等多態(tài)特性實(shí)現(xiàn)拓展應(yīng)用范圍的目的。例如針對不同求解器采用了如圖3所示的類設(shè)計(jì)，其中CFDSolver是所有基于CFD方法求解器的父類。為了方便程序的編寫和維護(hù)，N-S方程和RANS方程采用解耦計(jì)算策略，因此由父類CFDSolver派生出基于N-S方程解算的N-S求解器(NSSolver)和基于RANS方程解算的湍流求解器(TurbSolver)。為了滿足結(jié)構(gòu)/非結(jié)構(gòu)網(wǎng)格耦合計(jì)算的需要，再通過基類NSSolver和TurbSolver分別派生出基于結(jié)構(gòu)網(wǎng)格解算的求解器(NSSolverStruct與TurbSolverStr)和基于非結(jié)構(gòu)網(wǎng)格解算的求解器(NSSolverUnStruct與TurbSolverUnstr)。

圖3 求解器類結(jié)構(gòu)設(shè)計(jì)示意圖Fig.3 Design diagram of solver class structure

NNW-HyFLOW采用完全耦合的策略將N-S 方程與組分輸運(yùn)/振動(dòng)能量/電子能量方程聯(lián)合解算，因此在CFDSolver基類下不再獨(dú)立為非平衡求解器派生設(shè)計(jì)新類，非平衡求解功能仍基于NSSolver類進(jìn)行開發(fā)，相關(guān)代碼在NSSolverStruct和NSSolverUnStruct中擴(kuò)充和添加。另外，在Solver基類下派生了氣動(dòng)物理求解器類(AeroPHSolver)、有限元求解器類(FESolver)等與氣體輻射場、固體溫度場等計(jì)算相關(guān)的其他求解器類型，可作為非平衡流場與氣動(dòng)物理、結(jié)構(gòu)力學(xué)等解耦計(jì)算的預(yù)留接口。在CFDSolver基類下派生了磁流體求解器類(MHDSolver)，可以實(shí)現(xiàn)磁流體相關(guān)計(jì)算需求。

為適應(yīng)完全氣體、平衡氣體以及化學(xué)非平衡和熱化學(xué)非平衡氣體等多種計(jì)算氣體模型的一體化計(jì)算需求，NNW-HyFLOW采用了圖4所示的類結(jié)構(gòu)設(shè)計(jì)。在氣體模型基類(Gas)下派生出完全氣體(PerfectGas)和真實(shí)氣體(RealGas)兩個(gè)主要子類，其中完全氣體又衍生出量熱完全氣體(CaloricallyPerfectGas)和熱完全氣體(ThermalPerfectGas)，原有的完全氣體計(jì)算功能被量熱完全氣體類涵蓋，而熱完全氣體類則作為接口用于實(shí)現(xiàn)平衡氣體計(jì)算功能。另外，真實(shí)氣體類主要用于熱力學(xué)能量激發(fā)以及化學(xué)反應(yīng)氣體模型的計(jì)算，因此根據(jù)模擬大氣環(huán)境和化學(xué)反應(yīng)類型的差異，真實(shí)氣體類下又派生了地球反應(yīng)氣體(EarthReactionGas)、火星反應(yīng)氣體(MarsReactionGas)以及燃燒反應(yīng)氣體(CombustionGas)等子類，分別用于模擬地球大氣環(huán)境、火星大氣環(huán)境以及各種燃?xì)夥磻?yīng)體系。

圖4 氣體模型類結(jié)構(gòu)設(shè)計(jì)示意圖Fig.4 Design diagram of gas model class structure

上述類結(jié)構(gòu)設(shè)計(jì)為NNW-HyFLOW軟件功能擴(kuò)展提供了靈活方便的接口，充分滿足了多種計(jì)算氣體模型、多種物理化學(xué)效應(yīng)一體化計(jì)算的需要。NNW-HyFLOW與LAURA軟件在程序設(shè)計(jì)上的最大不同在于編程規(guī)范。采用C++面向?qū)ο缶幊趟枷胧沟肗NW-HyFLOW在功能擴(kuò)展方面更加靈活方便。LAURA軟件的最新版本使用Fortran 95語言規(guī)范，同樣也具備了面向?qū)ο蟮木幊趟枷?，能夠設(shè)計(jì)派生數(shù)據(jù)類型以適應(yīng)復(fù)雜數(shù)據(jù)結(jié)構(gòu)構(gòu)建的需要，但程序設(shè)計(jì)的主體仍然是面向過程，因此在底層代碼復(fù)用和功能拓展方面，它并不如NNW-HyFLOW易用和方便。

另外在變量存儲方面，NNW-HyFLOW僅對原有的數(shù)組空間進(jìn)行拓展，對其維度不做變更，可保持原有數(shù)組數(shù)據(jù)結(jié)構(gòu)的高效性。采用圖5所示的數(shù)組拓展方式設(shè)計(jì)底層數(shù)據(jù)結(jié)構(gòu)，兼顧了完全氣體、平衡氣體、單溫度化學(xué)非平衡氣體、兩溫度熱化學(xué)非平衡氣體、三溫度熱化學(xué)非平衡氣體以及多振動(dòng)溫度模型等多種計(jì)算氣體模型的一體化計(jì)算需要。相比于LAURA軟件，計(jì)算功能更加全面，適用范圍廣。

圖5 數(shù)組擴(kuò)展方式示意圖Fig.5 Diagram of array expansion method

1.3 耦合算法

基于圖3所示的類定義和設(shè)計(jì)，NNW-Hy-FLOW利用Solver/CFDSolver等父類對象的虛函數(shù)定義及多態(tài)特性實(shí)現(xiàn)不同類型求解器在同一控制流程下的合作計(jì)算，滿足結(jié)構(gòu)/非結(jié)構(gòu)網(wǎng)格耦合解算，以及非平衡求解器與其他求解器的耦合計(jì)算等要求。如圖6所示，在NNW-HyFLOW框架下不同求解器耦合計(jì)算的控制流程包含3層循環(huán)：

圖6 耦合計(jì)算控制流程示意圖Fig.6 Diagram of coupled computing control process

1) 第1層循環(huán)：遍歷每個(gè)迭代步，執(zhí)行CFD計(jì)算的一次時(shí)間推進(jìn)。

2) 第2層循環(huán)：遍歷不同求解器類型，執(zhí)行一類求解器的所有計(jì)算操作。不同求解器通過內(nèi)部接口定義實(shí)現(xiàn)相關(guān)數(shù)據(jù)共享。

3) 第3層循環(huán)：遍歷每個(gè)分區(qū)網(wǎng)格，執(zhí)行對應(yīng)求解器類型指定的計(jì)算流程。該層循環(huán)結(jié)束時(shí)，執(zhí)行并行通信，實(shí)現(xiàn)不同分區(qū)網(wǎng)格交接面的數(shù)據(jù)交換。

在求解器內(nèi)進(jìn)行并行通信的設(shè)計(jì)，使得開發(fā)人員只需關(guān)注某一類求解器的并行通信數(shù)據(jù)，卻不用考慮其他求解器的通信要求，從而給不同求解器耦合計(jì)算程序的開發(fā)提供便利。

由于LAURA軟件僅應(yīng)用于結(jié)構(gòu)網(wǎng)格，因此它并沒有上述耦合算法設(shè)計(jì)框架。NNW-HyFLOW的耦合算法設(shè)計(jì)不僅滿足了不同類型網(wǎng)格一體化計(jì)算的需要，同時(shí)給求解器層面的耦合計(jì)算提供了方便易用的接口，進(jìn)一步增強(qiáng)了軟件的可擴(kuò)展性。

1.4 并行計(jì)算

NNW-HyFLOW的并行算法設(shè)計(jì)如圖7所示，將求解器列表捆綁到分區(qū)上，每個(gè)分區(qū)對應(yīng)一個(gè)進(jìn)程，且每個(gè)分區(qū)里包含多塊網(wǎng)格，其并行計(jì)算過程同樣見圖6，當(dāng)一種求解器類型的計(jì)算執(zhí)行完成后，每個(gè)分區(qū)之間執(zhí)行一次并行通信，同時(shí)在分區(qū)內(nèi)部網(wǎng)格間執(zhí)行一次交接面數(shù)據(jù)交換。當(dāng)求解器列表遍歷完成后，也要進(jìn)行求解器之間的數(shù)據(jù)通信，以實(shí)現(xiàn)不同求解條件下流場或特殊邊界的數(shù)據(jù)傳遞和交換，滿足多物理場耦合計(jì)算的需要。

圖7 并行分區(qū)示意圖Fig.7 Diagram of parallel regions

在數(shù)據(jù)通信層面，NNW-HyFLOW直接采用了風(fēng)雷開源軟件提供的并行通信接口，其核心步驟包括數(shù)據(jù)壓縮、MPI通信和數(shù)據(jù)解壓[17]，如圖8所示。為適應(yīng)不同類型網(wǎng)格耦合并行計(jì)算的需要，NNW-HyFLOW同樣采用了風(fēng)雷開源軟件設(shè)計(jì)的“數(shù)據(jù)容器”[16](DataContainer)結(jié)構(gòu)，專門用于MPI層面的各種數(shù)據(jù)傳遞和交換，滿足差異化數(shù)據(jù)的統(tǒng)一存取和并行通信。這種并行通信接口設(shè)計(jì)大大簡化了并行計(jì)算的設(shè)計(jì)難度，提高了軟件開發(fā)的效率。

圖8 并行計(jì)算通信接口示意圖Fig.8 Diagram of parallel communication interface

目前，NNW-HyFLOW已經(jīng)完成了結(jié)構(gòu)網(wǎng)格CFD求解器的并行計(jì)算代碼實(shí)現(xiàn)，其根據(jù)圖5所示的數(shù)組退化或拓展方式對網(wǎng)格交接面的密度、速度、壓強(qiáng)、溫度以及組分質(zhì)量分?jǐn)?shù)等相關(guān)數(shù)據(jù)進(jìn)行封裝，然后調(diào)用并行通信接口函數(shù)進(jìn)行推送即完成并行計(jì)算代碼實(shí)現(xiàn)?；陲L(fēng)雷開源軟件提供的并行通信接口，NNW-HyFLOW軟件的并行設(shè)計(jì)具有簡單易用、擴(kuò)展性強(qiáng)等優(yōu)勢，這是LAURA軟件不具備的能力。

1.5 接口設(shè)計(jì)

與氣動(dòng)物理、流固耦合等NNW多學(xué)科軟件的耦合計(jì)算主要采用完全解耦策略。因此，NNW-HyFLOW軟件在后處理模塊中預(yù)留了相關(guān)輸出函數(shù)接口，將根據(jù)其他軟件的格式要求輸出流場或邊界面的相關(guān)參數(shù)。同時(shí)，在應(yīng)用層面同樣預(yù)留了相關(guān)求解器接口，當(dāng)其他軟件的求解器有基于風(fēng)雷開源平臺進(jìn)行程序開發(fā)的需要時(shí)，可以快速開展代碼移植和求解器耦合計(jì)算?？焖俅a移植和便利的耦合計(jì)算程序設(shè)計(jì)是NNW-HyFLOW區(qū)別于LAURA軟件的又一項(xiàng)優(yōu)勢。

NNW-HyFLOW軟件的最大功能特色是熱化學(xué)非平衡效應(yīng)模擬。針對復(fù)雜的化學(xué)模型計(jì)算，NNW-HyFLOW設(shè)計(jì)了專門的函數(shù)接口用于處理不同類型的化學(xué)反應(yīng)模型。類似“數(shù)據(jù)容器”的設(shè)計(jì)方式[16]，NNW-HyFLOW采用“化學(xué)模型容器”對內(nèi)置化學(xué)模型或其他類型的化學(xué)反應(yīng)模型進(jìn)行標(biāo)準(zhǔn)化統(tǒng)一管理。如圖9所示，“化學(xué)模型容器”包含3類主要信息：基本信息、化學(xué)反應(yīng)式、反應(yīng)速率模型?；拘畔ńM分總數(shù)、反應(yīng)式總數(shù)，以及各組分的化學(xué)名稱?；瘜W(xué)反應(yīng)式部分以3個(gè)二維數(shù)組的方式分別列舉和存儲各個(gè)組分在每個(gè)化學(xué)反應(yīng)式的前、后反應(yīng)中的化學(xué)計(jì)量系數(shù)，以及作為第三碰撞體的系數(shù)。反應(yīng)速率模型部分將列舉各反應(yīng)式前、后反應(yīng)式速率的擬合參數(shù)形式，這些參數(shù)同樣采用二維數(shù)組的方式存儲。

圖9 化學(xué)模型容器及接口設(shè)計(jì)示意圖Fig.9 Design diagram of chemical model container and its interfaces

當(dāng)前，NNW-HyFLOW軟件通過“化學(xué)模型容器”內(nèi)置和提供了Gupta、Dunn-Kang和Park等3類常用化學(xué)反應(yīng)模型，其接口函數(shù)能夠根據(jù)用戶輸入的氣體組分列表自動(dòng)匹配相關(guān)反應(yīng)式，同時(shí)構(gòu)建用于當(dāng)前計(jì)算的最小集合化學(xué)模型，實(shí)現(xiàn)化學(xué)反應(yīng)模型的自動(dòng)化設(shè)置。這一設(shè)計(jì)提高了軟件的計(jì)算效率，同時(shí)也大大增強(qiáng)了軟件的易用性和靈活性。

另外，用戶利用“化學(xué)模型容器”也可自定義化學(xué)反應(yīng)模型，其輸出接口函數(shù)可以自動(dòng)生成參考化學(xué)反應(yīng)模型的模板文件，用戶可利用界面程序?qū)ο嚓P(guān)參數(shù)進(jìn)行修改，或直接打開模板文件修改有關(guān)參數(shù)，計(jì)算時(shí)再利用輸入接口函數(shù)導(dǎo)入“化學(xué)模型容器”即可用于后續(xù)流場計(jì)算。LAURA軟件同樣具備化學(xué)反應(yīng)模型自定義功能，但在易用性、靈活性和通用性方面，NNW-HyFLOW的自定義操作更具優(yōu)勢。

2 理論方法

基于前述軟件框架的描述，以下僅從控制方程、物理化學(xué)模型以及核心數(shù)值方法等方面，對NNW-HyFLOW軟件的解算模型及其理論方法進(jìn)行簡單介紹。

2.1 控制方程

NNW-HyFLOW軟件主要基于Navier-Stokes控制方程組，通過對完全氣體層流控制方程的繼承、派生與擴(kuò)充，實(shí)現(xiàn)二維/三維熱化學(xué)非平衡流控制方程體系[1-2,5,11-12]的構(gòu)建。其最核心的控制方程形式為

(1)

式中：Q為守恒量向量；F、Fv分別為對流項(xiàng)通量和黏性項(xiàng)通量；S表示源項(xiàng)。

高溫氣體化學(xué)反應(yīng)和熱力學(xué)激發(fā)及其非平衡過程的表征，主要通過非平衡源項(xiàng)S，與流動(dòng)方程中各氣體組分質(zhì)量、動(dòng)量和能量的時(shí)間演化與空間輸運(yùn)過程相聯(lián)系。通過數(shù)據(jù)結(jié)構(gòu)和函數(shù)的重載，實(shí)現(xiàn)多種氣體模型和控制方程的統(tǒng)一。對于完全氣體，Q=[ρ,ρu,ρv,ρw,ρE]T，S=0；對于單溫度化學(xué)非平衡氣體，Q=[ρ,ρu,ρv,ρw,ρE,ρi]T,S=[0, 0, 0, 0, 0,Θi]T；對于兩溫度熱化學(xué)非平衡氣體，Q=[ρ,ρu,ρv,ρw,ρE,ρi,ρev]T,S=[0, 0, 0, 0, 0,Θi,Θv]T；對于三溫度熱化學(xué)非平衡氣體,Q=[ρ,ρu,ρv,ρw,ρE,ρi,ρev,ρee]T,S=[0, 0, 0, 0, 0,Θi,Θv,Θe]T；對于多振動(dòng)溫度熱化學(xué)非平衡模型[12],Q=[ρ,ρu,ρv,ρw,ρE,ρi,ρevi,ρee]T,S=[0, 0, 0, 0, 0,Θi,Θvi,Θe]T。這里ρ為氣體密度，u、v、w分別為直角坐標(biāo)系x、y、z方向速度；E、ev、ee分別為氣體的內(nèi)能、振動(dòng)能和電子能，對于不同氣體模型，其物理實(shí)質(zhì)、涵蓋范圍及表達(dá)形式存在差別；ρi和evi分別為第i個(gè)氣體組分的密度和振動(dòng)能；Θi、Θv、Θe分別為化學(xué)反應(yīng)生成源項(xiàng)、振動(dòng)非平衡源項(xiàng)和電子非平衡源項(xiàng)，Θvi為第i個(gè)組分的振動(dòng)非平衡源項(xiàng)。

控制方程中組分ρi的質(zhì)量連續(xù)性方程可以采用多種模式[18]：

1)i=1, 2,…,n。當(dāng)氣體組分不包含離子時(shí),n=ns-1(ns為總的氣體組分?jǐn)?shù))，結(jié)合ρ=∑ρi計(jì)算得到ρns；當(dāng)氣體組分包含離子時(shí),n=ns-2，進(jìn)一步結(jié)合電中性原則計(jì)算得到ρns-1。這樣處理的優(yōu)點(diǎn)在于可較好地繼承完全氣體方程形式及其軟件模塊，缺點(diǎn)在于相關(guān)非平衡偏導(dǎo)數(shù)形式較為復(fù)雜。

2)i=1, 2, … ,ns，不計(jì)算質(zhì)量連續(xù)性方程，由ρ=∑ρi計(jì)算得到。該方法對于完全氣體模塊的繼承性稍差。

3)i=1, 2, …,ns，同時(shí)計(jì)算ρ的連續(xù)性方程。其優(yōu)點(diǎn)在于較好的繼承完全氣體軟件模塊，同時(shí)非平衡偏導(dǎo)數(shù)形式相對簡單；缺點(diǎn)在于計(jì)算量較大，同時(shí)方程超定，需要額外進(jìn)行質(zhì)量歸一化和電中性修正。

由于NNW-HyFLOW軟件將兼顧多種物理效應(yīng)耦合的模擬需求，因此除核心的非平衡控制方程組外，還需耦合考慮湍流控制方程、氣體輻射方程、電磁麥克斯韋方程以及固體溫度場熱傳導(dǎo)方程等控制方程。為滿足多物理效應(yīng)模型的多樣化需求，這些控制方程采用各自獨(dú)立的構(gòu)架和耦合接口設(shè)計(jì)。

2.2 物理化學(xué)模型

在NNW-HyFLOW軟件框架上，完全氣體相關(guān)的熱力學(xué)特性參數(shù)在圖4所示的PerfectGas類中進(jìn)行包裝和實(shí)現(xiàn)，而熱化學(xué)非平衡計(jì)算相關(guān)的物理、化學(xué)模型及其相關(guān)參數(shù)則在RealGas類中進(jìn)行封裝和實(shí)現(xiàn)。以下分別對熱力學(xué)特性參數(shù)、化學(xué)反應(yīng)生成源項(xiàng)及化學(xué)反應(yīng)模型、振動(dòng)-電子能量方程源項(xiàng)的計(jì)算方法進(jìn)行簡要介紹。

完全氣體條件下的熱力學(xué)狀態(tài)參數(shù)采用理想氣體狀態(tài)方程計(jì)算。對于熱化學(xué)非平衡氣體，氣體壓強(qiáng)采用Dalton分壓定理[1]計(jì)算，各組分的熱力學(xué)狀態(tài)參數(shù)仍由理想氣體狀態(tài)方程計(jì)算。對于比熱、內(nèi)能以及焓等其他熱力學(xué)特性參數(shù)，熱力學(xué)溫度模型采用分子理論建模和計(jì)算，組分i的平-轉(zhuǎn)動(dòng)能etr,i、振動(dòng)-電子能eve,i、絕對焓hi等特性參數(shù)的計(jì)算方法參考文獻(xiàn)[1-2,12]。為簡化計(jì)算過程，單溫度模型通常采用CHEMKIN軟件[19-20]給出的擬合公式進(jìn)行計(jì)算。

為計(jì)算黏性系數(shù)，NNW-HyFLOW集成了文獻(xiàn)[11,20]提出的兩種擬合方法，其中Blotter方法作為默認(rèn)選項(xiàng)。Eucken關(guān)系式[2]用于計(jì)算組分熱傳導(dǎo)系數(shù)，混合氣體的總黏性系數(shù)和熱傳導(dǎo)系數(shù)均采用Wilke混合律計(jì)算。在假設(shè)Schmidt數(shù)為常數(shù)和考慮湍流對輸運(yùn)系數(shù)影響的基礎(chǔ)上，組分質(zhì)量擴(kuò)散系數(shù)采用文獻(xiàn)[2]給出的方法計(jì)算。

化學(xué)模型計(jì)算包含在源項(xiàng)的計(jì)算過程中，組分輸運(yùn)方程生成源項(xiàng)Θi的計(jì)算形式為

(2)

NNW-HyFLOW軟件目前可提供的化學(xué)模型包括Gupta模型[19]、Dunn-Kang模型[11,21]和Park模型[22]，應(yīng)用時(shí)可根據(jù)計(jì)算需要選取相應(yīng)的化學(xué)模型和組分個(gè)數(shù)構(gòu)建反應(yīng)式組合及其反應(yīng)速率系數(shù)。生成速率Rr項(xiàng)中的正、逆反應(yīng)速率系數(shù)kr,f和kr,b均采用Arrhenius經(jīng)驗(yàn)公式[18]進(jìn)行計(jì)算，Park模型的逆反應(yīng)速率系數(shù)已轉(zhuǎn)換為Arrhenius公式所需的參數(shù)形式，以方便程序的編寫以及保持計(jì)算形式的統(tǒng)一。

對于三溫度模型，振動(dòng)能量方程生成源項(xiàng)Θv和電子能量方程生成源項(xiàng)Θe的計(jì)算形式分別為

(3)

兩溫度模型方程通過振動(dòng)能量方程與電子能量方程迭加得到，因此兩溫度模型的能量方程生成源項(xiàng)由Θv和Θe相加得到。多振動(dòng)溫度模型源項(xiàng)及其相關(guān)參數(shù)模型方法參考文獻(xiàn)[12]。

2.3 數(shù)值方法

當(dāng)前，NNW-HyFLOW默認(rèn)條件下均采用LUSGS隱式時(shí)間推進(jìn)方法進(jìn)行流場數(shù)值迭代求解。LUSGS方法的統(tǒng)一形式[23]為

(4)

式中：Q為守恒量向量；R為右端項(xiàng)；D、L、U分別為通量項(xiàng)分裂得到的對角陣、下三角陣和上三角陣。熱化學(xué)非平衡計(jì)算條件下矩陣D包含源項(xiàng)S的雅可比矩陣Z。為此，NNW-HyFLOW軟件針對源項(xiàng)處理提供了全隱式和點(diǎn)隱式兩類計(jì)算方法。由于點(diǎn)隱式方法編程簡單，計(jì)算效率高，且穩(wěn)定性和收斂性也較好[5]，故默認(rèn)條件主要采用點(diǎn)隱式方法計(jì)算非平衡源項(xiàng)。

另外，底層數(shù)值算法庫可提供Roe、Steger和AUSM類等多種通量數(shù)值格式，以及Vanleer、Minmod、Vanalbada等常用限制器用于計(jì)算對流通量，黏性通量主要采用中心格式計(jì)算。

3 算例驗(yàn)證

完全氣體條件的算例測試與驗(yàn)證可參考文獻(xiàn)[17]。針對熱化學(xué)非平衡效應(yīng)，以下從激波位置、電子數(shù)密度、非平衡氣動(dòng)力和氣動(dòng)熱環(huán)境等非平衡流場特性的計(jì)算精準(zhǔn)度以及對復(fù)雜工程外形的適用性等方面對NNW-HyFLOW軟件已完成的求解器功能進(jìn)行校驗(yàn)和確認(rèn)。

3.1 HEG風(fēng)洞試驗(yàn)

針對HEG激波風(fēng)洞圓柱模型試驗(yàn)[24-25]開展數(shù)值模擬測試。計(jì)算采用單溫度、5組分Dunn-Kang化學(xué)模型，壁溫取Tw=300 K。計(jì)算網(wǎng)格如圖10所示，其對稱面網(wǎng)格量為129×91，第1層高度取Δh=1.0×10-6m。

圖10 HEG圓柱半模型及計(jì)算網(wǎng)格Fig.10 Half model of HEG cylinder and its computational grids

圖11給出了橫截面(z=0.01 m)的壓強(qiáng)等值線分布，同時(shí)圖12和圖13分別給出完全催化壁(FCW)和非催化壁(NCW)條件下計(jì)算得到的壁面壓力分布和熱流分布?？梢钥吹剑疚挠?jì)算得到的脫體激波位置與試驗(yàn)和文獻(xiàn)結(jié)果非常吻合，壁面壓力值分布與試驗(yàn)結(jié)果一致，同時(shí)兩種壁面條件計(jì)算得到的熱流分布趨勢符合實(shí)際物理機(jī)制和理論分析，證明了NNW-HyFLOW軟件具有可靠的非平衡流場模擬和激波捕獲能力，在高溫氣體非平衡氣動(dòng)力和氣動(dòng)熱預(yù)測方面具有較高的計(jì)算精度。

圖11 激波位置對比Fig.11 Comparison of shock wave positions

圖12 壁面壓力分布Fig.12 Pressure distribution on wall surface

圖13 壁面熱流分布Fig.13 Heat flux distribution on wall surface

3.2 RAM-C飛行試驗(yàn)

這里以RAM-C鈍錐模型[26-27]為模擬對象，考核NNW-HyFLOW軟件常用化學(xué)反應(yīng)模型的正確性，以及對等離子體分布預(yù)測的計(jì)算精度。測試飛行高度H=71 km，飛行速度U∞=7 650 m/s(Ma∞≈25.914 4)。計(jì)算采用單溫度模型以及7組分的Dunn-Kang、Gupta和Park等3種化學(xué)模型，來流氧氣(O2)和氮?dú)?N2)的質(zhì)量分?jǐn)?shù)分別為0.233和0.767，壁溫取Tw=1 500 K。計(jì)算網(wǎng)格如圖14所示，其對稱面網(wǎng)格量為161×64，第1層高度取Δh=1.0×10-5m。

圖14 RAM-C鈍錐模型及計(jì)算網(wǎng)格Fig.14 RAM-C blunt cone and its computational grids

圖15給出了軸向x/R=8.1處沿法向(y/R)的電子數(shù)密度分布?？梢钥吹?，3種化學(xué)模型分別計(jì)算得到的電子數(shù)密度曲線與AEROPH_Flow程序[12]的計(jì)算值非常吻合，與飛行試驗(yàn)測量值差異小，處于相同數(shù)量級。測試結(jié)果證明了NNW-HyFLOW軟件所含化學(xué)反應(yīng)模型的正確性和有效性，同時(shí)顯示了軟件在電子數(shù)密度分布預(yù)測方面具有較高的計(jì)算精度。

圖15 x/R=8.1處法向電子數(shù)密度分布Fig.15 Electron density distribution along normal direction at position of x/R=8.1

3.3 Electre鈍錐模型試驗(yàn)

這里選取Electre鈍錐模型[12,28]為研究對象，參考飛行試驗(yàn)和Muylaert等的數(shù)值分析結(jié)果[28]開展模擬測試，考核NNW-HyFLOW軟件對熱化學(xué)非平衡壁面熱流的模擬能力。計(jì)算采用兩溫度、7組分Gupta化學(xué)模型，來流O2和N2的質(zhì)量分?jǐn)?shù)分別為0.233和0.767，壁溫取Tw=343 K。計(jì)算網(wǎng)格采用類似RAM-C外形的網(wǎng)格劃分，半模型網(wǎng)格量約72萬，第1層網(wǎng)格高度取Δh=1.0×10-6m。

圖16和圖17分別給出完全催化壁(FCW)和非催化壁(NCW)條件下計(jì)算得到的駐點(diǎn)線組分質(zhì)量分?jǐn)?shù)分布和溫度分布，圖18給出了采用兩種壁面邊界條件計(jì)算得到的軸向壁面熱流分布。圖中：T為平動(dòng)溫度；Tv為振動(dòng)溫度?？梢钥吹?，本文計(jì)算的組分質(zhì)量分?jǐn)?shù)分布曲線與Muylaert數(shù)值結(jié)果吻合，溫度分布曲線與Muylaert數(shù)值結(jié)果一致。另外，兩種壁面條件計(jì)算的熱流值曲線與測量值相比分布趨勢合理，與AEROPH_Flow的計(jì)算值幾乎重合，與Muylaert的計(jì)算結(jié)果符合良好。測試結(jié)果證明了NNW-HyFLOW兩溫度模型求解器及化學(xué)反應(yīng)模型的正確性和可靠性，同時(shí)顯示了軟件在熱化學(xué)非平衡壁面熱流預(yù)測方面具有較高的計(jì)算精度。

圖17 駐點(diǎn)線溫度分布Fig.17 Temperature distribution on stagnation line

圖18 軸向壁面熱流分布Fig.18 Heat flux distribution on wall along axial direction

3.4 航天飛機(jī)STS飛行試驗(yàn)

選取“哥倫比亞”號(OV-102)航天飛機(jī)飛行試驗(yàn)條件[12,29]進(jìn)一步開展數(shù)值模擬測試。測試飛行高度H=57.2 km，飛行速度U∞=4 560 m/s(Ma∞≈14.25)，攻角α=40.56°。計(jì)算采用單溫度、5組分Park化學(xué)模型，來流O2和N2的質(zhì)量分?jǐn)?shù)分別為0.233和0.767，壁溫Tw=1 500 K。計(jì)算網(wǎng)格如圖19所示，半模型網(wǎng)格量約652萬，第1層網(wǎng)格高度取Δh=1.0×10-6m。

圖19 航天飛機(jī)OV-102模型及計(jì)算網(wǎng)格Fig.19 Space shuttle OV-102 model and its computational grids

圖20和圖21給出對稱面O2、NO的質(zhì)量分?jǐn)?shù)等值線分布，以及壁面的壓強(qiáng)、熱流分布云圖，同時(shí)圖22給出了采用非催化壁(NCW)條件計(jì)算得到的機(jī)身迎風(fēng)面中心線的熱流分布?？梢钥吹?，本文計(jì)算得到的壁面熱流值與STS-2、STS-5等飛行試驗(yàn)測量值[29]符合較好，與AEROPH_Flow計(jì)算值基本吻合，顯示NNW-HyFLOW軟件在工程型號復(fù)雜外形化學(xué)非平衡流動(dòng)及其影響的氣動(dòng)熱特性模擬方面具有較高的計(jì)算可信度，能夠適應(yīng)高超聲速復(fù)雜飛行器高溫非平衡流動(dòng)數(shù)值模擬的需求。

圖20 對稱面主要組分分布云圖Fig.20 Distribution contour of main species on symmetry plane

圖21 壁面主要?dú)鈩?dòng)參數(shù)分布云圖Fig.21 Distribution contour of main aerodynamic parameters on surface

圖22 機(jī)身迎風(fēng)面中心線熱流分布Fig.22 Heat flux distribution on centerline of windward side of fuselage

4 結(jié)束語

基于風(fēng)雷開源軟件的基礎(chǔ)框架和面向?qū)ο蟮木幊趟枷?，開發(fā)的NNW-HyFLOW高超聲速流動(dòng)模擬軟件具有底層代碼復(fù)用、功能兼容性好、擴(kuò)展能力強(qiáng)以及接口靈活等優(yōu)點(diǎn)。當(dāng)前測試版本初步實(shí)現(xiàn)完全氣體、化學(xué)非平衡氣體和熱化學(xué)非平衡氣體數(shù)據(jù)結(jié)構(gòu)的一體化設(shè)計(jì)，集成了常用空氣化學(xué)反應(yīng)5組分、7組分和11組分Park、Gupta和Dunn-Kang模型，主要構(gòu)建了熱力學(xué)一溫度、兩溫度和三溫度能量輸運(yùn)體系，初步具備了多種氣體模態(tài)熱化學(xué)非平衡效應(yīng)數(shù)值模擬能力。結(jié)合HEG風(fēng)洞圓柱模型試驗(yàn)、RAM-C飛行試驗(yàn)、Electre鈍錐模型試驗(yàn)和航天飛機(jī)STS飛行試驗(yàn)條件，開展了高溫氣體非平衡流數(shù)值模擬校驗(yàn)。軟件測試表明：該軟件測試版本在熱化學(xué)非平衡效應(yīng)及其影響的氣動(dòng)力特性、氣動(dòng)熱環(huán)境和等離子體分布特性預(yù)測與評估方面，具有較高的數(shù)值計(jì)算精準(zhǔn)度，初步滿足了高超聲速復(fù)雜飛行器高溫非平衡流動(dòng)數(shù)值模擬的需求。

為了進(jìn)一步提升NNW-HyFLOW軟件的綜合模擬能力，后續(xù)將在以下方面開展優(yōu)化：

1) 在網(wǎng)格架構(gòu)多樣性方面，基于現(xiàn)有數(shù)據(jù)結(jié)構(gòu)和功能模塊，實(shí)現(xiàn)結(jié)構(gòu)/非結(jié)構(gòu)混合和重疊網(wǎng)格擴(kuò)展，更好地滿足復(fù)雜外形飛行器局部細(xì)微結(jié)構(gòu)精細(xì)模擬能力。

2) 在功能的完整性方面，繼續(xù)集成國內(nèi)各類物理化學(xué)模型，包括但不限于新的地球大氣熱化學(xué)反應(yīng)體系、火星大氣反應(yīng)體系、多種燃?xì)夥磻?yīng)體系以及燒蝕反應(yīng)體系；進(jìn)一步完善多物理效應(yīng)耦合模擬模塊，包括但不限于湍流效應(yīng)、稀薄滑移、氣體輻射、表面有限催化及電磁效應(yīng)的耦合策略和接口設(shè)計(jì)。

3) 在代碼的可擴(kuò)展性和模塊的可移植性方面，繼續(xù)加強(qiáng)軟件各功能模塊的封裝，加強(qiáng)各級數(shù)據(jù)之間的分層管理，進(jìn)一步提升軟件移植性。

4) 在程序執(zhí)行的效率和穩(wěn)定性方面，持續(xù)優(yōu)化并行底層封裝和并行原理構(gòu)建，結(jié)合多種隱式耦合算法，提升軟件的計(jì)算效率和數(shù)值格式的魯棒性。

致 謝

本研究得到了中國空氣動(dòng)力研究與發(fā)展中心計(jì)算空氣動(dòng)力研究所風(fēng)雷課題組和氣動(dòng)物理課題組的大力支持，在此表示感謝。