橫流轉(zhuǎn)捩模型參數(shù)不確定度量化分析與應(yīng)用研究

向星皓，張毅鋒，陳堅(jiān)強(qiáng),2，袁先旭,2，陳樹生

(1. 中國(guó)空氣動(dòng)力研究與發(fā)展中心空氣動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，綿陽(yáng) 621000；2. 中國(guó)空氣動(dòng)力研究與發(fā)展中心計(jì)算空氣動(dòng)力研究所， 綿陽(yáng) 621000；3. 西北工業(yè)大學(xué)航空學(xué)院，西安 710072)

0 引 言

邊界層轉(zhuǎn)捩通常是指邊界層流動(dòng)由層流狀態(tài)發(fā)展為湍流狀態(tài)的過(guò)程，是一個(gè)多因素耦合影響、包含復(fù)雜轉(zhuǎn)捩機(jī)理的物理現(xiàn)象[1]。真實(shí)飛行器的三維邊界層轉(zhuǎn)捩往往由橫流轉(zhuǎn)捩主導(dǎo)[2]。作為轉(zhuǎn)捩預(yù)測(cè)手段之一的轉(zhuǎn)捩模型在構(gòu)造過(guò)程中涉及到較多模型參數(shù)，大部分參數(shù)是通過(guò)特定條件下的風(fēng)洞實(shí)驗(yàn)和理論分析來(lái)確定的，其適用范圍有限，在實(shí)際應(yīng)用中往往需要針對(duì)不同流動(dòng)類型進(jìn)行參數(shù)調(diào)整。

由于實(shí)驗(yàn)數(shù)據(jù)獲取困難以及缺乏對(duì)轉(zhuǎn)捩機(jī)理的深入理解，在轉(zhuǎn)捩模型的構(gòu)造過(guò)程中存在一定的認(rèn)知不確定性。雖然該不確定性會(huì)隨著人們對(duì)物理模型了解的深入和實(shí)驗(yàn)數(shù)據(jù)的豐富而逐步減小，但是目前仍會(huì)導(dǎo)致模型參數(shù)不確定性。這主要體現(xiàn)在給定的模型參數(shù)未必最優(yōu)，針對(duì)具體轉(zhuǎn)捩類型與流場(chǎng)特性，需要對(duì)參數(shù)進(jìn)行標(biāo)定與修正。

模型參數(shù)的標(biāo)定通常是確定性的，給出了參數(shù)具體取值或取值范圍。如文獻(xiàn)[3]用特定工況的理論分析和風(fēng)洞實(shí)驗(yàn)，結(jié)合建模中的部分假設(shè)直接確定模型參數(shù)取值。也有采用類似試錯(cuò)法進(jìn)行參數(shù)調(diào)整，如Lien等[4]在建模以后對(duì)參數(shù)進(jìn)行變參數(shù)調(diào)整優(yōu)化，在特定算例中對(duì)模型參數(shù)進(jìn)行試錯(cuò)調(diào)整，實(shí)現(xiàn)一個(gè)較小范圍的參數(shù)最優(yōu)。

不確定度量化分析相比以往確定性的參數(shù)標(biāo)定，能夠提供多參數(shù)對(duì)轉(zhuǎn)捩起始位置耦合作用的定量影響[5]、能夠定量給出計(jì)算結(jié)果不確定度范圍以及所需樣本量小[6]等諸多優(yōu)勢(shì)，對(duì)于參數(shù)標(biāo)定工作以及模型參數(shù)效能分析具有重要的指導(dǎo)意義。不確定度分析具有非侵入式以及對(duì)隨機(jī)輸入變量具有指數(shù)收斂性的特點(diǎn)[6]，能夠提供輸入條件對(duì)轉(zhuǎn)捩起始位置的定量影響[5]。

由于認(rèn)知不確定性造成的轉(zhuǎn)捩模型參數(shù)不確定性因素普遍存在，不確定度量化分析就顯得尤為必要。張涵信院士等[7]指出，CFD方法的不確定度或可靠性一直以來(lái)都是CFD中需要特別關(guān)心的問(wèn)題。不確定度量化分析對(duì)于轉(zhuǎn)捩模型參數(shù)研究具有重要意義，能夠提供參數(shù)敏感性以及不確定度的定量結(jié)論，進(jìn)而對(duì)參數(shù)效能研究以及模型參數(shù)的篩選調(diào)試提供具體有效的指導(dǎo)。在此基礎(chǔ)上，不確定度量化分析方法可以甄別對(duì)計(jì)算結(jié)果影響大的關(guān)鍵因素，為轉(zhuǎn)捩模型的改進(jìn)指明方向。

在轉(zhuǎn)捩/湍流模型的不確定度量化分析方面，目前關(guān)于模型不確定度研究主要集中在湍流模型方面，對(duì)于轉(zhuǎn)捩模型的不確定度研究則相對(duì)有限。Pecnik等[5]對(duì)由于自由來(lái)流不確定性與可壓縮修正所導(dǎo)致的轉(zhuǎn)捩模型不確定度進(jìn)行了分析，分析發(fā)現(xiàn)上述因素對(duì)轉(zhuǎn)捩起始位置與氣動(dòng)熱分布的不確定度造成顯著影響。Zhao等[8]對(duì)轉(zhuǎn)捩模型在高超聲速平板和尖錐中的不確定度和參數(shù)敏感性開展研究，采用非嵌入式多項(xiàng)式混沌(Non-intrusive polgnomial chaos,NIPC)方法對(duì)γ-Reθt和k-ω-γ轉(zhuǎn)捩模型進(jìn)行了基于來(lái)流不確定性的不確定度量化分析。

