陳利麗，宋筆鋒，宋文萍，楊文青

柔性撲翼氣動(dòng)結(jié)構(gòu)耦合特性數(shù)值研究

陳利麗1，宋筆鋒2，宋文萍2，楊文青2

(1.中國(guó)航空工業(yè)集團(tuán)第一飛機(jī)設(shè)計(jì)研究院，西安710089;
2.西北工業(yè)大學(xué)翼型葉柵空氣動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，西安710072)

微型撲翼體積小、重量輕，其柔性變形對(duì)氣動(dòng)特性有顯著的影響。通過(guò)求解雷諾平均N-S方程(Reynolds－Averaged Navier-Stokes，RANS)和結(jié)構(gòu)動(dòng)力學(xué)方程，對(duì)微型柔性撲翼飛行器的氣動(dòng)結(jié)構(gòu)耦合特性進(jìn)行了數(shù)值模擬研究。針對(duì)微型撲翼的大幅運(yùn)動(dòng)，發(fā)展了適用于撲翼的氣動(dòng)結(jié)構(gòu)耦合數(shù)值計(jì)算方法，研究了微型撲翼的氣動(dòng)結(jié)構(gòu)耦合特性。通過(guò)求解雷諾平均Navier-Stokes(RANS)方程得到微型撲翼的非定常氣動(dòng)特性;利用哈密頓原理(Hamilton Principle)推導(dǎo)了撲翼的結(jié)構(gòu)動(dòng)力學(xué)方程，采用結(jié)構(gòu)有限元方法對(duì)該動(dòng)力學(xué)方程進(jìn)行離散并求解，得到撲翼的動(dòng)態(tài)結(jié)構(gòu)特性;采用松耦合方法進(jìn)行迭代。計(jì)算結(jié)果與風(fēng)洞實(shí)驗(yàn)結(jié)果相比吻合良好，驗(yàn)證了所發(fā)展方法的有效性。在此基礎(chǔ)上研究了慣性力和關(guān)鍵運(yùn)動(dòng)參數(shù)對(duì)柔性撲翼氣動(dòng)及結(jié)構(gòu)特性的影響規(guī)律，有助于比較詳細(xì)、全面地了解微型撲翼的氣動(dòng)機(jī)理，為柔性撲翼的設(shè)計(jì)提供了參考依據(jù)。

微型撲翼;氣動(dòng)結(jié)構(gòu)耦合;雷諾平均Navier-Stokes;結(jié)構(gòu)動(dòng)力學(xué);非定常氣動(dòng)特性

0 引言

微型撲翼研究工作者們?cè)陂L(zhǎng)期研究中發(fā)現(xiàn)，柔性結(jié)構(gòu)變形對(duì)撲翼的氣動(dòng)特性有顯著影響［1-7］。因此，在微型撲翼飛行器的結(jié)構(gòu)設(shè)計(jì)中可以考慮根據(jù)氣動(dòng)方面的需要，控制柔性撲翼產(chǎn)生適當(dāng)?shù)慕Y(jié)構(gòu)變形來(lái)改善撲翼的氣動(dòng)特性。例如，在翼梢部分需要相對(duì)柔軟的結(jié)構(gòu)和更高的運(yùn)動(dòng)速度來(lái)提高推力特性;而翼根部分需要相對(duì)堅(jiān)固的結(jié)構(gòu)來(lái)抑制結(jié)構(gòu)變形，從而保證翼型剖面維持良好的升力特性［8］。

影響撲翼柔性結(jié)構(gòu)變形的因素主要有以下三方面:(1)氣動(dòng)載荷(撲翼自身的運(yùn)動(dòng)和飛行狀態(tài)產(chǎn)生);(2)慣性載荷(撲翼自身的運(yùn)動(dòng)和結(jié)構(gòu)重量產(chǎn)生);(3)結(jié)構(gòu)彈性(撲翼結(jié)構(gòu)剛度產(chǎn)生)。其中，氣動(dòng)載荷與結(jié)構(gòu)彈性是相互耦合的，氣動(dòng)載荷與結(jié)構(gòu)變形相互影響。因此，在對(duì)撲翼非定常氣動(dòng)特性進(jìn)行數(shù)值模擬時(shí)，有必要考慮結(jié)構(gòu)彈性的影響。慣性載荷取決于撲翼自身的運(yùn)動(dòng)以及結(jié)構(gòu)質(zhì)量分布。一旦確定以上條件，即可由結(jié)構(gòu)有限元方法(Finite Element Method，F(xiàn)EM)求得撲翼的結(jié)構(gòu)變形。

撲翼機(jī)構(gòu)的運(yùn)動(dòng)方式主要取決于兩個(gè)參數(shù):揮舞角和撲動(dòng)頻率。撲翼的運(yùn)動(dòng)方式影響慣性載荷以及空氣動(dòng)力特性，進(jìn)而影響撲翼的柔性結(jié)構(gòu)變形。

