孔祥志，韓宇峰*

天津大學(xué) 高速空氣動(dòng)力學(xué)研究室，天津 300072

一直以來，氣動(dòng)加熱在高超聲速飛行器的設(shè)計(jì)中占據(jù)著關(guān)鍵的地位，是工程實(shí)際及理論研究都關(guān)注的問題。飛行器表面的熱防護(hù)材料在經(jīng)過燒蝕后，會(huì)發(fā)生一定程度的變形，形成粗糙表面，對(duì)流動(dòng)和氣動(dòng)加熱產(chǎn)生不可忽略的影響。

在實(shí)際問題中，粗糙壁面形式多樣，開展對(duì)粗糙壁面峰值熱流特性的研究，首先要建立描述粗 糙 壁 面 特 征 的 模 型。Colebrook 和White［1］的研究發(fā)現(xiàn)，粗糙密度對(duì)于表征粗糙特征是至關(guān)重要的，后來Schlichting［2］量化表征了粗糙密度，用λ=Fr/F 定義粗糙密度，F(xiàn)r表示粗糙壁正面總投影面積，F(xiàn) 表示粗糙與壁面平行的投影面積；Nikuradse［3］基于粗糙密度提出等效砂礫高度的概念，通過試驗(yàn)建立起等效沙礫高度與粗糙壁面特征的關(guān)系；Sigal 和Danberg［4］擴(kuò)展 并改進(jìn)了這一映射關(guān)系，使用Sigal-Danberg 參數(shù)定義了粗糙密度，使等效砂礫高度獲得更為廣泛的應(yīng)用。此外，為了表征粗糙壁面特征，Napoli 等［5］提出有效斜率的概念衡量粗糙壁面的高度及密度，用于對(duì)完全粗糙壁面特征的表征，同時(shí)研究了有效斜率與粗糙壁面摩阻和熱流的關(guān)系。上述粗糙特征模型由于提出角度的不同，適用于不同壁面粗糙情況［6］。

對(duì)于高超聲速粗糙壁的氣動(dòng)熱規(guī)律，從20 世紀(jì)70 年代以來開展了大量的實(shí)驗(yàn)研究。Bertran等［7］在Jaeck［8］的試驗(yàn)報(bào)告基礎(chǔ)上，研究平板波紋型粗糙壁的傳熱特性，重點(diǎn)研究第1 波峰值熱流的演變規(guī)律，發(fā)現(xiàn)了粗糙高度、寬度、間隔等對(duì)于熱流的影響規(guī)律；美國(guó)PANT 計(jì)劃、文獻(xiàn)［8-14］利用實(shí)驗(yàn)得到的數(shù)據(jù)，擬合得到若干計(jì)算粗糙壁面熱流的經(jīng)驗(yàn)預(yù)測(cè)公式，受限于實(shí)驗(yàn)數(shù)據(jù)以及實(shí)驗(yàn)條件，這些公式不具有很大的通用性；Nestler［15］綜述了表面粗糙對(duì)傳熱影響的實(shí)驗(yàn)數(shù)據(jù)以及預(yù)測(cè)方法，這些表面粗糙的類型包括臺(tái)階、空穴、間隙和突起等，通過對(duì)該領(lǐng)域關(guān)鍵技術(shù)的評(píng)估，指出以往對(duì)于粗糙壁氣動(dòng)熱特性的研究缺乏對(duì)傳熱過程的物理認(rèn)知，建立的經(jīng)驗(yàn)相關(guān)性公式是不充分的，需要更完善的實(shí)驗(yàn)數(shù)據(jù)以及新的預(yù)測(cè)方法。

除了實(shí)驗(yàn)，學(xué)者們也采用各種數(shù)值方法對(duì)高超聲速粗糙壁的熱流進(jìn)行了研究。Polak 等［16］使用隱式有限差分技術(shù)研究不同馬赫數(shù)、雷諾數(shù)下的波紋型粗糙壁面的熱流演變規(guī)律，并對(duì)Jaeck等［8］提出的峰值熱流預(yù)測(cè)公式進(jìn)行評(píng)估，發(fā)現(xiàn)預(yù)測(cè)公式對(duì)試驗(yàn)結(jié)果存在嚴(yán)重依賴；Egorov 等［17］數(shù)值研究波紋粗糙平板壁面，得出壁面熱流的變化規(guī)律，指出在粗糙壁面峰值附近，熱流會(huì)達(dá)到最大；李俊紅等［18］用數(shù)值模擬和理論方法對(duì)高超聲速可壓縮湍流中的粗糙壁面熱增量進(jìn)行研究，討論了不同粗糙高度、密度等對(duì)熱流的影響規(guī)律，比較了幾種預(yù)測(cè)公式［12-14］的結(jié)果與數(shù)值結(jié)果的差異；李玉川和鮑麟［19］研究波紋粗糙壁面上高速氣流的流動(dòng)和傳熱特性，結(jié)果表明在波狀壁流動(dòng)中，壁面熱流系數(shù)與總阻力系數(shù)呈現(xiàn)出直接關(guān)聯(lián)。Kumar 和Reddy［20］在 平 板上布 置 三 維 凸起，通過實(shí)驗(yàn)研究凸起附近的熱流，使用不同工況下的實(shí)驗(yàn)數(shù)據(jù)擬合得到用于預(yù)測(cè)凸起附近熱流的關(guān)系式；Estruch-Samper［21］在平板上布置短壓縮斜坡粗糙元，研究粗糙偏轉(zhuǎn)角對(duì)于峰值熱流等的影響，評(píng)估幾種峰值熱流預(yù)測(cè)公式的有限適用性；Avallone 等［22］通過實(shí)驗(yàn)研究圓柱形粗糙元的流動(dòng)以及熱流峰值特征；Thakkar 等［23］使用多組不規(guī)則粗糙表面的模擬數(shù)據(jù)，開發(fā)了具有多個(gè)參數(shù)的粗糙壁面特征優(yōu)化模型，但是這一模型的參數(shù)較多，實(shí)際使用起來較為復(fù)雜。

