超大規(guī)模氣動彈性數(shù)值模擬軟件研制(2017)

王運濤, 孟德虹, 孫 巖, 洪俊武

(1. 中國空氣動力研究與發(fā)展中心 計算空氣動力研究所, 四川 綿陽 621000;2. 中國空氣動力研究與發(fā)展中心 空氣動力學(xué)國家重點實驗室, 四川 綿陽 621000)

0 引 言

高性能計算的峰值性能已于2008年跨越P級(PetaFlops，每秒千萬億次浮點運算)門檻，目前正在向E級(ExaFlops，每秒百億億次浮點運算)進發(fā)。自2010年11月以來，以“太湖之光”(峰值性能125.4P，持續(xù)性能93.0P)、“天河二號”(峰值性能54.9P，持續(xù)性能33.9P)為代表的超級計算機, 已經(jīng)多次在全球TOP 500超級計算機排行榜名列第一。截止2017年11月，中國超級計算機TOP 500的上榜總數(shù)達到歷史最高的202臺。隨著國內(nèi)超級計算機硬件技術(shù)的飛速發(fā)展，國內(nèi)超大規(guī)模的計算應(yīng)用也取得了突破性進展。2016年，中國科學(xué)院軟件研究所楊超等人[1]在美國鹽湖城舉行的全球超級計算大會上獲得了國際高性能計算應(yīng)用領(lǐng)域最高獎“戈登貝爾”獎，獲獎應(yīng)用名稱為“千萬核可擴展全球大氣動力學(xué)全隱式模擬”。2017年，清華大學(xué)付昊桓等人在美國丹佛舉行的全球超級計算大會上再次獲得“戈登貝爾”獎，獲獎應(yīng)用項目名稱為“非線性地震模擬”。隨著我國航空航天飛行器的自主創(chuàng)新發(fā)展，國內(nèi)飛行器設(shè)計與研制單位正在逐步裝備百萬億次量級的高性能計算機，其中中國空氣動力研究與發(fā)展中心裝備的高性能計算機已經(jīng)達到1.5P。以高保真氣動數(shù)據(jù)庫的高效生產(chǎn)和關(guān)鍵氣動技術(shù)攻關(guān)為牽引的計算空氣動力學(xué)(CFD)正在成為國產(chǎn)高性能計算機持續(xù)發(fā)展的原始驅(qū)動力之一。

高性能計算機的發(fā)展是飛行器多學(xué)科耦合設(shè)計的硬件基礎(chǔ)，而多學(xué)科耦合的數(shù)值模擬工具又是多學(xué)科耦合設(shè)計的軟件基礎(chǔ)，氣動彈性問題的數(shù)值模擬屬于典型的多學(xué)科耦合數(shù)值模擬問題。始于2010年、由AIAA組織的DPW(Drag Prediction Workshop)系列會議已經(jīng)成功舉辦了六屆，其主要目的是評估現(xiàn)代CFD技術(shù)模擬運輸機高速構(gòu)型的阻力預(yù)測能力，為CFD技術(shù)下一步的發(fā)展提出建議和意見。為考核現(xiàn)有流固耦合軟件的數(shù)值模擬能力，2016年6月召開的第六屆DPW組委會將CRM(Common Research Model)翼身組合體模型的流固耦合數(shù)值模擬(Case 5)作為可選工況。來自世界各地的25個研究機構(gòu)采用25種軟件提供了48組網(wǎng)格收斂性(Case 1)的數(shù)值模擬結(jié)果[2]，其中只有4家研究機構(gòu)提供了4組流固耦合的計算結(jié)果。由此可見相對于常規(guī)的CFD數(shù)值模擬，靜氣動彈性數(shù)值模擬技術(shù)尚未成熟[3]。始于2012年的由AIAA組織的AePW(Aeroelastic Prediction Workshop)系列會議已經(jīng)舉辦了兩屆，其主要目的是評估氣動彈性數(shù)值模擬技術(shù)現(xiàn)狀，為氣動彈性數(shù)值模擬技術(shù)和相應(yīng)風(fēng)洞試驗技術(shù)的發(fā)展提出意見和建議。2016年1月召開的第二屆AePW組委會將BSCW(Benchmark Supercritical Wing)模型的顫振模擬(Case2)作為必選工況[4]。來自世界各地的7個研究團隊提供了13組流固耦合的數(shù)值模擬結(jié)果，數(shù)值模擬得到的顫振起始動壓與試驗結(jié)果最大相差22%。這說明，顫振數(shù)值模擬技術(shù)與風(fēng)洞試驗技術(shù)本身，均需要進一步開展研究工作。

