董圣華，史愛(ài)明，葉正寅，田海濤

(西北工業(yè)大學(xué)航空學(xué)院，陜西 西安 710072)

超臨界翼型跨聲速抖振CFD計(jì)算和POD分析

董圣華，史愛(ài)明*，葉正寅，田海濤

(西北工業(yè)大學(xué)航空學(xué)院，陜西 西安 710072)

使用二階迎風(fēng)Roe格式、隱式時(shí)間推進(jìn)和S-A(Spalart-Allmaras)一方程湍流模型，通過(guò)求解基于格心格式有限體積法的RANS方程模擬了OAT15A超臨界翼型的跨聲速抖振流場(chǎng)。在模擬出激波/附面層相互干擾誘發(fā)的抖振現(xiàn)象的基礎(chǔ)上，對(duì)翼面激波運(yùn)動(dòng)過(guò)程中的氣流分離泡變化規(guī)律開展研究，揭示出激波變化相位略微超前于升力系數(shù)相位的非定?，F(xiàn)象。為剖析跨聲速抖振的物理成因，將抖振計(jì)算的非定常解作為快照，應(yīng)用本征正交分解POD(Proper Orthogonal Composition)方法提取POD模態(tài)，從流場(chǎng)相干結(jié)構(gòu)的角度對(duì)跨聲速抖振現(xiàn)象進(jìn)行了分析。

抖振;激波/附面層相互干擾;本征正交分解

0 引言

大型運(yùn)輸類飛機(jī)上多采用超臨界翼型，以達(dá)到增加阻力發(fā)散馬赫數(shù)的目的。但是在跨聲速范圍內(nèi)機(jī)翼有可能會(huì)產(chǎn)生由激波/附面層相互干擾而引起的跨聲速抖振現(xiàn)象，這將嚴(yán)重影響飛行器的飛行品質(zhì)與結(jié)構(gòu)壽命［1］，因此對(duì)于跨聲速抖振問(wèn)題的研究就顯得非常重要。

為了對(duì)由激波/附面層干擾而產(chǎn)生的跨聲速抖振現(xiàn)象進(jìn)行深入研究，法國(guó)國(guó)家航空航天科研局(ONERA)于2005年對(duì)OAT15A超臨界翼型的二維翼段進(jìn)行了風(fēng)洞實(shí)驗(yàn)［2-3］。由于 L.Jacquin［2-3］等使用OAT15A模型進(jìn)行的實(shí)驗(yàn)比以往同類實(shí)驗(yàn)提供了更加豐富詳盡的實(shí)驗(yàn)條件說(shuō)明與結(jié)果分析，所以很多研究者將其作為研究跨聲速激波抖振現(xiàn)象的標(biāo)準(zhǔn)算例［4-6］。

除了對(duì)于OAT15A翼型的研究之外，近期對(duì)于抖振現(xiàn)象的研究還包括新的抖振邊界預(yù)測(cè)方法［7-8］、抖振強(qiáng)度的控制方法［9-10］以及更接近于實(shí)際情況的抖振現(xiàn)象的研究［11］。

本征正交分解POD方法是一種將數(shù)據(jù)進(jìn)行壓縮的多變量統(tǒng)計(jì)分析方法，其基本思想是根據(jù)已有樣本數(shù)據(jù)計(jì)算得到一組最能夠代表這組數(shù)據(jù)的正交基函數(shù)，即使得樣本數(shù)據(jù)投影到此正交基上所產(chǎn)生的誤差最小。在構(gòu)造這組基函數(shù)時(shí)使得樣本數(shù)據(jù)在正交基上的投影分量按次序依次迅速衰減［12］，所以它可以將高維數(shù)據(jù)降階投影到低維空間，從而獲得數(shù)據(jù)的物理特征［13］。該方法已經(jīng)被廣泛應(yīng)用于建立流場(chǎng)的降階模型［14］、翼型的反設(shè)計(jì)［15］以及流動(dòng)結(jié)構(gòu)分析［16-18］等方面。

