孫蓬勃，周洲，李旭，王科雷

西北工業(yè)大學(xué) 航空學(xué)院，西安 710072

近年來，隨著環(huán)保要求和運(yùn)輸效率需求的提升，分布式電推進(jìn)技術(shù)受到了關(guān)注，并在飛行器的設(shè)計中逐步得到了應(yīng)用發(fā)展。以XV-24［1］和Lilium［2］為代表的分布式動力翼（Distribution-Propulsion-Wing，DPW）飛行器將涵道動力與機(jī)翼高度融合，形成了復(fù)雜的氣動-推進(jìn)耦合關(guān)系，已經(jīng)打破了傳統(tǒng)的飛發(fā)匹配概念。新型推進(jìn)方式的變革體現(xiàn)了飛行器設(shè)計從單一動力研究到飛行器總體、動力、氣動等技術(shù)集成研究的變化，為航空科技發(fā)展提供新的方向［3］。此類構(gòu)型的推進(jìn)單元本身成為升力產(chǎn)生面的一部分且動力裝置改變了機(jī)翼表面的流動，從而顯著影響整體的氣動力特性。這為分布式動力翼飛行器帶來了廣闊的設(shè)計空間，但也因為傳統(tǒng)的基于干凈來流的翼面設(shè)計方法失效，帶來了極大的設(shè)計困難。因此，深入開展分布式動力翼的氣動設(shè)計研究對于真正實現(xiàn)其優(yōu)勢具有重要的工程意義和學(xué)術(shù)價值，對此許多學(xué)者開展了研究。

Kerho［4-5］針對一種TeDP（Turboelectric Distributed Propulsion）驗證飛機(jī)開展了設(shè)計研究，對二維剖面外形的涵道唇口形狀和涵道安裝角進(jìn)行了設(shè)計研究，并進(jìn)一步開展了翼段三維設(shè)計與實驗研究。Schiltgen 等［6］和Papathakis［7］分 別對split-wing 這樣一種機(jī)翼完全由動力充當(dāng)?shù)臉?gòu)型進(jìn)行了數(shù)值模擬，結(jié)合進(jìn)出口寬度、唇口曲率等參數(shù)進(jìn)行了二維外形建模與設(shè)計分析；Lauer和Ansell［8-9］用貝塞爾曲線對分布式電推進(jìn)構(gòu)型的二維剖面進(jìn)行了參數(shù)化建模與優(yōu)化設(shè)計。夏濟(jì)宇等［10］結(jié)合二維理論模型與工程經(jīng)驗?zāi)Ｐ吞岢隽艘环N適用于矢量電推進(jìn)系統(tǒng)的氣動-推進(jìn)耦合模型。Guo 和Zhou［11］將分布式涵道單元內(nèi)外流分開考慮，基于超橢圓參數(shù)化方法對進(jìn)排氣道外形進(jìn)行了優(yōu)化推進(jìn)效率的設(shè)計。王科雷等［12］依次對動力翼二維剖面翼型、動力單元翼段、分布式動力翼整體的動力/氣動綜合特性進(jìn)行了數(shù)值研究，梳理了分布式動力翼內(nèi)外流耦合與其二維-單元-整體特性之間的內(nèi)在聯(lián)系。張星雨等［13］通過地面試驗結(jié)合數(shù)值模擬的方法對分布式電推進(jìn)技術(shù)驗證機(jī)的氣動性能及其氣動-推力耦合關(guān)系開展了研究。

飛行器設(shè)計通常需要滿足給定的升力和力矩特性。然而基于目前的研究可以看出，已有的動力翼研究多關(guān)注其自身的特性分析和優(yōu)化設(shè)計，但對如何將動力翼設(shè)計收斂到所需的目標(biāo)氣動特性的研究較少。如能在滿足氣動特性目標(biāo)值的前提下發(fā)揮動力翼的優(yōu)勢，則更具有實際工程意義。具體可從以下3 個方面展開討論：

1）從設(shè)計難度上看，相比常規(guī)機(jī)翼，動力翼的升力和俯仰力矩等氣動特性都有顯著改變；所受影響因素較多，不同情況下的氣動特性變化過程相對復(fù)雜。這增加了動力翼氣動設(shè)計的收斂難度，如何使動力翼設(shè)計方案能夠契合所給的設(shè)計目標(biāo)和約束成為需要研究的方向。

2）從設(shè)計思路上看，分布式動力翼具有典型的單元周期特性及動力壁面的類翼面特性。從文獻(xiàn)［11-12］中對動力翼特性規(guī)律的認(rèn)識以及設(shè)計邏輯的研究來看，可以根據(jù)內(nèi)外流與升力、推進(jìn)（阻力）特性之間關(guān)系對動力翼設(shè)計進(jìn)行拆分。取動力翼剖面形狀來進(jìn)行升力和力矩特性調(diào)控，取進(jìn)排氣道進(jìn)行推進(jìn)（阻力）特性調(diào)控，二者之間相互耦合，但又相對獨(dú)立。所以首先開展動力翼剖面二維設(shè)計研究是有重要作用的，其結(jié)果可直接指導(dǎo)動力翼外流面的形狀設(shè)計，并可進(jìn)一步作為進(jìn)排氣道三維形狀設(shè)計的基礎(chǔ)，以降低直接開展三維設(shè)計所需的計算迭代次數(shù)，從而提高動力翼的設(shè)計效率。此外由于其本身較快的計算和設(shè)計收斂速度，也可用于動力翼方案的快速對比評估。