由理論分析可知，對(duì)于做正弦規(guī)律運(yùn)動(dòng)的撲翼，其最大慣性載荷分別與撲動(dòng)頻率的平方和揮舞角呈正比關(guān)系。這說(shuō)明撲翼的柔性結(jié)構(gòu)變形對(duì)撲動(dòng)頻率的變化表現(xiàn)得更加“敏感”，并且增加撲動(dòng)頻率對(duì)機(jī)構(gòu)來(lái)說(shuō)更容易實(shí)現(xiàn)。因此，撲動(dòng)頻率可以用來(lái)作為控制柔性撲翼結(jié)構(gòu)變形的主要參數(shù)。固定撲動(dòng)頻率、增大揮舞角使撲翼獲得更高的瞬時(shí)撲動(dòng)速度，與對(duì)慣性力的影響相比，將會(huì)帶來(lái)更多的氣動(dòng)方面的影響。

21世紀(jì)以來(lái)，隨著實(shí)驗(yàn)技術(shù)和計(jì)算技術(shù)的發(fā)展，國(guó)內(nèi)外研究學(xué)者分別通過(guò)實(shí)驗(yàn)和計(jì)算兩種方法對(duì)柔性撲翼氣動(dòng)結(jié)構(gòu)耦合特性進(jìn)行了一些研究。Heathcote等人［9］通過(guò)實(shí)驗(yàn)方法對(duì)做沉浮運(yùn)動(dòng)的柔性撲翼剛度對(duì)推力特性的影響進(jìn)行了研究，結(jié)果表明，與剛性撲翼相比，柔性撲翼能夠產(chǎn)生更大的推力，并通過(guò)水洞實(shí)驗(yàn)對(duì)半展弦比為3、做沉浮運(yùn)動(dòng)的NACA0012矩形機(jī)翼氣動(dòng)特性進(jìn)行測(cè)量，發(fā)現(xiàn)展向柔性變形能夠顯著提高撲翼的推力。在數(shù)值計(jì)算方面，Wu等人［10］發(fā)展了一種多學(xué)科耦合實(shí)驗(yàn)方法，測(cè)量了六副不同結(jié)構(gòu)(每副機(jī)翼均為板條-薄膜結(jié)構(gòu))的微型撲翼做單自由度撲動(dòng)運(yùn)動(dòng)時(shí)的氣動(dòng)彈性響應(yīng)和氣動(dòng)特性。通過(guò)數(shù)字圖像關(guān)聯(lián)系統(tǒng)測(cè)量結(jié)構(gòu)變形，通過(guò)載荷傳感器記錄作用于撲翼上的氣動(dòng)載荷，利用數(shù)字粒子圖像測(cè)控技術(shù)(DPIV)測(cè)量流場(chǎng)特性。他們的研究結(jié)果表明，結(jié)構(gòu)變形對(duì)柔性撲翼氣動(dòng)特性有顯著影響，通過(guò)適當(dāng)?shù)娜嵝越Y(jié)構(gòu)變形，能夠獲得更好的氣動(dòng)特性。Singh等人［11-13］利用自行開(kāi)發(fā)的實(shí)驗(yàn)裝置和仿生撲動(dòng)機(jī)構(gòu)來(lái)測(cè)量機(jī)翼結(jié)構(gòu)的推力特性，結(jié)果表明，在其實(shí)驗(yàn)狀態(tài)范圍內(nèi)，慣性力對(duì)機(jī)翼結(jié)構(gòu)變形起主導(dǎo)作用。由此說(shuō)明，在相關(guān)氣動(dòng)結(jié)構(gòu)耦合研究中考慮慣性力的影響是非常必要的。Zhu等［14］發(fā)展了一套基于三維邊界元法和二維非線性薄鋼板結(jié)構(gòu)動(dòng)力學(xué)模型方法的氣動(dòng)結(jié)構(gòu)耦合分析求解器，并將數(shù)值模擬結(jié)果與Heathcote［9］做沉浮運(yùn)動(dòng)的柔性撲翼實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比，驗(yàn)證了算法的正確性。

西北工業(yè)大學(xué)針對(duì)微型撲翼氣動(dòng)結(jié)構(gòu)耦合特性開(kāi)展了一些研究工作:楊文青等人［15］發(fā)展了基于非定常雷諾平均N-S方程和靜態(tài)結(jié)構(gòu)有限元方法的氣動(dòng)結(jié)構(gòu)耦合求解器，初步研究了考慮柔性變形的微型撲翼氣動(dòng)特性。但是該方法沒(méi)有考慮到慣性力和運(yùn)動(dòng)加速度的影響，而對(duì)于碳桿-薄膜結(jié)構(gòu)撲翼飛行器，慣性力和運(yùn)動(dòng)加速度的影響是不應(yīng)忽視的。

本文發(fā)展了一套適用于撲翼氣動(dòng)結(jié)構(gòu)耦合特性的數(shù)值模擬方法:利用Hamilton原理，對(duì)動(dòng)能和應(yīng)變能進(jìn)行變分，得到撲翼動(dòng)力學(xué)方程，采用結(jié)構(gòu)有限元方法對(duì)該方程進(jìn)行離散并求解，得到微型撲翼的動(dòng)態(tài)結(jié)構(gòu)特性;通過(guò)求解三維非定常雷諾平均Navier-Stokes方程獲得精確的非定常氣動(dòng)特性［16］;采用松耦合方法進(jìn)行迭代;采用徑向基函數(shù)(radial basis function，RBF)方法進(jìn)行氣動(dòng)網(wǎng)格與結(jié)構(gòu)網(wǎng)格之間的數(shù)據(jù)傳遞［17］。采用上述方法對(duì)柔性撲翼氣動(dòng)特性進(jìn)行數(shù)值模擬，計(jì)算結(jié)果與風(fēng)洞實(shí)驗(yàn)結(jié)果吻合良好，驗(yàn)證了該方法的有效性。在此基礎(chǔ)上進(jìn)一步對(duì)柔性撲翼的氣動(dòng)結(jié)構(gòu)耦合特性展開(kāi)詳細(xì)的數(shù)值模擬研究，并對(duì)慣性載荷、氣動(dòng)載荷及其影響下的撲翼柔性結(jié)構(gòu)變形進(jìn)行定量分析，分別研究?jī)蓚€(gè)重要運(yùn)動(dòng)參數(shù)(撲動(dòng)頻率、揮舞角)對(duì)柔性撲翼氣動(dòng)結(jié)構(gòu)特性的影響。

