海鷗翼折轉運動的數值模擬及分析

常興華, 馬 戎, 張來平,*(1. 中國空氣動力研究與發(fā)展中心 空氣動力學國家重點實驗室, 四川 綿陽 61000;. 中國空氣動力研究與發(fā)展中心 計算空氣動力研究所, 四川 綿陽 61000)

0 引 言

撲翼飛行器一直是國內外很多研究機構的一個研究熱點，其設計靈感來源于自然界的昆蟲、鳥類等飛行生物。這些飛行生物經過了億萬年的進化，形成了非常出色的飛行能力。不同于固定翼飛行器，鳥翼即是升力產生機構，也是推力產生機構，其飛行效率更高、噪音更低，并且由于其氣動力產生方式更為靈活，因此機動性能更加出色。和固定翼相比，撲翼飛行的氣動力產生機制以及流動控制機制更加復雜，通過對自然界的鳥類、昆蟲等生物撲翼飛行進行觀察和研究，將有助于推動撲翼飛行器研制工作的進展。

中等體型鳥類的撲翼和昆蟲的撲翼存在較大差異。首先，兩者撲翼的減縮頻率不一致。昆蟲撲翼的減縮頻率很高，雷諾數較低，在翼的拍動過程中相對于翼面的來流迎角很大，因此會產生非常明顯的前緣渦以及后緣渦，其中前緣渦在翼的運動過程中保持不脫落，從而可以維持一個較高的升力[1]。而中等體型鳥類飛行速度較快，撲翼的減縮頻率較低，雷諾數較高，翼型本身的作用更加明顯[2-3]。其次，兩者撲翼的結構和運動方式差異較大。昆蟲翼一般由少量翅脈和大面積的柔性薄膜組成(如蜻蜓翼)，翼面作往復式旋轉，拍動過程中往往伴隨大角度的“8”字型扭轉，翼面相對于翅脈的柔性變形較小，且主要是被動彈性變形，因此早期的研究中將其視為剛性薄板進行簡化處理[4-5]；鳥翼的結構以及變形機制較為復雜，其由肌肉、骨骼、羽毛和多個關節(jié)組成，構造非常精細。其撲動過程具有較多的自由度，如翼的拍動、折轉、收縮、扭轉、飛羽的收縮-合攏等，主動以及被動變形非常劇烈。

由于鳥類撲翼的減縮頻率相對較低，雷諾數相對較高，翼型本身對升力的貢獻作用更為明顯，因此很多經典的研究工作正是基于翼型理論開展的，如Jones[6]、Angela[7]、Pennycuick[8]、Rayner[9]、Phlips[10]等人的工作。這些近似的分析工作定性的給出了鳥翼撲動過程中非定常升力、推力的產生機制。其中“面元法”作為一種簡單有效的數值模擬手段，在鳥類撲翼的研究中使用較多，如Smith[11]采用該方法并結合有限元模型對柔性撲翼運動進行了研究，國內的昂海松研究團隊[12]、余永亮[13]等也分別采用該方法開展了撲翼問題的研究。

雖然這些針對鳥類撲翼運動的研究已經初步揭示了鳥類的升力、推力產生機制，然而這些定性的認識還不能完全滿足撲翼飛行器設計上的需求。鳥類撲翼運動的許多細節(jié)問題如折轉收縮變形、展向扭轉變形、相位等對其力學特性有非常重要的影響，其中必定蘊含著重要的增升、減阻、減能耗機制需要流體力學工作者去不斷研究、挖掘。

計算流體力學的迅速發(fā)展為撲翼運動的精細化研究提供了條件，且已經在昆蟲的撲翼研究中得到了廣泛應用[14-18]。然而由于鳥翼撲動過程的復雜性，三維情況下復雜撲動過程的精細化數值研究還不多見，很多數值研究工作仍基于二維或者簡單的三維撲翼開展[19-20]。本文基于文獻資料構建了海鷗翼幾何模型并設計了簡化的撲動運動模型，采用動態(tài)混合網格技術及非定常數值模擬方法，對撲動過程進行了數值模擬，對拍動角、折轉角的影響進行了對比分析。

1 數值計算方法

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

1.1 動態(tài)混合網格生成技術

