考慮不確定性的數值模擬模型確認方法

任成坤, 熊芬芬, 王元豪, 李娜, 姜波, 郭志平

(1.西南技術工程研究所, 重慶 400039; 2.北京理工大學 宇航學院, 北京 100081)

對于基于仿真的優化設計,構建高保真、高可信的數值模擬模型是重中之重。使用與真實結果存在較大差異的數值模擬將造成預測結果不準確,對于具有高可靠性要求的飛行器而言,極有可能引入潛在風險,因此開展飛行器優化設計前必須對數值模擬進行模型確認和可信度評估。然而,隨著所模擬的物理過程日益復雜,一次性建立適合所有工況的精確仿真模型變得不太現實。比如,對于計算流體力學(computational fluid dynamics,CFD)數值模擬,現有研究表明商業或開源CFD軟件中模型參數默認值或文獻中給出的推薦值并不一定能使一般流動條件下的模型預測效果達到最優[1],而且機翼表面的壓力系數、激波在翼型上表面的位置和激波后壓力等對湍流模型封閉系數的變化非常敏感[2]。

傳統的模型確認以參數型修正為主,通過直接在確定性條件下使數值模擬與試驗數據之間的偏差最小來實現。該方式直接將模型參數作為修正對象,物理意義明確,是目前模型確認研究的主流方向[3]。王紀森等[4]針對油液流動的CFD數值模擬,采用參數遍歷法對Realizablek-ε兩方程湍流模型中的參數進行修正,使得壓力損失的仿真結果與試驗結果接近。張珍等[5]對雷諾平均方程(Reynolds-averaged Navier-Stokes,RANS)求解中的渦黏系數進行多次修正,提高了模型對流場各向異性特征預測的性能。

然而,實際的數值模擬中往往存在大量不確定性,如CFD數值模擬存在來流和邊界條件、湍流模型系數等不確定性,這對數值模擬輸出結果有著不可忽視的影響。例如,飛行器跨聲速飛行時,流場具有強非線性,幾何外形變化對氣動特性影響更加復雜[6]。實際中在不確定性的影響下,確定性情況下確認好的仿真模型極有可能產生很大的預測偏差,基于確定性模擬結果得出的結論有可能導致真實系統達不到預期的性能要求,引入潛在風險。因此,迫切需要考慮并高效量化不確定性對數值模擬結果的影響,開展不確定性條件下的模型修正。不確定性量化作為數值模擬模型確認的核心之一,目前圍繞其已經開展了大量研究工作并取得了豐碩的研究成果,為后續模型修正的開展提供了必要條件。Liu等[7]采用蒙特卡羅仿真(Monte Carlo simulation,MCS)、混沌多項式(polynomial chaos,PC)、徑向基函數和梯度增強Kriging等方法研究了跨聲速條件下機翼外形加工誤差對氣動特性的影響,研究表明,梯度增強的降階模型方法效率更高。Loeven等[8]通過對參數化建模后的機翼使用概率配置點法,研究了翼型的厚度、最大彎度位置和最大彎度等關鍵設計變量的不確定性對翼型氣動特性的影響。Geneva等[9]采用深度學習技術對RANS方程中的模型不確定性進行量化,并提供了諸如速度和壓力等輸出的概率區間。陳江濤等[10]使用混沌多項式方法研究了湍流模型系數的不確定度對RAE2822跨聲速翼型繞流模擬的影響,計算中關注了數值模擬的積分量(升力系數和阻力系數)和局部量(壁面壓力、摩擦因數和空間馬赫分布)的不確定度量化結果。

綜上可見,目前圍繞數值模擬的模型確認,大量的研究工作均未考慮模擬中的不確定性。雖然目前已開展了很多CFD不確定性量化的研究工作,這些工作為梳理各類不確定性對流動特性的影響規律及其大小提供了豐富的思路,給后續CFD數值模擬的模型修正工作的開展奠定了良好的基礎。但是,面對實際中的高維不確定性,不確定性量化面臨“維數災難”難題。同時,如何基于不確定性量化的結果進行高效的模型參數修正,構建高可信的預測模型,打破傳統依賴經驗或試湊的低效模型確認模式,形成不確定性下數值模擬模型確認的閉環流程,還較少有研究報道。因此,本文提出一套集不確定性量化、靈敏度分析、確認度量、參數修正為一體的考慮不確定性的數值模擬模型確認流程,并將提出的方法應用于NACA0012和ONERA M6翼型的CFD數值模擬模型確認。