總之，目前對(duì)于粗糙壁面的氣動(dòng)熱有了深入的認(rèn)知，但多數(shù)研究集中在再入飛行背景下的大尺度（毫米-厘米量級(jí)）燒蝕。隨著新型航天航空飛行器的發(fā)展以及其表面精細(xì)化的熱防護(hù)需求，復(fù)合燒蝕材料由于輕量以及高密度等優(yōu)勢(shì)，應(yīng)用越來越廣泛，相比于傳統(tǒng)燒蝕材料，復(fù)合燒蝕材料的線燒蝕率更低。例如在10～12 馬赫數(shù)范圍的高超聲速巡航或滑翔飛行器上，機(jī)翼結(jié)構(gòu)采用瓦狀非承重外部隔熱系統(tǒng)，對(duì)于低重量的隔熱瓦，會(huì)產(chǎn)生百微米量級(jí)的波紋型表面粗糙［24］。毫米以下的微粗糙和毫米、厘米量級(jí)的大粗糙的流場(chǎng)結(jié)構(gòu)是有顯著差異的。微粗糙尺度遠(yuǎn)小于邊界層厚度的，流場(chǎng)結(jié)構(gòu)較為簡(jiǎn)單，壁面熱流增量主要是由于壁面剪切耗散增大而導(dǎo)致［16］；大粗糙尺度與邊界層厚度相當(dāng)，流場(chǎng)結(jié)構(gòu)更為復(fù)雜，粗糙附近壓力波動(dòng)較大并出現(xiàn)流動(dòng)分離，熱流會(huì)受分離渦影響，熱增量機(jī)理更為復(fù)雜。目前已有的預(yù)測(cè)方法，都是針對(duì)大粗糙壁的基于實(shí)驗(yàn)擬合得到的半經(jīng)驗(yàn)平均熱流預(yù)測(cè)公式。對(duì)于微粗糙度的峰值熱流，這些方法是不適用的。此外在這種微燒蝕情況下的粗糙壁面峰值熱流的預(yù)測(cè)是很稀缺的［16］。

在實(shí)際中，飛行器表面因燒蝕而產(chǎn)生的粗糙形狀是復(fù)雜的，但是等效砂礫粗糙度概念的提出啟示我們可以利用一種標(biāo)準(zhǔn)粗糙元模型等效地模擬實(shí)際粗糙壁面形狀以獲得壁面熱流變化的規(guī)律。

本文重點(diǎn)關(guān)注微燒蝕條件下出現(xiàn)的百微米量級(jí)的波紋型粗糙壁，結(jié)合目前對(duì)于粗糙壁熱流影響因素的認(rèn)識(shí)以及預(yù)測(cè)方法，使用直接數(shù)值模擬計(jì)算不同高度、流動(dòng)參數(shù)下的流場(chǎng)及熱流。在此基礎(chǔ)上，進(jìn)行理論分析，總結(jié)數(shù)值結(jié)果的規(guī)律，提出定量預(yù)測(cè)粗糙壁面峰值熱流的模型。

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

1. 1 控制方程

通過直接數(shù)值求解直角坐標(biāo)系下的可壓縮守恒型Navier-Stokes（N-S）方程獲得平板邊界層的流場(chǎng)，其控制方程為

式中：t 為時(shí)間變量；x、y 分別為流向和法向的空間坐標(biāo)；U 為守恒型變量；E、F 分別為流向、法向的通量；Ev、Fv為相應(yīng)的黏性項(xiàng)，各項(xiàng)的具體表達(dá)式分別為

其中：u、v 分別為流向、法向的速度分量；p、ρ 和T分別表示壓力、密度和溫度；μ 為黏性系數(shù)；κ 為熱傳導(dǎo)系數(shù)；黏性系數(shù)使用Sutherland 公式計(jì)算：

式中：Tμ為參考溫度，取110.4 K。

熱傳導(dǎo)系數(shù)在普朗特?cái)?shù)以及定壓比熱為常數(shù)的條件下與黏性系數(shù)為線性關(guān)系：

式中：Tκ為參考溫度，取110.4 K。

在計(jì)算過程中需要選擇適當(dāng)?shù)奶卣鞒叨葘?duì)方程進(jìn)行無量綱化。上標(biāo)*號(hào)表示有量綱物理量，無上標(biāo)*號(hào)的量表示無量綱物理量。選取特征長(zhǎng)度L*對(duì)空間坐標(biāo)進(jìn)行無量綱化：

