薛臣，周洲,*，李旭，許曉平

1.西北工業(yè)大學(xué) 航空學(xué)院，西安 710072 2.西北工業(yè)大學(xué) 無人機特種技術(shù)重點實驗室，西安 710065

隨著航空技術(shù)的發(fā)展，飛行器的應(yīng)用范圍變得越來越廣泛，因此面臨的問題也更加復(fù)雜。近幾年，人們越來越重視綠色航空，太陽能無人機由于其自身獨有的特點越來越受到人們的重視，如美國的太陽神[1-3]、英國的西風(fēng)[4-5]。此外，谷歌和Facebook也都宣布了各自的太陽能無人機計劃。太陽能無人機因為受限于目前能源利用率較低的因素，一般展弦比很大，呈現(xiàn)超輕面密度的結(jié)構(gòu)特點，所以太陽能無人機在飛行中會面臨較嚴重的振動問題，研究這一問題對指導(dǎo)太陽能無人機設(shè)計有重要意義。

在翼型振動方面，目前多數(shù)學(xué)者的研究重點都在大迎角運動和增升機理方面。大迎角俯仰振蕩研究方面現(xiàn)在主要集中在翼型動態(tài)失速研究，Albertson等[6]、Helin和Walker[7]通過實驗發(fā)現(xiàn)翼型前緣加速側(cè)會出現(xiàn)動態(tài)失速渦(文獻[8-9]中也叫前緣渦)，在低雷諾數(shù)條件下，這個渦會顯著提高升力。文獻[10-11]也指出對于撲翼來說，這個動態(tài)失速渦會在增升方面扮演著重要的角色；Leknys等[12]采用粒子圖像測速(PIV)對NACA0021翼型的動態(tài)失速問題進行研究，發(fā)現(xiàn)翼型的轉(zhuǎn)動造成了翼型的升力有一個明顯的延遲現(xiàn)象，而且可以將翼型的整體分離流動出現(xiàn)的角度推遲到60°。

注意到在翼型的大幅振蕩研究中，已經(jīng)發(fā)現(xiàn)翼型的振動會推遲翼型的失速迎角，這種由翼型振動帶來的遲滯現(xiàn)象也吸引了很多學(xué)者的注意。實際上，這種非定常效應(yīng)引起的遲滯效應(yīng)與很多因素有關(guān)，例如振動頻率，振動軸位置和振動幅度等[13]，Jung和Park[14]發(fā)現(xiàn)在低雷諾數(shù)條件下，振動減縮頻率的增大會減小尾渦的脫落頻率； Yu等[15]研究了前緣渦對翼型的升力阻力和力矩系數(shù)的影響，發(fā)現(xiàn)隨著減縮頻率k的增大，翼型的最大升力系數(shù)和相應(yīng)迎角逐漸增大。Qin等[16]在低雷諾數(shù)情況下通過大渦模擬計算了不同減縮頻率下的SD7003翼型的氣動特性，發(fā)現(xiàn)減縮頻率k的增大會使升力系數(shù)隨時間變化(CL-t)曲線的相位變大；Yu和Bernal[17]發(fā)現(xiàn)在相同的減縮頻率情況下，振動軸位置的后移會延遲翼型周圍渦的發(fā)展，但是主要的流動結(jié)構(gòu)并沒有改變。

而在小振幅振動研究中，目前所見的研究并不多，高正紅[13]在無黏跨聲速條件下，研究了振動軸位置和減縮頻率對小振幅振動翼型氣動特性的影響，研究發(fā)現(xiàn)小振幅簡諧振動翼型的升力系數(shù)-迎角(CL-α)曲線會呈現(xiàn)一個橢圓狀的環(huán)型曲線，指出減縮頻率的變大會使得CL-α曲線的繞向從逆時針變?yōu)轫槙r針，而且當振動軸位置后移時，發(fā)生這種變向所對應(yīng)的減縮頻率k值越大；Ekaterinaris和Platzer[18]利用面元法計算了NACA0012翼型的強迫振動，迎角幅度為±10°(未失速)，發(fā)現(xiàn)翼型振動會帶來相位上的改變，而且減縮頻率k的增大會使得力矩系數(shù)遲滯環(huán)從逆時針變?yōu)轫槙r針。

