許常悅，鄭 靜，王 哲，王 彬

(南京航空航天大學(xué) 飛行器環(huán)境控制與生命保障工業(yè)和信息化部重點(diǎn)實(shí)驗(yàn)室, 南京 210016)

在工程中，人們會經(jīng)常遇到各類復(fù)雜流動現(xiàn)象，如附體邊界層流動、大范圍分離流、分離剪切層、湍流旋渦及尾跡等.鈍體繞流中包含這些復(fù)雜流動特征，具有較好的代表性，進(jìn)而引起人們的關(guān)注.根據(jù)邊界層分離位置的不同，鈍體可以分為兩類：一是具有移動分離位置；二是固定分離位置.在具有移動分離位置的鈍體中，圓柱屬于典型的規(guī)范構(gòu)型，相關(guān)研究較多，包括不可壓流和可壓縮流動[1-3].另一種具有固定分離位置的鈍體規(guī)范構(gòu)型是方柱.近年來，方柱不可壓繞流的研究報道較多[4-10]，然而缺乏關(guān)于可壓縮流動的研究.

由于方柱不可壓繞流具有豐富的實(shí)驗(yàn)和數(shù)值數(shù)據(jù)，故人們常把該流動問題選為新方法的驗(yàn)證算例.例如，Murakami等[4]針對雷諾數(shù)Re=22 000 的方柱不可壓繞流開展大渦模擬(LES)研究，評估了幾種亞格子模型對計算結(jié)果的影響.結(jié)果表明，基于Lagrange型動力學(xué)模型的計算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)吻合最好.鄧小兵等[5]以方柱不可壓繞流為研究對象，利用基于虛擬壓縮的LES方法進(jìn)行計算，并與實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對比，驗(yàn)證了該方法的可行性.唐鵬等[6]利用兩類RANS/LES(RANS為雷諾平均Navier-Stokes)混合方法模擬了Re=22 000 的方柱不可壓繞流，通過與實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對比，對混合方法的預(yù)測性能進(jìn)行了評估.Minguez等[7]運(yùn)用高階LES方法和實(shí)驗(yàn)手段對Re=21 400 的方柱不可壓繞流進(jìn)行深入研究.他們給出了豐富的定量數(shù)據(jù)，便于其他研究者用于新方法的廣泛驗(yàn)證工作.他們還探討了卡門渦街與自由剪切層Kelvin-Helmholtz (K-H)不穩(wěn)定性的關(guān)系.Trias等[8]采用直接數(shù)值模擬(DNS)方法對Re=22 000 的方柱不可壓繞流進(jìn)行精細(xì)模擬.DNS結(jié)果與已有實(shí)驗(yàn)數(shù)據(jù)吻合較好，能夠捕捉近壁面區(qū)域的轉(zhuǎn)捩現(xiàn)象.Bai等[9]利用LES和分離渦模擬(DES)方法對中等雷諾數(shù)下的方柱不可壓繞流進(jìn)行模擬.與實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對比發(fā)現(xiàn)，DES模型的預(yù)測性能優(yōu)于傳統(tǒng)的LES和RANS模型.李真子等[10]采用自主發(fā)展的LES程序?qū)e=22 000 的方柱繞流進(jìn)行精細(xì)模擬，并與DNS結(jié)果和實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對比，驗(yàn)證了LES方法的可靠性.

目前，關(guān)于方柱可壓縮繞流的文獻(xiàn)較少，且以實(shí)驗(yàn)研究為主.例如，Nakagawa[11]在跨聲速風(fēng)洞中對不同雷諾數(shù)和不同來流馬赫數(shù)(Ma)的方柱流場進(jìn)行了測量.雷諾數(shù)和來流馬赫數(shù)分別為0.696×105

