鈍頭體中的廣義雷諾比擬關(guān)系1)

陳星星 陳 皓 范晶晶 溫玉芬 張 正 馬友林

(中國運載火箭技術(shù)研究院，北京 100076)

引言

雷諾比擬描述了流體在壁面上的剪切和傳熱之間的關(guān)聯(lián)關(guān)系，最早由雷諾在1874 年提出[1]，在各類流體傳熱傳質(zhì)問題中受到廣泛關(guān)注和研究.在邊界層理論研究的早期，Blasius[2]在不可壓縮的半無窮長平板繞流中找到了邊界層方程的第一個理論解，從中可以計算得到平板壁面的雷諾比擬關(guān)系，結(jié)果表明摩阻系數(shù)和熱流系數(shù)的比值，即雷諾比擬系數(shù)，為常數(shù)1/2.

不可壓縮平板繞流中的雷諾比擬關(guān)系具有十分簡單的形式，直觀地揭示了邊界層中動量和能量輸運之間的關(guān)系.更進(jìn)一步的研究表明，即使考慮更加復(fù)雜的效應(yīng)，如可壓縮流動[3-4]、尖前緣的稀薄氣體效應(yīng)[5-6]、燃燒[7-8]等，平板繞流中的雷諾比擬系數(shù)仍然為常數(shù).

對于工程實踐中更常見的湍流邊界層流動，雷諾比擬也受到廣泛的關(guān)注和應(yīng)用.文獻(xiàn)[9-14]利用雷諾比擬關(guān)系構(gòu)建解析理論預(yù)示平板和尖錐表面的熱流分布，這些理論方法在工程中得到了廣泛的應(yīng)用.Chi 等[15]利用動量方程的湍流模型和雷諾比擬關(guān)系構(gòu)造了能量方程的湍流模型.作為目前湍流流動中少數(shù)的理論結(jié)果之一，雷諾比擬關(guān)系還被用來與各類數(shù)值計算、試驗測量結(jié)果進(jìn)行對比驗證[16-18].

近年來，陳星星等[19]進(jìn)一步對彎曲壁面上的雷諾比擬關(guān)系開展了研究，提出了廣義雷諾比擬以描述雷諾比擬系數(shù)與壁面外形之間的關(guān)系.傳統(tǒng)的雷諾比擬關(guān)系適用于平板、尖錐繞流等邊界層流向壓力梯度為零的流動，而對于存在流向壓力梯度的流動，則熱流和摩阻的比擬關(guān)系一般不是常數(shù)[20].以二維的圓柱繞流和軸對稱的圓球繞流為例，其熱流峰值處于駐點位置，沿壁面向下游單調(diào)下降，而摩阻則在駐點處為零，沿壁面向下游先增加后減小，二者的比值隨壁面位置變化.

陳星星等通過對鈍頭體壁面附近層流邊界層方程的理論分析，推導(dǎo)出在高超聲速條件下，壁面摩阻和熱流的雷諾比擬系數(shù)隨與壁面當(dāng)?shù)貎A斜角(壁面當(dāng)?shù)厍芯€與來流方向夾角的余角)呈正比關(guān)系，其比例系數(shù)與壁溫相關(guān)而與馬赫數(shù)、雷諾數(shù)無關(guān).理論分析結(jié)果與DSMC 數(shù)值仿真計算結(jié)果吻合[19].進(jìn)一步地，陳星星等[21]還將鈍頭體中的廣義雷諾比擬關(guān)系推廣到的包含化學(xué)反應(yīng)流動的情形.

廣義雷諾比擬關(guān)系理論與DSMC 數(shù)值模擬計算結(jié)果吻合較好.但是受DSMC 方法的限制[22-24]，主要適用于稀薄氣體或近連續(xù)區(qū)流動的計算，雷諾數(shù)較低，而對于工程中常見的較高雷諾數(shù)流動，乃至湍流流動，在當(dāng)前計算機(jī)硬件水平下，采用DSMC 方法很難實現(xiàn).

