劉旭博，王芬芬，于紀言

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

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

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

1.1 撲動簡介

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

圖1 撲翼運動的4種運動形式示意圖

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

1.2 剛性翅模型

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

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

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

圖2 鴿翼與剛性翅模型

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

1.3 柔性翅模型

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

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

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

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

圖4 計算域

1.4 撲動函數(shù)

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

(1)

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

(2)

2 計算前處理

2.1 計算方法

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

2.2 計算模型

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

1) 流體控制方程：

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

(3)

運動方程：

(4)

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

2) 固體控制方程：

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

ρd=▽·σ+F

(5)

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

3) 流固耦合方程：

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

(6)

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

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

2.3 計算平臺的搭建

ANSYS進行流固耦合計算時，在workbench中搭建如圖5所示的平臺。

圖5 計算平臺示意圖

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

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

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

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

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

3 計算過程

3.1 翅翼撲動升、阻力變化特點

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

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

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

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

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

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

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

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

3.2 波動時域特征

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

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

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

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

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

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

圖13 0.001 m支架柔性翅波動時域流場圖

圖14 0.002 m支架柔性翅波動時域流場圖

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

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

(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)動慣量，顯然0.002 m支架柔性翅的轉(zhuǎn)動慣量在同樣的變形情況下是0.001 m支架柔性翅轉(zhuǎn)動慣量的2倍。綜合考慮轉(zhuǎn)動慣量與內(nèi)力矩，可知0.002 m支架柔性翅在同樣的變形情況下的其恢復變形的角加速度是0.001 m支架柔性翅的4倍。因此也就出現(xiàn)了如圖12所呈現(xiàn)的變形與恢復變形的情況，0.001 m支架柔性翅的翅翼變形恢復時間遠長于0.002 m柔性翅翅翼。

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

4 試驗

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

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

圖15 鴿子的實際撲動與柔性翅仿真過程

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

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

5 結(jié)論

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

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

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