本文使用結(jié)構(gòu)化網(wǎng)格，針對(duì)風(fēng)洞實(shí)驗(yàn)狀態(tài)條件進(jìn)行了OAT15A超臨界二維翼型繞流的抖振問(wèn)題CFD數(shù)值方法研究。在比較了不同網(wǎng)格分布對(duì)于計(jì)算結(jié)果的影響之后，選擇其中一套網(wǎng)格的計(jì)算結(jié)果進(jìn)行了相應(yīng)的后處理并深入的研究了抖振現(xiàn)象的發(fā)展過(guò)程，發(fā)現(xiàn)激波運(yùn)動(dòng)的相位與升力系數(shù)波動(dòng)的相位并不完全相同。接下來(lái)提取了流場(chǎng)結(jié)果的POD模態(tài)，從物理流動(dòng)相干流場(chǎng)結(jié)構(gòu)角度對(duì)抖振現(xiàn)象進(jìn)行了分析，尋找抖振現(xiàn)象的主導(dǎo)模態(tài)。

1 跨聲速抖振流場(chǎng)的模擬及分析

1.1 算例描述及數(shù)值方法

L.Jacquin［2-3］等人對(duì)于OAT15A翼型在馬赫數(shù)為0.73，來(lái)流迎角為3.5°，雷諾數(shù)為3×106的實(shí)驗(yàn)條件下所獲得的結(jié)果進(jìn)行了比較全面的分析，因此選擇此實(shí)驗(yàn)條件作為數(shù)值模擬的條件。這個(gè)實(shí)驗(yàn)條件下的結(jié)果也是大多數(shù)研究者進(jìn)行數(shù)值模擬時(shí)所選用的對(duì)比結(jié)果［2-3，5-6］。

使用基于有限體積法的流場(chǎng)解算程序進(jìn)行數(shù)值模擬，空間離散格式為二階迎風(fēng)Roe格式［19］，時(shí)間項(xiàng)則采用全隱式雙時(shí)間格式［20］，物理時(shí)間步長(zhǎng)為6.6 ×10-5s，并使用當(dāng)?shù)貢r(shí)間步長(zhǎng)［21］來(lái)提高計(jì)算效率，湍流模型選用S-A一方程模型［22］。

1.2 網(wǎng)格收斂性練習(xí)

本文劃分了疏密不同的三套網(wǎng)格。三套網(wǎng)格均為O型網(wǎng)格，遠(yuǎn)場(chǎng)半徑為50倍參考弦長(zhǎng)，近壁面第一層網(wǎng)格滿足y+＜1，且翼型前緣節(jié)點(diǎn)分布也完全相同。這三套網(wǎng)格不同之處主要有兩點(diǎn):一是翼型鈍后緣處的網(wǎng)格節(jié)點(diǎn)數(shù)目依次為6、10、30;二是后緣之后的氣流分離區(qū)范圍之內(nèi)的網(wǎng)格的疏密不同，其中B網(wǎng)格的后緣網(wǎng)格最稀疏。第一套網(wǎng)格中翼型周圍的網(wǎng)格如圖1所示，三套網(wǎng)格的參數(shù)如表1。

圖1 A網(wǎng)格的網(wǎng)格拓?fù)浣Y(jié)構(gòu)Fig.1 Grid topology for Mesh A

表1 網(wǎng)格分布Table 1 Grid distributions of the meshs

使用這三套網(wǎng)格計(jì)算出的升力系數(shù)隨時(shí)間的變化如圖2所示。從圖中可以看到使用A與C網(wǎng)格計(jì)算出了周期性的升力系數(shù)波動(dòng)，其波動(dòng)的幅度相近，而B網(wǎng)格對(duì)應(yīng)的升力系數(shù)則最終收斂于一個(gè)定值。結(jié)合三套網(wǎng)格的不同之處來(lái)看，B網(wǎng)格無(wú)法計(jì)算出抖振現(xiàn)象的主要原因就是其在后緣之后弦向網(wǎng)格分布均勻的區(qū)域過(guò)小。

