賀紅林，周 翔，朱保利，余春錦

（南昌航空大學 航空制造工程學院，江西 南昌330063）

0 引 言

撲翼飛行器（FMAV）集舉升、懸停、推進功能于其機翼，具有尺寸小、質(zhì)量輕、機動性好、效率高等特點，在軍事與民用上具有廣泛的應用前景。因其在氣動與穩(wěn)定性等方面存在明顯優(yōu)勢，F(xiàn)MAV已成為微型飛行器的研究熱點。近些年，適于高性能FMAV研制的需要，對其氣動問題進行了不少研究，如Smith［1］和 Vest［2］采用面元法研究了飛鳥氣動特性；Sun［3］通過對N－S方程的數(shù)值求解并基于渦動力學理論闡釋了撲翼的非定常氣動特性；Joseph［4］采用非定常渦格法（UVLM）模擬柔性撲翼流場；Lin用MOC法求解Euler／Navier－Stokes方程以模擬撲翼運動［5］。在這些方法中，面元法簡單但計算結(jié)果較粗糙；CFD法的結(jié)果較準確，但計算開銷巨大，只適于剛性撲翼［6］；UVLM快速高效且可對流－構耦合態(tài)的柔性撲翼上的渦分布、尾渦變化等多種因素進行建模計算，并能獲得較準確氣動結(jié)果［7］?？偟膩砜矗延械难芯恐饕杏趽湟砹鲌龅哪M，對于撲翼的選型和設計支持不夠。事實上，實用中更關注的是撲翼結(jié)構及其撲動參數(shù)確定問題。據(jù)此，旨在為撲翼選型與設計提供一定依據(jù)，本文探索FMAV的氣動優(yōu)化，并將在對FMAV進行氣動分析及其UVLM建模的基礎上，引入GPU（Graphic Processing Unit）技術實現(xiàn)該算法，以滿足 優(yōu)化迭代中海量計算的需要，同時，考慮到UVLM模型不能提供導數(shù)信息，還將引入模式搜索（PS）法作為FMAV的尋優(yōu)工具。

1 非定常渦格法建模

1.1 翼面速度表示

為了表征撲翼的運動，在其上定義慣性坐標系O（X，Y，Z）和機翼坐標系o（x，y，z），前者固定不動，后者隨機翼運動，并設兩者在初始時刻位置重合。若在t（＞0）時刻，機翼坐標系在慣性系中位姿為R0＝（X0，Y0，Z0），D0＝（α，β，γ），則機翼系內(nèi)的翼面任意點A（x，y，z）在慣性系內(nèi)的速度

1.2 撲翼流場邊界條件

機翼撲動時，除前緣及翼尾的部分區(qū)域外，其它各處流場為不可壓無旋流，為此定義速度勢函數(shù)Φ（X，Y，Z）并滿足：▽2Φ＝0。因撲翼對遠處流場影響很小，因此其遠場邊界條件為Φ＝0。而近場邊界條件則滿足無穿透條件，即

式中，VS為翼面相對運動速度，V0為來流速度，VT、Ω分別為翼面平動和轉(zhuǎn)動速度，n為翼面單位法矢。

1.3 尾渦模型

撲翼的氣動狀態(tài)主要取決其尾渦。根據(jù)Kelven定律，當撲翼運行時其流場環(huán)量之和恒為常數(shù)。這意味著流場總環(huán)量的任何變化都會在尾流區(qū)引起等值、符號相反的渦量變化。該變化將發(fā)展成翼后緣沿角平分線處發(fā)散出的一段渦線。當翼面環(huán)量隨翼撲動而變化時，翼后緣將拖出一列連續(xù)渦帶。可按間距V（t）Δt對渦帶進行離散，且只需選取適當?shù)臅r間間隔Δt，就可用一組離散尾渦環(huán)模擬該渦帶。各離散渦環(huán)量應等于相應時段從翼后緣所脫離渦帶之渦強的積分，即

