二維彈性撲翼沉浮運動流動特性的數(shù)值研究

錢 靖，張正科，羅時鈞，劉 鋒

（1.翼型葉柵空氣動力學國防科技重點實驗室，西北工業(yè)大學，陜西 西安710072；2.加利福尼亞大學爾文分校機械與航空航天工程系，爾文，加利福尼亞，CA 92697－3975，美國）

0 引 言

由于微型飛行器（指尺寸小于20cm，巡航范圍大于10km，巡航時間大于20min的飛行器）有尺寸小等優(yōu)點，因而在諸如低空軍事偵察、戰(zhàn)場損傷評估、目標搜索、通訊中繼、生化探測等軍事領(lǐng)域和諸如交通監(jiān)控、邊境巡邏、野生勘察、森林防火等民用領(lǐng)域有廣闊的應(yīng)用前景，在上世紀末在國際上形成研究熱潮。同時近些年微電子技術(shù)和微型電腦技術(shù)的發(fā)展也使微型飛行器的研制和投入實際應(yīng)用成為了可能。

鳥類、昆蟲和魚類的撲翼推進方式激發(fā)了人們極大的興趣［1］。模擬鳥類、昆蟲和魚類的運動啟發(fā)了飛機的發(fā)明，但固定翼飛機只借鑒了鳥類滑翔的原理，卻放棄了撲翼驅(qū)動的方式。而當微型飛行器的尺度在15cm以下時，如仍采用固定翼的常規(guī)氣動布局，則由于雷諾數(shù)太低，達不到足夠的升阻比，而無法飛行。于是人們又回過頭來考慮撲翼驅(qū)動方式，不管是厘米級還是毫米級的昆蟲，都是采用撲翼的方式而具備了復(fù)雜多變的飛行能力。模擬鳥類、昆蟲和魚類運動方式的撲翼微型飛行器有更好的飛行性能，其機動性、靈活性及低能耗都優(yōu)于固定翼和旋轉(zhuǎn)翼微型飛行器。撲翼飛行方式成為微型飛行器首選的一種布局方式。

因而，對撲翼飛行器、撲翼的空氣動力效果的研究，包括實驗和數(shù)值模擬研究，也就成為一件很有理論價值和現(xiàn)實應(yīng)用價值的工作。Knoller［2］和Betz［3］分別在1909年和1912年注意到撲動的翼型會產(chǎn)生推力。1922年，Katzmayr［4］從實驗觀測中證實了此現(xiàn)象－－Knoller－Betz效應(yīng)。此后，對撲翼氣動特性的研究逐漸開展起來。早期的實驗和理論研究可見于文獻［5－7］。如今，用程序或商業(yè)軟件（比如Fluent）可以數(shù)值模擬剛性撲翼的非定常運動。所以對剛性撲翼的飛行性能已經(jīng)有了比較清楚的了解。而彈性撲翼與剛性撲翼有很大的差別。鳥類在每個拍翼周期中，可以通過主動地扭轉(zhuǎn)和彎曲翅膀，以改變翼展，從而獲得更大的推力、減小阻力。其翅膀隨流場被動的變化還可以防止流動的分離，提高升／阻比，從而有利于在起降等非定常運動情況的飛行［8］。M.S.Triantafyllou，G.S.Triantafyllou和 Yue［9］的研究顯示，魚體的彈性對魚類在水中的靈活性起到至關(guān)重要的作用。Heathcote等［10－11］用后緣帶柔性薄板的翼型在水槽中進行了沉浮運動實驗觀察，顯示具有適當彈性的薄板能使翼型產(chǎn)生推力。Tang等［12］用SIMPLE方法計算不可壓N－S方程求解流場，用Euler－Bernoulli梁振動微分方程求解薄板振動運動，并將流體流動、結(jié)構(gòu)變形耦合迭代，計算結(jié)果顯示，做沉浮運動的翼型因氣動載荷產(chǎn)生變形，改變了有效攻角，從而造成升力、阻力產(chǎn)生明顯改變，其位移隨時間變化規(guī)律與Heathcote和Gursul［11］的實驗結(jié)果一致。