由中國空氣動力研究與發(fā)展中心(CARDC)與大連理工大學(xué)、北京航空航天大學(xué)、中國商飛上海飛機設(shè)計研究院、成都飛機設(shè)計研究所等五家單位聯(lián)合申請的國家重點研發(fā)計劃項目“數(shù)值飛行器原型系統(tǒng)開發(fā)”于2016年8月正式立項。其中，研制具有完全自主知識產(chǎn)權(quán)的超大規(guī)模氣動彈性數(shù)值模擬軟件是項目的重要研究內(nèi)容之一。該軟件的研究目標(biāo)是：基于本項目開發(fā)的飛行器空氣動力學(xué)模擬軟件和結(jié)構(gòu)動力學(xué)分析軟件，突破百億量級超大規(guī)模流場計算網(wǎng)格高效魯棒變形、海量耦合數(shù)據(jù)高效高精度傳遞等技術(shù)瓶頸，實現(xiàn)百億量級網(wǎng)格規(guī)模、60萬處理器核的超大規(guī)模靜/動氣動彈性數(shù)值模擬，構(gòu)建飛行器氣動彈性機理研究和工程應(yīng)用的軟件平臺，提高飛行器氣動彈性問題數(shù)值模擬的效率和預(yù)測精度，滿足型號研制對氣動彈性模擬能力的需求。

本文介紹了超大規(guī)模氣動彈性數(shù)值模擬軟件的研制情況，主要包括超大規(guī)模結(jié)構(gòu)網(wǎng)格生成技術(shù)、基于結(jié)構(gòu)網(wǎng)格的動網(wǎng)格技術(shù)、氣動彈性數(shù)值模擬軟件主要功能模塊的并行效率和并行一致性測試、軟件的確認工作等內(nèi)容。

1 超大規(guī)模氣動彈性數(shù)值模擬軟件簡介

超大規(guī)模氣動彈性數(shù)值模擬軟件主要包括CFD求解模塊、CSM(Computational Structural Mechanics)求解模塊、耦合界面數(shù)據(jù)傳遞和動態(tài)網(wǎng)格變形4個主要功能模塊，通過主控程序有序組織上述功能模塊，從而實現(xiàn)復(fù)雜飛行器氣動彈性問題的數(shù)值模擬。其中，超大規(guī)模并行計算技術(shù)貫穿于軟件的各個功能模塊之中。目前，“數(shù)值飛行器原型系統(tǒng)開發(fā)”項目組已經(jīng)開發(fā)了超大規(guī)模氣動彈性數(shù)值模擬軟件的初級版本，各個功能模塊的主要算法簡單介紹如下。

(1) CFD求解模塊。采用CARDC自主研發(fā)的亞跨超CFD軟件平臺(Trisonic Platform，TRIP)——TRIP軟件,經(jīng)過了系統(tǒng)的驗證和確認工作[5-6]并廣泛應(yīng)用于多種飛行器的氣動設(shè)計與評估。該軟件采用有限體積方法離散任意坐標(biāo)下的雷諾平均Navier-Stokes方程(RANS)，RANS方程無黏項離散采用二階精度ROE迎風(fēng)格式[7]，黏性項離散采用二階中心格式，湍流模型包含Spalart-Allmaras SA一方程模型[8]、Menter’s SST兩方程模型[9]，定常問題離散方程組的求解采用LU-SGS方法[10]，非定常問題離散方程組的求解采用Jameson“雙時間步”方法[11]，采用低速預(yù)處理技術(shù)、多重網(wǎng)格技術(shù)和基于MPI或MPI/OpenMP的大規(guī)模并行技術(shù)加速收斂。

(2) CSM計算模塊。靜氣動彈性數(shù)值模擬采用柔度矩陣方法求解線性結(jié)構(gòu)靜力學(xué)方程，獲得氣動載荷作用下的物面變形；顫振等動氣動彈性數(shù)值模擬采用模態(tài)法求解結(jié)構(gòu)動力學(xué)方程。目前，結(jié)構(gòu)柔度矩陣以及模態(tài)的計算均基于飛行器有限元模型采用商業(yè)結(jié)構(gòu)動力學(xué)分析軟件獲得，適用于飛行器結(jié)構(gòu)線性變化范圍內(nèi)氣動彈性問題的模擬，對于結(jié)構(gòu)非線性氣動彈性問題的模擬需要進一步開展研究工作。

