李玉川，鮑 麟

(中國科學(xué)院大學(xué)工程科學(xué)學(xué)院， 北京 100049)(2018年5月9日收稿； 2018年5月28日收修改稿)

非光滑表面的高速流動傳熱規(guī)律既是一個(gè)基礎(chǔ)問題，也是航天工程領(lǐng)域廣泛關(guān)注和研究的課題。以再入式航天飛船為例，其再入行星大氣層時(shí)速度可高達(dá)20～30倍聲速，相應(yīng)產(chǎn)生高達(dá)數(shù)千K的高溫，為使飛行器免受高溫?fù)p壞，需采用燒蝕熱防護(hù)技術(shù)，即在其表面涂上反應(yīng)潛熱高的熱防護(hù)材料，依靠高溫?zé)g過程，帶走大量熱流，以保護(hù)飛行器。由于典型材料燒蝕過程是非均勻的，會形成粗糙壁面，研究其氣動加熱規(guī)律成為航天工程中的重要問題之一。從20世紀(jì)70年代以來，以美國PANT計(jì)劃[1]為代表，針對粗糙壁面熱流開展了大量研究，一是利用實(shí)驗(yàn)數(shù)據(jù)擬合得到若干經(jīng)驗(yàn)公式計(jì)算粗糙壁面的熱流，如Dirling[2]，Phinney[3]，Grabow[4]，二是用各種數(shù)值方法計(jì)算流場和粗糙壁面的熱流，如Yang等[5]，李俊紅等[6]。這些工作共同揭示影響粗糙壁面熱流主要參數(shù)有雷諾數(shù)、邊界層厚度以及等效粗糙元高度等，為工程設(shè)計(jì)提供了重要理論基礎(chǔ)。

由于前人的粗糙壁研究源于傳統(tǒng)的大鈍頭式再入飛行器，其波阻遠(yuǎn)大于摩擦阻力，因此這些研究大多只關(guān)心粗糙壁的熱流，而極少有對摩阻的討論。當(dāng)前正在研制的下一代航天飛行器將具備大氣層內(nèi)高速巡航能力，多采用扁平化的乘波體等氣動外形，其波阻大大減小，而表面摩阻所占比重極大地上升，最多可達(dá)30%以上[7]，因此，飛行器壁面流動的摩擦阻力研究越來越受到重視。另一方面，下一代飛行器將實(shí)現(xiàn)可重復(fù)使用功能，工程設(shè)計(jì)提出采用非燒蝕熱防護(hù)系統(tǒng)，使飛行器氣動外形保持不變。實(shí)際上，飛行器表面在極高溫作用下，不可能完全不燒蝕，尤其是多次長時(shí)間飛行后，其壁面應(yīng)當(dāng)看作存在微燒蝕，即是微粗糙的。綜合上述兩方面，微粗糙表面在高速流動中的壁面摩擦和傳熱是當(dāng)前值得研究的問題，對新型飛行器的氣動布局以及熱防護(hù)系統(tǒng)設(shè)計(jì)有應(yīng)用價(jià)值。

回顧光滑平板層流邊界層理論可知，壁面熱流與摩阻由臨近壁面的邊界層特性決定，且它們之間存在雷諾比擬關(guān)系[8]，即壁面熱流與摩阻是成正比的，沿著平板流動方向，其比值是一個(gè)常數(shù)，僅與來流馬赫數(shù)Ma∞、壁面溫度Tw、氣體的Pr數(shù)和比熱比γ有關(guān)。對微粗糙壁而言，高速氣體流經(jīng)非光滑表面時(shí)，不僅壁面摩阻與熱流發(fā)生變化，此時(shí)雷諾比擬是否仍成立存在疑問，而且與一般低速流動不同，壓強(qiáng)波動產(chǎn)生的波阻是超聲速時(shí)特有的現(xiàn)象，熱流又與波阻有什么關(guān)系，是需要討論的問題。

本文基于OpenFOAM程序數(shù)值模擬波狀粗糙壁面在高速來流中壁面阻力與熱流，及其流場特征，結(jié)果發(fā)現(xiàn)簡單雷諾比擬關(guān)系不適用于粗糙壁情況，應(yīng)代之以壁面熱流與阻力之間的比擬關(guān)系，并分析該比擬關(guān)系的主要影響參數(shù)。為了深入分析該問題中的基本規(guī)律，同時(shí)針對新型飛行器高空巡航的特點(diǎn)，本文計(jì)算僅考慮層流情況，流動介質(zhì)是空氣，且看作量熱完全氣體(即特定氣體常數(shù)R=287 J/(kg·K)、比熱比γ=1.4、Pr=0.72，以下不再說明)，但結(jié)果揭示出的流動規(guī)律和機(jī)理也可為進(jìn)一步研究各類工程實(shí)際情況提供參考。

