趙輝，胡星志，張健，陳江濤，馬明生

中國(guó)空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所，綿陽(yáng) 621000

近些年來(lái)，隨著網(wǎng)格生成、流場(chǎng)求解和高性能計(jì)算機(jī)等學(xué)科的快速發(fā)展，計(jì)算流體力學(xué)（Computational Fluid Dynamics，CFD）在工業(yè)領(lǐng)域發(fā)揮了越來(lái)越重要的作用。但是在復(fù)雜工程外形的數(shù)值模擬過(guò)程中，存在著多種來(lái)源的不確定性因素，包括幾何外形、網(wǎng)格分布、物理模型及數(shù)值計(jì)算方法等［1－2］，這使得CFD給出的最終結(jié)果也存在著復(fù)雜的不確定性。在工程外形的評(píng)估或者優(yōu)化問(wèn)題中，如果不考慮CFD結(jié)果的不確定性，往往不能準(zhǔn)確反映復(fù)雜物理?xiàng)l件下真實(shí)外形的氣動(dòng)性能。因此CFD的不確定度量化成為近期研究的熱點(diǎn)。

對(duì)于工程中的湍流問(wèn)題，通常都是通過(guò)湍流模型來(lái)?；字Z應(yīng)力，封閉雷諾平均方程。湍流模型是典型的認(rèn)知不確定性問(wèn)題。由于對(duì)湍流物理本質(zhì)的理解不完備，導(dǎo)致了湍流模型存在明顯的不確定性。比如模型中的系數(shù)很多都是通過(guò)試探方法和經(jīng)驗(yàn)給出，或者應(yīng)用漸進(jìn)性原則通過(guò)簡(jiǎn)單流動(dòng)的直接數(shù)值模擬或試驗(yàn)結(jié)果來(lái)確定的［3］。但是這樣得到的湍流模型并不能保證對(duì)所有流動(dòng)情況都適用。因此評(píng)估湍流模型系數(shù)的不確定性對(duì)數(shù)值模擬結(jié)果的影響對(duì)于CFD的確認(rèn)問(wèn)題或者工程外形的設(shè)計(jì)優(yōu)化都十分重要。

在常用的不確定度量化方法中，蒙特卡羅（Monte Carlo，MC）方法［4］由于構(gòu)造簡(jiǎn)單、易于實(shí)現(xiàn)被廣泛使用。但是MC方法收斂速度慢，需要大量的采樣點(diǎn)才能得到滿足精度要求的結(jié)果，因此限制了其在大維度多變量不確定性問(wèn)題中的應(yīng)用。隨機(jī)譜方法，包括本文要討論的多項(xiàng)式混沌（Polynomial Chaos，PC）方法［5－6］，是一種更有效率的不確定度量化方法。PC方法將整個(gè)數(shù)值模擬過(guò)程認(rèn)定為一個(gè)隨機(jī)過(guò)程，根據(jù)輸入?yún)?shù)的概率密度分布函數(shù)，選擇合適的正交基函數(shù)，然后將該隨機(jī)過(guò)程的輸出在基函數(shù)構(gòu)成的譜空間中展開(kāi)。通過(guò)嵌入式或者非嵌入式方法，確定展開(kāi)式中的自由度，從而得到輸出的完整的譜空間展開(kāi)和各種統(tǒng)計(jì)信息，比如平均值、均方根以及每個(gè)不確定輸入變量的全局敏感度指標(biāo)。使用PC方法的不確定度量化已經(jīng)有了很多成功的應(yīng)用［7－14］。

本文使用非嵌入式PC方法研究了Spalart-Allmaras（SA）模型系數(shù)的不確定度對(duì)RAE2822翼型數(shù)值模擬的影響。首先，介紹了非嵌入式PC方法和使用的流場(chǎng)解算工具，研究從單變量入手，考察了卡門常數(shù)的不確定度對(duì)升阻力系數(shù)、壓力分布、摩擦系數(shù)分布和馬赫數(shù)的影響。然后，同時(shí)考慮了SA模型中9個(gè)參數(shù)的不確定度對(duì)數(shù)值模擬的影響。最后，給出了研究結(jié)論。

