劉旭博，王芬芬，于紀(jì)言

(1.中國飛行試驗(yàn)研究院，西安 710089; 2.南京理工大學(xué)，南京 210094)

撲翼飛行器同時(shí)具有隱蔽性與機(jī)動(dòng)性兩種特征，符合現(xiàn)代戰(zhàn)爭隱蔽偵察與精確打擊的概念，具有良好的軍事用途前景[1]，因此以美國為首的多個(gè)國家[2]，近二十幾年來掀起了一股撲翼飛行器研究的浪潮[3]。隨著計(jì)算機(jī)技術(shù)的發(fā)展，研究撲翼飛行器的方法也愈來愈多。本文在前人對(duì)鳥類撲翼結(jié)構(gòu)研究的基礎(chǔ)上[4]，使用計(jì)算流體力學(xué)方法，建立三維鴿狀剛性翅與柔性翅翼展模型，著重通過對(duì)比不同翼展的氣動(dòng)力參數(shù)，給出柔性翅撲動(dòng)過程中的氣動(dòng)特性。

1 撲翼模型與撲動(dòng)函數(shù)

1.1 撲動(dòng)簡介

自然界鳥類與昆蟲的撲動(dòng)方式各異，產(chǎn)生的影響也各不相同，但經(jīng)過分解，撲動(dòng)過程的運(yùn)動(dòng)形式主要有以下4種：① 翅翼繞著翅根的上下?lián)鋭?dòng); ② 翅面的俯仰或翅翼展向的扭轉(zhuǎn); ③ 翅翼的前后向運(yùn)動(dòng)； ④ 翅膀沿展向彎曲折疊。其運(yùn)動(dòng)形式如圖1所示。

圖1 撲翼運(yùn)動(dòng)的4種運(yùn)動(dòng)形式示意圖

而鳥類的撲翼運(yùn)動(dòng)主要由翅翼的上下?lián)鋭?dòng)與翅面的俯仰運(yùn)動(dòng)組成，因此本文也以這兩種運(yùn)動(dòng)的特征作為數(shù)值計(jì)算的輸入條件。

1.2 剛性翅模型

剛性翅模型的建立主要參照鴿子的翅翼形狀及其特點(diǎn)，典型鴿子的飛行參數(shù)與翅翼特征在表1中列出。

表1 翅翼的特征參數(shù)

根據(jù)典型鴿子翅翼特征參數(shù)與具體鴿子的翅翼形狀，進(jìn)行翅翼建模，模型如圖2所示。

圖2 鴿翼與剛性翅模型

根據(jù)鴿翅的生物學(xué)構(gòu)造，鴿翅的前緣部分主要由肱骨支撐，量得肱骨直徑約為0.01 m。因此仿生翅的前緣半徑選擇為0.005 m。翅翼的后緣部分為羽毛尖端，厚度可以忽略不計(jì)，因此構(gòu)建模型時(shí)后緣厚度為0。鴿翅的厚度最大處位于肱骨與副羽的相接部分，約在弦向20%的位置，厚度為0.016 m。翼展前緣長度約0.13 m，后緣長度為0.3 m。鴿翅的下表面展開時(shí)沒有明顯的弧度，為平面型下表面。具體參數(shù)見表1所示。

1.3 柔性翅模型

建立柔性翅模型時(shí)考慮到后期與剛性翅氣動(dòng)特性的對(duì)比，保持柔性翅的外形與剛性翅完全相同，改變材料屬性即可。同時(shí)，為了探究柔性翅支架對(duì)其氣動(dòng)特性的影響，這里建立2種不同的柔性翅模型：① 橡膠翅翼，0.001 m厚鋼架支撐： ② 橡膠翅翼，0.002 m厚鋼架支撐：

其中模型(1)、模型(2)的外形與剛性翅并無差異，只是改變了其材料屬性，其外形如圖3所示，并且多出一個(gè)支架部分對(duì)翅翼起到支撐作用。模型(1)、模型(2)的差異則主要在于鋼架的厚度不同。

圖3 柔性翅翼及其支架示意圖

