趙 磊，高麗敏，俞一波，李瑞宇，王可鑫

（西北工業(yè)大學(xué)動力與能源學(xué)院，西安 710021）

0 引言

數(shù)字孿生技術(shù)是未來工業(yè)極具潛力和挑戰(zhàn)的發(fā)展方向，依托該技術(shù)所構(gòu)建形成的航空發(fā)動機(jī)數(shù)字孿生體[1]，對于解決航空發(fā)動機(jī)研制中面臨的多系統(tǒng)、多部件、多物理場耦合等挑戰(zhàn)性問題提供了新的思路。壓氣機(jī)作為航空發(fā)動機(jī)的3 大核心部件，其數(shù)字孿生體的研制對于形成航發(fā)整機(jī)的孿生體具有重要意義。

以壓氣機(jī)為物理實(shí)體構(gòu)建孿生模型，首先要保證壓氣機(jī)的數(shù)值模型能夠與現(xiàn)實(shí)環(huán)境中的物理實(shí)體一一對應(yīng)，即壓氣機(jī)數(shù)值模型具備完備性要求。“孿生”技術(shù)最早由美國國家航空航天局（NASA）在阿波羅項(xiàng)目[2]中提出；Grieves 等[3]指出孿生體是物理實(shí)體在虛擬空間的“數(shù)字復(fù)制品”；NASA 后期[4]發(fā)布了數(shù)字孿生的路線圖，其中數(shù)值仿真模擬為該技術(shù)的主要部分；Tuegel 等[5]、Stargel 等[6]、Cerrone 等[7]也指出數(shù)字孿生體需要對復(fù)雜系統(tǒng)中不同空間和時間的事件進(jìn)行準(zhǔn)確數(shù)值仿真預(yù)測；孫明霞等[8]、張志博等[9]、劉永泉等[10]以及王樂等[11]通過分析現(xiàn)代數(shù)字孿生技術(shù)，指出航空發(fā)動機(jī)的數(shù)字孿生技術(shù)中最核心的部分為孿生模型相關(guān)的數(shù)值模型的構(gòu)建?？梢姡瑪?shù)字孿生技術(shù)以數(shù)值仿真作為其核心技術(shù)手段。然而，在現(xiàn)有航空發(fā)動機(jī)研制過程中，壓氣機(jī)的數(shù)值仿真多以定常流動、均勻進(jìn)氣、單通道周期性假設(shè)等建立數(shù)值模型，所建立的數(shù)字孿生體不能完全反映壓氣機(jī)在真實(shí)環(huán)境下的工作性能等參數(shù)變化。因此，從計(jì)算流體力學(xué)角度出發(fā)構(gòu)建壓氣機(jī)數(shù)字孿生體，其孿生模型需要具備開展真實(shí)邊界條件下多級全環(huán)的非定常數(shù)值仿真能力。然而，壓氣機(jī)數(shù)字孿生模型中的全環(huán)非定常的計(jì)算需求巨大。Garbor 等[12]提出數(shù)字孿生的計(jì)算成本在于仿真所需的大量計(jì)算資源；Yamada 等[13]對某燃?xì)廨啓C(jī)7 級壓氣機(jī)開展全環(huán)計(jì)算，使用了近7300 多CPU核心的計(jì)算資源；Gourdain 等[14]指出未來采用大渦模擬開展多級壓氣機(jī)全環(huán)計(jì)算需要具有E級（運(yùn)算峰值性能達(dá)到百億億次，即1018flop/s）計(jì)算能力的高性能計(jì)算機(jī)；Slotnick 等[15]指出航空發(fā)動機(jī)整機(jī)在非設(shè)計(jì)工況流動的數(shù)值模擬（其中包含壓氣機(jī)多級全環(huán)非定常數(shù)值模擬）將是未來高性能并行計(jì)算領(lǐng)域的重大挑戰(zhàn)課題；Belzner等[16]曾提出數(shù)字孿生體中事件的模擬速度要盡可能比現(xiàn)實(shí)中的事件發(fā)生速度更快?？梢?，針對未來壓氣機(jī)數(shù)字孿生技術(shù)的發(fā)展必然要借助高性能并行計(jì)算技術(shù)，以提升數(shù)值仿真手段的計(jì)算能力和速度。

