仿鳥撲翼飛行器氣動力學建模精度測試

鄭和超，王建輝，胡紫陽，張忠海，何廣平，

1.北方工業大學 機械與材料工程學院，北京 100144

2.北京工業大學 材料與制造學部，北京 100124

3.北京航天測控技術有限公司，北京 100041

仿生學和空氣動力學研究結果表明，撲翼飛行器在利用渦流能量方面的效果突出［1］，在降低氣動噪聲、仿生隱藏、抵近偵察等方面優于固定翼和旋翼飛行器，近年來得到國內外研究人員的廣泛關注［2-3］。然而由于缺乏系統設計理論，目前試驗樣機［4-6］的設計制造還主要依靠形態仿生、直覺經驗、反復試錯和試驗修正。關于撲翼飛行器的研究主要包括：驅動機構設計［7］與運動學優化［8］、系統控制［9］及機翼形狀和結構參數對氣動力特性的影響［4，10-19］等。其中，對驅動機構和機翼結構特性測試和優化是重要的設計基礎。研究人員對此也進行了多種嘗試，例如，機翼展弦比和偏置角［10］、展向和弦向弱剛度［4，11-16］對撲翼空氣動力學的影響、撲翼氣動中心的測量方法［20］、撲翼周圍氣流煙霧的可視化［5］。這些工作主要集中在尋找機翼參數的設計原則或設計指標，相關經驗和結論可用于仿鳥撲翼飛行器的設計。

除了試驗，通過數值模擬研究撲翼的氣動性能也得到了較多關注。例如，Tong等［21］利用渦面場闡明了低展弦比撲翼繞流的渦動力學，并建立了用于估算尾跡中脫落渦面推力模型，但該模型僅適用于尾跡中具有離散渦面的撲翼。Yang等［1］通過數值模擬和試驗研究了仿鴿撲翼飛行器的受力機理，利用數字圖像動態變形數據分析氣動彈性效應和計算慣性力，最終的數值仿真結果與試驗結果一致性較好。Chen等［22］利用試驗與準穩態和數值模擬研究了仿昆蟲微型撲翼飛行器兩側機翼相位差對氣動性能的影響，數值模擬結果與試驗測得的力和渦量場吻合程度較高。雖然通過數值模擬方法定量預報撲翼產生的氣動力，可考慮機翼周圍流體狀態及機翼表面和空氣介質之間交互耦合的作用，但是數值模型的計算往往十分耗時，模擬結果也很難用于對撲翼飛行器的性能優化、閉環控制器設計和系統的運動穩定性分析。

針對以上問題，一些研究人員提出利用解析法預測撲翼飛行器的氣動力特性。Lane等［4］采用改進的葉素動量理論方法對撲翼飛行器推力預測，但由于該建模方法未考慮機翼的剛度和材料特性，并沒有正確表征機翼的運動學和幾何形狀，預測的平均推力僅在較高拍動頻率下與試驗結果基本一致。Nguyen等［10］基于齒輪齒條機構的實測推力，提出了一種預測推力經驗公式。該公式考慮了機翼面積、撲動頻率、撲動角和機翼展弦比等設計參數，但并未考慮機體姿態角和前向飛行速度等動態飛行參數，因此不能預報撲翼飛行器在實際飛行過程中的氣動性能。基于新型測試平臺，Mueller等［12］提出了一種預測撲翼飛行器升力和推力的半經驗氣動模型。該模型考慮了阻力和機翼柔度的影響，但預測結果僅在最低拍動頻率下與實測的瞬態氣動力相關性較好。年鵬等［23］提出了一種撲動軸瞬時氣動載荷半經驗高精度建模方法，最終模型確認系數大于0.89，但需基于風洞試驗數據建立該模型，不利于撲翼飛行器的廣泛研究。Jiao等［24］將一維細長體理論推廣，應用于二維平板驅動器的設計，提出了一種通用撲翼周期平均推力的解析估算方法，但這種方法在預報負推力時存在較大計算偏差。