1 物理模型與數(shù)值方法

在實(shí)際問題中，粗糙壁面形狀是極其復(fù)雜的，即使是微粗糙的壁面，也不易用較為簡單的形狀來描述。早期實(shí)驗(yàn)中為了對其進(jìn)行刻畫，通常采用在光滑表面上均布不同形狀和尺寸的粗糙元來模擬粗糙壁面，而為了使這些不同形狀大小和排布的粗糙壁面具有相同的度量特征，Nikuradse[9]提出等效砂礫粗糙度：排布著普通粗糙元的壁面熱流，和緊密排布著高度為ks的半球(三維)或半圓(二維)粗糙元的表面的熱流相同，那么ks即為此粗糙壁面的等砂礫粗糙度。Signal和Danberg[10]對普通粗糙元到等效砂礫粗糙度的映射關(guān)系做了改進(jìn)，在工程中得到廣泛應(yīng)用。

等效砂礫粗糙度概念的提出啟示我們，可以利用一種標(biāo)準(zhǔn)粗糙元模型等效地模擬實(shí)際粗糙壁面形狀，以獲得壁面熱流變化的規(guī)律，那么不失一般性，本文選取正弦波狀外形作為標(biāo)準(zhǔn)粗糙元模型。

圖1 正弦波狀粗糙壁面示意圖Fig.1 Schematic diagram of the sinusoidal rough wall

本文將利用數(shù)值計(jì)算方法求解上述粗糙壁模型的流動問題，主要探討壁面熱流、摩阻與壓阻之間的關(guān)系，為此，引入定義無量綱壁面壓強(qiáng)系數(shù)Cp、摩阻系數(shù)Cf與熱流率系數(shù)Ch[8]：

(1)

已知高速氣流的波狀板流動的控制方程為二維定常、可壓縮、黏性Navier-Stokes方程(以下簡稱N-S方程)：

(2)

有限體積法(finite volume methods，簡稱FVM)是計(jì)算流體力學(xué)中常用的數(shù)值方法，能夠較好地保證流動的物理守恒性，因此，通常認(rèn)為FVM類方法計(jì)算的壁面熱流或摩阻的精度較高，本文利用OpenFOAM開源FVM程序，選取二階Kurganov-Tadmor中心格式[11]求解上述N-S方程。具體格式詳見文獻(xiàn)[11]，在此不再贅述。重要的是該格式滿足TVD特性和最大值原理，同時(shí)模擬流場具有高解析度，且與時(shí)間步長大小無關(guān)，因此，是一種理想的高速流動的計(jì)算格式。此程序的具體算例驗(yàn)證將在下一節(jié)中給出。

2 數(shù)值程序的驗(yàn)證

高速粗糙壁流動中壓強(qiáng)梯度較大，這要求計(jì)算程序具有更高的流場解析度，為了驗(yàn)證本文采用的OpenFOAM程序的可靠性，首先，計(jì)算高馬赫數(shù)光滑平板邊界層前緣的黏性干擾區(qū)，并將壁面壓強(qiáng)、摩阻與熱流計(jì)算結(jié)果與已有的黏性干擾理論值作對比；其次，模擬超聲速正弦波狀板流動，在距離前緣充分遠(yuǎn)的區(qū)域，檢驗(yàn)流場的數(shù)值解是否與理論上的二維定常超聲速波形壁理論解一致；最后，對正弦波狀板的高速流動計(jì)算結(jié)果的時(shí)間步長和網(wǎng)格密度無關(guān)性做檢驗(yàn)。

首先，黏性干擾區(qū)算例的參數(shù)選取為Ma∞=8，單位雷諾數(shù)ReL=2×105m-1。依據(jù)Hayes和Probstein[12]得出的強(qiáng)干擾區(qū)摩阻和熱流的分布公式為：

(3)

圖2 平板黏性干擾區(qū)的Cf、Cp、Ch計(jì)算值與理論值的對比Fig.2 Comparison between the calculated and theoretical values of Cf, Cp， and Chalong a flat plate with viscous interference

