朱玉杰,琚亞平,張楚華

(西安交通大學能源與動力工程學院,710049,西安)

提高壓氣機葉柵氣動性能、擴大其穩(wěn)定高效運行范圍,始終是流體機械及工程、空氣動力學等學科,以及風機、航空等工業(yè)領(lǐng)域的關(guān)注重點。發(fā)展先進的葉柵葉型設(shè)計優(yōu)化方法是其重要基礎(chǔ),一般而言,葉柵葉型優(yōu)化設(shè)計可以分為正問題設(shè)計和反問題設(shè)計兩大類[1]。

正問題優(yōu)化設(shè)計方法是預先給定葉片幾何尺寸和來流條件,根據(jù)流場分析或性能試驗結(jié)果修正原始幾何尺寸,直到獲得性能最優(yōu)的葉型幾何,這種方法可以直接得到對應性能的葉型,但一般耗時較大。盡管近似模型和隨機性優(yōu)化方法的發(fā)展極大提高了正問題優(yōu)化求解的效率[2-4],但其大多局限于處理維數(shù)較低的壓氣機葉柵氣動優(yōu)化問題。當維數(shù)較高時,所需計算成本會急劇增大,面臨“維數(shù)災難”問題[5]。而因為伴隨變量的引入,基于最優(yōu)控制理論的伴隨方法,具有計算量與設(shè)計變量數(shù)目基本無關(guān)的優(yōu)點[6-10],為壓氣機葉柵高維氣動的正問題優(yōu)化設(shè)計提供了新的有效途徑。

反問題優(yōu)化設(shè)計通過預先給定葉片表面馬赫數(shù)或者壓力的分布規(guī)律等氣動控制參數(shù),定制出滿足特定氣動需求的先進葉型。反問題優(yōu)化求解可以實現(xiàn)葉柵葉型精細化設(shè)計,精確控制葉柵內(nèi)部和葉片表面氣動參數(shù)分布。通常的反問題設(shè)計是采用虛位移概念將目標控制參數(shù)與邊界條件聯(lián)系起來修改葉型[11-14],這種方法的邊界更新方式一般比較復雜,且易因更新過程中的過修改引起計算的不穩(wěn)定。相比而言,通過引入伴隨變量將反問題設(shè)計轉(zhuǎn)化成優(yōu)化問題的伴隨方法,因為設(shè)計變量的更新與目標函數(shù)的梯度方向相關(guān)而使得計算收斂更加穩(wěn)定。然而,由于目標氣動參數(shù)的給定往往與設(shè)計者的經(jīng)驗相關(guān),現(xiàn)有研究中,關(guān)于目標氣動參數(shù)的給定形式仍較為模糊,沒有簡單的規(guī)律可循[15],這從一定程度上限制了反問題的優(yōu)化求解方法在流體機械設(shè)計中的應用發(fā)展。

因此,本文基于連續(xù)型伴隨方法理論,分別建立了集葉片幾何參數(shù)化、網(wǎng)格生成、流場求解、伴隨場求解與優(yōu)化求解于一體的葉柵葉型正、反問題設(shè)計優(yōu)化求解方法?；趯θ~柵正問題優(yōu)化求解所得葉片表面壓力分布的分析,給出在更大進口氣流角情況下的葉柵反問題優(yōu)化求解所需目標壓力分布,旨在對葉柵內(nèi)部流動分離進行進一步控制,進而提高葉柵氣動性能。

1 流動方程及數(shù)值方法

在二維笛卡爾坐標系下,非定?？蓧嚎sNavier-Stokes方程可寫成如下形式

(1)

式中:Q、fi、fvi分別為守恒變量、對流通量、黏性通量;xi為笛卡爾坐標分量,i=1,2,t為時間。為使方程組封閉,采用S-A方程湍流模型。數(shù)值模擬通過已驗證有效的自編程序[7]來實現(xiàn),其中時間項采用五步Runge-Kutta法離散,對流項采用二階中心格式結(jié)合人工黏性方法離散,黏性項采用二階中心格式離散,同時采用多重網(wǎng)格、當?shù)貢r間步長和隱式殘差光順技術(shù)加速收斂。

2 正問題設(shè)計伴隨方程

根據(jù)伴隨方程推導過程及最終形式的不同,伴隨方法分為離散型伴隨方法和連續(xù)型伴隨方法。考慮到連續(xù)型伴隨方法具有物理意義明確、表達形式簡潔直觀等特點[16],本文基于該方法,建立與平面葉柵葉型正問題和反問題設(shè)計優(yōu)化目標相對應的伴隨方程及其邊界條件。

在平面葉柵葉型正問題優(yōu)化求解中,以出口總壓損失最小作為優(yōu)化目標,即

(2)

式中:p0為給定進口總壓;pt為總壓;ΔS為葉柵出口面積;ds為面積微元。在優(yōu)化循環(huán)中,忽略湍流黏性系數(shù)的改變。引入伴隨變量ψ={ψ1,ψ2,ψ3,ψ4}T,在流動區(qū)域內(nèi)各點計算其與流動方程變分的內(nèi)積,將結(jié)果代入目標函數(shù)在設(shè)計變量發(fā)生微小改變時的變分表達式,可得

