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

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

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

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

(1.中國航空工業(yè)集團第一飛機設(shè)計研究院，西安710089;
2.西北工業(yè)大學翼型葉柵空氣動力學國家重點實驗室，西安710072)

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

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

0 引言

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

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

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

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

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

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

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

1 流動控制方程

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

其中，

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

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

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

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

利用哈密頓原理，撲翼的結(jié)構(gòu)運動方程可以表示為:其中，t1，t2表示任意時刻，δU、δW、δT分別代表應(yī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。

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

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

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

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

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

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

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

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

其中，運動參數(shù)α、ψ均為時間t的函數(shù)。

圖4 撲翼坐標系Fig.4Flapping wing coordinate system

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

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

對位置坐標進行求導，可以得到任一點的速度表達式:

其中ρ為材料的密度。

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

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

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

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

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

3CFD/CSD耦合方法

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

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

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

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

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

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

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

4 算例驗證與分析

對柔性撲翼模型進行了數(shù)值模擬，并將計算結(jié)果與風洞實驗結(jié)果進行比較。

4.1 撲翼運動模式

撲翼機構(gòu)的運動為揮舞運動，如圖5所示。

圖5 撲翼運動示意圖Fig.5Flapping wing movement

撲翼運動規(guī)律表示如下:

其中，ψ的取值均以上反為正，Ψ1為揮舞運動的振幅。

4.2 算法驗證

4.2.1 網(wǎng)格和時間精度驗證

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

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

圖7 不同時間步長計算結(jié)果比較Fig.7Comparison of result with different time steps

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

4.2.2 算例驗證

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

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

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

計算結(jié)果與實驗測量結(jié)果的比較如圖11所示。

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

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

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

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

圖11 計算結(jié)果與實驗值比較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可以看出，柔性撲翼在撲動過程中均會產(chǎn)生柔性彎曲變形和扭轉(zhuǎn)變形。在撲動頻率為3Hz時，當撲翼從最高點往下?lián)鋾r，在最高點處，撲翼產(chǎn)生了上翹的柔性彎曲變形;從最高點撲動到平衡位置過程中，撲翼的彎曲變形幅度逐漸增大，在平衡位置附近，撲翼變形達到最大;從平衡位置撲動到最低點過程中，撲翼上翹變形逐漸反彈，在最低位置附近，撲翼變形達到最小;從最低點上撲到平衡位置過程中，撲翼產(chǎn)生向下的彎曲變形，變形量先增大后減小，在7/8T時刻達到最小，從7/8T時刻撲動到最高點過程中，撲翼產(chǎn)生上翹的變形。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

由此可知，增大揮舞角有利于柔性撲翼推力特性的產(chǎn)生，對柔性撲翼升力特性則幾乎不產(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)變形、氣動特性兩個方面研究了兩個關(guān)鍵運動參數(shù)(撲動頻率、揮舞角)對柔性撲翼氣動結(jié)構(gòu)耦合特性的影響、結(jié)構(gòu)變形對柔性撲翼氣動特性的影響以及柔性撲翼隨運動參數(shù)的變化趨勢，得到了柔性撲翼氣動特性方面的一些基本規(guī)律，研究結(jié)果在各節(jié)中均有小結(jié)。從研究結(jié)果可以得出，柔性撲翼產(chǎn)生的結(jié)構(gòu)變形對氣動特性有顯著影響，對柔性撲翼氣動特性進行數(shù)值模擬時，考慮結(jié)構(gòu)變形的影響是非常必要的。與揮舞角相比，撲動頻率對柔性撲翼的氣動結(jié)構(gòu)特性影響表現(xiàn)的更為明顯。

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

［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)曾銳，昂海松，梅源，等.撲翼柔性及其對氣動特性的影響［J］.計算力學學報，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è)大學學報，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)宋文萍，楊永，喬志德，等.雙時間法求解大迎角翼型繞流非定常N-S方程［J］.空氣動力學學報，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方法對微型撲翼翼型設(shè)計的研究［J］.西北工業(yè)大學學報，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)朱標，喬志德，高超.基于近似邊界條件的翼型跨聲速顫振數(shù)值模擬研究［J］.西北工業(yè)大學學報，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

中國博士后科學基金(20100481369)

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

陳利麗，宋筆鋒，宋文萍，等.柔性撲翼氣動結(jié)構(gòu)耦合特性數(shù)值研究［J］.空氣動力學學報，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.