王 威,王 軍,梁 鐘,李佳峻

(華中科技大學(xué) 能源與動(dòng)力工程學(xué)院,武漢 430074)

0 引 言

翼型廣泛應(yīng)用于航空航天、風(fēng)力發(fā)電等諸多流體機(jī)械領(lǐng)域,其性能對(duì)流體機(jī)械的工作性能起著至關(guān)重要的影響。高性能葉片和飛行器的設(shè)計(jì)中已不直接采用現(xiàn)有的翼型,而是對(duì)翼型進(jìn)行優(yōu)化設(shè)計(jì),尋找到功能符合設(shè)計(jì)要求的高性能翼型。傳統(tǒng)的翼型優(yōu)化[1-4]在研究其升阻力系數(shù)、效率等空氣動(dòng)力學(xué)性能、結(jié)構(gòu)特性等單一學(xué)科要求的基礎(chǔ)上,來(lái)進(jìn)行翼型的改進(jìn)、修型,但這些參數(shù)不能反映流場(chǎng)內(nèi)局部損失的具體來(lái)源和增加過(guò)程。如果能夠掌握流場(chǎng)中能量損失的大小及分布,可以針對(duì)局部損失較大的區(qū)域進(jìn)行優(yōu)化設(shè)計(jì),通過(guò)降低局部損失來(lái)改善壓氣機(jī)性能。

熵產(chǎn)理論基于熱力學(xué)第二定律,在不可逆現(xiàn)象比較集中的傳熱與流動(dòng)系統(tǒng)得到了廣泛深入的研究。由于熵產(chǎn)的大小揭示了過(guò)程的不可逆程度,反映了過(guò)程中能量轉(zhuǎn)換與利用的完善程度,國(guó)內(nèi)外學(xué)者創(chuàng)造性地將其用于翼型、風(fēng)機(jī)葉片以及飛行器的設(shè)計(jì)中。Shehata等[5]將熵產(chǎn)最小原理用于葉片性能分析中,分析了不同迎角下、不同翼型截面的熱力學(xué)第二效率。Mortazavi等[6]將熵分析方法和神經(jīng)網(wǎng)絡(luò)算法引入軸流風(fēng)機(jī)葉片翼型截面優(yōu)化中,分析了低雷諾數(shù)下NREL系列翼型截面的火用效率。Li和Figliola[7]將熵產(chǎn)作為優(yōu)化目標(biāo)對(duì)二維翼型進(jìn)行多目標(biāo)氣動(dòng)優(yōu)化設(shè)計(jì),并分析了熵產(chǎn)與翼型阻力的關(guān)系。王松齡等[8]基于熵產(chǎn)理論對(duì)離心風(fēng)機(jī)葉輪進(jìn)行參數(shù)化優(yōu)化設(shè)計(jì),優(yōu)化后葉輪和蝸殼熵產(chǎn)明顯降低,流動(dòng)得到改善。

本文將熵產(chǎn)方法引入到翼型優(yōu)化設(shè)計(jì)中,研究流場(chǎng)熵產(chǎn)在翼型氣動(dòng)優(yōu)化設(shè)計(jì)中的作用。通過(guò)建立參數(shù)化、網(wǎng)格生成、CFD計(jì)算、數(shù)值優(yōu)化方法一體化的翼型優(yōu)化平臺(tái),完成特定工況下以最大升阻比和最小熵產(chǎn)率為優(yōu)化目標(biāo)的多目標(biāo)翼型優(yōu)化設(shè)計(jì)。

1 熵分析優(yōu)化模型建立及求解

1.1 系統(tǒng)熵方程及熵產(chǎn)計(jì)算

熵產(chǎn)代表了流動(dòng)中能量損失的大小和系統(tǒng)的不可逆程度,根據(jù)熱力學(xué)第二定律,單位體積流體的熵產(chǎn)率[9-10]為:

其中Φ為耗散函數(shù),其張量形式表達(dá)式為:

式(1)表明引起流場(chǎng)熵產(chǎn)的因素是黏性引起的耗散和有限溫差引起的傳熱,其中耗散引起的熵產(chǎn)取決于流體的黏性和速度場(chǎng),傳熱引起的熵產(chǎn)取決于流體的導(dǎo)熱率和溫度場(chǎng)。

