基于氣動(dòng)反問(wèn)題的擴(kuò)壓葉柵優(yōu)化設(shè)計(jì)

姜俊鋒，陳 榴，胡 磊，戴 韌

1 前言

軸流風(fēng)扇和壓氣機(jī)葉型的基本構(gòu)造方法是在標(biāo)準(zhǔn)的中弧線上迭加標(biāo)準(zhǔn)厚度分布，例如NACA65系列［1］。這樣構(gòu)造出來(lái)的葉柵具有豐富的氣動(dòng)性能試驗(yàn)數(shù)據(jù)和詳細(xì)的設(shè)計(jì)應(yīng)用準(zhǔn)則［2,3］。選擇合適的葉柵稠度和葉片安裝角，即可滿足風(fēng)扇設(shè)計(jì)速度三角形的要求。隨著葉柵氣動(dòng)優(yōu)化方法的發(fā)展，風(fēng)扇葉型設(shè)計(jì)突破了標(biāo)準(zhǔn)葉型的束縛，依據(jù)具體工況要求，采用某種優(yōu)化算法，尋找與某個(gè)優(yōu)化目標(biāo)的最佳值所對(duì)應(yīng)“定制”葉型，提高了葉型的氣動(dòng)性能?？缫羲倭鲃?dòng)可控?cái)U(kuò)散葉型和低速流動(dòng)的自然層流葉型是2個(gè)重要成果［4～6］。

葉柵氣動(dòng)優(yōu)化方法需要大量的流場(chǎng)計(jì)算，雖然魯棒性好，但是耗費(fèi)的機(jī)時(shí)令人難以接受［7］。尤其重要的是，優(yōu)化結(jié)果的篩選完全取決于優(yōu)化算法的能力，往往陷入所謂的“局部”優(yōu)化，而且得到的優(yōu)化解往往缺乏合理的物理解釋。另一方面，經(jīng)過(guò)長(zhǎng)期的氣動(dòng)力學(xué)研究，我們已經(jīng)掌握了具有較好氣動(dòng)性能葉型的流動(dòng)分布特征。如果能預(yù)先給出某個(gè)氣動(dòng)目標(biāo)，再尋找對(duì)應(yīng)的葉型并做優(yōu)化，那么就可以節(jié)省大量的計(jì)算費(fèi)用，而且得到的解也具有物理上的合理性。

與葉柵流動(dòng)的正問(wèn)題相比，反問(wèn)題的理論還不夠完善。葉柵氣動(dòng)反問(wèn)題的第一個(gè)困難是如何建立氣動(dòng)目標(biāo)與葉型幾何之間的對(duì)應(yīng)關(guān)系，因?yàn)闆](méi)有流動(dòng)方程是以葉型為因變量，而是以速度或壓力分布為自變量而建立的。一個(gè)直觀可行的方法是從某個(gè)相近的初始葉型出發(fā)，比較其流動(dòng)與設(shè)計(jì)目標(biāo)的差距，逐步修改葉型逼近設(shè)計(jì)目標(biāo)［8～10］。其中 Kasra 采用虛擬速度法［8］，Dulikravich等人采用傅里葉級(jí)數(shù)［9］，姚征采用渦量修正的葉型修改方法［10］。戴韌將旋成面葉柵氣動(dòng)正命題的微分-積分解法與葉型修正的“噴氣模型”結(jié)合［11］，獲得反問(wèn)題的解，采用正反問(wèn)題結(jié)合的方法實(shí)現(xiàn)了葉片反設(shè)計(jì)。

葉柵反問(wèn)題求解的第二個(gè)難題是預(yù)定氣動(dòng)設(shè)計(jì)目標(biāo)的可行性是沒(méi)有保證的，因此，完全精確的反問(wèn)題解是不存在的。為了解決這個(gè)難題，以各種最優(yōu)化的方法，逼近氣動(dòng)設(shè)計(jì)目標(biāo)，獲得某種誤差范數(shù)的最小值所對(duì)應(yīng)的葉型，即為葉柵反問(wèn)題的解。Pehlivanoglu等采用遺傳算法［12］，Obayashi等采用進(jìn)化算法［13］，豐鎮(zhèn)平基于控制理論，采用伴隨梯度算法，都取得了不錯(cuò)的成果［14］。

