基于離散伴隨方法的透平葉柵壓力反設(shè)計(jì)研究

張朝磊，吳海燕

（1.江西理工大學(xué)機(jī)電工程學(xué)院，江西贛州341000；2.西安交通大學(xué)葉輪機(jī)械研究所，西安710049）

基于離散伴隨方法的透平葉柵壓力反設(shè)計(jì)研究

張朝磊1,2，吳海燕1

（1.江西理工大學(xué)機(jī)電工程學(xué)院，江西贛州341000；2.西安交通大學(xué)葉輪機(jī)械研究所，西安710049）

根據(jù)離散伴隨方法理論和自動(dòng)微分技術(shù)，由流場(chǎng)求解器源代碼構(gòu)造了相應(yīng)的離散伴隨場(chǎng)求解器.通過(guò)耦合參數(shù)化程序、網(wǎng)格生成程序、流場(chǎng)求解器和伴隨場(chǎng)求解器，建立了適用于葉輪機(jī)械葉柵的氣動(dòng)優(yōu)化設(shè)計(jì)系統(tǒng).利用該系統(tǒng)，對(duì)某二維跨音速透平葉柵在給定葉型表面目標(biāo)壓力分布的情況下，通過(guò)構(gòu)造合適的目標(biāo)函數(shù)將葉柵反設(shè)計(jì)問(wèn)題轉(zhuǎn)化為葉柵氣動(dòng)優(yōu)化設(shè)計(jì)問(wèn)題，成功進(jìn)行了氣動(dòng)反設(shè)計(jì).反設(shè)計(jì)結(jié)果表明，文中建立的葉柵氣動(dòng)優(yōu)化設(shè)計(jì)系統(tǒng)能夠根據(jù)給定的葉型表面壓力分布有效地進(jìn)行壓力反設(shè)計(jì)，驗(yàn)證了該方法的正確性與有效性.

離散伴隨方法；自動(dòng)微分；透平葉柵；反設(shè)計(jì)；靈敏度

0 前言

基于控制理論的伴隨方法能夠?qū)崿F(xiàn)快速精確的敏感度分析，并且敏感度計(jì)算量與設(shè)計(jì)變量數(shù)無(wú)關(guān)，有效克服了傳統(tǒng)敏感度分析方法計(jì)算量正比于設(shè)計(jì)變量數(shù)的問(wèn)題[1].基于伴隨系統(tǒng)建立方式的不同，伴隨方法可以分為連續(xù)伴隨方法和離散伴隨方法.與連續(xù)伴隨方法相比，離散伴隨方法具有目標(biāo)函數(shù)適用性廣，伴隨系統(tǒng)建立過(guò)程清晰、規(guī)范，存在成熟的自動(dòng)微分工具用于輔助開(kāi)發(fā)與考核伴隨系統(tǒng)等優(yōu)勢(shì).自上世紀(jì)九十年代以來(lái)該方法得到了相關(guān)研究者的關(guān)注[2-5].目前，離散伴隨方法在外流領(lǐng)域取得了一定的研究成果，在機(jī)翼乃至整機(jī)的氣動(dòng)優(yōu)化設(shè)計(jì)中得到了廣泛應(yīng)用.在內(nèi)流領(lǐng)域，離散伴隨方法自本世紀(jì)初以來(lái)逐步得到關(guān)注與發(fā)展，目前，研究人員已經(jīng)成功將其應(yīng)用于葉輪機(jī)械葉柵的氣動(dòng)優(yōu)化設(shè)計(jì)與反設(shè)計(jì)[6-11].

基于文獻(xiàn)[10-11]中的離散伴隨求解器實(shí)現(xiàn)方法，編程實(shí)現(xiàn)了無(wú)粘環(huán)境下的離散伴隨場(chǎng)求解器；建立了適用于透平葉柵氣動(dòng)反設(shè)計(jì)的設(shè)計(jì)系統(tǒng)，并利用該系統(tǒng)對(duì)某二維跨音速透平葉柵在給定葉型表面目標(biāo)壓力分布的情況下，通過(guò)構(gòu)造適當(dāng)?shù)哪繕?biāo)函數(shù)對(duì)該葉柵進(jìn)行了氣動(dòng)反設(shè)計(jì).