雖然已經(jīng)有不少學(xué)者研究了翼型振動帶來的一些特殊現(xiàn)象，但是這其中仍有很多問題并沒有得到很好的解釋，目前大多學(xué)者都把遲滯環(huán)方向的改變簡單歸結(jié)為非定常效應(yīng)引起的，張正秋等[19]指出，升力與扭矩的遲滯效應(yīng)隨各參數(shù)的變化趨勢相同，但是從文獻[13]的計算結(jié)果可以看出，升力系數(shù)遲滯環(huán)的變化和力矩系數(shù)遲滯環(huán)的變化是不相同的，并且文獻[20]中指出“非失速情況下，氣動力系數(shù)的遲滯環(huán)面積表征了遲滯效應(yīng)的大小”，注意到文獻[20]是在力矩系數(shù)變化情況下得到的結(jié)論，這一結(jié)論是否能推廣到升力系數(shù)遲滯環(huán)還有待研究。

可以看出目前專門針對翼型小振幅振動問題展開研究的并不是很多，而且多基于跨聲速情況，會有激波出現(xiàn)，對于低速情況，這種翼型小振幅振動會不會出現(xiàn)類似文獻[13]里的情況，有待研究。本文在前人的基礎(chǔ)上，結(jié)合本課題組在太陽能無人機實際飛行中遇到的特定問題，在低雷諾數(shù)條件下，從翼型小振幅振動入手研究振動軸位置和減縮頻率對振動翼型的氣動特性影響。

1 計算模型

本文采用的翼型為常用的低雷諾數(shù)翼型SD7037，網(wǎng)格劃分如圖1所示。Menter[21]指出，為準確地捕捉層流和轉(zhuǎn)捩邊界層，在低雷諾數(shù)流動計算中，需使y+接近于1。本文采用結(jié)構(gòu)化的O型二維網(wǎng)格，翼型弦長為c，計算域邊界半徑為25c，整個計算域及壁面附近網(wǎng)格如圖1所示。

圖1 振動翼型的計算網(wǎng)格Fig.1 Calculation grid of oscillation airfoil

2 數(shù)值方法

當弦長雷諾數(shù)在2×105～1×106的范圍內(nèi)變化時,基于雷諾時均的Navier-Stokes(RANS)方程也可以較好反映出翼型周圍的流場結(jié)構(gòu)[18,22]。本文采用滑移網(wǎng)格方法實現(xiàn)翼型的振動，采用k-kl-ω轉(zhuǎn)捩模型，因為S-A(Spalart-Allmaras)及SST(Shear Stress Transport)k-ω湍流模型主要考慮的是全湍區(qū)域，很難預(yù)測低雷諾數(shù)下的分離泡結(jié)構(gòu)，而k-kl-ω轉(zhuǎn)捩模型可以相對準確地捕捉觀察到翼型表面復(fù)雜流動現(xiàn)象及特征，轉(zhuǎn)捩位置及分離泡預(yù)測能力及精度相對較高，在低雷諾數(shù)轉(zhuǎn)捩流動數(shù)值求解中適用性較好[23]。

為了驗證本文所用轉(zhuǎn)捩模型的準確性，采用算例1的實驗數(shù)據(jù)[24]，驗證了k-kl-ω轉(zhuǎn)捩模型捕捉低雷諾數(shù)分離泡特征的能力，結(jié)果如圖2所示(圖中Cp為壓力系數(shù))。

算例1：翼型GA(W)-1，雷諾數(shù)Re=1.6×105，迎角α=12°。

本文中翼型的運動規(guī)律定義為

α=αmsin(ωt)

(1)

式中:αm為迎角運動幅度，本文中取為4°;ω為運動圓頻率，與減縮頻率k的關(guān)系為

(2)

其中：U∞為來流速度。

為了加快收斂，本文先利用定常計算得到一個較為穩(wěn)定的流場作為非定常運動計算的初始輸入，繼續(xù)進行非定常運動計算，每個時間步內(nèi)的子迭代步數(shù)選為30，時間步長按無量綱時間t/T=0.004 5選取，其中T為翼型運動周期。

為了驗證本文采用的計算方法對非定常運動計算的可靠性，采用了2個算例來驗證。

