昆蟲振翅飛行原理的數(shù)值模擬研究

1、<p>　　昆蟲振翅飛行原理的數(shù)值模擬研究</p><p>　　The Numerical Simulation of Insects Flight</p><p>　　北京大學 力學與工程科學系</p><p>　　理論與應用力學專業(yè) 98級 程暮林</p><p><b>　　摘要</b></p

2、><p>　　昆蟲經(jīng)過漫長的進化過程，獲得了令人驚嘆的飛行技巧，其高度的飛行機動性令現(xiàn)代飛行器黯然失色。隨著非定常流理論的不斷完善，人們逐漸認識到非定常效應在昆蟲飛行中起著重要作用。前人采用實驗和無粘假定下的理論分析的方法對此進行了一些研究。然而由于該問題本身的復雜性，使得理想化的理論分析和試驗的研究手段難于得到全面和深入的結(jié)果。本文采用流體力學數(shù)值模擬方法研究本問題，通過采取轉(zhuǎn)換參考系的方法對昆蟲的振翅飛行所擾動而

3、產(chǎn)生的流場進行了數(shù)值模擬，克服了直接對此問題進行數(shù)值模擬所必然遇到的困難－“動邊界問題”，縮短了計算時間。文中通過改變控制飛行的各個參數(shù)，從而比較詳細的研究了這些參數(shù)對昆蟲飛行的影響，為研究昆蟲的飛行控制機制提供了依據(jù)。</p><p><b>　　Abstract</b></p><p>　　Insects have gained very surprising a

4、bility of flight after thousand years of evolution. People have come to realize the important effect of unsteady flow in getting enough lift for insects with the improvement of Unsteady Flow Theory. We have got some use

5、ful results with experimental and theoretical methods, but now it is difficult to go futhure with such methods. In this article , I use the numerical simulation method (the standard k-εturbulence model) in a non-inertia

6、 reference linked to wi</p><p><b>　　前言。</b></p><p>　　經(jīng)過幾千萬年漫長的進化過程，昆蟲獲得了令人驚嘆的飛行本領，有的昆蟲能夠懸停于空中，有的能夠以很小的曲率半徑急速轉(zhuǎn)彎等等。這種高度靈活的飛行控制能力是人類現(xiàn)有的飛行機械所無法比擬的。因此，關于昆蟲振翅飛行原理的深入研究將有助于提高現(xiàn)有的飛行器的性能，推動微型飛

7、行器的開發(fā)和研制。</p><p>　　昆蟲學家們很早就對昆蟲的飛行發(fā)生了濃厚的興趣，隨著技術的發(fā)展，他們采用高速攝影，攝像技術對各種昆蟲的飛行進行了大量的拍攝，并且同時測得了一些與昆蟲飛行相關的參數(shù)（如振翅頻率，振翅的振幅等等）【1~4】，通過對這些數(shù)據(jù)進行數(shù)據(jù)擬合等處理后，他們發(fā)現(xiàn)這些參數(shù)間存在一定的關系，例如，翅的表面積正比于昆蟲質(zhì)量的2/3方，而振翅頻率正比于昆蟲質(zhì)量的－1/4方。雖然這些規(guī)律均是比較粗糙

8、的，但能夠證明這樣一個觀點：雖然經(jīng)過上千萬年的漫長進化過程，由于環(huán)境的差異導致各種昆蟲彼此的形態(tài)，大小，習性均有較大的差異，然而由于受到“振翅飛行”這一物理條件的限制，各種昆蟲的身體結(jié)構，翅的形態(tài)，翅的運動模式等飛行參數(shù)必然符合一定的力學規(guī)律。然而生物學家的研究受到一定的限制，未能對這個問題進行深入的數(shù)理分析，需要力學工作者的參與和合作。</p><p>　　傳統(tǒng)的空氣動力學只注意翼型在定常氣流中的升力，它不能解

9、釋昆蟲為什么能夠通過撲翅獲得足夠的升力而浮于空中：按傳統(tǒng)空氣動力學計算出來的升力遠小于昆蟲的重力【5~8】，因此產(chǎn)生了“大黃蜂飛不起來”的謬誤。隨著非定常流理論的完善，人們逐漸認識到非定常效應在昆蟲飛行升力的產(chǎn)生中起重要作用。T.Weis-Fogh的研究解開“昆蟲升力之謎”，在研究了大量的昆蟲飛行的資料以后，提出了振翅飛行的Clap－Fling機制【9】。所謂的Clap－Fling指的是：昆蟲的一對翅上抬至最高點，兩翅面保持對應翅脈平行

10、的合在一起（即Clap），然后從翅前緣（the leading-eadge）繞翅后緣(the trailing-eadge)快速的剝開。當兩翅面張開到一定程度時，兩翅面徹底分開，分別繼續(xù)作平動。出于理論分析和實驗上的考慮，這種運動引起的流動被簡化二維流動。右圖是這種運動的示意圖（fig.1）。</p><p>　　fig.1 開始階段兩翅面靠在一起（fig.1.a），然后兩翅面以某一角速度繞后緣勻速轉(zhuǎn)動（fig

11、.1.b），翅間的張角增大，空氣擠入兩翅面所張開的區(qū)域中。當兩翅面張開的角度增大到一定值時（通常認為是π/3）兩翅面脫離接觸(fig.1.c)，各自以常速度平動(fig.1.d)。（取自T.Weis-Fogh的論文）</p><p>　　在此之后，M.J.Lighthill【10】, R.H.Edwards和H.K.Cheng【11】采用無粘位勢流的假設從理論上應用復變函數(shù)論的方法，求得了此問題的理論解。實驗方面

12、，T.Maxworthy【12】 和G..R.Spedding【13】分別完成了相似的二維流動顯示實驗，其結(jié)果大致上驗證了理論上的解析解，同時也顯示了由于昆蟲飛行時的平均雷諾數(shù)（Re）處于100～2000這一范圍內(nèi)，所以空氣的粘性作用是不可忽略的。理論分析中的無粘假定導致了實驗與理論之間的差別。</p><p>　　近幾十年來，由于流體實驗技術的發(fā)展，使得對由于昆蟲振翅飛行所擾動而產(chǎn)生的流場的觀察和測量