在之前的研究工作中，作者所在的研究團隊建立了彈簧松弛法和背景網格映射法相結合的混合網格變形技術，具有較好的變形能力和動網格生成效率，并通過采用局部網格重構技術，提高了針對大變形、大位移、相對運動等復雜動邊界問題的適應能力[22-24]。在最近的研究工作中，進一步耦合了基于徑向基函數[25](RBF)的網格變形技術，其具有優(yōu)越的網格變形能力。標準的RBF方法在處理大規(guī)模網格時效率極差，為了提高其適用性，我們通過選擇有限的參考點來減少RBF算法中矩陣的規(guī)模，以提高計算效率。圖1所示為本文的RBF及相應的參考點選擇算法流程：計算準備階段，根據物面的初始位移，通過貪婪算法確定出參與矩陣求解的參考點；動網格生成的每一步，先對物面點進行插值，并進行誤差檢驗，如果物面誤差不滿足計算需求，則重新建立物面參考點序列，之后再對空間點進行插值。

圖1 RBF網格變形方法及參考點選擇流程Fig.1 Moving grid generation technique based on RBF approach and the selection of reference points

根據文獻中的觀測數據[31]，我們建立了海鷗翼的三維模型(如圖2所示)。其弦向截面為S1223翼型，根弦長c=0.2 m，展長L=0.5 m。本文采用混合網格離散計算域，翼前后緣附近采用非結構的四面體和三棱柱網格單元，其它區(qū)域采用六面體單元，并在尾流以及翼面上下等局部區(qū)域進行了網格加密處理，網格單元總數282萬(由于外形比較簡單，這里不再顯示具體網格分布)。圖2所示的紅點為通過貪婪算法選擇的RBF參考點，物面網格點數4萬，通過篩選只保留了1000個控制點，因此可以極大的提升動態(tài)網格生成效率。圖3所示為物面點的最大以及平均誤差隨參考點數目的變化情況，當選擇1000個控制點時，物面的最大誤差可保持在0.01 mm以下。

圖2 海鷗翼模型以及通過選擇得到的物面參考點Fig.2 Model of seagull wing and the selected reference points for RBF interpolation

圖3 物面最大誤差及平均誤差隨參考點數目變化Fig.3 The maximum and averaged errors of surface geometry with increasing of RBF reference points

1.2 非定常數值模擬方法

基于動態(tài)混合網格的非定常流場解算器采用了格心型的有限體積格式，時間離散采用二階的歐拉后插方法，為提高非定常計算效率，采用了雙時間步和BLU-SGS隱式計算方法。并采用了SA湍流模型模擬湍流流動。算法具體細節(jié)請參見文獻[24]。

在之前的研究中，我們基于誤差分析以及數值測試，對運動網格下的幾何守恒律問題進行了細致的研究，將現有的幾何守恒律算法[26-29]歸納為兩類：限制整體積分誤差的“體限制方法”和限制每個單元邊界面數值誤差的“面限制方法”。在此基礎上提出了一種簡便的滿足幾何守恒律的算法，其只需要約束單元邊界面的法向速度即可。以二階歐拉隱式格式為例，采用如下的法向速度求解方法：

(1)

2 翼型及運動方式

本文采用如下的坐標系定義：坐標原點位于翼根部翼型前緣頂點；x軸指向來流流動方向，y軸指向翼型上方，z軸指向展向方向。

真實海鷗的撲翼運動過程較為復雜，包括繞根部的拍動、沿展向的折轉、扭轉、收縮等。除此之外，還包括一些其它的提高力學特性的運動機制如飛羽的打開/合攏、翼的被動柔性變形等。因此完全考慮真實的撲翼過程是非常困難的。

本文將海鷗翼分為兩段，靠近根部的第一段長度0.2 m，剩下的為第二段，長度0.3 m。只考慮翼型的折疊和拍動，因此將撲翼運動簡化為如下兩種簡單運動的疊加：

1) 翼型的變形：第二段翼的折轉，折轉角定義為θ。

2) 翼型的剛性拍動：繞x軸的拍動(拍動角γ，下拍為正)，旋轉中心設置在翼根部。

圖4給出了翼型撲翼運動的示意圖。折轉角及拍動角的變化規(guī)律如下：

θ=0.5θmax[1-sin(φ+Δφ)]

γ=-γmaxcosφ

(2)

其中θmax、γmax分別為最大折轉角、最大拍動角。φ表示撲翼運動的相位，Δφ表示翼的折轉運動和拍動運動之間的相位差。相位φ的表達式如下：

