俯仰振蕩翼型流場的離散渦數(shù)值模擬

傅翔， 李高華， 王福新

(上海交通大學(xué) 航空航天學(xué)院， 上海 200240)

0 引言

拍動翼運動是自然界生物最常見的機動形式之一，魚類，鳥類以及昆蟲都是利用拍動翼來獲得升力或者推力。一個完整的拍動翼過程可以分解成沉浮振蕩和俯仰振蕩兩種最基本的運動?？茖W(xué)研究表明拍動翼的氣動效率遠(yuǎn)遠(yuǎn)高于噴氣飛機或者螺旋槳的氣動效率[1]，所以吸引了越來越多空氣動力學(xué)家的關(guān)注。

國內(nèi)外已經(jīng)獲得了一些針對拍動翼的研究結(jié)果。1973年，Weis-Fogh[2]進行了大量的昆蟲實驗，提出了一種理論，認(rèn)為升力的產(chǎn)生主要在于翅膀合攏后的快速打開過程。該發(fā)現(xiàn)在研究昆蟲飛行的非定常流動中起到了先導(dǎo)性的作用，并且被認(rèn)為是歷史上發(fā)現(xiàn)的第一個非定常效應(yīng)。Weis-Fogh的實驗是通過將昆蟲拴在風(fēng)洞中來進行觀察和測量的，這與自然界中的實際情況顯然是不符合的，不能夠真實地反映昆蟲的翅膀在拍動過程中的真實受力情況，但是，正是在Weis-Fogh等人的理論提出之后，人們開始從非定常流動角度去研究拍動翼的復(fù)雜流動機理。Triantafyllou[3,4]等人對NACA0012翼型進行了水動力試驗，翼型做沉浮俯仰組合運動，他們選取了一定范圍的參數(shù)，對推力系數(shù)和效率進行了測量。Anderson[5]利用PIV對二維拍動翼進行了研究，總結(jié)出了尾渦形態(tài)與振蕩頻率之間的關(guān)系。

而隨著數(shù)值模擬技術(shù)的興起和發(fā)展，越來越多的研究者選擇使用數(shù)值模擬的方法來研究振蕩翼型的運動機理。Guglielmini[6]等人利用數(shù)值模擬的方法研究了俯仰沉浮組合運動的振蕩翼型周圍的流場，認(rèn)為推力隨時間的變化主要受制于翼型的運動細(xì)節(jié)。Wang[7]等采用基于同位網(wǎng)格的SIMPLE算法對進行沉浮俯仰組合運動的二維振蕩翼型周圍的粘性流場進行了數(shù)值模擬，分析了翼型所受的力和尾渦結(jié)構(gòu)之間的關(guān)系，并提出一種用來確定渦的脫落頻率的方法。然而傳統(tǒng)的數(shù)值計算方法在研究拍動翼問題時還存在問題，一方面是因為振蕩翼型流場過于復(fù)雜，渦的尺度大小不一，導(dǎo)致計算精度難以保證；二是傳統(tǒng)數(shù)值計算方法在計算非定常問題時所需的計算量太大，常常需要幾天甚至幾個月的時間才能算完一個拍動周期。

與傳統(tǒng)的數(shù)值計算方法不同，離散渦方法作為一種拉格朗日求解方法，不需要設(shè)置網(wǎng)格，計算量小，在計算拍動翼問題時有著得天獨厚的優(yōu)勢。事實上，已經(jīng)有學(xué)者利用離散渦方法對拍動翼問題進行了研究，但由于原有的離散渦方法無法準(zhǔn)確模擬前緣渦的脫落，導(dǎo)致計算結(jié)果精度較差[8]。本文將利用邊界層渦量和分離點處脫落渦量的關(guān)系，推導(dǎo)出了高攻角情況下前緣渦脫落位置以及環(huán)量的計算公式，再對俯仰振蕩翼型流場進行離散渦數(shù)值模擬，并結(jié)合流場變化對翼型所受升力的演化過程進行分析。

