趙 磊,張樹海,袁俊杰,張忠海,周 林,何廣平

(1.北方工業(yè)大學(xué), 北京 100144; 2.北京航天測控技術(shù)有限公司, 北京 100041)

0 引言

微型撲翼飛行器具有隱蔽性好、機(jī)動性強(qiáng)等突出優(yōu)點,在軍事偵察、搜救等領(lǐng)域擁有廣泛的應(yīng)用前景,受到了各國國防部門的高度重視[1]。飛行昆蟲以其尺寸小、可懸飛等諸多優(yōu)點成為了微型撲翼飛行器設(shè)計的重要仿照對象[2]。為了設(shè)計出高效、可靠的仿昆蟲撲翼飛行器,學(xué)者們對中、大型昆蟲的空氣動力學(xué)進(jìn)行了大量研究,對于中、大型昆蟲的氣動力機(jī)制有了較為深入的理解。然而,目前對于體長<2 mm、翼展1 mm左右的微小型昆蟲的研究尚處于起步階段,對其氣動力機(jī)制的認(rèn)識還相對有限[3]。微小型昆蟲的振翅頻率通常>200 Hz,其基于弦長的雷諾數(shù)(reynolds number based on chord length,Rec)約為4～60[4-5],并且許多微小型昆蟲的翅膀上長有大量細(xì)長剛毛。剛毛翼的出現(xiàn)為飛行器的進(jìn)一步微型化提供了全新的思路。

拍合運動(合攏-打開,clap-fling)是小昆蟲常用的升力產(chǎn)生機(jī)制[6]。Clap是上行程末期昆蟲的雙翅在其背后相互靠近并繞前緣旋轉(zhuǎn)合攏的過程,fling是合攏的雙翅在下行程初期相互遠(yuǎn)離并繞后緣旋轉(zhuǎn)打開的過程。Miller推測剛毛翼可通過clap-fling在Rec～10時獲得較高的飛行效率[7]。Jones通過二維計算流體動力學(xué)(computational fluid dynamics,CFD)模擬發(fā)現(xiàn),剛毛間隙的流體泄漏可減小左右剛毛翼在fling階段所受阻力[8]。申研究了剛毛交叉對剛毛翼clap-fling過程中氣動特性的影響[9]。吳分析了fling階段剛毛附近的流場結(jié)構(gòu),并探討了剛毛間的相互作用[10]。

由以上研究可知,clap-fling是剛毛翼在懸飛模式下的重要運動方式。目前已證實提高clap-fling中拍動與俯仰重疊時間占clap-fling總時間的比重(拍動-俯仰重疊率)能夠有效提高實心撲翼的平均升力系數(shù)[7]。然而,拍動-俯仰重疊率對剛毛翼升力性能的影響目前尚不清楚。本研究中的主要工作是建立包含clap-fling運動的剛毛翼全撲動周期CFD模型,分析剛毛翼懸飛過程中的黏性效應(yīng),進(jìn)而探討不同拍動-俯仰重疊率對剛毛翼升力系數(shù)的影響。

1 剛毛翼幾何及運動學(xué)模型

1.1 幾何模型

剛毛翼昆蟲翼尖附近的剛毛數(shù)量和剛毛長度通常明顯高于其他區(qū)域[11],因此翼尖區(qū)域?qū)ιω暙I(xiàn)較大,通過截取翼尖附近的弦向截面(圖1(a))能夠近似反映剛毛翼的氣動特性[12]。為此,本研究中根據(jù)吳提供的幾何參數(shù)建立了剛毛翼的二維CFD簡化模型[10]。其中,剛毛翼的二維簡化幾何模型如圖1(b)所示。左右剛毛翼均由21根等間距的剛毛組成,其中每個圓代表一根剛毛的橫截面;d為剛毛直徑;g為相鄰剛毛中心之間的距離;c為弦長,表示前緣剛毛中心至后緣剛毛中心之間的距離;左右剛毛翼的弦線在運動開始前相互平行,初始翼間距為W。以上幾何參數(shù)的取值為:c=0.225 8 mm,d=0.005c,g=0.05c,W=0.2c。