(3)

式中k表示下拍時間和整個拍動周期的比值，ωdown、ωup分別為下拍、上拍過程的角頻率：

(4)

圖4 撲翼運動示意圖Fig.4 Sketch of the flapping motion

3 數值結果分析

3.1 拍動角的影響

以下分別針對拍動角γmax=5°、10°、20°、30°等四個狀態(tài)進行數值模擬，此時折轉角為0。該狀態(tài)下翼沒有變形，僅做單自由度的拍動。

圖5(a)所示為四個狀態(tài)下一個拍動周期內的升力系數、阻力系數以及能耗系數變化曲線。能耗系數定義為：

(5)

其中S為翼型表面，f為固壁上的氣動力，v表示固壁的運動速度，ρ∞為來流密度，U∞為來流速度，Lref為參考長度，本文取Lref=1.0 m。

一個周期內升力系數、阻力系數均出現了一個波峰和波谷，而能耗系數出現兩個波峰。升力系數峰值、阻力系數最小值均出現在下拍速度最快的時刻(t/T=0.35)；升力系數最小值、阻力系數最大值出現在上拍速度最快時刻(t/T=0.85)；而能耗系數峰值則分別出現在拍動速度最快的兩個時刻。隨著最大拍動角γmax的增加，升力系數、阻力系數以及能耗系數的峰值均單調增大。圖6給出了下拍最快時刻幾個典型狀態(tài)下翼上下表面的壓力云圖，可見第二段翼的壓力分布受γmax影響最為明顯：隨著γmax的增加，翼運動速度加快，第二段翼上表面的前緣出現較大的負壓區(qū)，下表面出現較強的高壓區(qū)，對升力、推力產生積極影響。

時均力學系數隨最大拍動角γmax的變化如圖5(b)所示。隨著γmax的增加，下拍過程(t/T=0～0.7)的時均升力系數單調增大，上拍過程的時均升力系數則單調減小，并在一定的范圍內變?yōu)樨撝?，而整個拍動周期的時均升力系數的變化則相對較小。下拍過程的時均阻力系數隨γmax的增加單調減小，并在一定γmax范圍內為負，說明起到了“推力”作用，而上拍過程則正相反，時均阻力系數隨γmax的增加而增大。由于下拍過程所占整個拍動周期的比重較大，因此總的時均阻力隨著γmax的增加而減小。上拍、下拍以及總的時均能耗系數均隨著γmax的增加而單調增大。

(a) 一周期內的氣動力系數(左：升力系數；中：阻力系數；右：能耗系數)

(b) 時均力學系數(左：升力系數；中：阻力系數；右：能耗系數)

(a) 上表面(Upper surface)

(b) 下表面(Lower surface)

對于上述單自由度剛性拍動的翼型，增加最大拍動角能夠較為明顯的改善下拍過程的升力特性，并能夠在下拍過程產生推力。這實際上是容易理解的，在相同的撲動頻率情況下，增加最大拍動角相當于增大撲動速度，由此導致撲翼各截面“等效迎角”變化加大，進而導致升阻力特性的變化。然而，由于上述計算由于沒有考慮翼型的彎曲、扭轉等變形機制，上拍過程中氣動力的不利影響較大，因此導致撲翼總的力學性能較差。

3.2 第二段翼折轉角的影響

翼面沿展向的折轉變形是鳥類撲翼過程中的一個顯著特征，通過展向的折轉，改變了翼型的有效迎風面積，因此可以改進拍動過程的升阻力特性。本節(jié)在3.1節(jié)的基礎上，考慮第二段翼的折轉，分析最大折疊角度、折轉角相位等對拍動過程升阻力特性的影響。

首先令公式(2)中的Δφ為0，表明下拍最快的時刻折轉角為0，上拍最快的時刻折轉角最大。最大拍動角γmax固定為30°，針對θmax=0°～50°等若干狀態(tài)進行數值模擬。

