劉錟韜,李瑞宇,高麗敏,*,趙磊
1.西北工業(yè)大學(xué) 動(dòng)力與能源學(xué)院,西安 710072
2.西北工業(yè)大學(xué) 翼型、葉柵空氣動(dòng)力學(xué)國(guó)家級(jí)重點(diǎn)實(shí)驗(yàn)室,西安 710072
3.西安交通大學(xué) 航天航空學(xué)院,西安 710049
對(duì)航空發(fā)動(dòng)機(jī)日益增長(zhǎng)的高可靠性、高性能、高經(jīng)濟(jì)性的需求不斷推動(dòng)著發(fā)動(dòng)機(jī)部件研制過(guò)程中流場(chǎng)預(yù)測(cè)技術(shù)的進(jìn)步。壓氣機(jī)是航空發(fā)動(dòng)機(jī)的壓縮部件,起著對(duì)來(lái)流氣體減速增壓的作用,對(duì)航空發(fā)動(dòng)機(jī)運(yùn)行的性能、安全性和穩(wěn)定性有重要影響。目前,壓氣機(jī)流場(chǎng)預(yù)測(cè)的主要方法有解析全尺度湍流的直接數(shù)值模擬(Direct Numerical Simulation,DNS)[1]、解析大尺度湍流的大渦模擬(Large Eddy Simulation,LES)[2-3]及基于雷諾平均 Navier-Stoke(Reynolds Averaged Navier-Stokes,RANS)方程的湍流模型[4]。DNS 和LES理論上較為完備,對(duì)湍流的刻畫更為準(zhǔn)確,可以獲得流場(chǎng)中精細(xì)的流動(dòng)結(jié)構(gòu),但其巨大的計(jì)算量及數(shù)據(jù)存儲(chǔ)量令其在工程領(lǐng)域的應(yīng)用極為有限,目前僅應(yīng)用在前沿學(xué)術(shù)研究中;在近期,湍流模型憑借其計(jì)算量小、魯棒性好、計(jì)算精度可以接受,仍將是工程領(lǐng)域中流場(chǎng)預(yù)測(cè)的主力軍。但是,基于RANS 的湍流模型大多基于湍流各項(xiàng)同性的渦黏假設(shè),在封閉雷諾應(yīng)力項(xiàng)的過(guò)程中應(yīng)用了大量的假設(shè)和經(jīng)驗(yàn)參數(shù),使湍流模型的計(jì)算結(jié)果都包含有較多的不確定性。這導(dǎo)致每種湍流模型都有一定的應(yīng)用范圍,且對(duì)于有強(qiáng)逆壓力梯度和流動(dòng)分離的復(fù)雜流場(chǎng),湍流模型預(yù)測(cè)雖然能夠反映一定的流動(dòng)特征,但預(yù)測(cè)準(zhǔn)確性相對(duì)較差。壓氣機(jī)內(nèi)包含逆壓力梯度、分離流及激波等復(fù)雜的流動(dòng)現(xiàn)象,在一些非設(shè)計(jì)工況下,預(yù)測(cè)精度的問(wèn)題比較明顯。因此,采用各種方式提高航空發(fā)動(dòng)機(jī)內(nèi)湍流預(yù)測(cè)精度十分有必要。
一般認(rèn)為,湍流模型計(jì)算結(jié)果的不確定主要來(lái)源于4 個(gè)方面[5]:① 雷諾平均過(guò)程;② 使用平均參數(shù)對(duì)雷諾應(yīng)力的表示;③ 特性擬合函數(shù)的選擇;④ 模型系數(shù)的確定。來(lái)源①雖然引入了不確定性,卻是構(gòu)造湍流模型的基礎(chǔ),一般難以修正。目前,通過(guò)改進(jìn)湍流模型提高預(yù)測(cè)精度的研究主要有3 類[6]:① 直接由DNS、LES 或者試驗(yàn)數(shù)據(jù)校正數(shù)值流場(chǎng)中的雷諾應(yīng)力;② 由湍流模型某一項(xiàng)進(jìn)行修正;③ 對(duì)模型系數(shù)進(jìn)行校準(zhǔn),分別對(duì)應(yīng)不確定度來(lái)源的②、③和④。這類研究的思路一般有2 類:① 根據(jù)湍流理論和相關(guān)經(jīng)驗(yàn)修正湍流模型;② 使用數(shù)據(jù)同化(Data Assimilation,DA)[7]方法利用高可信度的數(shù)據(jù)如DNS、LES 或試驗(yàn)測(cè)量數(shù)據(jù)校正RANS 流場(chǎng)和湍流模型。
DA 作為一種融合試驗(yàn)測(cè)量與模型預(yù)測(cè)結(jié)果的數(shù)學(xué)方法,最早起源于氣象領(lǐng)域[8],后被應(yīng)用到海洋科學(xué)[9]和地球物理[10]等學(xué)科,主要用于將衛(wèi)星或者地面站的觀測(cè)數(shù)據(jù)融入預(yù)測(cè)模型,進(jìn)行氣象或者環(huán)境預(yù)測(cè),近年來(lái)也被應(yīng)用到湍流流動(dòng)的增強(qiáng)預(yù)測(cè)[11]研究中。應(yīng)用到湍流預(yù)測(cè)的DA 算法主要有2 類:① 伴隨類算法;② 卡爾曼濾波類算法。伴隨類算法可以對(duì)全場(chǎng)流動(dòng)變量進(jìn)行反演,使預(yù)測(cè)結(jié)果更貼近物理狀態(tài),但需要求解伴隨方程以及改造CFD 求解器,在實(shí)現(xiàn)上比較復(fù)雜,主要用于學(xué)術(shù)研究;卡爾曼濾波類方法來(lái)源于經(jīng)典的卡爾曼濾波算法,既可用于校正模型系數(shù),也適合于全場(chǎng)參數(shù)重構(gòu),不侵入CFD 求解器,更適合工程應(yīng)用。
Foures[12]和He[13]等使用基于伴隨優(yōu)化的數(shù)據(jù)同化方法為RANS 流場(chǎng)添加了額外力項(xiàng)以構(gòu)造各項(xiàng)異性的湍流黏度。張亦知等[14]在S-A 湍流模型的生成項(xiàng)上添加了一個(gè)隨空間分布的修正因子β,利用DNS 計(jì)算結(jié)果及物理知識(shí)約束,通過(guò)伴隨方法確定了β的分布,提高了對(duì)方腔等流動(dòng)摩擦阻力的預(yù)測(cè)精度。 Singh 和Duraisamy[15]使用貝葉斯推斷結(jié)合離散伴隨方法,使用DNS、LES 及試驗(yàn)數(shù)據(jù)對(duì)槽道、彎管和翼型等流場(chǎng)中湍流模型的生成項(xiàng)進(jìn)行了校準(zhǔn),顯著提高了邊界層預(yù)測(cè)精度。但大部分的形式修正項(xiàng)僅為空間分布的函數(shù),僅針對(duì)參與校準(zhǔn)的流動(dòng)有效,缺乏普適性。He 等[16]為提高S-A 模型對(duì)壓氣機(jī)失速裕度的預(yù)測(cè)精度,對(duì)模型方程的剪切應(yīng)變力項(xiàng)乘以一個(gè)經(jīng)驗(yàn)系數(shù),該系數(shù)是壓力梯度和速度旋度的函數(shù),并在多個(gè)壓氣機(jī)算例上進(jìn)行了測(cè)試,具有一定普適性。
湍流模型系數(shù)的默認(rèn)值多數(shù)來(lái)源于平板、翼型等經(jīng)典的湍流案例,在處理復(fù)雜流動(dòng)時(shí)略顯預(yù)測(cè)精度不足。但實(shí)際上,湍流模型本質(zhì)上不能完整地刻畫湍流,所以很難找到一組具有完全普適性的湍流模型系數(shù),校正后的系數(shù)也具有流動(dòng)依賴性,但在一定范圍內(nèi)具有適用性。
王丹華等[17-18]調(diào)整了S-A 模型中Cb1參數(shù),改善了壓氣機(jī)葉柵角區(qū)分離的預(yù)測(cè)效果。馬力等[19]將壓力梯度引入到S-A 模型的剪切力系數(shù)中,亦改善了壓氣機(jī)葉柵角區(qū)分離的預(yù)測(cè)效果。但在根據(jù)經(jīng)驗(yàn)調(diào)整模型系數(shù)時(shí),試驗(yàn)數(shù)據(jù)僅用來(lái)驗(yàn)證效果,并未直接參與調(diào)整過(guò)程,參數(shù)選取和調(diào)整幅度主觀性較強(qiáng);而卡爾曼濾波類數(shù)據(jù)同化方法直接利用試驗(yàn)測(cè)量數(shù)據(jù)校正模型參數(shù),客觀性更強(qiáng),且可以同時(shí)校正多個(gè)參數(shù)。
集合卡爾曼濾波(Ensemble Kalman Filter,EnKF)方法是最常用的一種濾波類數(shù)據(jù)同化方法,由Evensen[20]從經(jīng)典的卡爾曼濾波發(fā)展而來(lái),EnKF 使用Monte-Carlo 方法構(gòu)造樣本協(xié)方差代替模型協(xié)方差矩陣,避免了復(fù)雜的公式推導(dǎo),特別適用于Navier-Stokes 方程等復(fù)雜非線性模型。為了適合不同的應(yīng)用場(chǎng)景和需求,EnKF 發(fā)展出了多個(gè)變種,如集合變換卡爾曼濾波[21](Ensemble Transformed Kalman Filter, ETKF),考慮正則項(xiàng)的正則集合卡爾曼方法[22-23](Regularized Ensemble Kalman Method, REKM)等。
Deng 等[24]使用EnKF 方法對(duì)4 種湍流模型參數(shù)在射流流場(chǎng)中進(jìn)行了修正,原本與試驗(yàn)數(shù)據(jù)符合并欠佳的k-ω和SST(Shear Stress Transport)模型在經(jīng)過(guò)系數(shù)校正后也可以較好地符合。房培勛等[25]使用EnKF 方法對(duì)蒸汽閥門流動(dòng)中SST 湍流模型的系數(shù)進(jìn)行了標(biāo)定,有效提高通流特性預(yù)測(cè)模型的精度,并在相近工況具有一定普適性。Kato 等[26]使用ETKF 對(duì)翼型流場(chǎng)中SST模型a1參數(shù)進(jìn)行了優(yōu)化,認(rèn)為a1=1.0 時(shí)比默認(rèn)值0.31 能更好地預(yù)測(cè)分離流和逆壓力梯度流動(dòng),并且認(rèn)為在參數(shù)優(yōu)化過(guò)程中應(yīng)該考慮風(fēng)洞壁面效應(yīng),將來(lái)流攻角和來(lái)流馬赫數(shù)納入校正過(guò)程中。
目前大部分基于試驗(yàn)數(shù)據(jù)的流場(chǎng)參數(shù)校準(zhǔn)研究主要針對(duì)平板、圓柱繞流、槽道和翼型等經(jīng)典流動(dòng),針對(duì)壓氣機(jī)流場(chǎng)的研究相對(duì)較少。壓氣機(jī)流動(dòng)有其自身特點(diǎn),如包含逆壓力梯度、明顯的流動(dòng)分離、激波邊界層作用等流動(dòng)現(xiàn)象,且工程上壓氣機(jī)試驗(yàn)流場(chǎng)中測(cè)點(diǎn)稀疏,開(kāi)展對(duì)壓氣機(jī)流動(dòng)的增強(qiáng)預(yù)測(cè)研究,提高流場(chǎng)預(yù)測(cè)精度,對(duì)發(fā)動(dòng)機(jī)部件設(shè)計(jì)和評(píng)估有重要意義。
由此本文基于集合卡爾曼濾波數(shù)據(jù)同化方法,開(kāi)展了針對(duì)壓氣機(jī)葉柵的試驗(yàn)數(shù)據(jù)驅(qū)動(dòng)的流場(chǎng)預(yù)測(cè)研究:首先,使用試驗(yàn)函數(shù)驗(yàn)證了EnKF 算法超參數(shù)的選取準(zhǔn)則;然后,使用兩步EnKF 數(shù)據(jù)同化方法,通過(guò)校正S-A 和SST 湍流模型系數(shù)及來(lái)流邊界條件,對(duì)壓氣機(jī)葉柵進(jìn)行由試驗(yàn)測(cè)量數(shù)據(jù)驅(qū)動(dòng)的流場(chǎng)預(yù)測(cè),實(shí)現(xiàn)了對(duì)試驗(yàn)工況流場(chǎng)的重構(gòu);最后,分析了葉柵的流動(dòng)校正機(jī)理及校正系數(shù)的分布規(guī)律,比較了數(shù)據(jù)同化前后的流動(dòng)差異。
EnKF 是一種序列濾波方法,原本用于處理非定常系統(tǒng),本文研究對(duì)象是定常流動(dòng)系統(tǒng),不同于文獻(xiàn)[26]的偽非定常方法,本文對(duì)EnKF 進(jìn)行簡(jiǎn)化,使其適用于定常系統(tǒng),計(jì)算量亦有所減小。定常系統(tǒng)預(yù)測(cè)和觀測(cè)過(guò)程可以描述為
式中:θ為來(lái)流馬赫數(shù)、攻角和湍流模型系數(shù)等組成的模型系數(shù)向量;x為狀態(tài)向量,在葉柵算例中x=[θT,MaIST]T,下角標(biāo)“f”表示該狀態(tài)向量為模型預(yù)測(cè)結(jié)果,MaIS為流場(chǎng)中對(duì)應(yīng)葉柵風(fēng)洞試驗(yàn)測(cè)量位置的葉片表面等熵馬赫數(shù);y為觀測(cè)向量;ε為觀測(cè)誤差;M 為模型算子,即使用湍流模型的求解器;H 為觀測(cè)算子,從狀態(tài)向量中提取出試驗(yàn)測(cè)量位置相應(yīng)的物理量。式(1)為預(yù)測(cè)過(guò)程;式(2)為觀測(cè)過(guò)程。
對(duì)來(lái)流邊界條件和湍流模型系數(shù)依試驗(yàn)工況和系數(shù)默認(rèn)值θ0進(jìn)行浮動(dòng),浮動(dòng)過(guò)程使用拉丁超立方方法進(jìn)行采樣,生成一個(gè)包含N個(gè)樣本的集合,記為θi,i=1,2,…,N,則每個(gè)樣本均可獲得一個(gè)預(yù)測(cè)結(jié)果:
構(gòu)造樣本協(xié)方差矩陣P、誤差協(xié)方差矩陣R和卡爾曼增益矩陣K:
式中:wi為每個(gè)樣本的系統(tǒng)噪聲,一般為一組零均值正態(tài)分布的隨機(jī)數(shù),在確定其方差時(shí),需同時(shí)考慮測(cè)量噪聲和物理過(guò)程與預(yù)測(cè)模型的物理差異產(chǎn)生的偏差,比如本研究中試驗(yàn)過(guò)程為三維葉柵流動(dòng),預(yù)測(cè)模型為二維葉柵流場(chǎng)數(shù)值計(jì)算,以及湍流模型不能完全準(zhǔn)確描述湍流流動(dòng),測(cè)量偏差的方差應(yīng)該適當(dāng)大于測(cè)量不確定度對(duì)應(yīng)值。
對(duì)每個(gè)集合樣本進(jìn)行EnKF 校正,即對(duì)應(yīng)卡爾曼濾波的分析步:
式中:yexp為試驗(yàn)測(cè)量數(shù)據(jù)。
在實(shí)際應(yīng)用過(guò)程中,可以多次生成噪聲,即對(duì)式(5)~式(8)多次計(jì)算,求取平均值,可以增加算法的穩(wěn)定性,從而避免了像文獻(xiàn)[24-25]中的多次內(nèi)迭代,進(jìn)一步減小計(jì)算量。
最后獲得校正后集合成員的均值即為數(shù)據(jù)同化后的狀態(tài)變量:
從中提取出新的入口邊界條件和湍流模型系數(shù)θnew,代入求解器再進(jìn)行一次CFD 計(jì)算獲得完整的預(yù)測(cè)流場(chǎng),其對(duì)應(yīng)的狀態(tài)變量xnew即為真實(shí)狀態(tài)變量xtrue的最優(yōu)估計(jì)。
以上為單步集合卡爾曼濾波的數(shù)據(jù)同化方法,輸入數(shù)據(jù)為模型參數(shù)的初始值,輸出參數(shù)為校正后的參數(shù)值。對(duì)于本文的研究對(duì)象,EnKF已經(jīng)滿足需求,故沒(méi)有采用正則化方法。
為了驗(yàn)證算法代碼編寫的準(zhǔn)確性、測(cè)試數(shù)據(jù)同化的效果以及探索超參數(shù)選取準(zhǔn)則,使用簡(jiǎn)單數(shù)學(xué)測(cè)試函數(shù)對(duì)以上算法進(jìn)行了驗(yàn)證,測(cè)試函數(shù)包含真實(shí)函數(shù)ytrue(x)和模型函數(shù)ym(x):真實(shí)函數(shù)用以模擬真實(shí)的物理過(guò)程,其結(jié)構(gòu)相對(duì)復(fù)雜;模型函數(shù)包含若干系數(shù),表示對(duì)真實(shí)物理過(guò)程的簡(jiǎn)化和建模,結(jié)構(gòu)相對(duì)簡(jiǎn)單。借鑒文獻(xiàn)[27]中的測(cè)試函數(shù),構(gòu)造測(cè)試函數(shù)為
式中:自變量x∈[0,1],為模擬觀測(cè)過(guò)程,每隔0.1 長(zhǎng)度取一個(gè)觀測(cè)點(diǎn),共11 個(gè)觀測(cè)點(diǎn);a、b、c、d為模型系數(shù);初值取θ0=[a,b,c,d]=[6,2,12,5]。真實(shí)函數(shù)和初始模型函數(shù)如圖1 所示??梢钥闯瞿P秃瘮?shù)與真實(shí)函數(shù)形狀相近,但數(shù)值上有較大差別。取狀態(tài)變量為x=[θT,yobsT]T,其中yobs表示觀測(cè)點(diǎn)處的預(yù)測(cè)值。
圖1 測(cè)試函數(shù)圖示Fig.1 Illustration of test functions
集合卡爾曼濾波的超參數(shù)主要有樣本數(shù)量、參數(shù)浮動(dòng)范圍、噪聲給定等??紤]到數(shù)據(jù)同化過(guò)程中的各參數(shù)中包含一定的隨機(jī)成分,即便參數(shù)設(shè)置相同,每次數(shù)據(jù)同化的結(jié)果也會(huì)有微小差異,為確保評(píng)估的準(zhǔn)確性、科學(xué)性,對(duì)同種參數(shù)設(shè)置進(jìn)行多次重復(fù)計(jì)算,統(tǒng)計(jì)不同超參數(shù)設(shè)置的數(shù)據(jù)同化效果,綜合考慮計(jì)算量和統(tǒng)計(jì)效果,設(shè)定每組超參數(shù)重復(fù)計(jì)算的次數(shù)為500,獲取預(yù)測(cè)結(jié)果及中間變量的均值及標(biāo)準(zhǔn)差。當(dāng)單一超參數(shù)改變時(shí),其他參數(shù)均采取默認(rèn)值,超參數(shù)默認(rèn)值見(jiàn)表1。圖2 顯示了在默認(rèn)參數(shù)下數(shù)據(jù)同化校正結(jié)果,可以看出數(shù)據(jù)同化后的結(jié)果明顯更接近真實(shí)函數(shù),從而驗(yàn)證了算法及代碼的準(zhǔn)確性,但也可以看出,當(dāng)x<0.5 時(shí)預(yù)測(cè)效果并沒(méi)有明顯改善,即單步的校正效果有限。
表1 超參數(shù)的默認(rèn)值Table 1 Default values of hyperparameters
圖2 默認(rèn)參數(shù)下數(shù)據(jù)同化校正結(jié)果Fig.2 Result of data assimilation under default parameters
圖3 顯示了分別給定不同樣本數(shù)N、參數(shù)浮動(dòng)范圍δ和噪聲方差σnoise時(shí)的觀測(cè)點(diǎn)及全局預(yù)測(cè)均方差(Mean Square Error, MSE)的均值和標(biāo)準(zhǔn)差分布,并顯示了卡爾曼增益矩陣K和協(xié)方差矩陣P的二范數(shù)隨超參數(shù)的變化情況,折線表示均值,誤差帶表示標(biāo)準(zhǔn)差。由于測(cè)試函數(shù)部分值接近0,計(jì)算MSE 時(shí)未使用相對(duì)值。
圖3 超參數(shù)選取對(duì)數(shù)據(jù)同化效果的影響Fig.3 Effects of hyperparameters on data assimilation
從圖3 可以看出3 種超參數(shù)設(shè)置對(duì)預(yù)測(cè)結(jié)果的影響,全局MSE 與觀測(cè)點(diǎn)MSE 隨超參數(shù)變化趨勢(shì)相同,當(dāng)觀測(cè)點(diǎn)基本覆蓋了整個(gè)定義域時(shí),觀測(cè)點(diǎn)的預(yù)測(cè)效果基本可以代表全局的預(yù)測(cè)效果。
圖3(a)顯示樣本數(shù)量越小,預(yù)測(cè)效果的標(biāo)準(zhǔn)差越大,即預(yù)測(cè)結(jié)果的隨機(jī)性越明顯,當(dāng)樣本數(shù)量>50 時(shí),預(yù)測(cè)效果的標(biāo)準(zhǔn)差隨樣本數(shù)量的變化已經(jīng)不明顯,但應(yīng)注意,本測(cè)試函數(shù)只包含4 個(gè)可調(diào)系數(shù),樣本數(shù)量應(yīng)該隨系數(shù)數(shù)量適當(dāng)調(diào)整;樣本數(shù)越多,||K||2越接近于1,其標(biāo)準(zhǔn)差也越小,||P||2均值基本保持不變,標(biāo)準(zhǔn)差逐漸減小,卡爾曼濾波矩陣K表征觀測(cè)值的權(quán)重,當(dāng)樣本數(shù)增加,權(quán)重比例基本保持不變。
圖3 (b)表明參數(shù)浮動(dòng)范圍δ過(guò)小會(huì)導(dǎo)致預(yù)測(cè)的誤差增大,此時(shí)模型預(yù)測(cè)結(jié)果的||P||2的數(shù)值過(guò)小,由式(7)可知,||K||2主要受噪聲矩陣R影響,也會(huì)較小,濾波過(guò)程中觀測(cè)值修正量權(quán)重降低,輸出結(jié)果更偏向模型預(yù)測(cè)值,預(yù)測(cè)結(jié)果不能得到有效校正;當(dāng)參數(shù)浮動(dòng)范圍δ過(guò)大時(shí),預(yù)測(cè)結(jié)果的MSE 增大,會(huì)導(dǎo)致數(shù)據(jù)同化結(jié)果不穩(wěn)定,對(duì)于本算例中的非線性測(cè)試函數(shù),參數(shù)范圍的過(guò)大調(diào)整會(huì)導(dǎo)致模型函數(shù)完全偏離真實(shí)結(jié)果。
圖3(c)顯示過(guò)大或者過(guò)小的給定觀測(cè)噪聲都會(huì)造成預(yù)測(cè)誤差和隨機(jī)性增加,噪聲的給定與樣本協(xié)方差矩陣無(wú)關(guān),故右圖沒(méi)有給出||P||2;當(dāng)噪聲過(guò)小時(shí),||K||2數(shù)值較大,觀測(cè)值權(quán)重大,但由于模型本身預(yù)測(cè)能力有限,會(huì)導(dǎo)致預(yù)測(cè)效果失真,當(dāng)噪聲過(guò)大時(shí),||K||2顯著減小,導(dǎo)致觀測(cè)結(jié)果對(duì)模型修正不足,本算例中觀測(cè)結(jié)果沒(méi)有添加額外噪聲,此處的觀測(cè)噪聲主要代表模型與真實(shí)函數(shù)之間的差異,故給定噪聲時(shí)應(yīng)該充分考慮模型的預(yù)測(cè)能力,即模型與真實(shí)過(guò)程的偏差,而不僅僅是測(cè)量的不確定性。
超參數(shù)的給定本身具有一定的經(jīng)驗(yàn)性,對(duì)于復(fù)雜問(wèn)題,受計(jì)算量限制,一般難以先進(jìn)行一遍超參數(shù)選取效果的研究,故很難選取到超參數(shù)的最優(yōu)值,但根據(jù)以上論述,經(jīng)過(guò)適當(dāng)調(diào)整,一般可以獲得一組較優(yōu)的超參數(shù)。
本文研究對(duì)象為中等負(fù)荷的可控?cái)U(kuò)散擴(kuò)壓葉型MAN GHH[28],具體參數(shù)如表2 所示,葉柵參數(shù)示意如圖4 所示。
表2 葉型設(shè)計(jì)參數(shù)Table 2 Design parameters of cascade
圖4 葉柵幾何構(gòu)型Fig.4 Geometry of cascade
為減小計(jì)算量,計(jì)算過(guò)程中使用二維幾何模型,網(wǎng)格結(jié)構(gòu)為O4H,網(wǎng)格總數(shù)約為7 萬(wàn),壁面第1 層網(wǎng)格y+約為1,網(wǎng)格劃分如圖5 所示,經(jīng)驗(yàn)證滿足網(wǎng)格無(wú)關(guān)性要求。流場(chǎng)數(shù)值計(jì)算使用ANSYS Fluent,采用有限體積法,重構(gòu)格式為二階迎風(fēng),計(jì)算格式為SIMPLEC 格式。計(jì)算過(guò)程中固定進(jìn)口總壓,通過(guò)用戶自定義函數(shù)進(jìn)行二次開(kāi)發(fā),自動(dòng)反饋調(diào)節(jié)出口靜壓,以獲得指定來(lái)流馬赫數(shù)。
圖5 網(wǎng)格劃分Fig.5 Illustration of mesh
在航空航天領(lǐng)域,最常使用的湍流模型為SA 模型[29]和SST 模型[30],為測(cè)試不同湍流模型數(shù)據(jù)同化的效果并對(duì)預(yù)測(cè)結(jié)果相互驗(yàn)證,采用這2 種湍流模型分別進(jìn)行計(jì)算。湍流模型的具體方程在此不進(jìn)行贅述,僅列出2 種模型進(jìn)行校正的參數(shù)名稱和默認(rèn)值,見(jiàn)表3 和表4,參數(shù)的具體含義參見(jiàn)文獻(xiàn)[29-30]。
表3 S-A 湍流模型系數(shù)默認(rèn)值Table 3 Default parameters of S-A models
表4 SST 湍流模型系數(shù)默認(rèn)值Table 4 Default parameters of SST models
考慮到葉柵風(fēng)洞壁面效應(yīng)可能會(huì)導(dǎo)致流場(chǎng)實(shí)際來(lái)流條件與標(biāo)稱的工況值有一定差別,來(lái)流攻角和馬赫數(shù)也納入到校正系數(shù)內(nèi)。邊界條件對(duì)流場(chǎng)影響顯著,單步校正效果有限,為了確保數(shù)據(jù)同化效果并考慮到湍流模型與邊界條件參數(shù)在數(shù)據(jù)同化過(guò)程中存在耦合,不同于以往文獻(xiàn)[24-25]的單步數(shù)據(jù)同化,在本文算法實(shí)施的過(guò)程中,分2 步進(jìn)行:① 以工況標(biāo)稱的來(lái)流攻角α0和來(lái)流馬赫數(shù)Main,0作為初值,初步校正流場(chǎng)的入口邊界條件;② 以默認(rèn)湍流模型參數(shù)θt,0和第1 步的結(jié)果α1和Main,1作為初值,校正湍流模型參數(shù)和來(lái)流邊界條件,獲得新的模型系數(shù)θnew,并以此計(jì)算新的流場(chǎng)。兩步EnKF 數(shù)據(jù)同化的過(guò)程如圖6 所示。
圖6 2 步數(shù)據(jù)同化流程圖Fig.6 Flowchart of two-step data assimilation
超參數(shù)的給定參照了1.2 節(jié)的結(jié)果及文獻(xiàn)[24-25],并進(jìn)行了試算??紤]到第1 步僅校正2 個(gè)參數(shù),第1 步數(shù)據(jù)同化集合樣本數(shù)目取32,第2步根據(jù)湍流模型參數(shù)數(shù)目,綜合考慮計(jì)算穩(wěn)定性和計(jì)算量,S-A 湍流模型樣本數(shù)目設(shè)置為150,SST 模型集合樣本數(shù)目設(shè)置為200。參數(shù)浮動(dòng)范圍為7.5%,系統(tǒng)噪聲為零均值高斯噪聲,標(biāo)準(zhǔn)差為0.01。在設(shè)計(jì)馬赫數(shù)(0.62)下,對(duì)4 種攻角(-4°、0°、4°、5°)共4 種工況進(jìn)行了數(shù)據(jù)同化研究。
數(shù)據(jù)同化結(jié)果如圖7、表5 和圖8 所示。圖7顯示了2 種湍流模型數(shù)據(jù)同化過(guò)程中預(yù)測(cè)結(jié)果的相對(duì)誤差,其中第0 步為初始默認(rèn)參數(shù)的計(jì)算結(jié)果??梢钥闯?,經(jīng)過(guò)2 步校正后,2 種湍流模型的預(yù)測(cè)相對(duì)誤差都有大幅下降:S-A 模型平均降低69.6%,SST 模型平均下降67.0%。但全過(guò)程中,SST 模型的預(yù)測(cè)誤差略高于S-A 模型,且第1 步校正誤差下降并不明顯。表5 顯示了2 種湍流模型在數(shù)據(jù)同化后的來(lái)流邊界條件參數(shù)??梢钥闯觯泄r的校正馬赫數(shù)都大于工況標(biāo)稱馬赫數(shù)(0.62),校正后都在馬赫數(shù)0.64 附近浮動(dòng),不同工況略有差別,但同一工況,2 種湍流模型的校正結(jié)果非常接近,最大偏差不超過(guò)0.002。攻角也有類似結(jié)果,除4°攻角外,其他工況校正攻角與工況的標(biāo)稱攻角偏差均<1°;但4°攻角工況,偏差超過(guò)了1.5°。一般的平面葉柵試驗(yàn)過(guò)程中,攻角并非直接測(cè)量得到,而是通過(guò)轉(zhuǎn)動(dòng)安裝葉柵的轉(zhuǎn)盤到特定的角度獲得指定攻角,但在風(fēng)洞中來(lái)流會(huì)受到風(fēng)洞壁面和葉片影響,并非完全平直流動(dòng),需要進(jìn)行進(jìn)口流場(chǎng)品質(zhì)調(diào)控[31-32],故攻角的不確性度主要取決于進(jìn)口氣流角的調(diào)控效果,一般來(lái)說(shuō)給定攻角和實(shí)際攻角都會(huì)有一定差別。推測(cè)4°工況在試驗(yàn)過(guò)程存在一些偶然因素,導(dǎo)致該工況攻角偏差遠(yuǎn)大于其他工況。對(duì)于同一工況,2 種湍流模型結(jié)果非常接近,最大偏差約為0.185°。2 種湍流模型結(jié)果進(jìn)行的相互驗(yàn)證可以進(jìn)一步表明數(shù)據(jù)同化算法的可靠性。
表5 數(shù)據(jù)同化校正后的來(lái)流邊界條件Table 5 Coming boundary conditions corrected by data assimilation
圖7 數(shù)據(jù)同化過(guò)程中預(yù)測(cè)結(jié)果的相對(duì)誤差Fig.7 Relative errors of prediction during data assimilation
圖8 葉片等熵馬赫數(shù)分布(BC 為僅校正邊界條件的預(yù)測(cè), DA 為數(shù)據(jù)同化后校正所有參數(shù)的預(yù)測(cè))Fig.8 Distribution of isentropic Mach number of blade(BC represents predictions with corrected boundary conditions and DA represents predictions with all parameters corrected by data assimilation)
雖然在本文中應(yīng)用的算例工況均為中等負(fù)荷葉柵的設(shè)計(jì)馬赫數(shù),但本文發(fā)展的基于數(shù)據(jù)同化的試驗(yàn)數(shù)據(jù)驅(qū)動(dòng)的葉柵流場(chǎng)預(yù)測(cè)方法對(duì)各種負(fù)荷的壓氣機(jī)葉柵流場(chǎng)及除堵塞工況外的各工況均可以適用。因?yàn)楫?dāng)流動(dòng)堵塞時(shí),入口馬赫數(shù)保持不變,無(wú)法再通過(guò)調(diào)整入口馬赫數(shù)來(lái)校正流場(chǎng)的邊界條件。
圖8 給出了4 種工況葉片表面等熵馬赫數(shù)的分布情況。左列為完全使用默認(rèn)參數(shù)與試驗(yàn)的對(duì)比度,預(yù)測(cè)結(jié)果和試驗(yàn)結(jié)果趨勢(shì)上總體相符,但在吸力面試驗(yàn)值明顯高于預(yù)測(cè)值,而壓力面數(shù)值上比較相近。為了分別顯示來(lái)流邊界條件參數(shù)和湍流模型系數(shù)的作用,中間列顯示了僅校正了來(lái)流邊界條件(Boundary Condition,BC)與試驗(yàn)的對(duì)比結(jié)果,可以看出前50%弦長(zhǎng)區(qū)域,預(yù)測(cè)效果明顯改善但是葉片尾緣區(qū)域的預(yù)測(cè)效果幾乎沒(méi)有改變,這是由于葉片尾緣存在流動(dòng)分離,僅修正邊界條件不能提高分離預(yù)測(cè)準(zhǔn)確度,SST模型預(yù)測(cè)的葉片尾緣等熵馬赫數(shù)明顯高于試驗(yàn)結(jié)果且偏差要大于S-A 模型。右列顯示了完整參數(shù)校正后的預(yù)測(cè)結(jié)果,可以看出預(yù)測(cè)出的等熵馬赫數(shù)分布與試驗(yàn)測(cè)量結(jié)果幾乎完全吻合,這說(shuō)明湍流模型參數(shù)的校正顯著提高了尾緣分離區(qū)的預(yù)測(cè)效果。
另一方面,圖8 可以解釋SST 模型第一步校正時(shí),誤差下降很小的原因。來(lái)流邊界條件校正主要提高了葉片前半部分的預(yù)測(cè)效果,但原始參數(shù)的SST 模型尾緣等熵馬赫數(shù)本身存在較大的過(guò)預(yù)測(cè),在校正來(lái)流條件時(shí)會(huì)使尾緣區(qū)域的過(guò)度預(yù)測(cè)更為明顯,為降低總體誤差,邊界條件調(diào)整幅值較小,即湍流模型系數(shù)與邊界條件參數(shù)在校正過(guò)程中存在耦合,導(dǎo)致SST 模型單純校正邊界條件效果不佳。而S-A 模型尾緣的預(yù)測(cè)偏差較小,對(duì)邊界條件校正影響小,即二者耦合程度低,故僅校正邊界條件預(yù)測(cè)誤差就有明顯減小。
表6 和表7 分別給出了2 種湍流模型各工況校正后的系數(shù),括號(hào)中的百分?jǐn)?shù)表示相對(duì)默認(rèn)值變化的百分比??梢钥闯霾糠窒禂?shù)有明顯的趨勢(shì)性變化,總體上相對(duì)默認(rèn)值有明顯的增大或者減小,如S-A 模型的Cb1、σv,SST 模型的、a1、σk1、σω2、κ等。文獻(xiàn)[17-18]為提高S-A 模型對(duì)分離流的預(yù)測(cè)能力,建議增大Cb1參數(shù),文獻(xiàn)[26]建議將SST 模型中的參數(shù)a1設(shè)置為1.0以提高分離流和逆壓力梯度流動(dòng)的預(yù)測(cè)能力。在調(diào)整方向上相關(guān)文獻(xiàn)與本研究校正結(jié)果相符,但由于本研究同時(shí)校正了多個(gè)參數(shù),而上述文獻(xiàn)僅調(diào)整了一個(gè)參數(shù),故數(shù)值上差別較大。部分參數(shù)如:S-A 模型Cb2、κ,SST 模型βi2、σω1變化幅度較小,其他參數(shù)變化幅度有明顯的工況相關(guān)性。值得注意的是,2 種模型中都有von Kármán 常數(shù)κ,但變化規(guī)律明顯不同,在S-A模型中基本保持不變而在SST 模型中有明顯的增加,原因在于2 種模型對(duì)該常數(shù)的使用不同。
表6 S-A 湍流模型校正系數(shù)Table 6 Corrected coefficients of S-A model
表7 SST 湍流模型校正系數(shù)Table 7 Corrected coefficients of SST model
校準(zhǔn)后湍流模型參數(shù)可以有效減小湍流模型對(duì)流動(dòng)分離的過(guò)度預(yù)測(cè),雖然在其他葉柵上的效果仍有待測(cè)試,但由于參數(shù)調(diào)整幅度不大,參數(shù)的調(diào)整并不會(huì)帶來(lái)明顯負(fù)面效果,當(dāng)其他壓氣機(jī)葉柵流場(chǎng)存在分離區(qū)過(guò)度預(yù)測(cè)時(shí),可以參照本文中近似工況的校準(zhǔn)參數(shù)進(jìn)行設(shè)置以降低分離區(qū)過(guò)度預(yù)測(cè),提高預(yù)測(cè)精度。
圖9 顯示了2 種湍流模型在0°和5°攻角預(yù)測(cè)結(jié)果的馬赫數(shù)云圖和流線圖,左列為僅校正邊界條件的預(yù)測(cè)結(jié)果,右列為校正邊界條件和湍流模型參數(shù)的預(yù)測(cè)結(jié)果。從流線圖可以看出,使用默認(rèn)的湍流參數(shù),吸力面尾緣的分離泡尺寸總體大于校正后的結(jié)果;馬赫數(shù)云圖顯示,湍流模型校正后的尾跡范圍明顯比校正前小。2 種湍流模型默認(rèn)參數(shù)條件下,預(yù)測(cè)分離區(qū)尺寸和起始分離位置有明顯差別,但數(shù)據(jù)同化后2 種模型預(yù)測(cè)結(jié)果非常接近。
圖9 2 種湍流模型數(shù)據(jù)同化過(guò)程的馬赫數(shù)云圖及流線圖Fig.9 Mach contours and streamlines of two models during data assimilation
圖9 也解釋了圖8 等熵馬赫數(shù)變化的原因:默認(rèn)參數(shù)的湍流模型預(yù)測(cè)的分離泡過(guò)大,導(dǎo)致葉片尾緣流道變窄,該流場(chǎng)為亞聲速流場(chǎng),流速及馬赫數(shù)增加,靜壓下降即等熵馬赫數(shù)增加,大于試驗(yàn)測(cè)量結(jié)果;參數(shù)校正后預(yù)測(cè)的分離泡尺寸減小,等熵馬赫數(shù)與試驗(yàn)測(cè)量結(jié)果更加接近。
為了進(jìn)一步分析吸力面尾緣分離區(qū)起始位置及葉片所受剪切力的情況,圖10 顯示了4 種工況葉片壓力面和吸力面剪切力系數(shù)Cf的分布情況。Cf的正方向?yàn)榱鲃?dòng)方向,當(dāng)Cf為負(fù)值時(shí)即發(fā)生了流動(dòng)分離??梢钥闯?,S-A 模型參數(shù)校正前后,葉片前50%的區(qū)域剪切力系數(shù)幾乎不變,變化主要發(fā)生在尾緣區(qū)域,分離點(diǎn)附近變化最為明顯。而SST 模型系數(shù)校正前后,整個(gè)葉片表面剪切力系數(shù)都有明顯變化,吸力面變化幅度大于壓力面。2 種模型在校正后都出現(xiàn)了分離起始點(diǎn)推遲,SST 模型推遲更大,分離點(diǎn)預(yù)測(cè)位置基本一致,在吸力面前約40% 弦長(zhǎng)、壓力面前20% 弦長(zhǎng),2 種模型校正后剪切力系數(shù)預(yù)測(cè)有一定差異,其余區(qū)域基本一致。剪切力系數(shù)在數(shù)據(jù)同化過(guò)程中沒(méi)有被直接校正,而是受到湍流模型系數(shù)的影響產(chǎn)生變化,SST 模型參數(shù)校正后前緣區(qū)域Cf的變化應(yīng)該是由于SST 模型邊界層對(duì)模型系數(shù)比較敏感。
圖10 2%~98%弦長(zhǎng)區(qū)域葉片表面剪切力系數(shù)分布Fig.10 Distributions of shear coefficient of blade surface along 2% to 98% chord
本研究使用的試驗(yàn)數(shù)據(jù)僅來(lái)自于葉片表面的靜壓孔,有必要分析數(shù)據(jù)同化后葉片下游流動(dòng)參數(shù)的變化情況。圖11 和圖12 顯示了葉片尾緣下游0.34 倍弦長(zhǎng)處氣流角和總壓損失系數(shù)ω沿截距方向的分布,xt表示距葉片尾緣的流向距離??梢钥闯?,經(jīng)數(shù)據(jù)同化校正后氣流角和總壓損失系數(shù)明顯減小,2 種湍流模型預(yù)測(cè)結(jié)果非常接近,葉片尾跡區(qū)的位置和參數(shù)變化幅值基本一致,尾跡范圍小于原始模型預(yù)測(cè)結(jié)果。葉片下游的氣流角和損失系數(shù)受吸力面分離區(qū)影響很大,校正后,分離區(qū)尺寸減小,從而降低了氣流角和總壓損失。
圖11 葉片下游氣流角(xt/C=0.34)Fig.11 Flow angles at blade downstream (xt/C=0.34)
圖12 葉片下游總壓損失系數(shù)Fig.12 Total pressure loss coefficients at blade downstream
為進(jìn)一步探究更遠(yuǎn)處下游流場(chǎng)情況,圖13顯示了2 種校正后的湍流模型在0°和5°攻角葉片尾緣0.34、1.34、2.34 倍弦長(zhǎng)處的總壓損失分布情況??梢钥闯觯幢阍谳^遠(yuǎn)的下游區(qū)域,2 種湍流模型預(yù)測(cè)結(jié)果仍然基本一致。這說(shuō)明使用集合卡爾曼濾波方法對(duì)壓氣機(jī)葉柵進(jìn)行試驗(yàn)數(shù)據(jù)驅(qū)動(dòng)的流場(chǎng)預(yù)測(cè)時(shí),全流場(chǎng)中大多數(shù)與葉柵性能直接相關(guān)的物理量的預(yù)測(cè)結(jié)果對(duì)湍流模型依賴性較弱。
圖13 葉片不同下游位置總壓損失系數(shù)分布Fig.13 Distributions of total pressure loss coefficient at several locations of blade downstream
本文使用集合卡爾曼濾波方法,使用S-A 和SST 這2 種湍流模型對(duì)MAN GHH 葉柵進(jìn)行了試驗(yàn)數(shù)據(jù)驅(qū)動(dòng)的流場(chǎng)預(yù)測(cè)。結(jié)論如下:
1)使用集合卡爾曼濾波方法進(jìn)行參數(shù)校正時(shí),應(yīng)該合理選擇超參數(shù),其中觀測(cè)噪聲的給定應(yīng)該同時(shí)包含測(cè)量不確定度和模型與真實(shí)物理過(guò)程之間的差異。
2)通過(guò)集合卡爾曼濾波的數(shù)據(jù)同化方法校正湍流模型系數(shù)和來(lái)流邊界條件,可以有效提高壓氣機(jī)葉柵的預(yù)測(cè)精度,可以使預(yù)測(cè)誤差降低將近70%,由于試驗(yàn)過(guò)程中風(fēng)洞壁面效應(yīng)的作用,校正來(lái)流邊界條件是必要的,校正后的湍流模型參數(shù)部分呈現(xiàn)一定共同規(guī)律性。
3)使用2 種湍流模型分別進(jìn)行兩步數(shù)據(jù)同化參數(shù)校正后,相比原始流場(chǎng),參數(shù)校正后的流場(chǎng)吸力面尾緣分離區(qū)尺寸減小,分離起始點(diǎn)延遲,其入口邊界條件及流場(chǎng)參數(shù)基本相同,即參數(shù)校正后的流場(chǎng)對(duì)湍流模型依賴性較弱。