蜻蜓非均勻柔性前翼撲動飛行的氣動性能研究1)

劉笑塵 何心怡 何國毅 王 琦 孫書美

(南昌航空大學飛行器工程學院,南昌 330063)

引言

撲翼飛行器相比于固定翼飛行器具有隱蔽性、機動性和靈活性等特征,滿足現在現代化戰爭的需求,可以更好的運用于軍事偵查和精確打擊[1].在以美國為引領的多個國家[2],撲翼飛行器的研究在近來二十幾年當中掀起了一股浪潮[3].雖然有關于撲翼的研究文獻越來越多,但是大部分的工作集中于剛性撲翼在非定常條件下氣動規律的研究,在真實情況下,撲翼無論其尺寸大小都具有一定的柔性因此柔性翼結構變形引起的流場和氣動參數變化的研究具有重要意義.

在自然界中,一些動物通過翅膀、鰭、尾部、軀干等結構的擺動來獲得推力以及升力是非常普遍的現象.Liu等[4]發現適當的組合柔性撲翼頻率和彎曲度可以顯著改變推力的產生和壓力分布,提供更高的氣動效率.Chen等[5]在得出改進后的UAM 非定常氣動模型后,計算得出柔性翼相比于剛性翼來說減小了負升力.文獻[6]在流體和結構相互作用領域的發現表明,柔性機翼可以比剛性機翼產生更多的升力.Addo-Akoto等[7]通過對比柔性翼與剛性翼氣動參數,結果表明柔性翼產生的力表現出了明顯的相位延遲.以上均為對于撲翼固定柔性的研究,但是柔性對于撲翼氣動參數的影響,還包括展向柔性[8-11]和弦向柔性兩個重要的領域[12-13].張云飛等[14]對于展向柔性的研究表明沿展向引入柔性改變了撲翼有效攻角,適當的展向柔性可以減小阻力,過大的展向柔性增加阻力.文獻[15]研究了在低雷諾數條件下,撲翼的展向柔性和弦柔性對撲翼性能的影響,引入一定的展向柔性可以增加推進效率.被動的柔性變形在很多的實驗中被證明是可以增加運動效率[16-19].值得一提的是,飛行類動物的翅膀柔性不是均勻分布,而且這種不均勻分布反而提高了推進效率[20].Wood[21]發現當水翼的前1/3 部分的柔性比較小,后2/3 部分的柔性比較大時,推進效率是最高的,也就是不均勻的柔性分布的推進效率要優于均勻柔性分布.

隨著撲翼發展,當需要尺寸在15 cm 以下的飛行器進行作戰任務時,微小型撲翼飛行器(flapping micro aerial vehicle,FMAV)就起到了至關重要的作用.卓越的飛行本領使得蜻蜓成為了研究微小型撲翼飛行器最適合的仿生對象.Hamamoto等[22-23]基于有限元分析方法處理蜻蜓懸停時的流固耦合問題,表明柔性翅膀的變形對氣動特性具有較小的影響.鄭光鵬[24]對于蜻蜓翅二維截面進行流固耦合計算,得出翼型越薄變形越大,撲動頻率越大產生的推力和升力越大,適合于大前飛速度.徐明等[25]通過主動控制二維蜻蜓翼截面的變形,計算了二等蜻蜓翼型截面周期性的動態柔性變形引起的氣動力響應,發現在下拍過程中,動態柔性變形顯著提高升力和推力.郝淑文[26]利用Matlab 編程建立了適合求解柔性翼氣動效能的程序,發現撲動平面與水平面的夾角在一定范圍內增大,可以提高推力和升力系數.孟令兵等[27]基于計算流體力學/計算結構力學(CFD/CSD)的雙向流固耦合方法,通過C++編程計算了前緣具有較大剛性的蜻蜓后翼平板模型在同一前飛速度下的氣動參數,發現了柔性翼降低了飛行阻力提高推力.蜻蜓翼飛行性能的研究大多數集中在均勻柔性分布下的二維柔性截面對于氣動性能影響的論述,但是蜻蜓的飛行是一個三維問題.本文就弦向柔性對于蜻蜓前翼氣動性能的影響展開研究.在前人對于撲翼流固耦合研究的基礎上,根據三維蜻蜓翼模型來探究均勻柔性分布和變柔性情況下對于蜻蜓前翼氣動性能參數的影響,對比柔性翼在兩種柔性分布方式下相對于剛性翼的優劣勢.