1 基于優質小樣本的模型確認方法

為提高數值模擬的可信度和預測能力,為后續開展分析和優化設計提供高保真的仿真模型,本文提出一種集不確定性量化、靈敏度分析、確認度量、參數修正為一體的數值模擬模型確認方法。圖1展示了所提出方法的實施流程,其具體實現步驟如下。

圖1 提出的模型確認方法流程圖

步驟1初始數值模擬模型構建,確定不確定性參數分布類型及其分布參數。

確定數值模擬的邊界條件等,確定數值模擬仿真模型及模型參數的不確定性分布類型和分布參數,如CFD數值模型中的湍流模型及其封閉系數。

步驟2不確定性量化,得到不確定性影響下數值模擬輸出響應的不確定性。

基于給定的不確定性分布信息,對數值模擬的輸出響應進行不確定性量化,針對低維問題(≤8維)推薦使用混沌多項式方法(見1.1.1小節),針對高維問題(>8維)采用深度學習方法(見1.1.2小節)。上述維度劃分是基于文獻[11-12]中的定義以及對混沌多項式與深度學習的特性和使用習慣經測試后得到的。

步驟3結合試驗數據進行數值模擬模型的確認度量(見1.2小節),若模型不滿足要求(即:確認度量指標>ε),則進入步驟4,否則模型確認程序結束,可基于當前數值模擬模型進行后續的仿真分析和優化設計。

步驟4全局靈敏度分析, 發掘重要的修正參數。

應用全局靈敏度分析方法(見1.3小節)量化模型參數對輸出響應的重要程度,并進行重要性排序,從中選取影響較大的參數進行后續修正。

步驟5基于優質小樣本的參數修正或模型重選。

通過結合試驗數據并應用優質小樣本方法(見1.4小節)選擇與試驗數據更加吻合的數值模擬樣本,并逆向反推得到相應的不確定性模型參數分布范圍或分布參數,實現參數修正。若參數修正迭代次數大于給定閾值K時仍無法獲取滿足確認度量準則的數值模擬,則可進行模型重選,例如對于CFD數值模擬,可選擇最佳的湍流模型。利用貝葉斯定理計算各個候選模型的后驗分布,選擇與試驗數據更吻合的模型。也可采用貝葉斯平均方法對多個候選模型進行融合,共同實現響應預測。

步驟6重復步驟2～5,直至確認度量指標滿足要求(EMR≤ε),則可利用最新的數值模擬模型進行分析和優化設計。

1.1 不確定性量化方法

從圖1中可以發現,不確定性傳播是必不可少的環節,它在整個模型確認流程中起著舉足輕重的作用。不確定性傳播的方法有很多,如蒙特卡羅仿真(Monte Carlo simulation,MCS)、混沌多項式(polynomial chaos,PC)和深度學習(deep learning,DL)等。PC方法理論完備且具有指數收斂速度,近年來在數值模擬的不確定性量化和優化設計中得到廣泛研究和應用[13]。但PC方法在非線性極強的高維問題上容易出現收斂緩慢甚至無法收斂的問題。而深度學習(deep learning,DL)技術由于其對復雜高維函數具有較強的表征能力,越來越多的復雜工程問題可通過深度學習得以解決[14]。深度學習技術的應用,相當于構建了高維函數的代理模型,不確定性量化則在代理模型上進行,為突破高維不確定性量化的“維數災難”提供了有效解決途徑。因此在本節中介紹應用最為廣泛的PC方法和深度學習方法。

1.1.1 基于混沌多項式的不確定性量化

p階PC模型可寫為[15-16]

(1)

式中:ξ為標準隨機變量;b為PC系數,正交多項式Φi(ξ)是各維隨機分布對應的一維正交多項式的乘積,即

(2)

