圓柱跨聲速繞流中的激波／湍流相互作用大渦模擬研究

許常悅，王從磊，孫建紅

（南京航空航天大學(xué) 航空宇航學(xué)院，江蘇 南京210016）

0 引 言

激波／湍流相互作用問題是可壓縮湍流研究領(lǐng)域中的挑戰(zhàn)性問題，與其相關(guān)的機(jī)理認(rèn)識(shí)已成為工業(yè)應(yīng)用上必須解決的問題，并得到了廣泛的研究［1－4］。該類問題的研究中常選用一些規(guī)范構(gòu)型作為研究對(duì)象，如平板、斜坡、鼓包、翼型和柱體，等等。研究的復(fù)雜性與構(gòu)型有關(guān)［3］，彎曲壁面會(huì)產(chǎn)生復(fù)雜的流場(chǎng)結(jié)構(gòu)。對(duì)于多數(shù)實(shí)際流動(dòng)而言，激波／湍流相互作用常發(fā)生在彎曲壁面上，這與該處的湍流邊界層中存在較大的壓力梯度有關(guān)［4］。研究彎曲壁面物體繞流中的激波／湍流相互作用問題受到了廣泛的關(guān)注，尤其是圓柱跨聲速繞流。

目前已有一些關(guān)于圓柱跨聲速繞流的實(shí)驗(yàn)和數(shù)值研究。Macha［5］、Murthy＆Rose［6］和 Rodriguez［7］實(shí)驗(yàn)研究了雷諾（Re）數(shù)約為105的圓柱跨聲速繞流問題，展示了此類流動(dòng)中的阻力上升現(xiàn)象和激波結(jié)構(gòu)。近年來(lái)，Riserda＆Leal［8］采用有限體積方法求解二維可壓縮Navier－Stokes方程，研究了來(lái)流馬赫數(shù)為0.8的圓柱繞流問題，主要分析了柱體的非定常受力和近場(chǎng)流場(chǎng)結(jié)構(gòu)。在之前的工作中［9－10］，我們已經(jīng)對(duì)來(lái)流馬赫數(shù)對(duì)圓柱跨聲速繞流的影響進(jìn)行了研究，并分析了不同流動(dòng)狀態(tài)下的局部超聲速區(qū)形成機(jī)制和剪切層的不穩(wěn)定性過(guò)程，然而關(guān)于大尺度的激波運(yùn)動(dòng)和非定常湍流剪切層的湍流特性等問題并未詳細(xì)研究，這正是本文的研究目的。

在當(dāng)前的計(jì)算流體力學(xué)中，雷諾平均Navier－Stokes方程（RANS）方法、直接數(shù)值模擬（DNS）方法和大渦模擬（LES）方法是主要的湍流研究手段。RANS方法在計(jì)算激波／湍流相互作用問題時(shí)無(wú)法給出可靠的結(jié)果，這是由于RANS方法不能正確捕捉流動(dòng)分離和再附過(guò)程，并且無(wú)法預(yù)測(cè)低頻的激波運(yùn)動(dòng)。DNS方法能夠再現(xiàn)湍流中各種尺度的渦結(jié)構(gòu)。在高Re數(shù)流動(dòng)中，湍流尺度有很寬的譜域，只有網(wǎng)格數(shù)達(dá)到Re3（包括時(shí)間維）的量級(jí)，才能有足夠的時(shí)空分辨能力。因此，目前的計(jì)算條件下DNS方法還只能研究Re數(shù)較低的湍流問題［11］。LES方法直接計(jì)算湍流中較大尺度的渦，而將動(dòng)力學(xué)意義不大的小渦用亞格子（SGS）模型來(lái)模擬，因此它可以反映較多的湍流信息，具有較好的普適性。近年來(lái)，LES方法已經(jīng)被成功應(yīng)用于激波／湍流相互作用問題的研究［3－4］，故本文采用 LES方法作為研究手段。

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

1.1 控制方程