1 離散伴隨系統(tǒng)原理

通過(guò)構(gòu)造適當(dāng)?shù)哪繕?biāo)函數(shù)，可以將葉輪機(jī)械葉柵的氣動(dòng)反設(shè)計(jì)問(wèn)題轉(zhuǎn)化為優(yōu)化設(shè)計(jì)問(wèn)題.通常，葉輪機(jī)械葉柵優(yōu)化設(shè)計(jì)對(duì)應(yīng)的最小值問(wèn)題可表示為：

其中，目標(biāo)函數(shù)J為設(shè)計(jì)變量α和狀態(tài)變量u的函數(shù)，同時(shí)，狀態(tài)變量u本身受狀態(tài)方程（Euler方程）R（α，u）=0約束.引入伴隨變量υ，定義新的目標(biāo)函數(shù)I：

線(xiàn)性化上述目標(biāo)函數(shù)，令目標(biāo)函數(shù)對(duì)狀態(tài)變量的變分系數(shù)為零，消除目標(biāo)函數(shù)對(duì)狀態(tài)變量敏感度的依賴(lài)，并引入時(shí)間相關(guān)項(xiàng)，可得到離散伴隨系統(tǒng)方程如下：

伴隨方程（3）中不存在對(duì)設(shè)計(jì)變量α的導(dǎo)數(shù)，因此伴隨方程的求解獨(dú)立于設(shè)計(jì)變量α的個(gè)數(shù).在求得伴隨變量υ后，目標(biāo)函數(shù)的梯度可按下式計(jì)算：

由公式（3）和（4）可知，為了求解目標(biāo)函數(shù)對(duì)設(shè)計(jì)變量的梯度，關(guān)鍵在于計(jì)算J/u、J/α、υT（R/α）以及υT（R/u），其中前三項(xiàng)的計(jì)算不需要迭代計(jì)算，而最后一項(xiàng)的求解必須迭代求解.

2 自動(dòng)微分

自動(dòng)微分是一種由某計(jì)算程序自動(dòng)生成對(duì)應(yīng)導(dǎo)數(shù)計(jì)算程序的技術(shù)，其理論依據(jù)為微積分的鏈導(dǎo)法則和程序各步操作可微性的假設(shè).自動(dòng)微分工具的實(shí)現(xiàn)方式主要包括兩種，即基于源代碼轉(zhuǎn)換實(shí)現(xiàn)和基于運(yùn)算操作符重載實(shí)現(xiàn).通常，基于運(yùn)算符重載的自動(dòng)微分工具需要計(jì)算機(jī)語(yǔ)言支持運(yùn)算符重載，自動(dòng)微分過(guò)程不重新生成代碼，一般情況下其導(dǎo)數(shù)程序計(jì)算效率較低；基于源代碼轉(zhuǎn)換的自動(dòng)微分工具通常必須根據(jù)原有計(jì)算程序另外生成導(dǎo)數(shù)計(jì)算程序，并且生成的自動(dòng)微分程序計(jì)算效率也比較高.

自動(dòng)微分可分為兩種模式，即前向模式和后向模式.通過(guò)自動(dòng)微分，對(duì)于一個(gè)標(biāo)量函數(shù)f和某一向量，前向模式給出了其方向?qū)?shù)Δf·，而其后向模式給出了其梯度Δf；而對(duì)于一個(gè)向量函數(shù)R和某一向量，前向模式給出了ΔR·，而后向模式給出了（ΔR）T·.因此，采用后向模式對(duì)計(jì)算流體力學(xué)源代碼中殘差子程序進(jìn)行自動(dòng)微分處理，即可獲得公式（3）和（4）中所需要的各項(xiàng).本研究中采用的自動(dòng)微分工具為T(mén)apenade[12]，該工具基于源代碼轉(zhuǎn)換實(shí)現(xiàn).