這樣，便得到圖1所示撲翼UVLM模型。若第i時段內(nèi)，環(huán)量Γ（t）的變化ΔΓi（t）＝Γi＋1（t）－Γi（t），則根據(jù)Kelven定律，從翼后緣脫離出的離散尾渦環(huán)量ΓWAKEi＝－ΔΓi（t）。為了確定該渦環(huán)位置，可假定在一定區(qū)域內(nèi)其渦強不變且不耗散。若用虛線表示第i時段內(nèi)翼后緣的運動軌跡，如圖2所示，則離散渦ΓWAKEi應位于該虛線上。需指出的是，確定渦環(huán)位置時，應避免將其渦心置于機翼后緣上，否則該渦的誘導速度將成為無窮大。本文將渦心定位于距翼后緣0.25V（t）Δt的位置。

圖1 尾渦離散示意圖Fig.1 Schematic of discretization of the wake vortex

當前時段前生成的離散渦均以當?shù)貧饬魉俣萔WAKEI運動，可通過對機翼及整個尾跡區(qū)在該點的誘導速度進行疊加求得。若設離散渦在第i時段相對于前一時段的位移ΔSi＝VWAKEiΔt，則其下一時刻渦環(huán)的位置為

依據(jù)上式可遞推出各離散尾渦的當前時刻位置，進而可得到各尾渦的運行狀態(tài)。

圖2 新生成尾渦的位置Fig.2 The position of the new wake vortex

1.4 FMAV 載荷計算

撲翼的周邊流場用非定常伯努力方程描述［6］

式中，pref為距翼面無窮遠處的壓強。與定常伯努力方程相比，上式的不同在于它增加了勢函數(shù)偏微項?Φ／?t。FMAV的承載能力隱含在壓強p之中并可通過數(shù)值解法求得。

1.5 UVLM數(shù)值求解

計算FMAV的氣動特性時，需先對它進行網(wǎng)格化，其方法與定常渦格法相似［8］，只是因FMAV的翼面形狀不斷變化，因此在各計算周期內(nèi)需對網(wǎng)格重新進行劃分。為此，定義翼弦向為i向，展向為j向，并將翼面劃分為M×N＝m個面元，且為各面元配置一個四邊形渦環(huán)。渦環(huán)前緣位于面元1／4弦線處，并取面元3／4弦線的中點為控制點，該點處氣流速度滿足式（1）。設面元兩對角線向量為A和B，則其單位法矢為

渦環(huán)環(huán)量方向由右手法則確定。這樣，利用m個渦環(huán)可模擬機翼對流場之作用。渦環(huán)（i，j）在任意位置產(chǎn)生的誘導速度為渦環(huán)四條邊共同誘導的結(jié)果，它可由B－S定律確定，即

dv為微量渦線段dl的誘導速度，r為渦線段至誘導點距離向量。通過對渦環(huán)上各渦線的誘導作用進行疊加，可確定渦環(huán)對指定點的誘導速度。面元的控制點處氣流速度是由翼面運動產(chǎn)生的相對氣流速度（u′，v′，w′）K、面元渦誘導速度、尾流誘導速度（uw，vw，ww）K構成。設第1面元單位強度渦環(huán)在第k面元控制點的誘導速度為（u，v，w），并引入氣動影響系數(shù)ak1表示其誘導作用，即

類似地，將其余面元的影響系數(shù)依次表示為aK2，…，akm。依據(jù)式（1），可將第 面元的無穿透邊界條件可完成

將上式中的已知項移至等號右側(cè)，可得

上式寫成矩陣方程形式，有

求解上式可得各面元渦量。再通過對式（4）進行差分處理，就可得到面元（i，j）的壓差

式中，τi、τj分別為面元弦向和展向切矢，Δcij和Δbij是面元弦向和展向?qū)挾?。通過疊加各面元壓差即可求得機翼的氣動力，即

F的垂直和水平分量即為撲翼升力和推力。

2 面向?qū)ο蟮腢LVM實現(xiàn)