(3) 耦合界面數(shù)據(jù)傳遞模塊。采用無限板樣條IPS(Infinite Plate Spline)方法或薄板樣條TPS(Thin Plate Spline)插值方法[12]構(gòu)建CFD計算模塊與CSM計算模塊之間的氣動載荷與結(jié)構(gòu)變形傳遞矩陣。

(4) 動態(tài)網(wǎng)格變形模塊。CFD計算網(wǎng)格的變形采用基于徑向基函數(shù)RBF(Radial Basis Functions)[13]與超限插值TFI(Transfinite Interpolation)[14]的復(fù)合型動態(tài)網(wǎng)格變形方法RBF_TFI[15]。

(5) 靜氣動彈性計算流程。采用松耦合方式建立靜氣動彈性模擬的一般流程(圖1)。第一步，CFD模塊計算出第n迭代步流場Un后，將壁面壓力Pn傳遞給CSM模塊；第二步，CSM模塊利用輸入的氣動載荷和邊界約束條件，計算出第n+1迭代步的結(jié)構(gòu)變形Wn+1；第三步，將結(jié)構(gòu)變形Wn+1傳遞給CFD模塊；第四步，CFD模塊利用變形后外形，更新計算網(wǎng)格，計算第n+1迭代步的流場Un+1；然后循環(huán)上述過程，不斷得到下一迭代步的位移和流場，直至變形位移和流場均達到收斂。

圖1 靜氣動彈性問題模擬流程Fig.1 Outline of static aero-elastic problem simulation

(6) 動氣動彈性計算流程。采用松耦合方式建立動氣動彈性模擬的一般流程(圖2)。第一步，CFD模塊計算出t0時間步流場Ut0后，將壁面壓力Pt0傳遞給CSM模塊；第二步，CSM模塊利用輸入的氣動載荷和邊界約束條件，基于模態(tài)迭加的結(jié)構(gòu)動力學(xué)方程，采用四步Runge-Kutta方法計算出t0+Δt時間步的結(jié)構(gòu)變形Wt0+Δt；第三步，將結(jié)構(gòu)變形Wt0+Δt傳遞給CFD模塊；第四步，CFD模塊利用變形后外形，更新流體計算網(wǎng)格，采用“雙時間步”方法計算t0+Δt時間步的流場Ut0+Δt。循環(huán)上述過程，不斷得到下一時刻的位移和流場，直至計算過程滿足收斂判據(jù)。由于CFD求解模塊和CSM求解模塊在數(shù)據(jù)交換時存在一個時間步的延遲，松耦合方式在時間方向只具有一階時間精度。下一步將采用預(yù)估-校正技術(shù)或?qū)SM求解模塊同樣構(gòu)造成含有子迭代的求解方式以提高動氣動彈性問題數(shù)值模擬的時間精度。

圖2 動氣動彈性計算流程Fig.2 Outline of dynamic aero-elastic problem simulation

2 超大規(guī)模氣動彈性數(shù)值模擬軟件研制

2.1 超大規(guī)模靜態(tài)結(jié)構(gòu)網(wǎng)格生成技術(shù)

超大規(guī)模靜態(tài)網(wǎng)格生成技術(shù)是實現(xiàn)復(fù)雜構(gòu)型超大規(guī)模氣動彈性問題數(shù)值模擬的前提。以包含機翼/機身/掛架/短艙/平尾/立尾的典型運輸機巡航構(gòu)型為例，目前模擬上述復(fù)雜構(gòu)型的結(jié)構(gòu)網(wǎng)格規(guī)模一般在千萬量級，在配置較高的微機或工作站上利用商業(yè)網(wǎng)格生成軟件即可完成。而構(gòu)造數(shù)億量級乃至百億量級的結(jié)構(gòu)網(wǎng)格已經(jīng)遠遠超出了當(dāng)前單機能力范圍，必須采用新的網(wǎng)格生成策略。超大規(guī)模靜態(tài)網(wǎng)格生成的關(guān)鍵技術(shù)主要包括物面網(wǎng)格的高保真度技術(shù)，高質(zhì)量、魯棒的空間網(wǎng)格構(gòu)造技術(shù)等方面。