以往相關研究工作表明，在撲翼飛行器氣動力預測模型中，需考慮機翼的運動學和幾何形狀、剛度及其彈性變形，以提高氣動力預測結果的準確性。考慮預測模型對樣機優化的實際意義，氣動力預測模型需具有足夠高的計算效率，并考慮撲翼機實際的飛行狀態。針對以上問題，提出了一種預測氣動力參數化解析計算模型。基于目前實驗室已有設備（高速攝像機和三維力傳感器），搭建了靜態和動態氣動力測試平臺，測得試驗樣機在試飛前的氣動性能，并明確氣動力模型中的機翼幾何形狀、運動學、變形角、俯仰角和前向飛行速度。與風洞試驗相比，該方法可以更低的成本，廣泛研究撲翼機的氣動性能。為進一步提高氣動力預測模型的準確性，基于靜態氣動力測試數據，修正了傳統偽穩態空氣動力學的力系數。并將該模型用于預測小型仿鳥撲翼飛行器產生的動態升力，取得了較好結果。

1 彈性平板翼仿鳥撲翼飛行器的氣動力學建模

1.1 撲翼機構的運動學模型

撲翼機構運動學如圖1所示，“曲柄-滑塊-搖桿”機構將曲柄旋轉運動轉化為搖桿往復擺動，左右搖桿通過齒輪嚙合傳遞動力，與以往的研究［1，4，6，8-13，15，22-24］相比，在結構上更加緊湊，兩機翼運動輸出具有更好的對稱性。并且搖桿輸出的拍動近似為余弦運動［16］。當機翼前緣沒有展向變形時，相應的拍動角φ0、角速度和角加速度可表示為

圖1 撲翼機構運動學示意圖Fig.1 Kinematic diagram of flapping-wing mechanism

本文研究的小型仿鳥撲翼飛行器拍動頻率較高、平板翼前緣剛度較小。實際拍動角φ除包含轉角φ0外，還需包含由機翼慣性力所導致的展向被動彎曲角ψs［11］。采用靜力學分析方法得到彈性平板翼的慣性力FI和展向被動彎曲變形角ψs的關系為

式中：ls為機翼氣動中心［9，20］距翼根的距離（見圖2）；ks為機翼前緣彎曲剛度；mw為機翼質量；L為機翼前緣長度。

圖2 彈性平板翼Fig.2 Elastic flat wing

當拍動振幅A=22°、偏置角ξ=12°、機翼質量mw=0.5 g、機翼前緣長度L=130 mm、平均拍動頻率=15.834 7 Hz、機翼前緣彎曲剛度ks=0.3 N·m/rad時，彈性平板翼的運動學規律如圖3所示。當考慮展向被動彎曲角ψs后，機翼實際拍動角φ的范圍接近51°。

圖3 彈性平板翼運動學（T =）Fig.3 Kinematics of elastic flat wing （T = ）

撲翼飛行器的空氣動力是由撲翼與流體之間相對運動產生，要建立氣動力計算模型，需分析機翼上任意一點的絕對速度［9］。彈性平板翼撲翼飛行器的坐標系如圖4所示。∑I代表固定在地球上的慣性系；∑B是固定在撲翼飛行器質心上的機體坐標系；∑W表示附著于機翼根關節的動坐標系；PIC表示機翼質心在慣性系∑I中的位置矢量；PIB是機體質心在慣性系∑I中的位置矢量；PWC是機翼質心在機翼坐標系∑W中的位置矢量；軸xB指向飛行器的前方，軸yB指向機體的右側，軸zB由左手定則確定；飛行器機體姿態定義為歐拉角α-β-γ；機翼的實際拍動角為φ；彈性平板翼在弦向方向上的扭曲變形角為ψc。為簡化坐標系變換，令軸xB和xW始終平行，而當機翼拍動角φ等于0時，軸yB和yW平行。

圖4 撲翼飛行器的坐標系Fig.4 Coordinate frames of flapping-wing vehicles

從∑W到∑B，以及從∑B到∑I的旋轉變換矩陣分別可寫為

機翼質心位置可以表示為

由于向量在機翼坐標系∑W中是不變的，即=0。線速度可表示為

1.2 氣動力學建模

