錢 靖，張正科，羅時(shí)鈞，劉 鋒

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

0 引 言

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

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

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

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

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

1 流場計(jì)算方法

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

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

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

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

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

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

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

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

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

采用二階隱式差分格式［18］，對振動(dòng)方程（2）進(jìn)行離散，離散方程可以寫為：

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

3 流－固耦合計(jì)算流程

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

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

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

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

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

4.1 彈性薄板的振動(dòng)及其對氣動(dòng)力的影響

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

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

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

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

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

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

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

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

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

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

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

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

5 結(jié) 論

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

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

［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)面法的機(jī)翼氣動(dòng)／結(jié)構(gòu)一體化優(yōu)化設(shè)計(jì)研究［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，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)治.計(jì)算流體力學(xué)基礎(chǔ)［M］.北京：北京航空航天大學(xué)出版社，1989：13.

［19］黃永強(qiáng)，陳樹勛.機(jī)械振動(dòng)理論［M］.北京：機(jī)械工業(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.