王 奇，朱寅鑫，牛培行，劉少寶

（1.南京航空航天大學(xué) 航空學(xué)院 機(jī)械結(jié)構(gòu)力學(xué)及控制國家重點(diǎn)實(shí)驗(yàn)室，南京 210016；2.南京航空航天大學(xué) 航空學(xué)院 多功能輕量化材料與結(jié)構(gòu)工信部重點(diǎn)實(shí)驗(yàn)室，南京 210016）

引言

盡管微型飛行器技術(shù)有了很大的發(fā)展，但昆蟲、鳥類等生物的飛行能力至今仍讓人造飛行器望塵莫及，因此研究撲翼氣動(dòng)性能具有重要意義[1-3].在以往的研究中，通常把撲翼翼型作為剛性翼型處理，忽略了其彈性變形的影響.剛性翼型在超過其相應(yīng)的臨界攻角時(shí)，會(huì)發(fā)生嚴(yán)重的流動(dòng)分離，從而影響飛行器的飛行狀態(tài).而昆蟲、鳥類等生物的翅膀在飛行中往往可以旋轉(zhuǎn)較大角度且仍能保持平穩(wěn)飛行.這不僅與其特殊的撲翼運(yùn)動(dòng)飛行方式有關(guān)，還與其柔性翼型結(jié)構(gòu)有關(guān).例如，Tang、Shyy 等[4]采用梁模型研究了低Reynolds 數(shù)下柔性平板翼型的氣動(dòng)性能，研究表明柔性結(jié)構(gòu)在空氣動(dòng)力載荷作用下發(fā)生變形，表現(xiàn)出類似撲翼的俯仰運(yùn)動(dòng)，產(chǎn)生了明顯的升力和推力.Visbal 等[5]采用高保真隱式大渦模擬方法，研究了表面微柔性薄膜的SD7003 翼型氣動(dòng)性能，研究表明柔性對(duì)膜翼型上流動(dòng)結(jié)構(gòu)使得升力增加和失速延遲.張興偉、周超英等[6]通過弱流固耦合（僅考慮結(jié)構(gòu)變形對(duì)流場(chǎng)的影響）方法分析了柔性撲翼變形的影響，發(fā)現(xiàn)適當(dāng)?shù)娜嵝宰冃慰梢蕴岣呱?、增大升阻?王姝歆、周建華等[7]以昆蟲的翅膀?yàn)榛A(chǔ)，對(duì)柔性翅模型進(jìn)行了合理的假設(shè)與簡化，建立了柔性翅模型并進(jìn)行了實(shí)驗(yàn)研究和分析，結(jié)果表明，柔性翅的柔性變形有助于改善飛行時(shí)的空氣動(dòng)力學(xué)特性.Kang 等[8]研究了前緣處局部柔性翼的繞流問題，發(fā)現(xiàn)局部柔性翼型對(duì)提高升力有著積極的影響.進(jìn)一步，陶真新等[9]數(shù)值分析了吸力面為三段柔性結(jié)構(gòu)的NACA0012 翼型的氣動(dòng)性能，發(fā)現(xiàn)翼面上表面的柔性起到提高升力、抑制流動(dòng)分離的作用.

盡管以上研究已涉及柔性翼型氣動(dòng)性能的研究，但他們往往僅關(guān)注展向柔性變形或局部柔性變形的情況，多簡化為薄膜、梁等簡單結(jié)構(gòu)進(jìn)行建模，而缺乏對(duì)全局柔性翼型氣動(dòng)性能的研究.本研究以柔性橢圓撲翼翼型為例，建立有限元模型，分析不同來流風(fēng)速、迎角下柔性橢圓翼型的空氣動(dòng)力學(xué)特性.

1 剛性撲翼翼型空氣動(dòng)力學(xué)仿真分析

1.1 計(jì)算模型與網(wǎng)格劃分

為了驗(yàn)證建模、仿真分析方法的有效性以及參數(shù)設(shè)置的合理性，首先對(duì)剛性橢圓翼型模型周圍的流場(chǎng)進(jìn)行仿真分析.利用ANSYS Workbench 開展雙向流固耦合計(jì)算，建立剛性翼型周圍的流場(chǎng)模型，采用Meshing模塊劃分網(wǎng)格，調(diào)用計(jì)算流體力學(xué)Fluent 模塊，基于有限差分格式進(jìn)行流場(chǎng)計(jì)算.