當(dāng)前采用的控制方程為守恒型的三維Favre濾波可壓縮Navier－Stokes方程。方程的無(wú)量綱化采用來(lái)流密度ρ∞、溫度T∞、速度U∞和圓柱直徑D作為特征量。無(wú)量綱的控制方程可以寫成如下形式：

其中，頂標(biāo)“－”和“～”分別表示空間濾波和Favre濾波，即＝／。變量ρ、ui、p和E分別為密度、速度分量、壓力和總能。擴(kuò)散項(xiàng)通過(guò)下面式子給出：

這里，μ為分子粘性系數(shù)，Cp為定壓比熱比，Pr為Prandtl數(shù)，Sij＝（?ui／?xj＋?uj／?xi）／2為應(yīng)變率張量。此外，理想氣體方程和計(jì)算μ的Sutherland公式被采用。

方程（2）和（3）中的SGS未封閉項(xiàng)的定義如下：

1.2 動(dòng)力SGS模型

控制方程經(jīng)過(guò)濾波后會(huì)產(chǎn)生一些SGS未封閉項(xiàng)，需要進(jìn)行模化。由于對(duì)能量方程的影響較小，該項(xiàng)可以忽略；粘性應(yīng)力非線性項(xiàng)和熱通量非線性項(xiàng)Q經(jīng)常被忽略。下面對(duì)、H和J未封閉項(xiàng)的?；^(guò)程進(jìn)行詳細(xì)介紹。

Yoshizawa［12］針對(duì)弱可壓縮湍流，提出了一種多尺度直接相互作用的近似模型：SGS應(yīng)力項(xiàng)τ反對(duì)稱部分τ－δijτ／3采用Smagorinsky模型進(jìn)行?；?，而對(duì)稱部分τ則單獨(dú)?；?/p>

其中，｜｜＝＝（ΔxΔyΔz）1／3為濾波寬度。CR和CI為模型系數(shù)，需要通過(guò)二次濾波動(dòng)態(tài)求解。Lilly［13］建議采用最小二乘法獲得模型系數(shù)：

其中，〈·〉表示局部光滑操作，即采用當(dāng)?shù)?7個(gè)網(wǎng)格單元的體積加權(quán)平均實(shí)現(xiàn)，該方法常被用來(lái)提高數(shù)值計(jì)算的穩(wěn)定性［14］。Lij、Mij和β的形式分別為：

SGS能量項(xiàng)H的?；缦拢?/p>

其中，＝2｜｜。同模型系數(shù)CR和CI，湍流Prandtl數(shù)的求法如下：

SGS湍流擴(kuò)散項(xiàng)的?；梃bKnight等人［15］的做法，即J＝。

1.3 離散格式和邊界條件

本文采用有限體積方法求解Favre濾波后的控制方程（1）－（3）。對(duì)流項(xiàng)的離散采用二階的中心－迎風(fēng)型混合格式。在該混合格式中，中心格式和Roe通量差分裂迎風(fēng)格式之間的切換，通過(guò)一個(gè)二進(jìn)制開關(guān)函數(shù)Φi＋1／2來(lái)實(shí)現(xiàn)：

程序中以自由來(lái)流為初始條件，遠(yuǎn)場(chǎng)條件采用基于當(dāng)?shù)匾痪SRiemann不變量的無(wú)反射特征邊界，壁面采用無(wú)滑移、無(wú)穿透的絕熱邊界條件，展向采用周期性邊界條件。

2 計(jì)算結(jié)果與分析

2.1 計(jì)算參數(shù)

基于已有的實(shí)驗(yàn)研究［5－7］，本文選取的來(lái)流馬赫數(shù)M∞為0.75，基于圓柱直徑D的雷諾數(shù)為2×105。采用O型計(jì)算網(wǎng)格，壁面和激波區(qū)附近的局部網(wǎng)格如圖1所示，壁面法向最小網(wǎng)格尺度為10－5D。經(jīng)過(guò)仔細(xì)測(cè)試［9－10］，流向、周向和展向的網(wǎng)格數(shù)分別取為513、513和121，徑向計(jì)算域直徑為50D，展向長(zhǎng)度為4D，時(shí)間步長(zhǎng)為0.00375D／U∞。為了表述方便，對(duì)本文的符號(hào)作如下說(shuō)明：〈〉表示時(shí)間和展向同時(shí)平均，｛｝表示密度加權(quán)平均，即｛φ｝＝φ〉／〉；脈動(dòng)密度和脈動(dòng)壓力為＝－〈〉和＝－〉，脈動(dòng)速度為＝－｝。

