求解滑移流區(qū)外掠圓柱體氣流流動與傳熱特性的數(shù)值模型

謝福壽,雷剛,邱一男,徐元元,王天祥,厲彥忠

(1.航天低溫推進劑技術(shù)國家重點實驗室,100028,北京;2.西安交通大學(xué)能源與動力工程學(xué)院,710049,西安)

近年來,關(guān)于稀薄氣體流動與傳熱問題日益受到人們的關(guān)注,最典型的4類案例為:①運載器穿過大氣層或再入大氣層的氣動加熱[1];②火星登陸器熱設(shè)計和熱管理[2];③微小精密熱線風(fēng)速儀研制[3];④高速真空隧道內(nèi)列車儀表的熱設(shè)計和熱管理[4]。對于運載火箭升空,隨著大氣密度的逐漸減小,氣體的稀薄效應(yīng)開始凸顯,連續(xù)介質(zhì)假設(shè)逐漸失效,此時氣體流態(tài)會依次經(jīng)歷連續(xù)流區(qū)(Kn<10-2)、滑移流區(qū)(10-210),如圖1所示。目前,由于稀薄氣體流動復(fù)雜,致使無法采用統(tǒng)一的算法或方程進行描述,故每個流態(tài)的求解方法不盡相同。

(a)不同特征長度下氣體流態(tài)劃分

(b)不同馬赫數(shù)和雷諾數(shù)下氣體流態(tài)劃分圖1 運載火箭穿越大氣層時氣體流態(tài)隨大氣高度的分布

滑移流區(qū)是指氣體與氣體之間的碰撞頻率,雖然其比氣體與物體表面間的碰撞頻率高相當(dāng)多,但稀薄效應(yīng)已不能再忽略,在靠近壁面處的非平衡效應(yīng)開始占主導(dǎo)作用,此時氣體流動與傳熱會呈現(xiàn)出兩個典型的特征:①速度滑移;②溫度跳躍。為了充分反映出該流域的典型特征,許多研究者采用壁面邊界條件修正的CFD技術(shù)進行流動與傳熱研究[5-10]。Kishore等通過速度滑移和定壁溫邊界條件修正對牛頓流體外掠球形顆粒的傳熱現(xiàn)象進行了數(shù)值仿真[5-6]。Choi等開發(fā)了嵌入Maxwell滑移邊界條件的非結(jié)構(gòu)化網(wǎng)格N-S方程求解模型,對圓柱體表面的速度滑移和溫度跳躍進行了研究[7]。索曉娜等采用連續(xù)介質(zhì)流動方程與壁面速度滑移和溫度跳躍邊界條件相結(jié)合描述微通道內(nèi)稀薄氣體二維可壓縮滑移流動[10]。Xie等對滑移流區(qū)稀薄空氣外掠圓柱體流動與傳熱特性展開了研究,其雷諾數(shù)范圍為0.001～20,克努森數(shù)范圍為0.01～0.1[11]。雖然許多學(xué)者針對滑移流區(qū)采用N-S方程進行了一定的研究,然而都沒有考慮壁面曲率對速度滑移和溫度跳躍的影響以及滑移邊界方程與CFD求解器之間坐標(biāo)系的轉(zhuǎn)換。Lockerby指出在目前稀薄氣體流動中,Maxwell最初滑移邊界條件被廣泛地誤用,其認為假如采用通用滑移邊界來模擬帶曲率的壁面問題時,一些關(guān)鍵的物理現(xiàn)象就會丟失:①通用滑移邊界條件計算出的切向滑移速度明顯低于Maxwell最初滑移條件得出的;②對于帶曲率的彎頭通道,采用Maxwell最初滑移邊界可以預(yù)測速度逆流現(xiàn)象,而通用滑移邊界條件則不能預(yù)測該流動行為;③對于外掠一個理想光滑球體,假定所有的入射分子都發(fā)生鏡面反射,即α=0,則采用Maxwell最初滑移邊界推導(dǎo)的Basset阻力方程預(yù)測表面摩擦阻力系數(shù)應(yīng)該為0,但采用通用滑移邊界條件推導(dǎo)的方程預(yù)測結(jié)果為負值,出現(xiàn)了非物理解[12]。