圖2 升力系數(shù)隨時(shí)間的變化Fig.2 Lift coefficient evolution versus time

截取A與C網(wǎng)格升力系數(shù)波動(dòng)的周期性部分并進(jìn)行傅里葉變換(FFT)，就得到了這兩套網(wǎng)格計(jì)算出的抖振頻率，分別為77.8 Hz與76 Hz，而實(shí)驗(yàn)值為69 Hz，C網(wǎng)格頻率的計(jì)算值更加接近于實(shí)驗(yàn)值。將計(jì)算得到的平均壓力系數(shù)、壓力均方根值與實(shí)驗(yàn)值以及文獻(xiàn)［6］中使用KKL湍流模型的計(jì)算值進(jìn)行對(duì)比，分別如圖3和圖4，可以看到A、C兩套網(wǎng)格的計(jì)算值幾乎重合。本文計(jì)算出的平均壓力系數(shù)在整個(gè)下表面以及上表面超聲速區(qū)與實(shí)驗(yàn)值吻合程度較好，而在激波運(yùn)動(dòng)區(qū)域與再壓縮區(qū)則與實(shí)驗(yàn)值則有一定偏差，文獻(xiàn)［5］使用KKL模型計(jì)算出的激波運(yùn)動(dòng)范圍與實(shí)驗(yàn)值相比整體靠后;均方根值的幅度高于實(shí)驗(yàn)值，但本文均方根值突增的起始位置與實(shí)驗(yàn)相同，使用KKL模型得到的波動(dòng)范圍較窄。

圖3 翼型表面平均壓力系數(shù)Fig.3 Mean surface pressure coefficient distributions

圖4 翼型上表面壓力均方根值分布Fig.4 Pressure rms distributions on the top side of OAT15A airfoil

激波的平均位置可以通過(guò)壓力脈動(dòng)偏斜因子［2］的分布得到，其定義如下:

式中，p'為壓力的脈動(dòng)量，即p'=p(t)－，p(t)為每一瞬時(shí)的壓力，為壓力的時(shí)均值。壓力脈動(dòng)偏斜因子考慮了壓力脈動(dòng)在一個(gè)周期之中的時(shí)間分布以及脈動(dòng)量的絕對(duì)大小。其值為正，說(shuō)明壓力在一個(gè)周期中的大部分時(shí)間內(nèi)要小于平均值;其值為負(fù)，說(shuō)明壓力在一個(gè)周期之中的大部分時(shí)間內(nèi)要大于壓力的平均值;這個(gè)值的絕對(duì)值越大，則表示在某個(gè)或某些時(shí)刻時(shí)壓力偏離壓力平均值的量越大。

上表面Sp的計(jì)算值與實(shí)驗(yàn)值沿弦向的分布如圖5。在激波運(yùn)動(dòng)區(qū)域之內(nèi)，Sp出現(xiàn)最大值，并在之后迅速衰減至最小值。這兩個(gè)值對(duì)應(yīng)于激波震蕩時(shí)上游與下游的極限位置，而其間Sp為零的點(diǎn)所對(duì)應(yīng)的弦向位置就是激波的平均位置。兩套網(wǎng)格的激波運(yùn)動(dòng)區(qū)域介于0.35c至0.6c之間，大于實(shí)驗(yàn)值測(cè)得的0.35c至0.5c的范圍。計(jì)算出的激波平均位置在0.5c處，比0.44c的實(shí)驗(yàn)值靠后。由于采用的URANS方法對(duì)于小脈動(dòng)量的模擬能力有所不足，造成壓力均方根和壓力脈動(dòng)偏斜因子的計(jì)算值與實(shí)驗(yàn)值有一定的差距。

