S1223翼型俯仰-沉浮運動的非定常氣動特性分析

常興華, 馬 戎, 張來平,*, 赫 新

(1. 中國空氣動力研究與發展中心 空氣動力學國家重點實驗室, 四川 綿陽 621000；2. 中國空氣動力研究與發展中心 計算空氣動力研究所, 四川 綿陽 621000)

S1223翼型俯仰-沉浮運動的非定常氣動特性分析

常興華1,2, 馬 戎2, 張來平1,2,*, 赫 新2

(1. 中國空氣動力研究與發展中心 空氣動力學國家重點實驗室, 四川 綿陽 621000；2. 中國空氣動力研究與發展中心 計算空氣動力研究所, 四川 綿陽 621000)

撲翼飛行器是未來航空領域的重要發展方向之一，而鳥類、昆蟲等自然界的飛行生物所具有的出色飛行能力，為人造撲翼飛行器的設計工作提供了很好的參照。本文以鳥類的撲翼飛行為研究背景，針對簡化的二維S1223翼型的剛性俯仰-沉浮運動，對其非定常氣動力特性進行了數值研究。研究采用了動態混合網格技術以及非定常數值計算方法。為了提高動態混合網格的變形能力，采用了基于徑向基函數的插值方法求解空間點的位移，翼型整個俯仰-沉浮運動周期內計算網格均維持了較好的質量，沒有發生網格重構。非定常算法方面，通過約束單元邊界面的法向速度從而滿足了運動網格下幾何守恒律的要求。空間離散采用了二階的有限體積格式，時間離散則采用了雙時間步和BLU-SGS相結合的隱式時間推進策略。計算得到了不同下拍時間、不同拍動角等條件下的升力、推力以及能耗，對其升力、推力產生機制進行了分析，并通過對氣動力以及流場進行對比，分析了各拍動參數的影響。計算結果表明，翼型自身的“靜態因素”是其產生升力的主要原因，非定常流動對增加升力起到了促進作用，而下拍時間、拍動角等運動參數對翼型的氣動性能影響較大。當下拍時間占到整個拍動周期的約65%-70%時，單位能耗下的時均升力最大，該結論和觀測數據較為一致。此外，通過對比分析得到了一組具有較好氣動特性的拍動角參數，為后續針對三維問題的研究提供了參考。

鳥類撲翼；動態混合網格；非定常流動；RBF動網格方法；幾何守恒律

0 引 言

自然界的鳥類經過億萬年的進化，形成了出色的飛行能力，其飛行效率、機動性都遠遠超出了目前所有的人造飛行器。具有較低能耗、噪聲，同時具有優越機動性能的小型撲翼飛行器是未來航空器發展的一個重要方向，為開展此類研制工作，首先需要對鳥類撲翼運動過程中的升力、推力產生機制以及流動控制機理進行深入細致的研究。

鳥類翅膀的構造非常精細，由肌肉、骨骼、羽毛和多個關節組成，撲翼過程中伴隨著翅膀的弦向轉動、展向的收縮-伸展以及扭轉[1]，其飛羽在上拍、下拍過程也會規律性的收縮-合攏。鳥類的升力主要源于翼型本身以及翼型的非對稱拍動，例如海鷗的翼型剖面和S1223翼型非常相似，在低速情況下具有很高的升阻比，而鳥翼在下拍過程中會完全打開，上拍過程中則會盡量收縮，且其上拍、下拍的時間也存在很大差異，這種非對稱拍動的目的也是為了增強升力。

目前國內外常見的關于鳥類撲翼飛行的研究工作大概可以分為兩類：一類通過觀測，采用數據統計和量綱分析的手段，分析鳥類形態和各飛行參數之間的關系，并得到一些經驗性的公式，如Greenewalt[2]、Tennekes[3]、Pennycuick[4]等人的工作。另一類則是采用試驗、理論或近似的數值分析方法，針對具體的撲翼過程開展研究，如Jones[5]、Angela[6]、Pennycuick[7]、Rayner[8]、Phlips[9]等人的工作。其中“面元法”作為一種簡單有效的數值模擬手段，在鳥類撲翼的研究中使用較多，如Smith[10]采用該方法并結合有限元模型對柔性撲翼運動進行了研究，國內的昂海松研究團隊[11]、余永亮[12]等也分別采用該方法對昆蟲及蝙蝠的撲翼進行了研究。

