撲翼柔性變形對懸停氣動特性影響的數值研究

張興偉，周超英，謝 鵬

(哈爾濱工業(yè)大學 深圳研究生院，518055廣東 深圳，cyzhou@hit.edu.cn)

撲翼的氣動特性研究是設計、優(yōu)化和控制微型撲動飛行器(Micro Air Vehicle，MAV)的重要環(huán)節(jié).近年來撲翼空氣動力學特性的研究有了很大的發(fā)展，但對于MAV所處的低雷諾數下的非定常氣動特性、撲翼柔性變形對氣動特性的影響等一些突出問題還有待進一步研究.特別是隨著材料、生物電池等技術的進步和發(fā)展，撲動飛行器進一步微型化、柔性化的發(fā)展趨勢日趨明顯，這些新型MAV的尺度、飛行雷諾數等參數都與小型昆蟲的相似，因此，研究昆蟲撲翅獲取高升力和推力的飛行機制能為發(fā)展MAV提供一些參考和靈感，其中研究低雷諾數下柔性翅的氣動特性尤為重要.近年來國內外學者對昆蟲的撲翼飛行方式及其飛行機制進行了深入研究，取得了顯著成果[1-5].

昆蟲在飛行中通過撲動翅膀來調整升力和推力從而改變其飛行速度和高度.近年來人們通過對多種飛行昆蟲的實驗觀測發(fā)現昆蟲撲翼在撲動飛行中會產生不同程度的彎曲和扭轉，不同昆蟲這種彎曲和扭轉的方式和程度也不同，這些變形對昆蟲撲動飛行性能有顯著的影響.此外，在MAV樣機試驗中也發(fā)現，采用仿生材料研制的撲動機翼在拍動過程中產生的柔性變形對撲翼的氣動特性有著重要影響.早期的數值研究大都假設翅翼為剛性無變形，隨著計算機技術和流體力學計算方法的進步和發(fā)展，越來越多的數值研究關注柔性變形對撲翼氣動特性的影響.本文采用流體力學數值計算方法有限體積法求解Navier-Stokes方程，對柔性撲翼在懸停時的氣動特性進行了研究，分析了撲翼的柔性變形對其懸停飛行性能的影響.一般而言昆蟲的撲動飛行是一個三維問題，在考慮展向尺度足夠長以及無展向扭轉變形的前提下，可忽略展向流影響從而將翼型簡化成二維模型來研究其撲動的非定常效應，對此Wang等[6]驗證了二維模型數值計算結果與三維實驗結果的可比性.為簡化問題，在此選取的柔性撲翼模型也為二維模型.

1 可變形撲翼物理模型及運動模式

建立弦長為c，最大厚度為弦長2%，截面為橢圓的二維撲翼模型，翼的懸停飛行及其柔性變形如圖1所示，其中:坐標xoy為固聯于地面的慣性坐標系;x'ory'為固定在翼表面的移動坐標系;坐標原點or固定于距翼前緣c/4處.在懸停撲動過程中，x'ory'坐標系相對于慣性坐標系xoy做往復平動和轉動.

圖1 懸停中撲翼往復平動和轉動示意圖

翼的懸停撲動由往復平動和轉動組成.在xoy坐標系下，平動位移H和轉動轉角α描述為

式中:A為撲動幅值;f為撲動頻率;T為撲動周期(T=1/f);α0為轉角幅值.通過實驗和數值研究表明，當A/c=3～5時，撲翼前緣處產生的渦即前緣渦不會從翼表面完全分離，從而產生較大的升力[6－7].因此，在計算中 A 取值為 4.5c，α0取值為45°.α(t)與瞬時攻角有90°的相位差，即瞬時攻角 =90°－|α(t)|.撲動平面為水平面.