超大規(guī)模結(jié)構(gòu)網(wǎng)格的生成步驟是：1) 采用商業(yè)軟件構(gòu)造千萬量級高質(zhì)量初始結(jié)構(gòu)網(wǎng)格；2) 利用初始結(jié)構(gòu)網(wǎng)格拓撲分塊信息，建立網(wǎng)格窗口的信息數(shù)據(jù)庫；3)基于網(wǎng)格窗口的信息數(shù)據(jù)庫創(chuàng)建不同網(wǎng)格窗口之間的關(guān)聯(lián)關(guān)系；4) 根據(jù)網(wǎng)格規(guī)模需求，對棱線網(wǎng)格點進行加密，重新計算棱線網(wǎng)格點分布及棱線網(wǎng)格點坐標(biāo)，通過超限插值(TFI)生成加密后的窗口網(wǎng)格；5) 采用線性投影方法，對邊界屬性標(biāo)記為物面的窗口網(wǎng)格重新向物面投影，獲得加密后的物面網(wǎng)格；6) 利用三維TFI生成加密后的體網(wǎng)格。上述超大規(guī)模網(wǎng)格生成策略與待模擬外形的復(fù)雜程度無關(guān)，所構(gòu)造的超大規(guī)模結(jié)構(gòu)網(wǎng)格拓撲結(jié)構(gòu)和網(wǎng)格單元質(zhì)量與初始網(wǎng)格完全相同，并可以對初始網(wǎng)格進行不同程度的粗化或細化。采用上述方法，可以快速獲得百億量級以上的高質(zhì)量結(jié)構(gòu)網(wǎng)格，同時可以構(gòu)造一系列的不同規(guī)模網(wǎng)格開展網(wǎng)格收斂性研究。

按照前述方法，對CRM-WBH構(gòu)型進行了驗證性網(wǎng)格重構(gòu)。初始網(wǎng)格規(guī)模為9 237萬網(wǎng)格單元，采用網(wǎng)格粗化和細化的方法獲得了105 219萬(10億量級)、591 187萬(60億量級)和1 154 662萬(115億量級)網(wǎng)格單元的不同規(guī)模網(wǎng)格。圖1給出了115億網(wǎng)格單元的CRM-WBH構(gòu)型空間網(wǎng)格拓撲。表1給出了不同規(guī)模網(wǎng)格的主要信息，其中Ratio表示棱線網(wǎng)格單元加密倍數(shù)，Nnode表示網(wǎng)格結(jié)點總數(shù)，NBL表示邊界層網(wǎng)格的法向網(wǎng)格點數(shù)目，y+表征第一層網(wǎng)格法向無量綱距離，Nblock表示總的網(wǎng)格塊數(shù)目，Size表示網(wǎng)格文件所占存儲量。

圖3 CRM-WBH構(gòu)型的網(wǎng)格拓撲(115億)Fig.3 Grid topology of CRM-WBH configuration (11.5 billion)

ParameterRatio=1.0Ratio=2.25Ratio=4.0Ratio=5.0Nnode92 372 9921 052 186 1125 911 871 48811 546 624 000NBL97217385481y+0.50.220.1250.1Nblock73737373Size/MB2 27725 556142 980278 829

2.2 超大規(guī)模結(jié)構(gòu)網(wǎng)格變形技術(shù)

高質(zhì)量的網(wǎng)格變形技術(shù)是CFD領(lǐng)域的研究熱點，也是氣動彈性數(shù)值模擬的關(guān)鍵技術(shù)之一。高質(zhì)量、魯棒高效的網(wǎng)格變形技術(shù)應(yīng)滿足以下要求：1)自動化。減少甚至避免網(wǎng)格變形過程中的人工干預(yù)是氣動彈性問題高效數(shù)值模擬的基本要求。2)并行化。網(wǎng)格變形過程的并行化是提高網(wǎng)格變形效率有效途徑，對于超大規(guī)模網(wǎng)格的動氣動彈性問題數(shù)值模擬，網(wǎng)格變形并行化的需求尤其迫切。3) 適用性。網(wǎng)格變形方法要能夠處理各種復(fù)雜外形、大變形，以及各種復(fù)雜變形(包括不規(guī)則變形)。4)通用性。網(wǎng)格變形方法不僅適用于結(jié)構(gòu)網(wǎng)格，還適用于各種類型的非結(jié)構(gòu)網(wǎng)格及混合網(wǎng)格。