圖3 無窮長波狀板二維超聲速流動的壓強(qiáng)分布的理論值與計(jì)算值對比Fig.3 Comparison between the calculated and theoretical values of the pressure distribution along a wavy plate in 2-D supersonic flow

圖4 不同計(jì)算參數(shù)下的波狀板計(jì)算Ch、CQ、Cp曲線Fig.4 The calculated curves of Ch, CQ, and Cpalong the wavy plate with different calculation parameters

3 算例及分析

3.1 沿波狀壁面的熱流-阻力關(guān)聯(lián)關(guān)系

(4)

(5)

CD=CDf+CDp，

(6)

(7)

圖5 h較小時(shí)CQ/CDf和CQ/CD隨的分布(Ma∞=4，Reλ=2×104)Fig.5 Distributions of CQ/CDfand CQ/CDalong at the small hvalues (Ma∞=4，Reλ=2×104)

圖6 不同h下的CQ/CD隨的變化曲線(Ma∞=4，Reλ=2×104)Fig.6 Distribution of CQ/CDalong at various hvalues(Ma∞=4，Reλ=2×104)

圖7 h=0.03和h=0.06時(shí)壁面附近的流場圖Fig.7 The flow field diagram near the wall when h=0.03 and 0.06

(8)

式中：Reθtr是以動量邊界層厚度θ作為特征尺度的雷諾數(shù)；k為粗糙元尺度，即為本文中2ε。定性上講，波狀板的粗糙元尺度越大，位置離前緣越遠(yuǎn)，流動越趨向轉(zhuǎn)捩乃至湍流。按本文設(shè)定的來流條件，設(shè)最大板長為x=1 m，即Rex約為105～106，則動量邊界層厚度θ量級約為10-3m[16]，代入式(8)得到臨界轉(zhuǎn)捩粗糙元高度約為k≈0.01 m，即h=0.1。因此，本文僅討論h<0.1時(shí)的波狀板流動，將其看作層流流動，并且不包含出現(xiàn)局部紊亂流場時(shí)的情況，下面將集中探討定常層流流動條件下的CQ/CD的具體變化規(guī)律。

3.2 熱流-阻力關(guān)聯(lián)關(guān)系的歸一化

3.2.1 未分離的波狀板流動

(9)

(10)

圖8 不同Ma∞和Reλ下波狀板附著流動CQ/CD隨的變化規(guī)律Fig.8 The law of the CQ/CDchange with at different Ma∞and Reλvalues along the wavy plate in attached flow

式(10)由數(shù)值計(jì)算結(jié)果擬合得出，并無可對比的實(shí)驗(yàn)結(jié)果。在此僅從數(shù)值驗(yàn)證角度估算其誤差范圍：由第2節(jié)，壁面熱流系數(shù)和摩阻系數(shù)最大偏差約為5%，而壓強(qiáng)系數(shù)的最大偏差為10%，因此，可大致估計(jì)CQ/CD最大偏差約為15%～20%，可以供工程預(yù)研或工程估算時(shí)參考。

3.2.2 存在分離的波狀板流動

(11)

(12)

圖9 不同Ma∞和Reλ下波狀板分離流動的CQ/CD隨等效位置的變化規(guī)律Fig.9 The law of the CQ/CDchange with the equivalent position at different Ma∞and Reλvalues along the wavy plate in separation flow

需要說明，雖然流動出現(xiàn)局部分離渦后，波狀板上CQ/CD具有相似性規(guī)律，但流動在何處產(chǎn)生分離仍是一個(gè)難題。定性上講，波狀板來流Ma∞越低，Reλ越高，流動越容易分離，并且波狀粗糙度越大，位置離前緣越遠(yuǎn)，流動越趨向于分離。然而暫時(shí)無法定量地給出流動開始分離的位置，這仍有待深入研究。

4 結(jié)論

CQ/CD=

無論從前人歷史研究的結(jié)論看，還是從波狀粗糙板流動本身的分析角度看，定量確定流動開始分離的位置仍然是一個(gè)非常困難的問題，本文僅討論壁面熱流-阻力的關(guān)系，分離條件仍有待下一步深入研究。

本文得到王智慧副教授的諸多指導(dǎo)和啟發(fā)，在具體工作中也得到博士生李賢冬、羅健的熱心幫助，在此謹(jǐn)表感謝。