本文利用邊界層積分解的Head方法建立葉型表面壓力系數(shù)分布與邊界層積分厚度的關(guān)系［15］，利用遺傳算法，獲得近似最優(yōu)的壓力系數(shù)分布。采用反問(wèn)題的“彈性模型”求解方法，逼近目標(biāo)壓力系數(shù)，并結(jié)合試驗(yàn)設(shè)計(jì)方法優(yōu)選反問(wèn)題的目標(biāo)壓力系數(shù)分布，從而建立基于反問(wèn)題的優(yōu)化設(shè)計(jì)方法。將此方法應(yīng)用到具體實(shí)例進(jìn)行驗(yàn)證，獲得符合氣動(dòng)設(shè)計(jì)目標(biāo)的葉型，改善葉柵流場(chǎng)，降低葉型的損失，提高氣動(dòng)性能。

2 葉柵氣動(dòng)反問(wèn)題設(shè)計(jì)方法

葉柵氣動(dòng)反問(wèn)題設(shè)計(jì)方法的計(jì)算流程如圖1所示，其中有3個(gè)重要步驟：確定目標(biāo)壓力系數(shù)分布、初始葉型的流動(dòng)分析和葉型的修正。葉柵反問(wèn)題的氣動(dòng)目標(biāo)與單機(jī)翼氣動(dòng)反問(wèn)題是不同的，因?yàn)槿~柵進(jìn)出口速度三角形限定了葉片繞流環(huán)量的大小，所給定的氣動(dòng)目標(biāo)需要滿足這個(gè)環(huán)量的約束，因此滿足氣動(dòng)目標(biāo)的解是有限的。葉型的修正有多種方法，但是修正后葉型的二階連續(xù)是至關(guān)重要的，因?yàn)樾途€曲率的間斷點(diǎn)是流動(dòng)發(fā)生階躍的位置，在擴(kuò)壓流動(dòng)中，會(huì)引起邊界層的快速增長(zhǎng)甚至是分離。

圖1 反問(wèn)題設(shè)計(jì)流程

2.1 目標(biāo)壓力分布的近似優(yōu)化

確定一個(gè)近似最優(yōu)的目標(biāo)壓力系數(shù)分布是反問(wèn)題設(shè)計(jì)中很重要的一步。本文根據(jù)Head邊界層計(jì)算方法，建立邊界層外緣速度ue與邊界層厚度以及壁面摩擦系數(shù)之間的關(guān)系。

二維不可壓的流體動(dòng)量積分方程為［15］：

引入Head裹入方程，與式（1）聯(lián)立求解邊界層動(dòng)量厚度θ沿流向的增長(zhǎng)。

其中 H=δ＊/θ,H1=（δ-δ＊）/θ

式中 H1——形狀因子

對(duì)于給定的速度分布ue，式（1）、（2）中包含3個(gè)未知數(shù)θ，Cf和H，還需補(bǔ)充若干關(guān)系式才能封閉求解［15］。設(shè)想沿葉型表面的速度分布ue（x）表述為由（n+1）點(diǎn)控制的n次參數(shù)化多項(xiàng)式：

其中Pi（t）是多項(xiàng)式基函數(shù)，ue,i為控制點(diǎn)。利用遺傳算法優(yōu)化{ue,i}，得出一個(gè)近似最優(yōu)的速度分布ue（x）。將此速度分布按照伯努利積分換算得出對(duì)應(yīng)的最優(yōu)壓力系數(shù)分布。

選取3種典型的速度分布，如圖2所示。經(jīng)過(guò)對(duì)比發(fā)現(xiàn)，內(nèi)凹的速度分布具有較小的動(dòng)量損失厚度，壁面摩擦系數(shù)也較小，如圖3，4所示。因此以?xún)?nèi)凹速度分布為初始速度分布進(jìn)行優(yōu)化，利用貝濟(jì)埃（Bezier）曲線對(duì)速度分布進(jìn)行參數(shù)化。

圖2 3種典型不同的速度分布

圖3 動(dòng)量損失厚度分布

圖4 壁面摩擦系數(shù)分布

優(yōu)化目標(biāo)：

Min：F=θout

Subject to：Cf≤ Cf.baseline

經(jīng)過(guò)優(yōu)化得到近似最優(yōu)的速度分布，并按照伯努利方程換算得出對(duì)應(yīng)的壓力系數(shù)分布如圖5所示。最優(yōu)壓力系數(shù)分布在x/Cax=0.5之前下降速度較快，曲線的切向與水平面的夾角從15°到55°之間進(jìn)行變化，x/Cax=0.5之后壓力系數(shù)變化較小，基本為水平線且?guī)в幸欢ǖ南掳肌?/p>

圖5 Head方法最優(yōu)壓力系數(shù)分布

