S1223翼型俯仰-沉浮運(yùn)動(dòng)的非定常氣動(dòng)特性分析

常興華, 馬 戎, 張來(lái)平,*, 赫 新

(1. 中國(guó)空氣動(dòng)力研究與發(fā)展中心 空氣動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室, 四川 綿陽(yáng) 621000；2. 中國(guó)空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所, 四川 綿陽(yáng) 621000)

S1223翼型俯仰-沉浮運(yùn)動(dòng)的非定常氣動(dòng)特性分析

常興華1,2, 馬 戎2, 張來(lái)平1,2,*, 赫 新2

(1. 中國(guó)空氣動(dòng)力研究與發(fā)展中心 空氣動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室, 四川 綿陽(yáng) 621000；2. 中國(guó)空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所, 四川 綿陽(yáng) 621000)

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

鳥(niǎo)類(lèi)撲翼；動(dòng)態(tài)混合網(wǎng)格；非定常流動(dòng)；RBF動(dòng)網(wǎng)格方法；幾何守恒律

0 引 言

自然界的鳥(niǎo)類(lèi)經(jīng)過(guò)億萬(wàn)年的進(jìn)化，形成了出色的飛行能力，其飛行效率、機(jī)動(dòng)性都遠(yuǎn)遠(yuǎn)超出了目前所有的人造飛行器。具有較低能耗、噪聲，同時(shí)具有優(yōu)越機(jī)動(dòng)性能的小型撲翼飛行器是未來(lái)航空器發(fā)展的一個(gè)重要方向，為開(kāi)展此類(lèi)研制工作，首先需要對(duì)鳥(niǎo)類(lèi)撲翼運(yùn)動(dòng)過(guò)程中的升力、推力產(chǎn)生機(jī)制以及流動(dòng)控制機(jī)理進(jìn)行深入細(xì)致的研究。

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

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

雖然這些針對(duì)鳥(niǎo)類(lèi)撲翼運(yùn)動(dòng)的研究已經(jīng)初步揭示了鳥(niǎo)類(lèi)的升力、推力產(chǎn)生機(jī)制，然而這些定性的認(rèn)識(shí)仍不能滿足撲翼飛行器設(shè)計(jì)上的需求。鳥(niǎo)類(lèi)撲翼運(yùn)動(dòng)的許多細(xì)節(jié)問(wèn)題如翼型、迎角、展向變形/扭轉(zhuǎn)、上/下拍動(dòng)時(shí)間比、拍動(dòng)頻率等對(duì)其力學(xué)特性有非常重要的影響，因此非常有必要針對(duì)這些細(xì)節(jié)問(wèn)題進(jìn)行深入研究。

計(jì)算流體力學(xué)的迅速發(fā)展為撲翼運(yùn)動(dòng)的精細(xì)化研究提供了條件，且已經(jīng)在昆蟲(chóng)的撲翼研究中得到了廣泛應(yīng)用[13-16]。然而由于鳥(niǎo)類(lèi)撲翼過(guò)程的復(fù)雜性，三維情況下復(fù)雜撲翼過(guò)程的精細(xì)化數(shù)值研究還不多見(jiàn)，很多數(shù)值研究工作仍基于二維或者簡(jiǎn)單的三維撲翼開(kāi)展[17-18]。

本文以海鷗的撲翼飛行為研究背景，選擇二維的S1223翼型作為研究對(duì)象，采用動(dòng)態(tài)混合網(wǎng)格技術(shù)及非定常數(shù)值模擬方法，對(duì)其剛性俯仰-沉浮運(yùn)動(dòng)進(jìn)行了研究，得到了下拍時(shí)間比、拍動(dòng)迎角等參數(shù)的影響規(guī)律。這些工作為后續(xù)的復(fù)雜三維問(wèn)題的數(shù)值研究奠定了基礎(chǔ)。

1 數(shù)值計(jì)算方法

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

1.1 動(dòng)態(tài)混合網(wǎng)格生成技術(shù)

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