雖然這些針對鳥類撲翼運動的研究已經初步揭示了鳥類的升力、推力產生機制，然而這些定性的認識仍不能滿足撲翼飛行器設計上的需求。鳥類撲翼運動的許多細節問題如翼型、迎角、展向變形/扭轉、上/下拍動時間比、拍動頻率等對其力學特性有非常重要的影響，因此非常有必要針對這些細節問題進行深入研究。

計算流體力學的迅速發展為撲翼運動的精細化研究提供了條件，且已經在昆蟲的撲翼研究中得到了廣泛應用[13-16]。然而由于鳥類撲翼過程的復雜性，三維情況下復雜撲翼過程的精細化數值研究還不多見，很多數值研究工作仍基于二維或者簡單的三維撲翼開展[17-18]。

本文以海鷗的撲翼飛行為研究背景，選擇二維的S1223翼型作為研究對象，采用動態混合網格技術及非定常數值模擬方法，對其剛性俯仰-沉浮運動進行了研究，得到了下拍時間比、拍動迎角等參數的影響規律。這些工作為后續的復雜三維問題的數值研究奠定了基礎。

1 數值計算方法

本文的數值模擬基于HyperFLOW軟件平臺[19]開展。該軟件平臺是中國空氣動力研究與發展中心研發的具有完全自主知識產權的大型CFD多學科通用求解平臺，具有優越的體系架構和生態系統，并已集成了結構/非結構NS方程流場解算器、動態混合網格生成技術、飛行力學/流體動力學一體化算法等，可進行完全氣體和化學非平衡氣體的定常/非定常計算。以下對其中的動態混合網格技術以及非定常算法進行簡要介紹。

1.1 動態混合網格生成技術

常用的網格變形技術有彈簧松弛法、求解偏微分方程的方法以及插值方法等。在之前的研究工作中，作者所在的課題組建立了彈簧松弛法和插值法相結合的混合網格變形技術，具有較好的變形能力和動網格生成效率，并通過結合局部網格重構技術，提高了針對大變形、大位移、相對運動等復雜動邊界問題的適應能力[20-22]。

在最近的研究工作中，我們進一步發展了基于徑向基函數(RBF)的網格變形技術，該技術具有優越的網格變形能力，是目前國內外CFD工作者爭相研究的一個熱點。而標準的RBF方法在處理大規模網格時效率極差，為了提高RBF網格變形技術的適用性，我們參考文獻[23]的做法，通過選擇有限的參考點來減少RBF算法中矩陣的規模，以提高計算效率。

1.2 非定常數值模擬方法

基于動態混合網格的NS方程解算器采用了格心型的有限體積格式，時間離散采用二階的歐拉后插方法，為提高非定常計算效率，采用了雙時間步算法和BLU-SGS隱式計算方法。算法具體細節請參見文獻[22]。

動網格下的NS方程離散算法必須滿足幾何守恒律。運動網格下的幾何守恒律(GCL)特指體積守恒，對于控制方程，假設流場均勻，則可得到GCL方程如下：

其中V為控制體Ω的體積，?Ω為控制體邊界，w表示控制體外邊界的運動速度，n為控制體邊界的外法向。導數符號“d”表示跟隨控制體的隨體導數，而不是跟隨流體的物質導數。

式(1)表示控制體單元的體積變化率等于邊界面法向運動速度的面積分，其在微分意義下是恒成立的，而計算所采用的離散格式也需要使式(1)得到滿足，才不至于在流場內引入額外的誤差。

滿足式(1)的數值格式很多，詳見Farhat[24]、Kallinderis[25]、劉君[26]、Wang[27]等的工作。這里我們采用了一種較為便捷的方法：直接對邊界面的法向速度進行約束。以一階歐拉隱式方法為例，式(1)的離散形式為：

式中v=w·n表示控制體邊界面的法向速度。對于單元的每一個離散的邊界面，限定其法向速度的求解方法：