雷諾數(shù)、馬赫數(shù)分別為

無量綱形式的黏性系數(shù)和熱傳導(dǎo)系數(shù)分別為

1. 2 數(shù)值格式和程序驗(yàn)證

壁面采用無滑移邊界條件，上邊界給定來流條件，入口采用自由來流條件，出口采用線性外推。采用基于有限差分的直接數(shù)值模擬方法開展流場(chǎng)計(jì)算，計(jì)算對(duì)流項(xiàng)采用三階WENO（Weighted Essentially Non-Oscillatory）格式，黏性項(xiàng)采用四階中心格式，通量分裂選擇Steger-Warming 分裂，采用四階龍格-庫塔方法在時(shí)間方向推進(jìn)求解。

為了驗(yàn)證數(shù)值程序的正確性，與Zhao 等［25］關(guān)于單個(gè)正弦粗糙元的高超聲速流動(dòng)研究中的工況進(jìn)行比對(duì)，圖1 展示了流場(chǎng)的壁面壓強(qiáng)分布和粗糙元附近的速度剖面的對(duì)比結(jié)果，可以看到，在粗糙元附近，計(jì)算結(jié)果基本吻合，驗(yàn)證了所使用程序在計(jì)算粗糙元時(shí)的正確性。

圖1 程序驗(yàn)證Fig.1 Program verification

1. 3 物理模型、網(wǎng)格劃分策略及計(jì)算參數(shù)

數(shù)值模擬的物理模型為全場(chǎng)均布正弦型波紋的平板，如圖2（a）所示，壁面連續(xù)分布有100 個(gè)周期的正弦凸起。壁面函數(shù)表達(dá)式為

圖2 正弦型壁面粗糙分布及熱流驗(yàn)證Fig.2 Sinusoidal wall roughness distribution and heat flux verification

式中：h0、kn和L 分別為粗糙元的高度、個(gè)數(shù)和寬度。近壁面使用貼體網(wǎng)格。同時(shí)為了使粗糙元以外的部分不受到網(wǎng)格彎曲帶來的影響，在距離粗糙元一定高度以上，仍采用平行網(wǎng)格，二者光滑過渡；為了準(zhǔn)確計(jì)算熱流，法向網(wǎng)格尺度使用基于分子平均自由程的網(wǎng)格準(zhǔn)則［26-27］得到，法向第1 層網(wǎng)格尺度在10-6m。在每個(gè)粗糙單元尺度上流向分布20 多個(gè)點(diǎn)，法向分布100 多個(gè)點(diǎn)，計(jì)算域總網(wǎng)格量為100 萬左右，具體網(wǎng)格分布如圖2（a）所示，圖中右上角為x、y 坐標(biāo)比例為1 的示意圖。不同法向第1 層長(zhǎng)度（10-7～10-5m）下熱流的計(jì)算結(jié)果見圖2（b），圖中y1表示法向第1 層網(wǎng)格長(zhǎng)度，可以看出，10-5m 的結(jié)果偏差很大，5×10-6m 的結(jié)果趨于收斂，但仍有細(xì)微的偏差。而在長(zhǎng)度為10-6、5×10-7、10-7m 時(shí)，熱流結(jié)果完全一致，說明在法向第1 層網(wǎng)格尺度10-6m下可保證熱流的計(jì)算是準(zhǔn)確的。后文計(jì)算中每個(gè)算例都經(jīng)過了網(wǎng)格無關(guān)性檢驗(yàn)，以保證計(jì)算結(jié)果的準(zhǔn)確性。

使用無量綱的熱流系數(shù)ch 來表征壁面熱流：

式中：chsm表示對(duì)應(yīng)光滑壁的熱流系數(shù)；計(jì)算的來流參數(shù)選自30 km 的大氣，特征長(zhǎng)度選取1 m，具體參數(shù)見表1。

表1 計(jì)算參數(shù)Table 1 Calculation parameters

2 數(shù)值結(jié)果及分析

粗糙壁的幾何特性會(huì)影響壁面熱流的演變規(guī)律，Sigal［4］、Bertram［7］、Polak［16］、李俊紅［18］等對(duì)粗糙壁高度、寬度、形狀等幾何特征對(duì)壁面熱流的影響進(jìn)行了研究，結(jié)果表明，高度和寬度是較為關(guān)鍵的影響因素，兩者通常以比值形式的關(guān)系［4］體現(xiàn)這種影響規(guī)律。

為了研究不同幾何特征以及來流參數(shù)對(duì)于粗糙壁峰值熱流的影響規(guī)律，數(shù)值計(jì)算了以下4種工況的流場(chǎng)和熱流，分別是：

1） 來流參數(shù)不變，粗糙壁寬度不變，不同高度下的粗糙壁。

2） 來流參數(shù)不變，粗糙壁高度不變，不同寬度下的粗糙壁。

3） 來流參數(shù)不變，高度寬度均變化，固定寬高比的粗糙壁。

4） 粗糙壁高度和寬度不變，不同馬赫數(shù)以及壁溫下的粗糙壁。