1 流動(dòng)控制方程

積分形式的三維非定?？蓧嚎s雷諾平均N-S方程可寫(xiě)成如下形式:

其中，

其中，q=(u，v，w)T為流體速度矢量，qb為柔性撲翼表面邊界的運(yùn)動(dòng)速度矢量，Ix，Iy，Iz分別為直角坐標(biāo)系單位坐標(biāo)向量。τ和β分別為粘性應(yīng)力和傳熱項(xiàng)。

采用有限體積法求解上述方程組［16，18］。時(shí)間推進(jìn)采用基于LU-SGS迭代的全隱式格式，空間離散采用Jameson中心格式，湍流模型采用k-ω SST模型。繞撲翼的網(wǎng)格為C-O拓?fù)浣Y(jié)構(gòu)的結(jié)構(gòu)化網(wǎng)格，采用基于廣義無(wú)限插值理論的代數(shù)網(wǎng)格方法自動(dòng)生成［19］，非常適合撲翼的撲動(dòng)運(yùn)動(dòng)，完全適用于柔性撲翼的氣動(dòng)結(jié)構(gòu)耦合計(jì)算，網(wǎng)格拓?fù)浣Y(jié)構(gòu)如圖1所示。

圖1 撲翼C-O結(jié)構(gòu)網(wǎng)格Fig.1Flapping wing C-O type grid

2 結(jié)構(gòu)動(dòng)力學(xué)方程

利用哈密頓原理，撲翼的結(jié)構(gòu)運(yùn)動(dòng)方程可以表示為:其中，t1，t2表示任意時(shí)刻，δU、δW、δT分別代表應(yīng)變能項(xiàng)、動(dòng)能項(xiàng)和外力所做的功。

柔性撲翼結(jié)構(gòu)如圖2所示，由碳纖維骨架和聚醚薄膜蒙皮結(jié)構(gòu)組成，考慮到聚醚薄膜蒙皮僅能承受張力，主要變形為碳纖維骨架所產(chǎn)生。碳纖維骨架作為主要的受力結(jié)構(gòu)，其承受大部分彎矩和所有的剪力。

圖2 柔性撲翼模型Fig.2Flexible flapping wing

撲翼模型的結(jié)構(gòu)幾何參數(shù)為:前梁2.0mm，斜梁1.5mm，翼肋1.2mm。

對(duì)碳纖維桿的材料特性進(jìn)行測(cè)量，結(jié)果如圖3所示。

圖3 碳纖維桿載荷-應(yīng)變曲線Fig.3Load-strain curve of carbon fiber

由圖可以看出，試驗(yàn)結(jié)果的重合性良好。經(jīng)過(guò)對(duì)試驗(yàn)測(cè)量結(jié)果的整理，得到了碳纖維桿試件的彈性模量，見(jiàn)表1。

表1 碳纖維桿復(fù)合材料彈性模量Table 1Elastic modulus of carbon fiber

采用基于鐵木辛哥梁結(jié)構(gòu)理論的梁?jiǎn)卧獊?lái)模擬碳纖維桿結(jié)構(gòu)。該單元模型適合于分析從細(xì)長(zhǎng)到中等粗短的梁結(jié)構(gòu)，并考慮了碳纖維桿剪切變形的影響。

首先，對(duì)應(yīng)變能項(xiàng)進(jìn)行變分。利用材料的本構(gòu)關(guān)系，可以得到應(yīng)變能的變分形式:

其中，［Q］為本構(gòu)關(guān)系矩陣。

對(duì)動(dòng)能項(xiàng)進(jìn)行變分。圖4為做俯仰揮舞運(yùn)動(dòng)的撲翼的參考坐標(biāo)系。其中慣性坐標(biāo)系Xi，Yi，Zi固連于翼根前緣位置處。揮舞角ψ表示揮舞坐標(biāo)系Xf、Yf、Zf繞Xp軸旋轉(zhuǎn)角度。俯仰角α為俯仰坐標(biāo)系Xp、Yp、Zp繞Xp軸旋轉(zhuǎn)角度。偏轉(zhuǎn)角θ為梁(前梁、翼肋、斜梁)局部坐標(biāo)系繞Zf軸旋轉(zhuǎn)角度。坐標(biāo)系之間的轉(zhuǎn)換關(guān)系可以表示為如下形式:

其中，運(yùn)動(dòng)參數(shù)α、ψ均為時(shí)間t的函數(shù)。