在偽穩態空氣動力學范疇［9，25］，由流固交互作用產生的氣流慣性力是仿鳥撲翼飛行器的主要動力來源，而由流固速度耦合引起的力（如離心力和哥氏力引起的摩擦力）則是次要氣動力。在三維空間中，彈性平板翼氣動力計算公式［9］為

式中：r為機翼單位長度，dr為沿翼展方向的位置，如圖2所示；為機翼在機體坐標系∑B中的絕對線速度矢量(r)和(r)分別為空氣動力在機翼坐標系∑W中的法向和切向分力，其方向總是與相反；W(r)為機翼在位置r處的弦寬；ρair=1.205 kg/m3為空氣密度；?為機器人的飛行迎角，定義為機翼平面與絕對線速度之間的夾角；(?)和(?)為機體坐標系∑B中的空氣動力學系數。

對于實際的飛鳥，可通過控制羽毛姿態以調節翅膀的空氣阻力系數［25］，因此即使對于恒定的迎角?，空氣動力學系數(?)和(?)通常是時變的。但是，類似昆蟲的機翼結構，參考圖5，對于彈性平板翼，由穩態經驗得到相應的空氣動力學系數［9］為

圖5 三維空氣動力學分析圖Fig.5 Diagram for aerodynamics analysis in 3D

式中：φ為機翼的實際拍動角；ψc為機翼在弦向方向上的扭曲變形角。

而式（7）計算出的切向力dfWT(r)為機翼坐標系中表示的實際力矢量。通過式（5）和式（7），機翼的總法向力和切向力沿翼展方向積分表示為

式中：AW為機翼面積為撲翼氣動中心在機體坐標系中的絕對線速度。而撲翼在機翼坐標系和慣性系中的總空氣動力分別表示為

式中：FB為撲翼飛行器在機體坐標系下的總氣動力為從機翼坐標系到機體坐標系的旋轉矩陣為從機體坐標系到慣性系的旋轉矩陣，由式（3）給出。

根據以往的研究工作［11］，發現在較高拍動頻率下，機翼的弦向扭曲變形角主要由柔性機翼上的空氣阻力引起。采用靜力學分析方法得到彈性平板翼上法向空氣動力FzW和被動扭轉角ψc之間的關系。相應公式可寫為

式中：lc為機翼氣動中心距機翼前緣的距離（如圖2所示）；kc為機翼的扭轉剛度；W為機翼的弦向最大寬度。

根據以往工作測試和計算方法［20］，機翼的扭轉剛度kc=0.007 8 N·m/rad。當W=95 mm、機體的俯仰角β=0°、平均前向飛行速度=0 m/s時，利用本節給出的建模方法，小型仿鳥撲翼飛行器兩側機翼在慣性系中所產生空氣動力的解析計算結果如圖6 所示。兩側機翼在一個拍動周期的平均升力≈-0.000 1 N，平均推力≈0.077 N。由于撲翼飛行器的對稱設計，其側向理論空氣動力為0。

圖6 小型仿鳥撲翼飛行器的運動學和氣動力（T =）Fig.6 Kinematics and aerodynamic forces of small bird-scale flapping-wing vehicles（T =）

關于小型仿鳥撲翼飛行器的氣動力預測和分析，本文僅對以下主要參數開展研究：空氣動力學系數（(?)和(?)）、機翼展向的被動彎曲角ψs和弦向的被動扭轉角ψc、機翼的平均拍動頻率、機體的俯仰角β及平均前向飛行速度

2 撲翼氣動力測試系統

2.1 撲翼飛行器及其彈性平板翼

小型仿鳥撲翼飛行器試驗樣機如圖7所示，相應的結構參數如表1所示。為使偽穩態空氣動力學模型預測的氣動力更具普適性，設計了3 種剛度不同的機翼，如表2 所示。在設計制造試驗機翼過程中，為避免不同機翼形狀的差異性和安裝精度對氣動力測試結果造成影響，通過在同一個翼面上增加或者減少支撐碳桿，來保證單純改變機翼結構的剛度［15］。另外，選用具有較高強度的10 μm厚聚亞酰胺薄膜為翼面，采用直徑0.6 mm碳纖維桿作為前緣和骨架。并且碳桿、翼膜和機翼前緣之間都使用20 μm厚透明膠帶高溫高壓粘結。最終，單個試驗機翼質量被控制在0.5 g以內，以減小翅膀慣性力對氣動力測試結果的影響。