課題組前期搭建了一套低密度風(fēng)洞試驗裝置,開展了大量的實驗研究,獲得了一些有價值的研究結(jié)果[13-14]。然而,這些實驗結(jié)果并不能形象地展示出流場和溫度場的變化規(guī)律。當(dāng)克努森數(shù)增加時,也不能獲得壓力系數(shù)和傳熱系數(shù)沿圓柱體壁面變化的規(guī)律,更不能揭示流體掠過圓柱體流動與傳熱的機理。同時,由于實驗條件的限制,對總結(jié)出的經(jīng)驗關(guān)聯(lián)式,其雷諾數(shù)使用范圍仍然比較窄。因此,有必要開發(fā)新的數(shù)值模型來深入研究滑移流區(qū)流體流動與傳熱問題。

為了精確預(yù)測滑移流區(qū)稀薄氣體流動與傳熱問題,本文擬基于有限容積法的ANSYS FLUENT 18.0平臺,開發(fā)一套可用于該流域氣體外掠圓柱體流動與傳熱仿真的數(shù)值模型,其考慮壁面曲率的影響,解決了CFD求解器中一些關(guān)鍵性技術(shù)問題,如坐標(biāo)變換、UDF開發(fā)、UDF植入、求解發(fā)散等。在此基礎(chǔ)上,與課題組之前實驗數(shù)據(jù)進行詳細的對比驗證,充分驗證數(shù)值模型的可靠性和精確度,為計算流體動力學(xué)實現(xiàn)滑移流區(qū)氣體流動與傳熱問題研究構(gòu)建一座橋梁。

1 數(shù)值模型構(gòu)建

本文所研究的問題為稀薄氣體外掠圓柱體,直徑為D,流態(tài)區(qū)域尺寸及坐標(biāo)如圖2所示。圓柱體壁面溫度定義為Tw,空氣與圓柱體之間的溫度跳躍和速度滑移定義為Ts和us。

圖2 計算區(qū)域與邊界條件示意圖

1.1 層流N-S控制方程

對于滑移流區(qū),氣體橫掠圓柱體在主流區(qū)的流動和傳熱現(xiàn)象可由N-S方程和能量方程來描述,其中二維、層流、不可壓縮流體的無量綱表達式如下

(1)

(2)

(3)

(4)

式中:u*、v*分別為x、y方向的無量綱速度分量;p*為無量綱壓力;T*為無量綱溫度。

無量綱變量定義為

式中:u、v分別為x、y方向速度分量;u∞為自由流速度;p為空氣壓力;ρ∞為氣體密度;T為氣體溫度;T∞為自由流溫度;Cp為壓力系數(shù);p∞為自由流靜壓力;Cf為表面摩擦系數(shù);τw為壁面剪切力。

1.2 邊界條件

滑移流區(qū)氣體橫掠圓柱體流動與傳熱問題的邊界條件如圖2所示,詳細設(shè)置如下。

在計算區(qū)域左側(cè),均勻的速度入口邊界條件和氣體溫度設(shè)置為

u*=1;v*=0;T*=0

在計算區(qū)域右側(cè),默認氣體在出口處為充分發(fā)展?fàn)顟B(tài),壓力出口邊界條件和空氣溫度設(shè)置為

沿著計算區(qū)域頂端和底端,壁面設(shè)置為移動壁面,速度等于入口速度,即

u*=1;v*=0;T*=0

在圓柱體壁面,速度滑移表達式以速度分量的形式植入到移動壁面中,其具體的表達式為[15]

(5)

式中:un為壁面法向的氣體速度;n為壁面法向的坐標(biāo);s為壁面切向的坐標(biāo);σm為動量協(xié)調(diào)系數(shù),具體取值見1.3小節(jié);Lm為氣體分子平均自由程。方程(5)右側(cè)第1項是一階等溫速度滑移,由圓柱體壁面處的切向應(yīng)力推導(dǎo)而來,而第2項是由于切向溫度梯度引起的一個熱蠕動。需要說明的是,在帶曲率的物體中應(yīng)將導(dǎo)數(shù)項?un/?s考慮進去,然而在大多數(shù)已公布的工作中,該項要么被忽略,要么被導(dǎo)數(shù)項?us/?s錯誤地取代。Lockerby研究發(fā)現(xiàn),將導(dǎo)數(shù)項?un/?s代入滑移邊界方程之中,對于越稀薄的氣體,計算越準(zhǔn)確[12]。

溫度跳躍的表達式以壁面溫度修正的形式植入到壁面中,其具體的表達式為

(6)

式中:a為熱協(xié)調(diào)系數(shù),具體取值見1.3小節(jié)。

1.3 壁面邊界條件滑移修正