此壓力系數(shù)分布是利用Head方法求解出的一個(gè)最優(yōu)的壓力系數(shù)分布，但是由于Head方法計(jì)算的時(shí)候沒(méi)有考慮曲率的影響，以及葉柵流動(dòng)之間的相互作用等因素，不能直接應(yīng)用到本文的反問(wèn)題設(shè)計(jì)中，還需要結(jié)合葉柵流動(dòng)的特點(diǎn)進(jìn)行適當(dāng)?shù)男薷?。因此以Head求解出的最優(yōu)壓力系數(shù)分布為初始?jí)毫ο禂?shù)，并利用貝濟(jì)埃（Bezier）曲線進(jìn)行參數(shù)化，如圖6所示，通過(guò)控制點(diǎn)使壓力系數(shù)在最優(yōu)壓力系數(shù)附近進(jìn)行變動(dòng)，并將此壓力系數(shù)分布作為目標(biāo)壓力系數(shù)分布?？刂泣c(diǎn)僅改變X3，X4和 X5，其余控制點(diǎn)保持不變。

圖6 壓力系數(shù)分布參數(shù)化

2.2 葉型氣動(dòng)反問(wèn)題的解法

葉型氣動(dòng)反問(wèn)題求解能否成功的關(guān)鍵在于，如何將葉型表面氣動(dòng)參數(shù)與氣動(dòng)目標(biāo)的差距，轉(zhuǎn)換為葉型修改的“動(dòng)力”。本文將葉片表面視為彈性膜，當(dāng)初始葉型表面壓力與目標(biāo)值不一致時(shí)，在兩者之間的壓差作用下，葉片表面發(fā)生局部“變形”，變形量的大小與當(dāng)?shù)貕翰畛烧?。?jù)此建立葉型修正模型。

式中 ΔS——葉型表面法向位移

δ——葉型的當(dāng)?shù)睾穸?/p>

k——變形率系數(shù)，相當(dāng)與彈性模型。

ΔCp—— 計(jì)算壓力系數(shù)與目標(biāo)壓力系數(shù)的差值

Cp——當(dāng)?shù)貕毫ο禂?shù)

在葉型表面選取m個(gè)反問(wèn)題目標(biāo)控制點(diǎn)，覆蓋部分葉型，余下葉型自由變動(dòng)，以滿足葉柵總環(huán)量的約束。m點(diǎn)葉型的變化后，需要與其余未修正的葉型光滑連接，并保證二階以上的連續(xù)。本文采用附加約束的最小二乘法進(jìn)行擬合得到新葉型［16］。

采用目標(biāo)壓力的收斂殘差Resp和葉型幾何殘差Resg考核氣動(dòng)反問(wèn)題的收斂性，兩者分別定義為：

式中 cp——計(jì)算壓力系數(shù)

cp*——目標(biāo)壓力系數(shù)

Si——原始葉型控制點(diǎn)

Si*——反設(shè)計(jì)葉型控制點(diǎn)

2.3 基于反問(wèn)題的優(yōu)化設(shè)計(jì)

應(yīng)用Head積分，獲得葉片吸力面一個(gè)最優(yōu)壓力系數(shù)分布，作為葉型優(yōu)化的初始反問(wèn)題目標(biāo)，并允許其有一定的“浮動(dòng)”空間。采用拉丁超立方的試驗(yàn)設(shè)計(jì)方法，選取n組不同的目標(biāo)壓力系數(shù)分布，針對(duì)每一組目標(biāo)壓力系數(shù)分布進(jìn)行反問(wèn)題求解，并對(duì)反問(wèn)題求解出的葉型進(jìn)行氣動(dòng)性能分析，對(duì)比葉柵靜壓升、總壓損失系數(shù)以及出口氣流角，從而得出較優(yōu)的反問(wèn)題葉型。建立以葉型壓力分布為自變量的優(yōu)化問(wèn)題，稱(chēng)為“氣動(dòng)反問(wèn)題優(yōu)化設(shè)計(jì)”。

3 反問(wèn)題設(shè)計(jì)方法的應(yīng)用

某風(fēng)扇中截面葉型采用NACA65［17］，葉柵稠度為1.4，葉片安裝角16°，進(jìn)口氣流相對(duì)葉片弦線的攻角為30度，擴(kuò)壓因子為0.56，超出了最佳氣動(dòng)特性范圍［2］。原葉柵流動(dòng)的CFD分析結(jié)果表明，在葉片吸力面的85%弦長(zhǎng)位置，流動(dòng)出現(xiàn)了分離，使得葉型流動(dòng)損失明顯增加，如圖7所示。為了提高該風(fēng)扇的氣動(dòng)性能，采用本文建立的反問(wèn)題優(yōu)化設(shè)計(jì)方法對(duì)該葉型進(jìn)行設(shè)計(jì)，獲得新的氣動(dòng)葉型。