另一方面，高超聲速鈍頭體繞流問題是高速氣體動力學(xué)中的經(jīng)典問題，在各類航天飛行器中廣泛存在，長期以來受到深入研究.文獻(xiàn)[25-26]分別研究了鈍頭體壁面和駐點熱流問題，給出了理論預(yù)示方法.國內(nèi)王智慧等[27-29]提出了駐點熱流受稀薄氣體效應(yīng)影響的理論判據(jù).在數(shù)值求解N-S 方程方面，文獻(xiàn)[30-32]研究了不同離散格式、不同數(shù)值網(wǎng)格對鈍頭體熱流計算結(jié)果的影響.Schwartzentruber 等[33]對比了采用DSMC 方法和數(shù)值求解N-S 方法計算得到的圓柱表面摩阻和熱流大小.

為了對鈍頭體中的廣義雷諾比擬關(guān)系開展進(jìn)一步研究，本文采用數(shù)值求解N-S 方程的方法計算了不同雷諾數(shù)、不同外形條件下的鈍頭體繞流問題，并將壁面上的雷諾比擬系數(shù)分布與理論預(yù)示結(jié)果進(jìn)行對比，為廣義雷諾比擬關(guān)系研究提供數(shù)據(jù)支撐.

1 圓柱壁面的廣義雷諾比擬關(guān)系

1.1 理論預(yù)示

分別定義摩阻系數(shù)Cf=和熱流系數(shù)Ch=，其中τw為壁面剪切力，qw為壁面熱流，為來流密度，為來流速度.

對于如圖1 所示的高超聲速繞鈍頭體流動，壁面摩阻和熱流系數(shù)的比值為

其中，Pr為氣體普朗特常數(shù)，ue為沿流向的邊界層外緣速度，分別為壁面上氣流無量綱速度和溫度的法向梯度.式中摩阻和熱流系數(shù)的比值Cf/Ch即為雷諾比擬系數(shù).

根據(jù)可壓縮流動的Bernoulli 方程得到的鈍頭體邊界層外緣速度分布公式[34]

式中，γ 為氣體常數(shù)，pw為壁面壓力，p0為駐點壓力.在來流馬赫數(shù)的條件下，上式可以近似為關(guān)于壁面當(dāng)?shù)貎A斜角θ 的線性分布[35]:，θ 的定義見圖1.此時可以得到線性的廣義雷諾比擬關(guān)系式[19]

其中，Cr為一個與壁面溫度相關(guān)的比例系數(shù)，根據(jù)前期研究結(jié)果,對于二維圓柱外形

對于軸對稱圓球外形

式中，Tw為壁面溫度，T0為來流氣體總溫.

圖1 鈍頭體繞流示意圖Fig.1 Hypersonic flow around a blunt-nosed body

1.2 N-S 方程數(shù)值計算

為了對更一般條件下的廣義雷諾比擬問題開展研究，本文采用數(shù)值求解N-S 方程的方法模擬圓柱繞流問題，計算壁面的摩阻分布和熱流分布.采用不同的數(shù)值方法進(jìn)行計算，并與理論預(yù)示以及文獻(xiàn)[32]中公開發(fā)表的數(shù)值計算結(jié)果進(jìn)行對比.計算工具為CFD++軟件(TVD 格式)[36].物理模型為以空氣為介質(zhì)的理想氣體定常層流流動.

為了方便對比，計算參數(shù)與文獻(xiàn)[32]保持一致：圓柱半徑8 cm，來流條件等效為海拔70 km 高度處大氣條件：來流靜壓p=4.85 Pa,來流密度=7.48×10?5kg/m3，來流馬赫數(shù)為3,6,12，其中不同馬赫數(shù)對應(yīng)的壁面溫度Tw依次分別為300 K,500 K,1000 K.

計算使用對稱半模型結(jié)構(gòu)化網(wǎng)格，如圖2 所示，在圓柱θ=90?的位置和計算出口邊界之間使用一段平直壁面作為過渡.網(wǎng)格數(shù)量為261(沿壁面)×125(壁面法向).沿壁面法向的第一層網(wǎng)格高度為1.5×10?5m.