根據(jù)以上4 種情況，在圖3 展示了具體計(jì)算工況的高度和寬度分布，圖3（a）列出了前3 種情況的粗糙壁高度和寬度分布，圖3（b）列出了第4 種情況下的工況分布。共計(jì)算了28 個(gè)不同工況下的粗糙壁流場(chǎng)和熱流，使用的網(wǎng)格流向與法向均已無關(guān)。同時(shí)，由于平板頭部存在黏性干擾，主要分析0.2～1 m 的流場(chǎng)。在圖3（c）與圖3（d）中給出了高度0.421 2 mm、馬赫數(shù)為6 時(shí)的計(jì)算結(jié)果的整體以及局部區(qū)域的速度云圖。

圖3 計(jì)算工況以及典型計(jì)算結(jié)果Fig.3 Calculation conditions and typical calculation results

2. 1 不同高度下的壁面熱流

首先分析第1 類工況下的壁面熱流，即來流參數(shù)不變，粗糙壁寬度不變，計(jì)算不同高度下粗糙壁的基本流場(chǎng)。具體的粗糙元配置：寬度均為10.53 mm，高 度 分 別 為0、0.105 3、0.126 36、0.210 6、0.294 84、0.379 08、0.421 2、0.526 5 mm。

圖4 給出了不同高度下，全場(chǎng)均布粗糙的壁面熱流和光滑平板的對(duì)比，同時(shí)將壁面形狀給出，方便與壁面熱流做比對(duì)。從圖4（a）可以發(fā)現(xiàn)，粗糙壁面的熱流隨壁面形式呈周期性變化，在正弦壁峰附近，熱流達(dá)到峰值；同一流向位置的壁面峰值熱流，隨著高度的增加而增大，這一熱增量不可忽略。由于粗糙元高度相比邊界層厚度很小，因此熱流的變化規(guī)律和壁面幾何形狀類似。圖4（b）給出了局部放大的結(jié)果。圖4（c）給出了較大高度（0.421 2 mm）下平均熱流與光滑壁熱流的比較，粗糙壁平均熱流是對(duì)波紋粗糙元一個(gè)周期內(nèi)的熱流進(jìn)行平均求得的，而光滑壁面的平均熱流指相同流向區(qū)域的熱流平均值。可見在微尺度的粗糙下，壁面平均熱流增量是很小的，即使在較大高度下，微粗糙下的平均熱流增量與光滑壁平均熱流的比值小于3%，因此本文重點(diǎn)關(guān)注峰值熱流。

圖4 不同高度下的壁面熱流Fig.4 Surface heat flux at different heights

為了分析熱流隨著高度的變化規(guī)律，將不同流向位置的壁面峰值熱流提取出來進(jìn)行分析比較，對(duì)于每一個(gè)工況可以得到100 個(gè)峰值點(diǎn)，不同高度下的壁面峰值熱流系數(shù)及相對(duì)熱流增量Qratio如圖5 所示。圖5（a）給出了不同高度下的壁面峰值熱流系數(shù)，可以發(fā)現(xiàn)隨著流向位置x 的增大，壁面峰值熱流減??；圖5（b）則給出了不同高度下的相對(duì)峰值熱流增量，很明顯，沿流向的相對(duì)熱流增量近似保持不變。同時(shí)，在同一流向位置，隨著高度增加，壁面峰值熱流和相對(duì)峰值熱流增量均增大。壁面峰值熱流沿流向位置的趨勢(shì)與反比例函數(shù)的形式近似，隨高度增大的趨勢(shì)與線性函數(shù)近似，相對(duì)峰值熱流增量與流向位置無關(guān)。為了確定壁面峰值熱流與流向位置的具體關(guān)系，圖6 給出了不同高度下壁面峰值熱流與x-0.5的關(guān)系。從圖中可以看出ch與x-0.5成線性相關(guān)關(guān)系，這一關(guān)系與光滑平板壁面熱流和x-0.5的線性規(guī)律一致。

圖5 不同高度下壁面峰值熱流及相對(duì)峰值熱流增量Fig.5 Surface peak heat flux and relative surface peak heat flux increment at different heights

圖6 不同高度下壁面峰值熱流與x-0.5 的關(guān)系Fig.6 Relationship between surface peak heat flux and x-0.5 at different heights

圖7給出了3 個(gè)不同流向位置下，壁面峰值熱流隨壁面高度變化的關(guān)系，很明顯，ch與h 具有較強(qiáng)的線性相關(guān)性。因此，對(duì)于粗糙壁面峰值熱流的計(jì)算，有理由做出以下假設(shè)：

圖7 不同流向位置壁面峰值熱流與高度的關(guān)系Fig.7 Relationship between surface peak heat flux and height at different flow direction positions

式中：第1 項(xiàng)表示由于粗糙引起的熱流增量，第2 項(xiàng)表示光滑平板的熱流。由式（6）得到相對(duì)峰值熱流增量為

式（8）與流向位置x 是無關(guān)的，符合圖5（b）的規(guī)律。

2. 2 不同寬度下的壁面熱流

對(duì)于第2 類工況，保持壁面高度不變，計(jì)算不同粗糙寬度的流場(chǎng)。使用的具體粗糙單元配置為：高度為0.421 2 mm，寬度分別為1.755、3.51、5.265、7.02、10.53 mm。對(duì)于小寬度，以0.4 m 為中心分布了100 個(gè)粗糙，在圖中對(duì)應(yīng)的曲線為這100 個(gè)粗糙的所在。