在最近的研究工作中，我們進(jìn)一步發(fā)展了基于徑向基函數(shù)(RBF)的網(wǎng)格變形技術(shù)，該技術(shù)具有優(yōu)越的網(wǎng)格變形能力，是目前國(guó)內(nèi)外CFD工作者爭(zhēng)相研究的一個(gè)熱點(diǎn)。而標(biāo)準(zhǔn)的RBF方法在處理大規(guī)模網(wǎng)格時(shí)效率極差，為了提高RBF網(wǎng)格變形技術(shù)的適用性，我們參考文獻(xiàn)[23]的做法，通過(guò)選擇有限的參考點(diǎn)來(lái)減少RBF算法中矩陣的規(guī)模，以提高計(jì)算效率。

1.2 非定常數(shù)值模擬方法

基于動(dòng)態(tài)混合網(wǎng)格的NS方程解算器采用了格心型的有限體積格式，時(shí)間離散采用二階的歐拉后插方法，為提高非定常計(jì)算效率，采用了雙時(shí)間步算法和BLU-SGS隱式計(jì)算方法。算法具體細(xì)節(jié)請(qǐng)參見(jiàn)文獻(xiàn)[22]。

動(dòng)網(wǎng)格下的NS方程離散算法必須滿足幾何守恒律。運(yùn)動(dòng)網(wǎng)格下的幾何守恒律(GCL)特指體積守恒，對(duì)于控制方程，假設(shè)流場(chǎng)均勻，則可得到GCL方程如下：

其中V為控制體Ω的體積，?Ω為控制體邊界，w表示控制體外邊界的運(yùn)動(dòng)速度，n為控制體邊界的外法向。導(dǎo)數(shù)符號(hào)“d”表示跟隨控制體的隨體導(dǎo)數(shù)，而不是跟隨流體的物質(zhì)導(dǎo)數(shù)。

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

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

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

在隨后的研究中，我們進(jìn)一步從誤差分析的角度對(duì)各種運(yùn)動(dòng)網(wǎng)格幾何守恒算法進(jìn)行了理論分析，發(fā)現(xiàn)可以將現(xiàn)有的幾何守恒算法歸為兩類(lèi)，即限制整體積分誤差的“體限制方法”和限制每個(gè)邊數(shù)值誤差的“面限制方法”，并分別可以針對(duì)不同的時(shí)間格式寫(xiě)成統(tǒng)一的形式。進(jìn)一步的對(duì)比發(fā)現(xiàn)，各種幾何守恒算法的計(jì)算結(jié)果大體一致，但是我們發(fā)展的幾何守恒算法在三維情況下更為簡(jiǎn)便。具體的理論分析詳見(jiàn)文獻(xiàn)[28]。

2 翼型及俯仰-沉浮運(yùn)動(dòng)

根據(jù)文獻(xiàn)描述[29]，本文選用二維S1223翼型作為研究對(duì)象，翼型弦長(zhǎng)c=0.2 m。將整個(gè)拍動(dòng)周期分為下拍和上拍兩個(gè)階段，其沉浮運(yùn)動(dòng)規(guī)律如下：

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

翼型繞其1/4弦長(zhǎng)處進(jìn)行俯仰，俯仰角(低頭為正)變化規(guī)律如下：

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

本文只針對(duì)二維問(wèn)題進(jìn)行數(shù)值模擬，其相當(dāng)于三維問(wèn)題的一個(gè)展向截面(如圖1所示)。

本文需要研究俯仰-沉浮運(yùn)動(dòng)的能耗特性，能耗系數(shù)CP的定義如下：

其中f表示流體作用在物體邊界上的力，w表示邊界元的運(yùn)動(dòng)速度，?Ω表示物體表面。

3 數(shù)值結(jié)果分析

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

3.1 下拍時(shí)間比對(duì)升力及能耗的影響

觀察發(fā)現(xiàn)[30]，海鷗等鳥(niǎo)類(lèi)巡航飛行時(shí)其下拍時(shí)間占到整個(gè)拍動(dòng)周期的60%～80%，說(shuō)明下拍時(shí)間比k與其力學(xué)特性關(guān)系密切。本節(jié)不考慮翼型的俯仰，僅考慮其上下沉浮運(yùn)動(dòng)，針對(duì)k=0.40～0.75等幾種情況進(jìn)行數(shù)值模擬。

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

