趙輝,胡星志,張健,陳江濤,馬明生
中國空氣動力研究與發(fā)展中心 計算空氣動力研究所,綿陽 621000
近些年來,隨著網(wǎng)格生成、流場求解和高性能計算機(jī)等學(xué)科的快速發(fā)展,計算流體力學(xué)(Computational Fluid Dynamics,CFD)在工業(yè)領(lǐng)域發(fā)揮了越來越重要的作用。但是在復(fù)雜工程外形的數(shù)值模擬過程中,存在著多種來源的不確定性因素,包括幾何外形、網(wǎng)格分布、物理模型及數(shù)值計算方法等[1-2],這使得CFD給出的最終結(jié)果也存在著復(fù)雜的不確定性。在工程外形的評估或者優(yōu)化問題中,如果不考慮CFD結(jié)果的不確定性,往往不能準(zhǔn)確反映復(fù)雜物理條件下真實外形的氣動性能。因此CFD的不確定度量化成為近期研究的熱點。
對于工程中的湍流問題,通常都是通過湍流模型來?;字Z應(yīng)力,封閉雷諾平均方程。湍流模型是典型的認(rèn)知不確定性問題。由于對湍流物理本質(zhì)的理解不完備,導(dǎo)致了湍流模型存在明顯的不確定性。比如模型中的系數(shù)很多都是通過試探方法和經(jīng)驗給出,或者應(yīng)用漸進(jìn)性原則通過簡單流動的直接數(shù)值模擬或試驗結(jié)果來確定的[3]。但是這樣得到的湍流模型并不能保證對所有流動情況都適用。因此評估湍流模型系數(shù)的不確定性對數(shù)值模擬結(jié)果的影響對于CFD的確認(rèn)問題或者工程外形的設(shè)計優(yōu)化都十分重要。
在常用的不確定度量化方法中,蒙特卡羅(Monte Carlo,MC)方法[4]由于構(gòu)造簡單、易于實現(xiàn)被廣泛使用。但是MC方法收斂速度慢,需要大量的采樣點才能得到滿足精度要求的結(jié)果,因此限制了其在大維度多變量不確定性問題中的應(yīng)用。隨機(jī)譜方法,包括本文要討論的多項式混沌(Polynomial Chaos,PC)方法[5-6],是一種更有效率的不確定度量化方法。PC方法將整個數(shù)值模擬過程認(rèn)定為一個隨機(jī)過程,根據(jù)輸入?yún)?shù)的概率密度分布函數(shù),選擇合適的正交基函數(shù),然后將該隨機(jī)過程的輸出在基函數(shù)構(gòu)成的譜空間中展開。通過嵌入式或者非嵌入式方法,確定展開式中的自由度,從而得到輸出的完整的譜空間展開和各種統(tǒng)計信息,比如平均值、均方根以及每個不確定輸入變量的全局敏感度指標(biāo)。使用PC方法的不確定度量化已經(jīng)有了很多成功的應(yīng)用[7-14]。
本文使用非嵌入式PC方法研究了Spalart-Allmaras(SA)模型系數(shù)的不確定度對RAE2822翼型數(shù)值模擬的影響。首先,介紹了非嵌入式PC方法和使用的流場解算工具,研究從單變量入手,考察了卡門常數(shù)的不確定度對升阻力系數(shù)、壓力分布、摩擦系數(shù)分布和馬赫數(shù)的影響。然后,同時考慮了SA模型中9個參數(shù)的不確定度對數(shù)值模擬的影響。最后,給出了研究結(jié)論。
PC方法源于Wiener提出的各向同性混沌理論[15],用于湍流問題中隨機(jī)變量的譜空間展開。Ghanem和Spanos將該方法應(yīng)用到有限元計算中,用來分析固體力學(xué)中的不確定性問題[16]。Xiu和Karniadakis將適合高斯隨機(jī)過程的 Hermite PC推廣到 Wiener-Askey PC,適用于滿足更一般的概率密度分布的隨機(jī)過程[6]。PC方法基于不確定量的譜空間展開,將隨機(jī)量分解為確定性的和隨機(jī)的兩部分。下面首先從單輸入變量入手,介紹PC方法。
假定隨機(jī)的輸入變量ξ,滿足某種概率密度分布f(ξ)。對于任意的隨機(jī)輸出變量y(可以是升力、阻力系數(shù)等積分量,也可以是流場中某處的壓力系數(shù)、摩擦力系數(shù)或者馬赫數(shù),或者分離區(qū)長度、再附點位置等流場特征量),可以在輸入變量ξ的正交基函數(shù)序列構(gòu)成的譜空間中展開,即
式中:αj為確定性分量,也可稱為自由度;ψj為隨機(jī)變量的基函數(shù)。理論上式(1)的展開有無窮多項,但是在實際應(yīng)用中通常將展開在某階次處截斷,即
正交基函數(shù)的選擇由隨機(jī)輸入變量的概率密度函數(shù)決定,滿足:
即 {ψj,j≥0}是帶權(quán)函數(shù)為f(ξ)的正交多項式序列。對于均勻分布的隨機(jī)變量,最優(yōu)的基函數(shù)序列是Legendre多項式;而對于正態(tài)分布的隨機(jī)變量,最優(yōu)的基函數(shù)序列是Hermite多項式。這里談到的最優(yōu),是指系統(tǒng)輸出的統(tǒng)計信息隨著多項式階數(shù)增加而收斂的速度與理論一致。
根據(jù)與求解器耦合方式的不同,可以將PC分為嵌入式多項式混沌(Intrusive Polynomial Chaos,IPC)法和非嵌入多項式混沌(Non-Intrusive Polynomial Chaos,NIPC)法。在IPC中,控制方程中的不確定變量和參數(shù)都用其PC展開替代,經(jīng)過投影之后,得到與原控制方程相似的隨機(jī)控制方程,自變量為PC展開中的自由度。IPC需要對現(xiàn)有代碼做較大的改動,不僅工作量大,而且很多時候現(xiàn)有代碼都是封閉的,這限制了IPC在工程問題中的應(yīng)用。而NIPC無需修改控制方程,將現(xiàn)有的求解器視為一個黑匣子。首先在輸入變量的隨機(jī)空間里通過恰當(dāng)?shù)某闃臃椒ǐ@得一系列的采樣點,傳遞給求解器,進(jìn)行確定性的計算,得到各采樣點對應(yīng)的輸出響應(yīng)值,然后通過下文介紹的投影法或回歸法,得到完整的PC展開。本文主要討論NIPC方法。
在NIPC中,主要有兩種方法求解PC展開式中的自由度,投影法和回歸法。投影法利用基函數(shù)的正交特性,在式(2)的左右兩側(cè)同時乘以ψj,并在ξ的支撐集內(nèi)積分,利用基函數(shù)的正交性質(zhì),可以得到
式中:右側(cè)積分〈yψj〉可以通過蒙特卡羅積分或者高斯積分得到。當(dāng)輸入變量的維數(shù)較大時,投影法需要大量的采樣點才能滿足精度要求,即使使用稀疏網(wǎng)格的高斯積分,計算量也是巨大的。
在回歸法中,通過確定性的CFD計算得到采樣點 {ξ1,ξ2,…,ξN}對 應(yīng) 的 響 應(yīng) 值 {y1,y2,…,yN}。通常情況下,采樣點的數(shù)目大于PC展開中自由度的個數(shù)。因此形成一個超定方程組Ψα=Y(jié),其中Ψ是測量矩陣,Ψij=ψj(ξi)。通過最小二乘法可以得到自由度組成的向量α。本文中自由度計算都是采用回歸法得到。在單變量的不確定度量化中,使用高斯積分法進(jìn)行結(jié)果的對比。
采樣點的數(shù)目N跟PC展開中的自由度個數(shù)有關(guān),即
式中:n為輸入隨機(jī)變量的維數(shù);p為混沌多項式的階數(shù);np定義為過采樣率。Hosder等[17]在計算中發(fā)現(xiàn),np=2時,PC展開對輸出變量統(tǒng)計信息的近似較好。因此在本文的計算中np均取為2。
得到輸出的PC展開之后,可以根據(jù)基函數(shù)的正交特性直接得到輸出的均值μ、均方根σ等統(tǒng)計信息,即
對于多輸入變量的情況,如果各個變量之間是相互獨立的,多維基函數(shù)可以直接取為單維基函數(shù)的乘積形式,聯(lián)合概率密度函數(shù)為各變量概率密度函數(shù)的乘積。
在多輸入變量的情況下,每個變量的不確定度對系統(tǒng)輸出的不確定度的貢獻(xiàn)大小不一,很多時候需要量化每個變量的相對貢獻(xiàn)大小,分析輸入變量的全局敏感性。Sobol指標(biāo)[18]提供了工具來量化每個輸入變量對系統(tǒng)輸出方差的貢獻(xiàn)程度。Sobol指標(biāo)通過Sobol分解得到,是一種基于方差分解的全局敏感性分析方法,其優(yōu)點是考慮了隨機(jī)輸入?yún)?shù)在整個取值范圍內(nèi)對輸出響應(yīng)的貢獻(xiàn)以及隨機(jī)輸入?yún)?shù)的交互作用[17]。Sobol指標(biāo)的詳細(xì)推導(dǎo)過程和具體形式可以參考文獻(xiàn)[18-19]。
本文考慮的算例是RAE2822翼型[20],這是一個典型的超臨界翼型,也是測試CFD程序跨聲速流動模擬能力的經(jīng)典算例,國內(nèi)外研究機(jī)構(gòu)針對該翼型開展了大量的風(fēng)洞試驗和數(shù)值模擬研究。本文的計算狀態(tài)是:馬赫數(shù)Ma=0.729,雷諾數(shù)Rec=6.5×106,迎角α=2.31°。使用的網(wǎng)格如圖1所示。
計算使用的程序是課題組自行研發(fā)的MFlow[21]。該程序采用基于格心的有限體積法,能夠處理多種網(wǎng)格單元類型(六面體、四面體、三棱柱、金字塔以及使用幾何多重過程中融合產(chǎn)生的多面體)。單元內(nèi)使用線性重構(gòu)達(dá)到二階精度。梯度計算使用基于格點的格林高斯方法[22]。使用Venkatakrishnan限制器[23]來抑制間斷附近的振蕩。使用Roe格式來計算無黏通量。計算使用一階歐拉后差來推進(jìn)到收斂狀態(tài),使用局部時間步長加速收斂過程。Jacobian通量使用一階迎風(fēng)格式得到。
圖1 RAE2822翼型網(wǎng)格分布Fig.1 Grid distribution for RAE2822airfoil
計算使用SA一方程模型[24],這是在工程界應(yīng)用非常廣泛的湍流模型。在使用中,假定流動為全湍流,忽略了原始模型中的轉(zhuǎn)捩項(Trip Term)。因此模型中有9個常數(shù),分別為:cb1、σ、cb2、κ、cw2、cw3、cv1、ct3、ct4。SA 模型的輸運方程形式和模型系數(shù)是通過量綱分析、伽利略不變性和部分試驗(包括二維混合層、尾流和平板邊界層等)結(jié)果給出的[25]。模型參數(shù)的取值存在著一定的不確定度。比較典型的是卡門常數(shù)κ,κ描述了邊界層內(nèi)平均流向速度U+隨離開壁面距離y+變化的對數(shù)關(guān)系[26]。在大氣表層試驗中,κ的值最大為Sheppard測得的0.46,最小為Kansas的結(jié)果0.35[27-28]。也有學(xué)者指出κ跟雷諾數(shù)和壓力梯度有關(guān)[29-30]。佘振蘇給出了規(guī)范壁湍流中普適卡門常數(shù)(κ=0.45)的定義[31-32]。本文參考了文獻(xiàn)[33]中8個變量的取值范圍,如表1所示,假定每個參數(shù)在各自的支撐集內(nèi)為均勻分布。
表1 SA模型的標(biāo)準(zhǔn)參數(shù)和支撐Table 1 SA model closure coefficients and associated support
在進(jìn)行不確定度量化之前需要指出的是本文側(cè)重于研究SA模型系數(shù)的不確定度對該算例模擬結(jié)果的不確定度影響,沒有考慮網(wǎng)格分布、數(shù)值格式等方面的不確定度帶來的影響。
從單變量入手,首先研究了卡門常數(shù)的不確定度對RAE2822翼型模擬的影響。由于事先無法知曉系統(tǒng)輸出對卡門常數(shù)的響應(yīng)關(guān)系,所以假定PC展開的截斷階次分別為p=1,2,3,4。圖2和圖3給出了翼型升力、阻力系數(shù)CL和CD的平均值和均方差值結(jié)果。可以清楚地看到,隨著PC展開的階次增加,升力、阻力系數(shù)的統(tǒng)計信息變化很小,特別是平均值只有0.01%以內(nèi)的差別。圖4給出了所有采樣點的升力、阻力系數(shù)隨著卡門常數(shù)變化的曲線,據(jù)此可以推斷升力、阻力系數(shù)對卡門常數(shù)是近似線性的響應(yīng)關(guān)系。圖5給出了升力、阻力系數(shù)的PC展開中各項自由度隨多項式階次的變化,二階和二階以上基函數(shù)對應(yīng)的自由度比一階基函數(shù)對應(yīng)的自由度量值小兩個量級以上,進(jìn)一步證實了上述結(jié)論。
圖2 升力系數(shù)的平均值和均方根值Fig.2 Mean values and standard deviations of lift coefficients
圖3 阻力系數(shù)的平均值和均方根值Fig.3 Mean values and standard deviations of drag coefficients
圖4 所有采樣點的升力、阻力系數(shù)Fig.4 Lift and drag coefficients for all samples
圖5 升力、阻力系數(shù)的自由度Fig.5 Freedoms of lift and drag coefficients
得到輸出的PC展開后,可以預(yù)測任意一點的輸入對應(yīng)的輸出響應(yīng)。表2給出了κ為模型標(biāo)準(zhǔn)取值0.41時,PC預(yù)測值和CFD計算的結(jié)果,兩者的誤差幾乎可以忽略不計。說明建立的PC展開能夠準(zhǔn)確地描述輸出對輸入的響應(yīng)。
在單變量情況下,利用高斯積分法和回歸法計算PC展開的自由度,計算量相當(dāng)。圖6給出了兩種方法得到的升力、阻力系數(shù)PC展開中各項自由度,兩者結(jié)果也基本一致。
在上面的分析中,系統(tǒng)輸出關(guān)注的是升力、阻力系數(shù)等積分量,下面進(jìn)一步討論流場中某處的壓力系數(shù)、摩擦力系數(shù)或者馬赫數(shù)。圖7和圖8給出了壁面壓力分布的不確定度分析結(jié)果。其中包括壓力分布的試驗值、平均值、最大值、最小值和方差值。κ的不確定度主要影響翼型上表面,特別是激波附近的壓力分布。減小κ會使得激波位置后移。
圖9給出了壁面摩擦系數(shù)Cf預(yù)測的均值和極值。κ的不確定度影響著翼型上、下表面的摩擦系數(shù)預(yù)測。當(dāng)κ在其支撐集內(nèi)取最小值時,預(yù)測的壁面摩擦系數(shù)最?。划?dāng)κ取最大值時,預(yù)測的壁面摩擦系數(shù)最大。對于激波位置和壁面摩擦系數(shù)隨著κ變化的原因,可以從湍流渦黏系數(shù)υt的變化來解釋。在簡單平面剪切湍流的等雷諾應(yīng)力層,渦黏系數(shù)υt與κ成正比[3],即υt=κuτd,其中uτ為壁面摩擦速度,d為到壁面的距離,因此減小κ會得到較小的渦黏系數(shù)。圖10給出了κ取最大值時計算的渦黏系數(shù)與κ取最小值時的渦黏系數(shù)的差值云圖。從圖可以清楚地看到,在流場的大部分區(qū)域,減小κ會減小渦黏系數(shù),從而導(dǎo)致預(yù)測的壁面摩擦系數(shù)減小,同時邊界層內(nèi)流體微團(tuán)動能損失較少,抵御逆壓梯度的能力更強(qiáng),預(yù)測的激波位置因此后移。
表2 PC展開與CFD得到的升力、阻力系數(shù)比較Table 2 Comparison of lift and drag coefficients obtained with PC and CFD
圖6 升力、阻力系數(shù)的自由度比較(回歸法和高斯積分法)Fig.6 Comparison of freedoms of lift and drag coefficients(regression method and Gauss integration method)
圖7 單變量下的壓力分布Fig.7 Pressure distribution under single variable
圖8 壓力分布的方差Fig.8 Variance of pressure distribution
圖11是馬赫數(shù)方差的云圖,κ的不確定度主要影響上表面激波位置的預(yù)測,從而也會影響激波附近馬赫數(shù)的預(yù)測。
圖9 單變量下壁面摩擦系數(shù)的平均值和極值Fig.9 Mean and extreme values of skin friction coefficient under single variable
圖10 渦黏系數(shù)差值云圖(υt|κ_max-υt|κ_min)
圖11 馬赫數(shù)方差云圖Fig.11 Contours of variance of Mach number
通過上面的分析可以知道κ的變化對激波位置和壁面摩擦預(yù)測的影響,也得到了系統(tǒng)輸出對κ的響應(yīng)。下面同時考慮SA模型中9個參數(shù)的不確定度對數(shù)值模擬的影響。計算中PC展開的階次為p=2,考慮了任意兩個參數(shù)的聯(lián)合影響。在9維空間中,采樣點通過拉丁超立方抽樣得到。
表3給出了升力、阻力系數(shù)預(yù)測中9個參數(shù)各自的Sobol指標(biāo),按照貢獻(xiàn)大小排序。與文獻(xiàn)[33]的結(jié)果存在明顯的差異,主要表現(xiàn)為在本文的計算中ct3、ct4這兩個參數(shù)對升力、阻力系數(shù)方差的貢獻(xiàn)很小。這是因為本文使用的SA模型忽略了原始模型中的轉(zhuǎn)捩項,由此很多學(xué)者指出包含ct3、ct4的ft2也是可以忽略的。同時cb2、cw3這兩個參數(shù)的貢獻(xiàn)也是可以忽略不計的。cb1、σ、κ、cw2、cv1這5個參數(shù)對于升力、阻力系數(shù)的不確定度影響都相對比較大,κ是其中貢獻(xiàn)最大的參數(shù),其他4個參數(shù)的貢獻(xiàn)大小排序并不完全一致。
表3 升力、阻力系數(shù)預(yù)測中9個參數(shù)的Sobol指標(biāo)Table 3 Sobol indices of 9closure coefficients for lift and drag coefficients
通過升力、阻力系數(shù)的響應(yīng),可以得到其平均值、均方差等統(tǒng)計信息,如表4所示。其中升力系數(shù) 的 最 大 值 取 在 {cb1_max、σ_max,κ_min,cw2_max,cv1_max},最 小 值 取 在 {cb1_min,σ_min,κ_max,cw2_min,cv1_min}。阻力 系 數(shù) 的 最 大 值 取 在 {cb1_max,σ_max,κ_max,cw2_max,cv1_min},最 小 值 取 在 {cb1_min,σ_min,κ_min,cw2_min,cv1_max}。ct3、ct4、cb2、cw3這 4 個 參 數(shù) 在支撐集內(nèi)變化對升力、阻力系數(shù)影響很小。
表4 升力、阻力系數(shù)的均值和均方差值Table 4 Mean value and standard deviation of lift and drag coefficients
得到輸出的PC展開后,可以預(yù)測任意一點的輸入對應(yīng)的輸出。表5給出了模型參數(shù)為標(biāo)準(zhǔn)取值時,PC預(yù)測值和CFD計算的對比,其中,CDp為壓差阻力,CDf為摩擦阻力,兩者誤差很小。說明建立的PC展開能夠準(zhǔn)確地描述輸出對多變量輸入的響應(yīng)。
表5 PC展開與CFD得到的升力、阻力系數(shù)比較Table 5 Comparison of lift and drag coefficients obtained through PC and CFD
通過確定的CFD計算可以得到壁面壓力分布響應(yīng)。圖12給出了壓力的試驗值、平均值和極值。在機(jī)翼的上表面,特別是激波附近,壓力的不確定度表現(xiàn)顯著。上表面壓力極小值取值在{cb1_max,σ_max,κ_min,cw2_max,cv1_max},極 大 值 取 值 在{cb1_min,σ_min,κ_max,cw2_min,cv1_min},與 升 力 系 數(shù) 的情況類似。圖13給出了9個輸入?yún)?shù)對壓力方差的貢獻(xiàn)度。與升力、阻力系數(shù)貢獻(xiàn)結(jié)論一致,ct3、ct4、cb2、cw3這4個參數(shù)對壓力的不確定度貢獻(xiàn)相對較小,其他5個參數(shù),特別是σ、κ貢獻(xiàn)較大。
圖14給出了壁面摩擦系數(shù)的平均值和極值。與壓力分布不同的是,ct3、ct4、cb2、cw3這4個參數(shù)在機(jī)翼的前緣對摩擦系數(shù)的不確定度也有著不可忽視的貢獻(xiàn),如圖15所示。
Fig.12 Pressure distribution under multi-variables
圖13 壓力分布的Sobol指標(biāo)分布Fig.13 Sobol index distributions for pressure distribution
圖14 多變量下壁面摩擦系數(shù)的平均值和極值Fig.14 Mean and extreme values of skin friction coefficient under multi-variables
圖15 壁面摩擦系數(shù)分布的Sobol指標(biāo)分布Fig.15 Sobol index distributions for skin friction coefficient distribution
SA模型中的9個參數(shù)綜合影響了湍流渦黏系數(shù)的生成、擴(kuò)散和耗散。比如參數(shù)cb1與渦黏系數(shù)的生成項和耗散項(部分)成正比,該參數(shù)的改變會同時增大或者減小生成和耗散,對渦黏系數(shù)的綜合作用尚且未知。從理論上很難得到參數(shù)變化對渦黏系數(shù)影響的規(guī)律,而且很可能在流場不同區(qū)域產(chǎn)生相反的結(jié)論。因此本文只對κ、ct3、ct4這3個參數(shù)引起的不確定度進(jìn)行簡單分析,其他參數(shù)的影響需要繼續(xù)深入研究。
參數(shù)ct3、ct4對模型的影響體現(xiàn)在生成衰減項ft2里。本文使用的SA模型忽略了原始模型中的轉(zhuǎn)捩項,由此ft2只在翼型前緣附近渦黏系數(shù)很小的區(qū)域發(fā)揮作用,在其他區(qū)域的邊界層內(nèi)基本沒有作用。因此ct3、ct4的不確定度對模擬的影響主要體現(xiàn)在前緣附近,特別是對摩擦系數(shù)的預(yù)測。
圖16給出了馬赫數(shù)的平均值和方差的云圖,在激波附近馬赫數(shù)表現(xiàn)出較大的不確定度,與設(shè)想的一致。
本文使用非嵌入式多項式混沌方法研究了湍流模型系數(shù)的不確定度對RAE2822跨聲速翼型繞流模擬的影響。
Fig.16 Contours of mean values and variances of Mach number
1)從單輸入變量入手,研究了卡門常數(shù)κ的不確定度對數(shù)值模擬的影響。研究發(fā)現(xiàn)升力、阻力系數(shù)對κ是近似線性的響應(yīng)關(guān)系。在單變量情況下,利用高斯積分法和回歸法計算PC展開的自由度,兩者結(jié)果也基本一致。κ的不確定度主要影響翼型上表面,特別是激波附近的壓力分布,同時還影響翼型上、下表面的摩擦預(yù)測結(jié)果和激波附近馬赫數(shù)的預(yù)測結(jié)果。κ減小時,預(yù)測的激波位置后移,摩擦系數(shù)預(yù)測值變小。
2)同時考慮模型中9個參數(shù)的不確定性帶來的影響。通過升力、阻力系數(shù)預(yù)測中9個參數(shù)各自的Sobol指標(biāo)可以發(fā)現(xiàn),ct3、ct4、cb2、cw3這4個參數(shù)對升力、阻力系數(shù)方差的貢獻(xiàn)很小,可以忽略不計。其他5個參數(shù)對于結(jié)果的不確定度影響都相對比較大,κ是其中貢獻(xiàn)最大的參數(shù)。建立的PC能夠給出壁面壓力分布、摩擦系數(shù)分布和馬赫數(shù)分布的不確定度。
3)通過PC展開和CFD計算的對比證實了建立的PC展開能夠準(zhǔn)確地描述輸出對輸入的響應(yīng)。這為量化計算結(jié)果的不確定度提供了切實可行、有較高效率的解決方案。需要指出的是本文的計算只考慮了RAE2822跨聲速翼型模擬的單一計算狀態(tài),影響規(guī)律是否可以推及其他工況和算例需要進(jìn)一步檢驗。