從上述研究可以看出，LES方法適用于研究方柱繞流問題.然而，為了使LES方法準(zhǔn)確地捕捉湍流邊界層，壁面附近仍需較密的網(wǎng)格.為了進(jìn)一步減少計算代價，人們構(gòu)造了多種RANS/LES混合方法，其中以DES方法的應(yīng)用最為廣泛.然而，早期的DES方法存在RANS和LES計算區(qū)域模糊問題，即所謂的“灰區(qū)”問題.不合適的網(wǎng)格設(shè)計，將會導(dǎo)致網(wǎng)格誘導(dǎo)分離現(xiàn)象.人們針對早期DES的不足，提出了各種改進(jìn)策略，如延遲分離DES(DDES)方法[16]及改進(jìn)的延遲分離DES(IDDES)方法[17]等.針對某些特定問題，這些改進(jìn)策略取得了一些成功的計算結(jié)果.然而，DDES和IDDES會引入額外的經(jīng)驗(yàn)參數(shù)，且會增加額外的計算量，這違背了DES發(fā)展的初衷.另外一種能夠?qū)崿F(xiàn)RANS/LES混合計算效果的方法是尺度自適應(yīng)模擬(SAS)方法.目前，應(yīng)用較多的是Menter等[18-19]發(fā)展的兩方程SAS方法.我們曾發(fā)展了一種一方程SAS方法[20-21]，具有計算量小且預(yù)測性能好的優(yōu)點(diǎn).本文將利用該SAS方法作為數(shù)值計算工具，進(jìn)一步評估該方法的預(yù)測能力，這也是本文的另一個研究目的.

1 數(shù)值計算方法

1.1 控制方程

(3)

1.2 尺度自適應(yīng)模型

(4)

(5)

式中：μT為湍流黏性系數(shù)，需要通過構(gòu)造尺度自適應(yīng)模型進(jìn)行求解；Sij=0.5(?ui/?xj+?uj/?xi)為速度應(yīng)變率張量，根據(jù)張量符號運(yùn)算規(guī)則，Skk=S11+S22+S33；PrT為湍流普朗特數(shù)，常近似取值0.92.為了構(gòu)造一方程的尺度自適應(yīng)模型，需要對Spalart-Allmaras (SA)湍流模型方程進(jìn)行修改.SA湍流模型輸運(yùn)方程可以寫成如下形式：

{Cb1[(1-ft2)fv2+ft2]/κ2-Cw1fw}

(6)

(7)

(8)

(9)

(10)

當(dāng)d>Lvk/κ時，SAS執(zhí)行類LES計算；當(dāng)d

1.3 離散格式

SAS的控制方程采用基于結(jié)構(gòu)網(wǎng)格的有限體積方法進(jìn)行求解.時間推進(jìn)采用近似因子分解(AF)方法.為了保證時間的二階精度，在AF方法中引入帶子迭代的偽時間項(xiàng).空間項(xiàng)的離散按照黏性項(xiàng)和對流項(xiàng)分別進(jìn)行處理.黏性項(xiàng)采用傳統(tǒng)的二階精度中心格式進(jìn)行離散.考慮到方柱跨聲速流場中會出現(xiàn)激波和湍流共存現(xiàn)象，對流項(xiàng)的離散格式需要同時滿足捕捉激波和湍流的要求：利用高耗散的格式(如迎風(fēng)格式)捕捉激波.湍流捕捉則采用低耗散格式，如中心格式.我們曾基于Roe通量差分裂格式發(fā)展了一類中心/迎風(fēng)混合格式[22]，中心和迎風(fēng)格式可以通過一個二進(jìn)制開關(guān)函數(shù)進(jìn)行自動切換.該數(shù)值策略已被成功應(yīng)用于多種可壓縮湍流的數(shù)值研究，如鈍柱體跨超聲速繞流[3, 22-23]、翼型跨聲速繞流[24]及超聲速的底部流動[21]等.因此，我們有理由認(rèn)為該數(shù)值策略能夠用于研究當(dāng)前問題.

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

2.1 計算細(xì)節(jié)及計算驗(yàn)證