圖1 剛毛翼二維簡化幾何模型的建立Fig.1 Establishment of the two-dimensional geometric model of a pair of bristled wings

1.2 運動學(xué)模型

如圖2所示,clap-fling中包含拍動和俯仰2種運動模式,常以俯仰運動的持續(xù)時間作為clap-fling的總時間。本文中采用以下考慮拍動-俯仰重疊率ξ的運動學(xué)模型[13]:

(1)

(2)

τtr,2=τrot,1+(1-ξ)Δτrot,1

(3)

τtr,6=τrot,4+(1-ξ)Δτrot,4

(4)

式(1)—式(4)中:ω為俯仰角速度;V為拍動線速度;τ為無量綱時間(τ=t/T,其中t為時間,T為撲動周期);i為俯仰角速度曲線的分段編號,i=1～6;j為拍動線速度曲線的分段編號,j=1～9;ωmax為俯仰角速度的最大值;Vmax為拍動線速度的最大值;ai為描述第i段俯仰角速度曲線的常數(shù);bj和cj為描述第j段拍動線速度曲線的常數(shù);τrot,i和Δτrot,i分別為第i段俯仰角速度曲線的起始時刻和所經(jīng)歷時間;τtr, j和Δτtr, j分別為第j段拍動線速度曲線的起始時刻和所經(jīng)歷時間;ωmax、Vmax、ai、bj、cj、τrot,i、Δτrot,i、τtr, j和Δτtr, j等參數(shù)的取值或表達(dá)式詳見文獻(xiàn)[13]。

圖2 Clap-fling過程的二維示意圖Fig.2 Two-dimensional diagrams of clap-fling

運動從下行程開始,首先進(jìn)行fling;初始時刻雙翼互相平行且翼間距最小。圖3為ξ=50%時的無量綱速度曲線。

圖3 ξ=50%時的全撲動周期無量綱速度曲線Fig.3 Dimensionless velocity profile in a complete flapping period at ξ=50%

2 剛毛翼CFD模型

2.1 控制方程與邊界條件

剛毛翼的雷諾數(shù)和馬赫數(shù)很小,屬于典型的不可壓縮層流問題,采用非定常不可壓Navier-Stokes方程進(jìn)行描述。

在懸飛模式下,剛毛翼流場的內(nèi)邊界為剛毛表面所構(gòu)成的運動壁面,其運動規(guī)律由式(1)—式(4)確定,壁面處滿足無滑移條件;外邊界為距離剛毛翼足夠遠(yuǎn)處的遠(yuǎn)場邊界。

2.2 重要無量綱參數(shù)定義

升力系數(shù)CL定義為:

(5)

式(5)中:FL為升力,以豎直向上為正方向;ρ為流體密度;c為弦長。

(6)

基于弦長的雷諾數(shù)Rec定義為:

(7)

式(7)中:ν為流體的運動學(xué)黏度。

2.3 重疊網(wǎng)格方案

重疊網(wǎng)格的基本思想是:對于各個部件附近的流場區(qū)域分別劃分部件網(wǎng)格,對整個流場區(qū)域劃分背景網(wǎng)格,通過允許部件網(wǎng)格與背景網(wǎng)格之間以及不同部件網(wǎng)格之間相互重疊來簡化網(wǎng)格劃分過程,以便在部件周圍生成高質(zhì)量的結(jié)構(gòu)網(wǎng)格和實現(xiàn)網(wǎng)格的參數(shù)化。重疊網(wǎng)格的基本流程如圖4所示:① 通過挖洞處理找出不參與計算的洞單元;② 通過將更多單元轉(zhuǎn)化為洞單元,將不同網(wǎng)格間的重疊區(qū)域最小化;③ 根據(jù)流場信息在重疊區(qū)域查找各套網(wǎng)格的插值接受單元與插值貢獻(xiàn)單元;④ 根據(jù)插值接受單元與貢獻(xiàn)單元之間的位置關(guān)系,插值計算出接受單元的流場信息,作為各套網(wǎng)格獨立求解時的重疊邊界條件;⑤ 根據(jù)重疊邊界條件和物理邊界條件對各套網(wǎng)格進(jìn)行獨立求解。