3）從設(shè)計方法上看，多數(shù)對于分布式動力翼飛機(jī)的設(shè)計研究均基于傳統(tǒng)的參數(shù)優(yōu)化方法，在考慮發(fā)動機(jī)影響的條件下進(jìn)行優(yōu)化并收斂到目標(biāo)，面對較多的變量參數(shù)和復(fù)雜流動機(jī)理，這種方式難以準(zhǔn)確收斂到所需的氣動特性目標(biāo)值，且計算次數(shù)較多，計算量極大。所以相比正向的直接優(yōu)化逐步接近目標(biāo)氣動特性的方式，動力翼的反設(shè)計成為值得探索的方法。

因此面向分布式涵道動力機(jī)翼的設(shè)計，基于上述思路，本文通過一種以弦向環(huán)量分布為目標(biāo)的動力翼二維反設(shè)計方法開展了動力翼二維形狀的反設(shè)計研究，分析了相同升力和俯仰力矩前提下，涵道壁位置與涵道壁弦長對反設(shè)計外形及其氣動性能的影響，并進(jìn)一步開展了分布式動力翼剖面二維外形優(yōu)化反設(shè)計，以期采用反設(shè)計方法實現(xiàn)給定氣動特性目標(biāo)下分布式動力翼剖面形狀的快速設(shè)計，為動力翼的進(jìn)一步三維設(shè)計提供基礎(chǔ)，并用于方案的快速對比評估和優(yōu)化。

1 動力翼二維反設(shè)計方法

反設(shè)計方法通過求解空氣動力學(xué)反問題的方式來使設(shè)計外形達(dá)到期望的氣動特性，具有計算效率高、針對性強(qiáng)等特點［14］。然而，傳統(tǒng)的翼型反設(shè)計方法主要針對均勻遠(yuǎn)場來流，當(dāng)設(shè)計對象是動力翼剖面這種內(nèi)外流深度耦合的構(gòu)型時，仍需發(fā)展相應(yīng)的氣動反問題求解方法。升力和力矩特性是典型的翼型設(shè)計目標(biāo)，而如圖1 所示，圖中x為弦向位置，γ為環(huán)量，翼型的弦向環(huán)量分布（即弦向升力分布）是對升力和力矩的表征，所以可將其分解到動力翼剖面的各部分作為反設(shè)計目標(biāo)；在此基礎(chǔ)上通過參數(shù)化表達(dá)目標(biāo)環(huán)量分布和其他構(gòu)型參數(shù)，可以進(jìn)一步實現(xiàn)動力翼反設(shè)計方案的優(yōu)化。

圖1 分布式動力翼二維反設(shè)計Fig.1 2D inverse design of distribution-propulsion-wing

為了給出分布式動力翼反設(shè)計的完整流程并提高反設(shè)計效率，主要介紹以下兩點的推導(dǎo)和方法：

1）考慮槳盤等其他部件誘導(dǎo)影響的翼型理論推導(dǎo)與初步反設(shè)計流程。

2）反設(shè)計修正方法。

在反設(shè)計迭代的過程中，雖然是二維流場計算，但因迭代次數(shù)一般較多，直接采用CFD 方法仍然是比較耗時的。因此在迭代初期使用快速高效的勢流方法，結(jié)合反設(shè)計方程推導(dǎo)，有助于快速收斂到目標(biāo)解附近，在迭代后期再基于CFD計算進(jìn)行進(jìn)一步修正。

1.1 初步反設(shè)計

1.1.1 翼型理論推導(dǎo)

根據(jù)文獻(xiàn)［15］，對于二維薄翼型，其氣動特性可通過中弧線上的點渦環(huán)量分布來描述，其任意一點x0處的邊界條件為

式中：左側(cè)為來流在中弧線上當(dāng)?shù)匚恢玫姆ㄏ蚍至?，右?cè)為弦向各ξ處點渦對當(dāng)?shù)胤ㄏ蛘T導(dǎo)速度之和；V∞為來流速度；α為翼型迎角；y′(x0)為當(dāng)?shù)刂谢【€斜率；c為翼型弦長。

對于動力翼二維剖面當(dāng)中的一個翼型，將槳盤及其他部件的誘導(dǎo)速度考慮進(jìn)來作為此翼型背景來流的一部分，在中弧線任意一點有疊加了槳盤等其他部件誘導(dǎo)速度后的局部流動迎角αlocal(x)及局部速度Vlocal(x)。

進(jìn)一步將弦向環(huán)量γ分布由傅氏級數(shù)表達(dá)：

式中：θ為弦向各處坐標(biāo)ξ的三角函數(shù)變換，相應(yīng)的坐標(biāo)x0對應(yīng)為Θ；A0、An為系數(shù)。

將式（2）代入式（1）右側(cè)，可得

通過式（3）可以看出，翼型沿弦向的環(huán)量γ分布可視為與中弧線斜率y′和當(dāng)?shù)貋砹魉俣萔local相關(guān)的函數(shù)。如果在翼型上，從前緣到后緣，布置N個控制點，通過給定在各控制點的目標(biāo)環(huán)量值γi(i=1，2，…，N)并計算流場得到各控制點處的當(dāng)?shù)厮俣萔local，將中弧線斜率值作為求解對象，可在任意控制點Θ處將式（3）代入式（2）：