當(dāng)?shù)退俨豢蓧呵也豢紤]傳熱時(shí),由傳熱項(xiàng)引起的流場(chǎng)熵產(chǎn)可忽略不計(jì)。因此,對(duì)于二維不可壓流場(chǎng)中單位體積熵產(chǎn)率可進(jìn)一步簡(jiǎn)化為:

當(dāng)流動(dòng)為湍流時(shí),本文采取雷諾平均法(RANS)對(duì)上式進(jìn)行時(shí)均運(yùn)算,由于T′作為分母出現(xiàn)在高階項(xiàng),故可忽略,可得:

式(4)表明,湍動(dòng)流動(dòng)過(guò)程中的耗散所引起的熵產(chǎn)可以分成兩項(xiàng),第一項(xiàng)是由于黏性耗散引起的,第二項(xiàng)是由于速度脈動(dòng)所引起的湍流耗散產(chǎn)生。

根據(jù)Boussinesq的渦黏性假設(shè),有

式中μt為湍流動(dòng)力黏性系數(shù),它是空間坐標(biāo)函數(shù),取決于流動(dòng)狀態(tài)而不是物性參數(shù)。

對(duì)式(6)進(jìn)行體積分,即可得到流場(chǎng)熵產(chǎn)率:

1.2 翼型參數(shù)化方法

本文采用Kulfan和Bussoletti[11-12]提出的CST(Class Shape Transformation)參數(shù)化方法,對(duì)于鈍前緣、尖尾緣翼型,翼型幾何型面的具體參數(shù)表達(dá)式如下:

式中,c為翼型弦長(zhǎng),x、y為翼型橫、縱坐標(biāo),Δzte/c為翼型后緣點(diǎn)縱坐標(biāo),Rle為翼型前緣半徑,β為后緣處的切線傾角,bi(i＝1,2,…,n－1)為引進(jìn)的權(quán)重因子。

為限制翼型上下表面曲線在后緣處封閉,故上下表面曲線的Δzte/c相等。因此,翼型的設(shè)計(jì)參數(shù)為翼型前緣半徑Rle,后緣點(diǎn)縱坐標(biāo)Δzte/c,上表面曲線在后緣處的切線傾角β和翼型后緣角βte,上、下表面控制參數(shù),(i＝1,2,…,n－1)。

本文分別采用7階Bernstein多項(xiàng)式對(duì)NACA0012和RAE2822進(jìn)行參數(shù)化擬合,并與翼型數(shù)據(jù)庫(kù)中翼型數(shù)據(jù)進(jìn)行對(duì)比,擬合效果如圖1所示,由圖可知,擬合出的翼型如與初始翼型符合得較好,充分體現(xiàn)了CST參數(shù)化方法的優(yōu)點(diǎn)。

圖1 NACA0012和RAE2822翼型CST參數(shù)化效果圖Fig.1 CST parametric representation of NACA0012 and RAE2822

1.3 流場(chǎng)數(shù)值計(jì)算

翼型的計(jì)算域和計(jì)算網(wǎng)格如圖2所示,翼型弦長(zhǎng)為c;入口邊界距機(jī)翼前緣11.5c,為半徑12.5c的半圓;出口邊界距機(jī)翼后緣,采用C型結(jié)構(gòu)網(wǎng)格,用四套網(wǎng)格:200×70、200×80、200×90和200×100進(jìn)行網(wǎng)格無(wú)關(guān)性驗(yàn)證,本文采用200×80網(wǎng)格進(jìn)行計(jì)算。

圖2 NACA0012翼型流場(chǎng)計(jì)算網(wǎng)格圖Fig.2 Computational mesh of NACA0012 flow field

對(duì)翼型進(jìn)行優(yōu)化前,需要對(duì)數(shù)值計(jì)算結(jié)果進(jìn)行驗(yàn)證,將其與實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對(duì)比,從而對(duì)計(jì)算的可信度和優(yōu)化質(zhì)量提供客觀的評(píng)價(jià)標(biāo)準(zhǔn)。本文計(jì)算工況為來(lái)流雷諾數(shù)2.88×106,進(jìn)出口溫度為299.15 K,計(jì)算迎角分別為0°、5°、10°和15°。湍流模型采用S-A單方程模型,邊界條件設(shè)置為速度入口、壓力出口,采用二階迎風(fēng)格式離散控制方程,選擇壓力基求解器,速度和壓力耦合采用SIMPLE算法,殘差控制在1×10－4數(shù)量級(jí)。