傳統(tǒng)不確定度分析方法如工程應(yīng)用較為普遍的蒙特-卡洛(Monte-Carlo)方法，在對(duì)計(jì)算流體力學(xué)氣動(dòng)力、熱不確定度分析時(shí)，樣本需求量大，計(jì)算資源消耗大[9]。相對(duì)于蒙特-卡洛方法而言，非嵌入式多項(xiàng)式混沌(NIPC)方法不需要修改求解器，具有計(jì)算量小、收斂快的特點(diǎn)，在復(fù)雜模型、系統(tǒng)的不確定度量化分析中得到廣泛應(yīng)用[10]。

本文首先基于課題組Chant 2.0計(jì)算平臺(tái)[11]實(shí)現(xiàn)了當(dāng)?shù)鼗瘷M流轉(zhuǎn)捩模型[12-13]，然后采用非嵌入式多項(xiàng)式混沌方法對(duì)NLF(2)- 0415后掠翼和S-K低速平板進(jìn)行了模型參數(shù)不確定度量化分析和參數(shù)敏感性分析，得到了當(dāng)?shù)鼗臋M流轉(zhuǎn)捩模型參數(shù)不確定性的定量分析結(jié)果。以分析結(jié)果為指導(dǎo)對(duì)橫流轉(zhuǎn)捩模型進(jìn)行了重新標(biāo)定與參數(shù)修正，最后將修正的轉(zhuǎn)捩模型在NLF(2)- 0415后掠翼、6:1標(biāo)準(zhǔn)橢球體和DLR-F4翼身組合體算例中進(jìn)行了驗(yàn)證。計(jì)算結(jié)果表明，重新標(biāo)定的模型對(duì)于亞聲速、跨聲速情況下的三維邊界層轉(zhuǎn)捩具有較好適用性。

1 計(jì)算方法

1.1 控制方程

本文采用有限體積法對(duì)雷諾平均Navier-Stokes(RANS)方程進(jìn)行求解，無(wú)黏通量計(jì)算格式為AUSMPW+，無(wú)黏項(xiàng)離散格式為二階精度NND格式，黏性項(xiàng)離散格式為二階中心差分格式，時(shí)間推進(jìn)采用LU-SGS方法，湍流模型采用k-ωSST模型，采用MPI技術(shù)進(jìn)行大規(guī)模并行計(jì)算。

1.2 橫流轉(zhuǎn)捩模型的不確定度輸入?yún)?shù)選取

本文在Chant 2.0計(jì)算平臺(tái)實(shí)現(xiàn)了文獻(xiàn)[12-13]提出的當(dāng)?shù)鼗瘷M流轉(zhuǎn)捩預(yù)測(cè)模型。由于計(jì)算平臺(tái)差異，模型參數(shù)需要重新標(biāo)定。對(duì)模型開展了轉(zhuǎn)捩相關(guān)量不確定度量化分析和參數(shù)研究，以指導(dǎo)參數(shù)標(biāo)定工作。橫流模型基于γ-Reθt轉(zhuǎn)捩模型進(jìn)行拓展，橫流效應(yīng)是通過(guò)在原始模型的Reθt輸運(yùn)方程中增加橫流源項(xiàng)來(lái)實(shí)現(xiàn)的。

(1)

(2)

(3)

cCF=0.6

(4)

在橫流源項(xiàng)DSCF中有兩個(gè)模型參數(shù)作為橫流源項(xiàng)的系數(shù)而存在，分別是cCF和cθt。上述參數(shù)大小會(huì)直接影響輸運(yùn)方程中的橫流源項(xiàng)大小，影響模型預(yù)測(cè)的橫流轉(zhuǎn)捩起始位置，是橫流轉(zhuǎn)捩不確定性研究的重點(diǎn)關(guān)注參數(shù)。

物面粗糙度作為物體表面固有屬性參數(shù)，在大多數(shù)風(fēng)洞實(shí)驗(yàn)、飛行試驗(yàn)以及CFD橫流模型預(yù)測(cè)中，這一重要參數(shù)往往是缺省的，或者測(cè)量存在一定的誤差。在本文所用的橫流轉(zhuǎn)捩模型中，表面粗糙度h是定常橫流轉(zhuǎn)捩判據(jù)ReSCF的重要參數(shù)之一，見式(5)。本文將粗糙度h作為輸入?yún)?shù)之一進(jìn)行不確定度和參數(shù)敏感性分析。

319.51+f(+ΔHCF)-f(-ΔHCF)

(5)

綜上所述，本文選取模型參數(shù)cCF和cθt以及表面粗糙度h進(jìn)行不確定度與參數(shù)敏感性分析，對(duì)模型參數(shù)進(jìn)行篩選與橫流標(biāo)定。橫流轉(zhuǎn)捩模型的其余參數(shù)則與γ-Reθt模型一致，見文獻(xiàn)[12]。

1.3 不確定度分析方法