將式（4）中的積分離散化，提取出y′的系數(shù)并將各控制點處構(gòu)造的方程組寫為矩陣形式：

式中：θi為各點的三角函數(shù)坐標(biāo)；(Δθ)j為第j處θ區(qū)間量Δθ的大小。

根據(jù)式（6）的推導(dǎo)，可以看到式（5）中矩陣A是奇異的（最后一行對應(yīng)后緣點處的等式，此時θi=π，所以此行值均為0），因而構(gòu)建方程組時應(yīng)去除掉過于接近后緣的幾個點，考慮到一般翼型彎度變化的連續(xù)性，這不會產(chǎn)生過大的誤差，反而帶來了求解上的便利。采用偽逆法求解式（6），便可以得到槳盤等其他部件影響下的薄翼型中弧線形狀。

1.1.2 翼型厚度影響分解

實際翼型設(shè)計過程中，厚度帶來的影響經(jīng)常是不可忽略的。文獻(xiàn)［16］中，假設(shè)流動是不可壓和有勢的，一個有厚度翼型所對應(yīng)的環(huán)量分布可分解為如圖2 所示的兩部分。Part1 表示具有對應(yīng)厚度分布的對稱翼型在該迎角下的環(huán)量分布，Part2 表示對應(yīng)彎度分布的薄翼型在0°迎角下的環(huán)量分布。所以需首先計算Part1 中對稱翼型的弦向環(huán)量分布，并在反求Part2 薄翼型彎度分布時，在目標(biāo)環(huán)量分布的給定中將這一部分去掉。

圖2 翼型厚度影響分解Fig.2 Decomposition of airfoil thickness influence

但是針對所研究的對象，由于存在槳盤等其他部件的誘導(dǎo)速度，翼型表面各處的背景來流是不同的。所以在計算Part1 對稱翼型弦向環(huán)量分布時，采用二維渦格法［16］把實際翼型表面各處的來流速度矢量轉(zhuǎn)移給對稱翼型上對應(yīng)表面位置的渦格，作為其速度邊界條件（圖3）。

圖3 對稱翼型的速度邊界條件Fig.3 Velocity boundary conditions of symmetrical airfoil

1.1.3 二維自由尾跡渦格法

為了能夠?qū)Ψ植际絼恿σ順?gòu)型進(jìn)行快速計算迭代并獲取槳盤等其他部件對所求翼面各處的誘導(dǎo)速度，采用自由尾跡渦格法求解。李旭等［17］參考Shollenberger［18-19］的求解思路，基于自由流線模型，利用渦面對射流進(jìn)行建模，結(jié)合一階渦格法，建立了二維自由尾跡渦格計算方法。如圖4 所示，以離散渦面加半無限長渦面表征二維槳盤噴流，其中離散段用來迭代受物面影響的自由尾跡；同時也將翼型離散為渦面。詳細(xì)的迭代計算過程可參考文獻(xiàn)［17］。

圖4 翼面與槳盤噴流的渦面Fig.4 Vortex surface of airfoil and disk jet

1.1.4 初步反設(shè)計流程

圖5 顯示了分布式動力翼二維剖面初步反設(shè)計流程，主要包括以下6 個步驟：

圖5 初步反設(shè)計流程Fig.5 Primary process of inverse design

1）給出槳盤的推力、半徑、弦向位置，翼型弦長，來流速度，并給出各翼面目標(biāo)環(huán)量弦向分布。初始的翼型中弧線形狀可以是無安裝角的直線，兩翼型之間保持一定高度差以容納槳盤。

2）初始的槳盤噴流可以是平直的直到無限遠(yuǎn)，通過渦格法計算翼型中弧線上各處所受的誘導(dǎo)速度，并通過式（5）計算得到翼型中弧線形狀。

3）基于給定的翼型厚度分布，得到翼型形狀。由于翼型弦向環(huán)量分布是由傅里葉級數(shù)擬合得到的，所以當(dāng)階數(shù)給定不合適的時候，容易出現(xiàn)過擬合的問題，表現(xiàn)為求解得到的翼型形狀的振蕩（圖6），這會嚴(yán)重影響到后續(xù)求解的收斂性；所以此處預(yù)先通過CST（Class Shape Function Transformation）擬合的方式［20-21］對翼型進(jìn)行了光順，從而提高設(shè)計程序的魯棒性；光順方法為分別光順上下表面，求均值得到新的中弧線，再通過給定的厚度分布生成新的翼型外形。進(jìn)一步垂直平移翼面，使得槳盤能被容納在給定的位置。

圖6 CST 光順翼型Fig.6 CST airfoil smoothing

4）結(jié)合噴口內(nèi)壁形狀更新槳盤噴流離散渦面的初始形狀，并通過自由尾跡渦格法迭代求解得到翼型表面各處的當(dāng)?shù)貋砹魉俣仁噶俊?/p>

5）由給定的厚度分布生成對稱翼型，將步驟4）中的翼型表面各處局部速度矢量賦給對稱翼型作為邊界條件，由渦格法計算對稱翼型表面壓力分布，進(jìn)而得到對稱翼型的弦向環(huán)量分布。