在隨后的研究中，我們進一步從誤差分析的角度對各種運動網格幾何守恒算法進行了理論分析，發現可以將現有的幾何守恒算法歸為兩類，即限制整體積分誤差的“體限制方法”和限制每個邊數值誤差的“面限制方法”，并分別可以針對不同的時間格式寫成統一的形式。進一步的對比發現，各種幾何守恒算法的計算結果大體一致，但是我們發展的幾何守恒算法在三維情況下更為簡便。具體的理論分析詳見文獻[28]。

2 翼型及俯仰-沉浮運動

根據文獻描述[29]，本文選用二維S1223翼型作為研究對象，翼型弦長c=0.2 m。將整個拍動周期分為下拍和上拍兩個階段，其沉浮運動規律如下：

其中dy為翼型的縱向位移，A為振幅，T為拍動周期，k表示下拍時間和整個拍動周期T的比值。本文令振幅A=0.2 m。

翼型繞其1/4弦長處進行俯仰，俯仰角(低頭為正)變化規律如下：

其中θ0表示上拍和下拍過程的起始拍動角，φ表示俯仰角變化與上拍-下拍運動之間的相位差。本文令φ=0、θ0=0，因此俯仰角振幅θdown、θup分別代表最大下拍角和最大上拍角，且最大拍動角出現在上拍或下拍過程的中間時刻。

本文只針對二維問題進行數值模擬，其相當于三維問題的一個展向截面(如圖1所示)。

本文需要研究俯仰-沉浮運動的能耗特性，能耗系數CP的定義如下：

其中f表示流體作用在物體邊界上的力，w表示邊界元的運動速度，?Ω表示物體表面。

3 數值結果分析

采用混合網格離散計算域，物面附近采用各向異性的四邊形單元。翼型外圍采用三角形網格單元，并在翼面附近以及尾跡內進行了加密處理，網格單元總數8.6萬(圖2)。計算采用SA湍流模型，不考慮轉捩因素的影響，采用二階空間和時間格式，一個拍動周期采用400個非定常時間步離散。

3.1 下拍時間比對升力及能耗的影響

觀察發現[30]，海鷗等鳥類巡航飛行時其下拍時間占到整個拍動周期的60%～80%，說明下拍時間比k與其力學特性關系密切。本節不考慮翼型的俯仰，僅考慮其上下沉浮運動，針對k=0.40～0.75等幾種情況進行數值模擬。

圖3(a)所示為不同k值情況下升力系數在一個周期內的變化曲線，橫坐標為拍動的相位角。下拍過程中(Phase<π)，在翼型的隨體坐標系下來流的等效迎角為正，因此產生正升力，但是隨著k值的增加，翼型下拍的速度減小，等效迎角減小，因此升力峰也隨之減弱；上拍過程中，翼型等效迎角為負，因此翼型產生負升力，且隨著k的增加，上拍速度增加，負升力峰增大。但是由于翼型本身具有較好的升阻比特性，在上拍和下拍的起始階段(Phase=0、π附近)仍然可以維持較大的升力系數，這對于增大整個拍動周期的平均升力是有益的。俯仰運動的瞬時能耗系數與翼型拍動速度的大小關系密切，在拍動速度最快的時刻(Phase=0.5π、1.5π)會出現能耗的峰值(如圖3b所示)，隨著k值的增加，下拍過程的能耗峰值減弱，而上拍過程能耗峰值則會增加。

圖4所示為拍動一個周期內幾個典型相位的渦量云圖，在Phase=0.5π、1.5π兩個相位下，不同的k值對應的流場結構基本類似，只是由于k值的變化影響了其下拍和上拍過程的等效迎角，因此背風區內的旋渦大小、形態等存在一些差異；而Phase=0、π兩個相位下，翼型的法向運動速度為0，但是加速度處在整個拍動周期的峰值，此時不同k值條件下流場的旋渦結構存在較大差異。以Phase=0為例，此時處于上拍結束的急劇減速以及下拍起始的急劇加速階段，翼型帶動周圍流體急劇加速/減速，受“附加質量效應”影響較大，且由于本文設計的運動規律在此時的加速度不連續，因此導致升力系數存在“不規律”的波動(圖3c)。然而，由于翼型拍動的減縮頻率較小，來流速度相對于拍動速度是一個較大的量，因此這些“不規則”的旋渦結構會很快消失在尾跡，升力系數也會很快進入單調增加的趨勢。

