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

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

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

引言

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

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

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

1 撲翼模型與函數(shù)

1.1 幾何模型建立

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

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

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

1.2 撲動函數(shù)

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

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

根據(jù)撲動角度函數(shù)先轉(zhuǎn)化為弧度制再進(jìn)行求導(dǎo),得到其撲動角速度函數(shù)

1.3 柔性分布函數(shù)

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

在變?nèi)嵝苑植贾卸x蜻蜓前翼前緣楊氏模量17 GPa[27],定義后緣楊氏模量3.8 GPa[30].通過Position 場函數(shù),調(diào)用前翼隨體坐標(biāo)系中y坐標(biāo).因此前后緣的y坐標(biāo)和楊氏模量構(gòu)成兩對坐標(biāo)點,得出(0.002,17)和(-0.006,3.8).線性函數(shù)表達(dá)式通過兩點即可得到,將二次函數(shù)頂點置于前緣處,并求解通過兩點的函數(shù),即可得到平方函數(shù)的表達(dá)式.整個蜻蜓前翼根據(jù)柔性分布函數(shù)發(fā)生柔性變化:

均勻柔性分布

線性函數(shù)柔性分布

平方函數(shù)柔性分布

2 研究方法

2.1 計算方法

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

剛性翼的撲動研究較較為簡單,主要采用動網(wǎng)格技術(shù)[30]實現(xiàn)對蜻蜓前翼撲動的仿真計算.柔性翼由于變形較大導(dǎo)致網(wǎng)格劇烈變化,因此在流固耦合的基礎(chǔ)上采用重疊網(wǎng)格和網(wǎng)格自適應(yīng)技術(shù).

將外流場和重疊區(qū)域均設(shè)置為流場域,流場域選擇隱式非穩(wěn)態(tài)求解,湍流模型選擇Spalar-Allmaras湍流模型[32].在流場域中,設(shè)置重疊保護(hù)和單元網(wǎng)格質(zhì)量校正,來提高撲動過程中的流域網(wǎng)格質(zhì)量,重疊區(qū)域設(shè)置為浮動,可以跟隨前翼固體模型一起運(yùn)動.

將前翼結(jié)構(gòu)設(shè)置為固體區(qū)域,在固體域中采用隱式非定常固體應(yīng)力求解模型,添加撲動函數(shù)和楊氏模量函數(shù)使蜻蜓前翼實現(xiàn)運(yùn)動及柔性分布.

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

2.2 控制方程

流體控制方程包含連續(xù)性方程和運(yùn)動方程,其中連續(xù)性方程如下

運(yùn)動方程

式中,ρ 為密度;u代表流體微團(tuán)的速度矢量;F表示作用在單位質(zhì)量上的質(zhì)量力;P為應(yīng)力張量.

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

式中,ρ 是固體密度;σ 是柯西應(yīng)力張量;F是體積力矢量;d是固體域當(dāng)?shù)丶铀俣仁噶?流固耦合需要流體域和固體域之間交換數(shù)據(jù),該過程必須滿足一些物理量的相等和守恒,即滿足方程如下

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

3 計算前處理

3.1 流場區(qū)域網(wǎng)格劃分

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

圖3 流體區(qū)域尺寸及網(wǎng)格Fig.3 The size of fluid domain and mesh

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

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

3.2 固體域網(wǎng)格劃分

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

圖4 蜻蜓前翼平板模型網(wǎng)格Fig.4 Mesh of forewing plat model

3.3 數(shù)值方法驗證

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

圖5 升力系數(shù)隨時間變化規(guī)律與文獻(xiàn)[34]對比結(jié)果Fig.5 The variation law of lift coefficient with time is compared to the results in Ref.[34]

3.4 網(wǎng)格無關(guān)性驗證

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

圖6 網(wǎng)格無關(guān)性驗證Fig.6 Verification of independence of grid number

4 不同柔性分布下柔性翼結(jié)構(gòu)變形結(jié)果與分析

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

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

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

在變?nèi)嵝苑植紬l件下,在撲動過程中蜻蜓前翼變形規(guī)律與楊氏模量較大時的翅翼變形規(guī)律保持一致,如圖7(b)所示.

4.2 不同柔性分布下柔性翼結(jié)構(gòu)變形數(shù)值監(jiān)測方法

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

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

4.3 不同柔性分布下柔性翼數(shù)值變形結(jié)果及流場分析

翼尖點F 的監(jiān)測變形量隨時間的變化如圖9 所示,可以發(fā)現(xiàn):(1)均勻柔性較大的前翼翼尖變形和均勻柔性較小、變?nèi)嵝苑植嫉囊砑庾冃蔚囊?guī)律相反,和前文的撲動變形過程相對應(yīng);(2)變?nèi)嵝苑植紩r,二次函數(shù)分布時翼尖點最大變形量較小,這是因為二次函數(shù)柔性分布時在靠近前緣部位,楊氏模量的減小速度較慢,因此到達(dá)翼尖位置時楊氏模量大于一次函數(shù)柔性分布.

圖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 均勻柔性分布的壓力云圖可以看出,上翼面壓強(qiáng)大于下翼面壓強(qiáng)是導(dǎo)致翼尖變形量在初始時刻為負(fù)位移的原因.

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