對于給定翼面，可將其離散成一系列四邊形渦環(huán)，各渦環(huán)又可分解為四條渦線段。將此分解思路進行抽象，便得到UVLM面向?qū)ο蟮膶崿F(xiàn)構架。為此，可先構建渦線段類以保存渦線段有關信息及其誘導速度算法；再使用該類構建渦環(huán)類，渦環(huán)類包含四個渦線段子類，保存面元法向量及控制點等信息；基于渦環(huán)類構建渦環(huán)網(wǎng)格類，網(wǎng)絡類將渦環(huán)按照給定的幾何參數(shù)組織成網(wǎng)格以描述翼面和尾跡，利用該類可派生出翼面類和尾跡類。最后，使用上述基本類構建機翼運動、線性方程求解、氣動力計算和輸出結(jié)果數(shù)據(jù)的輔助類及控制時間步進迭代類。這樣就可編寫出UVLM程序，圖3給出了程序框架。需要指出的是，因尾跡渦環(huán)量的更新計算量巨大，而該運算又特別適于并行處理，為提高算法效率，本文引入GPU求解器并利用其并行計算能力實現(xiàn)尾渦求解［10］，從而使ULVM求解速度達到CPU求解時的3倍。這不僅提高求解速度，更重要的是它使UVLM可作為FMAV優(yōu)化迭代的氣動計算工具。

圖3 UVLM的面向?qū)ο蟪绦蚩蚣蹻ig.3 The object oriented program framework for UVLM

3 模式搜索法尋優(yōu)

3.1 優(yōu)化方法確定

FWAV的氣動優(yōu)化基于UVLM的海量氣動計算。由于UVLM是一種不依賴于解析目標函數(shù)的離散數(shù)值法且它不提供導數(shù)信息，因此，該方法不適于需對目標函數(shù)進行求導的優(yōu)化場合（如牛頓迭代法）。眾所周知，PS法是一種直接搜索法，無需解析性目標函數(shù)且尋優(yōu)效率較高，因此，本文的撲翼優(yōu)化即選定它作為FMAV的氣動尋優(yōu)工具。

3.2 PS法的實現(xiàn)

式中，B∈Rn×n為基矩陣，在各迭代步中保持不變；Ck∈Zn×p為生成矩陣，形如

其中，Mk∈Zn×n和Lk∈Zn×（p－2n）中至少包含一列零向量。搜索方向即取Pk的某一列。

對于給定搜索步長Δk＞0，定義探測向量

式中，為Ck的第i列。第k步迭代時，需從，…中尋找滿足如下條件的，即

若 min｛f（xk＋y），y∈Δk［BMk－BMkBLk］｝＜f（xk），則f（xk＋）＜f（xk），可取xk＋1＝xk＋。

在2n個搜索方向中只要有一個能使f（x）下降，則搜索步長與該列相乘得到的就將使目標函數(shù)值下降。依據(jù)此思路且通過多次優(yōu)化迭代，就可尋得上述問題的最優(yōu)解。根據(jù)上述理論，編寫了模式搜索法程序模塊，并將UVLM嵌入到該優(yōu)化程序，從而實現(xiàn)了FMAV氣動計算與尋優(yōu)的整合。經(jīng)整合后的程序具備了FMAV氣動優(yōu)化能力。

4 氣動特性優(yōu)化

4.1 撲動模型

本文選定柔性撲翼作為優(yōu)化對象。根據(jù)FMAV運動特點，翼面撲動時其柔性變形主要發(fā)生在機翼的弦向，因此優(yōu)化所針對的撲動模型只考慮翼的弦向的柔性，并假定其為

式中，ηα是用于描述柔性形變的弦向扭轉(zhuǎn)角，它隨翼面形變沿展向和弦向增大，形如

式中，ηmax為最大扭轉(zhuǎn)角，φη為扭轉(zhuǎn)形變翼面撲動間的相位差，x和y分別為翼面坐標，c、b分別是翼的弦長和展長。

4.2 撲動相角優(yōu)化

為驗證上述優(yōu)化方法的有效性，首先研究了FMAV的一維氣動優(yōu)化。眾所周知，鳥飛時，翅翼不僅要上撲還需下?lián)?，并且撲動時還存在繞前緣的俯仰角變化，即其撲翅和俯仰間存在相位差Δφ＝φα－φβ。根據(jù)鳥撲動產(chǎn)生時氣動力的機理，易知該相位差對于鳥的氣動特別是其升力有很大影響。這就給我們一個啟示，要改善柔性撲翼的氣動條件，必須對其撲動與俯仰間的同步關系進行合理規(guī)劃。基于這種認識，選取某矩形直機翼進行優(yōu)化，該翼的展弦比為8，攻角為5°，最大撲動角αmax＝5°，βmax＝15°，減縮頻率k＝0.1。優(yōu)化時，以平均升力系數(shù)的最大化為優(yōu)化目標，以撲動－俯仰相位差為優(yōu)化變量，并沿翼面展向和弦向劃分16×6個面元，不考慮弦向柔性，只考慮撲動角和機翼繞前緣俯仰角變化并假定俯仰角沿翼展方向逐漸增大，優(yōu)化模型為