算例2：利用Piziali[25]的實驗結(jié)果進行驗證，翼型為NACA0015，馬赫數(shù)Ma=0.29，Re=1.95×106，運動規(guī)律為α=4°+4.2°sin(2πft)，其中頻率f=10 Hz。

圖2 GA(W)-1翼型壓力分布(α=12°)Fig.2 Pressure distribution of GA (W)-1 airfoil (α=12°)

算例3：翼型為NACA0012，Re=485 755，Ma=0.035 5，運動規(guī)律為α=10°+10°sin(2πft)，其中f=1.62 Hz，實驗條件及數(shù)據(jù)可以參考文獻[26]。

從圖3和圖4給出的升力系數(shù)遲滯環(huán)曲線中可以看出，本文采用的方法與實驗結(jié)果符合較好，說明本文采用的計算方法有足夠的可靠性。

圖3 NACA0015翼型升力系數(shù)遲滯環(huán)Fig.3 Hysteresis loop of lift coefficient of NACA0015 airfoil

圖4 NACA0012翼型升力系數(shù)遲滯環(huán)Fig.4 Hysteresis loop of lift coefficient of NACA0012 airfoil

3 結(jié)果及討論

通過計算可以發(fā)現(xiàn)，當翼型作俯仰振動時，CL-α曲線呈現(xiàn)為一個橢圓形，而且方向上也有變化。為了研究振動軸位置和減縮頻率k對CL-α曲線遲滯環(huán)的繞向的影響，本文首先在減縮頻率不變的情況下，研究振動軸位置對遲滯環(huán)的影響，隨后通過更改減縮頻率k，研究k對遲滯環(huán)繞向的影響，并分析原因。

3.1 振動軸的影響

本文中定義了2個坐標系，一個是隨著平板運動的體坐標系oxy，另一個是固定的慣性系OXY(風(fēng)軸系)，當翼型未振動時，這2個坐標系是重合的，如圖5所示。分別將振動軸定位于前緣處、1/4弦線、半弦長、后緣處，即x=0,0.25c,0.50c,c處(見圖5)。

這里僅展現(xiàn)k=0.3和k=1.0時的計算結(jié)果分別如圖6和圖7所示。從圖6(a)可以看出，當k=0.3，翼型繞前緣轉(zhuǎn)動時，此時的CL-α曲線“遲滯環(huán)”的繞向為順時針，對比圖6(d)可以發(fā)現(xiàn)此時的“遲滯環(huán)”方向變?yōu)槟鏁r針，即在k保持不變的情況下，當振動軸向后移動時，“遲滯環(huán)”方向有從順時針方向變?yōu)槟鏁r針的趨勢，而且在這個變化過程中，會經(jīng)歷一個“直線”過程，如圖6(c)所示。當k=1.0時，可以看出當振動軸從前緣移動到半弦長處時，對應(yīng)的遲滯環(huán)均為順時針，但是當振動軸移動到后緣處時，CL-α曲線此時也出現(xiàn)了近似直線的情況，如圖7(d)所示。

圖5 不同振動軸在坐標系中的示意圖Fig.5 Sketch of different pivot locations in coordinate system

為了弄清楚這種現(xiàn)象的原因，將升力系數(shù)隨時間變化曲線展開，如圖8所示。從圖中可以看出k=0.3，振動軸位置從前緣向后移動時，升力系數(shù)的峰(谷)值在時間上均出現(xiàn)了向后移動。

表1給出了不同雷諾數(shù)下對應(yīng)不同振動位置的升力系數(shù)曲線特性，表1中的Δt和φ是相對于α-t曲線計算得到的時間差和相位差，負值代表相位相對α-t曲線(時間)延遲，正值代表相位提前；CLmax和CLmin分別為升力系數(shù)的最大值和最小值，CLamp為振動時升力系數(shù)的幅值大小，CLmean為

圖6 k=0.3時不同振動位置的CL-α曲線(Re=7.89×105)Fig.6 CL-α curves at different pivot locations at k=0.3 (Re=7.89×105)

圖7 k=1.0時不同振動位置的CL-α曲線(Re=7.89×105)Fig.7 CL-α curves at different pivot locations at k=1.0 (Re=7.89×105)