圖4所示為拍動(dòng)一個(gè)周期內(nèi)幾個(gè)典型相位的渦量云圖，在Phase=0.5π、1.5π兩個(gè)相位下，不同的k值對(duì)應(yīng)的流場(chǎng)結(jié)構(gòu)基本類(lèi)似，只是由于k值的變化影響了其下拍和上拍過(guò)程的等效迎角，因此背風(fēng)區(qū)內(nèi)的旋渦大小、形態(tài)等存在一些差異；而Phase=0、π兩個(gè)相位下，翼型的法向運(yùn)動(dòng)速度為0，但是加速度處在整個(gè)拍動(dòng)周期的峰值，此時(shí)不同k值條件下流場(chǎng)的旋渦結(jié)構(gòu)存在較大差異。以Phase=0為例，此時(shí)處于上拍結(jié)束的急劇減速以及下拍起始的急劇加速階段，翼型帶動(dòng)周?chē)黧w急劇加速/減速，受“附加質(zhì)量效應(yīng)”影響較大，且由于本文設(shè)計(jì)的運(yùn)動(dòng)規(guī)律在此時(shí)的加速度不連續(xù)，因此導(dǎo)致升力系數(shù)存在“不規(guī)律”的波動(dòng)(圖3c)。然而，由于翼型拍動(dòng)的減縮頻率較小，來(lái)流速度相對(duì)于拍動(dòng)速度是一個(gè)較大的量，因此這些“不規(guī)則”的旋渦結(jié)構(gòu)會(huì)很快消失在尾跡，升力系數(shù)也會(huì)很快進(jìn)入單調(diào)增加的趨勢(shì)。

圖3 不同k值下一個(gè)拍動(dòng)周期內(nèi)的氣動(dòng)力系數(shù)
Fig.3 Aerodynamic coefficients in a period at different value of k

上述的數(shù)值模擬結(jié)果表明，翼型的下拍過(guò)程對(duì)升力的貢獻(xiàn)較大，而上拍過(guò)程則產(chǎn)生負(fù)升力。實(shí)際上，鳥(niǎo)類(lèi)在撲翼過(guò)程中正是通過(guò)翅膀展向的折轉(zhuǎn)來(lái)突出下拍的“有益”作用(下拍時(shí)翼面全部展開(kāi)、增大有效升力面)，減弱上拍的“不利”影響(上拍時(shí)適當(dāng)折疊撲翼，減小上拍迎風(fēng)阻力面)。為了近似反應(yīng)這種“三維效應(yīng)”，本文對(duì)升力、阻力、能耗等引入一個(gè)近似的加權(quán)系數(shù)：

λ在0.5和1.0之間周期性變化，φ表示加權(quán)系數(shù)和撲翼運(yùn)動(dòng)之間的相位差，本文令φ=0。將計(jì)算得到的瞬時(shí)氣動(dòng)力系數(shù)乘上加權(quán)系數(shù)λ，則本文的二維問(wèn)題可近似看作一個(gè)三維的翼型在做俯仰-沉浮運(yùn)動(dòng)的同時(shí)改變其展向長(zhǎng)度。