綜上所述，選擇C網(wǎng)格的瞬時(shí)計(jì)算結(jié)果作為分析數(shù)據(jù)。

圖5 翼型上表面Sp沿弦向的分布Fig.5 Spdistributions on the top side of OAT15A airfoil

1.3 抖振流場(chǎng)的非定常結(jié)果分析

壓力脈動(dòng)偏斜因子可以描繪出激波在一個(gè)周期之中的移動(dòng)范圍和平均位置，但是卻不能反映其隨時(shí)間的變化歷程。瞬時(shí)激波位置Xsh的變化曲線與升力系數(shù)波動(dòng)曲線對(duì)比如圖6。從圖中可以看到，使用上述方法計(jì)算出的激波位置運(yùn)動(dòng)的范圍與通過(guò)平均壓力系數(shù)和壓力脈動(dòng)因子所推測(cè)的范圍基本一致，能夠大致的描繪出激波在抖振中移動(dòng)規(guī)律。通過(guò)FFT變換計(jì)算得到激波位置移動(dòng)的頻率為76 Hz，與升力系數(shù)的波動(dòng)頻率一致，其相位角比升力系數(shù)波動(dòng)的相位角略微靠前。這表明激波的移動(dòng)與翼型升力系數(shù)的波動(dòng)有很強(qiáng)的關(guān)聯(lián)，但是這兩者之間有大約30°的相位差。

圖6 激波位置與升力系數(shù)在相同時(shí)間內(nèi)的波動(dòng)Fig.6 Shock location Xshand lift coefficient evolution in the same period

圖7 一個(gè)激波震蕩周期之內(nèi)不同時(shí)刻的馬赫數(shù)云圖Fig.7 Contours of Mach number at different instant time during one cycle of flow oscillation

圖8 一個(gè)激波震蕩周期之內(nèi)不同時(shí)刻的尾緣附近的流線Fig.8 Streamlines near the tailing edge at different instant time during one cycle of flow oscillation

圖9 馬赫數(shù)云圖對(duì)應(yīng)的時(shí)刻在一個(gè)升力波動(dòng)周期中的位置(t*=t/T，T為波動(dòng)周期)Fig.9 Instant times at which the Mach number contours are intercepted in one lift fluctuation cycle (t*=t/T，T is the fluctuation period)

分別截取一個(gè)周期之內(nèi)不同時(shí)刻的馬赫數(shù)云圖如圖7，對(duì)應(yīng)的翼型尾緣附近流線如圖8，這幾個(gè)時(shí)刻在一個(gè)升力系數(shù)波動(dòng)周期之中的位置如圖9所示。其中(a)對(duì)應(yīng)于激波位置在最下游的時(shí)刻，(b)對(duì)應(yīng)于升力系數(shù)波動(dòng)曲線的波峰時(shí)刻，(c)為從激波最下游位置到最上游位置之間的一個(gè)時(shí)刻，(d)對(duì)應(yīng)于激波位置處于最上游的時(shí)刻，(e)對(duì)應(yīng)于升力系數(shù)曲線波動(dòng)的波谷時(shí)刻，(f)為從波激波最上游位置到最下游位置之間的一個(gè)時(shí)刻。在(a)時(shí)刻，激波位于一個(gè)抖振周期之中最下游的位置，翼型上表面激波根部有一個(gè)小的分離氣泡，這使附面層增厚，進(jìn)而導(dǎo)致翼型的有效彎度減小，激波按著從(a)-(b)-(c)-(d)的順序大幅向上游移動(dòng)，同時(shí)分離區(qū)逐漸增大，直至到達(dá)激波在抖振周期中的最上游位置(d)，在此期間激波的強(qiáng)度也相應(yīng)的增大;在位置(d)處分離區(qū)的范圍最大，而從此時(shí)開始激波開始按照(d)-(e)-(f)-(a)的順序沿機(jī)翼表面向下游移動(dòng)，同時(shí)激波的強(qiáng)度減弱，其中按照(d)-(e)-(f)的順序分離區(qū)逐漸減小直至機(jī)翼上表面的氣流分離區(qū)完全消失，可以看到在(e)時(shí)刻，翼型上表面的激波根部與翼型后緣處同時(shí)存在兩個(gè)分離氣泡，這說(shuō)明在分離區(qū)消失的過(guò)程中，首先是整個(gè)分離區(qū)的范圍減小直到分成了激波根部與后緣處兩個(gè)分離區(qū)，然后這兩個(gè)分離區(qū)逐漸減小，上表面接近后緣的分離區(qū)完全消失，而從(f)-(a)的過(guò)程中，準(zhǔn)確的說(shuō)是在(a)之前的某一時(shí)刻則機(jī)翼上表面產(chǎn)生出了一個(gè)小的分離氣泡;另外，升力系數(shù)波峰與波谷所處的時(shí)刻(b)和(e)出現(xiàn)的時(shí)刻都略微滯后于激波在最下游與最上游的時(shí)刻(a)和(d)。