本文選擇了基于徑向基函數(shù)的超限插值方法(RBF_TFI)實現(xiàn)復(fù)雜構(gòu)型結(jié)構(gòu)網(wǎng)格的自動變形，利用刪減基點后的徑向基函數(shù)插值更新結(jié)構(gòu)網(wǎng)格塊的棱線坐標(biāo)，然后利用TFI插值重新生成計算網(wǎng)格。該方法的具體實施步驟是：1)根據(jù)輸入的多塊結(jié)構(gòu)網(wǎng)格信息，在物面網(wǎng)格處選擇徑向基函數(shù)基點，構(gòu)造可逆的插值矩陣；2)根據(jù)輸入的物面變形確定物面徑向基基點位移，利用可逆的插值矩陣確定徑向基插值系數(shù)；3)采用徑向基插值方法更新物面與空間網(wǎng)格的網(wǎng)格棱線坐標(biāo)；4)采用二維、三維TFI方法構(gòu)造各個網(wǎng)格塊的面網(wǎng)格與體網(wǎng)格。RBF_TFI方法憑借優(yōu)異的性能在很多程序和算例中得到了廣泛應(yīng)用，但在一些比較極端的情況下，邊界層網(wǎng)格第一層高度十分小的時候(這在高雷諾數(shù)流動問題和超大規(guī)模網(wǎng)格問題中經(jīng)常遇到)，RBF_TFI容易引起邊界層區(qū)域網(wǎng)格棱線的扭曲現(xiàn)象，并造成物面附近網(wǎng)格的交叉破壞問題，從而導(dǎo)致流場計算無法繼續(xù)進行?；诓煌憷龑υ摤F(xiàn)象進行了深入的分析，認為是物面網(wǎng)格點與空間網(wǎng)格點采用不同的坐標(biāo)更新方法導(dǎo)致的變形不一致引起了網(wǎng)格交錯破壞。通過現(xiàn)象、原因等方面的詳細分析，提出一種解決該問題的改進方法，改進后的網(wǎng)格更新流程見圖4。

圖4 RBF_TFI結(jié)構(gòu)網(wǎng)格變形工作流程Fig.4 Workflow of structured grid deformation with RBF-TFI method

在顫振問題的數(shù)值模擬中，為了保證動氣動彈性問題數(shù)值模擬時間上的二階精度，流體力學(xué)方程求解和結(jié)構(gòu)力學(xué)方程求解一般采用真實時間步加上子迭代的方式，在每一個子迭代步均需更新流場計算的網(wǎng)格。在這種情形下，CFD計算網(wǎng)格的更新效率至關(guān)重要，尤其是對于億量級或數(shù)十億量級網(wǎng)格的流場計算。本文采用基于MPI的并行處理技術(shù)，開展了結(jié)構(gòu)網(wǎng)格變形并行化技術(shù)研究，初步建立了結(jié)構(gòu)網(wǎng)格并行變形的一般流程。具體實現(xiàn)方法如下：首先主進程收集所有進程物面變形，計算全局插值系數(shù)，然后將插值系數(shù)播撒到所有子進程，每個核利用插值系數(shù)矩陣獨立采用RBF或RBF+TFI方法更新體網(wǎng)格。物面網(wǎng)格的更新方法與體網(wǎng)格的更新方法相同。

2.3 超大規(guī)模氣動/結(jié)構(gòu)數(shù)據(jù)并行傳遞技術(shù)

氣動彈性問題數(shù)值模擬時，需要將CSM模塊獲得氣動載荷作用下的結(jié)構(gòu)變形傳遞給動態(tài)網(wǎng)格處理模塊，進而開展下一個迭代步的CFD計算；將CFD模塊獲得的氣動載荷傳遞給CSM模塊以開展結(jié)構(gòu)分析，位移數(shù)據(jù)的傳遞與氣動載荷數(shù)據(jù)的傳遞均在物面進行。這涉及兩個方面的問題，一個是數(shù)據(jù)傳遞的精度與守恒性問題，另一個是針對超大規(guī)模氣動特性問題模擬的效率問題。第一個問題已經(jīng)在前期工作中得以解決[15-16]。針對第二個問題，為了使得數(shù)據(jù)傳遞的效率與氣動特性問題數(shù)值模擬的效率相匹配，采用了基于MPI的并行處理技術(shù)，并行化的耦合數(shù)據(jù)傳遞基于徑向基函數(shù)插值技術(shù)實現(xiàn)，與動態(tài)變形網(wǎng)格并行化的思路十分相似。具體實現(xiàn)方法是：首先在主結(jié)點實現(xiàn)插值系數(shù)方程組的求解，獲得全局耦合插值的系數(shù)，然后將插值系數(shù)播撒到各個CPU核，每個核獨立實現(xiàn)耦合變量的傳遞，最后組合實現(xiàn)全局耦合界面數(shù)據(jù)的高效傳遞。