6）將對稱翼型的弦向環(huán)量分布從目標(biāo)環(huán)量分布中去除，通過式（5）計算更新翼型中弧線形狀。返回步驟3），直到相近兩次迭代得到的翼型中弧線形狀均方誤差小于10-5，則認(rèn)為初步反設(shè)計收斂。

1.2 修正方法

1.2.1 CFD 方法驗證

使用商業(yè)數(shù)值模擬軟件FLUENT 進(jìn)行CFD計算，采用動量源方法（Momentum Source Method，MSM）模擬槳盤推力對來流的加速作用，采用k-ω剪切應(yīng)力傳輸（Shear Stress Transport，SST）湍流模型，求解RANS（Reynolds-Averaged Navier-Stokes）方程，數(shù)值計算過程中空間離散采用二階迎風(fēng)MUSCL（Monotone Upstreamcentered Scheme for Conservation Laws）插值的Roe 格式，時間離散與推進(jìn)則采用隱式AF（Approximate Factorization）方法。

首先針對NASA 涵道風(fēng)扇［22］采用動量源方法進(jìn)行計算分析并與實驗結(jié)果進(jìn)行對比。如圖7所示，將實驗中的槳盤拉力值以對應(yīng)大小的動量源的形式代替，最后比對30 m/s 來流工況下的總拉力T，邊界條件設(shè)置與實驗條件相一致，結(jié)果見表1，可見在不同網(wǎng)格數(shù)量下，MSM 計算結(jié)果與實驗值相對誤差η最大為3.37%，最小為1.28%。通過表1 結(jié)果對比，能夠排除網(wǎng)格數(shù)量對計算結(jié)果的影響，并且說明動量源方法的有效性。

表1 MSM 不同數(shù)量網(wǎng)格計算結(jié)果Table 1 Calculation results of MSM with different grid numbers

圖7 MSM 計算網(wǎng)格Fig.7 Mesh of MSM calculation

另外，也針對第2～4 節(jié)中的NLF（1）-0416 翼型經(jīng)過網(wǎng)格無關(guān)性驗證后給出二維CFD 計算結(jié)果，總網(wǎng)格數(shù)約為12 萬。計算工況為來流馬赫數(shù)0.1、雷諾數(shù)4×106、迎角1°［23］。翼型表面壓力分布如圖8 所示，可見與實驗結(jié)果基本一致。

圖8 NLF（1）-0416 二維計算網(wǎng)格與Cp分布Fig.8 2D mesh of NLF（1）-0416 and Cp distribution

通過對表1 和圖8 的分析，能夠說明所用CFD 方法的準(zhǔn)確性。

1.2.2 CFD 修正

由于在初步反設(shè)計流程中采用的是快速勢流方法，無法考慮黏性效應(yīng)且準(zhǔn)確性受渦核半徑等參數(shù)設(shè)置的影響，與實際流動存在一定偏差，所以還需通過CFD 計算來對反設(shè)計進(jìn)行校核修正。如圖9 所示，在初步的反設(shè)計迭代收斂后，對設(shè)計結(jié)果進(jìn)行二維CFD 計算，得到實際的弦向環(huán)量分布（圖9（a）中的初步反設(shè)計）并與目標(biāo)弦向環(huán)量分布進(jìn)行對比，通過將差量疊加修正進(jìn)初步反設(shè)計的設(shè)計目標(biāo)，反復(fù)進(jìn)行若干次初步反設(shè)計的方式（圖9（b）），使得初步反設(shè)計的實際結(jié)果最終與目標(biāo)弦向環(huán)量分布接近一致。

圖9 CFD 修正示例Fig.9 Example of CFD correction

1.2.3 Modified Garabedia-Mcfadden（MGM）修正

由于傅里葉級數(shù)擬合是對弦向環(huán)量分布的整體性描述，因此經(jīng)過1.1 節(jié)和1.2.2 節(jié)迭代得到的設(shè)計結(jié)果，某些時候在局部環(huán)量梯度較大的位置（尤其是槳盤附近和前緣）或所取階數(shù)n不合適時，仍可能與目標(biāo)環(huán)量分布存在少量差距，這時可使用MGM 對設(shè)計結(jié)果做進(jìn)一步微調(diào)。

Malone 等［24］基于線性小擾動理論，改進(jìn)得到了MGM 方法用于翼型反設(shè)計；其思想是將翼型表面壓力系數(shù)與氣動外形參數(shù)結(jié)合，得到計算、目標(biāo)壓力系數(shù)與翼型表面坐標(biāo)相關(guān)聯(lián)的MGM 方程。但由于設(shè)計對象是翼型中弧線形狀且目標(biāo)為環(huán)量分布，翼型表面目標(biāo)壓力分布是未知的，所以給出保持厚度分布、只修正中弧線形狀的MGM 方程（環(huán)量等效于翼型上下表面壓強(qiáng)差，所以相當(dāng)于上下表面各自的壓力分布修正疊加，因而公式形式與常規(guī)MGM 方程一致）：