圖8 k=0.3時不同振動位置的CL-t曲線(Re=7.89×105)Fig.8 CL-t curves at different pivot locations at k=0.3 (Re=7.89×105)

升力系數(shù)的平均值。從表1中可以看出，振動軸的位置對升力的平均值影響不大，而且振動軸的位置越靠近后緣，升力系數(shù)能達到的最大值越小；對于同一振動軸位置，隨著雷諾數(shù)的變大，時間和相位在數(shù)值上是逐漸增大的，雖然對于x=c位置來說，相位為負值，代表此時CL-t曲線相對于α-t曲線相位滯后，但是隨著雷諾數(shù)的增大，其滯后相位的絕對值也在變小，即當雷諾數(shù)變大時，會使得相位提前變多。

此外，從表1中還可以看出Re=7.89×105時，x=0位置處的相位較迎角曲線的相位要提前19.22°，x=0.25c位置相位提前10.43°，而對于x=c位置，相位滯后11.52°,這種相位上的差別就導(dǎo)致了遲滯環(huán)曲線的方向變化，這里可以從數(shù)學(xué)上加以證明。

表1 不同雷諾數(shù)下對應(yīng)不同振動位置的升力系數(shù)曲線特性(k=0.3)

迎角的運動形式為α(t)=αmsin(ωt)，穩(wěn)定后的升力系數(shù)曲線可以寫為：CL(t)=CLampsin(ωt+φ)，這點在3.2節(jié)可以被證明。

為了方便說明，以前1/4周期為例(其他情況類似)，t1(t1T/4)時刻剛好使得ωt2=θ2=π-θ1，此時α2=α1，有CL(t2)=CLampsin(ωt2+φ)=CLampsin(π-θ1+φ)，容易得到：φ>0°時，CL(t1)>CLamp(t2)，反映到CL-α曲線上，即為順時針方向；φ<0°時，CL(t1)

從數(shù)學(xué)關(guān)系上的證明可以看出，在一定范圍內(nèi)，相位φ的絕對值越大，t1、t2時刻對應(yīng)的升力系數(shù)差別越大，即遲滯環(huán)表現(xiàn)得會更加飽滿(無論是順時針還是逆時針)，而φ又與振動軸位置和k相關(guān)，所以這才是導(dǎo)致遲滯環(huán)變化的根本原因，而不是簡單的增大k就會導(dǎo)致遲滯環(huán)飽滿，對比圖7(d)和圖6(d)，可以看出增大k也可能使得遲滯環(huán)曲線變得扁平，甚至為直線；從圖6和圖7中還可以看出，向前移動振動軸位置也可能使得遲滯環(huán)變得更加飽滿。

圖9 相位差φ帶來的升力系數(shù)差Fig.9 Lift coefficient difference caused by phase difference φ

從以上分析可以看出，當振動軸位置向后緣移動時，CL-t曲線的相位會變小，即出現(xiàn)相位延遲，CL-α曲線也會從順時針向逆時針變化。對于此現(xiàn)象，可以有如下解釋：以平板為例，振動軸位置定在離前緣d處，如圖10所示。假設(shè)整個平板在來流速度為U∞的流場里以ω*的角速度繞振動軸順時針旋轉(zhuǎn)(即迎角增大)，則在體軸系里面觀察，前緣處會有一個vLE=ω*d的向上速度，如圖10所示，而當迎角為αoff時，此時在體軸系里觀察來流，在y方向上的速度為vy=U∞sinαoff，但是注意到此時前緣自身仍有一個vLE的向上速度，所以前緣處感受到的實際速度只有vLE-real=U∞sinαoff-vLE，也就是說實際迎角被減小了，而且只有當αoff=arcsin(ω*d/U∞)≈ω*d/U∞(小角度)才能剛好抵消前緣運動速度的影響，此時在體軸系觀察，來流沿x方向流過，迎角為0°，即相比較繞前緣振動，時間上滯后了。