圖7 原始葉型速度流線圖及壓力云圖

首先采用拉丁超立方的試驗(yàn)設(shè)計(jì)方法，給定控制點(diǎn)X3，X4，X5的上下限（0.8，1.2）、（1.6，2.6）、（0.12，0.28），如圖6所示，生成20組目標(biāo)壓力系數(shù)分布。對(duì)每一組目標(biāo)壓力系數(shù)分布進(jìn)行反問(wèn)題設(shè)計(jì)，并計(jì)算出反設(shè)計(jì)結(jié)果對(duì)應(yīng)的壓升δp、總壓損失系數(shù)cp以及出口氣流角β等參數(shù)。然后針對(duì)這20組樣本進(jìn)行分析，選取較優(yōu)的壓力系數(shù)分布作為最終的反問(wèn)題設(shè)計(jì)結(jié)果。

圖8是優(yōu)化設(shè)計(jì)過(guò)程中，反問(wèn)題設(shè)計(jì)葉型的表面壓力系數(shù)與目標(biāo)值之間的殘差收斂曲線。最終殘差穩(wěn)定在初始?xì)埐畹?%左右，殘差下降了2個(gè)數(shù)量級(jí)。圖9是反問(wèn)題葉型控制點(diǎn)殘差的收斂曲線，最終趨于穩(wěn)定，與壓力系數(shù)殘差趨勢(shì)一致。由圖9可知，第3個(gè)控制點(diǎn)的變化較大，其余2個(gè)葉型控制點(diǎn)的變化較小。

圖8 反問(wèn)題設(shè)計(jì)壓力系數(shù)殘差變化曲線

圖9 反問(wèn)題葉型控制點(diǎn)殘差變化曲線

目標(biāo)壓力系數(shù)分布、計(jì)算壓力系數(shù)分布以及原始?jí)毫ο禂?shù)分布，如圖10所示，計(jì)算的壓力系數(shù)分布與目標(biāo)值較為吻合。原始葉型吸力面，在x/Cax=0.4的軸向位置，壓力分布有一個(gè)“臺(tái)階”，在x/Cax=0.85的軸向位置其后流動(dòng)壓力幾乎保持不變，也說(shuō)明x/Cax=0.85以后，吸力面流動(dòng)是分離的。目標(biāo)壓力系數(shù)分布表明，x/Cax=0到x/Cax=0.15的吸力面壓力系數(shù)提高，x/Cax=0.15到x/Cax=0.7之間壓力系數(shù)減小。對(duì)比反設(shè)計(jì)葉型與原始葉型發(fā)現(xiàn)，x/Cax=0.2到x/Cax=0.8部分的葉型厚度較原始葉型有所減薄，如圖11所示。

圖10 原始葉型與反設(shè)計(jì)葉型壓力系數(shù)對(duì)比

圖11 原始葉型與反設(shè)計(jì)葉型對(duì)比

4 反問(wèn)題設(shè)計(jì)的驗(yàn)證

本文采用CFX14.0軟件，對(duì)反問(wèn)題設(shè)計(jì)結(jié)果進(jìn)行數(shù)值驗(yàn)證。其中湍流模型選擇帶 γ-θ 轉(zhuǎn)捩的SST模型，計(jì)算域采用結(jié)構(gòu)化網(wǎng)格，葉片周?chē)捎肙型網(wǎng)格，壁面的Y+≈1，所有網(wǎng)格單元的最小內(nèi)角正交性都在36°以上，拉伸比＜1000。

圖12為葉型優(yōu)化后，葉柵內(nèi)流動(dòng)的速度流線及壓力云圖。葉柵內(nèi)分離區(qū)較原始葉型有所縮小，葉片尾跡寬度減小。

圖12 反設(shè)計(jì)葉型速度流線及壓力云圖

圖13 為反設(shè)計(jì)葉型與原始葉型吸力面的曲率對(duì)比，由圖可知在x/Cax=0.2到x/Cax=0.78之間，反設(shè)計(jì)葉型要比原始葉型的曲率小，而在x/Cax=0.178到尾緣處則相反。由于葉型曲率的改變，從而改善了整個(gè)流場(chǎng)的流動(dòng)，提高了葉型的性能。曲率分布對(duì)流動(dòng)形態(tài)的作用，還有待后續(xù)研究。

圖13 反設(shè)計(jì)與原始葉型的吸力面曲率對(duì)比