彈性撲翼非定常流場數(shù)值模擬牽扯到結(jié)構(gòu)的彈性變形運動及與流場解的耦合，因而是一個非常復(fù)雜的問題。彈性薄板的變形非常明顯地增加了計算量和計算難度。N－S方程方法、Euler方程方法都要生成復(fù)雜的網(wǎng)格，對于剛性撲翼，非定常運動及網(wǎng)格變形已經(jīng)使計算量變得很大，如果再考慮彈性變形及彈性運動的收斂，則計算量將大幅增加。相比之下，勢流板塊法不需要生成網(wǎng)格，從而極大地提高了計算效率。在低速情況下，勢流方法的計算結(jié)果和實驗結(jié)果出 入 不 大。Anderson［13］，Anderson，Streitlien，Barrentt和 Triantafyllou［14］比較了NACA0012翼型水洞實驗與非線性不可壓勢流法計算結(jié)果，結(jié)果顯示當無明顯前緣渦脫落時，計算結(jié)果和實驗結(jié)果符合很好。而當翼型做高頻率小振幅振動時，前緣渦脫落對氣動力的影響會較?。?5］。因此，本文選用勢流板塊法來計算高頻率小振幅振動翼型的非定常氣動力，應(yīng)當具有足夠的計算精度，同時計算效率遠高于Euler方程方法和N－S方程方法。另外，對彈性薄板變形的準確描述較為困難。薄板的變形由慣性力、氣動力和彈性力決定。慣性力是給定的，而氣動力和彈性力又與薄板的變形和運動狀態(tài)相關(guān)。這無疑增加了薄板變形的計算難度。求解此氣動彈性問題的基本思想有兩種，即強耦合法（tight coupling method）和弱耦合法（loose coupling method）。在強耦合法中，按照未收斂的氣動力來計算變形，再把變形計入氣動力計算的迭代過程中，直至變形和氣動力都收斂。在弱耦合法中，只有在氣動力收斂后才計算變形。再計入變形，重新計算氣動力，直至變形收斂。

本文通過數(shù)值方法準確描述附加在翼型后緣的彈性薄板在氣動力作用下的彈性變形運動及回過頭來與流場產(chǎn)生的流－固耦合過程以及這種耦合產(chǎn)生的推力效應(yīng)。薄板變形根據(jù)Euler－Bernoulli梁振動微分方程進行數(shù)值求解。流場用非定常勢流板塊法計算。并且通過弱耦合［16］迭代算法，使得在變化的氣動載荷下，流場求解和懸臂梁微分方程求解的流－固耦合過程保持收斂，從而準確地描述薄板的彈性變形。

1 流場計算方法

通過非定常勢流板塊法［17］計算翼型的非定常氣動載荷。在文獻［17］中，假想每個時間步脫落的渦量并不直接變化為自由渦，而是如圖1所示，先分布在與翼型后緣相連的一個虛擬的直線尾跡板塊上，在下一個時間步，尾跡板塊上的全部渦量才會以自由渦的形式脫落到翼型尾跡中（如圖1所示）。

圖1 勢流板塊法示意圖Fig.1 Schematic diagram of the potential flow panel method

根據(jù)流動的邊界條件，物面上的n個板塊應(yīng)滿足板塊中點處流體法向速度等于物面法線速度的邊界條件，從而構(gòu)成關(guān)于單位長度源強和單位長度渦強的線性方程組。流動在后緣處還需滿足一個Kutta條件，即后緣處上下表面壓力連續(xù)，據(jù)此可由非定常伯努利方程得到一個后緣上下表面速度平方差與翼型環(huán)量變化率的關(guān)系式，成為一個附加方程。

根據(jù)Helmholtz渦量定理，在第k時間步，翼型總環(huán)量與尾跡板塊渦強度之和，等于上一個時間步中的翼型總環(huán)量：

其中，下標k指時間步，Γk為第k時間步翼型總環(huán)量，（γw）k為第k時間步尾跡板塊的單位長度渦強度，Δk為尾跡板塊長度。此外，尾跡板塊長度Δk及其與來流方向的夾角Θk與尾跡板塊中點處誘導(dǎo)速度的x，y分量（Uw）k，（Vw）k存在一定的關(guān)系。

解得翼型渦強、源強分布后，可以求出沿每一板塊的流體速度，進而根據(jù)非定常伯努利方程求得壓力，積分后可得到翼型的升、阻力。

2 結(jié)構(gòu)振動求解方法