對(duì)翅翼建模完成后，構(gòu)建外流域，外流域大小與流域流體的流速有關(guān)，由于是低速飛行，雷諾數(shù)較小，流場(chǎng)變化范圍不大，因此取外流域的長寬大約是翅翼的6倍即可。構(gòu)建后的外流域如圖4所示。

圖4 計(jì)算域

1.4 撲動(dòng)函數(shù)

確定撲動(dòng)函數(shù)時(shí)主要考慮的因素為撲翼飛行器的撲動(dòng)幅值以及撲動(dòng)頻率，由于翼型是根據(jù)鴿子的翅翼確定的，為了使結(jié)果具有實(shí)際參考價(jià)值，撲動(dòng)范圍以及撲動(dòng)頻率也依照鴿子的飛行特點(diǎn)進(jìn)行確定，鴿子的撲動(dòng)頻率在8 Hz左右，鴿子的主翅翼的撲扇幅度在45°左右，具體設(shè)置函數(shù)時(shí)，采用 8 Hz 的撲動(dòng)頻率，45°作為撲動(dòng)幅值。由于撲動(dòng)函數(shù)的確定過程較為簡單，方法成熟[5]，這里直接給出撲動(dòng)函數(shù)為：

(1)

求導(dǎo)得到其角速度函數(shù)：

(2)

2 計(jì)算前處理

2.1 計(jì)算方法

單純剛性翅的計(jì)算方法較為簡單，計(jì)算過程主要采用動(dòng)網(wǎng)格技術(shù)[6]，這里不再過多贅述，其具體過程可參看文獻(xiàn)[6]。柔性翅的計(jì)算則涉及流固耦合技術(shù)，流固耦合的核心在于計(jì)算時(shí)既要解得流場(chǎng)相關(guān)的物理量又要解得固體力學(xué)場(chǎng)的相關(guān)物理量，且不能忽略二者的相互作用[7]。簡單地說，就是在每一步求解中，都需要解兩類不同的方程并完成數(shù)據(jù)傳遞。因此流固耦合的計(jì)算量遠(yuǎn)遠(yuǎn)大于一般單項(xiàng)計(jì)算的情況，但其結(jié)果也是可觀的，其求解出的結(jié)果往往更接近真實(shí)的物理現(xiàn)象[8]。本文采用雙向耦合技術(shù)，即在翅翼在流場(chǎng)作用下發(fā)生形變，形變后翅翼的相關(guān)參數(shù)傳給流體域，流體域接收流體中翅翼形變數(shù)據(jù)，根據(jù)翅翼的新狀態(tài)進(jìn)行下一時(shí)間步的流場(chǎng)參數(shù)計(jì)算。

2.2 計(jì)算模型

在本文的計(jì)算過程中，剛性翅的計(jì)算模型較為簡單，其牽涉的計(jì)算模型只是柔性翅的一部分，而柔性翅計(jì)算時(shí)，既包含流體域計(jì)算，也包含固體力學(xué)計(jì)算，另外還涉及兩者之間的數(shù)據(jù)交換，下面對(duì)柔性翅的計(jì)算模型作以詳細(xì)說明：

1) 流體控制方程：

流體控制方程包括連續(xù)性方程與運(yùn)動(dòng)方程。連續(xù)性方程[9]如下。

(3)

運(yùn)動(dòng)方程：

(4)

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

2) 固體控制方程：

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

ρd=▽·σ+F

(5)

式中：ρ是固體密度；σ是柯西應(yīng)力張量；F是體積力矢量；d是固體域當(dāng)?shù)丶铀俣仁噶俊?/p>

3) 流固耦合方程：

流固耦合的過程即為數(shù)據(jù)傳遞與應(yīng)用的過程，該過程中則必須滿足流體域與固體域的一些重要物理量的相等與守恒，即滿足如下方程：

(6)

式中：τf、τs分別代表流體域與固體域力大??；nf、ns則接觸表面的法向單位向量；df、ds代表流體域與固體域微團(tuán)位移矢量；qf、qs分別代表流體域與固體域的熱流量；Tf、Ts分別代表流體域與固體域的溫度。

第一項(xiàng)則表示交界面上的流體域與固體域力相等，第二項(xiàng)表示位移相等，第三項(xiàng)表示熱流量相等，第四項(xiàng)表示溫度相等。