13、更為細致。C.P.Ellington【14】采用活的吊飛昆蟲進行流場顯示實驗，其結(jié)果表明：昆蟲向下拍翅過程中，在翅的前緣將產(chǎn)生一個渦，被稱為前緣渦（the leading eadge vortex）。由于這個渦的存在，導致翅面的上方產(chǎn)生一個低壓區(qū)，因而昆蟲獲得升力。由于在實驗中采取活的昆蟲作為實驗對象使測量受到很大的限制，例如無法測定翅根部的受力情況等等。人們開始制造微型振翅機械來代替活的昆蟲作為實驗對象。這樣做法的另一個好處是可以模擬

14、空中自由飛行的昆蟲的振翅模式，從而克服測量吊飛昆蟲中所無法避免的誤差。M.H.Dickinson等人按照這一思路完成了試驗【15】，取得了很大的成功（fig.2）。</p><p>　　fig.2 升力和阻力隨時間周期性變化的曲線。紅線是實驗測量值，藍線是翅平動產(chǎn)生的升力和阻力（采用機翼理論計算）。（取自M.H.Dickison的論文）</p><p>　　由于近幾年來計算流體力學的飛速

15、發(fā)展，數(shù)值模擬逐漸成為流體力學研究的另一重要手段。Z.JangWang通過利用渦流法模擬昆蟲的振翅飛行的二維流場【16】發(fā)現(xiàn)：由二維效應產(chǎn)生的升力(fig.3)足夠使昆蟲浮于空中。</p><p>　　fig.3 采用渦流法模擬昆蟲振翅飛行所得到的升力和阻力隨時間變化的曲線圖。可以看出與前面提到的模型實驗結(jié)果十分相似。（取自Z.JangWang的論文）</p><p>　　研究昆蟲飛行

16、的根本目的是建立一套完整的振翅飛行的理論，掌握昆蟲飛行和控制的機理，并便指導微型飛行器的研究和開發(fā)。以前的研究多集中于討論振翅飛行的機理性問題，如昆蟲是如何通過振動翅膀而獲得足夠的升力，但對于控制飛行的各個參數(shù)（如振翅頻率，振翅的振幅等等）的變化對升力和阻力的影響——進而控制飛行的方向，飛行的速度，飛行的姿態(tài)——這一問題討論得較少。實驗上，設計實驗的難度很大，并且試驗要求較高的測量精度和實驗技巧，這就極大的限制了實驗手段在這一問題上的應

17、用。數(shù)值模擬具有實驗手段無法比擬的優(yōu)點，因而比較適用于這一問題。然而如果直接進行數(shù)值模擬，將遇到“動邊界問題”，這在計算流體力學中至今仍是一個棘手的問題，即使可以計算，計算的時間也將比較長，對計算機的軟硬件的要求也會比較高。本文通過轉(zhuǎn)換坐標系的方法，克服了這些問題，大大縮短了計算周期，使得仔細研究控制飛行的各個參數(shù)對飛行的影響成為可能。文中通過改變控制飛行的各個參數(shù)：振幅，平動和轉(zhuǎn)動間的相位差等等，較為詳細的研究了各種情況下升力和阻力的

18、變化，而且發(fā)現(xiàn)不同的翅尖軌跡可以造成不同的升力和阻力等等一些現(xiàn)象。</p><p>　　基本假設和簡化的振翅模型。</p><p>　　昆蟲是一個十分復雜的生物系統(tǒng)，為了研究的方便將其進行了一些簡化，從而可以采用數(shù)值模擬的方法進行分析：</p><p>　　(1). 昆蟲振翅飛行所擾動的流場本質(zhì)上是一個三維流場，翼尖渦等三維效應是存在的。然而由于昆蟲的翅本身的特征長

19、度較小，這樣三維效應有可能并不十分顯著，另外近似的分析：取通過翅縱向方向中點的橫截面，在這一截面兩側(cè)距此截面距離相等的翅表面上所受到的升力，阻力的平均值將與該中截面近似相等。因此，本文中僅研究二維的昆蟲振翅飛行。</p><p>　　(2). 雖然昆蟲的軀體不具備良好的流線型，但由于其特征尺度與其到我們所研究的翅的中點橫截面的距離相比較小，所以在本文中忽略蟲體對流場的干擾作用。</p><p&

20、gt;　　(3). 不同種類的昆蟲振翅頻率相差較大，本文中以振翅頻率處于0～300Hz的昆蟲為主要研究對象。由此可以計算出相應的流場的特征速度，對應于振翅頻率為300Hz的情形，根據(jù)已有的生物學資料可估計出特征速度遠小于100m/s因此，本文采用空氣的不可壓縮模型進行計算。</p><p>　　(4). 生物學已有的資料表明：真實昆蟲的翅表面分布有剛毛，鱗毛等毛狀器官，這些結(jié)構的存在將破壞翅表面的邊界層，但是在本

21、文中認為翅表面是光滑的。此外，由于翅的柔性使得翅在拍動過程中受到空氣的作用后將發(fā)生變形（扭曲成螺旋狀），不過在本文中這也被忽略掉了。</p><p>　　綜合以上假設，我們進行數(shù)值模擬的假定是：昆蟲的翅被簡化為光滑剛性平板，周圍的空氣采取不可壓縮模型。</p><p>　　蟲振翅模型的數(shù)學描述。</p><p>　　如下圖所示，fig.4(a)是真實昆蟲振翅飛行時的

22、翅尖軌跡曲線（空間軌跡曲線和本體軌跡曲線）。由圖上可以看出，翅尖的本體軌跡基本上是一個很扁平的曲線，這條曲線的長軸方向所在平面被稱為平均拍翅平面。Fig.4(b)是簡化的拍翅過程示意圖。</p><p>　　Fig.4 昆蟲振翅飛行模型的數(shù)學描述。Fig.4(a) 通過高速攝像取得的昆蟲振翅時翅尖的軌跡。圖中所畫的兩條軌跡，一條是相對于地面的翅尖的空間軌跡，另一條是相對于昆蟲自身的本體軌跡，軌跡上的小短線代表