本文結(jié)合數(shù)值仿真方法給出了航發(fā)數(shù)字孿生技術(shù)中數(shù)值模型完備性的定義；針對某S 彎進(jìn)氣道中涵道風(fēng)扇的數(shù)字孿生體，對比分析了不同數(shù)值模型對孿生結(jié)果的影響規(guī)律；針對現(xiàn)代高性能計(jì)算平臺，給出了壓氣機(jī)數(shù)字孿生模型中的2 層次并行計(jì)算模型，并開展了多級壓氣機(jī)全環(huán)非定常計(jì)算驗(yàn)證。

1 航發(fā)孿生技術(shù)中數(shù)值模型完備性定義

航空發(fā)動機(jī)的數(shù)字孿生技術(shù)是集成多學(xué)科、多物理場、多尺度、概率性的仿真系統(tǒng)。航發(fā)孿生體是航空發(fā)動機(jī)的物理過程在數(shù)字虛擬世界中的完美映射，該孿生體由許多相互關(guān)聯(lián)的數(shù)值模型組成，各模型綜合起來用以描述發(fā)動機(jī)運(yùn)行過程中的氣動、傳熱、燃燒、強(qiáng)度等物理過程及其性能。要保證數(shù)字孿生體的正確性，首先需要滿足組成孿生體相關(guān)不同數(shù)值模型的完備性要求。由此，定義針對任意給定誤差δ≥0 和可描述航空發(fā)動機(jī)實(shí)時物理過程的控制方程Φ(x?)及其定解條件，若存在一個數(shù)值可解的計(jì)算模型φ(x?)使得

在給定解條件下的所有合法輸入x?均成立，那么認(rèn)為該數(shù)值可解的計(jì)算模型及其定解條件對于航發(fā)數(shù)字孿生體是計(jì)算完備的，即數(shù)值模型具備完備性要求。工程實(shí)踐中多采用計(jì)算模型的結(jié)果與真實(shí)環(huán)境下試驗(yàn)結(jié)果之間的誤差來定義δ。

壓氣機(jī)中氣動物理過程可用N-S方程及定解條件對其進(jìn)行描述，相應(yīng)的數(shù)值模型可結(jié)合現(xiàn)代計(jì)算流體力學(xué)方法（Computational Fluid Dynamics，CFD）給出。現(xiàn)有的CFD 中以有限體積方法居多，該方法將連續(xù)N-S方程及其定解條件離散化為物理空間中多個網(wǎng)格計(jì)算單元上的代數(shù)方程組，通過數(shù)值求解該方程組可近似獲得壓氣機(jī)中的氣動物理過程。因此，從計(jì)算流體力學(xué)角度出發(fā)，滿足壓氣機(jī)數(shù)值模型完備性的要求不僅要保證離散化連續(xù)方程過程中引入的離散誤差、舍入誤差等盡可能小，同時描述氣體流動湍流狀態(tài)等物理模型的誤差也要盡可能小。此外，幾何模型、邊界條件、初始條件等也要滿足一定的誤差要求。綜上，在滿足給定誤差δ的前提下，采用CFD 方法的壓氣機(jī)數(shù)值模型具備對壓氣機(jī)氣動物理過程的計(jì)算完備性要求。

2 某涵道風(fēng)扇/壓氣機(jī)數(shù)字孿生中數(shù)值模型完備性分析

選取某帶S 彎進(jìn)氣道的涵道風(fēng)扇/壓氣機(jī)（如圖1所示）作為研究對象，氣流過S 彎進(jìn)氣道后，在進(jìn)氣道與涵道風(fēng)扇/壓氣機(jī)的氣動交界面（Aerodynamic Inter?face Plane，AIP）上形成總壓畸變圖譜，如圖2 所示。對該涵道風(fēng)扇/壓氣機(jī)建立數(shù)字孿生模型可以采用3種不同數(shù)值模型方案（暫不考慮進(jìn)氣道和壓氣機(jī)的耦合效應(yīng)）：（1）將AIP 界面總壓進(jìn)行平均后取平均總壓作為進(jìn)氣邊界，即可采用現(xiàn)有的均勻進(jìn)氣周期性假設(shè)的單通道計(jì)算模型，如圖3 所示；（2）由于AIP 界面主要以總壓畸變?yōu)橹鳎乙灾芟蚩倝夯優(yōu)橹?，因此可從原始AIP 界面總壓畸變圖譜上將其簡化為周向總壓畸變，如圖4 所示。其中，從AIP 界面測得周向畸變角度為130.9°，總壓畸變強(qiáng)度為8%。選取簡化后的周向總壓畸變圖譜作為壓氣機(jī)進(jìn)氣邊界條件，并采用全通道幾何模型，如圖5 所示；（3）直接將AIP 界面的總壓畸變圖譜（圖2）給定為壓氣機(jī)進(jìn)口邊界條件，同時采用全通道計(jì)算模型（圖5）。以上3 種數(shù)值模型的邊界條件和幾何模型見表1。