方柱的來流參數(shù)依據(jù)已有實(shí)驗(yàn)[13]進(jìn)行選取，來流馬赫數(shù)和基于方柱邊長D的雷諾數(shù)分別為0.71和4×105.計算采用O型網(wǎng)格，并在近尾跡區(qū)和近壁處進(jìn)行網(wǎng)格的局部加密處理.O型計算域的直徑為50D，方柱的展向?qū)挾葹?D.為了評估網(wǎng)格分辨率對計算結(jié)果的影響，本文計算采用兩套網(wǎng)格：徑向網(wǎng)格數(shù)×周向網(wǎng)格數(shù)×展向網(wǎng)格數(shù)為197×257×81(網(wǎng)格1)和321×321×121(網(wǎng)格2).圖1給出了網(wǎng)格2在方柱橫截面的網(wǎng)格分布.兩套網(wǎng)格的徑向第一層網(wǎng)格厚度均為10-5D，對應(yīng)的利用壁面單位歸一化的網(wǎng)格間距值小于1，符合SAS模擬邊界層的要求.在程序中，基于網(wǎng)格1和網(wǎng)格2的計算時間步長分別為0.002D/U∞和0.001D/U∞(U∞為來流速度).為便于表述，基于網(wǎng)格1和網(wǎng)格2的SAS結(jié)果分別記為SAS-1和SAS-2.在本文的后處理數(shù)據(jù)中，符號〈 〉用于表示經(jīng)過時間和展向平均的物理量.

圖1 方柱的O型計算網(wǎng)格Fig.1 O-type computational mesh of square cylinder

表1和圖2給出了一些定量數(shù)據(jù)的比較，表中EXP-1和EXP-2分別為來自Layukallo等[13]和Nakagawa[11]的實(shí)驗(yàn)數(shù)據(jù).表1中〈CD〉為平均阻力系數(shù)，CLrms為升力系數(shù)的均方根.圖2中A和D分別為方柱前邊和后邊的中心點(diǎn)，B和C分別為方柱的兩個拐點(diǎn)，〈Cp〉為平均壓力系數(shù)，xs為壁面任意點(diǎn)距A點(diǎn)處的距離.可以看出，SAS-1和SAS-2的結(jié)果相差較小，可以認(rèn)為當(dāng)前計算滿足網(wǎng)格無關(guān)性的要求.然而，為了更好地捕捉精細(xì)的小尺度流動結(jié)構(gòu)，本文所討論和分析的數(shù)據(jù)為SAS-2結(jié)果.從表1中可以看出，SAS預(yù)測的平均阻力系數(shù)略低于實(shí)驗(yàn)值EXP-1和LES結(jié)果.渦脫落的Sr位于實(shí)驗(yàn)EXP-2的可信范圍內(nèi).從圖2中可以看出，在流動分離前，SAS結(jié)果與實(shí)驗(yàn)數(shù)據(jù)EXP-1及LES結(jié)果均吻合較好；在流動分離后，SAS的預(yù)測結(jié)果略高于實(shí)驗(yàn)數(shù)據(jù)EXP-1和LES結(jié)果.對于鈍體跨聲速流動而言，實(shí)驗(yàn)測量值常存在一定幅度的波動[2]，尤其是壁面物理量的測量.本文的SAS結(jié)果與實(shí)驗(yàn)數(shù)據(jù)及LES結(jié)果較為接近，可認(rèn)為當(dāng)前SAS方法具有較好的可靠性.

表1 一些物理量的對比Tab.1 Comparisons of some physical quantities

圖2 沿柱體表面的平均壓力系數(shù)Fig.2 Distributions of mean pressure coefficient on cylinder surface

考慮到流場的快速演化，計算程序選取自由來流參數(shù)作為初始條件.為了避免非物理波對流場信息的污染，遠(yuǎn)場采用基于一維Riemann不變量的無反射邊界條件.壁面速度取無滑移、無穿透條件.壁面溫度取絕熱條件.柱體展向采用周期性邊界條件.為了提高計算效率，程序采用分塊并行計算：周向分為8塊，展向分為3塊，共計24塊.

2.2 計算結(jié)果