表1對(duì)比了反設(shè)計(jì)后葉型阻力系數(shù)、升力系數(shù)、升阻比以及總壓損失系數(shù)等參數(shù)。新的葉型阻力系數(shù)減少16.4%了，升力系數(shù)幾乎不變，總壓損失系數(shù)減小了15.6%，總的性能有了較大的提升，證明了本文反問(wèn)題設(shè)計(jì)方法是可行的。

表1 原始葉型與反設(shè)計(jì)葉型性能參數(shù)對(duì)比

5 結(jié)論

（1） 本文采用的Head邊界層厚度計(jì)算方法，結(jié)合優(yōu)化算法求解最優(yōu)的壓力系數(shù)分布的方法具有可行性。成功地將邊界層厚度的計(jì)算應(yīng)用到葉型的反設(shè)計(jì)中，并結(jié)合優(yōu)化算法求解出近似最優(yōu)的壓力系數(shù)分布。

（2） 葉型的修改是反問(wèn)題的核心，本文采用的“彈性模型”修改葉型的方法，成功的設(shè)計(jì)出了滿足要求的葉型，表明本文的葉型修改方法是成功的。同時(shí)表明本文建立的反問(wèn)題求解程序不僅具有求解速度快，而且具有較好的實(shí)用性以及可靠性。

［1］ 王仲奇，秦仁，透平機(jī)械原理［M］. 北京：機(jī)械工業(yè)出版社，1988.

［2］ Herrig L J，Emery J C，Erwin J R. Systematic twodimensional cascade tests of NACA 65-Series compressor blades at low speeds［R］. NACA TN3916 Report Archive & Image Library，1951.

［3］ 軸流壓氣機(jī)氣動(dòng)設(shè)計(jì)［M］. 北京：國(guó)防工業(yè)出版社，1975.

［4］ Hobbs D E，Weingold H D. Development of controlled diffusion airfoils for multistage compressor application［J］. Journal of Engineering for Gas Turbines and Power，1984，106（2）：271-278.

［5］ Somers D M. Design and experimental results for a natural-laminar-flow airfoil for general aviation applications［R］. NASA Technocal Paper 1865，1981.

［6］ Gopalarathnam A，Selig M S. Low-speed naturallaminar-flow airfoils：Case study in inverse airfoil design［J］. Journal of aircraft，2001，38（1）：57-63.

［7］ Shahpar S. A Comparative Study of Optimisation Methods for Aerodynamic Design of Turbomachinery Blades［C］. ASME 2000-GT-523 - ASME Turbo Expo 2000：Power for Land，Sea，and Air. 2000：V001T03A087-V001T03A087.

［8］ Kasra Daneshkhah，Wahid S. Ghaly. An inverse blade design method for subsonic and transonic viscous flow in compressors and turbines［J］. Inverse Problems in Science & Engineering，2006，14（3）：211-231.

［9］ Dulikravich G S，Dennis B H. Inverse design and optimization using CFD［C］.European Congress on Computational Methods in Applied Sciences and Engineering，ECCOMAS，2000.

［10］ 姚征. 基于渦量修正的葉柵氣動(dòng)優(yōu)化設(shè)計(jì)［J］. 空氣動(dòng)力學(xué)學(xué)報(bào)，2001，19（4）：376-382.

［11］ 戴韌，劉高聯(lián). 旋成面葉柵氣動(dòng)反命題的新解法（Ⅱ）［J］. 航空動(dòng)力學(xué)報(bào)，1994（1）：4.

［12］ Pehlivanoglu Y V，Hacioglu A. Inverse design of 2-D airfoil via vibrational genetic algorithm［J］. Journal of Aeronautics and Space Technologies，2006，2（4）：7-14.

［13］ Obayashi S，Takanashi S. Genetic optimization of target pressure distributions for inverse design methods［J］.AIAA Journal，1996，34（5）：881-886.

［14］ 豐鎮(zhèn)平，厲海濤，宋立明，等. 基于控制理論的透平葉柵氣動(dòng)反設(shè)計(jì)優(yōu)化［J］. 中國(guó)科學(xué)：技術(shù)科學(xué)，2013（3）：257-273.

［15］ Cebeci T，Bradshaw P. Momentum transfer in boundary layers［M］. Washington，DC，Hemisphere Publishing Corp.；New York，McGraw-Hill Book Co.，1977.

［16］ 同濟(jì)大學(xué)計(jì)算數(shù)學(xué)教研室. 數(shù)值分析基礎(chǔ)［M］.上海：同濟(jì)大學(xué)出版社，1998.

［17］ 陳榴，劉雪驕，戴韌.彎曲動(dòng)葉對(duì)高負(fù)荷風(fēng)扇氣動(dòng)性能的影響［C］. 中國(guó)工程熱物理學(xué)會(huì)-流體機(jī)械.2013.