3 離散伴隨系統(tǒng)求解

3.1 靈敏度傳遞

在葉輪機(jī)械的數(shù)值計(jì)算中，從對(duì)應(yīng)于某幾何模型的一組設(shè)計(jì)變量α到某一確定的計(jì)算網(wǎng)格X，再到流場(chǎng)數(shù)值解U和目標(biāo)函數(shù)J，其間關(guān)系可表述如下：

3.2 伴隨場(chǎng)迭代求解

伴隨場(chǎng)的迭代求解策略與文獻(xiàn)[10,11]中相似，本研究中伴隨變量的迭代求解公式可寫(xiě)為：

為了獲得目標(biāo)函數(shù)對(duì)設(shè)計(jì)變量的導(dǎo)數(shù)，在求得伴隨場(chǎng)之后，利用公式（8）間接實(shí)現(xiàn)目標(biāo)函數(shù)對(duì)設(shè)計(jì)變量梯度的計(jì)算.即：

4 優(yōu)化流程

與文獻(xiàn)[10]中建立的葉輪機(jī)械葉柵氣動(dòng)優(yōu)化設(shè)計(jì)系統(tǒng)相似，文中建立的葉輪機(jī)械葉柵反設(shè)計(jì)系統(tǒng)包括參數(shù)化、網(wǎng)格生成、流場(chǎng)計(jì)算、伴隨場(chǎng)計(jì)算、梯度計(jì)算及優(yōu)化子程序六部分.其中葉柵參數(shù)化基于非均勻B-Spline實(shí)現(xiàn)；網(wǎng)格生成基于邊界規(guī)范化方法和Laplace網(wǎng)格光順技術(shù)實(shí)現(xiàn)；梯度計(jì)算過(guò)程中計(jì)算網(wǎng)格對(duì)設(shè)計(jì)變量的梯度基于中心差分方法實(shí)現(xiàn)，優(yōu)化子程序基于擬牛頓法實(shí)現(xiàn).整個(gè)葉柵反設(shè)計(jì)系統(tǒng)流程如圖1所示.

圖1 反設(shè)計(jì)系統(tǒng)組成及流程

5 反設(shè)計(jì)算例

基于建立的葉柵反設(shè)計(jì)系統(tǒng)，針對(duì)某跨音速葉柵在粘性環(huán)境下進(jìn)行了氣動(dòng)反設(shè)計(jì).該葉柵的氣動(dòng)參數(shù)如表1所示.

表1 葉柵氣動(dòng)參數(shù)

文中在對(duì)葉柵吸力面和壓力面幾何型線(xiàn)進(jìn)行參數(shù)化之前，首先根據(jù)設(shè)定的前緣和尾緣比例因子進(jìn)行葉柵前緣、尾緣和中間區(qū)域的分割，然后對(duì)吸力面和壓力面的中部區(qū)域進(jìn)行參數(shù)化.鑒于葉柵的優(yōu)化或反設(shè)計(jì)過(guò)程中葉柵前緣和尾緣的幾何型線(xiàn)通常是保持不變的.因此，采用文中的參數(shù)化方法，即僅針對(duì)中部區(qū)域進(jìn)行參數(shù)化不僅可以保證葉柵的前緣和尾緣形狀，減少參數(shù)化所需要的控制點(diǎn)數(shù)目，而且可以有效消除葉柵尾緣附近區(qū)域葉柵幾何型線(xiàn)的波動(dòng).文中的葉柵參數(shù)化基于上述方法，首先設(shè)定葉柵吸力面前緣和尾緣占吸力面型線(xiàn)長(zhǎng)度的5%為固定不變區(qū)域，其余部分為中部區(qū)域.分割后中部區(qū)域葉柵型線(xiàn)采用10個(gè)控制點(diǎn)進(jìn)行參數(shù)化，參數(shù)化的結(jié)果如圖2所示.在反設(shè)計(jì)過(guò)程中，控制點(diǎn)3rd-8th的坐標(biāo)選為設(shè)計(jì)變量，變量總數(shù)為12.