1 撲翼模型與函數

1.1 幾何模型建立

通常蜻蜓的兩側翼距離較大,相互干擾的作用較小,文獻[28]研究了蜻蜓懸停時的身體同側的前后翼并沒有很大的相互干擾作用,故本文選取了一個蜻蜓右前翼,模擬單獨的蜻蜓翼揮拍過程中的柔性變形,忽略蜻蜓翼之間的相互干擾.

蜻蜓前翼具有非常復雜的微觀結構,由于文獻[29]通過數值模擬方法,表明凹凸不平的褶皺結構并不影響撲翼整體的氣動響應,故而本次計算過程中將蜻蜓前翼簡化成具有相同輪廓的柔性平板模型,根據蜻蜓前翼的尺寸參數[30],建立蜻蜓前翼的幾何模型如下圖1 所示,幾何參數包括:平均弦長C為8 mm,展長L為40 mm,厚度 σ 為0.18 mm 以及幾何面積S為284.1 mm2.

圖1 蜻蜓前翼平板模型Fig.1 Plat model of dragonfly forewing

1.2 撲動函數

確定撲動函數需要確定兩個因素:撲動幅值和撲動頻率.撲動幅值以及撲動頻率依照蜻蜓的飛行特點進行確定,在實際飛行過程中,蜻蜓撲動幅值在10°～ 60°以及撲動頻率在30 Hz～ 50 Hz[31].在具體設置蜻蜓前翼的撲動函數時,根據文獻采用20°撲動幅值和43 Hz 撲動頻率[30].給出撲動示意圖如圖2以及撲動函數為

圖2 前翼撲動示意圖Fig.2 Schematic diagram of front wing flutter

根據撲動角度函數先轉化為弧度制再進行求導,得到其撲動角速度函數

1.3 柔性分布函數

蜻蜓前翼從前緣到后緣的柔性分布模式有兩種.(1)連續均勻分布:蜻蜓前翅的柔性從前緣端到后緣端都是相同的;(2)連續變化分布:本文將蜻蜓前翼的柔性從前緣端到后緣端按函數規律變化,包括線性函數,平方函數.

在變柔性分布中定義蜻蜓前翼前緣楊氏模量17 GPa[27],定義后緣楊氏模量3.8 GPa[30].通過Position 場函數,調用前翼隨體坐標系中y坐標.因此前后緣的y坐標和楊氏模量構成兩對坐標點,得出(0.002,17)和(-0.006,3.8).線性函數表達式通過兩點即可得到,將二次函數頂點置于前緣處,并求解通過兩點的函數,即可得到平方函數的表達式.整個蜻蜓前翼根據柔性分布函數發生柔性變化:

均勻柔性分布

線性函數柔性分布

平方函數柔性分布

2 研究方法

2.1 計算方法

采用STAR-CCM+計算流體力學軟件,對蜻蜓前翼撲動模型進行雙向流固耦合分析.

剛性翼的撲動研究較較為簡單,主要采用動網格技術[30]實現對蜻蜓前翼撲動的仿真計算.柔性翼由于變形較大導致網格劇烈變化,因此在流固耦合的基礎上采用重疊網格和網格自適應技術.

將外流場和重疊區域均設置為流場域,流場域選擇隱式非穩態求解,湍流模型選擇Spalar-Allmaras湍流模型[32].在流場域中,設置重疊保護和單元網格質量校正,來提高撲動過程中的流域網格質量,重疊區域設置為浮動,可以跟隨前翼固體模型一起運動.

將前翼結構設置為固體區域,在固體域中采用隱式非定常固體應力求解模型,添加撲動函數和楊氏模量函數使蜻蜓前翼實現運動及柔性分布.

將固體區域中的前翼表面和重疊區域中的前翼表面設置為映射接觸面,用于數據傳遞.在求解器中添加流體結構耦合(雙向),再在固體運動的基礎上添加固體變形位移即可實現在撲動過程中疊加變形,變形后的結構重新計算氣動力,數據在此界面上一直相互傳遞,由此來實現雙向流固耦合.

2.2 控制方程

流體控制方程包含連續性方程和運動方程,其中連續性方程如下

