許常悅，王從磊，孫建紅

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

0 引 言

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

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

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

1 數(shù)值計算方法

1.1 控制方程

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

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

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

方程（2）和（3）中的SGS未封閉項的定義如下：

1.2 動力SGS模型

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

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

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

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

SGS能量項H的?；缦拢?/p>

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

SGS湍流擴散項的模化借鑒Knight等人［15］的做法，即J＝。

1.3 離散格式和邊界條件

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

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

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

2.1 計算參數(shù)

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

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

2.2 計算結(jié)果驗證

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

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

圖2 沿柱體表面的平均壓力〈pw〉（a）、表面摩擦系數(shù)〈｜Cf｜〉（b）和壓力脈動均方根值（c）分布。這里，θ＝0對應(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 大尺度激波運動

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

圖3 利用密度梯度?！尅硎镜募げㄟ\動，時間間隔為0.9D／U∞Fig.3 Shock wave motion depicted by the magnitude of density gradient‖▽‖.The time interval is 0.9D／U∞

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

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

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

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

2.4 非定常湍流剪切層行為

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

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

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

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

圖7 P2點處的剪切應(yīng)力象限分解（a）和瞬時剪切應(yīng)力角概率密度函數(shù)分布（b），P2點的位置見圖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ù)也為兩個比例關(guān)系，即－｛u″w″｝／k和Ruw＝－｛u″w″｝／（σuσw）。這里，k＝（σ2u＋σ＋σ）／2為湍動能，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é) 論

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