式中：R為殘差，R=γt-γc，其中γc和γt分別為計算環(huán)量值與目標(biāo)環(huán)量值；x為翼型弦向坐標(biāo)；F0、F1、F2為常值系數(shù)，用以實現(xiàn)穩(wěn)定迭代，它們后面的項分別代表中弧線各點高度位置的修正量ΔY及其對x的一階和二階導(dǎo)數(shù)?；谖墨I(xiàn)［23］中的求解方法對中弧線進(jìn)行迭代微調(diào)，使動力翼環(huán)量分布與目標(biāo)分布的平均相對誤差（Average Relative Error，ARE，進(jìn)一步下降，直到收斂。第2～4 節(jié)在反設(shè)計研究中，權(quán)衡氣動特性的收斂精度和迭代計算次數(shù)，取收斂標(biāo)準(zhǔn)為ARE ＜0.06。

2 反設(shè)計示例

驗證第1 節(jié)介紹的反設(shè)計方法的有效性。為盡量保證所給目標(biāo)分布的合理性，同時也為研究相同升力和力矩特性下動力翼相比傳統(tǒng)翼型的優(yōu)勢并進(jìn)一步結(jié)合研究背景，以弦長1.5 m 的通航飛機(jī)典型研究翼型NLF（1）-0416（圖10）的弦向環(huán)量分布為目標(biāo)，反設(shè)計與其物面升力和力矩一致的動力翼二維外形。給出動力翼單元剖面中wing1（下翼面）、wing2（上涵道壁）、槳盤的弦向位置和長度，如圖11 所示（(xd-x1)/(x2-x1)=0.3）。動力翼剖面的兩翼型厚度分布形式與NACA 四位數(shù)翼型相一致，最大相對厚度均為8%。槳盤半徑r=0.15 m，推力T′=400 N（T′為二維等效推力，其與三維圓槳盤單元推力之間的關(guān)系詳見文獻(xiàn)［4］）。遠(yuǎn)場來流速度V∞=60 m/s，大氣參數(shù)按海拔1 000 m 給定，密度ρ=1.111 kg/m3，迎角α=0°。

圖10 NLF（1）-0416（弦長為1.5 m）Fig.10 NLF（1）-0416（chord is 1.5 m）

圖11 部件位置Fig.11 Part location

如無特別說明，氣動力系數(shù)的參考長度均為總弦長。此外，實際上動力翼俯仰力矩除受環(huán)量分布影響，也與參考點和兩翼型物面阻力矢量及槳盤推力矢量的位置關(guān)系有關(guān)，把評估俯仰力矩時的參考點放在槳盤中心線上，盡量減少力矩對比時槳盤推力線和壁面阻力對俯仰力矩的影響；參考點弦向位置在1/4 總弦長處。

對于如何給定動力翼剖面翼型的目標(biāo)環(huán)量分布形式，參考文獻(xiàn)［24］對翼型表面壓力分布的描述思路來進(jìn)行。文獻(xiàn)［24］將翼型上/下表面的壓力分布形式按激波階躍位置分為前加速段和后減速段。類似的，可將動力翼上部涵道壁的翼型環(huán)量分布按槳盤處的環(huán)量階躍分為前部槳盤抽吸段和后部恢復(fù)段兩部分。前部環(huán)量分布由于受翼型本身形狀變化和槳盤抽吸的影響，可能存在多個極值點，所以通過多項式的形式描述；而考慮到槳盤后的內(nèi)流面單側(cè)壓強(qiáng)較高，且翼型內(nèi)外壓強(qiáng)差（環(huán)量）向后緣一般單調(diào)恢復(fù)至0，所以參考文獻(xiàn)[24]指數(shù)函數(shù)的方式描述后部恢復(fù)段的環(huán)量分布形式。需要說明的是這種描述方法不一定能涵蓋所有環(huán)量分布形式，這里僅為一類較符合實際的形式，用以開展反設(shè)計方法的驗證工作。在實際設(shè)計時，可根據(jù)具體情況和經(jīng)驗制定不同的目標(biāo)環(huán)量分布形式。

結(jié)合上述討論，在反設(shè)計示例中，wing1 和wing2 的弦向目標(biāo)環(huán)量分布按如下方式給出：首先通過CFD 計算得到基準(zhǔn)翼型的弦向環(huán)量分布，然后將上涵道壁wing2 環(huán)量分布形式表示為

式中：γd為槳盤前處的wing2 當(dāng)?shù)丨h(huán)量值，考慮到槳盤抽吸一般會使wing2 前部下表面壓力低于外表面，所以γd一般為負(fù)數(shù)，取γd=-15 m2/s；Δγ=T/(2rρV∞)為槳盤前后壓力階躍所對應(yīng)的環(huán)量值階躍；n1、n2、n3為分布控制系數(shù)，在對比不同系數(shù)取值設(shè)計結(jié)果的基礎(chǔ)上，可得到相對合理的系數(shù)取值范圍，取n1=0.7、n2=5、n3=5。將目標(biāo)翼型與wing2 的目標(biāo)環(huán)量分布相減，得到下部翼面wing1 的目標(biāo)環(huán)量分布。目標(biāo)翼型及其分解為wing1 和wing2 的弦向環(huán)量分布如圖12 所示。

圖12 目標(biāo)環(huán)量分布Fig.12 Target circulation distributions