傳統(tǒng)的CFD數(shù)值模擬具有確定性，即針對(duì)一個(gè)確定的問(wèn)題通過(guò)一組確定性的輸入從而得到一個(gè)確定性的計(jì)算結(jié)果[10]。對(duì)于轉(zhuǎn)捩模型而言，由于轉(zhuǎn)捩機(jī)理的高度復(fù)雜性和在風(fēng)洞實(shí)驗(yàn)和飛行試驗(yàn)中測(cè)量的困難，以及對(duì)模型認(rèn)知存在的局限，導(dǎo)致構(gòu)建的轉(zhuǎn)捩模型的參數(shù)也存在較大不確定性[11]。因此有必要針對(duì)轉(zhuǎn)捩模型參數(shù)所導(dǎo)致的轉(zhuǎn)捩起始位置和轉(zhuǎn)捩特征量的預(yù)測(cè)不確定性進(jìn)行不確定度量化分析。

本文采用非嵌入式多項(xiàng)式混沌方法進(jìn)行不確定度量化分析。在CFD轉(zhuǎn)捩預(yù)測(cè)模型不確定度分析中，以壁面摩擦力系數(shù)Cf、轉(zhuǎn)捩起始位置等重要參數(shù)作為隨機(jī)輸出變量，表示為：

(6)

式中:α*是CFD直接計(jì)算結(jié)果，αj(x)是計(jì)算結(jié)果的確定部分耦合系數(shù)，ψj(ξ)是計(jì)算結(jié)果的隨機(jī)部分。ξ=(ξ1,…,ξn)是n維隨機(jī)變量，而隨機(jī)部分ψj(ξ)是以隨機(jī)變量ξ為自變量的隨機(jī)函數(shù)，為ξ的正交多項(xiàng)式。

隨機(jī)函數(shù)ψj(ξ)根據(jù)隨機(jī)變量的ξ分布，即本文輸入的轉(zhuǎn)捩模型參數(shù)的分布類型不同而有所區(qū)別，對(duì)應(yīng)選取不同形式的正交多項(xiàng)式。若模型參數(shù)滿足均勻分布時(shí)，正交多項(xiàng)式選擇Legendre正交多項(xiàng)式；若模型參數(shù)滿足正態(tài)分布時(shí)，正交多項(xiàng)式選取Hermite正交多項(xiàng)式[9]。

對(duì)多項(xiàng)式采取p階截?cái)?，設(shè)隨機(jī)參數(shù)的維數(shù)為n，則混沌多項(xiàng)式(PCE)項(xiàng)數(shù)可以表示為：

(7)

本文選取隨機(jī)響應(yīng)面法求解混沌多項(xiàng)式系數(shù)，開展不確定度分析。參數(shù)樣本量的設(shè)置，參考文獻(xiàn)[14]選用PCE系數(shù)兩倍的過(guò)采樣方法。根據(jù)Schaefer等[15]的比較結(jié)果，采用精度和收斂性均表現(xiàn)較好的拉丁超立方(LHD)抽樣方法，對(duì)模型參數(shù)進(jìn)行抽樣選取。

選取了Nt個(gè)隨機(jī)樣本點(diǎn)后，每一個(gè)樣本點(diǎn)對(duì)應(yīng)一個(gè)確定的CFD計(jì)算結(jié)果，由式(6)可得[9]：

(8)

采用最小二次回歸對(duì)PCE系數(shù)αi進(jìn)行求解[9]，平均值μ和方差D按如下計(jì)算：

μ=α0(x)

(9)

(10)

每一個(gè)輸入?yún)?shù)變量i對(duì)輸出變量不確定度貢獻(xiàn)的相對(duì)大小是通過(guò)敏感性指數(shù)來(lái)表征的。Sobol指數(shù)(STi)作為敏感性指數(shù)定義為部分方差與總方差的比值[16]：

(11)

其中部分方差與總方差分別為：

(12)

(13)

針對(duì)輸入?yún)?shù)i的Sobol指數(shù)(STi)則定義為包含變量i的所有部分Sobol指數(shù)之和：

(14)

2 不確定度與參數(shù)敏感性分析

本文選取S-K低速平板和NFL(2)- 0415后掠翼兩個(gè)算例進(jìn)行不確定度量化和參數(shù)敏感性分析。第一個(gè)算例的邊界層轉(zhuǎn)捩由T-S波主導(dǎo)，第二個(gè)算例是由橫流不穩(wěn)定性主導(dǎo)。通過(guò)上述兩個(gè)算例的NIPC分析，能夠得到模型參數(shù)在不同轉(zhuǎn)捩類型中的效能，從而進(jìn)行有針對(duì)性的參數(shù)修正工作。

選取三個(gè)模型參數(shù)作為輸入，選取兩個(gè)物理量作為輸出響應(yīng)，進(jìn)行不確定度量化分析。兩個(gè)輸出響應(yīng)為：轉(zhuǎn)捩起始位置xtr和壁面摩擦力系數(shù)Cf。三個(gè)輸入?yún)?shù)及依據(jù)為：1)動(dòng)量厚度雷諾數(shù)輸運(yùn)方程中橫流破壞項(xiàng)系數(shù)參數(shù)cCF,該參數(shù)作為橫流耗散項(xiàng)DSCF系數(shù)(見式(2))，對(duì)模型的橫流轉(zhuǎn)捩預(yù)測(cè)影響較大；2)表面粗糙度h，該參數(shù)構(gòu)成橫流轉(zhuǎn)捩判據(jù)ReSCF迭代式，必然對(duì)橫流轉(zhuǎn)捩產(chǎn)生影響；3)γ-Reθt轉(zhuǎn)捩模型和橫流模塊DSCF項(xiàng)共有的系數(shù)參數(shù)cθt。