圖8 顯示了不同寬度下，粗糙壁面的熱流結(jié)果與光滑壁面（x=0.4～0.41 m）的比較?？梢钥闯?，當(dāng)粗糙帶高度固定時(shí)，寬度越大，壁面峰值熱流越小。

圖8 不同寬度下的壁面熱流Fig.8 Surface heat flux at different widths

圖9 給出了粗糙壁峰值熱流及相對(duì)峰值熱流增量沿流向的演變規(guī)律。可以看到，隨著粗糙寬度的增加，壁面峰值熱流逐漸減小，但是整體并不成線性關(guān)系；同時(shí)，相對(duì)峰值熱流增量也隨寬度增加而減小。

圖9 不同寬度下的壁面峰值熱流及相對(duì)峰值熱流增量Fig.9 Surface peak heat flux and relative surface peak heat flux increment at different widths

2. 3 固定寬高比下的壁面熱流

在得到上述關(guān)于熱流與寬度成負(fù)相關(guān)的規(guī)律后，需要進(jìn)一步確定壁面熱流與寬度的定量關(guān)系。文獻(xiàn)中常用Sigal-Danberg 參數(shù)來表征粗糙帶特征，這一參數(shù)以高度與寬度的比值來定量刻畫粗糙度特征。本節(jié)使用第3 類工況，目的是研究高度和寬度對(duì)粗糙壁面峰值熱流的影響強(qiáng)度以及這兩個(gè)因素之間比值規(guī)律，具體的粗糙壁參數(shù)見表2。

表2 固定寬高比的計(jì)算參數(shù)Table 2 Fixed aspect ratio calculation parameters

所取工況的粗糙寬高比均固定為25，假如寬度和高度對(duì)壁面熱流有相同的影響強(qiáng)度，那么在固定寬高比的情況下，幾種不同配置的粗糙壁面峰值熱流在同一流向位置應(yīng)是相同的。圖10（a）給出了5 個(gè)工況下的壁面熱流?？梢钥吹?，在固定寬高比的情況下，粗糙壁面峰值熱流并不相同。圖10（b）給出了每個(gè)工況的壁面峰值熱流隨流向x 的變化，可以發(fā)現(xiàn)高度越大，壁面峰值熱流越大。這說明熱流受高度的影響更大。同時(shí)，圖10（c）進(jìn)一步顯示了相對(duì)峰值熱流增量與流向位置的無關(guān)性。通過上述分析可以認(rèn)識(shí)到，壁面峰值熱流受高度的影響大于受寬度的影響，熱流與寬度的冪次關(guān)系應(yīng)該在-1～0 這一范圍。

圖10 相同寬高比下的壁面熱流Fig.10 Surface heat flux with the same aspect ratio

2. 4 不同馬赫數(shù)下的壁面熱流

將粗糙壁面的寬度取為10.53 mm，高度取為0.379 08 mm，分析了不同馬赫數(shù)（4、5、6、7、8、10）下的壁面峰值熱流。圖11（a）給出了粗糙壁面峰值熱流結(jié)果的比較?？梢园l(fā)現(xiàn)，同一流向位置，在不同馬赫數(shù)下，壁面峰值熱流隨馬赫數(shù)的變化并非單調(diào)。這是由于光滑壁的熱流本身隨馬赫數(shù)的變化就是非單調(diào)的。圖11（b）給出了不同馬赫數(shù)下，相對(duì)峰值熱流增量的演變規(guī)律，可以發(fā)現(xiàn)同一流向位置，在不同馬赫數(shù)下，相對(duì)峰值熱流增量隨馬赫數(shù)增大而單調(diào)減小，但與流向位置x 無關(guān)，這一規(guī)律與前文相同。

圖11 不同馬赫數(shù)下壁面峰值熱流及其相對(duì)峰值熱流增量Fig.11 Surface peak heat flux and relative surface peak heat flux increment at different Mach numbers

2. 5 不同壁溫下的壁面熱流

同樣，研究壁溫對(duì)于粗糙壁面熱流的影響時(shí)，固定粗糙壁面的為寬度10.53 mm、高度0.379 08 mm，計(jì)算不同壁溫（350、400、500、600、650、700、800 K）下的粗糙壁熱流。

圖12 顯示了不同壁溫下，全場(chǎng)均布粗糙的壁面峰值熱流和相對(duì)峰值熱流增量，可以發(fā)現(xiàn)，在同一流向位置，壁溫越高，壁面峰值熱流和相對(duì)峰值熱流越?。幌鄬?duì)峰值熱流增量近似與流向位置無關(guān)。

圖12 不同壁溫下的壁面峰值熱流及相對(duì)峰值熱流增量Fig.12 Surface peak heat flux and relative surface peak heat flux increment at different wall temperatures

3 熱流預(yù)測(cè)模型

第2 節(jié)中不同高度、不同寬度、相同寬高比、不同馬赫數(shù)以及不同壁溫下，粗糙壁面峰值熱流的定性變化規(guī)律可以概括為