圖4 重疊網(wǎng)格基本流程

本研究中,采用商業(yè)軟件Fluent 2020R1進(jìn)行剛毛翼的重疊網(wǎng)格劃分和流場計算,其中重疊網(wǎng)格方案如圖5所示。背景網(wǎng)格為半徑為30c的圓形區(qū)域,內(nèi)部劃分結(jié)構(gòu)網(wǎng)格,用于捕捉較遠(yuǎn)處的流場信息;部件網(wǎng)格為半徑為2c的圓形區(qū)域,隨剛毛翼剛性運動,內(nèi)部劃分高質(zhì)量的結(jié)構(gòu)網(wǎng)格,并在剛毛表面附近進(jìn)行局部網(wǎng)格加密,用于捕捉剛毛附近的流場信息。針對以上網(wǎng)格系統(tǒng),通過建立部件網(wǎng)格和背景網(wǎng)格的重疊邊界,并對剛毛表面和剛毛翼部件網(wǎng)格施加運動速度,完成重疊網(wǎng)格的劃分。

圖5 重疊網(wǎng)格方案Fig.5 Overset grid scheme

2.4 數(shù)值求解算法及驗證

速度項采用二階迎風(fēng)格式進(jìn)行離散,時間格式采用一階全隱式時間積分,求解算法采用壓力-速度耦合求解算法。

為確保網(wǎng)格無關(guān)性,對網(wǎng)格總數(shù)為70萬、90萬和110萬時的瞬時升力系數(shù)模擬結(jié)果進(jìn)行了比較,比較結(jié)果如圖6所示。由圖6可知,3種網(wǎng)格數(shù)量下的升力系數(shù)非常相近,但考慮到總網(wǎng)格數(shù)量為110萬時渦量等值線的光滑性較好,本研究中采用網(wǎng)格總數(shù)為110萬的網(wǎng)格劃分方案。

圖6 網(wǎng)格無關(guān)性驗證Fig.6 Grid independence verification

為驗證所建立的剛毛翼二維CFD簡化模型的有效性,針對Rec=10時的fling階段,將本文中模擬得到的升力系數(shù)與吳模擬得到的升力系數(shù)[10]進(jìn)行了比較,比較結(jié)果如圖7所示。

圖7 CFD模型驗證Fig.7 Verification of CFD model

3 計算結(jié)果及分析

3.1 剛毛翼的黏性效應(yīng)

在剛毛翼可能面對的低雷諾數(shù)范圍內(nèi)(Rec=10～80),分析了其黏性效應(yīng)。圖8為ξ=50%時不同雷諾數(shù)下剛毛翼升力系數(shù)隨無量綱時間的變化。從圖8可以看出,當(dāng)Rec=10～80時,升力系數(shù)隨著Rec的增大而明顯減小。

為揭示雷諾數(shù)也即黏性效應(yīng)對剛毛翼升力的影響機(jī)制,分析對比了Rec=10和80時的剛毛受力分布和剛毛翼流場。由于clap階段、fling階段和純拍動階段的情況相似,以下僅展示純拍動階段τ=0.24時刻的比較結(jié)果。圖9-圖12分別為Rec=10和80時的各剛毛升力系數(shù)及其壓力與剪切力分量對比、速度場對比、渦量場對比和壓力場對比。

圖8 ξ=50%時不同雷諾數(shù)下升力系數(shù) 隨無量綱時間的變化Fig.8 Variation of lift coefficient with dimensionless time at ξ=50% for various Rec

圖9 Rec =10和80時純拍動階段τ=0.24時刻的 各剛毛升力系數(shù)及其壓力與剪切力分量對比Fig.9 Comparison of the lift coefficient and its pressure and shear components of each bristle for Rec=10 and 80 at τ=0.24 during translation