為了便于對(duì)比驗(yàn)證，撲翼翼型參考Pesavento 和Wang 所采用的果蠅翅膀尺寸[10]，典型果蠅質(zhì)量為1 mg，平均翼展長約為0.2 cm，橢圓翼型長軸長為0.068 cm，橢圓率為0.25（圖1）.流場(chǎng)尺寸為16 mm × 10 mm，橢圓位于流場(chǎng)中央，翼型迎角為27.5°.流體介質(zhì)為空氣，速度為2.94 m/s，方向垂直于入口邊界水平吹向橢圓翼型，左側(cè)為流場(chǎng)的速度入口，右側(cè)為出口，上下邊界為固定壁面.

圖1 橢圓翼型幾何模型：（a）流體域幾何尺寸；（b）橢圓翼型幾何尺寸Fig.1 The geometric model for the elliptical airfoil:(a)geometric sizes of the fluid domain;(b)geometric sizes of the elliptical airfoil

流場(chǎng)區(qū)域分為兩部分，橢圓壁后尾流區(qū)與外圍流場(chǎng)區(qū)，流場(chǎng)整體采用三角形非結(jié)構(gòu)網(wǎng)格進(jìn)行劃分，如圖2所示.為提高求解精度，同時(shí)也為了更好地呈現(xiàn)橢圓壁后流場(chǎng)，在尾跡區(qū)進(jìn)行網(wǎng)格加密處理；同時(shí)為了提高計(jì)算效率，外圍流場(chǎng)區(qū)采用較為稀疏的網(wǎng)格.考慮到低Reynolds 數(shù)下橢圓壁邊界層的黏性作用，橢圓邊界層使用結(jié)構(gòu)化網(wǎng)格進(jìn)行加密處理，邊界層網(wǎng)格膨脹率為1.1，共設(shè)置20 層，網(wǎng)格厚度從壁面向外逐層遞增，網(wǎng)格總數(shù)為19808.采用正交質(zhì)量檢查項(xiàng)進(jìn)行網(wǎng)格質(zhì)量檢查，網(wǎng)格質(zhì)量達(dá)到0.97.

圖2 流體網(wǎng)格劃分：（a）整體網(wǎng)格劃分；（b）壁面附近網(wǎng)格加密處理Fig.2 Meshing of the fluid domain:(a)meshing of the whole model;(b)mesh refinement near the wall

在Fluent 分析中，設(shè)置瞬態(tài)計(jì)算.由于果蠅的尺寸小，速度較低，在實(shí)際飛行中，其Reynolds 數(shù)范圍一般為102～104，因此在流場(chǎng)分析中主要分析其橢圓翼型低Reynolds 數(shù)下的空氣動(dòng)力學(xué)特性，計(jì)算模型采用k-ω SST 湍流模型[11-12].k-ω SST 湍流模型湍動(dòng)能k和其比耗散率ω 的關(guān)系為

變量k和ω 的控制方程分別為

方程（2a）和（2b）中右側(cè)前三項(xiàng)分別為生成項(xiàng)、耗散項(xiàng)和擴(kuò)散項(xiàng)，方程（2b）右側(cè)最后一項(xiàng)為交叉擴(kuò)散項(xiàng)，其中

式(1)～(4)中，ε為湍動(dòng)能耗散率；ρ為密度；ui,uj,uk為速度分量；μ為層流黏性；νt為湍流運(yùn)動(dòng)黏性；β?，β，γ，σk，σω和 σω2為封閉常數(shù).

采用SIMPLE 的二階格式求解方法對(duì)N-S 離散方程進(jìn)行求解，提高求解精度.時(shí)間步長為0.000 1 s，共計(jì)算100 個(gè)時(shí)間步，計(jì)算得到0.01 s 內(nèi)的結(jié)果.

1.2 模型驗(yàn)證

分別計(jì)算得到該模型流場(chǎng)的速度云圖和壓力云圖（圖3），橢圓翼型后產(chǎn)生了明顯的周期性尾渦脫落，這一特征與Pesavento 和Wang[10]數(shù)值計(jì)算的流場(chǎng)特征一致.

圖3 剛性翼型流場(chǎng)流速分布：（a）速度云圖；（b）壓力云圖Fig.3 Characteristics of the flow field around the rigid airfoil:(a)the velocity contour;(b)the pressure contour

為了便于與文獻(xiàn)結(jié)果定量對(duì)比驗(yàn)證，將時(shí)間無量綱化為