1） 壁面峰值熱流的變化與粗糙壁的高度成正比，與流向位置x 的0.5 次成反比；相對(duì)峰值熱流增量與流向位置近似無關(guān)。

2） 壁面峰值熱流的變化與粗糙壁的寬度成反比。

3） 壁面峰值熱流受高度的影響大于受寬度的影響。

4） 相對(duì)峰值熱流增量隨馬赫數(shù)增高而減小。

5） 相對(duì)峰值熱流增量隨壁溫增高而減小。

在以上定性規(guī)律認(rèn)知的基礎(chǔ)上，為了進(jìn)一步分析得到定量預(yù)測(cè)粗糙壁面峰值熱流的模型，首先從光滑平板壁面熱流的理論解出發(fā)給出熱流預(yù)測(cè)公式的形式。

對(duì)于可壓縮平板邊界層相似性解，存在理論解用來計(jì)算熱流，其具體推導(dǎo)過程見文獻(xiàn)［28］。對(duì)于熱流系數(shù)ch，理論解為

式（10）即光滑平板壁面熱流的理論解。在之前的數(shù)值分析中，對(duì)粗糙壁面峰值熱流做了的假設(shè)，當(dāng)式（7）中的h=0 時(shí)，即表示光滑壁面的熱流，應(yīng)與式（10）相同，那么由此可得：

對(duì)于式（7）中第1 項(xiàng)粗糙壁引起的峰值熱流增量，從第2 節(jié)的數(shù)值結(jié)果來看，顯然，它的系數(shù)c1應(yīng)與來流參數(shù)（馬赫數(shù)、壁溫比）和粗糙帶的特征有關(guān)。結(jié)合第2 節(jié)中對(duì)粗糙帶高度及寬度的數(shù)值分析結(jié)果，作出以下假設(shè)：

相應(yīng)地，相對(duì)峰值熱流增量為

式（13）中需要確定的參數(shù)有a、b 和C。使用2.2 節(jié)、2.4 節(jié)與2.5 節(jié)中得到的全場(chǎng)波紋壁上的峰值熱流數(shù)值實(shí)驗(yàn)數(shù)據(jù)，計(jì)算其相對(duì)峰值熱流增量（此量?jī)H于馬赫數(shù)、壁溫比、粗糙壁特征有關(guān)，與流向位置無關(guān)）。對(duì)流向取平均，使用冪次關(guān)系擬合這個(gè)值與寬度、馬赫數(shù)、壁溫比的關(guān)系（分別對(duì)應(yīng)2.2 節(jié)、2.4 節(jié)與2.5 節(jié)中的工況）。擬合的結(jié)果如圖13 所示，圖13（a）～圖13（c）分別給出了對(duì)寬度、馬赫數(shù)以及壁溫的擬合結(jié)果。由擬合公式可以最終確定，C 為0.5，a 為2.67，b 為1.38。

圖13 熱流預(yù)測(cè)公式冪次關(guān)系的修正Fig.13 Modification of exponential relationship of heating prediction formula

圖14 給出了參數(shù)確定后，利用熱流式（13）反求得到的不同寬度、馬赫數(shù)、壁溫情況下的Δ?？梢园l(fā)現(xiàn)這一常數(shù)整體集中在300～400 之間，根據(jù)本文的假設(shè)，即Δ 為表示粗糙類型的常數(shù)，通過將所有工況（均為正弦粗糙）先沿流向位置平均，計(jì)算得到的Δ 再取平均，最終確定正弦壁面的這一常數(shù)為322。

圖14 熱流預(yù)測(cè)公式中Δ的計(jì)算Fig.14 Calculation of Δ in heating prediction formula

這樣可以得到預(yù)測(cè)粗糙帶壁面峰值熱流的公式：

對(duì)于正弦波紋壁Δ可取322。那么，在已知粗糙壁的形狀、位置、來流參數(shù)的前提下，可以使用式（15）預(yù)測(cè)由粗糙壁導(dǎo)致的壁面峰值熱流。

4 預(yù)測(cè)模型驗(yàn)證

本文提出的壁面峰值熱流預(yù)測(cè)公式需要進(jìn)一步進(jìn)行以下驗(yàn)證：一是基于相似性解得到的熱流公式，與實(shí)際直接數(shù)值模擬得到的熱流是否一致；二是式（15）在實(shí)際運(yùn)用到粗糙壁的熱流預(yù)測(cè)時(shí)，其誤差及適用范圍有多大。

4. 1 光滑壁面的熱流驗(yàn)證

圖15 給出了光滑壁面，在不同馬赫數(shù)下基于相似性解通過變換得到的熱流和熱流預(yù)測(cè)公式（15）的比較，結(jié)果顯示，兩者得到的熱流結(jié)果是一致的，這驗(yàn)證了理論推導(dǎo)過程的正確性。

圖15 相似性解熱流結(jié)果與預(yù)測(cè)公式結(jié)果Fig.15 Similarity surface heat flux results and prediction formula results

在光滑平板上，圖16 給出了相似性解的熱流結(jié)果與實(shí)際DNS 數(shù)值結(jié)果的比較，可以看出，兩者的熱流計(jì)算結(jié)果差距小于3%，吻合的較好，這驗(yàn)證了使用相似性解來計(jì)算熱流的可行性。