由圖9可知,Rec=10和80時剛毛表面的剪切力和壓力均對升力有重要貢獻(xiàn),且各剛毛受力情況與其所在位置密切相關(guān);同時,Rec=10時剛毛表面的剪切力和壓差明顯大于Rec=80時。

圖10 Rec=10和80時純拍動階段τ=0.24時刻的 速度場對比Fig.10 Comparison of the velocity contour for Rec=10 and 80 at τ=0.24 during translation

圖11 Rec =10和80時純拍動階段τ=0.24時刻的 渦量場對比Fig.11 Comparison of the vorticity contour for Rec=10 and 80 at τ=0.24 during translation

由圖10和圖11可知,Rec=10和80時剛毛間隙處均存在一定的正、反間隙渦。Rec=10時,剛毛間隙處整體流速較高,間隙渦相對較弱,說明剛毛間隙處的流體泄漏量較小;Rec=80時,高流速區(qū)域緊貼剛毛表面,在剛毛間隙處存在明顯的低流速區(qū)域,流體與剛毛之間的相對運動速度較大,間隙渦明顯增強(qiáng),說明剛毛間隙處的流體泄漏量較大。由圖12可知,Rec=10時剛毛左下方與右上方之間的壓力差明顯大于Rec=80時。綜合以上分析可知,在剛毛翼最常見的懸飛雷諾數(shù)下(Rec=10),黏性效應(yīng)較強(qiáng),剛毛間隙泄漏明顯減弱,使得剛毛翼內(nèi)外側(cè)連通性較差、壓差較大,導(dǎo)致壓力對升力系數(shù)的貢獻(xiàn)明顯大于Rec=80時。

圖12 Rec =10和80時純拍動階段τ=0.24時刻的 壓力場對比Fig.12 Comparison of the pressure contour for Rec=10 and 80 at τ=0.24 during translation

為闡明Rec=10時剪切力對升力系數(shù)的貢獻(xiàn)明顯大于Rec=80時的原因,進(jìn)行了以下量綱分析:

(8)

式(8)中:Cshear為剛毛表面的無量綱剪切力;μ為流體的動力學(xué)黏度;Vleak為流體通過剛毛間隙泄漏時相對于剛毛的平均流速;V∞為剛毛的絕對運動速度;Cleak表示剛毛間隙處的流體泄漏率,為雷諾數(shù)的函數(shù),定義為Cleak=Vleak/V∞;Red為基于剛毛直徑的雷諾數(shù),定義為Red=ρV∞d/μ。

由式(8)可以看出,剛毛表面的無量綱剪切力Cshear總體上正比于Cleak/Red。在剛毛翼常見的剛毛間隙(g～10d)和可能面對的雷諾數(shù)(Rec=10～80,Red=0.043～0.35)下,根據(jù)兩平行圓柱繞流理論[14],盡管間隙泄漏率Cleak隨著Red的增大而增大,但其增加速率逐漸減小,因此Cleak/Red和Cshear隨著Rec的增大而不斷減小,Rec=10時剛毛表面的無量綱剪切力和剪切力對升力系數(shù)的貢獻(xiàn)明顯大于Rec=80時。

3.2 拍動-俯仰重疊率對升力系數(shù)的影響

為研究剛毛翼常見雷諾數(shù)下(Rec=10)拍動-俯仰重疊率ξ對升力系數(shù)的影響,在Rec=10下,分別模擬了ξ=0、25%、50%、75%、100%時的剛毛翼懸飛過程。其中,平均升力系數(shù)隨拍動-俯仰重疊率的變化和不同拍動-俯仰重疊率下瞬時升力系數(shù)隨無量綱時間的變化分別如圖13和圖14所示。

圖13 平均升力系數(shù)隨拍動-俯仰重疊率的變化

圖14 不同拍動-俯仰重疊率下升力系數(shù) 隨無量綱時間的變化Fig.14 Variation of lift coefficient with dimensionless time at various ξ