上述三個(gè)參數(shù)對(duì)橫流轉(zhuǎn)捩存在重要影響是顯而易見的，本文目的主要在于采用NIPC方法對(duì)參數(shù)影響進(jìn)行定量分析，得出不確定度和參數(shù)敏感性的定量結(jié)論，進(jìn)一步指導(dǎo)參數(shù)標(biāo)定工作。

輸入?yún)?shù)的不確定度范圍設(shè)置如下：橫流參數(shù)cCF基準(zhǔn)值根據(jù)文獻(xiàn)[13]設(shè)置為原始值0.6，正負(fù)偏差百分比根據(jù)模型標(biāo)定經(jīng)驗(yàn)[17]設(shè)置為±66.66%，對(duì)應(yīng)參數(shù)變化范圍[0.2，1.0]；參數(shù)cθt基準(zhǔn)值設(shè)置為0.03[12]，不確定度正負(fù)偏差百分比為±20%，對(duì)應(yīng)參數(shù)變化范圍[0.024，0.036]；表面粗糙度h根據(jù)機(jī)加工精度覆蓋從不可辨加工痕跡的超級(jí)加工光澤面至粗糙漆面[13,18-19]，粗糙度變化范圍為[0.5，10]，單位μm。

混沌多項(xiàng)式(PCE)采用2階截?cái)?，根?jù)式(7)采用過(guò)采樣方法，選取樣本數(shù)為20，采用拉丁超立方進(jìn)行抽樣，選擇轉(zhuǎn)捩起始位置和壁面摩阻系數(shù)作為輸出響應(yīng)。在95%的置信區(qū)間下，轉(zhuǎn)捩起始位置與摩阻系數(shù)Cf輸出響應(yīng)的不確定度為UQ%=100×1.96σR/μR。

2.1 S-K低速平板算例

S-K低速平板來(lái)流狀態(tài)為Ma=0.147,Re=3.34×106/m,來(lái)流湍流度Tu∞=0.18%。以摩阻系數(shù)Cf為輸出響應(yīng)進(jìn)行了基于三參數(shù)的不確定度與參數(shù)敏感性分析，通過(guò)不同站位上三個(gè)輸入?yún)?shù)對(duì)應(yīng)的敏感性指數(shù)比較，可以觀察輸入?yún)?shù)在不同流動(dòng)過(guò)程的影響程度。由圖1可知，在轉(zhuǎn)捩區(qū)(0.8

圖1 低速平板算例壁面Cf不確定度與參數(shù)敏感指數(shù)Fig.1 UQ% & sobol index of Cf in S-K flat plate

與摩阻系數(shù)類似，以轉(zhuǎn)捩起始位置為輸出響應(yīng)也可以給出三個(gè)參數(shù)的敏感性指數(shù)，計(jì)算結(jié)果見表1，模型參數(shù)cθt的敏感性指數(shù)大約是參數(shù)cCF與物面粗糙度h的敏感性指數(shù)的2倍，與圖1的結(jié)果一致。在相對(duì)幅值變化相同的條件下，參數(shù)cθt所引起的轉(zhuǎn)捩起始位置xtr的變化同樣更為顯著。

表1 低速平板算例轉(zhuǎn)捩起始位置的三參數(shù)敏感性指數(shù)Table 1 Sobol index in S-K flat plate case

該算例說(shuō)明，cθt是影響T-S不穩(wěn)定性主導(dǎo)轉(zhuǎn)捩計(jì)算的主要參數(shù)，調(diào)整模型參數(shù)cθt獲得的敏感性高于其他兩個(gè)參數(shù)，在對(duì)流向轉(zhuǎn)捩進(jìn)行標(biāo)定時(shí)，應(yīng)主要考慮參數(shù)cθt。

2.2 NFL(2)- 0415后掠翼算例

2.2.1模型參數(shù)不確定度量化

NFL(2)- 0415后掠翼上表面轉(zhuǎn)捩過(guò)程由橫流不穩(wěn)定性主導(dǎo)，選取上壁面摩阻系數(shù)Cf和轉(zhuǎn)捩起始位置xtr作為輸出響應(yīng)進(jìn)行分析。不確定度分析的NFL(2)- 0415后掠翼算例計(jì)算狀態(tài)為攻角-4°,Ma=0.2，Re=3.27×106/m，來(lái)流湍流度Tu∞=0.09%，邊界層遠(yuǎn)場(chǎng)湍流度自由衰減至Tu≤0.02%。

圖2是所有拉丁超立方抽樣樣本的上表面摩阻系數(shù)Cf分布圖，并給出NIPC均值進(jìn)行對(duì)照。在20個(gè)樣本曲線中，轉(zhuǎn)捩起始位置在x=0.43附近有5條Cf曲線基本重合，對(duì)應(yīng)樣本編號(hào)為7,12,13,19,20(見表2)。它們中的參數(shù)cCF與參數(shù)h有此消彼長(zhǎng)的關(guān)系，說(shuō)明兩參數(shù)對(duì)橫流轉(zhuǎn)捩的影響作用近似。

圖2 采樣樣本的上壁面Cf分布圖Fig.2 Cf distribution of all samples

表2則是采用拉丁超立方抽樣歸一化參數(shù)樣本(參數(shù)樣本：最大值→1，最小值→-1)。該不確定度分析手段具有“樣本量小、信息豐富”的特點(diǎn)。與常規(guī)變參分析進(jìn)行對(duì)比，若采用常規(guī)單參數(shù)調(diào)試的方法，假設(shè)每個(gè)參數(shù)在指定范圍內(nèi)進(jìn)行六個(gè)均勻采樣，則在樣本量相當(dāng)?shù)那闆r下，傳統(tǒng)方法只能得到單參數(shù)調(diào)試結(jié)果，無(wú)法得到參數(shù)耦合變化下的不確定度與參數(shù)敏感性分析結(jié)果(見圖3、圖4)。