用PS法求式（18）時，初始搜索點取為零點，初始步長設定為0.5，圖4給出了優(yōu)化搜索軌跡，PS法經(jīng)14次迭代后收斂于Δφ＝－1.227rad，即當撲動－俯仰相位差Δφ＝289.76°時，撲翼將獲得最大平均升力系數(shù)。該結(jié)果與文獻［11］中以15°為步長逐漸增大Δφ算得的平均升力變化曲線非常吻合，表明要使FMAV獲得盡可能大的升力系數(shù)，其撲動與俯仰運動的相位差有一個最佳值。它同時也表明，內(nèi)嵌了通過GPU實現(xiàn)UVLM的模式搜索法能有效地實現(xiàn)FMAV優(yōu)化。

圖4 模式搜索迭代軌跡Fig.4 Optimization of the phase difference between the flapping and pitching using Pattern－Search Method

4.3 翼型結(jié)構優(yōu)化

根據(jù)飛行器多學科優(yōu)化中的外形參數(shù)化設計思想，本文在矩直翼基礎上，進行了機翼氣動布局優(yōu)化計算。給定矩直翼參數(shù)為：翼展b＝0.6m，翼根與翼稍弦長均為c＝0.1m，翼面積S＝0.06m2。優(yōu)化條件是在保持翼面積不變前提下，將翼根和翼梢弦長減小，而在某一展向相對位置p∈［0，1］處將機翼的弦長增大為cmax∈［1.0，1.5］，如圖5所示。為了保證機翼面積不變，cmax與翼根弦長croot滿足關系

圖5 參數(shù)化外形優(yōu)化結(jié)果Fig.5 Optimal flapping wing shapes

優(yōu)化建模時，以式（17）為撲動模型，取αmax＝10°，βmax＝30°，減縮頻率ζ＝0.2，不考慮柔性形變引起的扭轉(zhuǎn)角，并分別以平均推力系數(shù)和平均升力系數(shù)為目標函數(shù)，以p和cmax作為設計變量。這樣，便得到撲翼的優(yōu)化模型

求解式（20）時，需先采用懲罰函數(shù)法將其轉(zhuǎn)為無約束優(yōu)化問題，再調(diào)用PS模塊進行優(yōu)化迭代，通過在迭代中進行大量UVLM計算尋找最佳氣動參數(shù)。圖6給出了迭代過程及其氣動狀態(tài)。表1給出了優(yōu)化結(jié)果。據(jù)此得到圖5所示兩種優(yōu)化翼形。可見，要提高平均推力，翼面宜呈倒梯形布局，且應使外翼段具有較大面積；反之，要想升力最大化，應使內(nèi)翼段面積取較大值。造成兩種翼形差異較大的原因，主要是因為翼在撲動過程中，其內(nèi)翼段與外翼段會形成幾何扭轉(zhuǎn)，而且當下?lián)鋾r，外翼段的扭轉(zhuǎn)使翼產(chǎn)生的合力在向前方向上有一個分量，該分量將形成推力和升力，但內(nèi)翼段卻主要產(chǎn)生升力，當翼面呈倒梯形時，其翼尖處的運動速度較快，造成翼剛度小、柔性形變大，因而可形成較大推力。但若翼采用正梯形布局，則內(nèi)翼段面積會增大，故可獲得較大升力。上述結(jié)果也意味著對于給定翼型，可能無法使其升力和推力同時最大化，它表明在翼面設計過程中需綜合多種因素進行權衡。

表1 氣動力優(yōu)化結(jié)果Table 1 Results of aerodynamic optimiation

4.4 撲動參數(shù)優(yōu)化