1 離散渦方法的基本理論

二維不可壓流動由連續(xù)性方程和Navier-Stokes方程控制[9]，即式(1)、(2)。

(1)

(2)

對式(2)取旋度，得式(3)。

(3)

(4)

如果略去粘性擴散效應(yīng)，則可得式(5)。

(5)

式(5)即等價于無粘歐拉方程。

離散渦方法，首先是離散，把平面上的連續(xù)渦量離散化為一系列的點渦。點渦的數(shù)量分布和位置坐標(biāo)事先確定，一種分布方法是根據(jù)需要的離散渦數(shù)將翼型按照弦長分段，在每一段的四分之一長處設(shè)置點渦。

離散后，利用點渦的速度環(huán)量值算出離散渦的誘導(dǎo)速度分布。每個點渦的速度由兩部分組成,一部分是流場本身的速度，另一部分是由其它離散渦共同誘導(dǎo)出的速度，該速度可以通過畢奧薩伐爾定律算出式(6)。

(6)

最后根據(jù)點渦前一時刻的速度和位置，可以計算下一時刻的新位置，采用歐拉近似，有式(7)。

(7)

由此可以模擬出點渦每一時刻的位置，進而模擬出整個渦量場，從而獲得流場每一時刻的特征量。

假設(shè)離散的渦元總數(shù)為N，則在時刻t場點r處渦量的近似值可以表示為式(8)。

(8)

其中δ(r)是狄拉克函數(shù)，χj(t)表示t時刻第j個渦元的位置，于是第i個渦元的運動可以描述如式(9)、(10)。

(9)

(10)

式(9)和式(10)即為經(jīng)典的點渦運動微分方程。

2 離散渦方法的數(shù)值化處理

2.1 非定常流場計算

假設(shè)一個點渦被放置在(x0,y0)處，那么根據(jù)畢奧薩伐爾定律，在xy平面上任意(x,y)處的速度勢為式(11)。

(11)

將其分解到x方向和y方向上有式(12)、(13)。

(12)

(13)

用矩陣的形式來表達的話，點渦Γj(xj,yj)在控制點i(xi,yi)處的誘導(dǎo)速度為式(14)。

(14)

式(14)中的r表示這兩個點之間的距離，為式(15)。

(15)

將Γj提取出來(亦即代入Γj=1)，可以得到單位點渦的誘導(dǎo)速度公式，將其與控制點i處的法向向量相點乘，可以得到位于j處的單位點渦在控制點i處所產(chǎn)生的誘導(dǎo)速度在當(dāng)?shù)胤ㄏ虻姆至浚Q其為影響系數(shù)aij，有式(16)。

aij=(ui,wi)·ni

(16)

注意到不可穿透邊界條件的實現(xiàn)需要滿足在翼型表面法向速度為零，對于定常情況，這里的法向速度由兩部分構(gòu)成，一部分是自由來流在當(dāng)?shù)胤ㄏ蚍至?，另一部分是所有點渦的誘導(dǎo)速度的法向分量，可以用下面的公式來表達為式(17)、(18)。

(17)

令Vi=-U∞·ni,有

(18)

式(18)即為點渦的邊界條件公式，對于控制點1，會受到所有的n個點渦的影響，則可以得到關(guān)于控制點1的邊界條件公式為式(19)。

a11Γ1+a12Γ2+a13Γ3+…+a1nΓn=V1

(19)

對所有的n個控制點都有這樣的方程，將其列在一起，可以得到一個矩陣形式的方程組為式(20)。

(20)

注意到式子中a矩陣由點渦和控制點的位置確定，是已知的，右邊的V矩陣由自由來流和各控制點法向量決定，也是已知的，未知的是布置在翼型表面的n個點渦的強度，通過求解這個矩陣方程可以得到初始的定常情況下的翼型表面的點渦強度。