圖13 給出了經(jīng)過8 輪CFD 修正初步反設(shè)計之后的環(huán)量分布，可見反設(shè)計得到的弦向環(huán)量分布與給定的目標(biāo)基本接近，但在槳盤和前緣附近仍有一些差距，這使得此時的環(huán)量分布平均相對誤差仍比較高。進(jìn)一步使用MGM 方法進(jìn)行修正，F(xiàn)0、F1、F2分別取值為4、2、2。經(jīng)過另外8 次結(jié)合CFD 計算的MGM 修正后，最終設(shè)計得到的環(huán)量分布與外形如圖14 所示，圖中CL為升力系數(shù)，Cm為俯仰力矩系數(shù)?？梢娊?jīng)過外形微調(diào)之后，最終設(shè)計的到的環(huán)量分布與目標(biāo)吻合良好，平均相對誤差為0.058 7，物面升力系數(shù)和俯仰力矩系數(shù)與目標(biāo)翼型一致。

圖13 CFD 修正設(shè)計后的環(huán)量分布Fig.13 Circulation distributions after CFD correction

圖14 最終設(shè)計結(jié)果Fig.14 Final design results

3 不同動力翼外形參數(shù)分析

在驗證了反設(shè)計方法的有效性之后，針對動力翼的兩個主要設(shè)計參數(shù)（上涵道壁的弦長Cduct=x2-x1和前后位置x2）對反設(shè)計結(jié)果的影響進(jìn)行分析；同時也與給定總環(huán)量分布時所參考的目標(biāo)翼型進(jìn)行對比，以研究動力翼相比傳統(tǒng)翼型的區(qū)別和優(yōu)勢。除所分析的參數(shù)外，來流和其余動力翼參數(shù)取值與第2 節(jié)一致。

3.1 涵道壁位置影響

計算分析了涵道壁后緣位置x2=1.5，1.3，1.1 時，不同動力翼反設(shè)計結(jié)果的升力、阻力和俯仰力矩隨迎角變化的情況及與目標(biāo)翼型NLF（1）-0416 的對比。各翼面目標(biāo)環(huán)量分布如圖15所示，反設(shè)計結(jié)果如圖16 所示（圖中俯仰力矩參考點取在（0.375，0）處，即1/4 弦長位置）。

圖15 不同x2 反設(shè)計目標(biāo)環(huán)量分布Fig.15 Target circulation distributions of inverse design with different x2

圖16 不同x2 反設(shè)計結(jié)果Fig.16 Inverse design results of different x2

可見在相同的總環(huán)量分布下，在0°迎角設(shè)計點時，各動力翼方案的升力系數(shù)CL和俯仰力矩系數(shù)Cm與目標(biāo)翼型基本一致；隨著上部涵道壁前移，其安裝角逐漸減小，下翼頭部也隨之下偏；相比翼型NLF（1）-0416，各動力翼方案的升力系數(shù)隨迎角變化的斜率CLα基本一致且均偏高，較大迎角時的抗分離能力提升；各動力翼方案物面阻力系數(shù)CD（動力翼剖面外形的壁面阻力之和，不包括槳盤推力）均小于目標(biāo)翼型，涵道壁靠前時物面阻力較小；隨著上部涵道壁前移，俯仰力矩系數(shù)隨迎角變化的斜率Cmα由負(fù)到正。

為研究圖15 中變化的原因，給出各動力翼方案在0°和12°迎角時的表面壓強(qiáng)分布、12°迎角時的流場速度分布以及6°迎角時的弦向總環(huán)量分布，如圖17～圖20 所示。

圖17 不同x2時0°迎角表面壓強(qiáng)分布Fig.17 Surface pressure distributions at α=0° with different x2

圖18 不同x2時12°迎角流場速度分布Fig.18 Velocity distributions of flow fields at α=12° with different x2

圖19 不同x2時12°迎角表面壓強(qiáng)分布Fig.19 Surface pressure distributions at α=12° with different x2

圖20 不同x2時6°迎角環(huán)量分布Fig.20 Circulation distributions at α=6° with different x2

從環(huán)量分布看，下翼面環(huán)量為正值，其對遠(yuǎn)場來流有順時針誘導(dǎo)偏折的效應(yīng)，當(dāng)上涵道壁越靠前緣，當(dāng)?shù)貋砹饔窍鄬υ酱螅湓谠O(shè)計點的目標(biāo)環(huán)量（升力）分布形式不變，因而其安裝角越小。從壓力分布看，槳盤對入流有抽吸加速的效果，當(dāng)槳盤位置跟隨上涵道壁前移，其對下翼前部附近的抽吸增升也會增加，但下翼前部的目標(biāo)環(huán)量分布不變，所以其前部會對應(yīng)下偏以抵消升力增量；在12°迎角時，目標(biāo)翼型的流動分離出現(xiàn)在上表面后緣，而對于動力翼，其下翼后緣處在噴流影響中，上涵道壁的后緣曲率也相對較小，所以動力翼大迎角抗分離能力相對目標(biāo)翼型有所提高；動力翼阻力系數(shù)的下降也是由于槳盤抽吸提高了動力翼尤其是下翼的前緣吸力，因而上涵道壁及槳盤越靠前時阻力越小。從6°迎角環(huán)量分布看，隨著迎角高于0°設(shè)計點，動力翼兩個翼面都有顯著的前緣吸力峰，且強(qiáng)度大于目標(biāo)翼型，因而CLα增加；隨著上部涵道壁的前移，其環(huán)量（升力）增量的作用位置也前移，使得動力翼整體焦點前移，引起Cmα由負(fù)到正的轉(zhuǎn)變。