圖3 不同k值下一個拍動周期內的氣動力系數
Fig.3 Aerodynamic coefficients in a period at different value of k

上述的數值模擬結果表明，翼型的下拍過程對升力的貢獻較大，而上拍過程則產生負升力。實際上，鳥類在撲翼過程中正是通過翅膀展向的折轉來突出下拍的“有益”作用(下拍時翼面全部展開、增大有效升力面)，減弱上拍的“不利”影響(上拍時適當折疊撲翼，減小上拍迎風阻力面)。為了近似反應這種“三維效應”，本文對升力、阻力、能耗等引入一個近似的加權系數：

λ在0.5和1.0之間周期性變化，φ表示加權系數和撲翼運動之間的相位差，本文令φ=0。將計算得到的瞬時氣動力系數乘上加權系數λ，則本文的二維問題可近似看作一個三維的翼型在做俯仰-沉浮運動的同時改變其展向長度。

考慮加權前后，撲翼運動的時均升力系數、時均能耗系數以及它們之間的比值隨k值的變化如圖5所示：與原始模式(Original Mode)相比，加權模式(Weighted Mode)下時均升力系數整體上有所增加，時均能耗系數則在整體上有明顯的降低。加權模式下，隨著k的增加，升力系數呈增加的趨勢，而時均升力/能耗比先增加后減少，并在k=0.65～0.7附近出現最大值。這個結論和自然界一些鳥類的觀測數據保持一致[29-30]。

圖4 典型相位下的渦量云圖(從上至下的依次為：k=0.4、0.5、0.7)
Fig.4 Vorticity contours at four typical phase angles (from upper to lower：k=0.4、0.5、0.7)

圖5 時均氣動力參數隨k值的變化
Fig.5 Time-averaged aerodynamic coefficients for original mode and weighted mode

3.2 拍動角的影響

本節取k=0.7，令θdown=θup，分別針對-5°～45°等若干拍動角的情況進行了數值模擬，分析拍動角對撲翼過程氣動力特性的影響。

一個拍動周期內升力系數、阻力系數的變化情況如圖6所示，圖中橫坐標為撲翼運動的相位。分別針對下拍過程和上拍過程進行單獨分析。下拍過程中(Phase=0～π)，隨著拍動角θdown的增加，相對于翼型而言來流迎角減小，因此升力系數的峰值減弱，且當θdown>25°時甚至會出現負的升力峰；就阻力系數而言，在一定的θdown范圍內為負值，說明產生了“推力”作用，但是隨著θdown的增加，“推力”的峰值逐漸減小，并最終導致阻力的產生。在上拍過程中(Phase=π～2π)，相對于翼型而言來流的等效迎角為負，因此當上拍角θup較小時，會出現負的升力峰值，但是隨著θup的增加，負升力峰逐漸減弱，并在一定條件下變為正值。

采用和3.1節相同的方法，為升力、阻力系數引入公式(9)所示的加權系數λ，得到拍動過程的時均力學系數如圖7所示，圖中橫坐標為最大拍動角。下拍過程的時均升力系數在θdown=-2°時出現最大值，之后隨著θdown的增加而單調減小，在θdown=30°左右變為負值；時均阻力系數的最小值出現在θdown=0°附近，之后隨著θdown的增加而單調增加，在θdown=21°左右變為正值。上拍過程的時均升力系數隨上拍角的增加而單調增加，并在θup=25°附近變為正值，說明上拍過程仍可能產生升力；就阻力而言，在各個θup條件下其時均值均為正，只是在一定的條件下具有較小的量，從圖7可以確定出這個阻力最小的拍動角約為θup=20°。

這些數值結果表明，翼型的拍動角度對其升阻力特性影響極大，因此采用合適的拍動角是其提高升力、增加推力的關鍵。下拍過程對升力、推力的產生都是有利的，且θdown=0°時具有較好的升阻力特性。由于本文所采用模型的局限性，在各個上拍角條件下，上拍過程中均無法產生推力，但是在合適的角度下能夠保證較小的阻力，且當上拍角大于一定數值時仍能夠產生升力。

3.3 優化的拍動方式

根據3.2節的分析，選擇θdown=0°，針對θup=20°的情況進行了數值模擬。