表1 涵道風(fēng)扇不同數(shù)字孿生模型對比

圖1 帶S彎進(jìn)氣道涵道風(fēng)扇/壓氣機(jī)

圖2 AIP界面總壓畸變圖譜

圖3 均勻進(jìn)氣與周期性假設(shè)的單通道數(shù)值模型

圖4 通過AIP界面獲得的簡化周向總壓畸變圖譜

圖5 涵道風(fēng)扇/壓氣機(jī)全通道計(jì)算模型

3 種孿生模型獲得的壓氣機(jī)特性如圖6 所示，其中數(shù)值方法以及網(wǎng)格無關(guān)性驗(yàn)證見文獻(xiàn)[17]。通過對比3 種不同數(shù)值模型可見，在相同轉(zhuǎn)速工況下，采用均勻進(jìn)氣和單通道假設(shè)的數(shù)值模型1 的壓比和效率要高于模型2、3 的，其中模型1 最高效率比模型2的高3.8%，比模型3 的高5.5%，同時相應(yīng)壓比也有不同程度的差異。此外從所獲得的壓氣機(jī)穩(wěn)定裕度可以看出，采用模型1 所得到的該涵道風(fēng)扇的穩(wěn)定裕度要遠(yuǎn)大于采用模型2、3 所獲得的穩(wěn)定裕度，其中采用模型2 所獲得的穩(wěn)定裕度僅為采用模型1 所獲得的50.5%，裕度范圍幾乎縮小一半。同樣的，對比模型2、3 可見，數(shù)值模型都采用全通道計(jì)算模型，僅是進(jìn)氣邊界的不同，采用模型3 所獲得的穩(wěn)定裕度范圍已經(jīng)縮小到采用模型2所獲得的29.4%。如果假設(shè)模型3（AIP 界面畸變圖譜為進(jìn)氣邊界條件，采用全通道計(jì)算模型）是該涵道風(fēng)扇的完備數(shù)值模型，通過對比可見，采用傳統(tǒng)單通道周期性和均勻進(jìn)氣假設(shè)所獲得的壓氣機(jī)性能整體參數(shù)將完全不適合作為該壓氣機(jī)的數(shù)字孿生模型。

圖6 3種孿生模型獲得的壓氣機(jī)特性

全通道模型不同進(jìn)氣邊界條件下涵道風(fēng)扇出口總壓分布如圖7 所示。二者均采用全通道計(jì)算模型，僅是進(jìn)口邊界條件不同。對比出口總壓分布云圖可見，不同進(jìn)口總壓畸變分布（兩者總壓畸變強(qiáng)度均相同），所得到的涵道風(fēng)扇出口的總壓分布存在較大差異，其中采用模型3 出口的高壓區(qū)的分布不均勻性明顯高于采用模型2 的，其壓力分布的不均勻更加明顯。由此可見，在數(shù)字孿生模型中，即使幾何模型一致，進(jìn)、出口邊界條件的差異（保持總壓畸變強(qiáng)度一定）也會使得進(jìn)氣畸變經(jīng)過風(fēng)扇轉(zhuǎn)子和靜子傳遞后，內(nèi)部流場特征存在較大差異。

圖7 全通道模型不同進(jìn)氣邊界條件下涵道風(fēng)扇出口總壓分布