2.3 計(jì)算平臺(tái)的搭建

ANSYS進(jìn)行流固耦合計(jì)算時(shí)，在workbench中搭建如圖5所示的平臺(tái)。

圖5 計(jì)算平臺(tái)示意圖

其中g(shù)eometry模塊用于翅翼模型的建立并共享給后面的兩個(gè)計(jì)算模塊。計(jì)算過程為瞬態(tài)過程，結(jié)構(gòu)體的力學(xué)計(jì)算使用Transient structure，流體部分CFX與Fluent都具有流固耦合的能力，但Fluent進(jìn)行動(dòng)網(wǎng)格計(jì)算時(shí)更有優(yōu)勢(shì)，因此這里流體部分使用Fluent模塊，最后采用耦合器system coupling。將結(jié)構(gòu)體與流體的setup部分與system coupling相聯(lián)接，計(jì)算過程中兩者在耦合器中進(jìn)行數(shù)據(jù)傳遞。

2.4 計(jì)算域網(wǎng)格劃分

平臺(tái)搭建完成后，需分別對(duì)流體與固體計(jì)算域進(jìn)行網(wǎng)格劃分。網(wǎng)格劃分時(shí)，需要根據(jù)不同網(wǎng)格的特點(diǎn)進(jìn)行考慮：六面體網(wǎng)格可以根據(jù)計(jì)算模型的結(jié)構(gòu)進(jìn)行塊劃分后生成，其計(jì)算速度快，準(zhǔn)確性高，但其不適用于形狀復(fù)雜的幾何模型且在動(dòng)網(wǎng)格中適應(yīng)性差，不宜收斂。四面體網(wǎng)格生成速度快，相對(duì)內(nèi)存占用較大，但適用于邊界條件有巨大變化的動(dòng)網(wǎng)格情況。因此這里采用四面體網(wǎng)格的劃分策略。

流體域的網(wǎng)格劃分可用Workbench中自帶的mesh模塊進(jìn)行過操作，也可聯(lián)接ICEM模塊進(jìn)行網(wǎng)格劃分，之后將網(wǎng)格直接導(dǎo)入mesh模塊。結(jié)構(gòu)力學(xué)部分的網(wǎng)格劃分過程相對(duì)簡單，選擇對(duì)應(yīng)的body，輸入控制網(wǎng)格大小參數(shù)即可，參數(shù)設(shè)置時(shí)盡量保持與流體接觸面網(wǎng)格大小相似，這樣能夠適當(dāng)減小插值次數(shù)，加快運(yùn)算速度。網(wǎng)格仍然選取四面體網(wǎng)格，因?yàn)槌嵋頃?huì)發(fā)生明顯變形，六面體網(wǎng)格不適合復(fù)雜形狀結(jié)構(gòu)體的運(yùn)算。具體劃分結(jié)果如圖6所示。

圖6 網(wǎng)格劃分情況

3 計(jì)算過程

3.1 翅翼撲動(dòng)升、阻力變化特點(diǎn)

在FLUENT中設(shè)置監(jiān)測(cè)時(shí)，一般得到的是升阻力系數(shù)cl、cd。無量綱量升阻力系數(shù)可以對(duì)比分析升阻力大小。

考慮到鴿子的飛行速度大致在11 m/s，且鳥類在實(shí)際飛行時(shí)很少具有絕對(duì)的平飛狀態(tài)，因此在FLUENT中設(shè)置氣體來流速度11 m/s，設(shè)置攻角為8°，可以監(jiān)測(cè)到剛性翅以及兩種柔性翅升阻力系數(shù)的變化情況。

從圖7、圖8可以看出，0.001 m支架柔性翅的升力、阻力系數(shù)變化失去了剛性翅時(shí)域上的規(guī)律性，但從周期性的變化規(guī)律來看，趨勢(shì)卻是與剛性翅類似的，無論是升力還是阻力系數(shù)，都符合剛性翅“高-低-高”的變化特征。另外在時(shí)域上，一個(gè)周期內(nèi)正升力部分所占的時(shí)域?qū)挾纫黠@大于負(fù)升力。阻力系數(shù)的分布于升力部分也是類似的。