可見柔性過大會導(dǎo)致翼尖會發(fā)生相反的柔性變形規(guī)律.

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

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

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

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

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

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

在變?nèi)嵝苑植紩r,翼尖部位的扭轉(zhuǎn)幅度大于翼中部位,說明翼尖部位后緣的上翹和下彎變形略微大于翼中部位的后緣.

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

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

從升力系數(shù)變化規(guī)律如圖13 所示,對比可以發(fā)現(xiàn):(1)在均勻柔性分布時,柔性較大的翅翼升力系數(shù)曲線變化趨勢滯后于剛性翼半周期,這是因為從結(jié)構(gòu)變形來看,將下?lián)溥\(yùn)動離散之后,每個瞬時時刻剛性翼相對于來流做向下?lián)鋭舆\(yùn)動,柔性較大的前翼相對于來流來說做向上變形運(yùn)動,從而導(dǎo)致了半個周期的剛性翼向下運(yùn)動,變?yōu)榱讼蛏献冃芜\(yùn)動;(2) 在變?nèi)嵝苑植紩r,柔性翼和剛性翼的變化趨勢相同.

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

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

為了探究柔性翼如何提高升力系數(shù)峰值,本文選取線性函數(shù)柔性分布的柔性翼和剛性翼下?lián)溥^程中流場的等渦量圖如圖14.

圖14 翅翼下?lián)溥^程中流場的等渦量圖(左列為柔性翼,右列為剛性翼)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 時刻,后緣渦在后緣遠(yuǎn)離邊緣地方逆時針方向向翼面上方內(nèi)卷,且向翼尖部位匯集.在匯集到翼尖后,上翼面渦達(dá)到最強(qiáng)時又逐漸向后緣脫落.此現(xiàn)象和鮑鋒等[35]對于單自由度撲翼渦PIV 實驗研究發(fā)現(xiàn)的規(guī)律較為一致.這是因為在撲動過程中翼尖的線速度最大,動壓大導(dǎo)致翼尖處靜壓最小,使得渦向翼尖匯集,形成堆積渦從而又影響了翅翼展向上速度分布,在下?lián)溥_(dá)到最大速度時柔性翼與剛性翼的速度分布呈現(xiàn)了堆積渦形狀如圖15.

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

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

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

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

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

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

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

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

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

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

在變?nèi)嵝苑植记闆r下,翅翼的阻力系數(shù)如圖18所示,相比于剛性翼來說:(1) 提高推力系數(shù)峰值;(2)明顯增加推力在一個周期內(nèi)的時間占比;(3)并不會較大增加阻力系數(shù)峰值.

圖18 變?nèi)嵝苑植枷氯嵝砸淼淖枇ο禂?shù)Fig.18 Drag coefficient of flexible wing with variable flexibility distribution

相比于均勻柔性翼來說,既獲得了3.8 GPa 均勻柔性分布時的高推力系數(shù)峰值,又綜合了17 GPa 均勻柔性分布時的推力在一個周期內(nèi)大的時間占比.因此表明剛?cè)岵⒋娴尿唑亚耙碓趽鋭訒r可以更好的提高推進(jìn)效率.

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

翅翼在下?lián)浜蜕蠐溥^程中規(guī)定渦量 ω 逆時針為正,順時針為負(fù),渦量最大(顏色最亮)為逆時針渦的渦心,渦量最小(顏色最深)為順時針渦的渦心.

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

圖19 柔性翼下?lián)渫屏ο禂?shù)峰值時刻渦量云圖Fig.19 Cloud diagram of vorticity at peak moment of thrust coefficient of flexible wing during upnward flapping

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

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

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

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

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

圖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 不同柔性分布下蜻蜓前翼獲得的動量增量及加速度

根據(jù)動量定理,在一定時間內(nèi),外力作用在質(zhì)點上的沖量,等于質(zhì)點在此時間內(nèi)動量的增量.積分形式表達(dá)式如下

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

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

根據(jù)動量定理和加速度表達(dá)式計算得出在不同楊氏模量分布下動量變化和加速度變化的結(jié)果如下表2 所示.(1)均勻柔性分布時,隨著楊氏模量的增加動量增量先增加后減小,加速度也先增加后減小;(2)變?nèi)嵝苑植嫉某嵋硐啾扔趧傂砸韥碚f,顯著提高了動量增量以及獲得了較高的加速度;(3)在翅翼做單自由度的撲動運(yùn)動情況下,線性函數(shù)變?nèi)嵝苑植紩r可以為蜻蜓提供較大的前飛動量增量以及加速度.

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

5.4 不同柔性分布下蜻蜓前翼在一個周期內(nèi)時均推力系數(shù)

一個周期內(nèi)的時均推力系數(shù)表達(dá)式如下

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

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

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

在變?nèi)嵝苑植紩r,時均推力系數(shù)相對于剛性翼來說得到了顯著提高.再次證明了剛?cè)岵⒋娴尿唑亚耙硎沟抿唑勋@得更好的氣動性能.

6 結(jié)論

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

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

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

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

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

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