翼型后緣連接的薄板相對于前部剛性翼型做彈性振動運動，設(shè)變形屬于小變形，則薄板的相對彈性運動滿足Euler－Bernoulli梁振動微分方程，無量綱梁振動微分方程為：

其中，＝I／c3＝／12，＝E／ρfU，＝ρs／ρf，＝b／c（x，t）＝q（x，t）／ρfU，c為翼型的弦長，ρf為流體密度，U∞為遠場自由來流速度，w為撓度，ρs為薄板密度，b為薄板厚度，E為材料楊氏模量，I為薄板截面的轉(zhuǎn)動慣量，q（x，t）代表薄板受到的分布載荷，包括慣性力和流體施加的壓強。

采用二階隱式差分格式［18］，對振動方程（2）進行離散，離散方程可以寫為：

薄板的邊界條件為：在薄板固定端，撓度和撓度的一階導(dǎo)數(shù)為零；在自由端，由于沒有集中載荷或力矩，所以自由端撓度的二階導(dǎo)數(shù)和三階導(dǎo)數(shù)為零。對邊界條件進行離散，與離散化的振動方程（3）進行聯(lián)立求解，可以得到薄板振動方程組。為了驗證上述格式及程序，取大小隨時間正弦變化的均布載荷，數(shù)值計算薄板自由振動頻率，得到與理論解［19］相一致的結(jié)果，說明差分方程（3）具有足夠計算精度，其結(jié)果是可信的。

3 流－固耦合計算流程

采用弱耦合迭代法［16］（loose coupling method），在氣動力收斂后計算變形，計入變形后再計算氣動力，如此反復(fù)，直至變形收斂。

每個時間步內(nèi)的程序流程為：（1）首先，輸入上一迭代步的壓力分布（第一迭代步采用上一時間步的壓力分布）和這個時間步的慣性力載荷；（2）然后，計算彈性薄板的變形；（3）接著，根據(jù)薄板的變形，計算翼型的外形，從而得到流場的新邊界條件；（4）進一步，根據(jù)新的邊界條件，求解流場的壓力分布；（5）最后，返回到第一步，再計算出彈性薄板新的變形，將此迭代步得到的變形和上一迭代步得到的變形進行比較，如果兩者相差較多，則繼續(xù)迭代，如果兩者接近，則停止迭代，程序進入下一個時間步。

4 計算結(jié)果與分析

在剛性NACA0012翼型后緣接一段長為25%弦長的彈性薄板，并用四階Bezier曲線對連接處進行光滑處理，作為本文算例使用的計算外形。該翼型沉浮運動規(guī)律為：

其中，h＝y(tǒng)為位移＝h／c無量綱位移，為無量綱振幅，k＝ωc／U∞為無量綱角頻率，ω為翼型簡諧振動角頻率。本文在所有的算例中，取ha＝0.018，k＝4.3，薄板無量綱密度為1500、無量綱厚度為0.001。

4.1 彈性薄板的振動及其對氣動力的影響

取薄板無量綱楊氏模量＝6×107，8×107，1×108，1.4×108，2.0×108等進行了計算，以薄板后緣的運動作為指標來考察薄板的運動。圖2給出了不同楊氏模量下薄板后緣相對于翼型前緣的位移隨時間的變化歷程。由圖可知，薄板后緣振動幅度隨薄板無量綱楊氏模量的增加先增加后減小。當＝1×108時，薄板振動幅度最大。并且，不同彈性模量薄板振動的相位也各不相同，隨著薄板值的增加，薄板振動位移達到最大的時刻逐漸提前。這種振動模式和 Heathcote等［10－11］觀察到的是一致的，和Tang等［12］計算的相對厚度為0.00056的結(jié)果也非常接近。

圖2 翼型前緣位移、翼型坐標下薄板后緣位移隨時間的變化Fig.2 Variation of elastic plate trailing edge displacement with time

圖3（a、b）分別給出了四個不同薄板彈性模量下計算的翼型阻力、升力系數(shù)隨時間的變化，可以看出，升力、阻力系數(shù)隨時間近似正弦變化。這個變化規(guī)律與Tang等［12］的計算結(jié)果很接近。

4.2 阻力系數(shù)修正與推進效率

勢流理論沒有考慮流體的粘性，因此，翼型有效推力（useful thrust）應(yīng)為勢流理論求得的推力減去粘性阻力。這里將翼型近似看作平板，分別用平板層流和湍流附面層單位長度表面摩擦阻力系數(shù)理論公式［20］：