1 PC方法展開(kāi)

PC方法源于Wiener提出的各向同性混沌理論［15］，用于湍流問(wèn)題中隨機(jī)變量的譜空間展開(kāi)。Ghanem和Spanos將該方法應(yīng)用到有限元計(jì)算中，用來(lái)分析固體力學(xué)中的不確定性問(wèn)題［16］。Xiu和Karniadakis將適合高斯隨機(jī)過(guò)程的 Hermite PC推廣到 Wiener-Askey PC，適用于滿足更一般的概率密度分布的隨機(jī)過(guò)程［6］。PC方法基于不確定量的譜空間展開(kāi)，將隨機(jī)量分解為確定性的和隨機(jī)的兩部分。下面首先從單輸入變量入手，介紹PC方法。

假定隨機(jī)的輸入變量ξ，滿足某種概率密度分布f（ξ）。對(duì)于任意的隨機(jī)輸出變量y（可以是升力、阻力系數(shù)等積分量，也可以是流場(chǎng)中某處的壓力系數(shù)、摩擦力系數(shù)或者馬赫數(shù)，或者分離區(qū)長(zhǎng)度、再附點(diǎn)位置等流場(chǎng)特征量），可以在輸入變量ξ的正交基函數(shù)序列構(gòu)成的譜空間中展開(kāi)，即

式中：αj為確定性分量，也可稱為自由度；ψj為隨機(jī)變量的基函數(shù)。理論上式（1）的展開(kāi)有無(wú)窮多項(xiàng)，但是在實(shí)際應(yīng)用中通常將展開(kāi)在某階次處截?cái)?，?/p>

正交基函數(shù)的選擇由隨機(jī)輸入變量的概率密度函數(shù)決定，滿足：

即 ｛ψj，j≥0｝是帶權(quán)函數(shù)為f（ξ）的正交多項(xiàng)式序列。對(duì)于均勻分布的隨機(jī)變量，最優(yōu)的基函數(shù)序列是Legendre多項(xiàng)式；而對(duì)于正態(tài)分布的隨機(jī)變量，最優(yōu)的基函數(shù)序列是Hermite多項(xiàng)式。這里談到的最優(yōu)，是指系統(tǒng)輸出的統(tǒng)計(jì)信息隨著多項(xiàng)式階數(shù)增加而收斂的速度與理論一致。

根據(jù)與求解器耦合方式的不同，可以將PC分為嵌入式多項(xiàng)式混沌（Intrusive Polynomial Chaos，IPC）法和非嵌入多項(xiàng)式混沌（Non-Intrusive Polynomial Chaos，NIPC）法。在IPC中，控制方程中的不確定變量和參數(shù)都用其PC展開(kāi)替代，經(jīng)過(guò)投影之后，得到與原控制方程相似的隨機(jī)控制方程，自變量為PC展開(kāi)中的自由度。IPC需要對(duì)現(xiàn)有代碼做較大的改動(dòng)，不僅工作量大，而且很多時(shí)候現(xiàn)有代碼都是封閉的，這限制了IPC在工程問(wèn)題中的應(yīng)用。而NIPC無(wú)需修改控制方程，將現(xiàn)有的求解器視為一個(gè)黑匣子。首先在輸入變量的隨機(jī)空間里通過(guò)恰當(dāng)?shù)某闃臃椒ǐ@得一系列的采樣點(diǎn)，傳遞給求解器，進(jìn)行確定性的計(jì)算，得到各采樣點(diǎn)對(duì)應(yīng)的輸出響應(yīng)值，然后通過(guò)下文介紹的投影法或回歸法，得到完整的PC展開(kāi)。本文主要討論NIPC方法。