圖4 撲翼坐標(biāo)系Fig.4Flapping wing coordinate system

固連于梁上任一點(diǎn)(ξ1，0)在局部坐標(biāo)系上的坐標(biāo)可以表示為:

其中u、w為ξ、t的函數(shù)，分別表示軸向變形位移和撓度變形位移。

對(duì)位置坐標(biāo)進(jìn)行求導(dǎo)，可以得到任一點(diǎn)的速度表達(dá)式:

其中ρ為材料的密度。

將應(yīng)變能、動(dòng)能的表達(dá)式離散到結(jié)構(gòu)單元內(nèi)，對(duì)動(dòng)能項(xiàng)進(jìn)行變分，從中提取出質(zhì)量矩陣、動(dòng)態(tài)剛度矩陣和單元體積力項(xiàng)，分別為:

其中N1、N2分別表示u和w的形狀函數(shù)。

利用單元節(jié)點(diǎn)之間的聯(lián)系，對(duì)單元質(zhì)量矩陣、靜態(tài)/動(dòng)態(tài)剛度矩陣和單元體積力矢量進(jìn)行組裝，可得到整體質(zhì)量、剛度矩陣和整體體積力矢量。由此，可以得到離散形式的結(jié)構(gòu)運(yùn)動(dòng)方程:

其中F(t)由氣動(dòng)力以及運(yùn)動(dòng)加速度產(chǎn)生的慣性力項(xiàng)組成。

對(duì)式(13)進(jìn)行離散時(shí)間離散并求解［20］，可以得到t時(shí)刻的位移、速度以及加速度。

3CFD/CSD耦合方法

本文的CFD/CSD耦合采用松耦合策略，具體過(guò)程如下:

(1)通過(guò)求解非定常Navier-Stokes方程得到非定常氣動(dòng)力;

(2)把氣動(dòng)網(wǎng)格的氣動(dòng)力轉(zhuǎn)換到結(jié)構(gòu)網(wǎng)格上，對(duì)撲翼結(jié)構(gòu)動(dòng)力學(xué)方程進(jìn)行求解得到相應(yīng)時(shí)刻撲翼的結(jié)構(gòu)變形位移;

(3)將結(jié)構(gòu)網(wǎng)格的變形位移轉(zhuǎn)換到氣動(dòng)網(wǎng)格上，更新柔性撲翼外形;

(4)檢查結(jié)構(gòu)變形是否收斂，如果收斂，則計(jì)算結(jié)束;

(5)更新結(jié)構(gòu)變形后回到第1步。

采用該耦合方法，可以在3-4個(gè)周期內(nèi)達(dá)到收斂。

4 算例驗(yàn)證與分析

對(duì)柔性撲翼模型進(jìn)行了數(shù)值模擬，并將計(jì)算結(jié)果與風(fēng)洞實(shí)驗(yàn)結(jié)果進(jìn)行比較。

4.1 撲翼運(yùn)動(dòng)模式

撲翼機(jī)構(gòu)的運(yùn)動(dòng)為揮舞運(yùn)動(dòng)，如圖5所示。

圖5 撲翼運(yùn)動(dòng)示意圖Fig.5Flapping wing movement

撲翼運(yùn)動(dòng)規(guī)律表示如下:

其中，ψ的取值均以上反為正，Ψ1為揮舞運(yùn)動(dòng)的振幅。

4.2 算法驗(yàn)證

4.2.1 網(wǎng)格和時(shí)間精度驗(yàn)證

首先對(duì)網(wǎng)格收斂性和時(shí)間精度進(jìn)行了驗(yàn)證。分別計(jì)算了在三種不同網(wǎng)格密度和時(shí)間步數(shù)下的周期性升力系數(shù)，結(jié)果如圖6和圖7所示。

圖6 不同網(wǎng)格密度計(jì)算結(jié)果比較Fig.6Comparison of result with different mesh sizes

圖7 不同時(shí)間步長(zhǎng)計(jì)算結(jié)果比較Fig.7Comparison of result with different time steps

從上圖可以看出，當(dāng)網(wǎng)格密度達(dá)到152×56×56，時(shí)間迭代步數(shù)達(dá)到80步時(shí)，網(wǎng)格精度和時(shí)間精度均趨于收斂。

4.2.2 算例驗(yàn)證

將撲翼模型放入西北工業(yè)大學(xué)低湍流度風(fēng)洞中進(jìn)行非定常氣動(dòng)力的測(cè)量。風(fēng)洞實(shí)驗(yàn)環(huán)境如圖8所示。計(jì)算和實(shí)驗(yàn)狀態(tài)為:風(fēng)速8m/s，撲動(dòng)頻率6Hz，撲動(dòng)幅度±35°。柔性機(jī)翼結(jié)構(gòu)參數(shù)見(jiàn)表1。

在數(shù)值計(jì)算中，由于CFD計(jì)算對(duì)網(wǎng)格拓?fù)浣Y(jié)構(gòu)有限制，需要對(duì)實(shí)際柔性撲翼模型做適當(dāng)簡(jiǎn)化。根據(jù)柔性翼翼肋直徑，選取最大厚度為1.2%的薄翼型作為撲翼的翼剖面，如圖9所示。用于CFD計(jì)算的網(wǎng)格單元數(shù)為152×56×56。撲翼結(jié)構(gòu)計(jì)算模型如圖10所示。