拍動一周期內的升力系數、阻力系數以及能耗系數變化曲線如圖7(a)所示。由于折轉角按照正弦曲線變化，下拍過程中完全展開，折轉角為0，因此下拍過程的氣動力受最大折轉角的影響較弱；而上拍過程中折轉角最大，對氣動特性的影響相對較大：隨著θmax的增加，上拍過程的負升力峰減弱，阻力、能耗的峰值均減小。圖8所示為上拍最快時刻(t/T=0.85)翼型的壓力云圖及渦量Q等值面，此時翼型折轉角最大，隨著θmax的增加，第二段翼前緣的壓力峰值逐漸減弱，翼梢前緣誘導的分離渦也呈減弱趨勢。從時均力學特性的變化情況可以看出(圖7b)，上拍過程受折轉角影響較大，隨著θmax增加，上拍過程的氣動性能得到較為明顯改善：升力增加，阻力、能耗均降低。而下拍過程所受影響相對較弱。圖9給出了不同折轉角情況下整個拍動周期內的升阻力系數和能耗系數的比值，可見隨著折轉角的增大，升力能耗比增加，阻力能耗比減少，說明增加折轉角對整個拍動周期是有益的。然而本文計算得到的能耗為翼型對流體所做的功，并不等同于翼型撲翼運動的總能耗。實際撲翼過程中存在一個從肌肉化學能到機械能再到有用功的轉換效率問題，隨著折轉角度增加，翼型的運動幅度增加，結構自身加速、減速運動必然會造成額外的能量消耗，可能會導致能量轉換效率的降低。此外，該定量的結論也僅針對本文簡化的撲翼過程。在海鷗的實際飛行過程中，除了翼的折轉運動之外，還存在明顯的展向扭轉運動和沿流向的收縮運動。在這些復合運動的作用下，上拍過程仍有可能產生升力，因此在實際的撲翼過程中折轉角并不是越大越好。因此，在下文關于折轉相位的分析中，根據文獻觀測數據[31]，折轉角度統(tǒng)一設定為θmax=30°。

(a) 一周期內的瞬時氣動力系數(左：升力系數；中：阻力系數；右：能耗系數)

(b) 時均氣動力系數(左：升力系數；中：阻力系數；右：能耗系數)

圖8 上拍最快時刻的物面壓力云圖及空間渦量Q等值面Fig.8 Pressure contours and iso-surface of Q at the time of maximum up-stroke velocity

折轉角相位差Δφ影響拍動過程中翼型展開-折疊的時機，其不僅影響到翼型有效迎風面積的變化，而且會影響到第二段翼的拍動速度。為分析其對翼型動態(tài)氣動力特性的影響，針對Δφ=-60°～90°等若干狀態(tài)進行了對比計算，計算中最大拍動角γmax固定為30°。

(a) 時均升力/時均能耗 (b) 時均阻力/時均能耗

Δφ>0時，翼型展開-折疊的相位提前，在到達下拍中點之前翼型已經完全打開，而在到達上拍中點之前已經完全折疊，Δφ<0時則相反。

圖10、圖11分別給出了拍動過程中的動態(tài)氣動力系數和時均氣動力系數。下拍過程中，隨著Δφ從-60°～90°變化，升力峰、阻力峰、能耗峰值均提前，其原因在于翼完全打開的時間提前。

翼型繞根部的拍動角速度為：

(6)

第二段翼的折轉角速度為：

(7)

上拍過程中，隨著Δφ由負變?yōu)檎?，負的升力系數、阻力系數、能耗系數的峰值均單調增加，原因在于Δφ改變了上拍過程的最快拍動速度。從時均力學系數的變化曲線(圖11b)來看，相對“滯后”的折轉角相位對上拍過程是有利的，可以減小阻力和“負升力”，并減少能耗。根據公式(7)可得，上拍過程中(φ=π～2π)，在φ=π～1.5π-Δφ區(qū)間內，第二段翼的折轉角速度為正，因此可以減弱第二段翼的上拍速度。而隨著Δφ的減小，這一區(qū)間是增大的，因此上拍的“不利”影響會得到減弱。

在本文所用模型和計算參數前提下，Δφ對于下拍和上拍過程的影響恰好相反，因此并沒有對全周期的力學性能(圖11c)產生明顯的有益影響。通過上述分析，折轉運動最理想的情況應該是：下拍過程中，第二段翼向下折疊運動，以增益下拍速度，且在下拍速度最快的時刻保證較小的折疊角，以增加有效迎風面積；上拍過程中，第二段翼向下折疊運動，以減弱上拍速度，且保證上拍速度最快時折疊角較大，以減小有效迎風面積。通過文獻[31]的觀測數據可以看出，翼的折轉角的變化規(guī)律不是簡單的正弦曲線(圖12所示，psi2為折轉角)，其折疊所用的時間遠大于展開所用時間：在整個下拍過程，以及上拍過程的前半段，翼都處在折疊運動的過程中，在上拍過程快結束時再迅速展開，這一觀測結果和本文的分析結果是相吻合的。