表1 小型仿鳥撲翼飛行器的結構參數Table 1 Structural parameters of small bird-scale flapping-wing vehicles

表2 機翼結構示意圖Table 2 Schematic diagram of wing structure

圖7 小型仿鳥撲翼飛行器試驗樣機Fig.7 Experimental prototype of small bird-scale flapping-wing vehicles

2.2 測試系統

靜態和動態氣動力測試平臺如圖8和圖9所示。高速攝像機（i-SPEED 3）捕捉的速率500 幀/s，三軸力傳感器的性能指標如表3所示。通過高速攝像機1和高速攝像機2在線捕捉機翼的瞬態變形角ψs和ψc。以機翼氣動中心為特征點，通過i-SPEED Suite便可獲得機翼前緣和后緣變形角的幅值和變化規律，用于偽穩態空氣動力學模型預測氣動力。為減小環境噪聲對氣動力測試結果的干擾，采用巴特沃斯低通濾波器對采集的力數據進行濾波。

表3 三維力傳感器的性能指標Table 3 Performance index of 3D force sensor

圖8 靜態氣動力測試平臺Fig.8 Test platform of static aerodynamic forces

圖9 動態氣動力測試平臺Fig.9 Test platform of dynamic aerodynamic forces

值得注意的是，靜態氣動力的測試分析可以考慮撲翼機零來流速度起飛時的實際情況，相應的推力特性對提高撲翼機飛行靈活性十分重要。而在動態升力的測試分析中，飛行器主動繞旋轉機構的被動轉軸旋轉，通過調整俯仰角調節機構，可以考慮飛行器在不同俯仰角β下的動態升力。基于高速攝像機3獲得的旋轉飛行周期Tr，當旋轉懸臂長度為Lr時，通過式（13）可計算飛行器的平均前向飛行速度

3 撲翼氣動力測試與解析模型修正

3.1 基于傳統偽穩態空氣動力學模型的氣動力分析和預報

基于靜態氣動力測試平臺的數據，使用第1節的氣動力學模型定量預測小型仿鳥撲翼飛行器在不同平均拍動頻率下的氣動力。3種不同剛度機翼實測和預測氣動力的比較如圖10～圖12所示，可看出預測氣動力與實測氣動力變化的趨勢一致。但是在行程開始和結束階段，實測與預測的瞬態氣動力數值及其變化曲率存在一定差異。文獻［5，25］研究結果也得到了佐證，這是由撲翼尾跡捕獲氣動力效應引起。在反轉運動時，機翼與前一行程中產生的流體尾跡相互作用。機翼后面的流體由于慣性，傾向于保持其速度。當機翼改變方向時，機翼和流體之間的相對速度大于絕對速度，從而在行程開始和結束階段會產生更大的力。而偽穩態空氣動力學系數只考慮了延遲失速和旋轉升力的非定常動力學效應［25］，因此在行程開始和結束階段所預測的氣動力數值和曲率存在誤差。除此之外，由于弦向剛度較小，機翼后緣會高階共振而降低后緣渦的強度，實測瞬態推力在行程中位面處的峰值不盡相同。這種現象在無骨架機翼的拍動行程中尤為明顯，而解析計算模型很難考慮撲翼后緣的異常波動現象［22］。通過定量分析圖10～圖12的氣動力數據，實測與預測的瞬態升力峰值的誤差均在-21.60%～0.87%，瞬態推力峰值的誤差均在-27.62%～-0.57%。

圖10 無骨架機翼的實測氣動力和預測氣動力比較（T =）Fig.10 Comparison of measured aerodynamic force and predicted aerodynamic force of frameless wings （T =）

圖11 少骨架機翼的實測氣動力和預測氣動力比較（T =）Fig.11 Comparison of measured aerodynamic force and predicted aerodynamic force of less-frame wings （T =）

圖12 多骨架機翼的實測氣動力和預測氣動力比較（T=）Fig.12 Comparison of measured aerodynamic force and predicted aerodynamic force of multi-frame wings （T=）