這里雖然只對前緣處作了分析，但是其他位置處的有效迎角是相似的，在振動軸前的部分有效迎角是減小的，振動軸后處的部分有效迎角在增大，這也可以看出當d較大時，即振動軸離前緣的位置越遠，有效迎角減小的部分更多，要達到和繞前緣振動相同狀態(tài)則需要更長的時間，反映在升力系數(shù)上，就是振動軸的后移延遲了相位，如圖8所示。此外，還可以看出振動軸帶來的相位滯后在k一定的情況下，主要和離前緣的距離d和當?shù)氐钠骄魉賃av的比值近似呈線性關(guān)系。

圖11中展現(xiàn)了翼型振動時沿X方向的速度vX云圖，從圖11中可以看出，當振動軸位于翼型后緣(x=c)時，基本上要比繞前緣處(x=0)振動滯后0.02 s，才會出現(xiàn)相似的流場結(jié)構(gòu)。

圖10 平板的速度示意圖Fig.10 Sketch of velocity of flat plate

圖11 繞前緣和后緣振動不同時刻的vX云圖對比(k=0.3)Fig.11 Comparison of vX contour at different times when pivot located at leading edge and trailing edge (k=0.3)

3.2 減縮頻率的影響

本節(jié)在相同振動軸位置，研究了不同k對遲滯環(huán)曲線的影響。圖12給出了x=0，0.25c，c時，k=0.05,0.15,0.30,1.00的計算結(jié)果。

從圖12中可以看到，在不同振動軸處，隨著k的變大，“遲滯環(huán)”的轉(zhuǎn)向會從逆時針向順時針轉(zhuǎn)變，中間會出現(xiàn)直線段過程，這并不如文獻[19]結(jié)論中指出的那樣，“隨著振動頻率的增加，氣動力系數(shù)遲滯曲線飽滿程度增加”?？梢钥闯鲈龃髃值類似于振動軸的前移過程，也可以提前相位。

從圖12中還可以看出，當k=1.00時，越靠近前緣處，升力系數(shù)遲滯環(huán)面積越大，這主要是因為當振動軸靠近前緣時，初始的負相位φ的絕對值較小，因此k增大到一個較小的值時，就會抵消這個負相位，此時升力系數(shù)遲滯環(huán)即變?yōu)橹本€型，繼續(xù)增大k這個相位φ就會變成正值，并不斷變大，按3.1節(jié)的解釋，遲滯環(huán)面積會變大，且此時的升力系數(shù)遲滯環(huán)為順時針轉(zhuǎn)向，而當振動軸位置靠近后緣時，這個初始的滯后相位的絕對值會較大，因此k需要增大到更大的水平才能抵消這個滯后相位，所以反映到遲滯環(huán)曲線上，如圖12，當k=1.00時，越靠近前緣的振動軸位置的曲線更加飽滿，當振動軸位于后緣時，CL-α曲線基本呈直線型，表明當振動軸位于后緣時，當k增大到1.00時，才能抵消相位滯后，而當振動軸位于前緣時，k只需到0.15，就可以抵消相位的滯后了。

圖13將不同k下的CL-t曲線與α-t曲線進行了對比，可以看出，隨著k的變大，CL-t曲線相對于此時對應(yīng)的α-t曲線的相位逐漸變大，即k的增大可以使得相位提前，此時對應(yīng)的CL-α曲線遲滯環(huán)也會從逆時針向順時針變化，直線段則對應(yīng)相位差為0°的情況。

圖12 不同振動軸位置隨k變化的CL-α曲線(Re=7.89×105)Fig.12 Variation of CL-α curves with k at different pivot locations (Re=7.89×105)

圖13 不同k的CL-t曲線與α-t曲線的對比(Re=7.89×105)Fig.13 Comparison between CL-t curves and α-t curves for different k (Re=7.89×105)

圖14給出了φ隨k的變化，可以清楚地看到對于同一振動軸，k的增大使得曲線的相位提前，而且曲線的相位是與減縮頻率k和振動軸有關(guān)的，此外可以看出相位與k之間存在一個函數(shù)關(guān)系，這點也在后文中進行了證明。