本文采用隨機響應面法(stochastic response surface method,SRSM)求解PC系數,將標準隨機樣本ξS和相應的真實函數響應值G代入(2)式中PC模型兩端,得

(3)

(3)式可簡寫為

Ab=G

(4)

基于(4)式,利用最小二乘法計算PC系數。

b=(ATA)-1ATG

(5)

當得到PC系數后,可用解析表達式(6)～(7)非常方便地計算輸出響應y的前兩階統計矩。

1.1.2 基于貝葉斯深度學習的不確定性量化

貝葉斯深度神經網絡(Bayesian neural networks,BNN)作為深度神經網絡的一個分支,具有穩健性高且能避免過擬合問題的優勢[17]。

對于一個具體問題來說,當給定d維輸入向量x,可定義一個具有L層的深度神經網絡

fθ(x)=wLφθ-(x)+bL

(8)

式中:wL和bL是網絡第L層(即輸出層)的權值和偏差;φθ-(x)是輸入變量x經前L-1層網絡運算后的迭代值;θ-表示前L-1層全部的網絡參數;fθ表示最后得到的網絡參數為θ的深度神經網絡。

fθ(x)=wφθ-(x)

(9)

式中,φθ-(x)是倒數第二層的輸出向量。由于神經網絡的計算速度極快,可通過直接在BNN上運行MCS的方法進行不確定性量化。

多分支網絡可以更好地應對存在相關性的多輸入和多輸出問題,其分支層可以根據問題的不同構建在輸入層或輸出層處。圖2展示了一個有2層隱藏層的多分支神經網絡結構圖,其中輸入層與前2層隱藏層為骨干層,輸出層為分支層。

圖2 多分支神經網絡結構

1.2 數值模擬模型的確認度量

通常情況下由于試驗數據獲取成本非常昂貴,往往僅能得到少量幾個試驗數據,此時可采用基于距離的確認度量指標,具體介紹如下。

(10)

如果存在若干個不同的工況需要同時進行模型確認,則可對多個工況下的EMR累加求和,作為模型確認度量的指標。

除此之外,若獲取到的試驗數據較多(ne較大),可有效構建輸出響應的經驗概率分布模型,則可采用面積法(Area Metric,AM)作為確認度量準則[20],其計算公式如下。

(11)

式中:y為數值模擬的輸出響應;Fa(y)為根據數值模擬所得輸出樣本所構建的y累積分布函數(cumulative distribution function,CDF);Fe(y)為根據試驗測量樣本ye所構建的y的經驗累積分布函數(empirical cumulative distribution function,ECDF)。

1.3 全局靈敏度分析

基于不確定性量化的結果,開展全局靈敏度分析。傳統的處理隨機不確定性的全局靈敏度分析方法大體可以分為基于回歸的方法和基于方差的方法[21],其中基于全局靈敏度指標的方差分析法因其簡單有效的特點得到了確認應用。對于一個有著d維輸入向量x=[x1,…,xd]的隨機響應函數y=g(x),其方差D可以表示為

(12)

式中,g0表示y的均值。

進一步,(12)式可以分解為

(13)

因此,全局指數定義為

(14)

基于(14)式,有以下特性

(15)

由于在不確定性量化中采用PC和深度學習方法,二者都相當于構建了數值模擬響應預測的一個隨機代理模型,計算非常快,因此靈敏度分析結合MCS可非常便捷地快速執行。考慮到混沌多項式方法將隨機輸出響應表示為一組正交多項式的加權組合,因此可解析地得到響應的方差(見1.1.1小節),從而非常方便地得到全局靈敏度指標的解析表達。相比于MCS,該方法可進一步減少計算資源消耗。

全局靈敏度指標可表示為

(16)

式中,DPC為總方差,可通過(7)式計算得到。

1.4 基于優質小樣本的模型參數修正

若步驟3中進行確認度量后,數值模擬模型不滿足要求,則基于當前參數得到數值模擬響應輸出的na個樣本,由于采用PC或深度學習進行不確定性傳播和量化,該過程不耗費任何函數調用,na可以設置得很大。從na個數值模擬輸出樣本中挑選出距離試驗數據ye最近的一定數目的樣本,將部分樣本稱作優質小樣本,設置一定的截斷比率r