拍動一周期內升力系數、阻力系數的變化曲線如圖8所示，橫坐標表示拍動的相位。拍動過程中幾個典型狀態下的瞬時流線及壓力云圖如圖9所示(為體現來流的“等效迎角”，速度場減去了翼型在y方向運動速度)。

圖8 優化拍動角情況下一周期內的氣動力系數
Fig.8 Aerodynamic coefficients in a period at the optimal flapping angle

圖9 拍動過程中幾個典型時刻的瞬時流線及壓力云圖
Fig.9 Streamlines and pressure contours at some typical times

下拍過程中(Phase=0～π)，翼型背風區內的流動分離較弱，僅在個別時刻在尾緣附近出現了較小的流動分離，因此整個下拍過程中翼型的升力系數均為正，且存在較大的升力峰值。從壓力云圖也可以看出，下拍過程中在翼型頭部存在較強的負壓區，可誘導出較大的升力以及正推力。

上拍過程的起始及結束階段，翼型的升力系數仍為正值，僅在下拍過程的中間階段出現負升力。這是由于上拍過程中翼型存在較大的上拍角度，相對于來流，翼型的等效“負”迎角減小，因此上拍的不利影響會得到減弱。

下拍-上拍的轉換階段(Phase=0、π)，由于本文所采用的運動規律在此時的俯仰角速度不連續，導致力學系數出現跳躍，但是這些瞬時的強非定?，F象對整個拍動周期的影響相對較弱。

采用與3.1節相同的做法，引入公式(9)所示的加權系數，將力學系數進行加權，得到的整個拍動過程的時均升力系數、阻力系數和能耗系數分別為1.10、-0.14、0.24。定義效率：

則撲翼的推進效率約為58.3%。

4 結 論

本文以海鷗等大型鳥類的巡航飛行狀態為研究背景，選擇S1223翼型以及相應的撲翼運動參數作為研究對象，希望通過這些初步的研究，探討翼型俯仰-沉浮運動過程的升力、推力的產生機理。

本文選用的撲翼方式具有較小的減縮頻率，和昆蟲等的撲翼方式不同，其非定常效應較弱，翼型的靜態因素對氣動力特性影響則較大。翼型通過沉浮運動改變來流的相對迎角，在下拍過程中，來流等效迎角為正，且由于非定常流動的遲滯作用，翼型背風區內的流動分離現象很弱，鑒于S1223翼型優秀的升阻比特性，氣動力合力指向翼型前上方，同時產生“升力”和“推力”(圖10所示)。

文獻[18]的數值模擬結果表明，三維柔性翼型在上拍過程仍有較大的推力產生。而在本文所選用的模型及運動方式條件下，上拍過程沒有能夠產生推力，其原因可能在于：(1)本文選用的翼型在負迎角條件下具有較差的力學特性；(2)本文采用剛性的俯仰-沉浮運動方式，沒有考慮到實際情況下翼型的柔性變形；(3)所選用的撲翼運動參數如頻率、拍動角等存在差異。

下拍時間、拍動角等撲翼參數對翼型的氣動力特性影響較大。在本文的數值模擬條件下，當下拍時間占整個拍動周期的65～70%左右時具有較好的時均升力/能耗比，而下拍角取為0°，上拍角取為20°時具有較好的升力、推力產生性能。

本文僅針對剛性的二維撲翼模型進行了數值模擬，初步分析了下拍時間、拍動角等的影響，對其推力、升力產生機理有了定性的認識。而三維情況下的撲翼運動具有更多的自由度，并且存在柔性變形如翼沿展向的彎曲、扭轉等。這些三維效應對升力、推力的影響很大，下一步我們將重點開展這方面的研究工作。

[1]Brown R H J. The flight of birds[J]. Journal of Experimental Biology, 1952, 25: 90-103.

[2]Greenewalt C H. Dimensional relationships for flying animals[J]. Smithson. Misc. Collect, 1962, 144: 1-46.

[3]Tennekes H. The simple science offlight: from insects to jumbo jets[M]. Boston, MA: MIT Press, 1996.

[4]Pennycuick C J. Speeds and wing-beat frequencies of migration birds compared with calculated benchmarks[J]. J. Exp. Biol., 2001, 204: 3283-3294.