圖15給出了SD7037翼型在不同雷諾數(shù)下的小迎角升力系數(shù)曲線，可以看出在氣動上，本文選取的翼型在小迎角范圍內(nèi)，升力系數(shù)和迎角近似呈線性關(guān)系CL=CLαα(CLα為升力線斜率)，現(xiàn)在將翼型及流場之間的共同作用看成一個系統(tǒng)，升力系數(shù)是輸出量，很明顯在小迎角范圍內(nèi)此時這個系統(tǒng)滿足線性系統(tǒng)的定義，本質(zhì)上是線性的。從控制的角度來看，一個系統(tǒng)的線性與否是其自身屬性，與輸入量沒有關(guān)系，雖然這里的輸入為迎角的變化：α=αmsin(ωt)，但是不會改變這個系統(tǒng)的線性特點，所以本文在這里將此系統(tǒng)看作是線性系統(tǒng)。另外從圖8中也可以看出，升力系數(shù)CL此時也是正弦變化形式，而且頻率保持與迎角變化頻率一致，也從側(cè)面驗證了這種假設(shè)的合理性。在線性系統(tǒng)假設(shè)的前提下，輸入(迎角)和輸出(升力系數(shù)、力矩系數(shù))的關(guān)系可證明如下。

圖14 CL-t曲線與α-t曲線之間的相位差φ隨k的變化Fig.14 Variation of phase shift φ between CL-t curve and α-t curve with k

圖15 SD7037翼型在不同雷諾數(shù)下的小迎角升力系數(shù)曲線Fig.15 Lift coefficient curves of SD7037 airfoil at small angles of attack with different Reynolds numbers

系統(tǒng)的輸入為α(t)=αmsin(ωt)，則輸出的拉氏變換可以寫為

CL(s)=Φ(s)αm(s)

(3)

式中：

Φ(s)為整個系統(tǒng)的傳遞函數(shù)；M(s)為Φ(s)的分子多項式；-s1,-s2,…,-sn為系統(tǒng)極點。為了方便討論，不妨認為系統(tǒng)穩(wěn)定，并假設(shè)特征根均為互異的負實根，利用部分分式展開法即有：

(4)

A2e-s2t+…+Ane-snt

(5)

當t趨于無窮時，輸出穩(wěn)態(tài)分量：

(6)

所以穩(wěn)態(tài)分量為

CL(t)|t→∞=αm|Φ(jω)|·

αm|Φ(jω)|sin(ωt+∠Φ(jω))=

CLamp(ω)sin(ωt+φ(ω))

(7)

從證明中可以看出，振幅CLamp和相位φ都是與ω有關(guān)的函數(shù)，這也驗證了圖14所示的規(guī)律。

結(jié)合以上證明，如果小振幅振動的翼型可以被看成一個線性系統(tǒng)，則可做一個合理的假設(shè)，即當雷諾數(shù)和振動軸位置確定后，運動的初始相位不會影響最終的CL-t曲線的相位。為了驗證這一假設(shè)，本文選取k=0.30、振動軸位置為x=0.25c和x=c處進行分析，從表1可以看出，振動軸位于0.25c處時，CL-t曲線相對α-t曲線的相位本身就是提前的，而振動軸位于后緣時，這個相位是負值，即相位滯后了。為了驗證初始相位對最終CL-t曲線相位的影響，在x=0.25c處預(yù)加-18°的初始相位，即α=αmsin(ωt-18°)，相應(yīng)的在x=c處預(yù)加+18°的初始相位，即α=αmsin(ωt+18°)。計算結(jié)果如表2所示，對比表1中的數(shù)據(jù)，可以發(fā)現(xiàn)初始相位的存在并沒有改變最終的相位特性，即CL-t曲線的相位不受迎角的初始相位影響。

表2 改變初始相位后的CL-t曲線相位特性

Table 2 Phase characteristics ofCL-tcurve after changing initial phase

在3.1節(jié)的分析中，已清楚了振動軸位置主要是通過影響有效迎角進而影響相位的變化,接下來研究k的變化對相位的影響。

注意到本文研究的是翼型的簡諧強迫振動，所以這里參考了二維Theodorsen非定常氣動力模型，具體表達式[27-28]為

將本文的翼型運動規(guī)律帶入化簡，可以得到升力系數(shù)表達式：

CL=πkαm(cos(ωt)+aksin(ωt))+…+

(9)

參考文獻[17,29]，本文在這里把翼型排開空氣的質(zhì)量稱為“附加質(zhì)量”，而且文獻[17,29]的研究結(jié)果都指出：“當k較大時，翼型排開的空氣質(zhì)量帶來的慣性反作用力會帶來額外的升力，從而提前相位”,而且可以看出當k較大時，帶來的額外升力也會比較明顯。