速度滑移邊界條件以一個移動壁面速度的形式加入,而溫度跳躍的邊界條件則以一個修正的定壁溫形式添加。該UDF采用C語言來編寫,以編譯的形式在Fluent求解器中執(zhí)行,其中所涉及到的速度梯度和溫度梯度分別從求解器中抽取,而速度滑移和溫度跳躍根據(jù)方程(5)(6)計算得到。

由于圓柱體表面存在曲率,所以自然坐標(biāo)系與求解器中默認的笛卡爾坐標(biāo)系是不統(tǒng)一的,正如圖3所示。為了表達自然坐標(biāo)系(s,n)中速度向量和所有的導(dǎo)數(shù)項,必須進行坐標(biāo)變換。笛卡爾坐標(biāo)系(x,y)中速度向量和所有的導(dǎo)數(shù)項可以從Fluent求解器中抽取。在坐標(biāo)變換之前,首先必須確定在壁面每個單元上兩個坐標(biāo)系之間的夾角θ。根據(jù)空間向量的性質(zhì),通過余弦函數(shù)則可以給出θ的表達式

(7)

式中:A=-An為邊界單元面積的法向向量,由Fluent?求解器提供;i、j為沿著x、y方向的單位向量。

圖3 在壁面處笛卡爾坐標(biāo)系與自然坐標(biāo)系及速度分量

使用鏈?zhǔn)椒▌t,可將方程(5)(6)的導(dǎo)數(shù)項表示為

(8)

(9)

(10)

(11)

自然坐標(biāo)系與笛卡爾坐標(biāo)系關(guān)系如下

(12)

在s、n方向?qū)、y分別求偏導(dǎo)數(shù)的表達式為

(13)

在笛卡爾坐標(biāo)系和自然坐標(biāo)系中,向量u可分別表達為

uxi+uyj

(14)

usk+ung

(15)

式中:k和g為s、n方向上的單位向量。

根據(jù)幾何關(guān)系,如4所示,可獲得向量(i,j)和(k,g)之間的關(guān)系如下

(16)

(a)自然坐標(biāo)系對i的表達 (b)自然坐標(biāo)系對j的表達圖4 單位向量變換

最后,將方程(16)代入方程(14),化解整理后,即可得到(ux,uy)和(us,un)之間的關(guān)系為

u=(uxcosθ+uysinθ)k+(-uxsinθ+uycosθ)g

(17)

于是

(18)

將方程(13)(18)代入方程(8)～方程(11),所有的導(dǎo)數(shù)項即可整理為

(19)

(20)

(21)

(22)

通過訪問ANSYS Fluent求解器中的數(shù)據(jù)內(nèi)存,可獲得速度和溫度的全部偏導(dǎo)數(shù)項(?ux/?x,?ux/?y,?uy/?x,?uy/?y,?T/?x,?T/?y)。對于角度偏導(dǎo)數(shù)(?θ/?x和?θ/?y),則需要通過極坐標(biāo)系來確定。

如圖5所示,笛卡爾坐標(biāo)系中(x,y)與極坐標(biāo)系中(r,θ)之間的關(guān)系可依據(jù)幾何關(guān)系來表示

(23)

式中r為圓的半徑。于是,角度偏導(dǎo)數(shù)可表示為

(24)

圖5 笛卡爾坐標(biāo)系與極坐標(biāo)系之間的變換

需要說明的是,動量協(xié)調(diào)系數(shù)σm和熱協(xié)調(diào)系數(shù)a一般取決于許多因素,如氣體與物體表面材料的種類、壁面光滑度、表面的清潔程度(氧化程度、是否有污染物覆蓋、是否吸附其他氣體)、溫度、壓力等以及一些尚未被人們知曉的因素,且兩個系數(shù)都在0(鏡面反射)到1(漫反射)之間變化。根據(jù)文獻[16-17]報道,對于大多數(shù)工程應(yīng)用來說,這兩個系數(shù)都非常接近于1,故在目前的工作中這兩個系數(shù)都取為1。

1.4 CFD數(shù)值計算方案確定

采用商業(yè)計算流體動力學(xué)代碼——ANYSYS Fluent展開數(shù)值仿真,即聯(lián)合邊界條件求解控制方程(1)～方程(4),最終獲得了稀薄氣體外掠圓柱體的壓力場、速度場和溫度場?？臻g離散梯度方案選擇基于單元最小二乘法梯度,壓力離散方案選擇二階格式,動量方程和能量方程離散方案選擇二階迎風(fēng)格式。在Fluent中,默認收斂標(biāo)準(zhǔn)是足以用來判斷收斂的,即連續(xù)方程和動量方程的全局殘差均降到了10-11,而能量方程的全局殘差降到了10-15,則認為迭代收斂。