在NIPC中，主要有兩種方法求解PC展開(kāi)式中的自由度，投影法和回歸法。投影法利用基函數(shù)的正交特性，在式（2）的左右兩側(cè)同時(shí)乘以ψj，并在ξ的支撐集內(nèi)積分，利用基函數(shù)的正交性質(zhì)，可以得到

式中：右側(cè)積分〈yψj〉可以通過(guò)蒙特卡羅積分或者高斯積分得到。當(dāng)輸入變量的維數(shù)較大時(shí)，投影法需要大量的采樣點(diǎn)才能滿足精度要求，即使使用稀疏網(wǎng)格的高斯積分，計(jì)算量也是巨大的。

在回歸法中，通過(guò)確定性的CFD計(jì)算得到采樣點(diǎn) ｛ξ1，ξ2，…，ξN｝對(duì) 應(yīng) 的 響 應(yīng) 值 ｛y1，y2，…，yN｝。通常情況下，采樣點(diǎn)的數(shù)目大于PC展開(kāi)中自由度的個(gè)數(shù)。因此形成一個(gè)超定方程組Ψα＝Y(jié)，其中Ψ是測(cè)量矩陣，Ψij＝ψj（ξi）。通過(guò)最小二乘法可以得到自由度組成的向量α。本文中自由度計(jì)算都是采用回歸法得到。在單變量的不確定度量化中，使用高斯積分法進(jìn)行結(jié)果的對(duì)比。

采樣點(diǎn)的數(shù)目N跟PC展開(kāi)中的自由度個(gè)數(shù)有關(guān)，即

式中：n為輸入隨機(jī)變量的維數(shù)；p為混沌多項(xiàng)式的階數(shù)；np定義為過(guò)采樣率。Hosder等［17］在計(jì)算中發(fā)現(xiàn)，np＝2時(shí)，PC展開(kāi)對(duì)輸出變量統(tǒng)計(jì)信息的近似較好。因此在本文的計(jì)算中np均取為2。

得到輸出的PC展開(kāi)之后，可以根據(jù)基函數(shù)的正交特性直接得到輸出的均值μ、均方根σ等統(tǒng)計(jì)信息，即

對(duì)于多輸入變量的情況，如果各個(gè)變量之間是相互獨(dú)立的，多維基函數(shù)可以直接取為單維基函數(shù)的乘積形式，聯(lián)合概率密度函數(shù)為各變量概率密度函數(shù)的乘積。

在多輸入變量的情況下，每個(gè)變量的不確定度對(duì)系統(tǒng)輸出的不確定度的貢獻(xiàn)大小不一，很多時(shí)候需要量化每個(gè)變量的相對(duì)貢獻(xiàn)大小，分析輸入變量的全局敏感性。Sobol指標(biāo)［18］提供了工具來(lái)量化每個(gè)輸入變量對(duì)系統(tǒng)輸出方差的貢獻(xiàn)程度。Sobol指標(biāo)通過(guò)Sobol分解得到，是一種基于方差分解的全局敏感性分析方法，其優(yōu)點(diǎn)是考慮了隨機(jī)輸入?yún)?shù)在整個(gè)取值范圍內(nèi)對(duì)輸出響應(yīng)的貢獻(xiàn)以及隨機(jī)輸入?yún)?shù)的交互作用［17］。Sobol指標(biāo)的詳細(xì)推導(dǎo)過(guò)程和具體形式可以參考文獻(xiàn)［18-19］。

2 計(jì)算外形和方法

本文考慮的算例是RAE2822翼型［20］，這是一個(gè)典型的超臨界翼型，也是測(cè)試CFD程序跨聲速流動(dòng)模擬能力的經(jīng)典算例，國(guó)內(nèi)外研究機(jī)構(gòu)針對(duì)該翼型開(kāi)展了大量的風(fēng)洞試驗(yàn)和數(shù)值模擬研究。本文的計(jì)算狀態(tài)是：馬赫數(shù)Ma＝0.729，雷諾數(shù)Rec＝6.5×106，迎角α＝2.31°。使用的網(wǎng)格如圖1所示。