圖2 計算域與網(wǎng)格示意圖Fig.2 Computation regime and mesh

不同馬赫數(shù)條件下圓柱壁面的雷諾比擬系數(shù)計算結(jié)果見圖3.根據(jù)Schwartzentruber 等[32]的計算結(jié)果，相比于DSMC 方法，采用N-S 方程計算得到的駐點熱流偏大3.5%～10.1%，表明在70 km 高度處，8 cm 半徑的圓柱繞流存在局部稀薄氣體效應(yīng)，處于近連續(xù)流動區(qū).對于摩阻和熱流的雷諾比擬系數(shù)Cf/Ch，采用DSMC 方法和N-S 方程計算得到的結(jié)果基本一致.

本文采用TVD 格式數(shù)值求解N-S 方程的結(jié)果也同時顯示在圖中.可見3 種數(shù)值方法計算得到的雷諾比擬關(guān)系較為一致.在高超聲速(=6,12)條件下，不同數(shù)值方法，包括DMSC 方法以及不同數(shù)值格式求解N-S 方程計算得到的雷諾比擬關(guān)系均與理論預(yù)示的結(jié)果吻合較好.

圖3 不同馬赫數(shù)條件下的雷諾比擬關(guān)系Fig.3 General Reynolds analogy relation in different Mach number

對不同來源、不同方法的計算數(shù)據(jù)對比分析發(fā)現(xiàn)，對于超聲速圓柱繞流中的雷諾比擬關(guān)系，采用基于粒子碰撞和統(tǒng)計抽樣的DSMC 方法和基于連續(xù)介質(zhì)假設(shè)求解N-S 控制方程方法計算的結(jié)果基本一致.在高超聲速條件下，數(shù)值計算結(jié)果與理論預(yù)示結(jié)果吻合較好，而在較低馬赫數(shù)(=3)條件下，數(shù)值計算結(jié)果與理論預(yù)示結(jié)果則存在一定區(qū)別，表明當(dāng)前的高超聲速廣義雷諾比擬關(guān)系的理論結(jié)果還需要進(jìn)一步改進(jìn)才可以適應(yīng)一般超聲速流動的情形.

為了對計算結(jié)果進(jìn)行確認(rèn)，檢查了本文采用的TVD 格式計算結(jié)果的收斂性以及不同網(wǎng)格數(shù)目對計算結(jié)果的影響.不同馬赫數(shù)條件下駐點熱流以及下游(θ=90?)位置的熱流計算結(jié)果收斂歷程見圖4.圖中縱坐標(biāo)為各時間步計算得到的熱流值與最后一步計算的熱流值之間的相對誤差：（ .可見，計算過程中熱流值穩(wěn)步收斂，最終收斂結(jié)果的波動幅度小于10?7.使用不同密度網(wǎng)格計算得到的雷諾比擬系數(shù)結(jié)果見圖5.三套網(wǎng)格對應(yīng)的壁面附近第一層網(wǎng)格高度由疏到密各自為2.2×10?5m,1.5×10?5m 和1.0×10?5m.可見，在不同的馬赫數(shù)條件下，網(wǎng)格數(shù)目對雷諾比擬系數(shù)的影響較小.

圖4 熱流系數(shù)計算相對誤差隨迭代步數(shù)的變化Fig.4 Heat Fluxes vary with iteration steps

圖5 不同網(wǎng)格數(shù)目下的雷諾比擬系數(shù)計算結(jié)果Fig.5 Reynolds analogy coefficients under different amount of mesh

1.3 雷諾數(shù)對廣義雷諾比擬關(guān)系的影響