葉柵網(wǎng)格生成采用H-O-H型網(wǎng)格拓?fù)浣Y(jié)構(gòu)，各塊網(wǎng)格生成采用邊界規(guī)范化方法，然后采用Laplace方法進(jìn)行光順，計(jì)算網(wǎng)格節(jié)點(diǎn)數(shù)均為9467.計(jì)算網(wǎng)格示意圖如圖3所示.

流場(chǎng)和伴隨場(chǎng)計(jì)算過(guò)程中邊界條件給定如下：進(jìn)口給定總壓、總溫和進(jìn)口氣流方向；出口給定平均靜壓；葉片表面給定絕熱無(wú)滑移壁面邊界條件.計(jì)算區(qū)域?yàn)閱瓮ǖ?，葉柵兩側(cè)面網(wǎng)格給定位移型周期性邊界條件.

圖2 葉柵參數(shù)化結(jié)果

圖3 計(jì)算網(wǎng)格

反設(shè)計(jì)過(guò)程中，通過(guò)給定葉柵表面目標(biāo)壓力分布，構(gòu)造目標(biāo)函數(shù)如下：

其中pd為目標(biāo)壓力分布.優(yōu)化設(shè)計(jì)采用擬牛頓算法，反設(shè)計(jì)過(guò)程中目標(biāo)函數(shù)收斂史如圖4所示，圖中目標(biāo)函數(shù)值以初始值為參考進(jìn)行了歸一化處理，經(jīng)過(guò)11次迭代后收斂，期間共計(jì)算了12次伴隨場(chǎng)，24次流場(chǎng).優(yōu)化前后目標(biāo)函數(shù)對(duì)設(shè)計(jì)變量的梯度對(duì)比如圖5所示，從圖可見(jiàn)，反設(shè)計(jì)結(jié)束時(shí)目標(biāo)函數(shù)對(duì)設(shè)計(jì)變量的梯度趨近于零.

圖4 目標(biāo)函數(shù)收斂史

圖5 設(shè)計(jì)前、后設(shè)計(jì)變量靈敏度對(duì)比

反設(shè)計(jì)前、后葉柵表面壓力對(duì)比如圖6所示.反設(shè)計(jì)前后葉柵型線(xiàn)對(duì)比如圖7所示.從圖可見(jiàn)，反設(shè)計(jì)后，葉柵表面壓力與給定的目標(biāo)壓力分布基本重合，驗(yàn)證了文中建立的葉輪機(jī)械葉柵反設(shè)計(jì)系統(tǒng)的正確性與有效性.

圖6 反設(shè)計(jì)前、后葉片壁面壓力對(duì)比

圖7 反設(shè)計(jì)前、后葉型對(duì)比

6 結(jié)論

文中建立了基于離散伴隨方法的葉輪機(jī)械葉柵氣動(dòng)優(yōu)化設(shè)計(jì)系統(tǒng)，并利用該系統(tǒng)在無(wú)粘環(huán)境下根據(jù)指定的目標(biāo)壓力分布對(duì)一跨音速葉柵進(jìn)行了壓力反設(shè)計(jì)，反設(shè)計(jì)結(jié)果驗(yàn)證了文中建立的葉柵氣動(dòng)優(yōu)化設(shè)計(jì)系統(tǒng)的正確性與有效性.文中提出的葉柵參數(shù)化方法和開(kāi)發(fā)的離散伴隨場(chǎng)求解器可直接應(yīng)用于三維葉輪機(jī)械葉柵的氣動(dòng)反設(shè)計(jì)，為下一步葉柵三維反設(shè)計(jì)優(yōu)化奠定了基礎(chǔ).