圖3(a)、3(b)分別給出了利用平均和瞬態(tài)當(dāng)?shù)伛R赫數(shù)展示的方柱跨聲速流動形態(tài)，圖中ωyD/U∞為展向渦量分量，P1～P4表示4個探測點(diǎn)，坐標(biāo)分別為： P1(-0.055 6, 0.571)，P2(0.721, 0.254)，P3(1.997, 1.485)，P4(-0.214, 0.948).后文將對這4個探測點(diǎn)處的壓力信號進(jìn)行分析與討論.與圓柱繞流類似，方柱的尾跡區(qū)會形成明顯的回流區(qū)，并有大尺度的渦脫落現(xiàn)象.與圓柱不同的是，方柱前緣拐點(diǎn)處會強(qiáng)迫流動發(fā)生分離，進(jìn)而在剪切層與柱體側(cè)壁之間形成回流區(qū).可壓縮效應(yīng)對流場存在明顯的影響，方柱兩側(cè)出現(xiàn)大范圍的局部超聲速區(qū)，且在剪切層下游形成運(yùn)動的激波結(jié)構(gòu)[15].此外，尾跡中也存在局部超聲速區(qū).然而，與圓柱跨聲速繞流[3,22]不同的是，方柱的分離點(diǎn)處未形成激波，這意味著方柱剪切層下游的壓縮波未能向分離點(diǎn)處匯聚.

圖3 方柱流動形態(tài)Fig.3 Flow pattern of square cylinder

對于鈍體繞流而言，其尾跡特性與剪切層的演化行為密切相關(guān).在可壓縮流動中，剪切層中的對流馬赫數(shù)Mc分布常用于分析剪切層的失穩(wěn)過程[25]，Mc的定義為Mc=(〈u〉1-〈u〉2)/(a1+a2)，下標(biāo)1和2分別表示剪切層兩側(cè)的高速和低速流動，〈u〉 和a分別為平均的流向速度和當(dāng)?shù)匾羲?剪切層的增長可以用剪切層的渦量厚度δω描述，δω的定義為δω=(〈u〉1-〈u〉2)/(?〈u〉/?z)max.

圖4 剪切層的Mc和δω沿流向分布Fig.4 Streamwise distributions of Mc and δω inside shear layer

圖5 斜壓項(xiàng)展向分量Fig.5 Isosurface of spanwise component of baroclinic term

圖6 不同截面處的Lamb矢量散度云圖Fig.6 Contours of Lamb vector divergence at different sectional planes

圖7 利用Q準(zhǔn)則繪制的三維渦結(jié)構(gòu)Fig.7 Three-dimensional vortex structures plotted by Q-criterion

為了展示尾跡中的大尺度渦脫落現(xiàn)象，圖7給出了利用Q準(zhǔn)則繪制的三維瞬態(tài)渦結(jié)構(gòu).在剪切層初始發(fā)展階段，剪切層中的渦量分布較為均勻，這與剪切層中的渦量厚度增長較慢有關(guān)(見圖4).隨著剪切層向下游繼續(xù)增長直至失穩(wěn)，尾跡中出現(xiàn)大尺度的發(fā)卡渦和條狀渦結(jié)構(gòu).從圖7可以看出，大尺度渦脫落區(qū)的流動速度仍以亞聲速為主.