在非定常情況下，會有新生渦的脫落，這時因為要考慮到新生渦與翼型表面渦之間的相互影響，以及新生渦對控制點處產(chǎn)生的誘導(dǎo)速度，所以需要對矩陣方程a進行修正。

當(dāng)翼型的攻角很小的時候，通常認(rèn)為沒有前緣渦的脫落，僅有尾緣渦脫落，尾緣渦在脫落以后強度不會再發(fā)生改變，所以對于每個脫渦的瞬時來說，有N+1個未知數(shù)，那么就需要N+1個方程才能求出結(jié)果，這里第N+1個方程使用環(huán)量守恒來構(gòu)建，即式(21)。

(21)

同時對矩陣a做出相應(yīng)的修改，將尾緣渦的誘導(dǎo)速度考慮進去，比較方便的處理方法是在方程右邊的V矩陣中添加尾緣渦項為式(22)。

V(t)i=-[U(t)∞+U(t-Δt)w]·ni

(22)

那么對于每一個控制點處的方程也相應(yīng)的進行修改為式(23)。

ai1Γ1+ai2Γ2+…+ainΓn+aiwΓW=V(t)i

(23)

注意這里的ΓW是新生的尾緣渦的強度，而在環(huán)量守恒中的尾緣渦強度Γ尾緣渦則是所有尾緣渦的強度之和?？梢缘玫叫碌木仃嚪匠虨槭?24)。

(24)

其中Γ(t-Δt)是上一時刻流場中的總渦量。

當(dāng)翼型的攻角較大或者攻角變化較大時，還要考慮前緣渦的脫落，關(guān)于前緣渦的脫落還沒有統(tǒng)一的處理方式，一種方案是由實驗或者經(jīng)驗確定分離點的位置，然后由其他的條件決定新生渦的強度，另一種方案是先通過邊界層理論決定初生渦的強度，然后由其他的條件來決定其位置。本文所采用的方法是，將前緣渦布置在翼型前緣附近的弦線上，這樣就確定了前緣渦脫落的位置，接下來就是確定前緣渦的強度。如圖1所示。

圖1 前緣渦脫落示意圖

可以將虛線方框內(nèi)包含的部分當(dāng)作是當(dāng)前時刻從分離點脫落的前緣渦，方框內(nèi)的渦量的增量可以通過公式得出式(25)。

(25)

在該分離點處，有ql≈0，而qu則可以從上一時間點求出，此時有式(26)。

(26)

另外由于前緣上的單位長度的渦量持續(xù)地脫落，可以得到式(27)。

(27)

其中ΓLE為前緣上的表面渦的強度，Δl為該渦所在微元段的長度，由此求得前緣渦的強度為式(28)、(29)。

ΓLEW=ΨΓ1

(28)

(29)

將式(28)并入矩陣a，可以得到新的矩陣方程為式(30)。

(30)

至此所有的求解點渦強度的方程已經(jīng)全部給出。

2.2 氣動力計算

以往離散渦方法中都是用非定常伯努利方程來求解氣動力，求解過程復(fù)雜且難以保證計算精度。這里采用吳鎮(zhèn)遠(yuǎn)提出的渦量矩定理[10]來計算物體所受的氣動力，即為式(30)。

?RsvdR

(31)

其中F為流體施加在浸沒在其中并相對于流體運動的固體上的作用力。Rs為固體所占據(jù)的區(qū)域，R∞為流體和固體所共同占據(jù)的無限大的空間，v是速度場，ω是渦量場，ρ是流體的密度，r是位置矢量。

將F正交分解到x,y方向上即可得出這兩個方向上的分力。

2.3 計算流程

在本離散渦算法的計算中，均采用了突然啟動的模式，即流動從靜止?fàn)顟B(tài)突然啟動而后保持設(shè)定好的速度勻速前進。