由圖13可以看出,平均升力系數(shù)隨著拍動-俯仰重疊率的增大而增大。結(jié)合圖14可知,造成這一現(xiàn)象的主要原因在于:當(dāng)clap和fling階段的拍動-俯仰重疊率較低時(ξ=0和25%),會在clap和fling中期形成負(fù)升力峰值;隨著拍動-俯仰重疊率的不斷增大,以上負(fù)升力峰值逐漸消失(ξ= 50%),并最終轉(zhuǎn)變?yōu)檎Ψ逯?ξ= 75%和100%)。

為進(jìn)一步闡明上述負(fù)升力峰值的轉(zhuǎn)變機(jī)制,選取fling階段和clap階段在ξ=0時的負(fù)升力峰值時刻τ=0.08和τ=0.9,分析對比了ξ=0和100%時的流場。由于流場左右對稱,以下僅選取左翼進(jìn)行分析。圖15—圖19分別為ξ=0和100%時在fling階段τ=0.08時刻的升力系數(shù)及其壓力與剪切力分量對比、速度場對比、渦量場對比、壓力場對比以及剪切力豎直分量的分布情況對比。圖20—圖24分別為ξ=0和100%時在clap階段τ=0.9時刻的升力系數(shù)及其壓力與剪切力分量對比、速度場對比、渦量場對比、壓力場對比以及剪切力豎直分量的分布情況對比。

圖15 ξ=0和100%時在fling階段τ=0.08時刻的 升力系數(shù)及其壓力與剪切力分量對比Fig.15 Comparison of the lift coefficient and its pressure and shear components of each bristle for ξ=0 and 100% at τ=0.08 during fling

圖16 ξ=0和100%時在fling階段τ=0.08時刻的 速度場對比Fig.16 Comparison of the velocity contour for ξ=0 and 100% at τ=0.08 during fling

圖17 ξ=0和100%時在fling階段τ=0.08時刻的 渦量場對比Fig.17 Comparison of the vorticity contour for ξ=0 and 100% at τ=0.08 during fling

由圖15—圖18可知,當(dāng)ξ=0時,流體通過剛毛間隙向左側(cè)泄漏時分別在剛毛左下方和右上方形成了局部低壓區(qū)和局部高壓區(qū),因此壓差產(chǎn)生的合力指向左下方,對升力的貢獻(xiàn)為負(fù);而當(dāng)ξ=100%時,流體通過剛毛間隙向右側(cè)泄漏時在剛毛左下方形成了局部高壓區(qū),而在右上方形成了局部低壓區(qū),壓差產(chǎn)生的合力指向右上方,對升力的貢獻(xiàn)為正。此外,結(jié)合圖15和圖19可知,對于fling中期,當(dāng)ξ=0時,剛毛左上方和右下方的剪切力豎直分量均為負(fù),導(dǎo)致剪切力對升力的貢獻(xiàn)為負(fù);而當(dāng)ξ=100%時,前緣附近剛毛剪切力對升力的貢獻(xiàn)為正,而其他位置剛毛剪切力對升力的貢獻(xiàn)為負(fù)。

圖18 ξ=0和100%時在fling階段τ=0.08時刻的 壓力場對比Fig.18 Comparison of the pressure contour for ξ=0 and 100% at τ=0.08 during fling

圖19 ξ=0和100%時在fling階段τ=0.08時刻的 剪切力豎直分量分布情況對比Fig.19 Comparison of the vertical wall shear distribution for ξ=0 and 100% at τ=0.08 during fling

對于clap階段中期(Rec=10),由圖20可知,當(dāng)ξ=0時,壓力和剪切力對升力的貢獻(xiàn)相當(dāng),且均為負(fù)值;隨著ξ的不斷增加,當(dāng)ξ=100%時,壓力仍對升力存在一定的負(fù)向貢獻(xiàn),但剪切力對升力的正向貢獻(xiàn)更為顯著,導(dǎo)致剛毛產(chǎn)生了明顯的正升力。

圖20 ξ=0和100%時在clap階段τ=0.9時刻的 各剛毛升力系數(shù)及其壓力與剪切力分量對比Fig.20 Comparison of the lift coefficient and its pressure and shear components of each bristle for ξ=0 and 100% at τ=0.9 during clap