從式(9)中還可以看出，這個附加質(zhì)量帶來的慣性反作用力只會對升力的幅值和相位有所影響，并不影響CL-t曲線的頻率特性，這與計算結(jié)果也是相符合的。

此外，從圖14中可以看出，當振動軸位置向后緣移動時，k增大帶來的相位增加并不明顯，即升力系數(shù)曲線的相位對k不再敏感了，這與文獻[13] 的結(jié)論相似，從本節(jié)的分析可以看出，當k變大時，翼型排開的附加質(zhì)量帶來的慣性反作用力會變大，使得翼型升力變大。從式(9)可以得到在周期開始時(瞬時迎角為0°)，由于反作用力帶來的升力系數(shù)增量為ΔCL=πkαm，所以當k=1.00 時，由于反作用力帶來的升力系數(shù)增量ΔCL1-t=0.219 3，k=0.30時，升力系數(shù)增量ΔCL2-t=0.065 8，即k=1.00相比于k=0.30，升力系數(shù)差值為ΔCL-t=0.153 5，從圖13中可以得到繞前緣振動時，k=1.00和k=0.30在周期開始時升力系數(shù)差值為ΔCL1-c=0.221 1，而繞后緣振動時，對應(yīng)的升力系數(shù)差值ΔCL2-c=0.116 6，此時有ΔCL1-c>ΔCL-t>ΔCL2-c，可以看出當振動軸后移時，升力增量的確是減小的，這主要是因為當振動軸在前緣時，翼型在迎角增大的過程中，有效迎角也在變大，如圖16(a)和圖16(b)所示，即此時升力的2個組成部分，附加質(zhì)量帶來的反作用力和常規(guī)的氣動力都是增加的，因此相位提前很顯著，而當振動軸后移時，特別地，對于振動軸在后緣時，翼型上的有效迎角是減小的，如圖16(c)和圖16(d) 所示，雖然由于附加質(zhì)量帶來的反作用力隨k的增大在增大，但是翼型上產(chǎn)生的氣動力由于有效迎角的變小在減少，而且振動軸后移本身會帶來負相位，這就大大抵消了慣性力帶來的額外升力增量，表現(xiàn)出來的就是k增大，相位增加緩慢。

圖16 繞前緣和后緣振動時翼型周圍流場隨k的變化(Re=7.89×105)Fig.16 Variations of flow field around airfoil with k when pivot located at leading edge and trailing edge (Re=7.89×105)

3.3 推論及應(yīng)用

從3.1節(jié)和3.2節(jié)的分析可知，利用已知的計算結(jié)果得到繞其他振動軸位置振動的升力系數(shù)曲線(幅值會略有差別)是可行的，例如，當k=0.30，Re=7.89×105時，從表1可以看出，繞x=0.50c處振動時CL-t曲線的相位偏移為4.87°，即時間上提前0.003 s，而繞x=c處振動時，相位和時間分別為-11.52°和-0.007 1 s(負號代表滯后)，根據(jù)本文的結(jié)論，在這兩個振動軸位置之間必然會出現(xiàn)相位為0的情況，即CL-α曲線呈現(xiàn)直線型，而且可以根據(jù)這兩個狀態(tài)推算出相位為0°時對應(yīng)的振動軸位置。首先，從這兩個狀態(tài)可以推算出x=0.50c到后緣處這一段翼型表面的平均流動速度為

Uav=(c-0.5c)/(0.003+0.007 1)=49.5cm/s

假設(shè)振動軸位置位于x0處時，相位差為0°，即時間差也為0，從3.1節(jié)的分析可以看出振動軸位置對相位的影響主要是線性的，所以有

x0-0.5c=Uav×0.003≈0.15c

即這個振動軸位置位于x=0.65c時，CL-α曲線近似呈現(xiàn)直線型。

圖17 k=0.30，x=0.65c時的CL-α曲線(Re=7.89×105)Fig.17 CL-α curve at k=0.30, x=0.65c (Re=7.89×105)

