基于非定常渦格法的撲翼飛行器氣動(dòng)特性?xún)?yōu)化

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

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

0 引 言

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

1 非定常渦格法建模

1.1 翼面速度表示

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

1.2 撲翼流場(chǎng)邊界條件

機(jī)翼?yè)鋭?dòng)時(shí)，除前緣及翼尾的部分區(qū)域外，其它各處流場(chǎng)為不可壓無(wú)旋流，為此定義速度勢(shì)函數(shù)Φ（X，Y，Z）并滿(mǎn)足：▽2Φ＝0。因撲翼對(duì)遠(yuǎn)處流場(chǎng)影響很小，因此其遠(yuǎn)場(chǎng)邊界條件為Φ＝0。而近場(chǎng)邊界條件則滿(mǎn)足無(wú)穿透條件，即

式中，VS為翼面相對(duì)運(yùn)動(dòng)速度，V0為來(lái)流速度，VT、Ω分別為翼面平動(dòng)和轉(zhuǎn)動(dòng)速度，n為翼面單位法矢。

1.3 尾渦模型

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

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

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

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

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

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

1.4 FMAV 載荷計(jì)算

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

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

1.5 UVLM數(shù)值求解

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

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

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

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

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

上式寫(xiě)成矩陣方程形式，有

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

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

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

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

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

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

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

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

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

3.2 PS法的實(shí)現(xiàn)

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

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

對(duì)于給定搜索步長(zhǎng)Δk＞0，定義探測(cè)向量

式中，為Ck的第i列。第k步迭代時(shí)，需從，…中尋找滿(mǎn)足如下條件的，即

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

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

4 氣動(dòng)特性?xún)?yōu)化

4.1 撲動(dòng)模型

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5 結(jié) 論

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

（2）為提高FMAV的升力，應(yīng)使撲翼的撲動(dòng)與俯仰運(yùn)動(dòng)間存在一定相位差，并且當(dāng)相位差接近289.76°時(shí)，其升力達(dá)到峰值。

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

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

［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］曾銳.仿鳥(niǎo)微型撲翼飛行器的氣動(dòng)特性研究［D］.［博士學(xué)位論文］.南京：南京航空航天大學(xué)，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］肖天航，昂海松.大幅運(yùn)動(dòng)復(fù)雜構(gòu)形撲翼動(dòng)態(tài)網(wǎng)格生成的一種新方法［J］.航空學(xué)報(bào)，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］曾銳，昂海松.撲翼柔性及其對(duì)氣動(dòng)特性的影響［J］.計(jì)算力學(xué)學(xué)報(bào).2005，26（6）：750－754.