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

許常悅，王從磊，孫建紅

（南京航空航天大學 航空宇航學院，江蘇 南京210016）

0 引 言

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

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

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

1 數值計算方法

1.1 控制方程

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

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

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

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

1.2 動力SGS模型

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

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

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

其中，〈·〉表示局部光滑操作，即采用當地27個網格單元的體積加權平均實現，該方法常被用來提高數值計算的穩定性［14］。Lij、Mij和β的形式分別為：

SGS能量項H的模化如下：

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

SGS湍流擴散項的?；梃bKnight等人［15］的做法，即J＝。

1.3 離散格式和邊界條件

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

程序中以自由來流為初始條件，遠場條件采用基于當地一維Riemann不變量的無反射特征邊界，壁面采用無滑移、無穿透的絕熱邊界條件，展向采用周期性邊界條件。

2 計算結果與分析

2.1 計算參數

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

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

2.2 計算結果驗證

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

表1 當前計算結果和實驗數據的對比Table 1 Comparison of the present calculated results with the experimental data

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

已有研究結果［7，9－10］表明，當M∞＝0.75時圓柱分離點處會形成一道斜激波。關于翼型表面的激波運動已得到廣泛的研究［17］，結果表明存在三種運動模式：正弦類運動（A類）、間歇式運動（B類）和向上游傳播類運動（C類）。然而，關于圓柱表面的激波運動模式研究較少。為了定性地分析圓柱表面的激波運動模式，圖3給出了不同時刻的瞬時流場圖?？梢钥闯觯瑘A柱分離點處的激波在柱體上下表面交替出現，并且向上游傳播，這與翼型表面的C類激波運動模式相似。為了獲得該激波運動的特征St數，在其運動區內的設置一個探測點P1，如圖3（a）所示。圖4給出了P1點處的能譜分布，可以看出激波運動的特征St數與渦脫瀉St數一致，均為0.19。此外，－5／3次方斜率關系表明當前計算結果達到了湍流的慣性子區［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分解的前四個模態空間分布Fig.5 Spatial distributions of first four POD modes based on the pressure fields

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

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

2.4 非定常湍流剪切層行為

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

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

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

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

圖7 P2點處的剪切應力象限分解（a）和瞬時剪切應力角概率密度函數分布（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）

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

圖8 正應力各向異性比例（a）和湍流結構參數（b）沿剪切層變化Fig.8 Evolution of the normal shear stress anisotropy ratio（a）and of turbulence structure parameter（b）along the shear layer

3 結 論

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

［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.