k=0.30,x=0.65c時的CL-α曲線如圖17所示，驗證了推算的準確性。相應(yīng)地，如果想推算已知振動軸位置的升力系數(shù)曲線，可以利用同樣的方法，推算出此振動位置的CL-t曲線與α-t曲線的相位差，如當k=0.30時，已知振動軸位置x=0.75c，利用x=0.50c和x=c的結(jié)果可以推算出其時間上滯后0.002 s，即相位差應(yīng)為-3.23°，則升力系數(shù)曲線為

CL-x=0.75c=CLmean+CLampsin(ωt-3.23°)

(10)

式中:CLmean和CLamp可以由表1插值得到，這里分別為0.362和0.297。

如圖18所示，推算的結(jié)果和計算的結(jié)果符合良好，可以看出本文得出的結(jié)論是正確可靠的，這對于工程實際應(yīng)用有較大的意義，例如可以通過計算或?qū)嶒灚@得容易測量位置的結(jié)果，以此推算出其他位置的結(jié)果。

此外，根據(jù)本文的分析，小振幅振動的翼型并沒有牽涉到渦的形成和發(fā)展，這也意味著傳統(tǒng)的勢流理論一定程度上也可以預(yù)測這種遲滯環(huán)變向的現(xiàn)象。

圖18 k=0.30, x=0.75c時的預(yù)測曲線和計算曲線對比Fig.18 Comparison between predicted curve and calculated curve at k=0.30, x=0.75c

本文利用面元法[30]計算了相應(yīng)的狀態(tài)，結(jié)果如圖19所示，可以看出當k<0.30時，面元法預(yù)測的結(jié)果和基于RANS方程計算的結(jié)果高度一致，但是當k=1.00時，升力系數(shù)計算結(jié)果不僅在幅值上有較大差別，而且相位預(yù)測也有差異。這主要是因為傳統(tǒng)的勢流理論只能處理附著流，而從圖16中可以看出，當k增大到一定程度時，振動翼型周圍流場特別是前緣處的有效迎角被改變，因此面元法計算結(jié)果也和CFD結(jié)果相差較大。但是，面元法在一定程度上還是提供了一種更加快速的研究小振幅振動翼型氣動特性的方法。

圖19 面元法與CFD動網(wǎng)格計算結(jié)果對比(Re=7.89×105)Fig.19 Comparison of calculated results of panel method and CFD sliding mesh (Re=7.89×105)

4 結(jié) 論

1) 遲滯環(huán)現(xiàn)象主要是升力系數(shù)、力矩系數(shù)隨時間變化的曲線在相位上與迎角變化規(guī)律曲線的差異，當升力系數(shù)曲線與迎角隨時間變化曲線的相位差φ>0°，即升力曲線相位提前時，升力系數(shù)隨迎角變化的遲滯環(huán)曲線呈現(xiàn)順時針變化規(guī)律，而當φ<0°，即升力曲線相位滯后時，遲滯環(huán)曲線呈現(xiàn)逆時針變化，當φ=0°時，此時的升力曲線恰好為直線型，而這個相位與振動軸位置和減縮頻率k有著密切關(guān)系。

2) 振動軸的影響：本文的計算結(jié)果顯示，振動軸位置后移會帶來相位上的延遲，而且這種影響近似呈線性關(guān)系。相位上的延遲，使得遲滯環(huán)向順時針變化就越困難，即相同的k下，振動軸的位置越靠后，升力系數(shù)遲滯環(huán)越可能呈現(xiàn)逆時針變化。

3) 減縮頻率k的影響：減縮頻率k對升力系數(shù)曲線的影響機制與振動軸不同，k的增大，帶動周圍空氣振動從而帶來的慣性反作用力增大，進而提高了升力，即提前了相位，但是并不改變CL-t曲線的頻率?？梢钥闯鲞@種變化不再類似振動軸位置的線性影響。

4) 雷諾數(shù)的影響：本文研究結(jié)果表明，當雷諾數(shù)變大時，流動黏性會相對變?nèi)?，因此相位滯后比低雷諾數(shù)時較小，但是不會使得遲滯環(huán)現(xiàn)象消失，文獻[13]在無黏情況下(可以看成雷諾數(shù)無窮大)也觀察到了遲滯環(huán)現(xiàn)象，與本文結(jié)論類似。