昆蟲翅膀主要由位于前緣的翅梁、翅脈和翅膜組成.翅膀前緣的剛度比較大，在飛行中形變非常小可近似成剛性;后緣剛度較小，懸停時會產生相對來說較大的變形.對所采用的橢圓截面翼型，考慮到翼的弦向變形尺度相對于弦長很小，因此，假設變形前后翼的弦長和厚度均保持不變;翼的前緣部分不發(fā)生變形，后部分的變形起點位于距前緣1/4弦長處，即x'ory'坐標中的原點or點.在懸停過程中，從or點到翼前緣這部分不發(fā)生變形，僅做平動和轉動;從or點到翼后緣這部分除了相對于坐標系xoy做平動和轉動外，同時還產生在x'ory'坐標系下相對于or點的柔性變形.圖2是翼型變形的示意圖，d/c=0.25，翼型的后緣點為最大變形處.取弦向的最大彎曲尺度為a0c，其中無量綱系數a0表征了后緣端點相對于弦長的變形程度.

圖2 懸停中翼變形的示意圖

本文主要研究翼型弦向的動態(tài)彎曲變形對懸停氣動特性的影響，因此計算中沒有考慮由于展向扭轉而引起的展向流影響.翼型在往復平動和轉動復合運動中產生的動態(tài)柔性形變可用Miao＆Ho[8]的模型進行描述，即在 x'ory'坐標系下，從or點到翼后緣部分翼型的弦向變形表述為

式中 φ 為相位角.Read[9]，Heathcote[10－11]等發(fā)現，柔性翼的最優(yōu)相位角為90°，因此，在此動態(tài)柔性模型的相位角φ取值為90°.b0為變形系數，當彎曲的起點or位于c/4處時，若翼型后緣點最大彎曲達0.1c，即取a0=0.1，則變形系數值為b0≈ 0.177 8.

利用坐標轉換關系式為

即可實現翼型變形變量從x'ory'坐標系到xoy坐標系的轉換.

2 計算方法

均勻不可壓縮流體運動的控制方程為Navier-Stokes方程

式中:ui為流速，二維情況下 i=1，2，j=1，2;P為動壓;雷諾數定義為Re=ρcuref/μ，其中ρ為流體密度，μ為流體粘性系數，uref為撲動平面內的平均平動速度.

采用有限體積法[12]，在控制體積中對上述控制方程進行積分，采用交錯網格[13]實施積分方程在有限單元上的數值離散，應用PISO壓力速度耦合法[14-15]迭代求解壓力場和速度場.時間項采用一階差分格式，對流項采用精度為二階的迎風差分格式，粘性項用二階中心差分格式.計算區(qū)域四周外邊界取在距翼型壁面30 c處，邊界條件給定為遠場無來流條件;翼型表面邊界條件為無滑移邊界條件.計算中選取了三角型網格.采用動網格技術實現了對撲動和變形動邊界的模擬，網格的運動變形滿足空間守恒定律[16].計算中首先對網格密度以及時間步長對計算結果精度的影響進行了大量試算和驗算，最終確定時間步長為Δt=0.001T，網格初始節(jié)點數為3.8 × 104.

翼型的升力系數CL和阻力系數CD分別定義為

式中:urms為平動速度的均方根值;FL和FD分別為瞬時升力和阻力.撲動周期內的平均升、阻力系數分別為

3 計算結果及分析

3.1 程序驗證

首先，通過計算剛性翼懸停時的氣動力及流場分布，檢驗了動網格計算的可靠性.圖3給出了Re=75時懸停剛性翼的氣動力計算結果，其翼型運動參數參照了文獻[6]的模型參數.結果顯示，所采用動網格方法的計算結果與文獻[6]的結果符合較好.

3.2 弦向彎曲變化對氣動力的影響

計算中假定了翼型不會因為流場變化而產生自適應的變形.為了研究柔性翼在懸停時不同變形尺度對氣動力產生的影響，分別計算了6種不同變形幅度翼型在Re=140時的氣動力及流場分布，其中 a0為 0.00，0.05，0.10，0.15，0.20 和0.30.