圖1 壁面和激波區(qū)附近的局部計(jì)算網(wǎng)格Fig.1 Local computational grid near the wall and shock region

2.2 計(jì)算結(jié)果驗(yàn)證

為了驗(yàn)證當(dāng)前計(jì)算結(jié)果的可信性，本文和已有實(shí)驗(yàn)數(shù)據(jù)［5－7］進(jìn)行比較，對(duì)比的物理量有壁 面壓力〈pw〉、表面摩擦系數(shù)〈｜Cf｜〉、壓力脈動(dòng)均方根值p、時(shí)均阻力系數(shù)〈CD〉t、升力系數(shù)均方根值CLrms和渦脫瀉Strouhal（St）數(shù)。圖2給出了沿柱體表面的〈pw〉、〈｜Cf｜〉和變化曲線。表1給出了LES計(jì)算結(jié)果和實(shí)驗(yàn)數(shù)據(jù)的對(duì)比情況。從圖2和表1中可以看出，當(dāng)前LES計(jì)算結(jié)果和已有的實(shí)驗(yàn)數(shù)據(jù)相符較好，這表明當(dāng)前的LES計(jì)算結(jié)果具有較好的可信性。

表1 當(dāng)前計(jì)算結(jié)果和實(shí)驗(yàn)數(shù)據(jù)的對(duì)比Table 1 Comparison of the present calculated results with the experimental data

圖2 沿柱體表面的平均壓力〈pw〉（a）、表面摩擦系數(shù)〈｜Cf｜〉（b）和壓力脈動(dòng)均方根值（c）分布。這里，θ＝0對(duì)應(yīng)柱體的前緣Fig.2 Distributions of the mean pressure〈pw〉（a），skin friction coefficient〈｜Cf｜〉（b）and root－mean－square value of pressure fluctuationon the cylinder surface（c）.Here，θ＝0increases from the front－point of the cylinder in the clockwise direction

2.3 大尺度激波運(yùn)動(dòng)

已有研究結(jié)果［7，9－10］表明，當(dāng)M∞＝0.75時(shí)圓柱分離點(diǎn)處會(huì)形成一道斜激波。關(guān)于翼型表面的激波運(yùn)動(dòng)已得到廣泛的研究［17］，結(jié)果表明存在三種運(yùn)動(dòng)模式：正弦類運(yùn)動(dòng)（A類）、間歇式運(yùn)動(dòng)（B類）和向上游傳播類運(yùn)動(dòng)（C類）。然而，關(guān)于圓柱表面的激波運(yùn)動(dòng)模式研究較少。為了定性地分析圓柱表面的激波運(yùn)動(dòng)模式，圖3給出了不同時(shí)刻的瞬時(shí)流場(chǎng)圖。可以看出，圓柱分離點(diǎn)處的激波在柱體上下表面交替出現(xiàn)，并且向上游傳播，這與翼型表面的C類激波運(yùn)動(dòng)模式相似。為了獲得該激波運(yùn)動(dòng)的特征St數(shù)，在其運(yùn)動(dòng)區(qū)內(nèi)的設(shè)置一個(gè)探測(cè)點(diǎn)P1，如圖3（a）所示。圖4給出了P1點(diǎn)處的能譜分布，可以看出激波運(yùn)動(dòng)的特征St數(shù)與渦脫瀉St數(shù)一致，均為0.19。此外，－5／3次方斜率關(guān)系表明當(dāng)前計(jì)算結(jié)果達(dá)到了湍流的慣性子區(qū)［18］。