對本文計算的勢流阻力系數(shù)進行修正。假定原NACA0012翼型的雷諾數(shù)Re為50000，則本文因使翼型加長了0.25倍，故應(yīng)取雷諾數(shù)為1.25×50000來計算摩擦阻力系數(shù)。

圖5（a、b）為分別經(jīng)層流、湍流附面層修正后的翼型平均阻力系數(shù)隨薄板無量綱楊氏模量的變化。由圖可看出，薄板為剛性時翼型的平均阻力系數(shù)（點劃線所示）為正，說明剛性撲翼不易產(chǎn)生有效推力；對于彈性薄板，存在一個阻力為負值的值區(qū)域，并且總可以找到一個值使翼型阻力系數(shù)達到最大負值。這與Heathcote等［10］的實驗結(jié)果是符合的。

推進效率可以定義為翼型產(chǎn)生的推進功率與輸入給翼型的總機械功率之比［21－22］。圖6（a、b）給出了分別由經(jīng)層流和湍流附面層修正后的推力所得的翼型推進效率隨薄板無量綱楊氏模量的變化規(guī)律。從圖中可看出，選取適當值的薄板，翼型可以達到最佳的推進效率。

4.3 尾跡脫落渦結(jié)構(gòu)

對不同楊氏模量＝8×107，9×107，1×108，1.4×108，2.4×108，2.8×108進行了翼型尾跡脫落點渦分布的計算觀察研究。圖7僅給出了典型楊氏模量值＝1×108下翼型尾跡點渦環(huán)量分布結(jié)果。圖中每個點渦環(huán)量的正負和大小通過每個點的顏色來表示。計算結(jié)果說明，所有楊氏模量下，環(huán)量為正（順時針旋轉(zhuǎn)），且環(huán)量值較大的點渦周期性地聚集在中心線（y＝0）附近偏下方，而環(huán)量為負（逆時針旋轉(zhuǎn)），且環(huán)量值較大的點渦周期性地聚集在中心線（y＝0）附近偏上方；而環(huán)量值較小的點渦則遠離中心線分布。因此，中心線附近流場的流速明顯大于自由來流流速。尾跡中點渦的這種環(huán)量分布規(guī)律與Jones，Dohring，Platzer的理論分析和實驗觀察［6］完全一致，也與Heathcote等的實驗結(jié)果［10－11］類似。并且＝108時的尾跡點渦最大環(huán)量值是所有值中最大的，這恰好是圖5（a）中推力系數(shù)最大的值，說明翼型尾跡中點渦環(huán)量最大值較大時，撲翼的平均推力系數(shù)也較大。

圖6 經(jīng)層流和湍流附面層摩阻修正后推進效率隨薄板無量綱楊氏模量的變化Fig.6 Airfoil propulsive efficiency dependency on nondimensional Young's modulus

圖7 翼型尾跡中的脫落點渦（＝108）Fig.7 Vortex distribution in airfoil wake

5 結(jié) 論

本文采用非定常勢流板塊法和Euler－Bernoulli梁振動微分方程數(shù)值模擬后緣帶彈性薄板的翼型在作高頻率小振幅沉浮運動時翼型繞流和薄板彈性振動之間的流固耦合問題及撲動運動的推力效應(yīng)。從計算結(jié)果可得如下結(jié)論：（1）翼型在做高頻小振幅振動時，非定常勢流板塊法具有較高的計算效率；（2）翼型尾部連接不同楊氏模量的彈性薄板，得到的氣動力的大小和分布會有明顯不同，翼型推力會隨薄板楊氏模量先增加后減小，調(diào)整彈性薄板的剛度，可以使翼型產(chǎn)生最佳的非定常推力；（3）翼型推進效率同樣隨彈性薄板的楊氏模量值先增加后減小，可以找到合適的楊氏模量值，使翼型的推進效率最佳；（4）彈性薄板的楊氏模量值同樣影響尾跡脫落點渦的強度，翼型推力最大時所對應(yīng)的尾跡脫落點渦的環(huán)量值域較寬。（5）翼型坐標下，薄板后緣振動幅度隨薄板楊氏模量的增加先增加后減小，并且振動相位也會發(fā)生變化，即隨著薄板楊氏模量值的增加，薄板振動位移達到最大的時刻逐漸提前。