在求解的過程中，首先設(shè)定好并初始化所有參數(shù)，然后根據(jù)需求在平板或者NACA0012翼型上布置好表面點渦(平板100個，NACA0012翼型300個)。再依據(jù)前文給出的誘導(dǎo)速度計算公式，計算出定常情況下的a矩陣，然后進入時間推進，每一個時間點根據(jù)渦脫落的設(shè)置確定每一個新生的前緣渦和尾緣渦的初始位置，接下來將定常情況下的n階a矩陣轉(zhuǎn)化為包含前緣渦尾緣渦脫落的非定常情況的n+2階矩陣。下一步就是根據(jù)翼型的運動情況和流體速度計算rhs矩陣，然后解矩陣方程即可獲得每一時刻的離散點渦的渦量大小。得到了渦量大小后就可以根據(jù)式(12)和式(13)計算出每一個點渦的運動速度，從而確定其下一時刻的位置，這里需要注意的是，翼型表面點渦滿足邊界條件所以始終保持在翼型上相同的位置，每一步位置會發(fā)生變化的是前緣渦和尾緣渦，但是表面點渦的渦量大小會收到脫落的前緣渦和尾緣渦的影響每一時刻都在變化。在每個瞬時獲得了每個點渦的渦量大小，渦量變化率和位置，速度這些參數(shù)之后，就可以利用渦量矩定理進行力的計算了。重復(fù)時間推進后的步驟直到時間達到設(shè)定值，程序運行結(jié)束。整個離散渦程序計算流程，如圖2所示。

圖2 離散渦程序計算流程

3 計算結(jié)果分析

3.1 算例驗證

首先利用離散渦方法計算出了平板翼型從0°攻角勻速抬起至45°攻角的流場。平板翼型弦長為1 m，來流速度為1 m/s，時間步長取為0.05 s。

平板翼型附近的點渦分布及演化過程，圖3所示。

可以看到在攻角較小的時候，流場類似定常算例中的流場，有明顯的啟動渦的卷起。當(dāng)攻角達到設(shè)定值時，前緣渦開始脫落，可以看到前緣渦的脫落對尾緣渦造成了影響，尾緣渦的流線出現(xiàn)了凸起。前緣渦在前緣不斷聚集，呈現(xiàn)順時針旋轉(zhuǎn)，而受到前緣渦影響的尾緣渦也在尾緣聚集，呈逆時針旋轉(zhuǎn)。當(dāng)前緣渦發(fā)展到一定大小時，逐漸被尾緣渦夾斷脫落，尾緣渦亦如是，如此往復(fù)，形成了如圖所示的卡門渦街。

同時程序還用渦量矩定理計算了這一過程中產(chǎn)生的升力，與Darrel K. Robertson[11]的實驗結(jié)果對比，誤差在可以接受的范圍內(nèi)，如圖4所示。

圖3 翼型附近點渦分布及演化過程

圖4 平板翼型抬起過程中的升力

3.2 俯仰振蕩翼型流場分析

在驗證了程序的可靠性后，對俯仰振蕩翼型的非定常流場進行模擬。這里所采用的翼型為弦長1 m的NACA0012翼型，來流速度均取為1 m/s，翼型上分布了300個離散點渦。

翼型采用的振蕩運動方式為正弦振蕩，這是自然界及工程領(lǐng)域最普遍存在的一種周期性振蕩方式，其運動方程可以用下面的式子來表示為式(31)。

αt=α0+αmaxsin(ωt)

(31)

式子中αt表示t時刻翼型的攻角，α0為翼型初始攻角，αmax為翼型攻角的振幅。ω為俯仰頻率，這里給出縮減頻率k，有式(32)。

(32)

本文主要關(guān)心的是俯仰運動的振幅對于升力的影響，因此對所有算例設(shè)定相同的參數(shù)α0=0，k=0.1，dt=0.02 s。翼型攻角的振幅則取了兩個不同的值：αmax=5°和αmax=30°。具體的數(shù)值模擬結(jié)果和分析如下。

(1) αmax=5°

首先從較小的攻角振幅開始模擬，從程序輸出的結(jié)果，如圖5所示。