圖3 利用密度梯度模‖▽‖表示的激波運(yùn)動(dòng)，時(shí)間間隔為0.9D／U∞Fig.3 Shock wave motion depicted by the magnitude of density gradient‖▽‖.The time interval is 0.9D／U∞

圖4 P1點(diǎn)處的能譜曲線。P1點(diǎn)的位置見圖3（a）Fig.4 Resolved energy spectrum at probe P1in Fig.3 （a）

圖5 基于壓力場(chǎng)POD分解的前四個(gè)模態(tài)空間分布Fig.5 Spatial distributions of first four POD modes based on the pressure fields

為了分析激波運(yùn)動(dòng)對(duì)流動(dòng)模態(tài)的影響，采用本征正交分解（POD）方法研究流場(chǎng)的主模態(tài)。對(duì)于一個(gè)給定的流動(dòng)變量f（x，t），POD分析可以確定一族正交函數(shù)φj（x），j＝1，2，…，則f在前n個(gè)函數(shù)上的投影為最

小誤差定義為〈‖f－‖2〉，這里分別表示時(shí)間平均和L2空間模，aj（t）代表第j個(gè)模態(tài)隨時(shí)間變化的系數(shù)［19］。圖5給出了基于壓力場(chǎng)POD分解的前四個(gè)模態(tài)，實(shí)線表示正值，虛線表示負(fù)值，30個(gè)等值線從－1.5ρ∞增加到1.5ρ∞U。可以看出，mode1和mode3呈現(xiàn)反對(duì)稱模態(tài)，這與柱體分離點(diǎn)處的大尺度激波運(yùn)動(dòng)有關(guān)；mode2和mode4呈現(xiàn)對(duì)稱模態(tài)，這與柱體尾跡中的渦脫瀉現(xiàn)象有關(guān)，相似的結(jié)論也見于翼型跨聲速繞流問題的研究［20］。

2.4 非定常湍流剪切層行為

已有的研究［9－10］主要關(guān)注剪切層的不穩(wěn)定性過(guò)程，本文將對(duì)剪切層中的壓力信號(hào)特性和剪切層的湍流特性進(jìn)行研究。

首先考察的是剪切層中的壓力信號(hào)特性。為此，在剪切層附近設(shè)置一個(gè)壓力探測(cè)點(diǎn)P2，如圖3（a）所示。P2點(diǎn)處的壓力信號(hào)功率譜（PSD）分析如圖6所示，可以看出該處的特征St數(shù)與渦脫瀉St數(shù)一致，并且PSD曲線中存在一些特征斜率關(guān)系。由于圓柱剪切層的擾動(dòng)源來(lái)自柱體分離點(diǎn)，故剪切層附近的壓力信號(hào)PSD曲線中呈現(xiàn)出的斜率關(guān)系與壁面附近的壓力信號(hào)一致：低頻區(qū)的0.4次方斜率關(guān)系與邊界層外層區(qū)域的湍流有關(guān)［21］，在低頻區(qū)向高頻區(qū)過(guò)渡區(qū)，邊界層對(duì)數(shù)律區(qū)的渦結(jié)構(gòu)會(huì)導(dǎo)致－1次方的斜率關(guān)系［18］，高頻區(qū)則呈現(xiàn)－5次方的斜率關(guān)系［22］。

圖6 P2點(diǎn)處的壓力信號(hào)功率譜（PSD）曲線，P2點(diǎn)的位置見圖3（a）Fig.6 Power spectral density（PSD）of pressure at P2in Fig.3 （a）