在文獻(xiàn)[32]提供的來流條件下，數(shù)值求解NS 方程以及DSMC 方法的雷諾比擬計算結(jié)果與理論預(yù)示結(jié)果均吻合較好.但是計算條件的來流密度和來流雷諾數(shù)較低.參考王智慧等[27,29]提出的駐點稀薄氣體流動判據(jù)Wr=，對于海拔高度H=70 km，參考長度8 cm 的計算條件，馬赫數(shù)=3，6，12 時，雷諾數(shù)分別為398，796，1592；Wr分別為0.013，0.018，0.026.對于計算的3 種狀態(tài)，Wr參數(shù)的大小表明流動處于近連續(xù)區(qū).而在常見的工程問題中，來流雷諾數(shù)往往大若干量級.為了在更一般的流動問題中考察廣義雷諾比擬關(guān)系的適應(yīng)性，在上述計算條件的基礎(chǔ)上，計算了更大雷諾數(shù)條件下的圓柱繞流問題，計算條件見表1.計算選擇高超聲速的=6 和12 狀態(tài)，在維持來流溫度以及壁面溫度條件下，通過提高來流壓力條件增大雷諾數(shù).

表1 計算條件Table 1 Simulation cases

圖7 中顯示了不同雷諾數(shù)條件下，駐點下游熱流和摩阻系數(shù)的分布曲線.圖中采用駐點熱流系數(shù)和最大摩阻系數(shù)對熱流系數(shù)和摩阻系數(shù)作歸一化處理.由圖可見，不同來流壓力條件下，沿圓柱壁面的熱流和摩阻分布規(guī)律基本一致，而在駐點下游較遠(yuǎn)的位置(θ >60?)，不同來流壓力條件下，歸一化的熱流和摩阻曲線出現(xiàn)了一定程度的差別.

圖6 駐點熱流系數(shù)計算結(jié)果Fig.6 Heat flux coefficients on stagnation point

圖7 歸一化的熱流和摩阻分布Fig.7 Normalized heat fluxes and skin frictions

歸一化的熱流和摩阻分布曲線預(yù)示不同雷諾數(shù)條件下的雷諾比擬關(guān)系規(guī)律近似一致.由圖8 可見，不同來流壓力p，也就是不同雷諾數(shù)條件下，廣義雷諾比擬關(guān)系具有如下規(guī)律：

(1)在θ ≤60?時，不同雷諾數(shù)條件下的Cf/Ch曲線分布一致，數(shù)值計算的雷諾比擬關(guān)系與理論吻合較好.在θ=60?時，數(shù)值計算與理論預(yù)示結(jié)果的相對誤差為2.1%～6.5%；

(2)在θ>60?時，Cf/Ch隨雷諾數(shù)增大而增大，且偏離理論預(yù)示的線性關(guān)系.在θ=90?時，數(shù)值計算與理論預(yù)示結(jié)果的相對誤差最大為31.6%.

圖8 不同來流條件下的廣義雷諾比擬關(guān)系Fig 8 General Reynolds analogy relation with p=4.83～4830

針對不同雷諾數(shù)條件的計算結(jié)果表明，對于圓柱外形，即使是較大雷諾數(shù)條件，理論推導(dǎo)的線性廣義雷諾比擬關(guān)系仍然可以較好地預(yù)示Cf/Ch分布.但是在比較靠近下游的區(qū)域，數(shù)值計算得到的廣義雷諾比擬關(guān)系則與理論預(yù)示存在差異，且二者之間的差異隨雷諾數(shù)增大而加大.

對于圓柱下游區(qū)域數(shù)值計算的雷諾比擬關(guān)系與線性理論預(yù)示結(jié)果的區(qū)別，可能是由于在靠近下游的部位，邊界層的自相似程度發(fā)生變化，導(dǎo)致線性的廣義雷諾比擬關(guān)系理論基礎(chǔ)與實際情況有所區(qū)別.另外，圓柱壁面的熱流分布隨θ 快速降低，在靠近下游的區(qū)域，數(shù)值計算得到的熱流相對誤差可能隨之加大，導(dǎo)致雷諾比擬關(guān)系的計算結(jié)果與實際情況出現(xiàn)差別.

2 冪次體壁面的廣義雷諾比擬關(guān)系

理論上，廣義雷諾比擬關(guān)系適用于一般的鈍頭體外形.除圓柱以外，本文另外選取一種冪次外形鈍頭體，采用數(shù)值求解N-S 的方法計算了不同來流條件下的壁面雷諾比擬關(guān)系.冪次體曲線方程為y=0.16x0.5,x=0～0.23 m.計算的來流條件與表1 一致.