圖5 αmax=5°時的點渦分布及演化

從圖5可以看到攻角振幅為5°時，流場與翼型靜止不動時的流場很類似。流動剛開始時，尾緣出現(xiàn)了很明顯的啟動渦，但是前緣渦沒有出現(xiàn)卷起，而是順著翼型表面向下游平順的運動。當(dāng)前緣渦運動到尾緣附近時，可以看到尾緣渦受到了前緣渦強烈的干擾，尾緣渦的脫落不再平順。當(dāng)前緣渦進入尾緣渦的范圍內(nèi)時，與尾緣渦交錯形成一個小的渦對，由于時間步長的設(shè)置，這里的圖像中每個小渦對里面的離散渦點數(shù)不是很多，看起來不夠明顯。同時可以看到在初始的時候，翼型正攻角，前緣渦從翼型上表面脫落，而當(dāng)翼型擺動到負(fù)攻角時，前緣渦從翼型下表面脫落。另外可以注意到，即使是從下表面脫落的前緣渦，依然保持著沿著翼型表面向下游流動的趨勢，沒有出現(xiàn)聚集和卷起。

(2) αmax=30°

攻角振蕩幅度增大到30°的流場與αmax=5°的流場有了明顯的不同。一個振蕩周期內(nèi)的點渦分布及演化，如圖6所示。

圖6 αmax=30°時的點渦分布及演化

從模擬出來的結(jié)果中可以看到，在攻角較小的時候，流場與αmax=5°的流場比較類似，尾緣渦逆時針卷起，而前緣渦則順著翼型上表面向下游運動。隨著攻角逐漸增大，前緣渦先于尾緣渦開始聚集，順時針卷起，在翼型表面上方形成了一個巨大的渦團，而后受到前緣渦的影響，尾緣渦不再直接向下游運動，而是在尾緣處逆時針卷起。尾緣渦聚集而成的渦團逐漸發(fā)展壯大，然后切斷了前緣渦的渦團。而后被切斷的前緣渦團又切斷了尾緣渦團，形成了一個蘋果型的渦對，如此交叉往復(fù)，就在流動中生成了一個個較大的渦對形成的渦街結(jié)構(gòu)?？梢钥吹?，在翼型攻角逐漸變小的過程中，渦對的大小明顯變小了，而且脫落的頻率更高。而進入后半個周期后，攻角變?yōu)樨?fù)值，前緣渦開始從翼型下表面脫落，并隨著負(fù)攻角的變大而出現(xiàn)了逆時針旋轉(zhuǎn)的渦團，相應(yīng)的，尾緣渦形成了順時針旋轉(zhuǎn)的渦團，與前半個周期的運動類似，前緣渦和尾緣渦的渦團互相夾斷，成對脫落。

3.3 俯仰振蕩翼型的升力分析

本文主要關(guān)心的是翼型的升力，利用渦量矩定理，對兩種不同攻角振幅下的算例進行了力的計算，首先是翼型的總升力系數(shù)。在較長周期的模擬中，總升力系數(shù)都表現(xiàn)出了周期性，這與翼型周期性振蕩的運動本質(zhì)是相符合的。翼型振蕩一個周期的升力系數(shù)時程曲線，如圖7所示。

圖7 翼型振蕩一個周期的升力系數(shù)時程曲線