3 氣動彈性數(shù)值模擬軟件的并行測試

2017年度，采用AGARD445.6機翼對超大規(guī)模氣動彈性數(shù)值模擬的動網(wǎng)格模塊開展了并行效率測試；采用CRM-WBH模型，對CFD求解模塊開展了并行一致性測試及并行效率測試；采用AGARD445.6機翼，對動氣動彈性模擬模塊開展了并行效率測試。

3.1 CRM-WBH模型與AGARD445.6顫振模型

CRM-WBH模型是2010年召開的第四屆AIAA阻力預(yù)測研討會(DPW IV)的基準研究模型[17]，設(shè)計馬赫數(shù)為0.85，設(shè)計升力系數(shù)為0.50。該模型先后在NASA Langley的NTF風(fēng)洞和NASA Ames的TWT風(fēng)洞中開展了多種構(gòu)型的風(fēng)洞試驗。風(fēng)洞試驗?zāi)Ｐ筒捎冒惭b于機身后體的葉片尾撐方式固定于風(fēng)洞迎角變換裝置，試驗結(jié)果包括了氣動特性、表面壓力分布等。DPW IV會議上數(shù)值模擬得到的氣動特性與相應(yīng)風(fēng)洞試驗結(jié)果的對比并不理想[18-19]，主要表現(xiàn)：相同迎角下，數(shù)值模擬得到的升力系數(shù)普遍大于試驗結(jié)果；相同升力系數(shù)下，數(shù)值模擬得到的低頭力矩系數(shù)普遍大于試驗結(jié)果。

AGARD445.6機翼是用于考核跨聲速顫振計算方法和軟件的標(biāo)準模型[20]，風(fēng)洞試驗是在美國NASA Langley 研究中心的跨聲速動態(tài)風(fēng)洞完成的。AGARD445.6機翼弦向翼剖面為NASA65A004翼型，試驗?zāi)Ｐ鸵圆馁|(zhì)均勻的桃花心木薄板制成。機翼根弦長0.5587 m、展長為0.762 m，機翼展弦比1.6525，根梢比1.5207，1/4弦長后掠角為45°。

3.2 動網(wǎng)格模塊并行效率測試

表2給出了AGARD445.6機翼2.05億結(jié)構(gòu)網(wǎng)格動網(wǎng)格模塊并行效率測試結(jié)果。采用128～4096個處理器測試了動網(wǎng)格模塊的并行效率，網(wǎng)格更新時間相對于流場單步子迭代的計算時間控制在14%以內(nèi)。在廣州超算中心的天河二號集群上開展了十億量級網(wǎng)格變形并行效率測試。圖5給出了AGARD445.6機翼模型16億多塊結(jié)構(gòu)網(wǎng)格變形并行測試的結(jié)果，圖中橫坐標(biāo)為采用的處理器核數(shù)，縱坐標(biāo)為網(wǎng)格更新時間與單步CFD流場計算時間的比值。對于16億的流場計算網(wǎng)格，從2 000核到25 000核，網(wǎng)格變形更新時間占單次CFD計算時間的比例基本穩(wěn)定在15%左右。上述測試表明，采用的動網(wǎng)格并行策略是可行的，并行效率不僅滿足工程問題的效率需求，也為數(shù)十億網(wǎng)格規(guī)模的氣動彈性問題數(shù)值模擬提供了一個可能的解決方案。

表2 2.05億結(jié)構(gòu)網(wǎng)格變形并行效率測試Table 2 Parallel efficiency test of 20.5 million structured grid deformation

圖5 16億結(jié)構(gòu)網(wǎng)格變形并行效率測試Fig.5 Parallel efficiency test of 1.6 billion structured grid deformation

3.3 靜氣動彈性模塊并行效率及一致性測試

