鈍頭體中的廣義雷諾比擬關系1)

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

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

引言

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

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

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

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

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

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

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

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

1 圓柱壁面的廣義雷諾比擬關系

1.1 理論預示

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

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

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

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

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

其中，Cr為一個與壁面溫度相關的比例系數，根據前期研究結果,對于二維圓柱外形

對于軸對稱圓球外形

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

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

1.2 N-S 方程數值計算

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

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

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

圖2 計算域與網格示意圖Fig.2 Computation regime and mesh

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

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

圖3 不同馬赫數條件下的雷諾比擬關系Fig.3 General Reynolds analogy relation in different Mach number

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

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

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

圖5 不同網格數目下的雷諾比擬系數計算結果Fig.5 Reynolds analogy coefficients under different amount of mesh

1.3 雷諾數對廣義雷諾比擬關系的影響

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

表1 計算條件Table 1 Simulation cases

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

圖6 駐點熱流系數計算結果Fig.6 Heat flux coefficients on stagnation point

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

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

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

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

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

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

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

2 冪次體壁面的廣義雷諾比擬關系

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

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

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

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

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

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

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

3 結論

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

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

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

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