其次，采用象限分解技術(shù)［23］研究非定常剪切層的湍流特性。流場(chǎng)中的瞬時(shí)流向脈動(dòng)速度u″和橫向脈動(dòng)速度w″通過(guò)采樣流場(chǎng)計(jì)算得到，并且u″和w″分別采用流向脈動(dòng)速度均方根值σu和橫向脈動(dòng)速度均方根值σw正則化。圖7（a）給出了剪切層區(qū)域內(nèi)P2點(diǎn)處的平均剪切應(yīng)力｛u″w″｝的象限分解。這里，剪切應(yīng)力值取為0.01。從圖中可以看出，剪切層區(qū)域的脈動(dòng)速度趨于有組織性。為了分析剪切應(yīng)力的演化特性，本文計(jì)算了瞬時(shí)的剪切應(yīng)力角ψ，其定義為ψ＝tan－1（w″／u″）。ψ＜0對(duì)應(yīng)象限分解中的第二、四象限；ψ＞0對(duì)應(yīng)象限分解中的第一、三象限；ψ＝0表示僅存在流向脈動(dòng)速度u″。圖7（b）給出了P2點(diǎn)處的瞬時(shí)剪切應(yīng)力角的概率密度函數(shù)（PDF）分布。這里，N表示采樣流場(chǎng)數(shù)目，NT表示總的采樣流場(chǎng)數(shù)目?？梢钥闯?，P2點(diǎn)處的剪切應(yīng)力角ψ趨于0，這表明該處的脈動(dòng)速度以流向脈動(dòng)速度為主。

圖7 P2點(diǎn)處的剪切應(yīng)力象限分解（a）和瞬時(shí)剪切應(yīng)力角概率密度函數(shù)分布（b），P2點(diǎn)的位置見圖3（a）Fig.7 Shear stress quadrant decomposition（a）and PDF of the shear－stress angle（b）at probe P2in Fig.3 （a）

最后，為了深入分析剪切層中的湍流特性，考察了正應(yīng)力各向異性比例和湍流結(jié)構(gòu)參數(shù)沿剪切層的變化，如圖8所示。這里，分析的正應(yīng)力各向異性比例有流向－展向比（σu／σw）2和橫向 －展向比（σv／σw）2。從圖8（a）中可以看出，流向－展向比的值沿剪切層從45迅速衰減，而橫向－展向比的值幾乎為0，這表明剪切層中的流向正應(yīng)力占主導(dǎo)地位，這與前文象限分解得出的結(jié)論一致。分析的湍流結(jié)構(gòu)參數(shù)也為兩個(gè)比例關(guān)系，即－｛u″w″｝／k和Ruw＝－｛u″w″｝／（σuσw）。這里，k＝（σ2u＋σ＋σ）／2為湍動(dòng)能，Ruw為剪切應(yīng)力相關(guān)系數(shù)。圖8（b）給出了－｛u″w″｝／k和Ruw沿剪切層的變化，可以看出－｛u″w″｝／k沿剪切層增加，而Ruw沿剪切層衰減，這與剪切層中的大尺度結(jié)構(gòu)組織性減小有關(guān)［24］。

圖8 正應(yīng)力各向異性比例（a）和湍流結(jié)構(gòu)參數(shù)（b）沿剪切層變化Fig.8 Evolution of the normal shear stress anisotropy ratio（a）and of turbulence structure parameter（b）along the shear layer

3 結(jié) 論

本文采用大渦模擬方法研究了圓柱跨聲速繞流中的激波／湍流相互作用問題，分析了大尺度的激波運(yùn)動(dòng)、剪切層中的壓力信號(hào)特征和湍流特性。來(lái)流馬赫數(shù)M∞取為0.75，基于圓柱直徑D的雷諾數(shù)為2×105。通過(guò)對(duì)計(jì)算結(jié)果的分析和討論有如下結(jié)論：（1）圓柱分離點(diǎn)處出現(xiàn)一道向上游傳播的斜激波，并且運(yùn)動(dòng)的St數(shù)與渦脫瀉St數(shù)一致；（2）分離點(diǎn)處的激波運(yùn)動(dòng)導(dǎo)致流場(chǎng)中出現(xiàn)反對(duì)稱的流動(dòng)模態(tài)；（3）圓柱剪切層中的壓力信號(hào)PSD曲線中存在0.4、－1和－5次方的斜率關(guān)系，剪切層中的脈動(dòng)速度以流向脈動(dòng)速度為主，并且沿剪切層的大尺度結(jié)構(gòu)組織性減小。

［1］DECK S.Numerical simulation of transonic buffet over a supercritical airfoil［J］.AIAAJ.，2005，43（7）：1556－1566.