圖10 結(jié)構(gòu)計(jì)算模型Fig.10Computational structural model

計(jì)算結(jié)果與實(shí)驗(yàn)測(cè)量結(jié)果的比較如圖11所示。

圖11中，實(shí)線表示CFD/CSD耦合計(jì)算結(jié)果，虛線表示不考慮柔性結(jié)構(gòu)變形的計(jì)算結(jié)果?？梢钥闯?，不考慮柔性結(jié)構(gòu)變形計(jì)算出的平均升力系數(shù)隨預(yù)安裝角變化基本呈線性，而考慮柔性結(jié)構(gòu)變形的計(jì)算結(jié)果呈現(xiàn)出與實(shí)驗(yàn)結(jié)果一致的非線性特性;平均推力系數(shù)的計(jì)算值均小于實(shí)驗(yàn)值，但趨勢(shì)與實(shí)驗(yàn)值一致，且考慮柔性結(jié)構(gòu)變形所得結(jié)果與實(shí)驗(yàn)值更為接近。對(duì)于計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果存在的差異，分析主要有以下原因:(1)實(shí)驗(yàn)中支架、洞壁等因素會(huì)對(duì)實(shí)驗(yàn)結(jié)果產(chǎn)生干擾;(2)計(jì)算模型的簡(jiǎn)化會(huì)對(duì)計(jì)算結(jié)果產(chǎn)生一定影響。從結(jié)果可以看出，對(duì)于研究的柔性撲翼模型，結(jié)構(gòu)變形對(duì)氣動(dòng)特性有很大影響，因此在數(shù)值模擬中考慮柔性結(jié)構(gòu)變形是十分必要的。

4.3 撲動(dòng)頻率對(duì)撲翼氣動(dòng)結(jié)構(gòu)耦合特性的影響

本節(jié)研究了撲動(dòng)頻率對(duì)柔性撲翼氣動(dòng)及結(jié)構(gòu)特性的影響，計(jì)算撲動(dòng)頻率從3Hz～12Hz范圍內(nèi)的變化情況，每隔3Hz計(jì)算一個(gè)結(jié)果。其余參數(shù)保持不變，分別為:風(fēng)速=6m/s，揮舞角ψ1=35°，預(yù)安裝角α0= 0°。撲翼的周期結(jié)構(gòu)變形計(jì)算結(jié)果如圖12所示。

翼梢前緣位置處柔性變形位移的大小反映了柔性撲翼結(jié)構(gòu)的彎曲剛度特性，但是從結(jié)構(gòu)變形位移云圖上不能對(duì)柔性結(jié)構(gòu)變形與氣動(dòng)特性之間的關(guān)系進(jìn)行進(jìn)一步直觀定量的分析。為此，根據(jù)撲翼柔性結(jié)構(gòu)變形的特點(diǎn)，將柔性變形分解為兩部分:彎曲變形(沿展向)和扭轉(zhuǎn)變形(沿弦向)，分別用翼梢前緣彎曲角度αbend和翼梢后緣扭轉(zhuǎn)角度αtwist表示。柔性撲翼產(chǎn)生的αbend和αtwist分別隨撲動(dòng)揮舞角的變化形成一個(gè)閉合的相位環(huán)，變化過(guò)程如圖13所示。

圖11 計(jì)算結(jié)果與實(shí)驗(yàn)值比較Fig.11Comparison between computational and experimental result

圖12 撲翼在3Hz、9Hz、12Hz的周期結(jié)構(gòu)變形Fig.12Periodic structural deformation in 3Hz、9Hz、12Hz

圖13 撲翼在不同頻率下的周期結(jié)構(gòu)變形相位圖Fig.13Periodic deformation phase in different frequencies

從圖12和圖13可以看出，柔性撲翼在撲動(dòng)過(guò)程中均會(huì)產(chǎn)生柔性彎曲變形和扭轉(zhuǎn)變形。在撲動(dòng)頻率為3Hz時(shí)，當(dāng)撲翼從最高點(diǎn)往下?lián)鋾r(shí)，在最高點(diǎn)處，撲翼產(chǎn)生了上翹的柔性彎曲變形;從最高點(diǎn)撲動(dòng)到平衡位置過(guò)程中，撲翼的彎曲變形幅度逐漸增大，在平衡位置附近，撲翼變形達(dá)到最大;從平衡位置撲動(dòng)到最低點(diǎn)過(guò)程中，撲翼上翹變形逐漸反彈，在最低位置附近，撲翼變形達(dá)到最小;從最低點(diǎn)上撲到平衡位置過(guò)程中，撲翼產(chǎn)生向下的彎曲變形，變形量先增大后減小，在7/8T時(shí)刻達(dá)到最小，從7/8T時(shí)刻撲動(dòng)到最高點(diǎn)過(guò)程中，撲翼產(chǎn)生上翹的變形。

撲動(dòng)頻率增大到9Hz時(shí)，撲翼一周期各時(shí)刻柔性變形的趨勢(shì)與3Hz時(shí)大致相同，主要區(qū)別有兩點(diǎn): (1)彎曲變形的幅度顯著增大;(2)撲翼最大彎曲變形和最大扭轉(zhuǎn)變形發(fā)生時(shí)刻提前，更靠近下?lián)浜蜕蠐涑跏嘉恢谩?/p>