考慮加權(quán)前后，撲翼運(yùn)動(dòng)的時(shí)均升力系數(shù)、時(shí)均能耗系數(shù)以及它們之間的比值隨k值的變化如圖5所示：與原始模式(Original Mode)相比，加權(quán)模式(Weighted Mode)下時(shí)均升力系數(shù)整體上有所增加，時(shí)均能耗系數(shù)則在整體上有明顯的降低。加權(quán)模式下，隨著k的增加，升力系數(shù)呈增加的趨勢(shì)，而時(shí)均升力/能耗比先增加后減少，并在k=0.65～0.7附近出現(xiàn)最大值。這個(gè)結(jié)論和自然界一些鳥(niǎo)類(lèi)的觀測(cè)數(shù)據(jù)保持一致[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 時(shí)均氣動(dòng)力參數(shù)隨k值的變化
Fig.5 Time-averaged aerodynamic coefficients for original mode and weighted mode

3.2 拍動(dòng)角的影響

本節(jié)取k=0.7，令θdown=θup，分別針對(duì)-5°～45°等若干拍動(dòng)角的情況進(jìn)行了數(shù)值模擬，分析拍動(dòng)角對(duì)撲翼過(guò)程氣動(dòng)力特性的影響。

一個(gè)拍動(dòng)周期內(nèi)升力系數(shù)、阻力系數(shù)的變化情況如圖6所示，圖中橫坐標(biāo)為撲翼運(yùn)動(dòng)的相位。分別針對(duì)下拍過(guò)程和上拍過(guò)程進(jìn)行單獨(dú)分析。下拍過(guò)程中(Phase=0～π)，隨著拍動(dòng)角θdown的增加，相對(duì)于翼型而言來(lái)流迎角減小，因此升力系數(shù)的峰值減弱，且當(dāng)θdown>25°時(shí)甚至?xí)霈F(xiàn)負(fù)的升力峰；就阻力系數(shù)而言，在一定的θdown范圍內(nèi)為負(fù)值，說(shuō)明產(chǎn)生了“推力”作用，但是隨著θdown的增加，“推力”的峰值逐漸減小，并最終導(dǎo)致阻力的產(chǎn)生。在上拍過(guò)程中(Phase=π～2π)，相對(duì)于翼型而言來(lái)流的等效迎角為負(fù)，因此當(dāng)上拍角θup較小時(shí)，會(huì)出現(xiàn)負(fù)的升力峰值，但是隨著θup的增加，負(fù)升力峰逐漸減弱，并在一定條件下變?yōu)檎怠?/p>

采用和3.1節(jié)相同的方法，為升力、阻力系數(shù)引入公式(9)所示的加權(quán)系數(shù)λ，得到拍動(dòng)過(guò)程的時(shí)均力學(xué)系數(shù)如圖7所示，圖中橫坐標(biāo)為最大拍動(dòng)角。下拍過(guò)程的時(shí)均升力系數(shù)在θdown=-2°時(shí)出現(xiàn)最大值，之后隨著θdown的增加而單調(diào)減小，在θdown=30°左右變?yōu)樨?fù)值；時(shí)均阻力系數(shù)的最小值出現(xiàn)在θdown=0°附近，之后隨著θdown的增加而單調(diào)增加，在θdown=21°左右變?yōu)檎怠Ｉ吓倪^(guò)程的時(shí)均升力系數(shù)隨上拍角的增加而單調(diào)增加，并在θup=25°附近變?yōu)檎?，說(shuō)明上拍過(guò)程仍可能產(chǎn)生升力；就阻力而言，在各個(gè)θup條件下其時(shí)均值均為正，只是在一定的條件下具有較小的量，從圖7可以確定出這個(gè)阻力最小的拍動(dòng)角約為θup=20°。

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

3.3 優(yōu)化的拍動(dòng)方式

根據(jù)3.2節(jié)的分析，選擇θdown=0°，針對(duì)θup=20°的情況進(jìn)行了數(shù)值模擬。

拍動(dòng)一周期內(nèi)升力系數(shù)、阻力系數(shù)的變化曲線如圖8所示，橫坐標(biāo)表示拍動(dòng)的相位。拍動(dòng)過(guò)程中幾個(gè)典型狀態(tài)下的瞬時(shí)流線及壓力云圖如圖9所示(為體現(xiàn)來(lái)流的“等效迎角”，速度場(chǎng)減去了翼型在y方向運(yùn)動(dòng)速度)。

圖8 優(yōu)化拍動(dòng)角情況下一周期內(nèi)的氣動(dòng)力系數(shù)
Fig.8 Aerodynamic coefficients in a period at the optimal flapping angle

圖9 拍動(dòng)過(guò)程中幾個(gè)典型時(shí)刻的瞬時(shí)流線及壓力云圖
Fig.9 Streamlines and pressure contours at some typical times

下拍過(guò)程中(Phase=0～π)，翼型背風(fēng)區(qū)內(nèi)的流動(dòng)分離較弱，僅在個(gè)別時(shí)刻在尾緣附近出現(xiàn)了較小的流動(dòng)分離，因此整個(gè)下拍過(guò)程中翼型的升力系數(shù)均為正，且存在較大的升力峰值。從壓力云圖也可以看出，下拍過(guò)程中在翼型頭部存在較強(qiáng)的負(fù)壓區(qū)，可誘導(dǎo)出較大的升力以及正推力。