從圖9、圖10可以看出，0.002 m支架柔性翅翅翼的升阻力系數(shù)的峰值得到了大幅度的提升，另一方面 2 mm柔性翅翅翼的升阻力系數(shù)與 1 mm柔性翅翅翼的升阻力系數(shù)有類似的特點(diǎn)，即一個(gè)周期內(nèi)整個(gè)變化曲線在時(shí)域上分布不均勻，依然出現(xiàn)了正的升阻力系數(shù)的時(shí)域范圍明顯大于負(fù)升阻力系數(shù)的時(shí)域。

圖7 剛性翅與0.001 m支架柔性翅升力系數(shù)曲線

圖8 剛性翅與0.001 m支架柔性翅阻力系數(shù)曲線

圖9 剛性翅與0.002 m支架柔性翅升力系數(shù)曲線

圖10 剛性翅與0.002 m支架柔性翅阻力系數(shù)曲線

3.2 波動(dòng)時(shí)域特征

剛性翅與柔性翅的升阻力系數(shù)都呈周期性變化，柔性翅的峰值到來的時(shí)間略晚于剛性翅，且峰值大于剛性翅。另外，周期內(nèi)正升阻力部分所占的時(shí)域?qū)挾纫黠@大于負(fù)升阻力。這種現(xiàn)象對(duì)于平均升力系數(shù)來說，明顯是有益的，下面對(duì)產(chǎn)生該種現(xiàn)象的原因作以分析，首先以兩種柔性翅的升力為例，對(duì)升力曲線的時(shí)域作以劃分，如圖11、圖12所示。

圖11 0.002 m支架柔性翅升力系數(shù)曲線

圖12 0.001 m支架柔性翅升力系數(shù)曲線

對(duì)比剛性翅與0.002 m支架柔性翅的時(shí)域可以發(fā)現(xiàn)，兩者的上升與下降趨。勢(shì)基本是同步的。而出現(xiàn)時(shí)域不對(duì)稱的原因，主要是由于柔性翅部分出現(xiàn)了波動(dòng)時(shí)域，波動(dòng)時(shí)域的長度一定程度上“占用了”柔性翅的負(fù)時(shí)域，從而出現(xiàn)了負(fù)時(shí)域時(shí)域變窄的情況。而在0.001 m柔性翅中也同樣出現(xiàn)了波動(dòng)時(shí)域，在圖12中可以看出，波動(dòng)時(shí)域的長度要大于0.002 m支架柔性翅的波動(dòng)時(shí)域長度，波動(dòng)時(shí)域不但減小了柔性翅負(fù)時(shí)域的長度，且波動(dòng)時(shí)域的延長直接導(dǎo)致了負(fù)升力的出現(xiàn)時(shí)刻發(fā)生了延遲。下面具體分析出現(xiàn)波動(dòng)時(shí)域的成因與原理。

首先，依據(jù)柔性翅的升力系數(shù)變化圖線，可知0.001 m支架柔性翅升力圖線波動(dòng)時(shí)域長度在0.06～0.08 s，0.002 m支架升力柔性翅圖線在0.05 s與0.10 s附近，下面給出兩種柔性翅翼在x=0.09 m處ZOY截面上的速度矢量圖。