23、翅面。可以看出翅尖的空間軌跡相似于正弦曲線。 Fig.4(b) 簡化的振翅模型圖。XOY是固連于地面的參考系，xoy是固連于翅面上的坐標系，msn是相對于XOY靜止的坐標系，在此坐標系下用（ζ，η）描述點的坐標。圖中實短線代表翅下拍過程中翅面的`位置，虛短線代表翅上抬過程中翅面的位置。</p><p>　　翅的運動可以分解為平動和轉(zhuǎn)動，翅的平動用翅尖的中心坐標（ζ，η）來描述，采用攻角來表征翅的轉(zhuǎn)動。綜合分析f

24、ig.4(a)并發(fā)展文獻【16】中的模型，給出如下方程描述翅的運動：</p><p><b>　　翅尖的運動方程：</b></p><p><b>　?。?）</b></p><p><b>　　（2）</b></p><p><b>　　攻角變化規(guī)律：</b&

25、gt;</p><p><b>　　（3）</b></p><p>　　其中，A0 為振幅，B0為橫向偏移量，為振翅頻率，為橫向振翅頻率，為啟動時的相位，為平動和轉(zhuǎn)動間的相位差。</p><p>　　微分（1）～（3）可以得到在msn坐標系下的速度和加速度：</p><p><b>　?。?）</b>

26、;</p><p><b>　?。?）</b></p><p><b>　?。?）</b></p><p><b>　?。?）</b></p><p>　　在msn坐標系下的角速度和加速度：</p><p><b>　?。?）</b>

27、;</p><p><b>　　（9）</b></p><p>　　由于坐標架msn與坐標架XOY相對靜止，所以上式也是在參考系XOY下的昆蟲的翅的角速度和角加速度的表達式。為了獲得在參考系XOY下昆蟲翅的速度和加速度，利用以下坐標變換公式：</p><p><b>　?。?0）</b></p><p&

28、gt;<b>　　（11）</b></p><p>　　其中，是平均拍翅平面的傾角。 是常數(shù)。</p><p><b>　　得到</b></p><p>　　相對于參考系XOY 的速度：</p><p><b>　　（12）</b></p><p>&l

29、t;b>　?。?3）</b></p><p>　　相對于參考系XOY 的加速度：</p><p><b>　?。?4）</b></p><p><b>　?。?5）</b></p><p>　　振翅飛行的昆蟲周圍流場所滿足的方程。</p><p>　　首先，

30、在固連于地面上的參考系XOY下考慮問題，</p><p>　　流體運動滿足連續(xù)性方程（不可壓縮條件）：</p><p><b>　?。?6）</b></p><p><b>　　其中，是速度場</b></p><p>　　流體運動滿足Navier-Stokes方程：</p><p

31、><b>　　（17）</b></p><p>　　其中，是空氣密度，由不可壓縮性條件，是常數(shù)。是壓力場，是動力學粘滯系數(shù)</p><p>　　初始時，假定昆蟲靜止，在t=0時刻開始起飛，因此流場滿足初條件：</p><p><b>　　（18）</b></p><p>　　其中，是無窮遠處來

32、流速度，如果昆蟲原地懸?；蛱幱诘躏w狀態(tài)，則此速度為零；如果昆蟲處于自由飛行的狀態(tài)，此速度與昆蟲的飛行速度相反。 是坐標系XOY的單位向量</p><p>　　在翅的表面上，滿足無滑移壁面條件：</p><p><b>　?。?9）</b></p><p>　　其中，C為翅面的邊界曲線， 為翅面上邊界點的瞬時速度。</p><

33、p>　　在無窮遠處，流場滿足：</p><p><b>　?。?0）</b></p><p>　　方程（16）～（20）即是在固連于地面上的慣性參考系XOY中此問題所滿足的流體運動方程。原則上，只要求解該方程組，即可解決此問題。然而，仔細觀察上述方程組可以發(fā)現(xiàn)，內(nèi)邊界（即翅面的邊界）是隨著翅的拍動而不斷變化的，這樣如果采取上述方程組直接進行數(shù)值模擬將十分困難。此

34、外，為了要達到研究不同的控制飛行參數(shù)對于昆蟲飛行的作用的目的，將導致計算所用的程序十分復雜，有可能計算時間也比較長。我們采取轉(zhuǎn)移參考系的辦法，來克服此困難。將參考系由固連于地面上的慣性參考系XOY 換成固連于昆蟲翅上的非慣性參考系xoy, 這樣動邊界問題消失了。原來的動邊界“變成”非定常的邊界條件出現(xiàn)在方程的邊條件中。采用以下變換：</p><p><b>　　（21）</b></p&

35、gt;<p><b>　?。?2）</b></p><p>　　其中，，是流體質(zhì)點在非慣性參考系xoy下的速度和加速度，，是翅尖在坐標系XOY下的速度和加速度（按照前面的假定翅是剛性不變形的，因此，翅尖的速度和加速度是已知的）。，是分別是速度和角速度。是矢徑。</p><p>　　代入上述方程組（16）～（20），得到非慣性參考系xoy下流場所滿足的方程

36、組：</p><p><b>　?。?3）</b></p><p><b>　?。?4）</b></p><p><b>　?。?5）</b></p><p><b>　　（26）</b></p><p><b>　?。?

37、7）</b></p><p>　　其中，慣性力項的表達式為：</p><p><b>　?。?8）</b></p><p>　　方程（23）～（28）就是非慣性參考系xoy下流場所滿足的方程。此方程組的分量形式為：</p><p><b>　　控制方程：</b></p>&

38、lt;p><b>　?。?9）</b></p><p><b>　　（30）</b></p><p><b>　?。?1）</b></p><p><b>　　初條件：</b></p><p><b>　　(32)</b><

39、;/p><p><b>　　(33)</b></p><p><b>　　其中，，，，</b></p><p><b>　　翅表面的邊條件：</b></p><p><b>　　（34）</b></p><p><b>　?。?/p>

40、35）</b></p><p>　　無窮遠處所滿足的邊界條件：</p><p><b>　?。?6）</b></p><p><b>　?。?7）</b></p><p>　　其中，F(xiàn)x，F(xiàn)y的形式為：</p><p><b>　?。?8）</b&

41、gt;</p><p><b>　?。?9）</b></p><p>　　實際進行數(shù)值模擬時，采用方程組（29）～（39）。</p><p>　　網(wǎng)格生成和數(shù)值求解方法。</p><p><b>　　網(wǎng)格的生成。</b></p><p>　　計算中假定翅為橢圓形長軸長為c＝1

42、cm，短軸b=c/10=0.1cm。計算區(qū)域為15倍于弦長c的圓形區(qū)域，計算區(qū)域的特征尺度遠遠大于翅的弦長，保證了在外邊界上近似地滿足上述無窮遠邊界條件（36），（37）。網(wǎng)格采用結(jié)構化網(wǎng)格，共有12960個結(jié)點，12800個控制體，最大控制體積1.004903E-4m3，最小控制體積為1.302605E-9m3。在翅表面，網(wǎng)格的特征尺度很小，小于翅表面附近的湍流邊界層的厚度?？紤]到翅尖附近流場變化較為劇烈，采用Gambit生成的網(wǎng)格如

43、下圖所示：</p><p>　　Fig.5 計算時所使用的網(wǎng)格。左圖是整個計算區(qū)域，可以看出相對于整個流場來說昆蟲的翅是很小的。因此，翅振動引起的擾動將由于空氣的粘性的耗散作用，而不能夠傳到外邊界，外邊界上的無窮遠邊界條件（36），（37）可以得到滿足。右圖是翅尖附近加密的網(wǎng)格。</p><p><b>　　數(shù)值求解方法。</b></p><p&g

44、t;　　昆蟲的翅是往復的振動，不但平動，而且同時轉(zhuǎn)動，因此，可以認為其擾動所形成的湍流是充分發(fā)展的湍流。采用標準的模型進行計算。將NS進行平均化處理得到雷諾應力方程：</p><p><b>　?。?0）</b></p><p><b>　　（41）</b></p><p>　　為了使上述方程組封閉，采用Boussines

45、q 假設：</p><p><b>　?。?2）</b></p><p><b>　　其中，k為湍流動能</b></p><p>　　是湍流粘滯系數(shù)，可以通過湍流動能,湍流耗散率計算：</p><p><b>　?。?3）</b></p><p>　　湍

46、流動能,湍流耗散率滿足的傳輸方程：</p><p><b>　　（44）</b></p><p><b>　?。?5）</b></p><p>　　方程（40）～（45）就組成了封閉的方程組，即是標準湍流模式，其中的常數(shù)經(jīng)過實驗確定為：</p><p><b>　　（46）</b&g

47、t;</p><p>　　實際計算中，采用Fluent完成計算。計算中發(fā)現(xiàn)大約經(jīng)過10～15個振翅周期升力和阻力開始呈穩(wěn)定的周期性變化，分析中著重于這一階段。</p><p><b>　　計算結(jié)果與討論。</b></p><p>　　在早期的研究中，沒有考慮到非定常效應對昆蟲獲得足夠的升力所起到的作用，因而產(chǎn)生了“大黃蜂飛不起來”的謬論。為了確

48、切的了解在昆蟲振翅飛行中非定常效應的作用，我首先采取定常的計算程序?qū)ハx的翅的空氣動力學性能進行計算（主要是升力和阻力系數(shù)）。另一方面，也可通過此計算了解計算程序的適用性和穩(wěn)定性，以及所采取的計算方法是否可行。通過計算不同攻角下翅面上（fig.6）的升力和阻力，并且采用以下公式計算得到升力和阻力系數(shù)：</p><p><b>　?。?7）</b></p><p>&l

49、t;b>　?。?8）</b></p><p>　　其中，V是來流速度，S是特征面積。</p><p>　　數(shù)值模擬的結(jié)果如下圖所示（fig.7）：</p><p>　　Fig.7 不同攻角下翅上的升力和阻力系數(shù). fig.6(a)是采用上述方法計算出來的結(jié)果；fig.6(b)是引自： M.H.Dickinson 的論文中的實驗結(jié)果。比較兩者，數(shù)值計

50、算與實驗結(jié)果十分相似，但計算結(jié)果稍小于實驗結(jié)果，可能是由于試驗和計算所采用的昆蟲翅的模型不是十分一致。</p><p>　　比較試驗和計算的結(jié)果可以發(fā)現(xiàn)：數(shù)值計算的升力和阻力系數(shù)的變化趨勢與試驗所得到的結(jié)果完全一致，但稍小于實驗的結(jié)果。這可能是由于試驗和數(shù)值計算所采用的昆蟲翅的幾何外形不是十分一致。</p><p>　　以上定常的昆蟲翅的空氣動力學性能的計算結(jié)果說明Fluent的程序有可能

51、滿足我們的需要。首先采用以下參數(shù)我們進行了計算，</p><p>　　振翅的振幅和橫向振幅分別為 ， </p><p><b>　　角度 ， ， </b></p><p><b>　　振翅的頻率為 ，</b></p><p>　　空氣的運動學粘滯系數(shù) </p><p>

52、;<b>　　空氣的密度 </b></p><p>　　計算大約10~15個周期后，升力和阻力穩(wěn)定的周期性變化，如下圖所示（fig.8）：</p><p>　　fig.8升力和阻力隨時間變化的曲線。從開始計算大約經(jīng)過10—15個周期升力和阻力開始隨時間穩(wěn)定的變化。圖中橫坐標表示周期數(shù)，縱坐標表示單位面積上作用的升力和阻力，因此，單位為N/m.</p>

53、<p>　　圖上顯示的結(jié)果和Z.Jane Wang 十分相似，表明計算程序的結(jié)果是可信的。正如前文提到的，計算是在固連于昆蟲翅上的參考系下進行的，所以不需要變動網(wǎng)格，僅僅邊界條件是隨時間變化的。然而在非慣性參考系下許多物理量，以及他們之間的關系是不同于慣性系下的關系，但是可以驗證在慣性參考系XOY下的渦量和在非慣性參考系xoy下的渦量僅僅相差一個翅的轉(zhuǎn)動角速度渦量，也就是：</p><p><b&

54、gt;　　（49）</b></p><p>　　計算所得的渦場隨時間變化的過程如下圖（fig.9）所示。雖然不同時間的渦場圖之間沒有可比性，因為不同時刻翅的轉(zhuǎn)動角速度不同，然而從圖上可以看出同一時刻渦量的分布。由圖上可以明顯看出，在平動速度最大時，即t=0 T，t=0.5 T時，翅兩端的渦量方向相反；在轉(zhuǎn)動角速度最大時，即t=0.25 T，t=0.75 T時，翅兩端的渦量方向相同。</p>

55、<p>　　Fig.9昆蟲振翅飛行時周圍空氣的渦場隨時間變化的過程（在非慣性參考系xoy下觀察）。慣性參考系XOY下的渦量與非慣性參考系下的渦量僅相差一個常數(shù)：翅的轉(zhuǎn)動角速度，而且這一關系在任一時刻對全流場的每一空間點都成立。左一列的兩張圖表示轉(zhuǎn)動的角速度最大時的渦場圖，在翅的兩端附近渦量具有相反的方向；右一列的兩張圖表示平動速度最大時的渦場圖，在翅的兩端附近渦量具有相同的方向</p><p>　　

56、正如在前言中提到的，在研究昆蟲飛行的早期，人們并未認識到非定常效應對昆蟲獲得升力的作用，為了確切分析非定常效應這一作用我們定義如下兩個量：</p><p><b>　　定常的升力：</b></p><p><b>　　（50）</b></p><p><b>　　定常的阻力</b></p>

57、<p><b>　　（51）</b></p><p>　　其中，是升力和阻力系數(shù)，可以采用上一節(jié)計算出來的值。是平動速度，由于翅的運動已知這一值可由（12）（13）式計算而得。下圖表示這樣處理的結(jié)果（fig.10）。</p><p>　　fig.10 一個多周期內(nèi)升力隨時間變化的曲線圖。圖中橫坐標軸代表周期數(shù)，縱坐標軸代表單位截面上升力的大小。其中，紅

58、線是數(shù)值計算的結(jié)果（定常效應＋非定常效應），藍線是按上述定義計算的定常效應產(chǎn)生的升力，黃線是兩者之差，即為非定常效應。非定常效應對阻力也有較大的影響。</p><p>　　從圖中可以明顯看出，非定常效應使升力的峰值明顯增大，通過實際的計算發(fā)現(xiàn)，非定常效應還使得平均升力增大20.07％，這樣大的差別足以解釋為什么采用定常流理論無法解釋昆蟲升力之謎。類似于升力，昆蟲翅上所受到的阻力也受到非定常效應的影響，雖然相對較小

59、一些（8.21%）。</p><p>　　在實際設計微型飛行器的工程中，關心的另外一個問題就是振翅飛行過程中所消耗的功率。由于昆蟲的翅的質(zhì)量相對于昆蟲的軀體的質(zhì)量來說是一個小量，所以以下在忽略昆蟲翅的質(zhì)量的假設下，討論飛行中功率的分配。昆蟲飛行過程中消耗的功率由三部分組成：產(chǎn)生升力消耗的功率，產(chǎn)生推力消耗的功率，克服扭矩轉(zhuǎn)動而消耗的功率。以上三部分功率可以表達成如下的形式：</p><p>

60、;<b>　　（52）</b></p><p><b>　?。?3）</b></p><p><b>　?。?4）</b></p><p>　　其中，分別是翅平動速度的X軸，Y軸分量，轉(zhuǎn)動角速度。</p><p>　　這樣，昆蟲飛行時消耗的總工率就可以表達為三者之和：</

61、p><p><b>　?。?5）</b></p><p>　　按公式（51）～（54），（55）計算三部分的即時功率，所得功率隨時間變化圖（fig.11）如下：</p><p>　　fig.11 飛行過程中消耗的即時功率隨時間變化的過程。圖中，橫坐標代表時間單位是妙；縱坐標代表翅的單位截面上消耗的功率。紅線代表升力消耗的功率，藍線代表推力消耗的功

62、率，黃線代表翅轉(zhuǎn)動所消耗的功率，淺紅色的曲線代表飛行過程中消耗的總工率。</p><p>　　由圖上可以清楚看到前半個振翅周期（即翅由最高點向下運動的過程）主要是升力在消耗功率，而推力消耗的功率很小；在后半個振翅周期（即翅由最低點向上運動的過程）主要是推力在消耗功率，而升力所消耗的功率很少。翅的轉(zhuǎn)動所消耗的功率在整個過程中都很小。對這些曲線在一個周期內(nèi)積分，可以得到平均功率：</p><p&g

63、t;<b>　?。?6）</b></p><p><b>　?。?7）</b></p><p><b>　　（58）</b></p><p><b>　?。?9）</b></p><p>　　在振翅過程中，推力和升力消耗的功率基本上分別消耗總功率的50%，

64、而轉(zhuǎn)動消耗的功率實際上可以忽略。對于飛行器設計來說這樣的功率分配是十分不合理的。因此，從這個意義上來說振翅飛行并不是一種有效的飛行方式。</p><p>　　正如在引言中提到的，計算是在固連于翅面上的參考系中進行的，所以不需要在計算過程中更改網(wǎng)格，動邊界條件轉(zhuǎn)化為隨時間變化的邊界條件，而且保持邊界不動。這樣對于不同的控制飛行的參數(shù)，僅僅需要改變?nèi)舾蓞?shù)即可實現(xiàn)對于不同條件的模擬。利用這一方法研究了若干控制飛行的參

65、數(shù)。</p><p>　　昆蟲的振翅過程可以看作翅的平動和轉(zhuǎn)動的疊加，因為平動最大速度發(fā)生于翅面處于水平位置，而對于的情形最大轉(zhuǎn)動角速度發(fā)生于上拍的最高點和下拍的最低點，這樣振幅其實控制了平動和轉(zhuǎn)動之間的空間距離。此外，顯然另一控制平動和轉(zhuǎn)動的參數(shù)就是平動和轉(zhuǎn)動間的相位差，它控制了平動和轉(zhuǎn)動之間的時間距離。增大振幅，升力曲線的最大值顯著提高，而在其余的位置上變化不大。利用下式（59）對平均升力進行無量綱化，<

66、;/p><p><b>　　（60）</b></p><p>　　這樣就得到了無量綱化的平均升力于振幅之間的關系（fig.12）：</p><p>　　fig.12 無量綱化的平均升力與振幅之間的關系。平均升力采用（59）式進行無量綱化，橫軸表示振幅，單位是cm。從圖上明顯可以看出當振幅增大到一定程度，無量綱化的平均升力趨于一個基本固定的常數(shù)（0

67、.55左右）。</p><p>　　無量綱化的平均升力隨著的振幅增大而增大，當振幅增大到一定值時，無量綱化的平均速度趨于一近似恒定的值（0.55）。這是因為，振幅增大使得轉(zhuǎn)動和平動的最大值發(fā)生的位置的空間距離加大，兩者所擾動的流場相互的干擾減小。另一方面，平動和轉(zhuǎn)動間的相位差也表現(xiàn)出相似的作用。下圖（fig.13）表示不同的相位差所對應的升力曲線（隨時間變化）。</p><p><b

68、>　　f</b></p><p>　　fig.13 不同平動和轉(zhuǎn)動間的相位差下升力隨時間變化的曲線，其中，相鄰的曲線間相差。圖中，時間軸是橫軸，縱軸表示單位長度的翅的截面上所受到的升力。</p><p>　　明顯可以看到，除了在幾個曲線的極大值和極小值的位置上，升力隨相位差的不同而發(fā)生變化外，在其他的時間上各條曲線基本上是重合在一起的。這一現(xiàn)象是由于相位差改變了最大轉(zhuǎn)動

69、速度發(fā)生的位置（最大平動速度發(fā)生的位置不隨相位差變化，參看（1），（2），（3）式）,通過平動和轉(zhuǎn)動所擾動的流場的作用，升力和阻力發(fā)生變化。通過積分一個周期內(nèi)的升力值，得到平均升力隨平動與轉(zhuǎn)動間相位差變化的曲線圖（fig.14）</p><p>　　fig.14 平均升力雖平動和轉(zhuǎn)動間相位差變化的曲線圖。圖中，橫軸表示相位差，縱軸表示單位長度上的升力。圖中，紅色的叉子表示計算所得到的平均升力，相鄰兩個計算值得相

70、位差相差。當相位差減小時，平均升力也隨之減??；反之平均升力隨相位差的增大而增大。</p><p>　　圖上每一個紅色的叉子代表一個計算值，相鄰兩點之間相位差相差。相位差增大時，平均升力也隨之增大，反之，平均升力減小。相位差發(fā)生較小的變化就可以引起升力和阻力的較大變化，這一現(xiàn)象在試驗中已經(jīng)被發(fā)現(xiàn)，并被稱為“升力和阻力對相位差的敏感性”。這一現(xiàn)象的發(fā)現(xiàn)為解釋昆蟲如何在飛行中十分靈活的控制飛行的方向，飛行的姿態(tài)提供了依

71、據(jù)。一種簡單的假設是：昆蟲在飛行過程中可以通過控制左右兩個翅的振動的平動和轉(zhuǎn)動的相位差，使左右兩個相位差不同，從而使作用于左右兩個翅上的升力和阻力不同，這樣就可以形成偏航力矩，俯仰力矩等等。由于控制這種相位差相對于控制振翅的振幅來說，需要的消耗能量少的多，而且敏感性比較強，所以在微型飛行體的設計中是一個可以考慮的控制飛行的方案。</p><p>　　在第二節(jié)的“基本假設”中，昆蟲的翅被簡化為具有光滑表面的平板，而

72、忽略了翅表面上的各種微觀結(jié)構對周圍流場可能的影響?？梢酝ㄟ^另外一個途徑，近似地了解這種微觀結(jié)構的作用。由于微觀結(jié)構主要是影響翅面附近流場的邊界層，并且主要是通過于流體的粘性摩擦作用耗散流場的動能，而在邊界層中粘性力和慣性力相當，粘性力與流體的粘滯系數(shù)成正比，因此微觀結(jié)構對邊界層的這種作用可以通過改變流體的粘滯系數(shù)加以模擬。流體具有較小的粘滯系數(shù)就相當于昆蟲的翅表面比較光滑；反之，具有較大粘滯系數(shù)的流體的流場代表昆蟲的翅表面各種微觀結(jié)構比

73、較多。利用數(shù)值計算的方法可以比較容易的實現(xiàn)各種情況，下圖（fig.15）就是計算的結(jié)果：</p><p>　　fig.15 平均升力隨空氣的動力學粘滯系數(shù)變化的曲線，紅色的點是計算的平均升力的值，相鄰的點對應的空氣的動力學粘滯系數(shù)相差一個量級。值得注意的是，橫坐標是空氣的動力學粘滯系數(shù)，此軸為對數(shù)軸。在空氣的動力學粘滯系數(shù)處于一定范圍內(nèi)時，平均升力變化很小。</p><p>　　由圖上可

74、以看出，當空氣的動力學粘性系數(shù)在正?？諝獾膭恿W粘滯系數(shù)左右變化時平均升力變化不大。當空氣的動力學粘滯系數(shù)增大到的量級上時，平均升力急劇下降，說明如果昆蟲翅上的微觀結(jié)構的作用還是比較大的。</p><p>　　前面的討論都是在假定昆蟲是靜止不動（相當于昆蟲處于吊飛狀態(tài)）的前提下進行的。這種假定情形與處于飛行的昆蟲的情形有一個十分明顯的力學上的差別：處于水平飛行的昆蟲，從固連于昆蟲上的參考系來看，其無窮遠處的邊界條

75、件上來流速度不為零?？紤]到這種差別以后，翅上所獲得的升力和阻力有相應的變化，實際計算的結(jié)果如下圖所示（fig.16）：（假定昆蟲處于水平勻速飛行狀態(tài)，無窮遠來流速度與飛行速度相反，參看（20），（36），（37）式）</p><p>　　fig.16 平均升力隨速度（用2m/s無量綱化）變化的曲線。</p><p>　　除了上面提到的無窮遠邊界條件不同外，吊飛昆蟲和自由飛行的昆蟲在其他方

76、面也存在著差異。根據(jù)已有的生物學資料：處于吊飛狀態(tài)下的昆蟲在振翅過程中其翅尖軌跡是近似于橢圓的一條十分扁平的曲線；而處于自由飛行的昆蟲的翅尖軌跡是“八字形”或“雙八字形”；或是更為復雜的曲線，如圖所示（fig.17）。</p><p>　　fig.17 昆蟲翅尖軌跡圖。(a)橢圓形軌跡代表昆蟲處于吊飛狀態(tài)下，(b),(c)“八字形”軌跡和“雙八字形軌跡”對應于處于自由飛行的昆蟲振翅模式。軌跡上的小箭頭表示翅尖的

77、運動方向。</p><p>　　以上的討論都是假定翅尖沿某一傾斜的直線往復運動，即橫向振幅。關于這一問題的研究似乎未引起太多的注意，但計算的結(jié)果顯示翅尖軌跡對昆蟲振翅所獲得升力和阻力有較大的影響。在前面提到的昆蟲振翅的數(shù)學模型（1），（2），（3）中，取不同的的值就可以實現(xiàn)不同的翅尖軌跡。具體計算中，取該比值為1，2，3分別模擬了橢圓形軌跡，八字形軌跡，雙八字形軌跡，橫向振幅，，。計算結(jié)果如下圖所示（fig.18

78、）: </p><p>　　fig.18(a) 橢圓形翅尖軌跡時，升力隨時間變化的過程。</p><p>　　fig.18(b) “八字形”翅尖軌跡時，升力隨時間變化的過程。</p><p>　　fig.18(c) “雙八字形”翅尖軌跡時，升力隨時間變化的過程。</p><p>　　Fig.8 不同的軌跡下，振翅飛行的昆蟲所獲得的瞬時升

79、力隨時間變化的過程。fig.18(a) 橢圓形翅尖軌跡時，升力隨時間變化的過程; fig.18(b) “八字形”翅尖軌跡時，升力隨時間變化的過程。fig.18(c) “雙八字形”翅尖軌跡時，升力隨時間變化的過程。橫軸表示時間，縱軸表示瞬時升力。明顯可以看出，在(a)中三條曲線基本上重合在一起，而在(b),(c)中在極大值的位置上三條曲線明顯分開，尤其是在(c)中這一現(xiàn)象更加明顯。</p><p>　　由圖上

80、可以看到，當翅尖軌跡是橢圓形時，增大橫向振幅，升力曲線基本上沒有變化。但當翅尖軌跡是“八字形”時，增大橫向振幅，升力曲線在原極大值的位置上繼續(xù)增大，當翅尖軌跡是“雙八字形”，這一效應更加明顯了。同樣地，對瞬時的升力和阻力在一個振翅周期內(nèi)積分，可以得到平均的升力和阻力，這一計算結(jié)果如下表(tab.1)所示：</p><p>　　tab.1 不同翅尖軌跡對應的平均升力和平均推力值，箭頭表示橫向增大時，平均推力或平均

81、升力變化的方向。</p><p>　　雖然目前為止還找不到一個合適的方式解釋這一問題，但可以提供一個可能解釋此問題的想法：仔細觀察不同的翅尖軌跡，可以發(fā)現(xiàn)在橢圓形軌跡，“八字形”軌跡，“雙八字形”軌跡，三者中橢圓形軌跡所圍的區(qū)域最大，“雙八字形”軌跡所圍的區(qū)域（三個小塊區(qū)域）最小。因此，橢圓形軌跡對流體的“限制作用”最小，“八字形軌跡”其次，“雙八字形”軌跡對流體的“限制作用”最強，因為它圍成的區(qū)域最小。“限制作

82、用”越強所獲得的升力的峰值也就越大，但同時升力變化得也就越劇烈。這一解釋方式明顯是不十分嚴謹?shù)?，需要進一步的研究。然而這一現(xiàn)象對于微型飛行器的設計卻有比較重要的意義：一方面，設計中可以利用這一方式提高可以獲得的升力大小，還可以作為一種產(chǎn)生控制力，控制力矩（俯仰力矩，偏航力矩等等）的方案。但另一方面，必須防止由這一現(xiàn)象產(chǎn)生的升力或推力劇烈變化，導致飛行器喪失飛行的穩(wěn)定性。</p><p>　　利用這一振翅模型和計算

83、方法，我們還計算了起飛段的情形。簡單的，起飛段的模型可以從（1），（2），（3）式稍加變化的來。考慮下式：</p><p><b>　?。?1）</b></p><p>　　其中，是穩(wěn)定時的振幅，n控制達到穩(wěn)定所需要的時間，T為振翅周期。</p><p>　　將（60）式帶入（1），（2），（3）式即可得到起飛段昆蟲振翅過程的數(shù)學描述。采用同樣

84、的計算方法可以得到升力隨時間的變化如下圖所示：</p><p>　　fig.19 起飛段升力隨時間的變化過程。振翅模式按照（61），（1），（2），（3）決定。</p><p><b>　　七．結(jié)論。</b></p><p>　　以上的計算結(jié)果表明轉(zhuǎn)換參考系處理這一問題的方法是可行的，所得到的結(jié)果與已有的試驗的結(jié)果可以比較好的符合。雖然從理論

85、上可以有許多參數(shù)（如振幅，頻率等等等，頻率，最大攻角的影響比較復雜，但我們也做了相關的數(shù)值計算，但需要進行進一步的研究。）可以控制飛行得速度，方向，飛行的姿態(tài)等，但是對于生物體——昆蟲來說，控制某些參數(shù)（翅的幾何外形等）是不可能的，或者控制某些參數(shù)不是“十分經(jīng)濟的”。這樣根據(jù)以上的計算結(jié)果，可以確定一些參數(shù)可能被昆蟲所使用，來控制飛行。例如平動和轉(zhuǎn)動間的相位差，以及翅尖軌跡等等。最后，值得一提的是，三維情形下的振翅的數(shù)值模擬相對于二維有

86、更多的困難。轉(zhuǎn)換參考系的方法原則上可以同樣地應用于三維情形，但是此時就必須考慮昆蟲的軀體對流場的影響，因此可能需要其他的方法加以克服。</p><p><b>　　致謝</b></p><p>　　首先，我要感謝李政道和他夫人提供這一基金，使得我有機會與處于科學研究第一線的學者接觸，直接聆聽他們的教誨，學習科學研究的方法。</p><p>　　

87、本文是在陳濱教授的指導下完成的，在論文工作的各個階段，陳濱教授總是給我很多的啟示。正是在陳教授的指導，鼓勵，幫助下，這篇論文才能得以完成。在我完成“君政”學者計劃的兩年期間，陳教授不僅在專業(yè)上引導我進入學術研究的殿堂，教給我科學的研究方法，而且更重要的是交給我做人的道理。陳濱教授的淵博的學術知識，敏銳的抓住問題核心的能力，和對力學發(fā)展方向的洞察力都是我學習的榜樣。陳濱教授兩年間對我的學習，生活上的指導，幫助將使我終身受益。</p&

88、gt;<p>　　在我的后一年的研究時間里，我到陳耀松教授的實驗室，完成數(shù)值計算的工作。在此期間，陳耀松教授不但無私的提供給我必要的設備外，還熱忱地對我的科研工作給于了具體的指導和幫助。對我的具體的研究方向和論文的撰寫進行了細致的關心。陳耀松教授淵博的學術知識，敏銳的抓住問題核心的能力，簡單清晰的物理建立模型的能力，“活到老學到老”的學習精神和對力學發(fā)展方向的洞察力都是我學習的榜樣。陳耀松教授對力學無盡的癡迷之情，對事業(yè)的

89、執(zhí)著追求，對他人無私的幫助，以及堂堂正正的生活態(tài)度給我以深刻的印象，值得我永遠回味和學習。</p><p>　　系里很多教授同樣給我以熱忱的幫助。與魏慶鼎教授，陸來教授，翟大熙教授，陳十一教授，陳國謙教授，吳介之教授，李植老師，譚文長教授，蘇衛(wèi)東教授，蔡慶東教授的討論使我受益匪淺。特別是，李植老師一直關心我的研究進展，并且提供了許多好的建議。他還對我的文章進行了十分認真的修改，達到字斟句酌的程度，他的認真的態(tài)度令

90、我感動。陳國謙教授，陳十一教授，吳介之教授在我的研究過程中一直關心我的研究進展情況，他們一直積極的引導我將研究推向更加深入的境界。</p><p>　　同時，實驗室里濃厚的學術氣氛，輕松開放的學術氛圍，都是我得以完成這一論文的必要條件。在這兩年期間，安亦然給我以十分巨大的幫助，很多問題是在與他進行反復討論的基礎上完成的，在與他的討論中我不僅可以了解到許多力學發(fā)展的前沿，同時我也看到了自己的差距。此外，他肯于鉆研的

91、精神值得我學習。此外，實驗室里的其他同學，張立柱，李勇，李向群，李正紅，陳彬，都給我以很大的地幫助。我的很多想法都是在于他們的熱烈的討論中獲得的，十分感謝他們的啟發(fā)。從張立柱同學那里我學到了腳踏實地的研究作風，從李向群同學那我了解到許多有用的計算機，計算機網(wǎng)絡的知識，等等。在此我向他們以及其他所有關心和幫助過我的人表示衷心的感謝。</p><p>　　最后，我要感謝我的父母，沒有他們的理解和鼓勵，我是不可能克服研

92、究過程中所遇到的各種困難，最終完成這篇論文。</p><p><b>　　參考文獻</b></p><p>　　【1】R.J.Wooton, Annu.Rev.Entomol, 37:113-40(1992)</p><p>　　【2】A.F.G..Dixon and P.Kindlmann, Ecology, 80(5), pp.1670-1

93、690(1998)</p><p>　　【3】J.H.Marden, J.Exp.Biol, 130, 235-258(1987)</p><p>　　【4】G.Ruppell, J.Exp.Biol, 144, 13-42(1989)</p><p>　　【5】T.Weis-Fogh and M.Jensen, Proc.R.Soc.London B, 239,

94、415(1956)</p><p>　　【6】G.R.Spedding, Contemp.Math, 141, 401(1993)</p><p>　　【7】R.Dudley, The Biomechanics of Insect Flight, Priceton University Press , Princeton 2000</p><p>　　【8】C.P.

95、Ellington, Philos.Trans.R.Soc.Lodon B, 305, 1(1984)</p><p>　　【9】T.Weis-Fogh, J.Exp.Biol, 59, 169-230(1973)</p><p>　　【10】M.J.Lighthill, J.Fluid Mech.vol.60, part.1 pp.1-17(1973)</p><p&

96、gt;　　【11】R.H.Edwards and H.K.Cheng, J.Fluid Mech, vol.120, pp.463-473(1982)</p><p>　　【12】T.Maxworthy, J.Fluid Mech, vol,93, part.1, pp.47-63(1979)</p><p>　　【13】G.R.Spedding and T.Maxworthy, J.Fl

97、uid Mech, 165, 247(1986)</p><p>　　【14】C.P.Ellington, C.V.D.Berg, A.P.Willmott and A.L.R.Thomas, Nature, 384, 626(1996)</p><p>　　【15】M.H.Dickinson, F.O.Lehmann and S.P.Sane, Sci, 284, 1954(1999)&

98、lt;/p><p>　　【16】Z.Jane Wang, Physical Review Letters, vol.85, num.10, 2216-2219(2000)</p><p>　　作者簡介：程暮林，男，1979年9月22日生于大連。高中畢業(yè)于大連市第一高中，1998年9月進入北京大學力學與工程科學系理論與應用力學專業(yè)學習。1999年師從陳濱教授開始“君政”專題研究計劃。在大學期間，