表2 拉丁超立方抽樣輸入?yún)?shù)及轉(zhuǎn)捩起始位置Table 2 LHD sampling and transition location

圖3是基于參數(shù)抽樣結(jié)果的壁面摩阻系數(shù)均值及誤差帶分布。其意義在于僅采用了少量樣本，即提供了變參數(shù)橫流轉(zhuǎn)捩模型輸出量的整個(gè)空間分布及其變化范圍。如采用傳統(tǒng)變參分析，其參數(shù)樣本量與計(jì)算量將遠(yuǎn)大于此。

由圖3可知，在x<0.4的層流區(qū)和x>0.7以后的充分發(fā)展湍流區(qū)，誤差帶范圍明顯低于在0.4

圖3 NIPC方法壁面摩擦系數(shù)均值及誤差帶分布Fig.3 Surface Cf distribution with NIPC method

圖4 后掠翼算例Cf不確定度與參數(shù)敏感指數(shù)分布Fig.4 UQ & Sobol Index of Cf on back-swept wing

將上述參數(shù)采樣和對(duì)應(yīng)計(jì)算結(jié)果作為輸入與輸出，采用NIPC方法進(jìn)行不確定度分析。圖4為摩阻Cf作為輸出響應(yīng)的不確定度與參數(shù)敏感性指數(shù)沿弦向分布，敏感性在第2.2.2節(jié)進(jìn)行分析。不確定度的分布與圖3誤差帶區(qū)間分布具有一致性，在x<0.4與x>0.7的區(qū)間近似為0，在轉(zhuǎn)捩區(qū)間呈現(xiàn)出“兩端小、中間大”的分布規(guī)律。轉(zhuǎn)捩區(qū)域中x/c=0.43時(shí)Cf不確定度達(dá)到最大，大小約為189%。

2.2.2模型參數(shù)敏感性分析

對(duì)轉(zhuǎn)捩起始位置xtr和壁面摩擦系數(shù)Cf分布分別進(jìn)行了參數(shù)敏感性的點(diǎn)分析和線分析。二者反映不同物理量對(duì)參數(shù)變化的敏感程度，轉(zhuǎn)捩起始位置三參數(shù)點(diǎn)分析的敏感性指數(shù)(見表3)，不能通過(guò)壁面摩擦系數(shù)敏感性指數(shù)分布(見圖4)直接得到。

以轉(zhuǎn)捩起始位置xtr為輸出響應(yīng)，進(jìn)行參數(shù)敏感性的點(diǎn)分析。表3給出了經(jīng)不確定度點(diǎn)分析的轉(zhuǎn)捩起始位置敏感性指數(shù)，三參數(shù)的Sobol指數(shù)分別為0.56153，0.43323，0.12676，該指數(shù)表征參數(shù)敏感程度，參數(shù)cCF與h的敏感性指數(shù)約為參數(shù)cθt敏感性指數(shù)的3～4倍，其對(duì)橫流轉(zhuǎn)捩起始位置的影響要遠(yuǎn)大于參數(shù)cθt。

表3 后掠翼轉(zhuǎn)捩起始位置的三參數(shù)敏感性指數(shù)Table 3 Sobol index in back-swept wing case

以摩阻Cf為輸出響應(yīng)，進(jìn)行參數(shù)敏感性的線分析見圖4。從圖4可以看出，對(duì)轉(zhuǎn)捩區(qū)間Cf分布影響較大的參數(shù)是cCF與表面粗糙度h，其敏感性指數(shù)在轉(zhuǎn)捩區(qū)域內(nèi)(0.40.8的區(qū)間，cθt的敏感性指數(shù)趨于0，即cθt幾乎不影響轉(zhuǎn)捩前的層流區(qū)域和轉(zhuǎn)捩后的湍流區(qū)域的Cf分布。但cCF,h與cθt的敏感性指數(shù)規(guī)律不同，它們僅在x<0.4的層流區(qū)域近似為0，在x>0.8轉(zhuǎn)捩完成以后的湍流區(qū)域，cCF敏感性指數(shù)約為0.3，h的敏感性指數(shù)約為0.4。說(shuō)明橫流參數(shù)cCF與表面粗糙度h不僅影響轉(zhuǎn)捩起始位置，還會(huì)對(duì)轉(zhuǎn)捩之后的Cf分布造成一定影響。與該觀點(diǎn)相互印證的是圖3中誤差帶分析結(jié)果，x∈[0.8,1]區(qū)域內(nèi)誤差帶不為零，該區(qū)域內(nèi)的不確定度大約在4%～8%的范圍內(nèi)。調(diào)整參數(shù)cCF和h，能使轉(zhuǎn)捩起始位置相對(duì)于流向轉(zhuǎn)捩提前發(fā)生，模擬橫流轉(zhuǎn)捩效果。

2.3 不確定度與參數(shù)敏感性分析小結(jié)

相比常規(guī)變參分析方法，不確定度量化分析方法具有如下顯著優(yōu)勢(shì)：

1) 不確定度量化分析能夠給出研究者所關(guān)心的輸出物理量的不確定度定量范圍，具有重要的工程指導(dǎo)意義。