圖3為NACA0012翼型在迎角0°、10°和15°下的壓力系數(shù)分布與文獻(xiàn)[13]中實(shí)驗(yàn)數(shù)據(jù)的對(duì)比。其中實(shí)驗(yàn)來(lái)流速度為55 m/s,模型弦長(zhǎng)為0.76 m,雷諾數(shù)為2.88×106。表1為迎角0°、5°和10°下數(shù)值計(jì)算NACA0012翼型升、阻力系數(shù)和文獻(xiàn)[14]實(shí)驗(yàn)數(shù)據(jù)對(duì)比,可以看出:數(shù)值計(jì)算得到的壓力分布與風(fēng)洞實(shí)驗(yàn)數(shù)據(jù)非常接近,升力系數(shù)的數(shù)值計(jì)算結(jié)果與實(shí)驗(yàn)值相當(dāng),阻力計(jì)算存在一定誤差,數(shù)值結(jié)果表明此時(shí)流場(chǎng)網(wǎng)格劃分和數(shù)值計(jì)算的設(shè)置適用于該工況下翼型流場(chǎng)計(jì)算。

圖3 0°、10°和15°迎角下壓力系數(shù)分布與實(shí)驗(yàn)值比較Fig.3 Comparison of the pressure coefficient distribution obtained from calculation and experiment

表1 不同迎角下升、阻力系數(shù)與實(shí)驗(yàn)值對(duì)比Table 1 Comparison of the lift and drag coefficients obtained from calculation and experiment

1.4 優(yōu)化模型及流程

優(yōu)化的數(shù)學(xué)表達(dá)式為:

其中,D為優(yōu)化設(shè)計(jì)空間,d為優(yōu)化設(shè)計(jì)變量的向量,gi(X)為約束條件。

將優(yōu)化算法、CST參數(shù)化建模、翼型網(wǎng)格生成、流場(chǎng)CFD計(jì)算進(jìn)行耦合,生成適用于多目標(biāo)翼型優(yōu)化的計(jì)算程序。根據(jù)CST翼型參數(shù)化方法,輸入設(shè)計(jì)變量的值,生成翼型坐標(biāo)點(diǎn)數(shù)據(jù)文件,繼而完成翼型網(wǎng)格生成。通過(guò)翼型流場(chǎng)CFD計(jì)算,輸出翼型氣動(dòng)性能數(shù)據(jù)文件。

在翼型的優(yōu)化設(shè)計(jì)中,設(shè)計(jì)變量與目標(biāo)函數(shù)之間存在復(fù)雜的非線性關(guān)系,本文采用遺傳算法進(jìn)行優(yōu)化計(jì)算。對(duì)于多目標(biāo)優(yōu)化問(wèn)題,采用帶精英策略的非支配排序遺傳算法(NSGA2)[15]。它應(yīng)用Pareto最優(yōu)解理論來(lái)處理多目標(biāo)優(yōu)化問(wèn)題,由于各個(gè)目標(biāo)沒(méi)有人為地設(shè)定權(quán)重分配,所以所有Pareto最優(yōu)解的集合構(gòu)成了該問(wèn)題的Pareto最優(yōu)解集。此外NSGA2算法采用快速非支配排序方法和精英策略,降低了算法的計(jì)算復(fù)雜度,提高了計(jì)算效率。NSGA2優(yōu)化算法將翼型氣動(dòng)性能數(shù)據(jù)作為輸入文件,根據(jù)遺傳算法既定的輸入?yún)?shù),計(jì)算各個(gè)翼型的適應(yīng)度來(lái)決定下一代種群的個(gè)體(翼型),從而形成一個(gè)優(yōu)化循環(huán)。

2 優(yōu)化結(jié)果與分析

根據(jù)上述優(yōu)化流程,本文對(duì)NACA0012翼型進(jìn)行單目標(biāo)和多目標(biāo)優(yōu)化設(shè)計(jì)。算例一為基于熵產(chǎn)的多目標(biāo)優(yōu)化設(shè)計(jì)方法,優(yōu)化目標(biāo)為升阻比最大以及熵產(chǎn)最小,約束條件為升力系數(shù)以及設(shè)計(jì)變量的幾何約束。算例二為傳統(tǒng)的翼型氣動(dòng)優(yōu)化設(shè)計(jì)方法,優(yōu)化目標(biāo)為阻力系數(shù)最小,約束條件與算例一一致。其數(shù)學(xué)描述分別為:

翼型的設(shè)計(jì)狀態(tài)為來(lái)流雷諾數(shù)為2.88×106,飛行迎角5°,流場(chǎng)網(wǎng)格和CFD參數(shù)設(shè)置選擇上述經(jīng)驗(yàn)證的設(shè)置。翼型的參數(shù)化方法選擇CST方法,設(shè)計(jì)參數(shù)和相應(yīng)的約束范圍如表2所示。

表2 設(shè)計(jì)變量名稱(chēng)、取值范圍Table2 Range of design variables

優(yōu)化設(shè)計(jì)的各優(yōu)化控制參數(shù)選取為:初始種群規(guī)模為80,進(jìn)化代數(shù)為50,雜交概率為0.70,變異概率為0.12,變量采用實(shí)數(shù)編碼。

兩種算例優(yōu)化過(guò)程中的種群分布如圖4、圖5所示?？芍?初始種群分布混亂,會(huì)出現(xiàn)升阻比小的個(gè)體,經(jīng)過(guò)數(shù)代進(jìn)化后,通過(guò)雜交、變異,淘汰了部分劣勢(shì)個(gè)體,種群的分布朝好的方向進(jìn)化,但仍呈現(xiàn)不規(guī)律性。第10代以后有收斂的趨勢(shì),在以后的優(yōu)化中,收斂趨勢(shì)更加明顯,不斷優(yōu)化出升阻比更大、熵產(chǎn)率或阻力系數(shù)更小的個(gè)體,說(shuō)明優(yōu)化的全局收斂性。

圖4 算例一優(yōu)化過(guò)程種群分布情況Fig.4 Population distribution in optimization procedure in case 1

圖6(a、b)為兩種算例下優(yōu)化后翼型目標(biāo)函數(shù)的Pareto分布,可以看出采用基于熵產(chǎn)率作為優(yōu)化目標(biāo)的算例一得到的優(yōu)化翼型具有相對(duì)較大的升阻比。圖中Pareto分布中每個(gè)數(shù)據(jù)點(diǎn)代表一種可行的優(yōu)化翼型,它是兩種互為矛盾的優(yōu)化目標(biāo)不斷調(diào)和的結(jié)果,因此可以根據(jù)設(shè)計(jì)需要選擇合適的最優(yōu)翼型。由Pareto分布圖我們可以看出,翼型的升阻比越大,對(duì)應(yīng)流場(chǎng)的熵產(chǎn)率越大,變化趨勢(shì)與文獻(xiàn)[7]一致。因此提高翼型的升阻比,是以較大的流動(dòng)損失作為代價(jià)的。

圖5 算例二優(yōu)化過(guò)程種群分布情況Fig.5 Population distribution in optimization procedure in case 2

同時(shí),圖6也給出了四種優(yōu)化翼型的壓力云圖和熵產(chǎn)率云圖,四種優(yōu)化翼型分別為兩種算例中Pareto分布中對(duì)應(yīng)的最大、最小升阻比對(duì)應(yīng)翼型(分別命名為Opt_1、Opt_2、Opt_3和Opt_4,如圖所示)。由壓力云圖可知四種翼型升阻比越大,翼型下表面高壓區(qū)越大;由熵產(chǎn)率云圖可知翼型的流動(dòng)損失主要發(fā)生在翼型前緣、近壁面以及尾緣處,熵產(chǎn)云圖可以準(zhǔn)確反映出翼型流動(dòng)損失較大的部位,四種翼型熵產(chǎn)率越小,翼型尾緣處的流動(dòng)損失越小,而前緣和壁面處的損失幾乎不變。

圖6 優(yōu)化翼型目標(biāo)函數(shù)的Pareto分布Fig.6 Pareto curve of the airfoil shapes optimization

表3給出了對(duì)應(yīng)四種優(yōu)化翼型的氣動(dòng)性能對(duì)比,可知翼型優(yōu)化后升力系數(shù)得到很大程度提高,阻力系數(shù)亦有所增加。但升阻比相比于基準(zhǔn)翼型仍有很大提高。本文優(yōu)化結(jié)果與文獻(xiàn)[16]對(duì)RAE2822翼型進(jìn)行氣動(dòng)外形優(yōu)化設(shè)計(jì)后氣動(dòng)性能提高具有相似之處?；陟禺a(chǎn)理論的優(yōu)化翼型相比于傳統(tǒng)的基于最大升阻比和最小阻力系數(shù)的優(yōu)化翼型,能夠在控制流動(dòng)損失的前提下,得到升阻比盡可能大的翼型。優(yōu)化翼型Opt_1和Opt_2為基于熵產(chǎn)理論優(yōu)化得出,可以看出相比于基準(zhǔn)翼型,其升阻比分別提高了52.28%和24.13%。