其中，v為來流速度，v=2.94 m/s；t為時(shí)間；c為橢圓翼型的弦長，c=0.68 mm.

數(shù)值仿真計(jì)算所得橢圓翼型升力隨時(shí)間的變化，開始階段由于求解過程未收斂出現(xiàn)較大幅度擾動(dòng)，隨后升力值在0.4～0.6 mg 之間呈現(xiàn)出周期性變化，周期約為4（圖4）.這是因?yàn)橐硇秃缶壋霈F(xiàn)周期性尾渦脫落.升力平均值約為0.5 mg，這是因?yàn)閱沃怀岚虺袚?dān)了果蠅重量.計(jì)算所得升力曲線的這些特征（擾動(dòng)周期、平均升力），與Pesavento 和Wang[10]的數(shù)值計(jì)算結(jié)果保持一致.表明所建立的模型、計(jì)算方法及參數(shù)設(shè)置等是有效的，可供后續(xù)研究.采用三種網(wǎng)格數(shù)量（1.23×104，1.98×104，5.23×104）進(jìn)行了網(wǎng)格無關(guān)性驗(yàn)證，稀疏網(wǎng)格對(duì)升力峰值的監(jiān)測(cè)不夠，中等網(wǎng)格和密網(wǎng)格結(jié)果相近.因此，后續(xù)研究采用中等網(wǎng)格（即網(wǎng)格數(shù)1.98×104）即可滿足計(jì)算精度.

圖4 剛性翼型升力與Pesavento 和Wang[10]的數(shù)值計(jì)算結(jié)果對(duì)比Fig.4 Comparison of lift forces of the rigid airfoil to the numerical results of Pesavento and Wang[10]

1.3 不同風(fēng)速、迎角下的流場(chǎng)及氣動(dòng)性能

選取翼型迎角為0°，6°，12°，16°，20°，24°，28°，32°，風(fēng)速為2 m/s，5 m/s，10 m/s，分別計(jì)算得到橢圓翼型周圍的流場(chǎng)（圖5），壓力場(chǎng)（圖6），升、阻力及升阻比（圖7）.

圖5 不同迎角下的剛性翼型流場(chǎng)速度云圖Fig.5 Velocity contours of the flow field around the rigid airfoil at different attack angles

從速度云圖（圖5）可以看出，隨著攻角增大，橢圓翼型后流動(dòng)分離范圍逐漸增大.當(dāng)迎角較小時(shí)，流動(dòng)為穩(wěn)定的層流；當(dāng)迎角達(dá)到24°，翼型后的流動(dòng)開始出現(xiàn)擾動(dòng)；當(dāng)迎角達(dá)到28°以上時(shí)，翼型后出現(xiàn)周期性的渦脫落.從不同迎角下的壓力云圖（圖6）可以看出，當(dāng)迎角增大到一定程度，由于脫落渦的出現(xiàn)，出現(xiàn)了脫離翼型的負(fù)壓中心.

圖6 不同迎角下的剛性翼型流場(chǎng)壓力云圖Fig.6 Pressure contours of the flow field around the rigid airfoil at different attack angles

從升力曲線圖(圖7(a))可以看出，升力隨來流風(fēng)速的增大而增大；當(dāng)迎角較小時(shí)，升力隨迎角線性增大；隨著迎角進(jìn)一步增大，由于橢圓翼型后流動(dòng)分離加劇，升力隨迎角的增大呈非線性變化，迎角越大升力增大的速率越慢.從阻力曲線圖(圖7(b))可以看出，當(dāng)迎角較小時(shí)，阻力隨著迎角緩慢增大；當(dāng)迎角增大達(dá)到16°左右時(shí)，阻力隨迎角迅速增大.而升阻比隨迎角變化成非單調(diào)趨勢(shì)(圖7(c))，隨迎角的增大先增大再減小，當(dāng)迎角為10°～15°時(shí)達(dá)到最大值.

圖7 不同風(fēng)速下剛性翼型的升力、阻力及升阻比：（a）升力；（b）阻力；（c）升阻比Fig.7 Lift forces,drag forces and lift-drag ratios of the rigid airfoil at different wind speeds:(a)lift forces;(b)drag forces;(c)lift-drag ratios

2 柔性翼型流固耦合仿真分析

2.1 流固耦合模型