3.2 涵道壁弦長影響

計算分析了涵道壁弦長Cduct=1.4，1.0，0.6時，不同動力翼反設(shè)計結(jié)果的升力、阻力和俯仰力矩隨迎角變化的情況及與目標(biāo)翼型NLF（1）-0416的對比。反設(shè)計目標(biāo)環(huán)量分布如圖21 所示，結(jié)果曲線如圖22 所示。可見在相同總環(huán)量目標(biāo)弦向分布下，在0°迎角設(shè)計點時，各動力翼方案的升力系數(shù)CL和俯仰力矩系數(shù)Cm與目標(biāo)翼型基本一致；隨著上部涵道壁弦長增加，其安裝角下降，并且因為相對厚度參數(shù)不變，其整體厚度增加；相比翼型NLF（1）-0416，各動力翼方案的升力系數(shù)隨迎角變化斜率CLα仍然均偏高，其中Cduct=1.4 時 最高；中小迎角時，各動力翼方案物面阻力系數(shù)CD仍均小于目標(biāo)翼型，涵道壁弦長較長時，其厚度增加，但物面阻力仍相對較小；隨著上部涵道壁弦長增加，俯仰力矩系數(shù)曲線斜率Cmα整體仍由負(fù)到正變化。

圖21 不同Cduct 反設(shè)計目標(biāo)環(huán)量分布Fig.21 Target circulation distributions of inverse design with different Cduct

圖22 不同Cduct 反設(shè)計結(jié)果Fig.22 Inversed design results of different Cduct

為研究圖22 中現(xiàn)象出現(xiàn)的原因，同樣給出不同動力翼方案在0°和12°迎角時的表面壓力分布、12°迎角時的流場速度分布及6°迎角時的總環(huán)量弦向分布，如圖23～圖26 所示。

圖23 不同Cduct時0°迎角表面壓強(qiáng)分布Fig.23 Surface pressure distributions at α=0° with different Cduct

圖24 不同Cduct時12°迎角流場速度分布Fig.24 Velocity distributions of flow fields at α=12° with different Cduct

圖25 不同Cduct時12°迎角表面壓強(qiáng)分布Fig.25 Surface pressure distributions at α=12° with different Cduct

圖26 不同Cduct時6°迎角環(huán)量分布Fig.26 Circulation distributions at α=6° with different Cduct

從目標(biāo)環(huán)量分布、流場與壓力分布看，上涵道壁安裝角變化和下翼頭部下偏原因與3.1 節(jié)分析類似。此外在較大迎角時，目標(biāo)翼型的流動分離出現(xiàn)在上表面后緣，而對于動力翼來說，其下翼后緣處在噴流影響中，上涵道壁則是在靠近后部且弦長較短時受下翼影響其局部迎角較小，弦長較長時前緣半徑較大且安裝角較小，所以動力翼大迎角抗分離能力相對目標(biāo)翼型均有所提高；動力翼阻力系數(shù)的下降也是由于槳盤抽吸提高了動力翼尤其是下翼前緣吸力并且強(qiáng)于翼面變厚增加的形阻。從6°迎角環(huán)量分布看，隨著迎角高于0°設(shè)計點，動力翼兩個翼面都有顯著的前緣吸力峰，且強(qiáng)度大于目標(biāo)翼型，因而CLα增加；隨著上部涵道壁弦長增加，因其翼面尺度變大，其前緣的吸力峰強(qiáng)度和范圍顯著提高，使得動力翼整體升力系數(shù)曲線斜率變大；另外其Cmα由負(fù)到正的轉(zhuǎn)變原因也與3.1 節(jié)類似。

綜合3.2 節(jié)和3.1 節(jié)反設(shè)計所得到的外形來看，其下翼形狀有時會出現(xiàn)中部的凹陷，與常規(guī)翼型有所區(qū)別。這是因為進(jìn)行的是分布涵道機(jī)翼的剖面二維外形的反設(shè)計，為滿足所給特定的目標(biāo)環(huán)量分布，有時會出現(xiàn)這種外形結(jié)果；并且對于這種涵道風(fēng)扇與翼面耦合的構(gòu)型，相當(dāng)于涵道進(jìn)口附近的下表面較低，這種情況是正常的；根據(jù)后續(xù)多迎角的分析，其也基本沒有產(chǎn)生流動分離等不良影響。實際設(shè)計過程中，可從多學(xué)科的角度對設(shè)計結(jié)果進(jìn)行分析取舍，并對應(yīng)調(diào)整和限制目標(biāo)環(huán)量分布的形式、涵道壁相對位置和長度等參數(shù)，使最終的剖面形狀符合氣動、結(jié)構(gòu)等各學(xué)科的期望。

4 動力翼優(yōu)化設(shè)計

在驗證了反設(shè)計方法的有效性及對兩個涵道壁參數(shù)的影響進(jìn)行對比分析之后，繼續(xù)開展固定總弦向環(huán)量分布下的動力翼二維優(yōu)化設(shè)計，以期實現(xiàn)具有最小物面阻力系數(shù)（即除槳盤外其他部件的阻力之和）CD的設(shè)計結(jié)果；此外，考慮到實際轉(zhuǎn)子槳盤應(yīng)工作在盡可能均勻的來流條件下以降低畸變帶來的功率損失及涵道進(jìn)氣段應(yīng)使總壓損失盡可能小，所以在二維優(yōu)化設(shè)計的過程中，使槳盤前20 mm 處的總壓恢復(fù)系數(shù)σ和速度均勻指數(shù)U最小，也同樣作為優(yōu)化目標(biāo)之一。