2) 不確定度量化分析所給出的參數(shù)敏感性指數(shù)，可以直接量化參數(shù)在特定空間位置對(duì)特定物理量的影響大小，而常規(guī)變參分析只能給出影響大小的宏觀定性排序。前者對(duì)于參數(shù)篩選與標(biāo)定工作更具指導(dǎo)意義。

3) 基于NIPC的不確定度分析，采用拉丁超立方抽樣需要的樣本量小，得到的信息豐富，效費(fèi)比高，能夠綜合考慮多參數(shù)耦合變化時(shí)的參數(shù)作用效應(yīng)。常規(guī)變參數(shù)分析不具備上述特點(diǎn)。

三輸入?yún)?shù)以不同的作用方式都對(duì)轉(zhuǎn)捩模擬結(jié)果有影響，不同的轉(zhuǎn)捩方式需采用不同的參數(shù)進(jìn)行標(biāo)定：

1) 參數(shù)cθt對(duì)流向轉(zhuǎn)捩相關(guān)量的不確定度影響較大，敏感性較高，且對(duì)于橫流轉(zhuǎn)捩影響較小。流向轉(zhuǎn)捩模型的標(biāo)定應(yīng)主要針對(duì)參數(shù)cθt開展。

2) 參數(shù)cCF對(duì)橫流轉(zhuǎn)捩相關(guān)量的不確定度影響較大，敏感性較高，且cCF對(duì)于流向轉(zhuǎn)捩影響較小，對(duì)橫流模型的標(biāo)定應(yīng)主要針對(duì)參數(shù)cCF開展。

3) 表面粗糙度h是反映模型加工表面光潔程度的固有參數(shù)，對(duì)橫流轉(zhuǎn)捩位置xtr以及Cf分布均存在顯著影響。不同風(fēng)洞實(shí)驗(yàn)、飛行試驗(yàn)的模型表面由于加工精度不同，粗糙度差別較大，需要對(duì)其精確測(cè)量。

3 橫流轉(zhuǎn)捩模型參數(shù)標(biāo)定與校驗(yàn)

本文主要研究橫流轉(zhuǎn)捩模型參數(shù)，根據(jù)不確定度和參數(shù)敏感性分析結(jié)果，對(duì)模型參數(shù)進(jìn)行篩選，給出參數(shù)cCF的標(biāo)定和校驗(yàn)結(jié)果。

3.1 橫流模型參數(shù)標(biāo)定

針對(duì)不同計(jì)算平臺(tái)的模型參數(shù)標(biāo)定的廣泛需求，以不確定度和參數(shù)敏感性分析結(jié)果為指導(dǎo)，進(jìn)行模型參數(shù)篩選和標(biāo)定。

NFL(2)- 0415后掠翼算例來(lái)流狀態(tài)為Ma=0.2,Re=1.9～3.3×106/m (1.93×106, 2.19×106, 2.40×106, 2.73×106, 3.27×106)。在Chant 2.0計(jì)算平臺(tái)上采用基準(zhǔn)參數(shù)的橫流模型，由于計(jì)算平臺(tái)差異，采用原始基準(zhǔn)參數(shù)的模型預(yù)測(cè)的NFL- 0415(2)后掠翼轉(zhuǎn)捩起始位置與實(shí)驗(yàn)結(jié)果偏差較大，如圖5所示。在雷諾數(shù)為3.27×106算例中轉(zhuǎn)捩起始位置偏差大約在30%以上。根據(jù)現(xiàn)有的NFL- 0415后掠翼轉(zhuǎn)捩實(shí)驗(yàn)數(shù)據(jù)[19]，采用改變模型參數(shù)值的方法調(diào)整橫流轉(zhuǎn)捩起始位置，對(duì)橫流轉(zhuǎn)捩模型進(jìn)行重新標(biāo)定。表面粗糙度設(shè)置為實(shí)驗(yàn)測(cè)量結(jié)果[13,19](h=3.3 μm)。

圖5 基準(zhǔn)參數(shù)模型預(yù)測(cè)轉(zhuǎn)捩起始位置與實(shí)驗(yàn)對(duì)比Fig.5 Transition location of modeling and experiment

由第2節(jié)不確定度與參數(shù)敏感性分析結(jié)果可知，在橫流不穩(wěn)定性主導(dǎo)的轉(zhuǎn)捩中，模型參數(shù)cCF和cθt均影響轉(zhuǎn)捩位置，但cCF的參數(shù)敏感性指數(shù)更高(見表1和圖4)。同時(shí)cCF僅是動(dòng)量厚度輸運(yùn)方程橫流源項(xiàng)系數(shù)，該系數(shù)的調(diào)整與橫流轉(zhuǎn)捩位置變化具有單調(diào)一致性。因此選取參數(shù)cCF進(jìn)行橫流標(biāo)定。

雖然常規(guī)變參分析可以選取橫流耗散項(xiàng)DSCF系數(shù)cCF進(jìn)行直接標(biāo)定，但參數(shù)選取及變參數(shù)調(diào)試缺乏定量的誤差帶和敏感性數(shù)據(jù)支持，經(jīng)驗(yàn)性較大。

圖6 變參數(shù)的后掠翼的計(jì)算和實(shí)驗(yàn)結(jié)果[19]Fig.6 Experiment & modeling results with varying parameter values[19]