采用CRM-WBH構(gòu)型、1.8億網(wǎng)格，開展了CFD模塊的并行效率測試。表3給出了采用單重網(wǎng)格、960個處理器核到9600個處理器核并行計算相對于96個處理器核的加速比和并行計算效率。在同構(gòu)MPI并行模式下，9600個處理器核的并行效率達到70%左右。并行一致性是大規(guī)模并行計算的另一個關(guān)鍵問題，采用CRM-WBH構(gòu)型開展了CFD軟件的并行一致性測試。計算來流條件為：Ma=0.85，Re=5.0×106，α=3.0°。表4給出了96～9600處理器核并行計算得到的CRM-WBH構(gòu)型的氣動特性。當(dāng)處理器規(guī)模由96個增加到9600個時，升力系數(shù)的變化量為0.12%、阻力系數(shù)的變化量為0.34%、俯仰力矩系數(shù)的變化量為1.1%。計算結(jié)果一致性良好。

表3 1.8億結(jié)構(gòu)網(wǎng)格CFD并行效率測試Table 3 Parallel efficiency test of CFD with 18.0 million structured grid

表4 1.8億結(jié)構(gòu)網(wǎng)格CRM-WBH構(gòu)型并行一致性測試Table 4 Parallel consistency test of CRM-WBH configuration with 18.0 million structured grid

3.4 動氣動彈性模塊并行效率測試

采用AGARD445.6機翼16億規(guī)模的多塊對接結(jié)構(gòu)網(wǎng)格，在廣州超算中心天河二號超級集群上進行了動氣動彈性模塊的大規(guī)模并行效率測試，并行規(guī)模為1536～24 576個處理器核。圖6給出了并行效率和加速比，由圖可以看出，采用MPI并行方式，相對于1536個處理器核，并行規(guī)模為24 576個處理器核時，并行效率達到了77.63%。

圖6 AGARD445.6機翼顫振計算并行效率Fig.6 Parallel efficiency of flutter simulation for AGARD 445.6 wing

4 氣動彈性數(shù)值模擬軟件的確認

在前期工作中，開展了HIRENASD機翼模型、DLR-F6模型[16]，CRM-WB模型[21]靜氣動彈性模擬測試，初步確認了氣動彈性模擬軟件的靜氣動彈性模擬精度。本文采用CRM-WBH構(gòu)型進一步開展了靜氣動彈性數(shù)值模擬精度確認，采用AGARD445.6機翼開展了動氣動彈性數(shù)值模擬精度確認。

4.1 CRM-WBH構(gòu)型靜氣動彈性數(shù)值模擬

采用CRM-WBH模型、支撐裝置數(shù)模及風(fēng)洞試驗?zāi)Ｐ偷挠邢拊Ｐ烷_展了CRM-WBH構(gòu)型的靜氣動彈性模擬。以下為敘述方便，將包含尾撐裝置的CRM翼/身/平尾組合體構(gòu)型簡稱為CRM-WBHS。圖7給出了包含支撐裝置的CRM-WBHS構(gòu)型多塊對接結(jié)構(gòu)網(wǎng)格拓撲及表面網(wǎng)格，網(wǎng)格規(guī)模6061萬。圖8給出了CRM-WBH風(fēng)洞試驗?zāi)Ｐ偷挠邢拊Ｐ汀?/p>

圖7 CRM-WBHS構(gòu)型網(wǎng)格拓撲及表面網(wǎng)格Fig.7 Grid topology and surface grid for CRM-WBHS configuration

圖8 CRM-WBH風(fēng)洞試驗有限元模型Fig.8 Finite element model CRM-WBH wind tunnel test

圖9為不同迎角下靜彈性計算最大結(jié)構(gòu)變形隨迭代步數(shù)的變化曲線，不同迎角下的計算結(jié)果都可以看出,在迭代到第5步后，最大變形量基本不變。

圖10給出了CRM-WBH構(gòu)型氣動特性隨迎角的變化，計算來流狀態(tài)：Ma=0.85，Re=5.0×106，α=0°～4.00°；流固耦合計算時，速壓q=61.062kPa，載荷因子q/E=3.342×10-7,E為彈性模量。其中，CRM-WBHS_CFD代表采用CFD方法得到的氣動特性數(shù)值模擬結(jié)果，CRM-WBHS_FSC代表采用流固耦合方法得到的氣動特性數(shù)值模擬結(jié)果，同時圖中給出了ETW風(fēng)洞試驗的測力結(jié)果。由圖可見，采用流固耦合方法，考慮模型的靜氣動彈性變形影響導(dǎo)致相同迎角下升力系數(shù)、阻力系數(shù)下降，俯仰力矩系數(shù)增加。CRM-WBHS_FSC的模擬結(jié)果更接近試驗結(jié)果。氣動特性數(shù)值模擬結(jié)果與試驗結(jié)果的差異需要從洞壁干擾、轉(zhuǎn)捩模擬等方面進一步開展工作。