計(jì)算使用的程序是課題組自行研發(fā)的MFlow［21］。該程序采用基于格心的有限體積法，能夠處理多種網(wǎng)格單元類型（六面體、四面體、三棱柱、金字塔以及使用幾何多重過(guò)程中融合產(chǎn)生的多面體）。單元內(nèi)使用線性重構(gòu)達(dá)到二階精度。梯度計(jì)算使用基于格點(diǎn)的格林高斯方法［22］。使用Venkatakrishnan限制器［23］來(lái)抑制間斷附近的振蕩。使用Roe格式來(lái)計(jì)算無(wú)黏通量。計(jì)算使用一階歐拉后差來(lái)推進(jìn)到收斂狀態(tài)，使用局部時(shí)間步長(zhǎng)加速收斂過(guò)程。Jacobian通量使用一階迎風(fēng)格式得到。

圖1 RAE2822翼型網(wǎng)格分布Fig.1 Grid distribution for RAE2822airfoil

計(jì)算使用SA一方程模型［24］，這是在工程界應(yīng)用非常廣泛的湍流模型。在使用中，假定流動(dòng)為全湍流，忽略了原始模型中的轉(zhuǎn)捩項(xiàng)（Trip Term）。因此模型中有9個(gè)常數(shù)，分別為：cb1、σ、cb2、κ、cw2、cw3、cv1、ct3、ct4。SA 模型的輸運(yùn)方程形式和模型系數(shù)是通過(guò)量綱分析、伽利略不變性和部分試驗(yàn)（包括二維混合層、尾流和平板邊界層等）結(jié)果給出的［25］。模型參數(shù)的取值存在著一定的不確定度。比較典型的是卡門常數(shù)κ，κ描述了邊界層內(nèi)平均流向速度U＋隨離開(kāi)壁面距離y＋變化的對(duì)數(shù)關(guān)系［26］。在大氣表層試驗(yàn)中，κ的值最大為Sheppard測(cè)得的0.46，最小為Kansas的結(jié)果0.35［27－28］。也有學(xué)者指出κ跟雷諾數(shù)和壓力梯度有關(guān)［29－30］。佘振蘇給出了規(guī)范壁湍流中普適卡門常數(shù)（κ＝0.45）的定義［31－32］。本文參考了文獻(xiàn)［33］中8個(gè)變量的取值范圍，如表1所示，假定每個(gè)參數(shù)在各自的支撐集內(nèi)為均勻分布。

表1 SA模型的標(biāo)準(zhǔn)參數(shù)和支撐Table 1 SA model closure coefficients and associated support

在進(jìn)行不確定度量化之前需要指出的是本文側(cè)重于研究SA模型系數(shù)的不確定度對(duì)該算例模擬結(jié)果的不確定度影響，沒(méi)有考慮網(wǎng)格分布、數(shù)值格式等方面的不確定度帶來(lái)的影響。

3 不確定度分析

3.1 單變量