2 抖振流動(dòng)POD分析

2.1 POD簡(jiǎn)介及快照方法

以提取靜壓的POD模態(tài)為例，按照相同的時(shí)間間隔tm得到M份快照:

式中:m=1，2，…，M。

壓力場(chǎng)可以分解為平均流動(dòng)分量和瞬時(shí)波動(dòng)分量:

使用POD方法就可以將p'(x，t)化為如下形式:

式中:pi(x)是i階POD基函數(shù)，時(shí)間變化特性通過(guò)第i階模態(tài)的時(shí)間系數(shù)ai(t)表現(xiàn)出來(lái)。

定義M×M階自相關(guān)矩陣C，其中每一個(gè)元素Cmn如下:

式中:m，n=1，2，…，M，下標(biāo)Ω表示在整個(gè)的空間域進(jìn)行積分。

為了獲得出第i個(gè)模態(tài)的時(shí)間系數(shù)向量a(i)，需要求解如下的特征值問(wèn)題

式中:特征值 λi代表了流場(chǎng)投影在第 i階基函數(shù)pi(x)上的平均“能量”。相關(guān)系數(shù)矩陣的特征向量就是時(shí)間系數(shù)向量。

按照正實(shí)特征值由大到小的順序?yàn)槟B(tài)進(jìn)行排序，即λ1≥λ2≥λ3≥… ＞0;特征值為零的模態(tài)對(duì)速度場(chǎng)無(wú)影響，所以將其去掉。由于自相關(guān)矩陣C是自伴隨矩陣并且是半正定的，所以時(shí)間系數(shù)相互正交?？梢缘玫絇OD模態(tài)如下:

使用快照方法在時(shí)域中進(jìn)行POD分解，則其產(chǎn)生的相關(guān)系數(shù)矩陣的規(guī)模為M×M;而如果在空間域進(jìn)行POD分解，則其關(guān)系數(shù)矩陣的規(guī)模為Ng× Ng，Ng為網(wǎng)格點(diǎn)的數(shù)目?？梢?jiàn)，快照方法更適用于處理大量的CFD計(jì)算數(shù)據(jù)。

2.2 OAT15A翼型抖振流場(chǎng)的POD分析

為了提取抖振流場(chǎng)靜壓分布的POD模態(tài)，使用C網(wǎng)格的280個(gè)瞬時(shí)流場(chǎng)結(jié)果作為“快照”，確保其包含的時(shí)間長(zhǎng)度大于一個(gè)周期，并將計(jì)算出的壓強(qiáng)除以的遠(yuǎn)場(chǎng)壓強(qiáng)得到無(wú)量綱壓強(qiáng)。提取的前六階POD模態(tài)每個(gè)模態(tài)的特征值占全部特征值的百分比以及前n階特征值之和占全部特征值的百分比如表2，前50階模態(tài)特征值的衰減如圖10所示?？梢钥吹降谝浑A與第二階模態(tài)的特征值占全部模態(tài)特征值的比重最大，這兩者的大小相近，其它階模態(tài)的特征值相對(duì)較小，且從第三階模態(tài)開始特征值衰減的速率就開始減緩。如前所述，模態(tài)特征值的大小表征了此階模態(tài)“能量”所占系統(tǒng)總“能量”的大小，所以可以認(rèn)為第一階模態(tài)與第二階模態(tài)在抖振現(xiàn)象中占主導(dǎo)地位。