［2］PIROZZOLI S，GRASSO F.Direct numerical simulation of impinging shock wave／turbulent boundary layer interaction atM＝2.25［J］.Phys.Fluids，2006，18（6）：065113.

［3］LOGINOV M S，ADAMS N A，ZHELTOVODOV AA. Large－eddy simulation of shock－wave／turbulentboundary－layer interaction［J］.J.FluidMech.，2006，565：135－169.

［4］SANDHAM H D，YAO Y F，LAWAL A A.Large－eddy simulation of transonic turbulent flow over a bump［J］.Int.J.HeatFluidFlow，2003，24（4）：584－595.

［5］MACHA J M.Drag of circular cylinders at transonic Mach numbers［J］.J.Aircraft，1977，14（6）：605－607.

［6］MURTHY V S，ROSE W C.Detailed measurements on a circular cylinder in cross flow［J］.AIAAJ.，1978，16（6）：549－550.

［7］RODRIGUEZ O.The circular cylinder in subsonic and transonic flow［J］.AIAAJ.，1984，22（12）：1713－1718.

［8］MISERDA R F B，LEAL R G.Numerical simulation the unsteady aerodynamic forces over a circular cylinder in transonic flow［R］.AIAA Paper 2006－1408，2006.

［9］XU C Y，CHEN L W，LU X Y.Effect of Mach number on transonic flow past a circular cylinder［J］.ChineseSci.Bull.，2009，54（11）：1886－1893.

［10］XU C Y，CHEN L W，LU X Y.Numerical investigation of shock wave and turbulence interaction over a circular cylinder［J］.Mode.Phys.Lett.B，2009，23（3）：233－236.

［11］MOIN P，MAHESH K.Direct numerical simulation：a tool in turbulence research［J］.Annu.Rev.Fluid Mech.，1998，30：539－578.

［12］YOSHIZAWA A.Statistical theory for compressible turbulent shear flows，with the application to subgrid modeling［J］.Phys.Fluids，1986，29（7）：2152－2164.

［13］LILLY D K.A proposed modification of the Germano subgrid scale closure method［J］.Phys.Fluids，1992，4（3）：633－635.

［14］MOIN P，SQUIRES K，CABOT W，et al.A dynamic subgrid－scale model for compressible turbulence and scalar transport［J］.Phys.Fluids，1991，3（11）：2746－2757.

［15］KNIGHT D，ZHOU G，OKONG＇O N，et al.Compressible large eddy simulation using unstructured grids［R］.AIAA paper 1998－0535，1998.

［16］SWANSON R C，TURKEL E.On central－difference and upwind schemes［J］.J.Comput.Phys.，1992，101（2）：292－306.

［17］TIJDEMAN H，SEEBASS R.Transonic flow past oscillating airfoils［J］.Ann.Rev.FluidMech.，1980，12：181－222.

［18］NA Y，MOIN P.The structure of wall－pressure fluctuations in turbulent boundary layers with adverse pressure gradient and separation［J］.J.FluidMech.，1998，377：347－373.

［19］BERKOOZ G，HOLMES P，LUMLEY J L.The proper orthogonal decomposition in the analysis of turbulent flows［J］.Annu.Rev.FluidMech.，1993，25：539－575.

［20］BOURGUET R，BRAZA M.Reduced－order modeling for unsteady transonic flows around an airfoil［J］.Phys.Fluids，2007，19（11）：111701.

［21］PANTON R L，LINEBARGER J H.Wall pressure spectra calculations for equilibrium boundary layers［J］.J.FluidMech.，1974，65：261－287.

［22］BLAKE W K.Mechanics of flow－induced sound and vibration［M］.New York：Academic Press，1986.

［23］WALLACE J M，ECKELMANN H，BRODKEY R S.The wall region in turbulent shear flow［J］.J.Fluid Mech.，1972，54：39－48.

［24］HERRIN J L，DUTTON J C.The turbulence structure of a reattaching axisymmetric compressible free shear layer［J］.Phys.Fluids，1997，9（11）：3502－3512.