撲動(dòng)頻率繼續(xù)增大到12Hz時(shí)，撲翼的變形幅度持續(xù)增大，最大彎曲變形和最大扭轉(zhuǎn)變形發(fā)生在下?lián)浜蜕蠐涞某跏紩r(shí)刻。

進(jìn)一步研究撲翼氣動(dòng)特性隨撲動(dòng)頻率的變化趨勢(shì)以及柔性結(jié)構(gòu)變形對(duì)撲翼氣動(dòng)特性的影響。周期性升力系數(shù)和推力系數(shù)變化曲線如圖14所示。

圖14 不同頻率下周期性升力系數(shù)和推力系數(shù)的變化曲線Fig.14Periodic lift and thrust coefficient in different frequencies

從圖14中可以觀察到，對(duì)于剛性撲翼，隨著撲動(dòng)頻率的增加，升力系數(shù)和推力系數(shù)曲線峰值均顯著增加，并且增加幅度越來(lái)越大。對(duì)于柔性撲翼，在撲動(dòng)頻率較低(低于6Hz)時(shí)，與剛性撲翼相比，柔性撲翼周期性升力系數(shù)曲線峰值有所降低;當(dāng)撲動(dòng)頻率較高(高于9Hz)時(shí)，與剛性撲翼相比，柔性撲翼周期性升力系數(shù)曲線峰值有所增長(zhǎng)，并且隨著撲動(dòng)頻率的增大，增長(zhǎng)幅度越來(lái)越大。

平均升力系數(shù)和平均推力系數(shù)隨撲動(dòng)頻率的變化如圖15所示。

可以觀察到，對(duì)于剛性撲翼，平均升力系數(shù)隨撲動(dòng)頻率的增加而增大，平均推力系數(shù)隨撲動(dòng)頻率的增加略微有所增大;對(duì)于柔性撲翼，平均升力系數(shù)與剛性撲翼相比有所降低，并且隨撲動(dòng)頻率的增加呈減小趨勢(shì)，平均推力系數(shù)與剛性撲翼相比有顯著增大，頻率越大，增大趨勢(shì)越明顯。

由以上結(jié)果分析可知，由于結(jié)構(gòu)變形的影響，柔性撲翼氣動(dòng)特性隨撲動(dòng)頻率變化的趨勢(shì)與剛性撲動(dòng)情況是不同的，并且隨著撲動(dòng)頻率的增大，柔性變形對(duì)氣動(dòng)特性的影響變得越來(lái)越顯著。采用柔性結(jié)構(gòu)能夠顯著提高推力特性，但對(duì)升力特性則會(huì)產(chǎn)生稍許不利影響;增大撲動(dòng)頻率有利于柔性撲翼推力特性的產(chǎn)生，但對(duì)柔性撲翼升力特性則有不利影響。

圖15 平均升力系數(shù)和推力系數(shù)隨著撲動(dòng)頻率的變化曲線Fig.15Variation of average lift and thrust coefficient with frequencies

4.4 揮舞角對(duì)撲翼氣動(dòng)結(jié)構(gòu)耦合特性的影響

本節(jié)研究揮舞角對(duì)柔性撲翼氣動(dòng)及結(jié)構(gòu)特性的影響，計(jì)算了揮舞角從15°～35°范圍內(nèi)變化的情況，每隔5°計(jì)算一個(gè)結(jié)果。其余參數(shù)保持相同，分別為:風(fēng)速=6m/s，撲動(dòng)頻率=6Hz，預(yù)安裝角α0=0°。

圖16為撲翼在揮舞角為15°、25°和35°時(shí)的周期結(jié)構(gòu)變形。柔性撲翼產(chǎn)生的αbend和αtwist隨撲動(dòng)過(guò)程中揮舞角變化過(guò)程如圖17所示。

從圖16和圖17中可以看出，與隨撲動(dòng)頻率變化趨勢(shì)有所不同，柔性撲翼的彎曲變形和扭轉(zhuǎn)變形幅度隨揮舞角的增大呈線性增長(zhǎng)，而最大彎曲變形和最大扭轉(zhuǎn)變形發(fā)生時(shí)刻則基本不發(fā)生變化。由此可知，揮舞角變化只改變結(jié)構(gòu)變形的峰值，不改變?nèi)嵝宰冃蔚南辔弧?/p>

為研究撲翼氣動(dòng)特性隨揮舞角的變化趨勢(shì)，分別計(jì)算了不同揮舞角下剛性和柔性撲翼的氣動(dòng)特性。周期性升力系數(shù)和推力系數(shù)變化曲線如圖18所示。

從圖18中可以看到，對(duì)于剛性撲翼，隨著揮舞角的增加，升力系數(shù)和推力系數(shù)曲線峰值均有所增加。柔性撲翼升力系數(shù)曲線峰值與剛性撲翼相比略微有所降低，但推力系數(shù)曲線峰值則有明顯升高，并且隨著揮舞角的增大，升高幅度越來(lái)越大。

圖16撲翼在15°、25°、35°的周期結(jié)構(gòu)變形Fig.16Periodic structural deformation in15°、25°、35°

圖17 撲翼在不同揮舞角下的周期結(jié)構(gòu)變形相位圖Fig.17Periodic deformation phase in different flapping angles