撲動參數(shù)優(yōu)化并非什么新問題，但此方面已有的研究主要針對單個參數(shù)的優(yōu)化。實際上，如同鳥飛時需同時協(xié)調(diào)好多個撲動參數(shù)一樣，要想FMAV產(chǎn)生理想的氣動狀態(tài)，同樣需考慮多個撲動參數(shù)的共同影響。為此，本文選取NACA0012矩直翼（展長為0.8，展弦比A＝8）為對象，以平均推力系數(shù)最大化為優(yōu)化目標，同時對撲動頻率與柔性扭轉(zhuǎn)角進行了優(yōu)化。根據(jù)撲翼飛行的實情，假定撲動頻率f∈［4Hz，10Hz］，弦向扭轉(zhuǎn)角φη∈［10°，45°］，來流速度Vo＝5m／s，撲動角βmax＝40°，俯仰角αmax＝0°，將翼面劃分為16×4個面元。這樣，其優(yōu)化問題描述為

表2列出了撲翼的多參數(shù)優(yōu)化結(jié)果?？梢姡蛊骄屏ο禂?shù)最大的撲動頻率為頻率范圍內(nèi)的最大值，而柔性形變扭轉(zhuǎn)角為24.85°。根據(jù)文獻［11］，當撲動頻率固定為6Hz而不斷增大扭轉(zhuǎn)角時，將在扭轉(zhuǎn)角φη＝30°處使平均推力達到最大值。但由本文的優(yōu)化結(jié)果可知，隨著撲動頻率增加，對應于最大推力系數(shù)的扭轉(zhuǎn)角將發(fā)生一定的變化，并且隨著頻率增加，最大值扭轉(zhuǎn)角應相應地減小，即翼面柔性呈降低趨勢才能保證在該撲動頻率下的推力最大。

表2 撲動頻率與弦向扭轉(zhuǎn)角優(yōu)化結(jié)果Table 2 Optimiation of twist angle and the flapping frequency

5 結(jié) 論

（1）基于UVLM并通過GPU實現(xiàn)其氣動計算的模式搜索優(yōu)化法，運算速度快、計算精度高，可望成為FMAV氣動優(yōu)化的重要工具。

（2）為提高FMAV的升力，應使撲翼的撲動與俯仰運動間存在一定相位差，并且當相位差接近289.76°時，其升力達到峰值。

（3）為實現(xiàn)FMAV推力最大化，應在增大撲動頻率時適度地減小機翼的柔性扭轉(zhuǎn)角。

（4）要使FMAV產(chǎn)生較大氣動推力，其翼面宜呈倒梯形結(jié)構；反之，要獲得大的推力，翼面結(jié)構宜呈梯形布局。

［1］SMITH M J C，WILKIN P J，WILLIAMS M H.The advantages of an unsteady panel method in modeling the aerodynamic forces on a rigid flapping wing［J］.J.Exp.Biol.，1996，199：1073－1083.

［2］VEST M S，KATZ J.Unsteady aerodynamic model of flapping wings［J］.AIAAJournal，1996，34（5）：1435－1440.

［3］SUN M，TANG J.Aerodynamic force generation and power requirements in forward flight in a fruit fly with modeled wing motion［J］.J.Exp.Biol.，2003，206（12）：3065－3083.

［4］JOSEPH K，ALLEN P.Low speed aerodynamics：from wing theory to panel method［M］.New York：McGraw －Hill Book Co.1991.

［5］LIN S Y，HU J J.Aerodynamics performance study of flapping－wing flow fields［C］.23rd AIAA Applied Aerodynamics Conference 6－9June 2005：4123－1427.

［6］曾銳.仿鳥微型撲翼飛行器的氣動特性研究［D］.［博士學位論文］.南京：南京航空航天大學，2005.

［7］FRITZ，LONG L N.Object－oriented unsteady vortex lattice method for flapping flight［J］.JournalofAircraft，2004，41（6）：1275－1290.

［8］DANIEL T L，COMBES S A.Flexing wings and fins：bending by inertial or fluid－dynamic forces［J］.IntegrativeComparativeBiology，2002，42（10）：1044－1049.

［9］肖天航，昂海松.大幅運動復雜構形撲翼動態(tài)網(wǎng)格生成的一種新方法［J］.航空學報，2008，29（1）：41－48.

［10］LEWIS R M，TORCZON V.Pattern search method for bound constrained minimization［J］.JournalonOptimization，1999，28（4）：1365－1372.

［11］曾銳，昂海松.撲翼柔性及其對氣動特性的影響［J］.計算力學學報.2005，26（6）：750－754.