致謝：本文第一作者感謝蔡晉生教授的耐心指點和劉亞師兄、詹磊師兄、馮春娟師姐和徐嘉師兄的幫助。

［1］SPEDDING G R，LISSAMAN P B S.Technical aspects of microscale flight systems［J］.JournalofAvianBiology，1998，29（4）：458－468.

［2］KNOLLER R.Die gesetze des luftwiderstandes［J］.Flug－undMotortechnik（Wein），1909，3（21）：1－7.

［3］BETZ A.Ein beitrag zur erkl?erung des segelfluges［J］.ZeitschriftfuerFlugtecnikundMotorluftschiffahrt，1912，3（21）：269－272.

［4］KATZMAYR R.Effect of periodic changes of angle of attack on behaviour of airfoils［R］.NACA Rept.147，Oct.1922.

［5］KOOCHESFAHANI M M.Vortical patterns in the wake of an oscillating airfoil［J］.AIAAJournal，1989，27（9）：1200－1205.

［6］JONES K D，DOHRING C M，PLATZER M F.Experimental and computational investigation of the Knoller－Betz effect［J］.AIAAJournal，1998，36（7）：1240－1246.

［7］LAI J C S，PLATZER M F.Jet characteristics of a plunging airfoil［J］.AIAAJournal，1999，37（12）：1529－1537.

［8］SHYY W，BERG M，LJUNGQVIST D.Flapping and flexible wings for biological and micro air vehicles［J］.ProgressinAerospaceSciences，1999，35：455－505.

［9］TRIANTAFYLLOU M S，TRIANTAFYLLOU G S，YUE D K P.Hydrodynamics of fishlike swimming［J］.AnnualReviewofFluidMechanics，2000，32：33－53.

［10］HEATHCOTE S，MARTIN D，GURSUL I.Flexible flapping airfoil propulsion at zero freestream velocity［J］.AIAAJournal，2004，42（11）：2196－2204.

［11］HEATHCOTE S，GURSUL I.Flexible flapping airfoil propulsion at low reynolds number［A］.In：43rd AIAA Aerospace Sciences Meeting and Exhibit［C］.Reno，Nevada.AIAA－2005－1405，Jan 2005.

［12］TANG J，VIIERU D，SHYY W.A Study of aerodynamics of low reynolds number flexible airfoils［A］.In：37th AIAA Fluid Dynamics Conference and Exhibit［C］.Miami，F(xiàn)L.AIAA－2007－4212，June 2007.

［13］ANDERSON J M.Vorticity control for effcient propulsion［D］.Massachusetts Institute of Technology and Woods Hole Oceanographic Institution，1996.

［14］ANDERSON J M，STREITLIEN K，BARRENTT D S，et al.Oscillating foils of high propulsive efficiency［J］.FluidMech.，1998，360：41－72.

［15］YOUNG J，LAI J C S.Oscillation frequency and amplitude effects on the wake of a plunging airofil［J］.AIAA Journal，2004，42：2042－2052.

［16］劉金輝，喬志德，楊旭東，等.基于響應(yīng)面法的機翼氣動／結(jié)構(gòu)一體化優(yōu)化設(shè)計研究［J］.空氣動力學學報，2006，24（3）：300－306.

［17］TENG N－H.The development of a computer code for the numerical solution of unsteady，inviscid and incompressible flow over an airfoil［D］.Naval Postgraduate School，Monterey，California，1987：

［18］劉導(dǎo)治.計算流體力學基礎(chǔ)［M］.北京：北京航空航天大學出版社，1989：13.

［19］黃永強，陳樹勛.機械振動理論［M］.北京：機械工業(yè)出版社，1996：123.

［20］ANDERSON Jr J D.Fundamentals of aerodynamics［M］.Boston：McGraw Hill Companies，Inc，2007 ：883－925.

［21］YANG S，LUO S，LIU F.Computation of the flows over flapping airfoil by the euler equations［A］.In：43rd AIAA Aerospace Sciences Meeting［C］.Reno，NV，AIAA－2005－1407，Jan 2005.

［22］LAI A，LIU F.Computation of a rigid airfoil in plunging motion with a flexible trailing beam［A］.In：47th AIAA Aerospace Sciences Meeting and Exhibit［C］.Jan.2009，Orland，F(xiàn)lorida，AIAA Paper 2009－0726.