[1]劉靜，丁凌蓉.基于Matlab的低矮式破碎機(jī)V帶傳動(dòng)模糊優(yōu)化設(shè)計(jì)[J].江西理工大學(xué)學(xué)報(bào)，2009，30(2)：24-27.

[2]Elliott J，Peraire J.Practical 3D aerodynamic design and optimization using unstructured meshes[J].AIAA Journal，1997，35(9)：1479-1485.

[3]Mohammadi B.Optimal shape design，reverse mode of automatic differentiation and turbulence[C].AIAA Paper，1997，97-99.

[4]Giles M B，Pierce N A.Adjoint equations in CFD：duality，boundary conditions and solution behavior[C].AIAA Paper，1997，97-1850.

[5]Marta A C，Mader C A，Martins J R R A，et al.A methodology for the developmentofdiscreteadjointsolversusingautomaticdifferentiation tools[J].InternationalJournalofComputationalFluidDynamics，2007，21(9-10)：307-327.

[6]FloreaR，HallKC.Sensitivityanalysisofunsteadyinviscidflowthrough turbomachinery cascades[J].AIAA Journal，2001，39(6)：1047-1056.

[7]Campobasso M S，Duta M C，Giles M B.Adjoint methods for turbomachinery design[C].ISABE Conference，2001.

[8]DutaMC，GilesMB，CampobassoM S.The harmonic adjoint approach to unsteady turbomachinery design[J].International Journal for Numerical MethodsinFluids，2002，3(3-4)：323-332.

[9]Campobasso M S，Duta M C，Giles M B.Adjoint calculation of sensitivities of turbomachinery objective functions[J].AIAA Journal of Propulsion and Power，2003，19(4)：693-703.

[10]張朝磊，厲海濤，豐鎮(zhèn)平.基于離散伴隨方法的透平葉柵氣動(dòng)優(yōu)化設(shè)計(jì)[J].工程熱物理學(xué)報(bào)，2012，33(1)：47-50.

[11]Zhang C L，F(xiàn)eng Z P.Aerodynamic shape design optimization for turbomachinery cascade based on discrete adjoint method[C].ASME Paper GT2011-45805.

[12]Hascoёt L，Pascual V.Tapenade 2.1 user’s guide.Technical Report 0300[R]，INRIA，2004.

Study on pressure inverse design of turbomachinery cascade based on discrete adjoint method

ZHANG Chao-lei1,2,WU Hai-yan1
(1.School of Mechanical and Electrical Engineering,Jiangxi University of Science and Technology,Ganzhou 341000,China;2.Institute of Turbomachinery,Xi'an Jiaotong University,Xi'an 710049,China)

A discrete adjoint solver is constructed according to the discrete adjoint theory and the automatic differentiation technology based on an in-house flow solver code.An aerodynamic optimization design system for turbomachinery cascades is established by coupling a parameterization program,a structured multi-block grid generator,a flow solver and the discrete adjoint solver established in this work.Based on the optimization system,the pressure inverse design problem could be translated to an optimization problem by constructing the objective function appropriately.An inverse design for a typical 2D transonic turbomachinery cascade is preformed under the inviscid flow environment according to the target pressure distribution specified beforehand.The inverse design result indicates that the optimization system established in this article could implement the pressure inverse design efficiently.The validity and efficiency of the optimization design system are proved.

discrete adjoint method;automatic differentiation;turbomachinery cascade;inverse design;sensitivity

V228.7

A

2012-02-13

國(guó)家自然科學(xué)基金資助項(xiàng)目（51076121）；江西理工大學(xué)科研基金項(xiàng)目計(jì)劃（JXXJ11038）

張朝磊（1979-），男，在讀博士，講師，主要從事葉輪機(jī)械氣動(dòng)優(yōu)化設(shè)計(jì)等方面的研究，E-mail：zhangchaolei@yahoo.cn.

2095-3046（2012）03-0047-04