平均升力系數(shù)和平均推力系數(shù)隨揮舞角的變化如圖19所示。

由圖19中可以觀察到，對(duì)于剛性撲翼，隨著揮舞角的增加，平均升力系數(shù)基本保持不變，平均推力系數(shù)有略微增長(zhǎng)的趨勢(shì);對(duì)于柔性撲翼，平均升力系數(shù)隨揮舞角的增加基本不發(fā)生變化，但與剛性撲翼相比有所降低，平均推力系數(shù)與剛性撲翼相比顯著增大，隨揮舞角的增加一直增大。

由此可知，增大揮舞角有利于柔性撲翼推力特性的產(chǎn)生，對(duì)柔性撲翼升力特性則幾乎不產(chǎn)生影響。

圖18 不同揮舞角下周期性升力系數(shù)和推力系數(shù)的變化曲線Fig.18Periodic lift and thrust coefficient in different flapping angles

圖19 平均升力系數(shù)和推力系數(shù)隨著揮舞角的變化曲線Fig.19Variation of average lift and thrust coefficient with flapping angle

5 結(jié)論

本文從結(jié)構(gòu)變形、氣動(dòng)特性兩個(gè)方面研究了兩個(gè)關(guān)鍵運(yùn)動(dòng)參數(shù)(撲動(dòng)頻率、揮舞角)對(duì)柔性撲翼氣動(dòng)結(jié)構(gòu)耦合特性的影響、結(jié)構(gòu)變形對(duì)柔性撲翼氣動(dòng)特性的影響以及柔性撲翼隨運(yùn)動(dòng)參數(shù)的變化趨勢(shì)，得到了柔性撲翼氣動(dòng)特性方面的一些基本規(guī)律，研究結(jié)果在各節(jié)中均有小結(jié)。從研究結(jié)果可以得出，柔性撲翼產(chǎn)生的結(jié)構(gòu)變形對(duì)氣動(dòng)特性有顯著影響，對(duì)柔性撲翼氣動(dòng)特性進(jìn)行數(shù)值模擬時(shí)，考慮結(jié)構(gòu)變形的影響是非常必要的。與揮舞角相比，撲動(dòng)頻率對(duì)柔性撲翼的氣動(dòng)結(jié)構(gòu)特性影響表現(xiàn)的更為明顯。

在氣動(dòng)特性方面，柔性撲翼結(jié)構(gòu)對(duì)推力特性有顯著地提升作用，并且隨著撲動(dòng)頻率和揮舞角的增大這種作用越來(lái)越明顯，而對(duì)升力特性卻有不利的影響。

［1］Zeng R，Ang H S，Mei Y，et al.Flexibility of flapping wing and its effect on aerodynamic characteristic［J］.Chinese Journal of Computational Mechanics，2005，22(6):750-754.(in Chinese)曾銳，昂海松，梅源，等.撲翼柔性及其對(duì)氣動(dòng)特性的影響［J］.計(jì)算力學(xué)學(xué)報(bào)，2005，22(6):750-754.

［2］Unger R，Haupt M C，Horst P.Structural design and aeroelastic analysis of an oscillating airfoil for flapping wing propulsion［C］//In 46th AIAA Aerospace Sciences Meeting and Exhibit，Reno，Nevada，January 2008.AIAA 2008-306.

［3］Liani E，Guo S，Allegri G.Aerodynamics and aeroelasticity of flexible flapping wings［R］.2007.

［4］Kim D K，Lee J S，Lee J Y，et al.An aeroelastic analysis of a flexible flapping wing using modified strip theory［J］.Active and Passive Smart Structures and Integrated Systems，2008，6928:1-8.

［5］Hamamoto M，Ohta Y，Hara K，et al.Application of fluid-structure in-teraction analysis to flapping flight of insects with deformable wings［J］.Advanced Robotics，2007，21(1-2):1-21.

［6］Liani E，Guo S，Allegri G.Aeroelastic effect on flapping wing perfor-mance［C］//In 48th AIAA/ASME/ASCE/AHS/ASC Structures，Structural Dynamics，and Materials Conference，Honololu，Hawaii，April 2007.AIAA 2007-2412.

［7］Ennos A R.The importance of torsion in the design of inset wings［J］.Journal of Experimental Biology，1988，140:137-160.

［8］Pin Wu.Experimental characterization，design，analysis and optimization of flexible flapping wings for micro air vehicles［D］.［PhD Dissertation］.University of Florida，2010.

［9］Heathcote S，Wang Z，Gursul I.Effect of spanwise flexibility on flapping wing propulsion［J］.Journal of Fluids and Structures，2008，24(2):183-199.

［10］Wu P，Ifju P，Stanford B，et al.A multidisciplinary experimental study of flapping wing aeroelasticity in thrust production［C］//In 50th AIAA/ASME/ASCE/AHS/ASC Structures，Structural Dynamics，and Materials Conference，Palm Springs，CA，May 2009.AIAA 2009-2413.

［11］Singh B.Dynamics and aeroelasticity of hover capable flapping wings:experiments and analysis［D］.［PhD Dissertation］.Department of Aerospace Engineering，University of Maryland，College Park，Maryland，2006.