(3)

式中:D為計算區(qū)域;dv為體積微元。為將式(3)中流動變量變分有關(guān)項和設(shè)計變量變分有關(guān)項分離,引入流動變量空間導數(shù)的變分關(guān)系式[7],并結(jié)合高斯積分公式,式(3)中右端第4項轉(zhuǎn)換為

(4)

(5)

式(4)中右端第4項轉(zhuǎn)換為

(6)

在計算區(qū)域D內(nèi),將式(5)、(6)代入式(3)即可實現(xiàn)目標函數(shù)變分表達式中流動變量變分項和幾何變量變分項的分離。為消除目標函數(shù)變分對流動變量變分項的依賴,并引入虛擬時間導數(shù)項以便數(shù)值求解,可得如下伴隨方程

(7)

同時,由式(3)(4)(5)可知,在計算區(qū)域邊界,影響目標函數(shù)變分的流動變量變分項為

(8)

由式(8)可得伴隨方程的進出口邊界條件

(9)

在固壁邊界,采用絕對無滑移絕熱邊界條件,則

(10)

由式(2)可得伴隨方程的固壁邊界條件

(11)

式(7)(9)(11)加之周期邊界共同組成了以出口總壓損失最小為目標的葉柵葉型氣動優(yōu)化伴隨方程及其邊界條件。在此基礎(chǔ)上,考慮進出口網(wǎng)格節(jié)點固定以及出口面積不變的特點,即δxi=δΔS=0,目標函數(shù)變分的最終表達式為

(12)

3 反問題設(shè)計伴隨方程

在平面葉柵葉型反問題優(yōu)化求解中,以葉片表面與目標壓力分布差值最小為目標函數(shù),即

(13)

式中:ptar為目標壓力。類似上述推導[7],可得伴隨方程

(14)

進出口及固壁邊界條件

(15)

目標函數(shù)變分的最終表達式

(16)

4 葉柵葉型設(shè)計優(yōu)化算例

基于流場和伴隨場數(shù)值求解方法,完成流場和伴隨場數(shù)值求解并獲取目標函數(shù)梯度,最后采用BFGS擬牛頓法、一維尋優(yōu)法確定尋優(yōu)方向和尋優(yōu)步長[7]。

4.1 正問題設(shè)計算例

選取NACA 65-(12)10葉型為原始葉型,流場計算條件是:進口總溫為300 K,進口總壓為101 325 Pa,出口靜壓為92 563 Pa,葉片安裝角為25.9°,進口氣流角為42.0°,弦長l為0.127 m,相對柵距為1.0。對葉柵葉型進行正問題設(shè)計優(yōu)化求解,以期減少葉柵總壓損失。

分別采用6次Bezier曲線對原始葉型壓力面和吸力面進行擬合,選取Bezier曲線控制頂點的10個縱坐標為設(shè)計變量。在優(yōu)化過程中,為避免葉片幾何型線不合理現(xiàn)象并約束葉片厚度,控制吸力面和壓力面的第2、第6個Bezier控制點的變化范圍在±5%以內(nèi),其他控制點的變化范圍在±15%以內(nèi);為保證葉柵做功能力以及葉型幾何在前緣的一階連續(xù),固定葉片前尾緣不動,吸力面的第2個和第6個Bezier控制點分別限制在S1S2及S6S7兩條線上移動,壓力面的第2個和第6個Bezier控制點分別限制在P1P2及P6P7兩條線上移動,如圖1所示。經(jīng)過9個優(yōu)化迭代步之后,計算收斂。為評估正問題優(yōu)化設(shè)計結(jié)果,分別通過表1和圖2給出了優(yōu)化前后葉柵總壓損失系數(shù)和葉柵流場葉片尾緣速度矢量的局部放大圖。由表1可知,葉柵總壓損失系數(shù)下降至0.002 32,相比于優(yōu)化前的0.002 46,下降約5.69%。由圖2可知,優(yōu)化前葉柵流道在葉片尾緣附近存在著明顯的流動分離,而優(yōu)化后流動分離得到有效抑制,流動損失降低,葉柵氣動性能得到改善。

圖1 Bezier曲線及其控制點

葉柵總壓損失系數(shù)靜壓比原始葉型0.002 461.069 6正問題解0.002 321.071 0變化量/%-5.690.13

(a)優(yōu)化前

(b)優(yōu)化后圖2 正問題設(shè)計的葉片尾緣速度矢量分布

圖3 正問題設(shè)計優(yōu)化前后葉型對比

為進一步了解葉型對葉柵氣動性能的影響,對比了正問題設(shè)計伴隨方法優(yōu)化前、后的葉型,如圖3所示。由圖3可知,優(yōu)化所得葉型尤其是吸力面葉型曲線較初始葉型沿弦線方向變化平緩,曲率減小,對應葉型彎度變小。由氣動角度看,這一幾何特征有利于降低葉片中后部逆壓梯度,減小流動分離,與優(yōu)化結(jié)果一致,如圖4所示。本文以葉柵出口總壓損失最小為目標的正問題設(shè)計優(yōu)化方法,可以有效降低葉柵流道損失,提高葉柵效率。