從圖7中可以看到5°的攻角振動幅度下，升力系數(shù)很小，而30°的攻角振動幅度下，升力系數(shù)出現(xiàn)了明顯的起伏，基本呈現(xiàn)在前半個周期大于零，后半個周期小于零的形式，而且曲線類似于正弦函數(shù)的形狀，升力系數(shù)的極值出現(xiàn)在大約四分之一周期和四分之三周期處，此時翼型的攻角處于最大值，但是角速度為最小值零。30°的攻角振動幅度下，升力系數(shù)先達到一個最大值，然后突然變小。這個最大值點出現(xiàn)在t=4s左右，回顧圖6，發(fā)現(xiàn)此時前緣渦在翼型表面聚集形成了巨大的前緣附著渦團。隨后，前緣渦團與翼型表面發(fā)生了分離，導(dǎo)致升力突然急劇下降。下降的最低點對應(yīng)的時間點，流動逐漸從翼型表面分離，流場中尾緣渦卷起，并有切斷前緣渦的趨勢，翼型上表面出現(xiàn)向前運動的尾緣點渦，與自由流動方向相反，此時翼型上表面存在逆流，之后流動再次附著到翼型表面，升力開始增長。故認(rèn)為是翼型表面的前緣渦脫落和尾緣渦的逆向流動導(dǎo)致了升力下降，而前緣渦的再次附著又導(dǎo)致了升力重新增加。

4 總結(jié)

利用邊界層渦量和分離點處渦量的對應(yīng)關(guān)系，推導(dǎo)出了前緣渦脫落位置及環(huán)量的計算公式，對離散渦方法進行了修正。用修正后的離散渦方法對俯仰振蕩翼型的非定常流場進行了數(shù)值模擬，分析了高低兩種攻角振幅對應(yīng)的前緣渦和尾緣渦生長及脫落過程，結(jié)合流場對翼型所受升力的演化過程進行了定量化研究，主要結(jié)論有：

(1) 攻角振幅越大，升力系數(shù)的極值越大，通常出現(xiàn)在攻角最大值附近。

(2) 在一個振蕩周期內(nèi)，對于前后緣渦來說，前半個周期前緣渦的影響較大，后半個周期尾緣渦的影響較大。

(3) 大攻角下流動分離更加劇烈，導(dǎo)致升力系數(shù)會出現(xiàn)驟降。

[1] Shyy W, Aono H, Chimakurthi S K, et al. Recent progress in flapping wing aerodynamics and aeroelasticity[J]. Progress in Aerospace Sciences, 2010, 46(7): 284-327.

[2] Weis-Fogh T. Quick estimates of flight fitness in hovering animals, including novel mechanisms for lift production[J]. Journal of Experimental Biology, 1973, 59(1): 169-230.

[3] Triantafyllou M S, Triantafyllou G S, Gopalkrishnan R. Wake mechanics for thrust generation in oscillating foils[J]. Physics of Fluids A: Fluid Dynamics, 1991, 3(12): 2835-2837.

[4] Triantafyllou G S, Triantafyllou M S, Grosenbaugh M A. Optimal thrust development in oscillating foils with application to fish propulsion[J]. Journal of Fluids and Structures, 1993, 7(2): 205-224.

[5] Anderson J M. Vortex control for efficient propulsion[D]. Massachusetts Institute of Technology, Cambridge, MA, 1996.

[6] Guglielmini L, Blondeaux P. Numerical experiments on the transient motions of a flapping foil[J]. European Journal of Mechanics-B/Fluids, 2009, 28(1): 136-145.

[7] Zhi-Dong W, Xiao-Qing Z, Yu-Min S, et al. Analysis of the Caudal Vortices Evolvement around Flapping Foil[J]. Journal of Bionics Engineering, 2(04): 195-201.

[8] Lin H, Vezza M, Mcd. Galbraith R A. Discrete Vortex Method for Simulating Unsteady Flow[J]. AIAA Journal, 1997, 35(3): 494-499.

[9] Ginevsky A S, Zhelannikov A I. Discrete Vortex Method[J]. Vortex Wakes of Aircrafts, 2009: 13-32.

[10] Wu J C. Theory for aerodynamic force and moment in viscous flows[J]. AIAA Journal, 1981, 19(4): 432-441.

[11] Robertson D, Reich G, Joo J. Vortex particle aerodynamic modelling of perching manoeuvres with micro air vehicles[C]//51st AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference 18th AIAA/ASME/AHS Adaptive Structures Conference, 2010: 2825.