總壓恢復(fù)系數(shù)σ和速度均勻指數(shù)U定義為［1，25］

式中：為來流總壓；為各網(wǎng)格處的總壓；Ai為各網(wǎng)格面積，二維情況下以截面上的網(wǎng)格高度表征；qi為流量函數(shù)。

式中：Vi為各網(wǎng)格上的速度大?。粸闃P處平均速度大小。

優(yōu)化案例同樣以第2 節(jié)中的弦向總環(huán)量分布為優(yōu)化過程中每個方案的反設(shè)計目標(biāo)，以第2節(jié)的動力翼形狀為優(yōu)化基準(zhǔn)，同時固定相同的來流狀態(tài)和槳盤推力；優(yōu)化過程中的變量為涵道壁后緣弦向位置x2以及翼型厚度分布參數(shù)，其他參數(shù)固定為第2 節(jié)中的值。翼型厚度分布Z(x)通過6 階CST 方法進(jìn)行參數(shù)化：

式中：dz為后緣厚度，厚度分布涉及a0～a6共7 個控制變量。優(yōu)化工作的數(shù)學(xué)表達(dá)式為

3個優(yōu)化目標(biāo)的權(quán)值比例為0.50∶0.25∶0.25。

優(yōu)化流程如圖27 所示：首先使用最優(yōu)拉丁超立方方法生成100 個初始樣本點，建立Kriging 代理模型以加快迭代速度；然后使用多島遺傳算法進(jìn)行尋優(yōu)，使用尋優(yōu)結(jié)果變量進(jìn)行反設(shè)計校核，并更新代理模型數(shù)據(jù)庫，以此循環(huán)；最終得到全局優(yōu)化設(shè)計結(jié)果。優(yōu)化算法主要參數(shù)設(shè)置為：①亞種群規(guī)模為20；②島數(shù)為10；③進(jìn)化代數(shù)為10；④交叉、變異和遷移概率分別為1.0、0.01、0.01；⑤ 遷移間隔為5。

圖27 優(yōu)化框架Fig.27 Optimization framework

圖28 顯示了初始方案與優(yōu)化之后的方案對比。相比初始方案，優(yōu)化后的動力翼剖面形狀在設(shè)計點升力和力矩特性不變的前提下，物面厚度有所增加，上涵道壁前移，物面阻力系數(shù)下降了160%；同時槳盤前入流的總壓損失很小，總壓恢復(fù)系數(shù)基本維持不變；此外速度分布的均勻性進(jìn)一步提升。阻力進(jìn)一步下降的主要原因與第3 節(jié)分析類似，圖29 為動力翼表面壓力分布對比，可見上涵道壁帶動槳盤位置前移，增強(qiáng)了對動力翼前部壁面的抽吸，優(yōu)化后的上涵道壁前部整體的壓力水平較低，進(jìn)一步降低了前后壓差阻力。而通過調(diào)整厚度分布，使槳盤入流的均勻性和總壓恢復(fù)得到了保持，較厚的動力翼剖面形狀也有利于結(jié)構(gòu)布置。

圖28 優(yōu)化結(jié)果Fig.28 Optimization results

圖29 優(yōu)化結(jié)果Cp分布Fig.29 Cp distributions of optimization results

5 結(jié)論

基于一種以弦向環(huán)量分布為目標(biāo)的動力翼二維反設(shè)計方法，對升力和俯仰力矩相同的動力翼反設(shè)計方案進(jìn)行了涵道壁參數(shù)影響分析和進(jìn)一步的減阻優(yōu)化設(shè)計，主要有以下結(jié)論：

1）所使用的動力翼二維反設(shè)計方法通過基于翼型理論推導(dǎo)的初步反設(shè)計、CFD 修正和MGM 修正方法，能夠快速收斂到目標(biāo)環(huán)量分布附近，平均相對誤差為0.058 7。

2）隨著動力翼上涵道壁位置前移，上涵道壁安裝角降低，槳盤跟隨前移并對下翼前部的抽吸影響增加，使物面阻力下降、下翼前部跟隨下偏；此外動力翼焦點前移，Cmα由負(fù)到正變化；迎角大于0°設(shè)計點時，動力翼前緣吸力峰比目標(biāo)翼型更高，所 以CLα較高，并且在12°迎角時仍未出現(xiàn)分離。

3）隨著動力翼上涵道壁弦長增加，上涵道壁安裝角降低，同樣有槳盤跟隨前移使其對下翼前部抽吸影響增加，使下翼前部下偏且槳盤抽吸的影響強(qiáng)于涵道壁厚度增加，使整體阻力下降；另外動力翼焦點前移，Cmα由負(fù)到正變化；動力翼相比目標(biāo)翼型的CLα均較高，上涵道壁弦長較長時其前緣吸力峰強(qiáng)度和范圍較大，CLα最高。

4）相比動力翼的初始形狀，優(yōu)化后的動力翼涵道壁前移，厚度有所增加，阻力系數(shù)下降了160%，同時槳盤入流總壓基本沒有損失，速度分布均勻性進(jìn)一步增加。