就參數(shù)cCF對(duì)輸運(yùn)方程的影響進(jìn)行定性分析：該系數(shù)越小，Reθt輸運(yùn)方程橫流破壞項(xiàng)占比越低，由間歇因子輸運(yùn)方程啟動(dòng)項(xiàng)對(duì)Reθt的單調(diào)性，減小該系數(shù)使轉(zhuǎn)捩起始位置后移。極限情況下cCF參數(shù)為0，橫流模型恢復(fù)至原始γ-Reθt轉(zhuǎn)捩模型。

根據(jù)上述分析以及圖3的采樣樣本結(jié)果，以表1樣本17, 3, 8為參考，設(shè)置參數(shù)cCF的標(biāo)定范圍為[0，0.6]，選取cCF=0, 0.05, 0.1, 0.15, 0.2, 0.4, 0.6七個(gè)參數(shù)值在雷諾數(shù)Re=3.27×106算例中進(jìn)行精細(xì)化標(biāo)定。參數(shù)值取0.2時(shí)與實(shí)驗(yàn)數(shù)據(jù)擬合最為準(zhǔn)確，進(jìn)一步在較大雷諾數(shù)范圍進(jìn)行測(cè)試，表面粗糙度h=3.3 μm，Re=2.37×106與Re=3.27×106的轉(zhuǎn)捩起始位置偏差都得到了有效修正，轉(zhuǎn)捩起始位置與風(fēng)洞實(shí)驗(yàn)符合較好(見圖6)。

3.2 橫流轉(zhuǎn)捩算例

前文基于固定粗糙度的后掠翼算例對(duì)橫流模型參數(shù)進(jìn)行了標(biāo)定，本節(jié)對(duì)三種粗糙度、不同雷諾數(shù)條件下的后掠翼以及15°攻角橢球體和DLR-F4翼身組合體進(jìn)行計(jì)算，以驗(yàn)證橫流模型參數(shù)的適用性。

3.2.1NLF(2)- 0415后掠翼

NFL(2)- 0415無(wú)限展長(zhǎng)后掠翼在-4°攻角下，機(jī)翼上表面的轉(zhuǎn)捩由橫流不穩(wěn)定性主導(dǎo)，是低速橫流轉(zhuǎn)捩典型算例。選取該算例，采用標(biāo)定后的模型對(duì)變粗糙度的風(fēng)洞實(shí)驗(yàn)[19]轉(zhuǎn)捩位置進(jìn)行預(yù)測(cè)，驗(yàn)證模型在不同粗糙度下橫流轉(zhuǎn)捩模擬能力。攻角為-4°，后掠角45°，雷諾數(shù)范圍為Re=1.9×106～3.8×106，網(wǎng)格量約50萬(wàn)，物面第一層網(wǎng)格法向間距y+<1，壁面粗糙度(h)根據(jù)風(fēng)洞模型加工精度設(shè)置為0.5 μm,3.3 μm,9 μm，來(lái)流湍流度Tu∞≈0.1%，νt/ν=5，前緣邊界層遠(yuǎn)場(chǎng)湍流度Tu≤0.02%。

圖7是三種粗糙度、多雷諾數(shù)條件下的計(jì)算結(jié)果，預(yù)測(cè)的轉(zhuǎn)捩起始位置隨粗糙度和雷諾數(shù)變化規(guī)律與風(fēng)洞實(shí)驗(yàn)相同，粗糙度促進(jìn)轉(zhuǎn)捩，增大雷諾數(shù)轉(zhuǎn)捩靠前。高雷諾數(shù)下模型預(yù)測(cè)的轉(zhuǎn)捩起始位置和實(shí)驗(yàn)符合較好，低雷諾數(shù)下預(yù)測(cè)的轉(zhuǎn)捩起始位置相比實(shí)驗(yàn)值靠前。

圖7 不同粗糙度下NFL- 0415后掠翼的實(shí)驗(yàn) 和計(jì)算結(jié)果對(duì)比[13,19]Fig.7 Modeling and experimental results of different surface roughness values[13,19]

3.2.2帶傾角的標(biāo)準(zhǔn)橢球體

采用6:1標(biāo)準(zhǔn)橢球體[20]對(duì)重新標(biāo)定的橫流模型在大攻角條件下的橫流轉(zhuǎn)捩預(yù)測(cè)能力進(jìn)行驗(yàn)證。算例選取攻角15°、雷諾數(shù)6.5×106的典型狀態(tài)進(jìn)行測(cè)試。網(wǎng)格量為200萬(wàn)，物面第一層網(wǎng)格法向間距y+<1。來(lái)流湍流度按衰減公式[12]進(jìn)行設(shè)定，保證靠近橢球體的來(lái)流Tu≈0.1%，黏性比νt/ν=5，壁面粗糙度取默認(rèn)值[13]h=3.3 μm。

圖8為實(shí)驗(yàn)測(cè)量以及模型預(yù)測(cè)的橢球體表面摩阻系數(shù)分布云圖。圖8由上至下依次是實(shí)驗(yàn)結(jié)果、γ-Reθt計(jì)算結(jié)果、Chant平臺(tái)橫流模型計(jì)算結(jié)果以及文獻(xiàn)計(jì)算結(jié)果[13]。在高雷諾數(shù)、大攻角且橫流不穩(wěn)定性占主導(dǎo)的情況下，原始γ-Reθt模型的預(yù)測(cè)轉(zhuǎn)捩陣面與實(shí)驗(yàn)存在較大差別，而Chant平臺(tái)標(biāo)定后的橫流模型能較準(zhǔn)確預(yù)測(cè)橢球表面的橫流轉(zhuǎn)捩陣面，轉(zhuǎn)捩區(qū)Cf值比實(shí)驗(yàn)結(jié)果偏低，與文獻(xiàn)橫流模型[13]預(yù)測(cè)能力相當(dāng)。圖9為沿長(zhǎng)軸周向展開后的摩阻系數(shù)分布，轉(zhuǎn)捩模型預(yù)測(cè)橫流轉(zhuǎn)捩陣面與穩(wěn)定性分析和風(fēng)洞實(shí)驗(yàn)符合較好。