圖10 一周期內的氣動力系數(左：升力系數；中：阻力系數；右：能耗系數)Fig.10 Aerodynamic coefficients in a flapping cycle(Left: Lift coefficient; Middle: Drag coefficient; Right: Power coefficient)

(a) 下拍過程(down-stroke)

(b) 上拍過程(Up-stroke)

(c) 全周期(the full cycle)

圖12 海鷗撲翼運動參數(摘自文獻[31])Fig.12 Typical kinematic parameters of a seagull wing (From Ref.[31])

4 結 論

本文通過對簡化的海鷗撲翼過程進行數值模擬，分析了拍動角、折轉角對非定常氣動力以及流場結構的影響。拍動是撲翼的一個主要特征，隨著拍動角度的增加，平均升力、推力均單調增加，同時能耗也明顯增大；折轉角則體現了翼的主動柔性變形，通過折轉，改變上拍-下拍過程中翼的有效迎風面積以及第二段翼面的運動速度，對氣動力影響較大。本文的數值模擬結果表明，折轉角度主要影響上拍過程，尤其對能耗影響較大，而折轉角相位對上拍及下拍過程的平均氣動力均會產生明顯影響。為了增益撲動過程的平均氣動力特性，應當延長翼的折疊時間，減少翼的展開時間，這和實驗觀測的結果一致。

然而本文所得到的結論僅針對本文簡化的撲翼過程，對于實際的撲翼過程，各種運動機制如拍動、扭轉、收縮、折轉是協同運作的，相互之間肯定存在明顯的影響，單純的將各個因素獨立出來分析必然存在較大的近似。希望本文的工作能夠為下一步針對更為真實撲翼過程的研究提供一些指引。

[1]Sun M. Aerodynamics of insect flight[J]. Advances in Mechanics, 2015, 45: 201501: 1-27. (in Chinese)孫茂. 昆蟲飛行的空氣動力學[J]. 力學進展, 2015, 45: 201501: 1-27

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

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

[4]Tong B G, Sun M, Yin X Z. A brief review on domestic research developments in biofluid dynamics of animal flying and swimming[J]. Chinese Journal of Nature, 2005, 27(4): 191-198. (in Chinese)童秉綱, 孫茂, 尹協振. 飛行和游動生物流體力學的國內眼角進展概述[J]. 自然雜志, 2005, 27(4): 191-198

[5]Tong B G, Lu X Y. A review on biomechanics of animal flight and swimming[J]. Advances in Mechanics, 2004, 34(1): 1-8. (in Chinese)童秉綱, 陸夕云. 關于飛行和游動的生物力學研究[J]. 力學進展, 2004, 34(1): 1-8

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

[7]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

[8]Pennycuick C J. Power requirements for horizontal flight in the pigeon[J]. J. Exp. Biol., 1968, 49: 527-555

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

[10]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

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

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

[13]Guan Z W, Yu Y L. Morphing models of a bat wing in flapping flight[C]//The 4th International Conference of Bionic Engineering, 2013, Nanjing

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

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

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

[17]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

[18]Xiao T H, Ang H S, Zhou X C. Numerical method for unsteady flows of flexible flapping wings[J]. Acta Aeronautica et Astronautica Sinica, 2009, 30(6): 990-999. (in Chinese)肖天航, 昂海松, 周新春. 柔性撲翼非定常流場的數值計算方法[J]. 航空學報, 2009, 30(6): 990-999

[19]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

[20]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

[21]He X, Zhang LP, Zhao Z, et al. Validation of the structured-unstructured hybrid CFD software HyperFLOW[J]. Acta Aerodynamica Sinica, 2016, 34(2): 267-275

[22]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]. Acta Aerodynamica Sinica, 2009, 27(1): 26-32. (in Chinese)張來平, 段旭鵬, 常興華, 等. 基于Delaunay背景網格插值方法和局部網格重構的動態(tài)混合網格生成技術[J]. 空氣動力學學報, 2009, 27(1): 26-32

[23]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

[24]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

[25]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

[26]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

[27]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

[28]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

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

[30]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

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

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