圖4 正問題設(shè)計的葉片表面壓力

4.2 反問題設(shè)計算例

文獻[7]中驗證了本文反問題優(yōu)化求解方法可以有效獲得表面壓力分布逼近目標壓力分布的葉型。本算例擬開展更大進口氣流角條件下的葉柵葉型反問題優(yōu)化求解設(shè)計,初始葉型采用NACA 65-(12)10翼型,改變流場進口氣流角為43.5°,其他計算條件均與4.1節(jié)算例相同。初始葉柵葉片尾緣速度矢量圖如圖5所示,可知在葉片尾緣附近存在明顯的流動分離現(xiàn)象。

本算例擬從減弱流動分離的角度出發(fā),通過給定合適的葉片表面目標壓力分布,實現(xiàn)對葉柵葉型反問題設(shè)計的優(yōu)化求解。由載荷分布角度可知,優(yōu)化后為前加載,優(yōu)化后的葉片載荷較優(yōu)化前在前緣附近有明顯的增加,在葉片20%～80%弦長范圍內(nèi)有所減小;由壓力梯度的角度可知,優(yōu)化后從葉片中前部開始,逆壓梯度減小,具體表現(xiàn)為優(yōu)化前吸力面壓力在2%～30%弦長范圍內(nèi),上升趨勢平緩,在30%～80%弦長范圍內(nèi),上升趨勢陡峭,而優(yōu)化后,在30%～80%弦長范圍內(nèi)吸力面壓力上升曲線平緩?；谏鲜隽鲌龇治?為減小43.5°進口氣流角下的流動分離,根據(jù)分析所得正問題優(yōu)化后相對于優(yōu)化前的葉片表面壓力分布變化趨勢,對反問題初始葉型表面的壓力分布做相應修正,從而給出該進口氣流角條件下的理想目標壓力分布,如圖6所示。

經(jīng)過9步迭代,計算收斂。為了便于對比分析,圖6同時給出了反問題設(shè)計優(yōu)化求解后的壓力分布。由圖6可知:反問題設(shè)計優(yōu)化求解所得壓力分布很好地逼近了目標壓力分布;相比初始壓力分布,優(yōu)化后壓力分布為前加載,在吸力面30%～80%弦長范圍內(nèi)逆壓梯度減小,壓力上升較為平緩。

圖5 43.5°進口氣流角下葉片尾緣速度矢量分布

圖6 反問題設(shè)計的葉片表面壓力

圖7 反問題設(shè)計的葉片尾緣速度矢量分布

圖8 反問題設(shè)計優(yōu)化前后葉型對比

為評估反問題優(yōu)化設(shè)計結(jié)果,驗證所給目標壓力分布有效性,優(yōu)化后葉柵流場葉片尾緣速度矢量的局部放大圖如圖7所示,優(yōu)化前后葉柵總壓損失系數(shù)和靜壓比對比如表2所示。相比于圖5，圖7所示的葉片尾緣附近的流動分離得到明顯減弱。由表2可知,優(yōu)化后靜壓比基本不變,而總壓損失系數(shù)減小了4.50%,葉柵氣動性能得到提升,高效穩(wěn)定運行范圍有所拓寬,從而驗證了所給定的目標壓力分布是合理有效的。為進一步了解葉片表面壓力分布對葉型的控制,對比了優(yōu)化前、后的葉型幾何,如圖8所示。由圖8可知,優(yōu)化后葉片吸力面型線較優(yōu)化前更為平緩,變化趨勢與圖3中正問題優(yōu)化求解結(jié)果類似,說明所給定的目標壓力分布能夠有效控制葉型。

5 結(jié) 論

(1)發(fā)展了針對葉柵葉型正反問題設(shè)計的伴隨優(yōu)化方法和數(shù)值模型,實現(xiàn)了對葉柵葉型正問題和反問題設(shè)計的優(yōu)化求解。新設(shè)計出的葉型葉柵氣動性能均優(yōu)于原始葉柵,驗證了該優(yōu)化求解方法和模型的有效性,展示了該方法和模型在提高葉柵葉型氣動性能方面的優(yōu)化能力。

(2)對葉柵葉型正問題設(shè)計優(yōu)化求解所得葉片表面壓力分布進行分析,可以給出較優(yōu)的葉柵葉型反問題設(shè)計目標壓力分布,較好地解決反問題中難以給定氣動控制參數(shù)的問題。

(3)對NACA 65-(12)10葉型的優(yōu)化結(jié)果表明,在所研究的進氣條件下,葉柵葉型吸力面曲率減小、葉片彎度降低的幾何特征以及葉片前加載方案、吸力面中后部逆壓梯度減小的氣動特征,有利于抑制或消除葉柵葉片尾緣流動分離現(xiàn)象,提高葉柵葉型氣動性能。