涂國華, 陳堅強,,*, 袁先旭, 楊 強, 張毅鋒

(1.空氣動力學(xué)國家重點實驗室, 四川 綿陽 621000; 2.中國空氣動力研究與發(fā)展中心 計算空氣動力研究所, 四川 綿陽 621000)

0 引 言

高超聲速情況下湍流的摩阻和熱流通常是層流的3-5倍[1]，湍流/層流轉(zhuǎn)捩問題嚴(yán)重影響高超聲速飛行器的氣動特性和熱防護系統(tǒng)。高超聲速情況下邊界層流動容易經(jīng)歷層流/湍流轉(zhuǎn)捩[2]，發(fā)展高超聲速邊界層轉(zhuǎn)捩控制方法，推遲邊界層流動從層流轉(zhuǎn)捩成湍流，對減阻和熱防護具有重要意義。但是高超聲速邊界層轉(zhuǎn)捩現(xiàn)象非常復(fù)雜，影響因素眾多，轉(zhuǎn)捩控制極為困難。在轉(zhuǎn)捩被動控制方法中，國外的研究表明，多孔表面、微槽道和超聲波吸聲材料或許可用來推遲轉(zhuǎn)捩[3-7]。但是，這類控制方法的風(fēng)洞實驗難度大成本高，數(shù)值模擬也面臨著算法、網(wǎng)格和計算量等方面的諸多困難。因此，通過簡化模型、采用線性穩(wěn)定性理論(LST)從轉(zhuǎn)捩機理的角度對這類控制方法開展研究是目前少數(shù)可用的快速且低成本的技術(shù)途徑之一。Fedorov[3]和Kozlov[4]等通過聲波在無窮長管道中的傳播關(guān)系對多孔表面的擾動邊界條件進行了簡化，Sandham[5]和Wartemann[6]等通過時間模式的直接數(shù)值模擬驗證了這種簡化方法是可行的。我們采用Fedorov[3]和Kozlov[4]多孔表面邊界處理方法，分析了開孔率和孔半徑對第二模態(tài)擾動波的影響，得到了抑制擾動波幅值增長的最優(yōu)開孔率和孔半徑。

1 線性穩(wěn)定性理論

在邊界層穩(wěn)定性和轉(zhuǎn)捩研究中，LST被廣泛用于研究邊界層中小擾動演化情況(比如幅值增長率)，然后可通過半經(jīng)驗eN方法預(yù)測轉(zhuǎn)捩。

設(shè)流動可用基本流疊加擾動構(gòu)成，如

(1)

其中φ0=[ρ0,u0,v0,w0,T0]T為基本流，φ′=[ρ′,u′,v′,w′,T′]T為小擾動。把疊加了擾動的流場代入Navier-Stokes方程，減去基本流滿足的部分，略去二階及高階,得到線化擾動方程：

(2)

并假定基本流是平行流(即v=0)，并進一步假設(shè)擾動可表示成如下行波的線性疊加。

(3)

式中c.c.表示對應(yīng)的復(fù)共軛。把式(3)代入線化擾動方程，可得到常微分方程:

(4)

其中系數(shù)矩陣

C0=D-iωΓ+iαA+iβC+α2Vxx+

β2Vzz+αβVxz

C1=B-iαVxy-iβVyz

C2=-Vyy

可見C0～C2與流向波數(shù)α、展向波數(shù)β、圓頻率ω、和式(2)的線性系數(shù)矩陣Γ、A、B等有關(guān)。當(dāng)給定基本流場后，式(2)的線性系數(shù)矩陣就已知了，式(4)所表達的是一個特征值問題，只有當(dāng)α、β和ω滿足一定關(guān)系時才有特征解，這個關(guān)系被稱為色散關(guān)系式。穩(wěn)定性分析時可給定其中兩個參數(shù)，只需求解剩下的一個。實際求解通常分時間模式和空間模式兩種情況：對于時間模式的擾動，α和β為實數(shù)，ω=ωr+iωi為復(fù)數(shù)；對于空間模式，ω為實數(shù)，α和β為復(fù)數(shù)。由于時間模式比空間模式求解更方便，本文僅考慮了時間模式，ωi表示擾動波的幅值隨著時間按eωi指數(shù)型增長/衰減情況。

通過合適的邊界條件(下文討論)，便可求解式(4)。求解時可以把它轉(zhuǎn)化成單變量的8次方程，也可直接求解，通常采用高階精度方法。本文采用Chebychev點配置法對式(4)進行離散。

2 邊界條件

對于光滑壁面，壁面取無滑移邊界和等溫邊界，即：