［12］Singh B，Chopra I.Insect-based flapping wings for micro hovering air ve-hicles:experimental investigationsp［C］//American Helicopter Society International Specialists Meeting on Unmanned Rotorcraft，Arizona，January 2004.

［13］Singh B，Chopra I.Dynamics of insect-based flapping wings:loads and validation［C］//In 47t/AIAA/ASME/ASCE/AHS/ASC Structures，Structural Dynamics，and Materials Conference，Newport，Rhode Island，May 2006.AIAA 2006-1663.

［14］Zhu Q.Numerical simulation of a flapping foil with chordwise or spanwise flexibility［J］.AIAA Journal，2007，45(10):2448-2457.

［15］Yang W Q，Song B F，Song W P.Fluid-structure coupling research for micro flapping-wing［C］//the 27th International Congress of Aeronautic Sciences，Nice，F(xiàn)rance.Sep，19-24，2010.

［16］Xie H，Song W P，Song B F.Numerical solution of navier-stokes equations for flow over a flapping wing［J］.Journal of Northwestern Polytechnical University，2008，26(1):104-109.(in Chinese)謝輝，宋文萍，宋筆鋒.微型撲翼繞流的N-S方程數(shù)值模擬［J］.西北工業(yè)大學(xué)學(xué)報(bào)，2008，26(1):104-109.

［17］Bechert A，Wendland H.Multivariate interpolation for fluid-structure-interaction problems using radial basis functions［J］.Aerospace Science and Technology，2001，5(2):125-134.

［18］Song W P，Yang Y，Qiao Z D，et al.Unsteady N-S calculations using a dual-time stepping method for airfoils oscillating at high angle of attack［J］.Acta Aerodynamica Sinica，1999，17(4):466-471.(in Chinese)宋文萍，楊永，喬志德，等.雙時(shí)間法求解大迎角翼型繞流非定常N-S方程［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，1999，17(4):466-471.

［19］Xie H，Song W P，Song B F.Airfoil design of a micro-flapping wing based on CFD［J］.Journal of Northwestern Polytechnical University，2008，27(2):227-233.(in Chinese)謝輝，宋文萍，宋筆鋒.基于CFD方法對(duì)微型撲翼翼型設(shè)計(jì)的研究［J］.西北工業(yè)大學(xué)學(xué)報(bào)，2008，27(2):227-233.

［20］Zhu B，Qiao Z D，Gao C.Numerical simulation of transonic flutter for airfoil using approximate boundary conditions［J］.Journal of Northwestern Polytechnical University，2008，26(2):254-258.(in Chinese)朱標(biāo)，喬志德，高超.基于近似邊界條件的翼型跨聲速顫振數(shù)值模擬研究［J］.西北工業(yè)大學(xué)學(xué)報(bào)，2008，26(2):254-258.

Numerical aerodynamic-structural coupling research for flexible flapping wing

Chen Lili1，Song Bifeng2，Song Wenping2，Yang Wenqing2
(1.The AVIC First Aircraft Institute，Xi'an710089 China; 2.National Key Laboratory of Science and Technology on Aerodynamic Design and Research，School of Aeronautics，Northwestern Polytechnical University，Xi'an710072，China)

Due to the small size and light weight，flexible deformation plays an important part in flapping wing aerodynamics.the aero-elastic performances of Flapping-wing Micro Air Vehicle(FMAV) are researched by solving the Reynolds-Averaged Navier-Stokes(RANS)equations and structural dynamic equations.An aerodynamic-structural coupling computational framework is developed，which is able to simulate aerodynamic-structural coupling characteristics of flexible flapping wings.The flapping wing's unsteady aerodynamic characteristics is obtained by solving Reynolds Average Navier-Stokes(RANS)equations.Structural dynamic equations capable of describing the flapping wing's movement are derived by using Hamilton principle，then discreted through finite element method and solved by Newmark solution to get structural characteristics.Loose coupling method is used.On the basis of the above work，CFD/ CSD(Computational Fluid Dynamics/Computational Structural Dynamics)grid data exchange method，dynamic grid technology as well as coupling method are further investigated，and finally a complete set of CFD/CSD coupling solver is developed.Computational results show good agreement with experimental results，which prove that the method developed in this paper are valid and suit for simulation of flexible flapping wing.Effect of inertia loads and kinetic parameters arealso investigated，which can help to understand the mechanisms of flexible flapping wing's aeroelasticity and give a guidance in the design of flexible flapping wing.

micro flapping wing;aerodynamic-structural coupling;Reynolds-averaged Navier-Stokes;structural dynamics;unsteady aerodynamic characteristics

V211.3

Adoi:10.7638/kqdlxxb-2012.0210

0258-1825(2015)01-0125-08

2012-12-19;

2014-12-23

中國(guó)博士后科學(xué)基金(20100481369)

陳利麗(1985-)，女，安徽合肥人，博士研究生，主要從事微型飛行器方向研究.E-mail:lilichen@mail.nwpu.edu.cn

陳利麗，宋筆鋒，宋文萍，等.柔性撲翼氣動(dòng)結(jié)構(gòu)耦合特性數(shù)值研究［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2015，33(1):125-133.

10.7638/kqdlxxb-2012.0210.Chen L L，Song B F，Song W P，et al.Numericalaerodynamic-structural coupling research for flexible flapping wing［J］.Acta Aerodynamica Sinica，2015，33(1):125-133.