運動方程

式中,ρ 為密度;u代表流體微團的速度矢量;F表示作用在單位質量上的質量力;P為應力張量.

固體控制方程一般由牛頓第二定律導出,如下

式中,ρ 是固體密度;σ 是柯西應力張量;F是體積力矢量;d是固體域當地加速度矢量;流固耦合需要流體域和固體域之間交換數據,該過程必須滿足一些物理量的相等和守恒,即滿足方程如下

式中,τs和 τf分別代表固體域與流體域力大小;ns和nf則接觸表面的法向單位向量;ds和df代表固體域與流體域微團位移矢量;第一項表示交界面上的流體域與固體域力相等,第二項表示位移相等.

3 計算前處理

3.1 流場區域網格劃分

建立一個長方體外流域,其尺寸大小為30C×30C×25C,如圖3(a)所示.將靠近蜻蜓前翼翼根的面的邊界條件設置為對稱面,將左面和上下面設置為速度入口邊界條件,上下邊界使用速度入口條件可以減少邊界對流場的影響[33],右面的邊界條件設置為壓力出口.

圖3 流體區域尺寸及網格Fig.3 The size of fluid domain and mesh

流域網格劃分采用多面體網格生成器和網格重構方法來提高網格質量.在流體域中增加體控制設置為加密區域,使得蜻蜓前翅在撲動過程中始終處于加密區域,采用多面體網格方法合理設置加密區域的尺寸大小和網格尺寸,保證了撲動時計算的精確性.

重疊區域采用多面體網格和網格重構方法劃分網格,區域大小如圖3(b)所示,翅翼的邊界層網格采用4 層棱柱層網格.加密區域尺寸及重疊區域在加密區域中的相對網格尺寸如圖3(c)所示.

3.2 固體域網格劃分

CSD 的數值方法之一是有限元法,將三維蜻蜓前翼的翼根為固定面,翼面為交界面,采用四面體網格和薄體網格的方法離散化蜻蜓前翼,即可得到固體求解時的有限元網格,固體域網格如圖4 所示.

圖4 蜻蜓前翼平板模型網格Fig.4 Mesh of forewing plat model

3.3 數值方法驗證

為了驗證本文方法的數值計算撲翼運動的正確性,選取了文獻[34]撲翼運動模型,并將計算結果同文獻結果對比.對比結果如圖5 所示,升力系數曲線隨時間的變化的趨勢極為近似,數值大小較為接近,因此該方法可以用于數值分析撲翼運動的氣動特性.

圖5 升力系數隨時間變化規律與文獻[34]對比結果Fig.5 The variation law of lift coefficient with time is compared to the results in Ref.[34]

3.4 網格無關性驗證

本文在來流速度為3 m/s、撲動頻率為43 Hz、撲動幅值為20°情況下,采用三套網格對蜻蜓前翼升力系數進行監測如圖6 所示,Case A 網格數為160 萬,Case B 網格數為230 萬,Case C 網格數為280 萬.可以看出來升力系數曲線圖基本一致,考慮到計算效率以及計算結果準確性,本文選用230 萬網格數進行計算.

圖6 網格無關性驗證Fig.6 Verification of independence of grid number

4 不同柔性分布下柔性翼結構變形結果與分析

4.1 不同柔性分布下柔性翼撲動過程

在均勻柔性分布條件下,蜻蜓前翼雖然都是在上撲結束和下撲結束時刻到達變形最大值,但是楊氏模量的變化導致前翼出現了不一樣的變形規律.例如在撲動過程中:(1)當楊氏模量較小時,如圖7(a)所示,翼尖在0T時刻達到最大下彎變形,在T/2時刻達到最大上翹變形,0T時刻雖然翼根部位處于最大上撲動角度狀態,但是翅翼整體處于最大下彎變形狀態;(2)當楊氏模量較大時,如圖7(b)所示,翼尖在0T時刻達到最大上翹變形且在T/2 時刻達到最大下彎變形,0T時刻翼根處于最大上撲動角度狀態,翅翼整體處于最大上翹變形狀態.

圖7 柔性翼撲動變形過程Fig.7 The deformation of flexible wing during flapping

在變柔性分布條件下,在撲動過程中蜻蜓前翼變形規律與楊氏模量較大時的翅翼變形規律保持一致,如圖7(b)所示.