綜合以上分析可知，涵道風(fēng)扇孿生模型的不同對“現(xiàn)實(shí)環(huán)境”中的風(fēng)扇壓氣機(jī)總體性能參數(shù)影響非常大，孿生模型即壓氣機(jī)數(shù)值計(jì)算模型的不同會對壓氣機(jī)本身數(shù)字孿生體的內(nèi)部流場以及宏觀性能參數(shù)會產(chǎn)生巨大影響。針對進(jìn)氣畸變的孿生場景研究中發(fā)現(xiàn)，幾何模型（單通道周期性假設(shè)）的不同對于所獲得的性能總體參數(shù)的影響最大。

3 數(shù)字孿生模型中的高效計(jì)算技術(shù)

航發(fā)孿生技術(shù)對高性能計(jì)算技術(shù)有著較高要求，如何充分挖掘現(xiàn)有高性能計(jì)算潛能以推動數(shù)值孿生技術(shù)的發(fā)展，是未來航發(fā)數(shù)字孿生需要突破的瓶頸之一。一方面，計(jì)算機(jī)技術(shù)的發(fā)展正在促使高性能的浮點(diǎn)計(jì)算性能朝著E 級（百億億次浮點(diǎn)計(jì)算能力）計(jì)算發(fā)展，甚或是未來量子計(jì)算成為可能，都有望促使基于3 維的物理過程仿真與現(xiàn)實(shí)世界同步；另一方面，對諸如壓氣機(jī)等航發(fā)部件實(shí)現(xiàn)部件級物理過程仿真（例如實(shí)現(xiàn)壓氣機(jī)整機(jī)全環(huán)非定常數(shù)值模擬），對于航發(fā)孿生體不同部件在不同工作環(huán)境下數(shù)字孿生體的準(zhǔn)確構(gòu)建與驗(yàn)證具有重要意義。因此，針對壓氣機(jī)，未來壓氣機(jī)數(shù)字孿生模型不僅需要具備多級全環(huán)非定常的并行計(jì)算能力，同時相應(yīng)仿真計(jì)算能力的構(gòu)建需要能夠充分挖掘高性能計(jì)算的浮點(diǎn)計(jì)算潛力。

現(xiàn)有高性能計(jì)算平臺均采用分布式計(jì)算系統(tǒng)，大都具有機(jī)柜、節(jié)點(diǎn)、多CPU、多計(jì)算核心等不同層次的硬件架構(gòu)特征。本文結(jié)合該計(jì)算架構(gòu)特征，發(fā)展了針對壓氣機(jī)多級計(jì)算的兩層次并行計(jì)算模型，2 層次并行計(jì)算模型如圖8 所示。結(jié)合計(jì)算系統(tǒng)中計(jì)算進(jìn)程、計(jì)算線程在不同層次架構(gòu)中的計(jì)算特征，可將上述硬件結(jié)構(gòu)特征抽象為2 層計(jì)算模型，上層主要以計(jì)算進(jìn)程為主，下層主要以計(jì)算線程為主。針對該2 層計(jì)算特征，可分別采用不同的并行編程模型。例如針對計(jì)算進(jìn)程特征層，一般采用MPI 編程模型；針對計(jì)算線程特征層，可結(jié)合具體的硬件架構(gòu)特征，采用OPENMP、CUDA 等編程模型。計(jì)算進(jìn)程可采用區(qū)域分解算法對計(jì)算對象進(jìn)行任務(wù)級分解，可充分發(fā)揮不同節(jié)點(diǎn)、不同CPU 上的計(jì)算潛力。計(jì)算線程層次主要針對現(xiàn)有的異構(gòu)加速硬件，針對任務(wù)級中的計(jì)算密集型任務(wù)實(shí)現(xiàn)高效加速。

圖8 2層次并行計(jì)算模型

結(jié)合上述2 層次高性能計(jì)算模型，本文自主發(fā)展了相應(yīng)的高性能計(jì)算軟件，對某15 級壓氣機(jī)開展了全環(huán)非定常計(jì)算驗(yàn)證。15 級壓氣機(jī)全通道計(jì)算網(wǎng)格如圖9所示。針對該15級壓氣機(jī)的全環(huán)模型，轉(zhuǎn)靜交界面的處理采用滑移網(wǎng)格法，計(jì)算時網(wǎng)格量約為1.5×109，每排葉排網(wǎng)格數(shù)約為4.6×107，單個葉片通道的網(wǎng)格量大于106，y+為1。共使用256 個計(jì)算節(jié)點(diǎn)進(jìn)行并行計(jì)算，每個計(jì)算節(jié)點(diǎn)為2 路CPU，每顆CPU 均具有48個計(jì)算核心。該多級壓氣機(jī)全環(huán)非定常計(jì)算采用3072 個計(jì)算進(jìn)程，每個計(jì)算進(jìn)程采用4 個計(jì)算線程實(shí)現(xiàn)加速計(jì)算。