圖8 橢球體表面摩擦系數(shù)分布云圖Fig.8 Cf contours on spheroid

圖9 Chant平臺(tái)橫流模型模擬表面摩擦力云圖Fig.9 Cf contour on spheroid: Chant modeling, stability theory and experimental result

3.2.3DLR-F4翼身組合體

DLR-F4翼身組合體是在歐洲跨聲速風(fēng)洞中進(jìn)行的實(shí)驗(yàn)。風(fēng)洞實(shí)驗(yàn)[21]采用溫敏漆技術(shù)顯示層/湍流區(qū)域以及轉(zhuǎn)捩位置。圖10、圖11中實(shí)驗(yàn)照片的明暗界線即為轉(zhuǎn)捩位置。風(fēng)洞實(shí)驗(yàn)Ma=0.785，機(jī)翼上表面為跨聲速流動(dòng)區(qū)域。翼面有轉(zhuǎn)捩發(fā)生，上翼面靠近翼根部分的轉(zhuǎn)捩由橫流不穩(wěn)定性主導(dǎo)，靠近翼尖區(qū)域的轉(zhuǎn)捩由T-S不穩(wěn)定性主導(dǎo)。選取兩個(gè)攻角狀態(tài)：-2.59°，-0.87°，雷諾數(shù)Re=6×106，網(wǎng)格量約為320萬(wàn)，物面第一層網(wǎng)格法向間距y+<1，來(lái)流湍流度Tu∞≈0.1%，保證邊界層遠(yuǎn)場(chǎng)湍流度自由衰減至Tu≤0.05%，黏性比νt/ν=1，表面粗糙度設(shè)置為0.15 μm[5]。

圖10 -2.59°攻角下的轉(zhuǎn)捩位置示意圖Fig.10 Transition location on DLR-F4 with AoA of -2.59°

圖11 -0.87°攻角下的轉(zhuǎn)捩位置示意圖Fig.11 Transition location on DLR-F4 with AoA of -0.87°

圖10、圖11是不同攻角下DLR-F4風(fēng)洞實(shí)驗(yàn)和數(shù)值計(jì)算的轉(zhuǎn)捩位置對(duì)比。轉(zhuǎn)捩模型能夠較準(zhǔn)確預(yù)測(cè)靠近翼根處的橫流不穩(wěn)定性占主導(dǎo)的轉(zhuǎn)捩區(qū)域。圖11所示機(jī)翼中段部位，數(shù)值計(jì)算同實(shí)驗(yàn)存在差距，現(xiàn)有模型預(yù)測(cè)的T-S波轉(zhuǎn)捩過(guò)早發(fā)生。

4 結(jié) 論

本文橫流轉(zhuǎn)捩模型不確定度量化分析與參數(shù)研究結(jié)論如下：

1) 參數(shù)調(diào)整對(duì)轉(zhuǎn)捩不確定度貢獻(xiàn)明顯。在參數(shù)cCF,h與cθt給定范圍后，轉(zhuǎn)捩區(qū)Cf不確定度較大，層流區(qū)為0。低速平板轉(zhuǎn)捩區(qū)Cf不確定度最大達(dá)90%，后掠翼轉(zhuǎn)捩區(qū)不確定度最大達(dá)189%。

2) 根據(jù)參數(shù)靈敏度指數(shù)分析結(jié)果，在橫流轉(zhuǎn)捩中，參數(shù)cCF對(duì)Cf和轉(zhuǎn)捩起始位置的影響要大于參數(shù)cθt的影響。應(yīng)選取橫流參數(shù)cCF對(duì)模型進(jìn)行標(biāo)定。粗糙度h作為固有參數(shù)也會(huì)影響橫流轉(zhuǎn)捩，在風(fēng)洞實(shí)驗(yàn)與飛行試驗(yàn)中應(yīng)進(jìn)行精密測(cè)量。

3) 以不確定度及敏感性分析結(jié)果為指導(dǎo)進(jìn)行重新標(biāo)定后的橫流模型對(duì)變粗糙度NFL- 0415翼型、帶傾角的6:1標(biāo)準(zhǔn)橢球體和DLR-F4翼身組合體的橫流轉(zhuǎn)捩能夠進(jìn)行較好地預(yù)測(cè)。

基于參數(shù)不確定性的橫流轉(zhuǎn)捩模型不確定度與參數(shù)敏感性研究，能夠量化模型參數(shù)對(duì)轉(zhuǎn)捩的影響，提供轉(zhuǎn)捩結(jié)果誤差帶定量分布，具有工程指導(dǎo)價(jià)值。參數(shù)導(dǎo)致的模型不確定度和敏感性定量結(jié)論，對(duì)研究者進(jìn)行模型使用、標(biāo)定與修改具有參考價(jià)值。研究方法具有通用性，可用于各類轉(zhuǎn)捩/湍流模型的不確定度分析，甄別模型的關(guān)鍵參數(shù)，指導(dǎo)模型的精細(xì)化調(diào)試，對(duì)模型的改進(jìn)與工程應(yīng)用有所裨益。