上拍過(guò)程的起始及結(jié)束階段，翼型的升力系數(shù)仍為正值，僅在下拍過(guò)程的中間階段出現(xiàn)負(fù)升力。這是由于上拍過(guò)程中翼型存在較大的上拍角度，相對(duì)于來(lái)流，翼型的等效“負(fù)”迎角減小，因此上拍的不利影響會(huì)得到減弱。

下拍-上拍的轉(zhuǎn)換階段(Phase=0、π)，由于本文所采用的運(yùn)動(dòng)規(guī)律在此時(shí)的俯仰角速度不連續(xù)，導(dǎo)致力學(xué)系數(shù)出現(xiàn)跳躍，但是這些瞬時(shí)的強(qiáng)非定?，F(xiàn)象對(duì)整個(gè)拍動(dòng)周期的影響相對(duì)較弱。

采用與3.1節(jié)相同的做法，引入公式(9)所示的加權(quán)系數(shù)，將力學(xué)系數(shù)進(jìn)行加權(quán)，得到的整個(gè)拍動(dòng)過(guò)程的時(shí)均升力系數(shù)、阻力系數(shù)和能耗系數(shù)分別為1.10、-0.14、0.24。定義效率：

則撲翼的推進(jìn)效率約為58.3%。

4 結(jié) 論

本文以海鷗等大型鳥(niǎo)類(lèi)的巡航飛行狀態(tài)為研究背景，選擇S1223翼型以及相應(yīng)的撲翼運(yùn)動(dòng)參數(shù)作為研究對(duì)象，希望通過(guò)這些初步的研究，探討翼型俯仰-沉浮運(yùn)動(dòng)過(guò)程的升力、推力的產(chǎn)生機(jī)理。

本文選用的撲翼方式具有較小的減縮頻率，和昆蟲(chóng)等的撲翼方式不同，其非定常效應(yīng)較弱，翼型的靜態(tài)因素對(duì)氣動(dòng)力特性影響則較大。翼型通過(guò)沉浮運(yùn)動(dòng)改變來(lái)流的相對(duì)迎角，在下拍過(guò)程中，來(lái)流等效迎角為正，且由于非定常流動(dòng)的遲滯作用，翼型背風(fēng)區(qū)內(nèi)的流動(dòng)分離現(xiàn)象很弱，鑒于S1223翼型優(yōu)秀的升阻比特性，氣動(dòng)力合力指向翼型前上方，同時(shí)產(chǎn)生“升力”和“推力”(圖10所示)。

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

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

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

[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). 曾銳. 仿鳥(niǎo)微型撲翼飛行器的氣動(dòng)特性研究[D]. 南京航空航天大學(xué)博士學(xué)位論文, 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). 余永亮. 昆蟲(chóng)前飛拍翼非定?？諝鈩?dòng)力學(xué)的理論?；芯縖D]. 中國(guó)科學(xué)技術(shù)大學(xué)博士學(xué)位論文, 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). 張來(lái)平, 段旭鵬, 常興華, 等. 基于Delaunay背景網(wǎng)格插值方法和局部網(wǎng)格重構(gòu)的動(dòng)態(tài)混合網(wǎng)格生成技術(shù)[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 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). 劉君, 白曉征, 張涵信, 等. 關(guān)于變形網(wǎng)格“幾何守恒律”概念的討論[J]. 航空計(jì)算技術(shù), 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

國(guó)家自然科學(xué)基金(11532016、11672324);國(guó)家重點(diǎn)研發(fā)計(jì)劃(2016YFB0200701)

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

張來(lái)平*(1968-)，男，湖北監(jiān)利人，研究員，博士，研究方向：非結(jié)構(gòu)網(wǎng)格技術(shù)及算法，大規(guī)模計(jì)算. E-mail: zhanglp_cardc@126.com

常興華, 馬戎, 張來(lái)平, 等. S1223翼型俯仰-沉浮運(yùn)動(dòng)的非定常氣動(dòng)特性分析[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 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