圖16 相似性解熱流結(jié)果與DNS 數(shù)值模擬結(jié)果Fig.16 Similarity solution surface heat flux results and DNS numerical simulation results

4. 2 全場(chǎng)均布波紋粗糙壁面的熱流驗(yàn)證

完成對(duì)于式（15）理論依據(jù)的驗(yàn)證后，對(duì)其在預(yù)測(cè)粗糙壁面峰值熱流時(shí)的誤差進(jìn)行分析。首先使用控制變量法，對(duì)高度、寬度、壁溫、馬赫數(shù)分別進(jìn)行驗(yàn)證。使用相對(duì)誤差以及絕對(duì)誤差來衡量其大小。理論公式預(yù)測(cè)的壁面峰值熱流與實(shí)際DNS 數(shù)值結(jié)果的相對(duì)誤差：

式中：chDNS為直接數(shù)值模擬得到的熱流系數(shù)。絕對(duì)誤差：

式中：Qw 為有量綱的熱流；Qwpredict為預(yù)測(cè)公式得到的有量綱的熱流。

圖17～圖21 給出了在不同高度、不同寬度、相同寬高比、不同馬赫數(shù)以及不同壁溫下的壁面峰值熱流的絕對(duì)誤差和同一流向位置實(shí)際與預(yù)測(cè)結(jié)果的比較。同時(shí)將其相對(duì)誤差在圖中標(biāo)出。可以看出，在大多數(shù)工況，除了在平板頭部由于黏性干擾導(dǎo)致誤差較大，下游的相對(duì)誤差都在5%以下。值得注意的是，在不同寬度下的相對(duì)誤差較大（＜10%），其原因可能為模型使用的用于衡量粗糙特征的冪次表達(dá)有一定的局限性，但絕對(duì)誤差也可保持在一個(gè)較小的范圍內(nèi)。相同寬高比、不同馬赫數(shù)以及不同壁溫下的預(yù)測(cè)結(jié)果誤差都在3%之內(nèi)。表明在已有工況下，式（15）預(yù)測(cè)效果不錯(cuò)。

圖17 不同高度下全場(chǎng)波紋壁預(yù)測(cè)公式與DNS 結(jié)果壁面峰值熱流誤差Fig.17 Surface peak heat flux error between full-field corrugated wall prediction formula and DNS results at different heights

圖18 不同寬度下全場(chǎng)波紋壁預(yù)測(cè)公式與DNS 結(jié)果壁面峰值熱流誤差Fig.18 Surface peak heat flux error between full-field corrugated wall prediction formula and DNS results at different widths

圖19 相同寬高比下的全場(chǎng)波紋壁預(yù)測(cè)公式與DNS 結(jié)果壁面峰值熱流誤差Fig.19 Surface peak heat flux error between full-field corrugated wall prediction formula and DNS results with same aspect ratio

圖20 不同馬赫數(shù)下全場(chǎng)波紋壁預(yù)測(cè)公式與DNS 結(jié)果壁面峰值熱流誤差Fig.20 Surface peak heat flux error between full-field corrugated wall prediction formula and DNS results at different Mach numbers

圖21 不同壁溫下全場(chǎng)波紋壁預(yù)測(cè)公式與DNS 結(jié)果壁面峰值熱流誤差Fig.21 Surface peak heat flux error between full-field corrugated wall prediction formula and DNS results at different wall temperatures

4. 3 測(cè)試工況的驗(yàn)證

為了進(jìn)一步驗(yàn)證式（15）對(duì)于更為復(fù)雜的工況，即多變量均變化情況下的適用性，對(duì)更多的測(cè)試工況進(jìn)行了驗(yàn)證。進(jìn)行了高度、壁溫、馬赫數(shù)均變化的工況的驗(yàn)證，計(jì)算工況見表3。

表3 測(cè)試工況的計(jì)算參數(shù)Table 3 Calculation parameters of test conditions

圖22 給出3 個(gè)驗(yàn)證工況下的壁面峰值熱流以及相對(duì)誤差結(jié)果，可以看出，在多變量變化的情況下，所提出的預(yù)測(cè)模型的誤差可以維持在較低的范圍內(nèi)（不超過3%），這也證明了式（15）對(duì)于多變量變化的復(fù)雜情況的適用性。

圖22 測(cè)試工況驗(yàn)證結(jié)果Fig.22 Test case verification results

4. 4 極限工況的驗(yàn)證

對(duì)建立預(yù)測(cè)模型工況的壁溫、馬赫數(shù)的最大和最小值，即極限的幾個(gè)工況，模型的適用性進(jìn)行進(jìn)一步驗(yàn)證的。使用的工況見表4。

表4 極限工況的計(jì)算參數(shù)Table 4 Calculation parameters of limited conditions

圖23 給出壁面峰值熱流比較以及模型預(yù)測(cè)的相對(duì)誤差?？梢园l(fā)現(xiàn)，實(shí)際數(shù)值模擬的壁面峰值熱流與預(yù)測(cè)公式的結(jié)果基本吻合；其中，馬赫數(shù)為4、壁溫為800 K 時(shí)的相對(duì)誤差較大，原因在于此時(shí)壁溫接近絕熱壁溫（841.51 K），熱流較小，相對(duì)偏差會(huì)更大一些。