圖10 前50階模態(tài)特征值的變化曲線Fig.10 The eigenvalue attenuation of the first fifty modes

表2 壓力前六階POD模態(tài)每個(gè)模態(tài)特征值所占比重以及前n階模態(tài)特征值之和所占比重Table 2 Eigenvalue of the six first POD modes and the sum of the n first POD modes expressed as a percentage of the sum of all eigenvalues

計(jì)算得到的流場(chǎng)靜壓的第一階、第二階模態(tài)如圖11。可以看到第一階模態(tài)包含了激波運(yùn)動(dòng)與由其產(chǎn)生的氣流分離，其中尤以激波根部與附面層之間相互作用的區(qū)域最為顯著;第二階模態(tài)則主要與氣流分離有關(guān)。

根據(jù)公式(5)可知，流場(chǎng)中壓力的波動(dòng)值為計(jì)算所得的每一階的模態(tài)值乘上同一時(shí)刻的時(shí)間系數(shù)再求和，時(shí)間系數(shù)為無(wú)量綱值。計(jì)算得到的前六階模態(tài)的時(shí)間系數(shù)隨時(shí)間的變化如圖12，可以看到第一階模態(tài)與第二階模態(tài)的時(shí)間系數(shù)的波動(dòng)幅度最大，與其它模態(tài)的時(shí)間系數(shù)相比占主導(dǎo)地位;且隨著模態(tài)階數(shù)的提高，模態(tài)的時(shí)間系數(shù)的幅值越來(lái)越小、頻率越來(lái)越高并更加不規(guī)則;第一階模態(tài)與第二階模態(tài)的時(shí)間系數(shù)頻率均為76Hz，與升力系數(shù)波動(dòng)的頻率一致，可以認(rèn)為抖振現(xiàn)象主要為第一階模態(tài)和第二階模態(tài)合成的結(jié)果。另外，第一階模態(tài)的時(shí)間系數(shù)波動(dòng)的波峰值與第二階模態(tài)的波峰值之間大概相差四分之一周期，也就是說(shuō)第一階模態(tài)滯后于第二階模態(tài)90°的相位，這是由POD模態(tài)的正交性所決定的。

圖11 第一階與第二階POD模態(tài)云圖Fig.11 Contours of the first and the second POD modes

圖12 前六階模態(tài)的時(shí)間系數(shù)隨時(shí)間的變化Fig.12 Evolution of temporal coefficients of the first six mode

將第一階、第二階模態(tài)的時(shí)間系數(shù)與相同時(shí)間內(nèi)升力系數(shù)的波動(dòng)在同一副圖中比較，如圖13，其中第一階模態(tài)的時(shí)間系數(shù)的相位角要比升力系數(shù)的相位角超前大約30°，與激波位置的相位角十分接近，說(shuō)明第一階模態(tài)主要表征了激波運(yùn)動(dòng)對(duì)壓力分布的影響。由于這些模態(tài)都對(duì)應(yīng)于實(shí)際的空間結(jié)構(gòu)，那么這也從另一個(gè)方面表明第二階模態(tài)主要表征了翼型上表面的氣流分離。

圖13 前兩階模態(tài)的時(shí)間系數(shù)與相同時(shí)間內(nèi)的升力系數(shù)Fig.13 Temporal coefficients of the first two modes and lift coefficient in the same period

