基于熵產(chǎn)理論設(shè)計方法的多目標翼型優(yōu)化

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

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

0 引 言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1.3 流場數(shù)值計算

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié) 論

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

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

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

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

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