圖9 15級壓氣機(jī)全通道計(jì)算網(wǎng)格

15 級壓氣機(jī)軸向速度瞬時流場分布如圖10 所示。為了深入對比不同時刻多級壓氣機(jī)內(nèi)部流場的非定常變化，前4級壓氣機(jī)50%葉不同時刻的壓力分布如圖11 所示。從圖中可見，受到上下游葉片排的相互影響，在葉片排軸向間隙中的流場變化隨著動葉的旋轉(zhuǎn)都在變化。第3級動葉排出口靜壓分布隨時間的變化如圖12 所示。從圖中第3 級動葉排出口的壓力分布可見，不同時刻第3 級動葉排出口的壓力分布的周向不均勻性變化明顯，且在不同時刻的壓力分布的不均勻性均存在。第10~12 級壓氣機(jī)50%葉高在不同時刻的壓力分布和第11 級動葉排出口靜壓分布隨時間的變化如圖13、14 所示。對比前面級的仿真結(jié)果可見，多級計(jì)算中不同葉排間、單葉排不同葉道在不同時刻的流場結(jié)果均呈現(xiàn)出差異，若僅是采用單通道定常假設(shè)的數(shù)值模型是無法獲得該流場的細(xì)節(jié)差異。

圖10 15級壓氣機(jī)軸向速度瞬時流場分布

圖11 前4級壓氣機(jī)50%葉不同時刻的壓力分布

圖12 第3級動葉排出口靜壓分布隨時間的變化

圖13 第10~12級壓氣機(jī)50%葉高在不同時刻的壓力分布

圖14 第11級動葉排出口靜壓分布隨時間的變化

綜合上述分析可見，壓氣機(jī)數(shù)字孿生模型要具備多級全環(huán)非定常計(jì)算能力離不開高性能計(jì)算技術(shù)支撐。同時結(jié)合現(xiàn)有高性能計(jì)算平臺發(fā)展適合于壓氣機(jī)整機(jī)非定常計(jì)算的并行計(jì)算模型，對于發(fā)掘壓氣機(jī)內(nèi)部流動的非定常流動特性極具必要性。

4 總結(jié)

（1）針對S 彎進(jìn)氣道中涵道風(fēng)扇建立了數(shù)字孿生的不同數(shù)值模型，主要從進(jìn)氣邊界模型和幾何模型2方面進(jìn)行了分析。通過研究發(fā)現(xiàn)，采用均勻進(jìn)氣邊界和單通道周期性假設(shè)計(jì)算模型與采用AIP 界面原始畸變圖譜、全通道計(jì)算模型所獲得的壓氣機(jī)性能參數(shù)差別巨大，尤其是壓氣機(jī)穩(wěn)定裕度范圍完全不同?？梢姡瑥臄?shù)值模型完備性角度出發(fā)，壓氣機(jī)數(shù)字孿生模型需要同時在均勻進(jìn)氣和畸變進(jìn)氣環(huán)境下獲得準(zhǔn)確結(jié)果時，采用全環(huán)計(jì)算模型構(gòu)建壓氣機(jī)數(shù)字孿生體較為合適。

（2）壓氣機(jī)數(shù)字孿生模型需要具備多級全環(huán)非定常計(jì)算能力，其對計(jì)算量的需求巨大。本文結(jié)合現(xiàn)代高性能計(jì)算技術(shù)發(fā)展了2 層次大規(guī)模并行計(jì)算模型，實(shí)現(xiàn)了15 級壓氣機(jī)全環(huán)非定常模擬仿真，計(jì)算規(guī)模使用到了12288 個計(jì)算核心，獲得了多級壓氣機(jī)的整機(jī)非定常流場特征，為后續(xù)形成壓氣機(jī)孿生模型奠定了基礎(chǔ)。