u′=v′=w′=T′=0

(5)

對于遠場，施加速度及溫度擾動等于零的條件。

對于多孔介質(zhì)表面，遠場邊界條件與光滑壁面相同，但是壁面條件需要考慮孔壁對擾動波的干擾作用。假設(shè)多孔介質(zhì)的孔徑遠遠小于邊界層厚度，對邊界層宏觀流動(基本流場)無影響，但對聲波、T-S波等的影響不可忽略。同時假設(shè)孔的深度遠大于半徑，于是擾動在孔中的傳播可以簡化為聲波在無窮長管道中的傳播，根據(jù)Fedorov[3]和Kozlov[4]等的研究，壁面上的邊界條件為:

(6)

其中系數(shù)A和B通過聲波在無窮長管道中的傳播關(guān)系來確定，

(7)

其中,

上面各式中：n為開孔率，定義為控制區(qū)域內(nèi)孔面積所占總面積的比例(由于孔的間距不能小于孔的直徑，所以最大開孔率為π/4)；r為孔的半徑，在給定開孔率的情況下，孔半徑越小，表示單位面積上孔的個數(shù)越多；h為孔的深度；J0和J2為第一類Bessel方程；Me為邊界層外緣馬赫數(shù)；Re為雷諾數(shù)；Pr為普朗特數(shù)；γ為比熱比；下標(biāo)“w”表示物面。