本文采用POD模態(tài)分析方法相對(duì)于直接對(duì)非定常結(jié)果進(jìn)行分析而言，可以更加明確的辨別主導(dǎo)跨聲速抖振現(xiàn)象的主要流動(dòng)相干結(jié)構(gòu)，但是仍然需要進(jìn)行非定常流場(chǎng)的計(jì)算。

3 結(jié)論

(1)在抖振過(guò)程中激波運(yùn)動(dòng)的相位并不與升力系數(shù)的波動(dòng)相位完全重合，這兩者之間有一定的相位差，要明確他們之間的關(guān)系需要更進(jìn)一步的研究。

(2)在根據(jù)OAT15A的非定常計(jì)算結(jié)果所提取的POD模態(tài)中，第一階模態(tài)與第二階模態(tài)的頻率與抖振頻率相同，是主導(dǎo)抖振現(xiàn)象的主要模態(tài)。

［1］ Mou R K，Yang Y N.Advances of studies for the buffet problem of aircraft［J］.Chinese Journal of Applied Mechanics，2001，18: 142-150.(in Chinese)

牟讓科，楊永年.飛機(jī)抖振問(wèn)題研究進(jìn)展［J］.應(yīng)用力學(xué)學(xué)報(bào)，2001，18(SI):142-150.

［2］ Jacquin L，Molton P，Deck S，et al.An experimental study of shock oscillation over a transonic supercritical profile［C］//Toronto: American Institute of Aeronautics and Astronautics，Inc.，2005.

［3］ Jacquin L，Molton P，Deck S，et al.Experimental study of shock oscillation over a transonic supercritical profile［J］.AIAA Journal，2009，47(9):1985-1994.

［4］ Sebastien Deck.Numerical simulation of transonic buffet over a supercritical airfoil［J］.AIAA J.，2005，l(43):1156-1566.

［5］ Mylene Thiery，Eric Coustols.Numerical prediction of shock induced oscillations over a 2D airfoil:Influence of turbulence modelling and test section walls［J］.International Journal of Heat and Fluid Flow，2006，27:661-670.

［6］ Sebastian Illi，Thorsten Lutz，Ewald Kramer.On the capability of unsteady RANS to predict transonic buffet［C］//Raunschweig，Germany:Institute of Aerody-namics and Gas Dynamics，2012.

［7］ Crouch J D，Garbaruk A，Magidov，et al.Global structure of buffeting flow on transonic airfoils［C］//In IUTAM Symposium on Unsteady Separated Flows and Their Control，Springer:2009.

［8］ Crouch J D，Garbaruk A，Magidov，et al.Origin of transonic buffet on aerofoils［J］.J.Fluid Mech.，2009，628:357-369.

［9］ Xiong Juntao，Li Feng，Luo Shijun.Computation of NACA0012 airfoil transonic buffet phenomenon with unsteady Navier-Stokes equations［C］//Nashville，Tennessee:University of California Irvine，2012.

［10］Xiong Juntao，Liu Feng.Numerical simulation of transonic buffet on swept wing of supercritical airfoils［C］//San Diego:University of California Irvine，2013.

［11］Molton P，Bur R，Lepage A，et al.Control of buffet phenomenon on a transonic swept wing［C］//Applied Aerodynamics Symposium，2012.

［12］Zhang W，Chen C，Sun D J.Numerical simulation of flow around two side-by-side circular cylinders at low Reynolds numbers by a POD-Galerkin spectral method［J］.J.of Hydrodinamics，2009，24 (1)，83-88.(in Chinese)

張偉，陳誠(chéng)，孫德軍.低Reynolds數(shù)橫向排列雙圓柱繞流的POD—Galerkin譜方法數(shù)值模擬［J］.水動(dòng)力學(xué)研究與進(jìn)展，2009，24(1):83-88.

［13］Wang A X，Ma Y C，F(xiàn)u Y.Proper orthogonal decomposition for the nonstationary Navier-Stokes equations based on two-grid method［J］.Basic Sciences Journal of Textile Universtiles，2009，22(1): 76-81.(in Chinese)