4.2 不同柔性分布下柔性翼結構變形數值監測方法

采取在監測點監測變形量的方法,更容易準確得出柔性翅的翼尖變形量和扭轉變形量即扭轉角度.在翼尖處放置一個監測點F,用來監測翼尖在z 軸的變形量.分別在距離翼根20 mm (翼中部位)和36 mm 處(翼尖部位)的前,后緣放置4 個監測點分別命名為A,B,C,D監測前后緣在z軸的變形量,監測點位置如圖8 所示.

圖8 蜻蜓前翼的監測點位置Fig.8 Location of monitoring points on the forewing

4.3 不同柔性分布下柔性翼數值變形結果及流場分析

翼尖點F 的監測變形量隨時間的變化如圖9 所示,可以發現:(1)均勻柔性較大的前翼翼尖變形和均勻柔性較小、變柔性分布的翼尖變形的規律相反,和前文的撲動變形過程相對應;(2)變柔性分布時,二次函數分布時翼尖點最大變形量較小,這是因為二次函數柔性分布時在靠近前緣部位,楊氏模量的減小速度較慢,因此到達翼尖位置時楊氏模量大于一次函數柔性分布.

圖9 不同柔性分布下柔性翼的翼尖變形Fig.9 Tip deformation of flexible wing with different flexibility distributions

選取楊氏模量為3.8 GPa 均勻柔性分布和楊氏模量為10.4 GPa 均勻柔性分布的柔性翼在一個周期的初始時刻翼尖部位的壓力云圖,如圖10 和圖11所示.

圖10 楊氏模量為3.8 GPa 均勻柔性分布壓力云圖Fig.10 Pressure contours of uniform flexible distribution with Young's modulus of 3.8 GPa

圖11 楊氏模量為10.4 GPa 均勻柔性分布壓力云圖Fig.11 Pressure contours of uniform flexible distribution with Young's modulus of 10.4 GPa

從楊氏模量為3.8 GPa 均勻柔性分布的壓力云圖可以看出,上翼面壓強大于下翼面壓強是導致翼尖變形量在初始時刻為負位移的原因.

從楊氏模量為10.4 GPa 均勻柔性分布的壓力云圖可以看出,上翼面壓強小于下翼面壓強從而導致翼尖變形量在初始時刻為正位移變形.

可見柔性過大會導致翼尖會發生相反的柔性變形規律.

根據監測點A,B,C,D位移可以得到后緣減前緣即(B-A和D-C)的差值,示意圖如圖12 所示.當差值為正時,表示后緣點在前緣點上方,此時扭轉角度為負,當差值為負時,表示后緣點在前緣點下方,此時扭轉角度為正.

圖12 后緣點與前緣點差值示意圖Fig.12 Schematic diagram of difference between trailing edge point and leading edge point

再根據監測點之間的差值代入扭轉角公式可以得出蜻蜓前翼在翅翼中部和翼尖部位的扭轉角度.整個翅翼伴隨著撲動過程也一直在按照余弦函數規律做扭轉振蕩運動,差值最大值時刻對應的也就是扭轉角度曲線最大幅值時刻.

翅翼扭轉角度曲線的周期和撲動周期一致,因此確定幅值后也就得到了扭轉曲線,初始時刻的最大幅值如下表1 所示.其中β1代表的是在翼中部位20 mm 處的扭轉角,β2代表的是在翼尖部位36 mm處的扭轉角.扭轉角為正則此刻迎角為正,扭轉角為負則此刻迎角為負.

表1 蜻蜓前翼0T 時刻的扭轉角幅值Table 1 Torsional angle amplitude of dragonfly forewing at 0T

在均勻柔性分布時:(1)扭轉振蕩幅值隨著楊氏模量的增加逐漸減小,楊氏模量較小時在初始時刻的扭轉角度為正,與其他柔性分布相反,也是因為撲動過程的變形規律相反造成;(2)翼中和翼尖部位扭轉振蕩幅度不同,隨著楊氏模量的增加,又基本保持一致.

在變柔性分布時,翼尖部位的扭轉幅度大于翼中部位,說明翼尖部位后緣的上翹和下彎變形略微大于翼中部位的后緣.

5 不同柔性分布下蜻蜓前翼氣動特性分析