[5]Jones K D, Center K B. Wake structures behind plunging airfoils: a comparison of numerical and experimental results[R]. AIAA paper, 96-0078, 1996.

[6]Angela M B, Andrew A B. Wing and body kinetics of takeoff and landing flight in the pigeon[J]. J. Exp. Biol. 2010, 213: 1651-1658[7]Pennycuick C J. Power requirements for horizontal flight in the pigeon[J]. J. Exp. Biol. 1968, 49: 527-555.

[8]Rayner J M V. A new approach to animal flight mechanics[J]. J. Exp. Biol. 1979, 80: 17-54.

[9]Phlips P J, East R A, Pratt N H. An unsteady lifting line theory of flapping wings with application to the forward flight of birds[J]. J. Fluid Mech. 1981, 112: 97-125.

[10]Smith M J C. Simulating moth wing aerodynamics: Towards the development of flapping-wing technology[J]. AIAA J, 1996, 34: 1348-1355.

[11]Zeng R. Aerodynamic characteristics of flapping-wing MAV simulating bird flight[D]. Nanjing University of Aeronautics and Astronautics, 2004 (in Chinese). 曾銳. 仿鳥微型撲翼飛行器的氣動特性研究[D]. 南京航空航天大學博士學位論文, 2004.

[12]Yu Y L. A theoretical modeling study on unsteady aerodynamics of wing flapping during insect forward flight[D]. University of Science and Technology of China, 2004 (in Chinese). 余永亮. 昆蟲前飛拍翼非定常空氣動力學的理論模化研究[D]. 中國科學技術大學博士學位論文, 2004.

[13]Ito Y, Nakahashi K. Flow simulation of flapping wings of an insect using overset unstructured grid[R]. AIAA paper, 2001-2619.

[14]Miller L A, Peskin C S. When vortices stick: an aerodynamic transition in tiny insect flight[J]. J. Exp. Biol. 2004, 207: 3073-3088.

[15]Liu H, Ellingtona C P. Computational fluid dynamic study of hawkmoth hovering[J]. J. Exp. Biol., 1998, 201: 461-477.

[16]Liang B, Sun M. Aerodynamic interactions between contralateral wings and between wings and body of a model insect at hovering and small speed motions[J]. Chinese Journal of Aeronautics, 2011, 24(4): 396-409.

[17]Ashraf M A, Young J, Lai J C S. Reynolds number, thickness and camber effects on flapping airfoil propulsion[J]. Journal of Fluids and Structures, 2011, 27: 145-160.

[18]Unger R, Haupt M C, Horst P, et al. Fluid-structure analysis of a flexible flapping airfoil at low Reynolds number flow[J]. Journal of Fluids and Structures, 2012, 28: 72-88.

[19]He X, Zhang L P, Zhao Z. The research and development of structured-unstructured hybrid CFD software[C]//The ninth Asian CFD conference, Nanjing, 2012.

[20]Zhang L P, Duan X P, Chang X H, et al. A hybrid dynamic grid generation technique for morphing bodies based on Delaunay graph and local remeshing[J]. ATCA Aerodynamica Sinica, 2009, 27(1): 26-32 (in Chinese). 張來平, 段旭鵬, 常興華, 等. 基于Delaunay背景網格插值方法和局部網格重構的動態混合網格生成技術[J]. 空氣動力學學報, 2009, 27(1): 26-32.

[21]Zhang L P, Chang X H, Duan X P, et al. Applications of dynamic hybrid grid method for three-dimensional moving/deforming boundary problems[J]. Computers & Fluids, 2012, 62: 45-63.

[22]Zhang L P, Wang Z J. A block LU-SGS implicit dual time-stepping algorithm for hybrid dynamic meshes[J]. Computers & Fluids, 2004, 33: 891-916.

[23]Rendall T C S. Allen C B. Reduced surface point selection options for efficient mesh deformation using radial basis functions[J]. Journal of Computational Physics, 2010, 229: 2810-2820.

[24]Lesoinne M, Farhat C. Geometric conservation laws for flow problems with moving boundaries and deformable meshes, and their impact on aero-elastic computations[J]. Computer Methods in Applied Mechanics and Engineering, 1996, 134: 71-90.

[25]Ahn H T, Kallinderis Y. Strongly coupled flow/structure interactions with a geometrically conservative ALE scheme on general hybrid meshes[J]. Journal of Computational Physics, 2006, 219: 671-696.