圖23 極限工況驗(yàn)證結(jié)果Fig.23 Limited case verification results

4. 5 其他粗糙壁面形式的適用性

對(duì)于建立的熱流預(yù)測(cè)模型式（15），詳細(xì)驗(yàn)證了其對(duì)于正弦型連續(xù)波紋壁面的可用性，確定了其適用范圍，為了進(jìn)一步研究此模型對(duì)于其他壁面形式的適用性，數(shù)值計(jì)算分散的正弦型壁面以及橢圓型壁面，壁面形式如圖24 所示。圖24（a）為分散的正弦壁面，在0～1 m 的范圍內(nèi)，布置了50 個(gè)高度為0.210 6 mm、寬度為10.53 mm 的正弦型突起，相鄰粗糙元之間的間距為10.53 mm，在粗糙元附近進(jìn)行了局部加密；圖24（b）為分散的橢圓壁面，在0～1 m 的范圍內(nèi)，布置了40 個(gè)高度為0.210 6 mm、寬度為7.02 mm 的橢圓型凸起，相鄰粗糙元之間間距為14.04 mm，粗糙元附近進(jìn)行了局部加密。兩套網(wǎng)格流向分布2 501 個(gè)點(diǎn)，法向分布401 個(gè)點(diǎn)。

圖24 壁面形狀Fig.24 Wall shape

使用上述配置進(jìn)行數(shù)值計(jì)算，在計(jì)算分散的正弦壁面時(shí)，Δ取322，與之前的正弦波紋壁是一致的，這是符合Δ表征粗糙帶類型的常數(shù)這一假設(shè)的；而計(jì)算橢圓時(shí)，首先需要確定Δ，使用的方法與第4 節(jié)中方法相同，使用數(shù)值得到的相對(duì)熱流增量求得Δ，發(fā)現(xiàn)其分布為1 條與流向無關(guān)的直線，符合對(duì)于Δ的定義，對(duì)所有流向位置得到的Δ取平均，得到Δ為300。圖25 給出了不同壁面形式下粗糙附近的熱流分布，不同壁面形式下，粗糙附近的熱流變化是不同的。圖26 給出了2 種形式下壁面峰值熱流的變化規(guī)律，同時(shí)比對(duì)了其與預(yù)測(cè)值的差異，發(fā)現(xiàn)兩者是基本吻合的。圖27 則給出了壁面峰值熱流預(yù)測(cè)模型的相對(duì)誤差，可以看出，相對(duì)誤差在2%以下，說明預(yù)測(cè)公式對(duì)于其他2 種壁面形狀也具有一定的適用性。

圖25 不同壁面形式下的壁面熱流Fig.25 Surface heat flux with different wall shapes

圖26 不同壁面形式下的壁面峰值熱流及其預(yù)測(cè)Fig.26 Surface peak heat flux and its prediction with different wall shapes

圖27 不同壁面形式下的壁面峰值熱流預(yù)測(cè)誤差Fig.27 Prediction error of surface peak heat flux with different wall shapes

5 結(jié) 論

針對(duì)高超聲速平板層流邊界層，研究了表面均布高度小于600 μm 正弦型波紋凸起的粗糙導(dǎo)致的壁面熱流的變化規(guī)律，并提出了一種可用于定量預(yù)測(cè)由波紋粗糙壁引起的壁面峰值熱流的公式，詳細(xì)驗(yàn)證了其適用性。

使用直接數(shù)值模擬研究粗糙壁在不同馬赫數(shù)、高度、寬度以及壁溫下的壁面峰值熱流，發(fā)現(xiàn)壁面峰值熱流的變化與粗糙壁的高度成正比，而與粗糙壁的寬度成反比；壁面峰值熱流受高度的影響大于受寬度的影響；壁面相對(duì)峰值熱流增量隨馬赫數(shù)增大而減小；壁面相對(duì)峰值熱流增量隨壁溫增大而減小。

基于數(shù)值結(jié)果得到的定性規(guī)律，從光滑壁面熱流的理論解出發(fā)，假設(shè)熱流增量與粗糙高度、寬度、馬赫數(shù)和壁溫成冪次關(guān)系，得到預(yù)測(cè)粗糙壁面峰值熱流的公式，并根據(jù)數(shù)據(jù)擬合得到了各冪次關(guān)系中的系數(shù)。首先通過控制變量法，對(duì)于寬度、高度、馬赫數(shù)、壁溫等單一變量進(jìn)行了驗(yàn)證。然后驗(yàn)證了這些參數(shù)均變化時(shí)公式的有效性，并在極限工況下測(cè)試了其適用范圍。最后在2 種不同的壁面粗糙形式下進(jìn)一步探究了公式的可用性。整體上，預(yù)測(cè)模型的誤差維持在8%以下。

本文關(guān)注的粗糙元都是緩變的，對(duì)于高度與寬度大小相當(dāng)?shù)拇植谠?，本文的定性結(jié)論應(yīng)該能夠適用，但定量預(yù)測(cè)公式是否適用，還有待進(jìn)一步研究。

致 謝

感謝中國(guó)空氣動(dòng)力研究與發(fā)展中心國(guó)義軍研究員有啟發(fā)性的幫助。