5.1 不同柔性分布下蜻蜓前翼的升力系數

從升力系數變化規律如圖13 所示,對比可以發現:(1)在均勻柔性分布時,柔性較大的翅翼升力系數曲線變化趨勢滯后于剛性翼半周期,這是因為從結構變形來看,將下撲運動離散之后,每個瞬時時刻剛性翼相對于來流做向下撲動運動,柔性較大的前翼相對于來流來說做向上變形運動,從而導致了半個周期的剛性翼向下運動,變為了向上變形運動;(2) 在變柔性分布時,柔性翼和剛性翼的變化趨勢相同.

圖13 不同柔性分布下柔性翼的升力系數Fig.13 Lift coefficients of flexible wings with different flexibility distributions

兩種柔性分布情況下,因為運動為單自由度的撲動運動,撲動過程上下對稱,所以只提高了升力系數峰值,并未改變正升力系數在一個周期內的時間占比,并沒有提高平均升力系數.

為了探究柔性翼如何提高升力系數峰值,本文選取線性函數柔性分布的柔性翼和剛性翼下撲過程中流場的等渦量圖如圖14.

圖14 翅翼下撲過程中流場的等渦量圖(左列為柔性翼,右列為剛性翼)Fig.14 Iso-vorticity diagram of the flow field in the plunging process(flexible wing on the left and rigid wing on the right)

從0T到T/2 時刻,后緣渦在后緣遠離邊緣地方逆時針方向向翼面上方內卷,且向翼尖部位匯集.在匯集到翼尖后,上翼面渦達到最強時又逐漸向后緣脫落.此現象和鮑鋒等[35]對于單自由度撲翼渦PIV 實驗研究發現的規律較為一致.這是因為在撲動過程中翼尖的線速度最大,動壓大導致翼尖處靜壓最小,使得渦向翼尖匯集,形成堆積渦從而又影響了翅翼展向上速度分布,在下撲達到最大速度時柔性翼與剛性翼的速度分布呈現了堆積渦形狀如圖15.

圖15 T/4 時刻速度云圖(上為剛性翼,下為柔性翼)Fig.15 Velocity cloud chart at T/4 (rigid wing up and flexible wing down)

從柔性翼流場的等渦量圖可以看出:柔性翼堆積渦現象顯著,從而導致上下翼面速度差增大,根據伯努利原理得出上下翼面的壓強差增大,從而導致升力的明顯增加,但是前緣渦脫落較快,這種現象對于升力具有消極的作用[33],導致柔性翼升力系數下降速度相比于剛性翼較快.

從剛性翼流場的等渦量圖可以看出:在下撲過程中翼根部位始終有渦的存在,此現象和文獻[36]對于三維剛性撲翼在單自由度下的二維截面渦量云圖顯示一致,由于在實際情況中蜻蜓身體的存在,因此在下撲時翼根處向上卷起的渦會給予軀體一個小的升力,但是對于翅翼的升力基本沒有影響.

文獻[37]對于飛蛾懸停實驗的研究表明,撲翼飛行中大約有2/3 的升力取決于前緣渦.因此本文選取楊氏模量為10.4 GPa 均勻分布的柔性翼和剛性翼在升力系數峰值時刻,距離翼根不同位置處(6 個等距截面)的渦量云圖來觀察前緣渦.如圖16 所示:柔性翼的前緣渦從翼根到翼尖逐漸擴散且擴散程度大于剛性翼,且前緣渦的強度在靠近翼根的三分之二位置以內(即靠近翼根的前4 個截面)明顯強于剛性翼.

圖16 T/4 時刻翅翼二維截面渦量云圖(上為剛性翼,下為柔性翼)Fig.16 2D section vorticity cloud at T/4 (rigid wing up and flexible wing down)

根據扭轉角做余弦函數的變化過程可以看出:沿展向柔性翼會發生弦向的周期性扭轉振蕩并改變了原有振蕩的形式,這一振蕩過程會使得前緣渦增強[38].因此導致在撲動速度最大時刻(即升力系數達到峰值時刻)做周期性扭轉振蕩的柔性翼的前緣渦要強于剛性翼,從而相比于剛性翼來說提高了升力系數峰值.

5.2 不同柔性分布下蜻蜓前翼的阻力系數