從單變量入手，首先研究了卡門常數(shù)的不確定度對(duì)RAE2822翼型模擬的影響。由于事先無(wú)法知曉系統(tǒng)輸出對(duì)卡門常數(shù)的響應(yīng)關(guān)系，所以假定PC展開(kāi)的截?cái)嚯A次分別為p＝1，2，3，4。圖2和圖3給出了翼型升力、阻力系數(shù)CL和CD的平均值和均方差值結(jié)果?？梢郧宄乜吹剑S著PC展開(kāi)的階次增加，升力、阻力系數(shù)的統(tǒng)計(jì)信息變化很小，特別是平均值只有0.01%以內(nèi)的差別。圖4給出了所有采樣點(diǎn)的升力、阻力系數(shù)隨著卡門常數(shù)變化的曲線，據(jù)此可以推斷升力、阻力系數(shù)對(duì)卡門常數(shù)是近似線性的響應(yīng)關(guān)系。圖5給出了升力、阻力系數(shù)的PC展開(kāi)中各項(xiàng)自由度隨多項(xiàng)式階次的變化，二階和二階以上基函數(shù)對(duì)應(yīng)的自由度比一階基函數(shù)對(duì)應(yīng)的自由度量值小兩個(gè)量級(jí)以上，進(jìn)一步證實(shí)了上述結(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 所有采樣點(diǎn)的升力、阻力系數(shù)Fig.4 Lift and drag coefficients for all samples

圖5 升力、阻力系數(shù)的自由度Fig.5 Freedoms of lift and drag coefficients

得到輸出的PC展開(kāi)后，可以預(yù)測(cè)任意一點(diǎn)的輸入對(duì)應(yīng)的輸出響應(yīng)。表2給出了κ為模型標(biāo)準(zhǔn)取值0.41時(shí)，PC預(yù)測(cè)值和CFD計(jì)算的結(jié)果，兩者的誤差幾乎可以忽略不計(jì)。說(shuō)明建立的PC展開(kāi)能夠準(zhǔn)確地描述輸出對(duì)輸入的響應(yīng)。

在單變量情況下，利用高斯積分法和回歸法計(jì)算PC展開(kāi)的自由度，計(jì)算量相當(dāng)。圖6給出了兩種方法得到的升力、阻力系數(shù)PC展開(kāi)中各項(xiàng)自由度，兩者結(jié)果也基本一致。

在上面的分析中，系統(tǒng)輸出關(guān)注的是升力、阻力系數(shù)等積分量，下面進(jìn)一步討論流場(chǎng)中某處的壓力系數(shù)、摩擦力系數(shù)或者馬赫數(shù)。圖7和圖8給出了壁面壓力分布的不確定度分析結(jié)果。其中包括壓力分布的試驗(yàn)值、平均值、最大值、最小值和方差值。κ的不確定度主要影響翼型上表面，特別是激波附近的壓力分布。減小κ會(huì)使得激波位置后移。

圖9給出了壁面摩擦系數(shù)Cf預(yù)測(cè)的均值和極值。κ的不確定度影響著翼型上、下表面的摩擦系數(shù)預(yù)測(cè)。當(dāng)κ在其支撐集內(nèi)取最小值時(shí)，預(yù)測(cè)的壁面摩擦系數(shù)最?。划?dāng)κ取最大值時(shí)，預(yù)測(cè)的壁面摩擦系數(shù)最大。對(duì)于激波位置和壁面摩擦系數(shù)隨著κ變化的原因，可以從湍流渦黏系數(shù)υt的變化來(lái)解釋。在簡(jiǎn)單平面剪切湍流的等雷諾應(yīng)力層，渦黏系數(shù)υt與κ成正比［3］，即υt＝κuτd，其中uτ為壁面摩擦速度，d為到壁面的距離，因此減小κ會(huì)得到較小的渦黏系數(shù)。圖10給出了κ取最大值時(shí)計(jì)算的渦黏系數(shù)與κ取最小值時(shí)的渦黏系數(shù)的差值云圖。從圖可以清楚地看到，在流場(chǎng)的大部分區(qū)域，減小κ會(huì)減小渦黏系數(shù)，從而導(dǎo)致預(yù)測(cè)的壁面摩擦系數(shù)減小，同時(shí)邊界層內(nèi)流體微團(tuán)動(dòng)能損失較少，抵御逆壓梯度的能力更強(qiáng)，預(yù)測(cè)的激波位置因此后移。

表2 PC展開(kāi)與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ù)測(cè)，從而也會(huì)影響激波附近馬赫數(shù)的預(yù)測(cè)。

圖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

3.2 多變量

通過(guò)上面的分析可以知道κ的變化對(duì)激波位置和壁面摩擦預(yù)測(cè)的影響，也得到了系統(tǒng)輸出對(duì)κ的響應(yīng)。下面同時(shí)考慮SA模型中9個(gè)參數(shù)的不確定度對(duì)數(shù)值模擬的影響。計(jì)算中PC展開(kāi)的階次為p＝2，考慮了任意兩個(gè)參數(shù)的聯(lián)合影響。在9維空間中，采樣點(diǎn)通過(guò)拉丁超立方抽樣得到。