ω既是時間模式的解，又是邊界條件的輸入?yún)?shù)，采用迭代方法求解ω，當(dāng)?shù)`差小于1×10-8時，認為結(jié)果已經(jīng)收斂。

由于高超聲速平板邊界層中第二模態(tài)擾動波占主導(dǎo)作用，且第二模態(tài)擾動波中的二維擾動增長最快，所以我們僅分析二維擾動波，即展向波數(shù)β=0。

3 多孔表面對馬赫數(shù)6平板邊界層穩(wěn)定性的影響

3.1 基本流場和方法考核

邊界層基本流動的參數(shù)與Sandham等[5]相同，馬赫數(shù)為6.0，來流溫度216.65 K，壁面溫度1522.44 K。但是，本文獲取基本流的方法與Sandham等的方法不同。Sandham等在邊界層相似解的假設(shè)下采用Blasius解，而本文采用基于WCNS格式[8]的高階精度CFD方法[9-10]求解Navier-Stokes方程。本文方法的優(yōu)點是更接近實際情況，因為實際流動中，平板前緣因邊界層開始生成，會出現(xiàn)一道微弱的斜激波(圖1)，而Blasius相似解忽略了這道斜激波。

圖1 邊界層法向速度云圖Fig.1 Wall-normal velocity of the boundary layer

為了便于與文獻比較，取排移厚度等于2 mm的邊界層剖面進行穩(wěn)定性分析，并以此尺度作為全文計算的參考長度。由于文獻沒有考慮法向速度，本文首先強制法向速度為零，通過與文獻比較考察本文方法的正確性。圖2給出了流向波數(shù)為α=2π/3的擾動波形狀函數(shù)，可見與文獻[5]非常一致。圖3在光滑壁面上比較了本文計算的擾動增長率和文獻[6](文獻[5]和文獻[6]為同一邊界層)的結(jié)果，可見二者基本一致。圖4比較了光滑表面和多孔介質(zhì)表面上的擾動波幅值增長率，還同時給出了文獻[6]的計算結(jié)果，可見本文結(jié)果與文獻結(jié)果一致，都表明：多孔介質(zhì)能大大降低第二模態(tài)的最大增長率；多孔介質(zhì)對第一模態(tài)有輕微的放大作用。

圖2 α=2π/3擾動波的形狀函數(shù)Fig.2 Shape function of the α=2π/3 wave

圖3 光滑表面上擾動波的增長率Fig.3 Increase rates of disturbances in theboundary layer of smooth wall

仔細觀察圖3和圖4可以發(fā)現(xiàn)，當(dāng)流向波數(shù)較小時，本文結(jié)果與文獻結(jié)果幾乎完全一致；當(dāng)流向波數(shù)較大時(比如α>2)時，本文結(jié)果與文獻有輕微差別，這種差別可能是由基本流速度剖面的細微差別引起的。文獻中的速度剖面采用的是Blasius相似解，而本文采用的是DNS結(jié)果(但在LST分析時法向速度強制賦零)。由于法向速度很小，所以對波長較長(即波數(shù)較小)的擾動波的影響非常小，僅對波長較短(即波數(shù)較大)的擾動波有輕微影響。下文還將進一步討論基本流的法向速度對多孔介質(zhì)控制效果的影響。

3.2 最優(yōu)開孔率和孔半徑

文獻[6]研究表明，當(dāng)孔深超過0.8倍邊界層排移厚度時，控制效果基本上不再隨著孔的深度變化。本文不研究孔深的影響，而是給定孔深為2 mm(等于x=235 mm處的邊界層排移厚度)，研究開孔率和孔半徑的影響。

圖4 多孔介質(zhì)表面與光滑表面的比較Fig.4 Comparison between the porous wall and smooth wall

我們早期的考察結(jié)果表明[11]，并不是開孔率越大越好，最不穩(wěn)定擾動波的增長率并不是隨著開孔率單調(diào)變化，而是存在最優(yōu)開孔率。早期考察還表明，控制效果隨孔半徑的變化趨勢也不是單調(diào)的，而是存在最優(yōu)孔半徑。為了尋找最優(yōu)的控制策略，我們在較為廣泛的控制參數(shù)范圍內(nèi)全局尋優(yōu)，其中開孔率n的范圍為0.05～0.75；孔半徑r的范圍為0.01～0.3 mm。圖5給出了距離平板前緣235 mm處最不穩(wěn)定第二模態(tài)擾動增長率ωi隨開孔率n和孔半徑r的變化情況?？梢娮顑?yōu)控制的組合參數(shù)是n≈0.35，r≈0.20 mm。此時最不穩(wěn)定時間模式擾動波的增長率僅為光滑平板情況下的38%，失穩(wěn)擾動波的增長率大幅減小，有利于推遲轉(zhuǎn)捩。

圖5 x=235 mm處最不穩(wěn)定第二模態(tài)擾動波增長率ωi與開孔率n和孔半徑r的關(guān)系Fig.5 Increase rates of the most unstable mode-2 wavesin the porosity n and pore radius r coordinates

圖6比較了多孔介質(zhì)表面和光滑表面的最不穩(wěn)定第二模態(tài)擾動波的幅值增長率。通過采用多孔介質(zhì)，擾動波開始增長的位置(中性點)從距離平板前緣24 mm推遲到了50 mm的地方，推遲比例約1倍。不僅如此，擾動波的幅值增長率也大幅減小。在所考察的1000 mm長的平板范圍內(nèi)，擾動波的幅值增長率不超過光滑平板的40%。圖6還表明，除了頭部區(qū)域，多孔平板情況下最不穩(wěn)定擾動波增長率沿流向單調(diào)降低，與光滑平板類似。在頭部區(qū)域，擾動波增長率被抑制的比例更大，這說明把多孔介質(zhì)安裝在擾動波增長最快的區(qū)域可達到最佳控制效果。

圖6 最不穩(wěn)定第二模態(tài)擾動波增長率沿流向比較Fig.6 Increase rates of the most unstable mode-2waves along the wall

圖7給出了不同流向位置上的最優(yōu)開孔率和孔半徑。在x<600 mm的地方，開孔率沿流向并不是單調(diào)變化，在x=450 mm的地方達到最大值0.66左右。最優(yōu)孔半徑的變化范圍比較小，最大值為0.22，最小值為0.19，大部分都在0.21 mm左右，這對設(shè)計多孔壁非常有利。

在最優(yōu)多孔控制參數(shù)情況下，最不穩(wěn)定擾動波的圓頻率有所降低，接近光滑壁第二模態(tài)中性曲線的下支(圖8)。從圖9可以看出，多孔壁使最不穩(wěn)定第二模態(tài)擾動波的相速度和流向波數(shù)輕微減小。

圖10給出了擾動波長與最優(yōu)孔半徑的關(guān)系，可見最不穩(wěn)定擾動波長為最優(yōu)孔半徑的20倍以上，且隨著下游發(fā)展，擾動波長與最優(yōu)孔半徑之比持續(xù)增加，這可能是由邊界層增厚引起的。邊界層越厚，最不穩(wěn)定擾動波的波長越長。定義如下衡量擾動波波長與孔半徑的無量綱參數(shù):

(8)

其中,λ為擾動波波長，r為孔半徑，Rex為當(dāng)?shù)乩字Z數(shù)(參考量長度為當(dāng)?shù)刂燎把氐木嚯x，其他參考量為邊界層外緣參數(shù))。從圖10可以看出，除頭部區(qū)域外，τ為0.02左右。

圖8 第二模態(tài)擾動波圓頻率沿流向分布Fig.8 Circular frequencies of mode-2 waves along the wall

圖9 最不穩(wěn)定第二模態(tài)的相速度c和流向波數(shù)αFig.9 Phase speeds and stream-wise wavenumbers ofmost unstable mode-2 waves

圖10 最優(yōu)孔半徑與最不穩(wěn)定擾動波長之間的關(guān)系Fig.10 Optimal pore radius vs.stream-wise wavelength

3.3 當(dāng)?shù)鼗玖鞣ㄏ蛩俣鹊挠绊?/h3>

邊界層流動大多不是嚴(yán)格的平行流，但是上述分析采用了平行流假設(shè)(v0=0)，這會帶來一定的誤差；另外，多孔壁的擾動邊界條件(6)中法向擾動不等于零，所以有必要分析基本流非平行性的影響。我們保留高階精度CFD計算得到的基本流的法向速度(v0≠0)。即對基本流的速度進行如下擴展:

(9)

從嚴(yán)格意義上理解，在非平行流假設(shè)下并不能把小擾動表示成式(3)的行波解形式。幸運的是大部分邊界層的法向速度相對于流向速度是一個小量，式(3)仍然近似可行，且在實驗中已經(jīng)普遍證明了邊界層確實存在行波(比如第一、二模態(tài)擾動波)。這說明仍然可在LST框架下保留法向速度。如果忽略形狀函數(shù)和空間波數(shù)的流向變化，則仍然可以把擾動方程化為與式(4)類似的常微分方程，其邊界條件和求解方法與傳統(tǒng)LST一樣。這樣處理的優(yōu)點是可以直接借用傳統(tǒng)LST的分析方法和經(jīng)驗，缺點是僅能考慮當(dāng)?shù)胤瞧叫行缘挠绊懀荒荏w現(xiàn)上下游的關(guān)系，是介于傳統(tǒng)LST和拋物化穩(wěn)定性分析方法之間的一種分析方法。拋物化穩(wěn)定性分析方法保留了法向速度和流向1/Re量級的變化量。

圖11比較了法向速度對控制效果的影響。法向速度對最優(yōu)開孔率和孔半徑的影響非常小，除了頭部區(qū)域，對最優(yōu)開孔率的影響小于10%，對最優(yōu)孔半徑的影響小于5%。法向速度對控制效果有一定的影響，在頭部區(qū)域，擾動波增長率可減小超過30%，但是隨著向下游發(fā)展，考慮法向速度與不考慮法向速度的差異逐漸減小?？傊?，基于平行流假設(shè)得到的最優(yōu)開孔率和孔半徑仍然有效，對控制效果的估計偏保守。

圖11 基本流法向速度對控制效果的影響(圖中線條表示不考慮法向速度的結(jié)果，符號表示考慮法向速度的結(jié)果)Fig.11 Comparison between the results with andwithout wall-normal velocity of base flow (The symbols are results with wall-normal velocity; the lines are results without wall-normal velocity.)

4 結(jié)果討論

當(dāng)參數(shù)合適時，多孔介質(zhì)表面對第二模態(tài)擾動波的幅值增長率具有抑制作用。多孔介質(zhì)的關(guān)鍵參數(shù)有三個：開孔率、孔半徑和孔深。文獻中采用時間模式DNS對馬赫數(shù)6平板邊界層的分析表明，當(dāng)孔深超過邊界層排移厚度的0.8倍時，控制效果基本不再隨孔深變化。但是文獻還未研究開孔率和孔半徑對控制效果的影響。我們采用時間模式的線性穩(wěn)定性理論，在孔深為2 mm的情況下，分析了開孔率和孔半徑對控制效果的影響。在分析時，通過聲波在無限長管道中的傳播關(guān)系來確定LST的多孔壁面條件。

研究發(fā)現(xiàn)，采用多孔介質(zhì)表面可以大大推遲第二模態(tài)擾動波的中性點，對于本文的馬赫數(shù)6平板邊界層，中性點距離平板前緣的位置從24 mm推遲到了50 mm。最優(yōu)開孔率與流向位置緊密相關(guān)，大致變化范圍是0.31至0.66之間。最優(yōu)孔半徑與流向位置關(guān)系不大，基本在0.2 mm左右波動。通過最不穩(wěn)定擾動波長、最優(yōu)孔半徑和當(dāng)?shù)乩字Z數(shù)組成的無量綱參數(shù)τ約等于0.02。

我們還考察了基本流的法向速度對擾動波的影響。與平行流假設(shè)相比，法向速度對最優(yōu)開孔率和孔半徑的影響較小，對控制效果的影響僅在頭部區(qū)域比較明顯，但是隨著下游發(fā)展，影響趨勢逐漸減弱。

上述結(jié)果是在多孔表面對基本流的影響可以忽略的假設(shè)下采用時間模式LST得到的。當(dāng)開孔率較大時，多孔表面對基本流的影響將變得重要起來；另外，當(dāng)孔半徑較大時，多孔介質(zhì)本身還會誘導(dǎo)聲波；當(dāng)孔半徑較小時，需要考慮稀薄氣體效應(yīng)，文獻[12]發(fā)現(xiàn)稀薄效應(yīng)可以增強多孔介質(zhì)的控制效果，但本文還未開展相關(guān)研究。由于空間模式LST與真實流動更一致，建議后續(xù)采用空間模式LST開展研究。

致謝:在研究過程中得到了天津大學(xué)曹偉教授、劉建新博士和趙磊博士的指點和幫助，對他們表示感謝。

參考文獻：

[1]Li F, Xie S F, Bi Z X.Experimental study of several on aerodynamicproblems on hypersonic vehicles[J].Modern Defence Technology, 2014, 42(5): 1-7.(in Chinese)

李鋒, 解少飛, 畢志獻, 等.高超聲速飛行器中若干氣動難題的實驗研究[J].現(xiàn)代防御技術(shù), 2014, 42(5): 1-7.

[2]Chen J Q, Tu G H, Zhang Y F, et al.Hypersnonic boundary layer transition: what we know, where shall we go[J].Acta Aerodynamics Sinica, 2017, 6, 35(3): 311-337.(in Chinese)陳堅強, 涂國華, 張毅鋒, 等.高超聲速邊界層轉(zhuǎn)捩研究現(xiàn)狀與發(fā)展趨勢[J].空氣動力學(xué)學(xué)報, 2017, 35(3): 311-337.

[3]Alexander V Fedorov.Stabilization of hypersonic boundary layers by porous coatings[J].AIAA J, 2001, 39(4): 605-610.

[4]Kozlov V F, Fedorov A V, Malmuth N D.Acoustic properties of rarefied gases inside pores of simple geometries[J].J.Acoustical Society of America, 2005, 117(6): 3402-3412.

[5]Neil D Sandham, Heinrich Lüdeke.Numerical study of Mach 6 boundary-layer stabilization by means of a porous surface[J].AIAA J, 2009, 47(9): 2243-2252.

[6]Wartemann W, Lüdeke H, Sandham N D.Numerical investigation of hypersonic boundary-layer stabilization by porous surfaces[J].AIAA Journal, 2012, 50(6): 1281-1290.

[7]Nicola DeTullio, Neil D Sandham.Direct numerical simulation of breakdown to turbulence in a Mach 6 boundary layer over a porous surface[J].Physics of Fluids, 2010, 22: 094105.

[8]Deng X G, Mao M L, Tu G H, et al.Extending weighted compact nonlinear schemes to complex grids with characteristic-based interface conditions[J].AIAA Journal, 2010, 48(12): 2840-2851.

[9]Tu G H, Deng X G, Min Y B, et al.Method for evaluating spatial accuracy order of CFD and applications to WCNS scheme on four typically distorted meshes[J].Acta Aerodynamica Sinica, 2014, 32(04): 425-432.(in Chinese)涂國華, 鄧小剛, 閔耀兵, 等.CFD空間精度分析方法及4種典型畸形網(wǎng)格中WCNS格式精度測試[J].空氣動力學(xué)學(xué)報, 2014, 32(4): 425-432.

[10]Tu G H, Deng X G, Mao M L.Spectral property comparison of fifth-order nonlinear WCNS and WENO difference schemes[J].Acta Aerodynamica Sinica, 2012, 30(6): 709-712.(in Chinese)涂國華, 鄧小剛, 毛枚良.5階非線性WCNS和WENO差分格式頻譜特性比較[J].空氣動力學(xué)學(xué)報, 2012, 30(6): 709-712.

[11]Tu G H, Chen J Q, Yuan X X, et al.Stabilizing a Mach 6 flat plate boundary layer by porous wall[C]//The 17th national conference on CFD, 2017, Hangzhou, China, CARS-2017-03-165, 141-146.(in Chinese)涂國華, 陳堅強, 袁先旭, 等.多孔介質(zhì)對馬赫6平板邊界層穩(wěn)定性的影響[C]//第十七屆全國計算流體力學(xué)會議, 2017, 杭州: CARS-2017-03-165, 141-146.

[12]Fedorov A, Kozlov V, Shiplyuk A, et al.Stability of hypersonic boundary layer on porous wall with regular microstructure[J].AIAA J, 2006, 44(8): 1866-1871.