在均勻柔性分布的情況下,翅翼的阻力系數在不同楊氏模量的變化曲線如圖17 所示.阻力系數曲線隨著楊氏模量的增加:(1)最大阻力系數先減小后增加;(2)最大推力系數逐漸減小;(3)曲線的上下波動的程度逐漸減小,隨后保持基本一樣的波動幅度;(4)推力產生時長先小于剛性翼,隨后逐漸增大.

圖17 均勻柔性分布下柔性翼的阻力系數Fig.17 Drag coefficient of flexible wing with uniform flexible distribution

總體來說,在一定剛度范圍內,剛度的增加可以減小阻力,增加推力及產生推力時長(例如17 GPa,30.2 GPa 和33.6 GPa)這與高鵬聘等[39]對于柔性翼在不同剛度下的水動力特性試驗研究發現的結果較為一致.

在變柔性分布情況下,翅翼的阻力系數如圖18所示,相比于剛性翼來說:(1) 提高推力系數峰值;(2)明顯增加推力在一個周期內的時間占比;(3)并不會較大增加阻力系數峰值.

圖18 變柔性分布下柔性翼的阻力系數Fig.18 Drag coefficient of flexible wing with variable flexibility distribution

相比于均勻柔性翼來說,既獲得了3.8 GPa 均勻柔性分布時的高推力系數峰值,又綜合了17 GPa 均勻柔性分布時的推力在一個周期內大的時間占比.因此表明剛柔并存的蜻蜓前翼在撲動時可以更好的提高推進效率.

為了探究阻力曲線變化中兩次推力峰值產生的原因,將在推力系數峰值時刻對柔性翅翼尾渦運動進行分析.因為三維翅翼在撲動過程中,靠近翅翼后緣的第一對脫落渦距離翅翼最近且比后方的渦量大得多[35],也就對翅翼推力影響較大,所以本文主要研究靠近翅翼后緣的第一對脫落渦.雖然渦的結構為三維,但在一個近軸平面上,渦的旋轉速度主要在流向平面[35],因此距離翼根36 mm的YZ平面上的脫落渦在一定程度上反映了整體渦系變化.

翅翼在下撲和上撲過程中規定渦量 ω 逆時針為正,順時針為負,渦量最大(顏色最亮)為逆時針渦的渦心,渦量最小(顏色最深)為順時針渦的渦心.

從圖19 可以看出,在下撲過程中脫落的渦②與上一周期上撲過程中脫落的渦①在翅翼后緣形成一對反卡門渦街.從圖20 看出,在上撲過程中脫落的渦③與下撲過程脫落的渦②再次形成一對反卡門渦街.在渦③完全脫落時渦①基本消散,這是因為渦在移動過程中的摩擦和碰撞導致渦量消散的較快[30].

圖19 柔性翼下撲推力系數峰值時刻渦量云圖Fig.19 Cloud diagram of vorticity at peak moment of thrust coefficient of flexible wing during upnward flapping

圖20 柔性翼上撲推力系數峰值時刻渦量云圖Fig.20 Cloud diagram of vorticity at peak moment of thrust coefficient of flexible wing during upnward flapping

因此每半個周期的脫落渦均可以和上半周期的脫落渦形成一對反卡門渦街來為蜻蜓前飛提供推力,這也導致了阻力系數曲線中會出現兩次推力系數峰值.

從圖21 中可以看出,在上撲推力系數峰值時刻剛性翼尾緣會產生上下兩對反卡門渦街,此現象和文獻[40]的仿蝌蚪剛性游動的尾跡出現上下兩條反卡門渦街現象較為一致,但兩對反卡門渦街的渦量強度較低,產生的推力相比于柔性翼來說較小.

圖21 剛性翼上撲推力系數峰值時刻渦量云圖Fig.21 Cloud diagram of vorticity at peak moment of thrust coefficient of rigid wing during upnward flapping

從圖22 和圖23在T/4 時刻的流場帶有速度顯示的等渦量圖可以看出來,無論翅翼是柔性還是剛性,均會出現連續的上下交替排列的后緣尾渦向來流方向發展且逐漸減弱.根據張志君等[41]對于5 種翼型(NACA0014,NACA2414,NACA4414,NACA8414,S1223-RTL)的仿生撲翼氣動性能的研究,發現NACA8414 和S1223-RTL 產生了明顯的后緣尾渦,有助于推力的產生.由此可得,因為柔性翼的后緣渦較強,即此刻(推力系數峰值時刻)產生的推力相較于剛性翼來說較大.