冪次體外形與圓柱外形計算得到的流場馬赫數(shù)云圖如圖9 所示.冪次體壁面附近的流動圖像仍然是由脫體弓形激波主導(dǎo)，但由于外形相對更加尖銳，激波脫體距離比圓柱外形明顯減小.

圖9 圓柱和冪次體繞流馬赫數(shù)云圖Fig.9 Mach number counters on flows around a cylinder and power-law body

圖10 冪次體壁面的雷諾比擬關(guān)系Fig.10 Reynolds analogy relations on a power-law body

不同來流馬赫數(shù)和壓力條件下計算得到的壁面雷諾比擬系數(shù)見圖10.結(jié)果表明，在不同的雷諾數(shù)條件下冪次體壁面的雷諾比擬系數(shù)基本一致，且相比于圓柱外形，比擬系數(shù)更接近線性分布.另一方面，在冪次體壁面，雷諾比擬系數(shù)Cf/Ch沿θ 分布的斜率小于圓柱外形上分布.在=6,Tw=500 K 的條件下，式(2)預(yù)示的圓柱壁面的雷諾比擬關(guān)系為=1.942，而冪次體外形的計算結(jié)果則近似為≈1.48；=12,Tw=1000 K 條件下，圓柱壁面的結(jié)果為=1.522，而冪次體外形的計算結(jié)果近似為≈1.24.在θ=80?時，數(shù)值計算與理論預(yù)示結(jié)果的相對誤差最大為10.2%.

利用CFD 方法計算鈍體壁面壓力分布，根據(jù)式(1)分別計算得到的圓柱體和冪次體ue速度分布見圖11.結(jié)果表明不同外形邊界層外緣的速度分布存在一定差異，可能是導(dǎo)致不同鈍頭外形下雷諾比擬系數(shù)分布規(guī)律有所區(qū)別的原因.具體的影響機(jī)理還有待后續(xù)進(jìn)一步的理論研究.

圖11 邊界層外緣速度分布Fig.11 Velocities along boundary layer edge

3 結(jié)論

本文采用數(shù)值求解N-S 方程的方法研究了鈍頭體壁面的廣義雷諾比擬關(guān)系.采用不同的數(shù)值方法計算了不同來流條件和鈍頭外形下的壁面的熱流、摩阻以及雷諾比擬系數(shù)的分布.數(shù)值計算結(jié)果與理論預(yù)示以及文獻(xiàn)中公布的數(shù)據(jù)吻合.

通過改變來流壓力條件，計算了不同雷諾數(shù)條件下圓柱壁面的熱流和摩阻分布，數(shù)值計算得到的駐點熱流與經(jīng)典理論預(yù)示結(jié)果一致.不同來流條件下的歸一化熱流分布和摩阻分布在距離駐點下游較遠(yuǎn)的位置隨雷諾數(shù)增加存在一定差異.受其影響，在不同的雷諾數(shù)條件下，壁面雷諾比擬關(guān)系在較下游的位置也存在一定區(qū)別，較高雷諾數(shù)下計算得到的Cf/Ch隨θ 分布曲線偏離理論預(yù)示的線性分布規(guī)律.

本文同時計算了拋物線形鈍頭體壁面的雷諾比擬關(guān)系.在相同的來流條件下，相比于圓柱外形，冪次體壁面的雷諾比擬關(guān)系與線性規(guī)律吻合較好，且在不同雷諾數(shù)條件下的區(qū)別較小.受到鈍頭體外形對邊界層外緣以及邊界層內(nèi)速度和溫度的影響，不同外形壁面的雷諾比擬系數(shù)存在一定區(qū)別.

本文的數(shù)值計算和對比研究結(jié)果表明，理論預(yù)示的廣義雷諾比擬關(guān)系存在于較大范圍的雷諾數(shù)區(qū)間以及不同類型的鈍頭體外形中.在實際工程應(yīng)用中，如果針對具體的流動條件和鈍頭外形進(jìn)行修正，則可以進(jìn)一步提高預(yù)示精度.