表3給出了升力、阻力系數(shù)預(yù)測(cè)中9個(gè)參數(shù)各自的Sobol指標(biāo)，按照貢獻(xiàn)大小排序。與文獻(xiàn)［33］的結(jié)果存在明顯的差異，主要表現(xiàn)為在本文的計(jì)算中ct3、ct4這兩個(gè)參數(shù)對(duì)升力、阻力系數(shù)方差的貢獻(xiàn)很小。這是因?yàn)楸疚氖褂玫腟A模型忽略了原始模型中的轉(zhuǎn)捩項(xiàng)，由此很多學(xué)者指出包含ct3、ct4的ft2也是可以忽略的。同時(shí)cb2、cw3這兩個(gè)參數(shù)的貢獻(xiàn)也是可以忽略不計(jì)的。cb1、σ、κ、cw2、cv1這5個(gè)參數(shù)對(duì)于升力、阻力系數(shù)的不確定度影響都相對(duì)比較大，κ是其中貢獻(xiàn)最大的參數(shù)，其他4個(gè)參數(shù)的貢獻(xiàn)大小排序并不完全一致。

表3 升力、阻力系數(shù)預(yù)測(cè)中9個(gè)參數(shù)的Sobol指標(biāo)Table 3 Sobol indices of 9closure coefficients for lift and drag coefficients

通過(guò)升力、阻力系數(shù)的響應(yīng)，可以得到其平均值、均方差等統(tǒng)計(jì)信息，如表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 個(gè) 參 數(shù) 在支撐集內(nèi)變化對(duì)升力、阻力系數(shù)影響很小。

表4 升力、阻力系數(shù)的均值和均方差值Table 4 Mean value and standard deviation of lift and drag coefficients

得到輸出的PC展開(kāi)后，可以預(yù)測(cè)任意一點(diǎn)的輸入對(duì)應(yīng)的輸出。表5給出了模型參數(shù)為標(biāo)準(zhǔn)取值時(shí)，PC預(yù)測(cè)值和CFD計(jì)算的對(duì)比，其中，CDp為壓差阻力，CDf為摩擦阻力，兩者誤差很小。說(shuō)明建立的PC展開(kāi)能夠準(zhǔn)確地描述輸出對(duì)多變量輸入的響應(yīng)。

表5 PC展開(kāi)與CFD得到的升力、阻力系數(shù)比較Table 5 Comparison of lift and drag coefficients obtained through PC and CFD

通過(guò)確定的CFD計(jì)算可以得到壁面壓力分布響應(yīng)。圖12給出了壓力的試驗(yàn)值、平均值和極值。在機(jī)翼的上表面，特別是激波附近，壓力的不確定度表現(xiàn)顯著。上表面壓力極小值取值在｛cb1＿max，σ＿max，κ＿min，cw2＿max，cv1＿max｝，極 大 值 取 值 在｛cb1＿min，σ＿min，κ＿max，cw2＿min，cv1＿min｝，與 升 力 系 數(shù) 的情況類似。圖13給出了9個(gè)輸入?yún)?shù)對(duì)壓力方差的貢獻(xiàn)度。與升力、阻力系數(shù)貢獻(xiàn)結(jié)論一致，ct3、ct4、cb2、cw3這4個(gè)參數(shù)對(duì)壓力的不確定度貢獻(xiàn)相對(duì)較小，其他5個(gè)參數(shù)，特別是σ、κ貢獻(xiàn)較大。