圖22 T/4 時刻柔性翼流場的等渦量圖Fig.22 Iso-vorticity diagram of flow field of flexible wing at T/4

圖23 T/4 時刻剛性翼流場的等渦量圖Fig.23 Iso-vorticity diagram of flow field of rigid wing at T/4

5.3 不同柔性分布下蜻蜓前翼獲得的動量增量及加速度

根據動量定理,在一定時間內,外力作用在質點上的沖量,等于質點在此時間內動量的增量.積分形式表達式如下

通過Matlab 編程擬合并求解阻力曲線在一個周期內的積分,積分為負時表示推力作用在翅翼上的沖量大于阻力作用,獲得的動量增量為負,加速度方向和前飛方向一致,加速度為正.翅翼獲得的動量增量以及加速度表達式如下

式中,FD為阻力;T為撲動周期;m代表蜻蜓前翼質量.

根據動量定理和加速度表達式計算得出在不同楊氏模量分布下動量變化和加速度變化的結果如下表2 所示.(1)均勻柔性分布時,隨著楊氏模量的增加動量增量先增加后減小,加速度也先增加后減小;(2)變柔性分布的翅翼相比于剛性翼來說,顯著提高了動量增量以及獲得了較高的加速度;(3)在翅翼做單自由度的撲動運動情況下,線性函數變柔性分布時可以為蜻蜓提供較大的前飛動量增量以及加速度.

表2 蜻蜓前翼周期內的動量增量及加速度Table 2 Momentum increment and acceleration in the period of dragonfly's forewing

5.4 不同柔性分布下蜻蜓前翼在一個周期內時均推力系數

一個周期內的時均推力系數表達式如下

通過Matlab 編程計算得出在不同楊氏模量分布下的時均推力系數如表3 所示.其中時均阻力系數為負值時表示一個周期產生的凈阻力為推力反之為阻力.

表3 蜻蜓前翼不同柔性下的時均推力系數Table 3 Average thrust coefficient of dragonfly forewing under different flexibility

在均勻柔性分布時,楊氏模量較小時,并不會幫助蜻蜓獲得更好的推力,反而增加了阻力.隨著楊氏模量的增加,時均推力系數先增加后減小.這種變化規律與錢靖等[42]對于翼型后緣薄板不同楊氏模量下產生的平均推力系數仿真結果規律一致.

在變柔性分布時,時均推力系數相對于剛性翼來說得到了顯著提高.再次證明了剛柔并存的蜻蜓前翼使得蜻蜓獲得更好的氣動性能.

6 結論

本文基于流固耦合探究了兩種柔性分布方式對于蜻蜓前翼氣動性能的影響.

均勻柔性分布時,楊氏模量較小時,翅翼翼尖的變形趨勢會與楊氏模量較大的均勻柔性分布相反且導致升阻力系數變化趨勢滯后剛性翼半周期.

兩種柔性分布方式下的柔性翼相比于剛性翼來說均可以顯著提高升力系數峰值.其一,由于柔性翼在撲動過程中產生的堆積渦使得上下翼面壓差增大.其二,由于柔性翼撲動過程中會產生扭轉振蕩使得前緣渦增強.

兩種柔性分布方式下的柔性翼相比于剛性翼來說均可以顯著提高推力系數峰值.其一,柔性翼撲動在后緣產生的反卡門渦街的渦量強度要強于剛性翼.其二,柔性翼較強的尾緣渦有助于推力的產生.

在一個周期內,柔性翼在均勻柔性分布下,隨著楊氏模量的增加,給予蜻蜓的動量增量增量先增加后減小,加速度先為減速后為加速且加速大小先增加后減小.在變柔性分布下,柔性翼相對于剛性翼來說,給予蜻蜓的動量增量較大,可以提供給蜻蜓較大的前飛加速度.

時均推力系數在均勻柔性分布下,隨著楊氏模量的增加,先增加后減小(其中楊氏模量較小時的時均推力系數體現為阻力).在變柔性分布下,時均推力系數相比于剛性翼顯著增大,其中線性函數分布時,時均推力系數顯著提高.