2 結(jié)果討論與分析

2.1 連續(xù)流低壓狀態(tài)下定克努森數(shù)結(jié)果比較

圖6展示了連續(xù)流區(qū)低壓狀態(tài)下無滑移邊界數(shù)值預(yù)測值與滑移邊界數(shù)值預(yù)測值的比較結(jié)果,其中低壓狀態(tài)的壓力為50 kPa,所對應(yīng)的克努森數(shù)為0.006。從圖中對比可知,只要氣體流態(tài)處于連續(xù)流區(qū)時,無論是常壓狀態(tài)還是低壓狀態(tài),無滑移邊界與滑移邊界的計算結(jié)果都非常相近,兩者相對誤差在1.5%以內(nèi)。

圖6 連續(xù)流區(qū)數(shù)值計算結(jié)果

2.2 滑移流區(qū)低壓狀態(tài)下定克努森數(shù)結(jié)果比較

圖7呈現(xiàn)了在滑移流區(qū)狀態(tài)下無滑移邊界數(shù)值預(yù)測值、滑移邊界數(shù)值預(yù)測值和實驗結(jié)果的對比?；屏鲄^(qū)低壓狀態(tài)分別是指氣壓壓力為6 kPa和3 kPa,所對應(yīng)的克努森數(shù)分別為0.054和0.108。從圖中可看出,當(dāng)氣體流態(tài)由連續(xù)流區(qū)進入到滑移流區(qū)后,采用無滑移邊界條件的計算結(jié)果與實驗結(jié)果的誤差已達到±15%,尤其是壓力越低,克努森數(shù)越大,計算誤差越大。即使在同一壓力下,采用無滑移邊界條件的計算結(jié)果與實驗結(jié)果的誤差也會隨著雷諾數(shù)增大而增大。采用滑移邊界條件的計算結(jié)果可將誤差控制在±2.5%以內(nèi),所以在稀薄氣體滑移區(qū)的流動換熱問題中,本文修正的滑移流區(qū)數(shù)值模型可保障數(shù)值計算的準(zhǔn)確性。

(a)6 kPa

(b)3 kPa圖7 滑移流區(qū)數(shù)值計算結(jié)果

2.3 過渡流區(qū)低壓狀態(tài)下定克努森數(shù)結(jié)果比較

圖8給出了過渡流區(qū)狀態(tài)下無滑移邊界數(shù)值預(yù)測值、滑移邊界數(shù)值預(yù)測值和實驗結(jié)果的對比。過渡流區(qū)低壓狀態(tài)是指氣體壓力為1 kPa,所對應(yīng)的克努森數(shù)為0.323。從圖中可看出,當(dāng)?shù)蛪簹怏w的流態(tài)處于過渡流區(qū)時,采用無滑移邊界條件的計算結(jié)果與實驗結(jié)果的相對誤差已達到±60%,即使采用滑移邊界條件修正,也只能將相對誤差縮小到±20%以內(nèi)。

圖8 過渡流區(qū)數(shù)值計算結(jié)果

2.4 跨流域變克努森數(shù)結(jié)果比較

為了進一步驗證該數(shù)值模型的可靠性和精確度,維持氣體速度(3 m/s)不變,通過改變氣體壓力,使壓力在100 Pa～100 kPa之間變化,獲得了數(shù)值預(yù)測值、實驗值,并與傳統(tǒng)經(jīng)驗關(guān)聯(lián)式(Fand[18],Hilpert[19],McAdams[20],Collis[21])預(yù)測值進行了比較,如圖9所示。從圖中可清楚地看到,在連續(xù)流區(qū),數(shù)值計算結(jié)果與本文實驗值、Collis的經(jīng)驗關(guān)聯(lián)式預(yù)測值和Fand的經(jīng)驗關(guān)聯(lián)式預(yù)測值都吻合良好,說明數(shù)值計算過程都是合理的,包括邊界條件設(shè)置、求解方法選擇、物性設(shè)置、區(qū)域無關(guān)性研究和網(wǎng)格無關(guān)性研究等。

圖9 跨流域變克努森數(shù)結(jié)果比較

在滑移流區(qū),采用滑移邊界條件修正的數(shù)值計算結(jié)果與本文實驗值吻合良好,其最大相對偏差處于±2.5%以內(nèi),而采用無滑移邊界條件的數(shù)值計算結(jié)果逐漸偏離實驗值,最終達到10%的相對偏差。也就是說,通過滑移邊界條件修正,數(shù)值計算結(jié)果的精度有了顯著提升。