在以上剛性翼計(jì)算模型的基礎(chǔ)上，將橢圓翼型考慮為柔性翼.果蠅翅膀上分布著沿展向翅脈，且翅脈的剛度遠(yuǎn)大于翼面.為了簡化模型，在柔性撲翼橢圓翼型模型中心的空心圓柱（直徑為0.02 mm）面，采用固定邊界，以模擬剛性翅脈（圖8）.流體作用力與柔性翼型固體結(jié)構(gòu)相互耦合，一方面流體在固體結(jié)構(gòu)表面產(chǎn)生作用力，使固體結(jié)構(gòu)產(chǎn)生變形；結(jié)構(gòu)形狀的改變，引起流場(chǎng)結(jié)構(gòu)的改變.流固耦合仿真分析在ANSYS Workbench 平臺(tái)中用Fluent 計(jì)算流體模塊進(jìn)行流場(chǎng)分析，用Transient Structural 瞬態(tài)結(jié)構(gòu)模塊進(jìn)行結(jié)構(gòu)分析，使用System Coupling耦合求解器將兩個(gè)模塊關(guān)聯(lián)起來，雙向傳遞流體的計(jì)算數(shù)據(jù)和固體結(jié)構(gòu)的計(jì)算數(shù)據(jù)，將流場(chǎng)與結(jié)構(gòu)變形進(jìn)行耦合.在流固耦合分析時(shí)，需要建立流固界面進(jìn)行數(shù)據(jù)傳遞，因此柔性翼型模型應(yīng)為三維模型.

圖8 柔性橢圓翼型有限元模型：（a）網(wǎng)格劃分；（b）流固耦合面與固定面Fig.8 The finite element model for the flexible elliptical airfoil:(a)the meshing;(b)the fluid-solid coupling and the fixed surface

柔性翼型固體結(jié)構(gòu)分析的動(dòng)力學(xué)控制方程[9]可寫為

式中M,C,K分別為質(zhì)量矩陣、阻尼矩陣、剛度矩陣；F(t)為流體作用在固體結(jié)構(gòu)表面上的瞬態(tài)耦合力，其對(duì)結(jié)構(gòu)的綜合作用可以是力，也可以是力矩.柔性材料密度為1 000 kg/ m3，材料本構(gòu)模型采用線彈性模型，Poisson 比為0.3，彈性模量可調(diào).

2.2 網(wǎng)格劃分

流固耦合計(jì)算在Fluent 中采用動(dòng)網(wǎng)格技術(shù)，采用非結(jié)構(gòu)網(wǎng)格進(jìn)行劃分.在流固耦合計(jì)算中，當(dāng)固體結(jié)構(gòu)發(fā)生變形后，流場(chǎng)需要接收固體變形的位移數(shù)據(jù)，流場(chǎng)邊界也相應(yīng)的改變，因此需要使用動(dòng)網(wǎng)格進(jìn)行網(wǎng)格重構(gòu).具體為在Dynamic Mesh 中設(shè)置Smoothing 與Remeshing，將橢圓壁面設(shè)為System Coupling，端部壁面設(shè)為Deforming，其余設(shè)置與剛性翼型一致.固體結(jié)構(gòu)也采用非結(jié)構(gòu)網(wǎng)格進(jìn)行劃分.

2.3 模型驗(yàn)證

在驗(yàn)證模型中，橢圓翼型迎角設(shè)置為27.5°，來流速度為2.94 m/s.柔性翼型彈性模量為0.01 MPa，0.1 MPa，1000 MPa，分析流固耦合模型的流場(chǎng)和固體結(jié)構(gòu)變形.

從翼型周圍流場(chǎng)速度云圖（圖9）可以看出，隨著橢圓翼型彈性模量的減小，翼型周圍的最大流速減小，但是翼型對(duì)周圍流場(chǎng)的擾動(dòng)范圍擴(kuò)大.類似地，從翼型周圍流場(chǎng)壓力云圖（圖10）可以看出，隨著橢圓翼型彈性模量的減小，翼型周圍產(chǎn)生的最大負(fù)壓值減小，但是翼型周圍的負(fù)壓區(qū)擴(kuò)大.從柔性翼型結(jié)構(gòu)位移云圖（圖11）可以看出，隨著橢圓翼型彈性模量的減小，柔性翼型結(jié)構(gòu)位移增大.

圖9 不同彈性模量的柔性翼型流場(chǎng)速度云圖Fig.9 Velocity contours of the flow field around the flexile airfoil with different Young’s moduli

圖10 不同彈性模量的柔性翼型流場(chǎng)壓力云圖Fig.10 Pressure contours of the flow field around the flexile airfoil with different Young’s moduli