(17)

式中,n′為優質小樣本的樣本數目。

截斷比率r需要在模型確認開始之前定義。它主要用于控制模型確認迭代過程的收斂速度,其大小取決于初始數值模擬數據與試驗數據的差異程度。一般來說,r取值越大,每次截取的優質小樣本的容量越大,則修正過程收斂越慢,從而需要更多的迭代步數;r取值越小,則修正過程收斂速度越快,但有可能導致最后修正精度不足。同時,收斂速度也與待修正參數的數量有關,待修正的參數越多,收斂速度越慢。基于作者的經驗,建議在一次修正過程中同時修正的參數不超過3個,截斷比率r取0.05～0.3。基于優質小樣本的模型參數修正方法其原理如圖3所示,通過反推找出每個優質小樣本相對應的輸入參數樣本,基于所有優質小樣本對應的輸入,結合模型參數的分布類型更新分布參數,如均勻分布的上下界、正態分布的均值和方差等。度量指標ε也是該方法的重要參數之一,ε的值取決于用戶對模型確認精度的需求,更高的精度需要更長的收斂過程。基于作者經驗,建議度量指標ε<0.05。

2 仿真驗證

本節將考慮湍流模型系數存在不確定性,將提出的數值模擬模型確認閉環流程應用于NACA0012翼型的CFD模型確認,提升CFD氣動預測數據與試驗數據的一致性。在本節中,模型確認問題中的截斷比率取r=0.1,確認度量采用MRE方法,度量指標取ε=0.01。

本算例選取SA湍流模型[22]進行CFD數值模擬,考慮其中的6個封閉系數存在不確定性并進行修正,認為其在變化范圍內服從均勻分布。利用ICEM軟件生成計算用網格,在空氣動力學分析中,使用Fluent 19.0作為CFD求解器,電腦配置為16 G內存,顯卡為NVIDIA 1660,CPU為Intel(R) Core(TM) i5-9400F。

考慮進行模型確認的流動條件為Ma=0.15和Re=6×106,且存在攻角α為0°,2°,4°,6°,8°,10°和11°多種工況。模型確認中進行確認度量和參數修正所用到的升力系數試驗數據來源于NASA網站的公開數據[23]。在給定模型封閉系數初始分布時,若封閉系數的初始區間設置較大,在這個區間對封閉系數取值進行CFD數值模擬,得到的輸出響應上下邊界范圍內一定可以包含試驗數據,也就是說當前流動情況下的模型系數最佳值一定包含在初始區間內。但很多情況下,給定的初始區間很可能較小且極有可能偏離最佳參數范圍。為了較為全面地驗證所提出的模型確認方法的有效性,考慮2種情況進行測試:

情況1封閉系數的初始范圍較大,封閉系數在其變化區間內服從均勻分布,其變化范圍見表1。

情況2封閉系數的初始范圍較小,假設湍流模型系數在默認值處服從±10%的均勻分布,默認值見表1。

針對以上2種情況所采取的模型確認策略也有所不同。對于情況1,由于封閉系數的初始不確定性極大,在模型確認過程中需要降低其不確定性,即維持其分布類型不變并減小封閉系數的變化區間;對于情況2,封閉系數的初始不確定性較小,但是極有可能處在不合適的區域,因此模型確認需要調整其范圍,即維持其分布類型不變并修正均值。

首先,對NACA0012基準翼型在名義流動條件下(α=10°,Ma=0.15,Re=6×106)進行了網格收斂測試,避免網格帶來的氣動預測不確定性,在測試中湍流模型選擇SA,模型系數取表2所示的默認值。網格收斂性測試的結果見表2,從中可以看出,當網格密度增加到一定數值時(遠場:150;翼型邊界:300),升力系數幾乎沒有變化。因此,在本小節的CFD模擬中,網格遠場的節點數被設定為150,而翼型邊界的節點數被設定為300。同時,對比該工況下的升力系數試驗數據(1.080 9),可以發現基于CFD仿真得到的升力系數存在7%左右的誤差,也證明了對CFD進行模型確認的必要性。

表2 網格收斂性測試