由圖21和圖22可知,對于clap階段中期,當(dāng)ξ=0時,流場中的渦以剛毛間隙渦為主,流體通過剛毛間隙向右泄漏;當(dāng)ξ=100%時,在拍動和俯仰運動的耦合作用下,左右翼加速閉合,在兩翼之間產(chǎn)生的增壓作用較為顯著,使得流場中同時存在明顯的逆時針前緣渦、順時針后緣渦和明顯強(qiáng)于ξ=0時的剛毛間隙渦,流體通過剛毛間隙向左泄漏。

圖21 ξ=0和100%時在clap階段τ=0.9時刻的 速度場對比Fig.21 Comparison of the velocity contour for ξ=0 and 100% at τ=0.9 during clap

圖22 ξ=0和100%時在clap階段τ=0.9時刻的 渦量場對比Fig.22 Comparison of the vorticity contour for ξ=0 and 100% at τ=0.9 during clap

由圖23可知,ξ=0時,剛毛左上方的高壓區(qū)導(dǎo)致了負(fù)升力的產(chǎn)生;ξ=100%時,流體通過剛毛間隙泄漏時在剛毛右側(cè)形成了局部高壓區(qū),且高壓區(qū)隨著剛毛位置遠(yuǎn)離前緣而向上移動,導(dǎo)致壓力對大部分剛毛的升力貢獻(xiàn)為負(fù)值。與此同時,由圖24可知,對于clap中期,ξ=0時剛毛表面大部分區(qū)域的剪切力豎直分量為負(fù)值,導(dǎo)致剪切力對升力的貢獻(xiàn)為負(fù);而ξ=100%時剛毛表面大部分區(qū)域的剪切力豎直分量為正值,導(dǎo)致剪切力對升力的貢獻(xiàn)為正。

圖23 ξ=0和100%時在clap階段τ=0.9時刻的 壓力場對比Fig.23 Comparison of the pressure contour for ξ=0 and 100% at τ=0.9 during clap

圖24 ξ=0和100%時在clap階段τ=0.9時刻的 剪切力豎直分量分布情況對比Fig.24 Comparison of the vertical wall shear distribution for ξ=0 and 100% at τ=0.9 during clap

4 結(jié)論

1)Rec=10～80時,隨著雷諾數(shù)的降低,黏性效應(yīng)明顯增強(qiáng),剛毛間隙泄漏和剛毛翼內(nèi)外側(cè)連通性明顯減弱,導(dǎo)致剛毛翼內(nèi)外側(cè)壓差增大,壓力對升力系數(shù)的貢獻(xiàn)明顯增強(qiáng);與此同時,剛毛表面的無量綱剪切力同剛毛間隙泄漏率與雷諾數(shù)之比成正比,與雷諾數(shù)呈負(fù)相關(guān),導(dǎo)致無量綱剪切力和剪切力對升力系數(shù)的貢獻(xiàn)隨著雷諾數(shù)的降低而明顯增大。

2) 當(dāng)拍動和俯仰運動無重疊時,clap和fling階段中期的渦以成對出現(xiàn)且方向相反的剛毛間隙渦為主,流體泄漏方向與剛毛翼運動方向相同,泄漏過程中產(chǎn)生的剪切力和壓力對升力的貢獻(xiàn)為負(fù),從而形成了明顯的負(fù)升力峰值。

3) 當(dāng)拍動-俯仰重疊率為100%時,clap和fling階段中期的流體泄漏方向與剛毛翼運動方向相反,流場中同時存在明顯的前緣渦、后緣渦和剛毛間隙渦;其中,clap階段在流體泄漏過程中產(chǎn)生了明顯的豎直向上剪切力合力,而fling階段則在流體泄漏過程中產(chǎn)生了明顯的豎直向上壓差,導(dǎo)致clap和fling階段的負(fù)升力峰值轉(zhuǎn)變?yōu)檎Ψ逯?平均升力系數(shù)明顯增大。