在過渡流區(qū),采用滑移邊界條件修正的數(shù)值計算結(jié)果與本文實驗結(jié)果相比,計算偏差逐漸增加,而采用無滑移邊界條件的數(shù)值計算結(jié)果則完全偏離了本文實驗結(jié)果,因此在該流態(tài)下,數(shù)值計算模型需要進一步修正。然而,本文側(cè)重于滑移流區(qū)低壓空氣外掠圓柱體流動與傳熱特性研究,所以在此不考慮過渡流區(qū)。

2.5 滑移效應(yīng)對低壓氣體外掠圓柱體的流動影響

圖10展現(xiàn)了在滑移流區(qū)當(dāng)雷諾數(shù)為1時無量綱滑移速度沿圓柱體壁面的分布規(guī)律。從圖中可看出,當(dāng)克努森數(shù)保持不變時,無量綱滑移速度在前止滯點處為0,隨著角度α的增加,其逐漸增加到某一最大值,然后開始減小,直到在后止滯點處變?yōu)?。因為在前止滯點處,空氣動壓轉(zhuǎn)變?yōu)殪o壓,使速度變?yōu)?。隨后,靜壓會逐漸再次轉(zhuǎn)變?yōu)閯訅?使壁面的切向速度最終達到一個最大值。由于稀薄效應(yīng)的影響,最大無量綱滑移速度會隨著克努森數(shù)的增加而增大,當(dāng)克努森數(shù)達到0.1時,最大無量綱滑移速度可達到自由來流速度的16%。

圖10 無量綱滑移速度沿圓柱體壁面的分布

圖11描述了在滑移流區(qū)雷諾數(shù)為1時壓力系數(shù)沿圓柱體壁面的分布規(guī)律。當(dāng)克努森數(shù)保持不變時,壓力系數(shù)最大值位于圓柱體前止滯點處,且沿著圓柱體壁面逐漸減小。然而,由于滑移效應(yīng)的影響,圓柱體前止滯點處的壓力系數(shù)會隨著克努森數(shù)的增加而逐漸減小,不過克努森數(shù)對壓力系數(shù)的影響其實并不大。

圖12 表面摩擦系數(shù)沿圓柱體壁面的分布

圖12描述了在滑移流區(qū)雷諾數(shù)為1時表面摩擦系數(shù)沿圓柱體壁面的分布規(guī)律。從圖中可知,當(dāng)克努森數(shù)一定時,由于雷諾數(shù)較低,黏性力起主導(dǎo)作用,從而抑制邊界層分離,使圓柱體前止滯點與后止滯點速度為0,故表面摩擦系數(shù)沿著圓柱體壁面先從0增加至最大值,然后逐漸減小至0。由于稀薄效應(yīng)影響,最大表面摩擦系數(shù)隨著克努森數(shù)增加而逐漸減小,當(dāng)克努森數(shù)達到0.1時,最大表面摩擦系數(shù)比連續(xù)流區(qū)最大表面摩擦系數(shù)減小了38%。

3 結(jié) 論

針對滑移流區(qū)氣體外掠圓柱體問題,構(gòu)建了一套預(yù)測低壓氣體流動與傳熱的數(shù)值模型,考慮了軸向溫度梯度引起的熱蠕動影響和壁面曲率的影響,通過與實驗數(shù)據(jù)對比,得出以下結(jié)論:

(1)當(dāng)氣體流態(tài)處于連續(xù)流區(qū)時,滑移邊界修正模型的計算結(jié)果與無滑移邊界的計算結(jié)果相近;

(2)當(dāng)氣體流態(tài)處于滑移流區(qū)時,預(yù)測外掠圓柱體傳熱的傳統(tǒng)經(jīng)驗關(guān)聯(lián)式已不再適用,此時采用滑移邊界修正的數(shù)值模型,可實現(xiàn)對稀薄氣體滑移流區(qū)流動與換熱的精確預(yù)測,且數(shù)值計算結(jié)果明顯優(yōu)于傳統(tǒng)經(jīng)驗關(guān)聯(lián)式預(yù)測結(jié)果;

(3)在滑移流區(qū),稀薄效應(yīng)對低壓氣體外掠圓柱體的流體動力學(xué)特性影響顯著,隨著克努森數(shù)增加,最大無量綱滑移速度呈線性增加,而最大表面摩擦系數(shù)呈現(xiàn)線性減小。