王阿霞，馬逸塵，付英.基于雙重網(wǎng)格法的非定常N-S方程POD數(shù)值模擬［J］.紡織高校基礎(chǔ)科學(xué)學(xué)報(bào)，2009，22(1):76-81.

［14］Couplet M，Basdevant C，Sagaut P.Calibrated reduced-order PODGalerkin system for fluid flow modelling［J］.Journal of Computational Physics，2005，207(1):192-220.

［15］Bui-Thanh T，Damodaran M，Willcox K.Aerodynamic data reconstruction and inverse design using proper orthogonal decomposition［J］.AIAA Journal，2004，42(8):1501-1516.

［16］Lumley J L.The structure of inhomogeneous turbulence［C］//Moscow:Nauka，1967.

［17］Aubry N，Holmes P，Lumley J L，et al.The dynamics of coherent structures in the wall region of a turbulent boundary layer［J］.J.Fluid Mech.，1988，192(1):115-173.

［18］Moehlis J，Smith T R，Holmes P，et al.Models for turbulent plane Couette flow using the proper orthogonal decomposition［J］.Phys.Fluids，2002，14(7):2493-2507.

［19］Roe P L.Characteristic based schemes for the euler equations［J］.Annual Review of Fluid Mechanics，1986，l(18):337-365.

［20］Jameson A.Time dependent calculations using multi-grid，with applications to unsteady flows past airfoils and wings［R］.AIAA-91-1596，1991.

［21］Golias N A，Tsiboukis T D.An approach to refin-ing three dimensional tetrahedral meshes based on delaunay transformations［J］.Int.J.Num.Meth.Eng.，1994，l(37):793-812.

［22］Spalart P R，Allmaras S R.A one-equation turbulence model for aerodynamic flows［R］.AIAA-92-0439，1992.

CFD computation and POD analysis for transonic buffet on a supercritical airfoil

Dong Shenghua，Shi Aiming*，Ye Zhengyin，Tian Haitao
(College of Aeronautics，Northwestern Polytechnical University，Xi’an 710072，China)

The transonic buffet flow on the OAT15A supercritical airfoil is simulated by the solution of the cell-centered finite-volume method(FVM)based RANS(Reynolds-averaged Navier-Stokes)equations with the S-A(Spalart-Allmaras)one-equation turbulence model，as well as the implicit time-stepping scheme and the Roe scheme.On the basis of successful simulation of transonic buffet phenomenon caused by the shock-boundary layer interaction，studying the evolution of flow separation and it reveals a hysteresis phenomenon that the phase of shock movement is slightly ahead of the phase of lift fluctuation.To view the physical nature of transonic buffet，apply POD(Proper Orthogonal Composition)method to extract POD modes with the transient results as snapshots.An POD analysis from the perspective of coherent structure has been implemented.

buffet;shock-boundary layer interaction;proper orthogonal composition

V211.3

A

10.7638/kqdlxxb-2013.0100

0258-1825(2015)04-0481-07

2013-11-15;

2013-12-02

國(guó)家自然科學(xué)基金(10602046)

董圣華(1988-)，男，朝鮮族，黑龍江牡丹江人，碩士研究生，主要研究領(lǐng)域:跨聲速抖振.E-mail:1024424029@qq.com

史愛(ài)民*(1977-)，男，江蘇金壇人，副教授/博士，主要從事空氣動(dòng)力學(xué)和流固耦合力學(xué)研究.E-mail:sam@nwpu.edu.cn

董圣華，史愛(ài)明，葉正寅，等.超臨界翼型跨聲速抖振CFD計(jì)算和POD分析［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2015，33(4):481-487.

10.7638/kqdlxxb-2013.0100 Dong S H，Shi A M，Ye Z Y，et al.CFD computation and POD analysis for transonic buffet on a supercritical airfoil［J］.Acta Aerodynamica Sinica，2015，33(4):481-487.