圖11 不同彈性模量的柔性翼型位移云圖Fig.11 Deformation contours of the flexile airfoil with different Young’s moduli

從升力曲線圖（圖12（a））可以看出，柔性橢圓翼型升力開始階段由于求解過程未收斂出現(xiàn)較大幅度波動(dòng)，隨后升力值穩(wěn)定在0.4～0.6 mg 之間呈現(xiàn)出周期性變化，總體變化趨勢(shì)與剛性翼型一致.但隨著彈性模量的降低，升力值穩(wěn)定擾動(dòng)振蕩的周期增大.而隨著彈性模量的降低，擾動(dòng)振蕩的振幅呈先減小再增大的趨勢(shì).當(dāng)彈性模量降低到0.1 MPa 時(shí)，升力值擾動(dòng)振蕩幅度幾乎消失，升力值達(dá)到一穩(wěn)定值（約0.5 mg）；而當(dāng)彈性模量進(jìn)一步降低時(shí)，升力值又呈現(xiàn)出周期性的震蕩變化.由此可見，較剛性翼型，柔性翼型對(duì)于穩(wěn)定升力擾動(dòng)有著積極的作用.類似地，從阻力曲線圖（圖12（b））中可以看出，較剛性翼型，柔性翼型的阻力略大于剛性翼.且隨著彈性模量的降低，柔性翼型的阻力逐漸出現(xiàn)了周期性振蕩.

圖12 柔性翼型的升力和阻力：（a）升力；（b）阻力Fig.12 Lift forces and drag forces of the flexible airfoil:(a)lift forces;(b)drag forces

2.4 不同迎角下流場(chǎng)及氣動(dòng)性能

柔性翼型的迎角分別選取0°，6°，12°，16°，20°，24°，28°，32°進(jìn)行仿真分析.風(fēng)速為5 m/s，柔性翼型彈性模量為0.1 MPa.

由翼型周圍流場(chǎng)（圖13）可以看出，隨著迎角逐漸增大，剛性翼型后的流場(chǎng)逐漸出現(xiàn)擾動(dòng)，翼型后伴隨有尾渦脫落.而柔性翼型后的流場(chǎng)與剛性翼型有很大不同，流體流過柔性翼型后盡管出現(xiàn)了流動(dòng)分離，但流場(chǎng)并未出現(xiàn)明顯的擾動(dòng).由此可見，較剛性翼型，柔性翼型延緩了尾渦脫落時(shí)間，顯著抑制了尾流流場(chǎng)的擾動(dòng)范圍及幅度.

圖13 不同迎角下的剛性翼型（左）與柔性翼型（右）的周圍流場(chǎng)Fig.13 Velocity contours of the fluid field around the rigid(left)and the flexible(right)airfoil with different attack angles

由柔性翼型的升力曲線（圖14（a））可以看出，升力隨著迎角增大而增大，且柔性翼型的升力略大于剛性翼型.由柔性翼型的阻力曲線（圖14（b））可以看出，阻力隨迎角的增大而增大，且柔性翼型的阻力略大于剛性翼型.但從升阻比（圖14（c））看，柔性翼型與剛性翼型并無相顯著差異.也就是說，柔性翼型并未損失飛行效率.

圖14 柔性翼型的升力、阻力及升阻比隨迎角變化：（a）升力；（b）阻力；（c）升阻比Fig.14 Lift forces,drag forces and lift-drag ratios of the flexible airfoil vs.the attack angle:(a)lift forces;(b)drag forces;(c)lift-drag ratios

3 結(jié)論

本文從仿生學(xué)的角度出發(fā)，通過流固耦合數(shù)值仿真分析了柔性翼型在不同風(fēng)速、迎角下的空氣動(dòng)力學(xué)特性.研究發(fā)現(xiàn)，較剛性翼型，柔性翼型延緩了尾渦脫落時(shí)間，有效降低升力擾動(dòng)振蕩頻率；柔性翼型顯著抑制了尾流流場(chǎng)的擾動(dòng)，降低升力擾動(dòng)振蕩幅值，合適的彈性模量翼型使得擾動(dòng)振蕩完全消除，這將有利于改善飛行器運(yùn)動(dòng)平順性、降低氣動(dòng)噪聲.

致謝本文作者衷心感謝江蘇省仿生功能材料重點(diǎn)實(shí)驗(yàn)基金對(duì)本文的資助.