圖9 CRM 構(gòu)型機翼最大變形量收斂歷程Fig.9 Convergence process of wing maximal deformation for CRM configuration

(a) CL～α

(b) CD～α

(c) Cm～α

4.2 AGARD445.6機翼顫振數(shù)值模擬

采用松耦合方式開展AGARD445.6機翼的顫振特性數(shù)值模擬。CFD計算網(wǎng)格規(guī)模為2566萬，CFD網(wǎng)格的變形采用并行處理方式；結(jié)構(gòu)動力學(xué)方程求解取AGARD445.6機翼前四階模態(tài)；CFD計算網(wǎng)格、各階模態(tài)的固有頻率如圖11所示。圖12給出了AGARD445.6機翼顫振計算得到的廣義位移響應(yīng)曲線，來流狀態(tài)為Ma=0.96，q/qe=0.8650時，計算得到的顫振速度系數(shù)為Vf=0.2817，其中qe為試驗動壓。

圖11 AGARD445.6機翼CFD計算網(wǎng)格與結(jié)構(gòu)模態(tài)Fig.11 CFD grid and structural modal for AGARD445.6 wing

圖12 AGARD445.6機翼廣義位移隨時間變化歷程Fig.12 Time history of generalized displacement for AGARD445.6 wing

圖13為AGARD445.6機翼顫振計算得到的各馬赫數(shù)臨界顫振速度系數(shù)與試驗的比較，馬赫數(shù)范圍0.499～1.072，同時給出了320萬網(wǎng)格規(guī)模的數(shù)值模擬結(jié)果。采用320萬網(wǎng)格規(guī)模進行顫振模擬，采用720個FT1500處理器核，迭代5000個時間步需要65.9 h；采用2566萬網(wǎng)格規(guī)模進行顫振模擬，采用1600個FT1500處理器核，迭代5000個時間步需要255.6 h。由圖13可見，計算結(jié)果準確預(yù)測了跨聲速顫振的“凹坑”現(xiàn)象；在亞跨聲速范圍內(nèi)，計算結(jié)果與試驗結(jié)果吻合良好；在超聲速區(qū)域，計算結(jié)果明顯高于試驗結(jié)果，超聲速區(qū)域計算結(jié)果與其他文獻的計算結(jié)果類似[22-23]。這個具體算例的數(shù)值結(jié)果一方面驗證了本文超大規(guī)模顫振問題的數(shù)值模擬能力；另一方面，網(wǎng)格規(guī)模的增加并沒有使得計算結(jié)果與試驗結(jié)果更加吻合。計算結(jié)果與試驗結(jié)果的差異需要從風(fēng)洞試驗以及數(shù)值模擬方法兩個方面進一步開展工作。

圖13 AGARD445.6機翼顫振特性數(shù)值模擬Fig.13 Flutter numerical simulation for AGARD445.6 wing

5 結(jié) 論

建立了靜氣動彈性、顫振模擬的一般流程并開發(fā)了超大規(guī)模氣動彈性數(shù)值模擬軟件的初級版本，實現(xiàn)了百億量級超大規(guī)模靜態(tài)結(jié)構(gòu)網(wǎng)格生成和十億量級結(jié)構(gòu)網(wǎng)格自動變形，開展了動網(wǎng)格模塊、CFD求解模塊、顫振模擬的并行效率和并行一致性測試，采用CRM-WBH構(gòu)型和AGARD445.6機翼模型進一步開展了氣動彈性數(shù)值模擬軟件的確認工作。

下一步工作主要包括：采用MPI+OpenMP并行策略，進一步提高氣動彈性數(shù)值模擬的并行效率；開展時域非定常氣動彈性計算方法研究，提高非定常氣動彈性問題數(shù)值模擬精度；利用高保真度靜氣動彈性、顫振風(fēng)洞試驗數(shù)據(jù)，進一步開展氣動彈性模擬軟件的確認工作。

致謝:感謝參加本課題的張書俊、李偉、楊小川等三位同志在超大規(guī)模網(wǎng)格生成、結(jié)構(gòu)有限元分析和計算結(jié)果處理等方面的工作。