圖4給出了a0=0.05時柔性翼在懸停運動狀態(tài)下升、阻力系數的時程圖，并與相同雷諾數下剛性翼的結果進行了比較.結果顯示，柔性翼升、阻力系數的瞬時峰值分別由剛性翼的3.5和4.6提高到4.8和5.7，顯然懸停飛行中柔性彎曲對瞬時氣動力峰值的影響非常顯著.另外，柔性彎曲變形沒有改變瞬時氣動力的分布，即在一個撲動周期內，剛性翼和柔性翼升、阻力峰值發(fā)生的時刻基本相同，顯示旋轉環(huán)流、尾跡捕捉等非定常效應與翼型的變形機制相比仍然是決定非定常氣動力分布的主要因素.

圖3 剛性翼懸停時升、阻力系數時程曲線(Re=75)

為研究柔性變形對氣動特性的影響，對升、阻力系數的周期時均值及升阻比進行了比較分析.圖5 為Re=140 時變形尺度系數a0為0.00，0.05，0.10，0.15，0.20和0.30時升、阻力系數時均值隨a0的變化曲線.結果表明，當后緣產生一定程度彎曲時如a0=0.05時，雖然升、阻力系數時均值的變化幅度同瞬時值的變化幅度相比較小，但柔性變形還是使翼翅的時均氣動性能有了一定的改變，升、阻力系數時均值略有增大，因此適當的柔性彎曲可以改善翼的懸停特性，這一結果與金曉怡等[17]試驗觀察的現象一致;但是，隨著后緣點彎曲幅度a0繼續(xù)增大，升、阻力系數時均值快速減小，Zhao等[18]通過樣機試驗發(fā)現了相同的現象，他們還發(fā)現當翼翅保持某一固定彎曲幅值時，阻力系數時均值隨功角轉動幅值增大而增大.圖6為升阻比隨變形尺度a0的變化規(guī)律，結果顯示升阻比變化趨勢與升、阻力系數時均值的變化趨勢基本相同，柔性翼在a0=0.05，0.10，0.15 時得到了高過剛性翼的升阻比，但隨著彎曲度的進一步增大，升阻比迅速減小，氣動性能惡化.

圖4 撲翼懸停時升、阻力系數時程曲線(Re=140)

圖5 柔性翼時均升、阻力系數隨變形幅度的變化曲線

圖6 升阻比隨變形幅度的變化曲線

3.3 弦向彎曲程度變化對渦量分布影響

圖7和圖8分別給出了剛性翼和不同變形幅度柔性翼在1/2撲動周期內的幾個典型時刻的渦量分布圖.當翼型位于最左端位置時(圖7和圖8中t/T=0.5)，上半個周期產生的前緣渦位于翼型右半面的前緣附近.隨著撲翼向右運動，前緣渦向后緣發(fā)展，最后脫落產生較長的尾渦;同時，翼左半面的上緣處逐漸產生新的前緣渦(圖7和圖8 中 t/T=0.625).

圖7 剛性翼在不同時刻的渦量分布圖

對比圖7，圖8可知，翼的柔性變形使得撲動過程中產生的前緣渦渦強變弱(圖7中t/T=0.875和圖8(a)、圖8(b)中 t/T=0.875).因此，相應的升、阻力時均值與剛性翼的相比較小.這一現象與文獻[18]通過樣機實驗測量的研究結果相一致.柔性翼流場與剛性翼流場相比，后緣處的渦強變弱，且后緣處在半個周期內產生的尾跡在x方 向變長(圖8(a)、(b)中t/T=1.0)，顯示柔性彎曲不但改變了前緣、后緣渦的強度也改變了尾跡脫落渦的傳輸方向.

圖8 柔性翼在不同時刻的渦量分布圖

4 結論