結(jié)合1.3節(jié)中的撲翼角度函數(shù)，分析兩種柔性翅t(yī)=0.05 s時(shí)的狀態(tài)。t=0.05 s時(shí)翅翼剛剛經(jīng)過下極限位置，整個(gè)翅翼正在向水平位置復(fù)位，此時(shí)展項(xiàng)中前端由于受到慣性力與氣動(dòng)力的影響，翼展的大部分區(qū)域還處在大變形的狀態(tài)。首先分析0.001 m支架柔性翅的流場(chǎng)情況(見圖13)，0.001 m支架柔性翅在0.05 s時(shí)，盡管此時(shí)翅翼的翼根部已經(jīng)過了極限位置并正在向水平位置復(fù)位，但結(jié)合撲動(dòng)函數(shù)易知，翼根部位在復(fù)位前期，角速度極低，是整個(gè)撲動(dòng)過程中，角速度最接近0的位置。而此時(shí)翅翼的前端與中部受到慣性力與氣動(dòng)力的影響，繼續(xù)向下運(yùn)動(dòng)，到0.06 s時(shí)能夠明顯觀察到翅的變形。從流場(chǎng)速度矢量可以明顯地看出，在該截面中大部分氣流在ZOY平面上的投影都指向了水平方向，只有少部分的氣流在翼根處指向y軸的正方向，而此時(shí)翼根處一方面角速度較低，另一方面距離旋轉(zhuǎn)中心較近，因此實(shí)際翼根處的旋轉(zhuǎn)速度非常低。因此此時(shí)翅翼產(chǎn)生的力，大部分是水平力，如果翅翼是對(duì)稱分布，那么水平方向的力是完全能夠抵消掉，不會(huì)對(duì)飛行過程造成影響。且0.06～0.08 s期間，翅翼持續(xù)變形，速度矢量在ZOY平面的投影大部分都持續(xù)指向水平方向。這種持續(xù)變形的現(xiàn)象反映到升力系數(shù)變化圖線上，即為波動(dòng)時(shí)域的在周期內(nèi)的時(shí)間占比較長。

類似地，分析0.002 m支架柔性翅的波動(dòng)時(shí)域?qū)?yīng)的流場(chǎng)速度矢量圖(見圖14)。產(chǎn)生波動(dòng)時(shí)域的原因與0.001 m支架柔性翅的原因相同，即柔性翅翅翼中前端的持續(xù)運(yùn)動(dòng)導(dǎo)致了一定時(shí)間內(nèi)，翅翼主要產(chǎn)生了水平力，而升力則基本在0點(diǎn)附近徘徊。但與0.001 m支架柔性翅不同的是，0.002 m支架柔性翅出現(xiàn)變形的時(shí)間更早，在0.05 s附近，翅翼已經(jīng)明顯變形，且恢復(fù)變形的速度也更快，在0.06 s時(shí)，流場(chǎng)速度矢量的已經(jīng)指向了Y軸，或者說流場(chǎng)速度矢量在Y軸的分量在0.06 s時(shí)已經(jīng)不容忽視，翅翼產(chǎn)生的主要力已經(jīng)由水平力轉(zhuǎn)變成了升力。0.10 s時(shí)，翅翼的變形過程、波動(dòng)時(shí)域產(chǎn)生的原因與0.05 s時(shí)均是同理，這里便不再贅述。

圖13 0.001 m支架柔性翅波動(dòng)時(shí)域流場(chǎng)圖

圖14 0.002 m支架柔性翅波動(dòng)時(shí)域流場(chǎng)圖

造成兩者變形過程不同的因素是顯然的，二者流場(chǎng)相同、撲翼函數(shù)相同、翼型相同，只有支架的厚度相差了0.001 m。而支架影響翅翼恢復(fù)形狀的原因則主要是由于彈性變形的回彈力矩產(chǎn)生了巨大差異。下面根據(jù)彈塑性力學(xué)理論，作以簡要分析。

由于該支架的厚度只有0.001～0.002 m，其厚度遠(yuǎn)遠(yuǎn)小于整體支架的長度，因此可以作為薄板處理，根據(jù)彈性力學(xué)的薄板理論[10]，薄板的內(nèi)力矩計(jì)算公式如下：

(7)

式中：Mx、My、Mz是作用在中面上每單位長度的彎矩與扭矩；h為薄板厚度；ω為薄板撓度；D為板件的抗彎剛度；κx、κy為變形后中面沿x、y方向的曲率；κxy為扭曲率；v為泊松比且有

(8)

在式(7)中，D為板件的抗彎剛度。且有：

(9)

其中：h為板件的厚度；E為彈性模量；v為泊松比。根據(jù)以上公式可知，若0.001 m支架柔性翅與0.002 m支架柔性翅發(fā)生同樣程度的變形，那么0.002 m支架柔性翅產(chǎn)生的內(nèi)力矩將是0.001 m柔性翅內(nèi)力矩的8倍。

再有根據(jù)公式：

(10)

(11)