利用混沌多項式(PC)方法進行不確定性量化,構建以6個湍流模型封閉系數為輸入、升力系數為輸出的2階PC模型,并利用最小二次回歸方法計算PC系數。

2.1 情況1

基于湍流系數的分布類型和初始分布參數進行第一次UQ,獲得CFD數值模擬輸出響應的取值范圍,即升力系數的不確定性區間。UQ的結果如圖4所示,從中可以看出初始階段CFD預測的升力系數的不確定度(圖中藍色區域面積)較大,尤其在大攻角范圍(>10°)非常明顯,同時可以計算得到此時的最大EMR為0.43,明顯不滿足精度要求且遠超容忍值(ε=5%),因此需要開展參數修正。

圖4 初始階段CFD不確定性量化的結果(情況1)

為了降低參數修正的難度,提高收斂速度,首先采用1.3節介紹的方法對升力系數進行全局靈敏度分析,查找對升力系數影響最大的湍流模型系數,各個封閉系數的全局指數如表3所示。從表中可以看出,對于這7種工況,Cb1對升力系數的影響最大,且都遠大于其他幾個封閉系數。因此,僅對Cb1進行參數修正,此時其他系數固定在默認值處。

表3 湍流系數全局指數 %

CFD參數修正過程中的升力系數的不確定性區間見圖5中藍色部分。隨著修正過程的不斷迭代,對封閉系數不斷修正,相當于封閉系數取值的認知不確定性逐步減少,則升力系數的不確定度也隨之逐漸降低(即圖片中藍色區域逐漸縮小),最終CFD仿真的結果與試驗數據十分吻合,完成模型確認。

圖5 參數修正CFD升力系數預測區間和試驗數據(情況1)

Cb1的確認結果展示在表4中,從中可以發現經模型確認后Cb1的變化范圍中并非包含默認值(0.135 5),對Cb1的變化區間取中值,顯然與默認值略有不同。這也印證了文獻[1]中的觀點,按照湍流模型給定的默認系數,CFD仿真結果并非與試驗數據高度吻合。

表4 迭代過程中Cb1區間的變化(情況1)

2.2 情況2

首先進行初始狀態的UQ,結果展示于圖6a)。

圖6 參數修正CFD升力系數預測區間和試驗數據(情況2)

從圖中可以看出,初始狀態下的CFD預測的升力系數不確定性區間并沒有很好地覆蓋試驗數據。基于表3得知,對升力系數影響最大的都是Cb1,因此僅對Cb1進行修正,其他系數固定在默認值處。

CFD參數修正過程中的升力系數不確定性區間見圖6中藍色部分。隨著修正過程的不斷迭代,CFD不確定性區間逐漸向上偏移,最終能完全覆蓋試驗數據。

Cb1的確認結果展示于表5中,可以發現經確認后Cb1的取值與表4中的結果高度吻合,這也再一次驗證了所提出方法的有效性。

表5 迭代過程中Cb1區間的變化(情況2)

3 結 論

本文針對數值模擬模型確認,提出了一種考慮不確定性的數值模擬模型確認方法,并通過NACA0012翼型CFD數值模擬的模型確認應用,驗證了提出的數值模擬模型確認閉環流程的有效性。研究表明,該方法有如下優勢:

1) 建立一套涵括不確定性量化、全局靈敏度分析、參數修正、確認度量的模型確認閉環流程,可以合理有效地降低數值模擬中的不確定性;

2) 利用全局靈敏度分析篩選對輸出影響較大的模型參數,有效降低了模型確認的計算量;

3) 建立一種基于優質小樣本的模型參數修正策略,選取部分優質小樣本實現反向不確定性量化,從而快速有效地實現模型參數修正;

4) 提出的方法豐富了目前不確定性下的數值模擬模型確認理論和規范,為建立高可信度的學科或系統仿真分析模型提供了科學系統的方法,為后續開展飛行器穩健優化設計提供高保真的分析模型;

5) 提出的方法通過將模型看作一個黑箱問題,僅關注待修正參數和輸出的關系,極大地降低了模型確認的難度,并能夠推廣至工程應用中的各種模型確認問題。