1)撲翼的弦向彎曲變形能改變瞬時氣動力的峰值.以后緣最大變形尺度為0.05倍弦長的柔性撲翼為例，升力和阻力最大峰值的增加分別可達37%和24%.

2)撲翼弦向彎曲變形對氣動力時均值的影響相對瞬時值較小.適度的彎曲變形可以有效增大升力從而獲得較高的升阻比，但過大的柔性彎曲會導致翼型的氣動特性急劇惡化，升、阻力系數時均值及升阻比迅速減小.

3)撲翼撲動產生的非定常效應同翼的彎曲變形機制相比仍是決定氣動力分布的主要因素.但柔性變形能在一定程度上改變前緣和后緣渦的渦強，以及脫落渦的傳輸方向.

[1]SANE S P.The aerodynamics of insect flight[J].Journal of Experimental Biology，2003，206(23):4191-4208.

[2]WANG Z J.Dissecting insect flight[J].Annual Review of Fluid Mechanics，2004，37(1):183-210.

[3]鮑麟，童秉綱.模型翼拍動中動態(tài)柔性變形效應的數值研究[J].中國科學院研究生院學報，2005，22(6):676-684.

[4]WALKER S M，THOMAS A L R，TAYLOR G K.Deformable wing kinematics in free-flying hoverflies[J].Journal of The Royal Society Interface，2010，7(42):131-143.

[5]周超英，林玉峰.撲翼氣動力特性的數值研究[J].哈爾濱工業(yè)大學學報，2006，38(9):1403-1405.

[6]WANG Z J，BIRCH J M，DICKINSON M H.Unsteady forces and flows in low Reynolds number hovering flight:Two-dimensional computations vs robotic wing experiments[J].Journal of Experimental Biology，2004，207(3):449-460.

[7]DICKINSON M H，G?TZ K G.Unsteady aerodynamic performance of model wings at low Reynolds numbers[J].Journal of Experimental Biology，1993，174:45-64.

[8]MIAO J M，HO M H.Effect of flexure on aerodynamic propulsive efficiency of flapping flexible airfoil[J].Journal of Fluids and Structures，2006，33(3):401-419.

[9]READ D A，HOVER F S，TRIANTAFYLLOU M S.Forces on oscillating foils for propulsion and maneuvering[J].Journal of Fluids and Structures，2003，17(1):163-183.

[10]HEATHCOTE S，MARTIN D，GURSUL I.Flexible flapping airfoil propulsion at zero freestream velocity[J].AIAA Journal，2004，42(11):2196-2204.

[11]HEATHCOTE S，GURSUL I.Flexible flapping airfoil propulsion at low Reynolds numbers[J].AIAA Journal，2007，45(5):1066-1079.

[12]VERSTEEG H K，MALALASEKERA W.An Introduction to Computational Fluid Dynamics:The Finite Volume Method[M].New York:Wiley，1995.

[13]王福軍.計算流體動力學分析:CFD軟件原理與應用[M].北京:清華大學出版社，2004.

[14]ISSA R I.Solution of the implicitly discretized fluid flow equations by operator-splitting[J].Journal of Computational Physics，1985，62(1):40-65.

[15]ISSA R I，AHMADI-BEFRUI B，BESHAY K R，et al.Solution of the implicitly discretized reacting flow equations by operator-splitting[J].Journal of Computational Physics，1991，93:388-410.

[16]LESOINNE M，FARHAT C.Geometric conservation laws for flow problems with moving boundaries and deformable meshes and their impact on aeroelastic computations[J].Computer Methods in Applied Mechanics and Engineering，1996，134(1/2):71-90.

[17]金曉怡，顏景平，周建華.仿生撲翼飛行機器人柔性翅試驗研究及其理論解釋[J].中國機械工程，2007，18(9):1028-1031.

[18]ZHAO L，HUANG Q，DENG X，et al.Aerodynamic effects of flexibility in flapping wings[J].Journal of the Royal Society Interface，2010，7(44):485-497.