式中：β為角加速度；I為轉(zhuǎn)動(dòng)慣量，顯然0.002 m支架柔性翅的轉(zhuǎn)動(dòng)慣量在同樣的變形情況下是0.001 m支架柔性翅轉(zhuǎn)動(dòng)慣量的2倍。綜合考慮轉(zhuǎn)動(dòng)慣量與內(nèi)力矩，可知0.002 m支架柔性翅在同樣的變形情況下的其恢復(fù)變形的角加速度是0.001 m支架柔性翅的4倍。因此也就出現(xiàn)了如圖12所呈現(xiàn)的變形與恢復(fù)變形的情況，0.001 m支架柔性翅的翅翼變形恢復(fù)時(shí)間遠(yuǎn)長于0.002 m柔性翅翅翼。

綜上所述，波動(dòng)時(shí)域產(chǎn)生的原因是：柔性翅在極限位置附近發(fā)生了變形，導(dǎo)致翅翼的中上段延遲到達(dá)原本的既定位置，在該過程中，翅翼產(chǎn)生的主要力為水平力，而在雙翅翼的情況下，水平力相互抵消。升力部分則主要由翼型與翼根部分對(duì)氣流的作用決定，且其作用力很小。而柔性翅支架的形狀與材料，決定了柔性翅的變形恢復(fù)速度，也就決定了波動(dòng)時(shí)域的長短。

4 試驗(yàn)

本文撲翼研究始終以鴿子翼為仿生對(duì)象，展開了一系列的數(shù)值計(jì)算，得出了相關(guān)剛性翅與柔性翅的相關(guān)結(jié)論。在翅翼形變的基礎(chǔ)上提出了波動(dòng)時(shí)域理論。為了驗(yàn)證波動(dòng)時(shí)域理論在鳥類飛行過程中是否具有實(shí)際意義，設(shè)計(jì)了高速攝影實(shí)驗(yàn)以拍攝鴿子的撲動(dòng)過程。

由于鴿子的撲動(dòng)過程并非一個(gè)固定的函數(shù)，為了方便對(duì)比，將鴿子的撲動(dòng)過程以周期T來表征，結(jié)果如圖15所示。

圖15 鴿子的實(shí)際撲動(dòng)與柔性翅仿真過程

從圖13可以看出，在下?lián)潆A段的變形過程中，鴿子翅翼的變形過程與數(shù)值計(jì)算的結(jié)果十分相似，但在翅翼上撲階段，即0.4t～0.8t之間，鴿子翅翼的變形過程與數(shù)值計(jì)算結(jié)果不同。這是由于鴿子在翅翼的上撲過程中，中前段翅翼一直盡量保持在豎直狀態(tài)，也就是說，鴿子翅翼在上撲過程中，鴿子盡量放緩了翅翼恢復(fù)成一個(gè)平面的過程，使得翅翼在上撲過程中盡量不產(chǎn)生負(fù)升力。如果聯(lián)系0.002 m支架柔性翅升力系數(shù)曲線圖，可以認(rèn)為，鴿子在上撲過程中利用身體肌肉與骨骼結(jié)構(gòu)延長了第一段波動(dòng)時(shí)域在整個(gè)周期的時(shí)間占比，縮短了負(fù)時(shí)域的長度。

因此，實(shí)際鳥類在飛行中不但運(yùn)用了柔性翅的波動(dòng)時(shí)域特征，而且其自身的骨骼結(jié)構(gòu)提升了波動(dòng)時(shí)域的利用率，更利于撲翼飛行。

5 結(jié)論

1) 柔性翅的升力特性優(yōu)于剛性翅，柔性翅產(chǎn)生大升力的主要原因在于波動(dòng)時(shí)域的存在，波動(dòng)時(shí)域延遲了負(fù)升力的到來時(shí)間，縮短了負(fù)升力在周期內(nèi)的時(shí)間占比。

2) 柔性翅的支架對(duì)柔性翅的波動(dòng)時(shí)域與峰值特性影響較大，主要由其結(jié)構(gòu)、材料特性決定。

3) 波動(dòng)時(shí)域在自然界鳥類飛行起著重要作用，并且由于鳥類的自身骨骼結(jié)構(gòu)影響，鳥類對(duì)波動(dòng)時(shí)域的利用率更高。