方柱中激波運(yùn)動、剪切層失穩(wěn)及渦脫落均會使流場信息發(fā)生規(guī)律變化.為了獲取相關(guān)的特征頻率，我們在方柱的關(guān)鍵位置布置了4個探測點(diǎn)，如圖3(a)所示.借助功率譜密度(PSD)分析，可以獲得這4個探測點(diǎn)處壓力信號的特征頻率，如圖8所示，圖中k為斜率.PSD曲線的峰值對應(yīng)流動的特征頻率，Sr≈0.135屬于柱體的渦脫落特征頻率，而Sr≈0.27則是渦脫落特征頻率的倍頻.圖8(a)中PSD曲線高頻區(qū)呈現(xiàn) -2 次方斜率關(guān)系，這與柱體側(cè)壁與剪切層之間的回流區(qū)存在流動再附現(xiàn)象有關(guān)[15, 28].圖8(b)、8(c)中PSD曲線的高頻區(qū)呈現(xiàn) -5/3 次方斜率關(guān)系，這說明SAS方法能夠模擬湍流的慣性子區(qū)[28]，進(jìn)一步驗(yàn)證了SAS方法的可靠性.圖8(d)中PSD曲線的高頻區(qū)未存在快速衰減，這與P4點(diǎn)處的流動仍以層流流動為主有關(guān).圖8(a)、8(b) 和8(d)中的PSD曲線出現(xiàn)了明顯的倍頻關(guān)系，甚至3倍頻關(guān)系，這與剪切層中存在的渦合并現(xiàn)象有關(guān)[29].同時，這也說明剪切層中的渦合并現(xiàn)象會導(dǎo)致流場中出現(xiàn)向周圍傳播的壓縮波，進(jìn)而使得剪切層外側(cè)附近流場和剪切層內(nèi)側(cè)的回流區(qū)受到影響.然而，圖8(c)中則未出現(xiàn)倍頻關(guān)系.從圖3(a)可以看出，P3探測點(diǎn)位于方柱外側(cè)局部超聲速區(qū)的下游，這說明剪切層激發(fā)的壓縮波未能穿透該局部超聲速區(qū)向下游傳播.

在方柱剪切層的失穩(wěn)過程中，剪切層中發(fā)生渦合并現(xiàn)象，導(dǎo)致激發(fā)的壓縮波向外傳播，剪切層的失穩(wěn)導(dǎo)致尾跡中出現(xiàn)大尺度的渦脫落現(xiàn)象.可以看出，剪切層的失穩(wěn)行為直接影響方柱可壓縮流動的模態(tài).為了分析方柱可壓縮剪切層的失穩(wěn)對流動模態(tài)的影響，借助本征正交分解(POD)技術(shù)對方柱的壓力場進(jìn)行分解.關(guān)于POD技術(shù)的完整描述可以參考文獻(xiàn)[30]，這里僅對POD的定義作簡單說明.對于一個待分解的物理量f(x,t)，POD分析可以確定若干正交函數(shù)φj(x)，j=1, 2,….該物理量f(x,t)在前n個函數(shù)上的投影可以寫成如下形式：

(11)

圖8 4個探測點(diǎn)處的壓力信號PSD曲線Fig.8 PSD curves of pressure signals at four probes

圖9 壓力場的前兩個POD模態(tài) Fig.9 First two POD modes for pressure field

3 結(jié)論

當(dāng)前工作采用SAS方法研究了來流馬赫數(shù)為0.71、雷諾數(shù)為4×105的方柱跨聲速繞流問題.通過對計算結(jié)果進(jìn)行分析與討論，可以得到如下的研究結(jié)論：

(1) 通過與LES結(jié)果和實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對比，可以看出當(dāng)前SAS方法具有較好的可靠性，能夠較好地模擬具有固定分離點(diǎn)的鈍體可壓縮繞流問題.

(2) 方柱剪切層中的對流馬赫數(shù)約為0.6，這意味著Kelvin-Helmholtz不穩(wěn)定性主導(dǎo)分離剪切層的演化過程.在剪切層的初始階段，剪切層中的擾動渦展向呈現(xiàn)滾筒狀；Lamb矢量散度呈現(xiàn)負(fù)值包裹正值的“夾心三層”結(jié)構(gòu)，抑制了渦量厚度的增長.

(3) 通過對方柱流場中布置一些壓力信號探測點(diǎn)，借助功率譜密度分析方法，發(fā)現(xiàn)剪切層外側(cè)附近和方柱的回流區(qū)均出現(xiàn)倍頻現(xiàn)象，這與剪切層中存在明顯的渦合并現(xiàn)象有關(guān).

(4) 壓力場的POD分析表明，方柱跨聲速流場中的主導(dǎo)流動模態(tài)為反對稱模態(tài)，這與尾跡中的渦脫落現(xiàn)象和剪切層激發(fā)的壓縮波傳播有關(guān).