圖14給出了壁面摩擦系數(shù)的平均值和極值。與壓力分布不同的是，ct3、ct4、cb2、cw3這4個(gè)參數(shù)在機(jī)翼的前緣對(duì)摩擦系數(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個(gè)參數(shù)綜合影響了湍流渦黏系數(shù)的生成、擴(kuò)散和耗散。比如參數(shù)cb1與渦黏系數(shù)的生成項(xiàng)和耗散項(xiàng)（部分）成正比，該參數(shù)的改變會(huì)同時(shí)增大或者減小生成和耗散，對(duì)渦黏系數(shù)的綜合作用尚且未知。從理論上很難得到參數(shù)變化對(duì)渦黏系數(shù)影響的規(guī)律，而且很可能在流場(chǎng)不同區(qū)域產(chǎn)生相反的結(jié)論。因此本文只對(duì)κ、ct3、ct4這3個(gè)參數(shù)引起的不確定度進(jìn)行簡(jiǎn)單分析，其他參數(shù)的影響需要繼續(xù)深入研究。

參數(shù)ct3、ct4對(duì)模型的影響體現(xiàn)在生成衰減項(xiàng)ft2里。本文使用的SA模型忽略了原始模型中的轉(zhuǎn)捩項(xiàng)，由此ft2只在翼型前緣附近渦黏系數(shù)很小的區(qū)域發(fā)揮作用，在其他區(qū)域的邊界層內(nèi)基本沒(méi)有作用。因此ct3、ct4的不確定度對(duì)模擬的影響主要體現(xiàn)在前緣附近，特別是對(duì)摩擦系數(shù)的預(yù)測(cè)。

圖16給出了馬赫數(shù)的平均值和方差的云圖，在激波附近馬赫數(shù)表現(xiàn)出較大的不確定度，與設(shè)想的一致。

4 結(jié) 論

本文使用非嵌入式多項(xiàng)式混沌方法研究了湍流模型系數(shù)的不確定度對(duì)RAE2822跨聲速翼型繞流模擬的影響。

Fig.16 Contours of mean values and variances of Mach number

1）從單輸入變量入手，研究了卡門常數(shù)κ的不確定度對(duì)數(shù)值模擬的影響。研究發(fā)現(xiàn)升力、阻力系數(shù)對(duì)κ是近似線性的響應(yīng)關(guān)系。在單變量情況下，利用高斯積分法和回歸法計(jì)算PC展開(kāi)的自由度，兩者結(jié)果也基本一致。κ的不確定度主要影響翼型上表面，特別是激波附近的壓力分布，同時(shí)還影響翼型上、下表面的摩擦預(yù)測(cè)結(jié)果和激波附近馬赫數(shù)的預(yù)測(cè)結(jié)果。κ減小時(shí)，預(yù)測(cè)的激波位置后移，摩擦系數(shù)預(yù)測(cè)值變小。

2）同時(shí)考慮模型中9個(gè)參數(shù)的不確定性帶來(lái)的影響。通過(guò)升力、阻力系數(shù)預(yù)測(cè)中9個(gè)參數(shù)各自的Sobol指標(biāo)可以發(fā)現(xiàn)，ct3、ct4、cb2、cw3這4個(gè)參數(shù)對(duì)升力、阻力系數(shù)方差的貢獻(xiàn)很小，可以忽略不計(jì)。其他5個(gè)參數(shù)對(duì)于結(jié)果的不確定度影響都相對(duì)比較大，κ是其中貢獻(xiàn)最大的參數(shù)。建立的PC能夠給出壁面壓力分布、摩擦系數(shù)分布和馬赫數(shù)分布的不確定度。

3）通過(guò)PC展開(kāi)和CFD計(jì)算的對(duì)比證實(shí)了建立的PC展開(kāi)能夠準(zhǔn)確地描述輸出對(duì)輸入的響應(yīng)。這為量化計(jì)算結(jié)果的不確定度提供了切實(shí)可行、有較高效率的解決方案。需要指出的是本文的計(jì)算只考慮了RAE2822跨聲速翼型模擬的單一計(jì)算狀態(tài)，影響規(guī)律是否可以推及其他工況和算例需要進(jìn)一步檢驗(yàn)。