王運(yùn)濤, 孟德虹, 孫 巖, 洪俊武

(1. 中國空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所, 四川 綿陽 621000;2. 中國空氣動(dòng)力研究與發(fā)展中心 空氣動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室, 四川 綿陽 621000)

0 引 言

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

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

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

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

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

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

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

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

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

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

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

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

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

圖2 動(dòng)氣動(dòng)彈性計(jì)算流程Fig.2 Outline of dynamic aero-elastic problem simulation

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

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

超大規(guī)模靜態(tài)網(wǎng)格生成技術(shù)是實(shí)現(xiàn)復(fù)雜構(gòu)型超大規(guī)模氣動(dòng)彈性問題數(shù)值模擬的前提。以包含機(jī)翼/機(jī)身/掛架/短艙/平尾/立尾的典型運(yùn)輸機(jī)巡航構(gòu)型為例，目前模擬上述復(fù)雜構(gòu)型的結(jié)構(gòu)網(wǎng)格規(guī)模一般在千萬量級(jí)，在配置較高的微機(jī)或工作站上利用商業(yè)網(wǎng)格生成軟件即可完成。而構(gòu)造數(shù)億量級(jí)乃至百億量級(jí)的結(jié)構(gòu)網(wǎng)格已經(jīng)遠(yuǎn)遠(yuǎn)超出了當(dāng)前單機(jī)能力范圍，必須采用新的網(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)造千萬量級(jí)高質(zhì)量初始結(jié)構(gòu)網(wǎng)格；2) 利用初始結(jié)構(gòu)網(wǎng)格拓?fù)浞謮K信息，建立網(wǎng)格窗口的信息數(shù)據(jù)庫；3)基于網(wǎng)格窗口的信息數(shù)據(jù)庫創(chuàng)建不同網(wǎng)格窗口之間的關(guān)聯(lián)關(guān)系；4) 根據(jù)網(wǎng)格規(guī)模需求，對(duì)棱線網(wǎng)格點(diǎn)進(jìn)行加密，重新計(jì)算棱線網(wǎng)格點(diǎn)分布及棱線網(wǎng)格點(diǎn)坐標(biāo)，通過超限插值(TFI)生成加密后的窗口網(wǎng)格；5) 采用線性投影方法，對(duì)邊界屬性標(biāo)記為物面的窗口網(wǎng)格重新向物面投影，獲得加密后的物面網(wǎng)格；6) 利用三維TFI生成加密后的體網(wǎng)格。上述超大規(guī)模網(wǎng)格生成策略與待模擬外形的復(fù)雜程度無關(guān)，所構(gòu)造的超大規(guī)模結(jié)構(gòu)網(wǎng)格拓?fù)浣Y(jié)構(gòu)和網(wǎng)格單元質(zhì)量與初始網(wǎng)格完全相同，并可以對(duì)初始網(wǎng)格進(jìn)行不同程度的粗化或細(xì)化。采用上述方法，可以快速獲得百億量級(jí)以上的高質(zhì)量結(jié)構(gòu)網(wǎng)格，同時(shí)可以構(gòu)造一系列的不同規(guī)模網(wǎng)格開展網(wǎng)格收斂性研究。

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

圖3 CRM-WBH構(gòu)型的網(wǎng)格拓?fù)?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)域的研究熱點(diǎn)，也是氣動(dòng)彈性數(shù)值模擬的關(guān)鍵技術(shù)之一。高質(zhì)量、魯棒高效的網(wǎng)格變形技術(shù)應(yīng)滿足以下要求：1)自動(dòng)化。減少甚至避免網(wǎng)格變形過程中的人工干預(yù)是氣動(dòng)彈性問題高效數(shù)值模擬的基本要求。2)并行化。網(wǎng)格變形過程的并行化是提高網(wǎng)格變形效率有效途徑，對(duì)于超大規(guī)模網(wǎng)格的動(dòng)氣動(dò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)實(shí)現(xiàn)復(fù)雜構(gòu)型結(jié)構(gòu)網(wǎng)格的自動(dòng)變形，利用刪減基點(diǎn)后的徑向基函數(shù)插值更新結(jié)構(gòu)網(wǎng)格塊的棱線坐標(biāo)，然后利用TFI插值重新生成計(jì)算網(wǎng)格。該方法的具體實(shí)施步驟是：1)根據(jù)輸入的多塊結(jié)構(gòu)網(wǎng)格信息，在物面網(wǎng)格處選擇徑向基函數(shù)基點(diǎn)，構(gòu)造可逆的插值矩陣；2)根據(jù)輸入的物面變形確定物面徑向基基點(diǎn)位移，利用可逆的插值矩陣確定徑向基插值系數(shù)；3)采用徑向基插值方法更新物面與空間網(wǎng)格的網(wǎng)格棱線坐標(biāo)；4)采用二維、三維TFI方法構(gòu)造各個(gè)網(wǎng)格塊的面網(wǎng)格與體網(wǎng)格。RBF_TFI方法憑借優(yōu)異的性能在很多程序和算例中得到了廣泛應(yīng)用，但在一些比較極端的情況下，邊界層網(wǎng)格第一層高度十分小的時(shí)候(這在高雷諾數(shù)流動(dòng)問題和超大規(guī)模網(wǎng)格問題中經(jīng)常遇到)，RBF_TFI容易引起邊界層區(qū)域網(wǎng)格棱線的扭曲現(xiàn)象，并造成物面附近網(wǎng)格的交叉破壞問題，從而導(dǎo)致流場計(jì)算無法繼續(xù)進(jìn)行?；诓煌憷龑?duì)該現(xiàn)象進(jìn)行了深入的分析，認(rèn)為是物面網(wǎng)格點(diǎn)與空間網(wǎng)格點(diǎn)采用不同的坐標(biāo)更新方法導(dǎo)致的變形不一致引起了網(wǎng)格交錯(cuò)破壞。通過現(xiàn)象、原因等方面的詳細(xì)分析，提出一種解決該問題的改進(jìn)方法，改進(jìn)后的網(wǎng)格更新流程見圖4。

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

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

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

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

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

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

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

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

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

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

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

表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 靜氣動(dòng)彈性模塊并行效率及一致性測試

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

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

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

4 氣動(dòng)彈性數(shù)值模擬軟件的確認(rèn)

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

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

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

圖7 CRM-WBHS構(gòu)型網(wǎng)格拓?fù)浼氨砻婢W(wǎng)格Fig.7 Grid topology and surface grid for CRM-WBHS configuration

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

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

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

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

(a) CL～α

(b) CD～α

(c) Cm～α

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

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

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

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

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

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

5 結(jié) 論

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

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

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