為直觀反映靜態氣動力與機翼剛度之間的

關系，定量預測試驗樣機在不同平均拍動頻率下的平均氣動力，結果如圖13所示。與實際測得的平均升力相比，可明顯看到預測的平均升力也幾乎為0 N［1，5，17］。這符合靜態測試中，飛行器的俯仰角為0°、平均前向飛行速度為0 m/s時，翅膀在上下沖程運動中產生升力相互抵消的實際情況。其次，與圖14所示變形角的變化呈現相反的規律，多骨架機翼的平均推力大于少骨架機翼，而無骨架機翼的數值最小。根據第1節的理論分析發現氣動力會受到展向和弦向變形的共同影響。多骨架機翼因為具有較小的展向變形和較穩定的弦向變形，使得機翼的結構變形能更穩定地向后緣方向偏轉，從而釋放出強大而且穩定的后緣渦，進而產生較大的推力［5，11，15，18-19］。隨著平均拍動頻率的增大，盡管預測的平均推力也隨之增大，但是相應的誤差卻在-10.42%～22.87%。因此，基于1.2節介紹的各向同性氣動力系數，只能定性預測撲翼機的氣動力，對開展樣機改進和優化，利用價值有限。

圖13 測量和預測的平均氣動力比較Fig.13 Comparison of measured and predicted data of average aerodynamic forces

圖14 3種機翼的變形角Fig.14 Deformation angles of three wings

3.2 偽穩態空氣動力學模型的修正和氣動力預報

為提高氣動力解析計算模型的預測精度，需要對模型修正。由以往的研究工作［12-13］可知，撲翼氣動力系數通常取決于相應的結構形狀和流體雷諾數。對于小型低速雙平板翼撲翼飛行器，各向同性氣動力學系數取值方式表現出了明顯的局限性。由于所研究機翼的迎角較大，由機翼表面摩擦引起的切向力分量只能起到很小的作用［25］，在保證運動學參數相同的情況下，基于圖13和圖14的測試數據，對式（14）所示的法向氣動力系數(?（)i=x，y，z）修正。而機翼表面的實際法向力主要來源于zW方向的法向氣動力分量，因此本文僅對相應的氣動力系數(?)修正，并假定其與機翼的弦向剛度kc有關。新的氣動力系數表現出各向異性特征，并呈現相同的變化規律，如式（15）所示，確定系數R2＞0.985 3，標準誤差RMSE＜0.030 9，可知建立的模型具有較高的準確性和參考性。

通過采用修正后的氣動力系數，3種不同剛度薄膜機翼實測與預測氣動力的比較結果如圖15～圖17所示。與傳統偽穩態空氣動力學模型相比，預測的氣動力峰值精度得到了改善：升力峰值的誤差從-26.60%～0.87%降到了-9.83%～9.07%；推力峰值的誤差從-27.62%～-0.57%降到了-9.05%～5.71%。由圖18可以發現，預測的平均升力近似為0 N。相較于傳統偽穩態空氣動力學模型的預測結果，理論預測的平均推力誤差從-10.42%～22.87%降到了-3.37%～7.36%，顯示出良好的預測精度。因此，在采用修正后的氣動力系數，解析計算模型的精度能夠得到有效提升，可以更加準確地預測撲翼產生的靜態氣動力。

圖15 無骨架機翼的實測氣動力和預測氣動力比較（修正）Fig.15 Comparison of measured aerodynamic force and predicted aerodynamic force of frameless wings （revised）

圖16 少骨架機翼的實測氣動力和預測氣動力比較（修正）Fig.16 Comparison of measured aerodynamic force and predicted aerodynamic force of less-frame wings （revised）

圖17 多骨架機翼的實測氣動力和預測氣動力比較（修正）Fig.17 Comparison of measured aerodynamic force and predicted aerodynamic force of multi-frame wings （revised）

圖18 測量和預測的平均氣動力比較（修正）Fig.18 Comparison of measured and predicted data of average aerodynamic forces （revised）