圖7(a、b)分別為上述兩組優(yōu)化算例的四種優(yōu)化翼型壓力系數(shù)分布和幾何外形與初始翼型的對(duì)比。圖7(a)為兩組算例Pareto解中最大升阻比對(duì)應(yīng)的翼型Opt_1和Opt_3?？梢钥闯?優(yōu)化后的翼型彎度增加明顯,前緣半徑亦有所增加,該變化情況與對(duì)應(yīng)的升力系數(shù)增加明顯、阻力系數(shù)增大結(jié)果與文獻(xiàn)[17]具有一致性,翼型升阻比得到提高。翼型上表面厚度增加,有利于增加上表面流動(dòng)速度,從而降低上表面壓力,彎度增大,翼型下表面向內(nèi)凹進(jìn),便于產(chǎn)生更大升力。由壓力系數(shù)分布亦可以看出,上、下表面壓力系數(shù)明顯降低,使得升力系數(shù)明顯提高。雖然前緣半徑增大使得阻力系數(shù)增大,但升力系數(shù)增幅大于阻力系數(shù)增幅,使得升阻比相對(duì)于基準(zhǔn)翼型仍有所增大;圖7(b)為兩組算例Pareto解中最小升阻比對(duì)應(yīng)的翼型Opt_2和Opt_4,兩組翼型的厚度明顯變小,翼型下表面內(nèi)凹,使得下表面壓力系數(shù)變大。相比于NACA0012,兩組優(yōu)化翼型的前緣半徑減少。

表3 優(yōu)化翼型氣動(dòng)性能對(duì)比表Table 3 Comparison of aerodynamic characteristics of optimized airfoil

3 結(jié) 論

本文基于熵產(chǎn)理論對(duì)翼型的優(yōu)化設(shè)計(jì)進(jìn)行了一些探索,并與傳統(tǒng)的翼型優(yōu)化設(shè)計(jì)方法進(jìn)行了對(duì)比。其研究目的并不是要取代傳統(tǒng)的翼型優(yōu)化設(shè)計(jì)方法,而是基于熱力學(xué)第二定律提出一種新方法,作為傳統(tǒng)方法的一種補(bǔ)充。

通過(guò)參數(shù)化方法完成翼型參數(shù)化建模,并采用遺傳算法,最終對(duì)來(lái)流雷諾數(shù)為2.88×106,飛行迎角5°的翼型進(jìn)行多目標(biāo)優(yōu)化,得出如下結(jié)論:

(1)將優(yōu)化算法、CFD求解器、參數(shù)化方法和網(wǎng)格變形方法進(jìn)行了耦合,實(shí)現(xiàn)了二維翼型的多目標(biāo)優(yōu)化設(shè)計(jì),NSGA2優(yōu)化算法在小種群優(yōu)化中表現(xiàn)出良好的收斂性,能夠滿足復(fù)雜的非線性?xún)?yōu)化設(shè)計(jì)要求。

(2)CST翼型參數(shù)化方法有效減少了翼型設(shè)計(jì)變量的個(gè)數(shù),提高了優(yōu)化設(shè)計(jì)的效率,為未來(lái)三維復(fù)雜氣動(dòng)外形優(yōu)化設(shè)計(jì)提供了良好的參數(shù)化設(shè)計(jì)基礎(chǔ)。(3)相比于傳統(tǒng)的翼型氣動(dòng)優(yōu)化設(shè)計(jì)方法,基于熵產(chǎn)理論的氣動(dòng)優(yōu)化設(shè)計(jì)方法能夠具體掌握能量損失的大小及分布,得到氣動(dòng)性能更佳的翼型。

(4)通過(guò)對(duì)NACA0012翼型的多目標(biāo)優(yōu)化設(shè)計(jì)結(jié)果表明:與初始種群相比,最終優(yōu)化翼型在升阻比提高的條件下,流場(chǎng)熵產(chǎn)率減少,能量效率提高。