[26]Liu J, Bai X Z, Zhang H X, et al. Discussion about GCL for deforming grids[J]. Aeronautical Computing Technique, 2009, 39(4): 1-5 (in Chinese). 劉君, 白曉征, 張涵信, 等. 關于變形網格“幾何守恒律”概念的討論[J]. 航空計算技術, 2009, 39(4): 1-5.

[27]Wang Z J, Yang H Q. Unsteady flow simulation using a zonal multi-grid approach with moving boundaries[R]. AIAA Paper, 1994-0057.

[28]Chang X H, Ma R, Zhang L P, et al. Further study on the geometric conservation law for finite volume method on dynamic unstructured mesh[J]. Computers & Fluids, 2015, 120(5): 98-110.

[29]Liu T S. Avian wing geometry and kinematics[J]. AIAA J., 2006, 44: 954-963.

[30]McCroskey W J, McAlister K W. Dynamic stall on advanced airfoil sections[J]. Journal of the American Helicoper Society, 1981, 26(3): 40-50.

Numerical study of the phunging-pitching motion of S1223 airfoil

Chang Xinghua1,2, Ma Rong2, Zhang Laiping1,2,*, He Xin2

(1.StateKeyLaboratoryofAerodynamicsofChinaAerodynamicsResearchandDevelopmentCenter,Mianyang621000,China;2.InstituteofComputationalAerodynamicsofChinaAerodynamicsResearchandDevelopmentCenter,Mianyang621000,China)

Flapping-wing micro air vehicle (FMAV) is one of the most important research directions of future aviation. The bird or insect in nature possess the advanced ability of flying because of thousands of years’ evolution, so the study of their flight kinematics as well as flow mechanism can greatly contribute to the development of FMAV. To investigate the unsteady aerodynamic mechanism of flapping wing, simplified plunging/pitching motion of a rigid 2-D S1223 airfoil is firstly taken into consideration. Numerical method is based on the hybrid dynamic mesh technique and unsteady flow solver. To improve the moving ability of dynamic mesh, the radius basis function is used to calculate the displacement of volume nodes. For unsteady flow solver on moving mesh, geometric conservation law is satisfied by constraining the normal velocity of cell surface. Second-order finite volume method is used. Dual-time stepping method and BLU-SGS implicit scheme are adopted for time marching. The unsteady lifts, thrusts as well as power consumptions with different flapping parameters are abtained. The flow mechanism for lift and thrust is studied. The influence of down-stroke time ratio and stroke angle is analyzed,which would greatly influence the aerodynamic performance and the power consumption. Numerical results show that the ‘static effect’ of the airfoil plays the main role for lift generation.Meanwhile, the unsteady flow induced by the plunging/pitching motion can enhance the lift by weakening the flow separation and increasing the equivalent angle of attack. The maximum lift per power consumption can be got when down-stroke time takes about 65%-70% of the whole stroke cycle, which agrees well with observation data. Further, a set of appropriate flapping angles are specified based on these analyses, which would greatly contribute to the realistic 3-D simulations in the future.

flapping wing; dynamic hybrid mesh; unsteady flow; radius basis function; geometric conservation law

0258-1825(2017)01-0062-09

2015-10-08;

2015-11-05

國家自然科學基金(11532016、11672324);國家重點研發計劃(2016YFB0200701)

常興華(1982-)，男，河南焦作人，副研究員，博士，研究方向：動態混合網格技術，非定常數值模擬技術. E-mail：cxh_cardc@126.com

張來平*(1968-)，男，湖北監利人，研究員，博士，研究方向：非結構網格技術及算法，大規模計算. E-mail: zhanglp_cardc@126.com

常興華, 馬戎, 張來平, 等. S1223翼型俯仰-沉浮運動的非定常氣動特性分析[J]. 空氣動力學學報, 2017, 35(1): 62-70.

10.7638/kqdlxxb-2015.0183 Chang X H, Ma R, Zhang L P, et al. Numerical study of the phunging-pitching motion of S1223 airfoil[J]. Acta Aerodynamica Sinica, 2017, 35(1): 62-70.

V211.3

A doi: 10.7638/kqdlxxb-2015.0183