為了使修正后的氣動力系數更具普適性，如圖19所示，通過觀察3種不同剛度薄膜機翼的變形角和預測氣動力，可以發現與傳統偽穩態空氣動力學模型的預測規律相同，多骨架薄膜機翼因為具有較小的展向變形和較為穩定的弦向變形，使得相應的推力峰值和平均值最大，少骨架薄膜機翼次之，無骨架薄膜機翼的推力性能最差。并且3種機翼的升力峰值也表現出相同的規律。相應的平均升力因為在上下沖程中相互抵消，與展向變形和弦向變形的關系并不明顯。

圖19 預測的氣動力與展向變形和弦向變形的關系Fig.19 Relationship between predicted aerodynamic force and wingspan deformation and chord deformation

4 動態氣動升力的測試和預報

根據小型仿鳥撲翼飛行器靜態飛行氣動力測試和分析預測結果，可以發現雖然無骨架機翼和少骨架機翼的質量較小，但相應比剛度卻小于多骨架機翼，這并不符合偽穩態空氣動力學中的“運動剛體”對象假設，使得相應的氣動性能難以達到實際飛行需求。因此基于圖9測試平臺，測量了裝配多骨架機翼實驗樣機的升力，并用于驗證式（14）～式（15）所示氣動力系數的有效性。與傳統偽穩態空氣動力學模型相比，修正后的解析計算模型所預測的升力更接近測量值，結果如圖20～圖22所示。當俯仰角分別為5°、10°和15°時，升力峰值的誤差均有所減小，分別從-23.38%～-7.48%、-19.75%～-2.55%、-30%～5.38%降到了-13.98%～-8.56%、-13.71%～-0.63%、-13.65%～1.98%。

圖20 俯仰角為5°時多骨架機翼的實測升力和預測升力的比較Fig.20 Comparison of measured and predicted lift of multi-frame wings at pitch angle of 5°

圖21 俯仰角為10°時多骨架機翼的實測升力和預測升力的比較Fig.21 Comparison of measured and predicted lift of multi-frame wings at pitch angle of 10°

圖22 俯仰角為15°時多骨架機翼的實測升力和預測升力的比較Fig.22 Comparison of measured and predicted lift of multi-frame wings at pitch angle of 15°

圖23直觀反映動態升力與俯仰角和前向飛行速度之間的關系，定量預測了撲翼機的平均升力。隨著俯仰角的增大，前向飛行速度即使在一定范圍內減小，機翼產生的平均升力也會增大，說明俯仰角對升力有明顯的影響，這也與自然界鳥類的飛行規律相同。與實測平均升力的規律類似，隨著平均拍動頻率的增加，預測的平均升力也呈現近似線性增長的趨勢。但是相較于傳統偽穩態空氣動力學模型的預測結果，由式（14）～式（15）所預測的平均升力更接近實測平均值：理論預測值與測量值之間的誤差從-27.92%～-17.43%降到了-15.61%～-4.82%。因此，基于修正的氣動力系數提出的解析計算模型能夠較好地預測小型低速雙平板翼撲翼飛行器產生的氣動力，對開展仿鳥撲翼飛行器的優化設計和加快實驗樣機改進，有著重要的實際意義。

圖23 小型仿鳥撲翼飛行器的平均前向飛行速度和平均升力Fig.23 Average forward flight speed and average lift of small bird-scale flapping-wing vehicles

5 結論

1） 基于飛行器機體位置固定條件下的靜態撲翼氣動力測試結果，對傳統偽穩態空氣動力學模型的氣動力系數修正，表現出各向異性特征，并且修正后的公式具有相同的結構和變化規律。通過對2種模型預測結果進行對比，修正后的氣動力預測模型所預測的瞬態氣動力峰值和平均值更接近實測氣動力：靜態升力和靜態推力的峰值誤差由±22%和±28%均降到±10%，平均推力誤差由±23%降到±8%。

2） 基于多骨架機翼動態測試數據，對比了2種模型的動態升力預測精度。測試和分析結果表明，修正后的氣動力預測模型精度更高。升力峰值的誤差分別從±30%降到±14%；平均升力的誤差從±28%降到±16%。因此，改進后的解析計算模型能較好地預測小型低速雙平板翼撲翼飛行器產生的氣動力，對開展仿鳥撲翼飛行器優化設計和加快實驗樣機改進有著重要意義。