謝福壽,雷剛,邱一男,徐元元,王天祥,厲彥忠
(1.航天低溫推進(jìn)劑技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室,100028,北京;2.西安交通大學(xué)能源與動(dòng)力工程學(xué)院,710049,西安)
近年來(lái),關(guān)于稀薄氣體流動(dòng)與傳熱問(wèn)題日益受到人們的關(guān)注,最典型的4類案例為:①運(yùn)載器穿過(guò)大氣層或再入大氣層的氣動(dòng)加熱[1];②火星登陸器熱設(shè)計(jì)和熱管理[2];③微小精密熱線風(fēng)速儀研制[3];④高速真空隧道內(nèi)列車儀表的熱設(shè)計(jì)和熱管理[4]。對(duì)于運(yùn)載火箭升空,隨著大氣密度的逐漸減小,氣體的稀薄效應(yīng)開始凸顯,連續(xù)介質(zhì)假設(shè)逐漸失效,此時(shí)氣體流態(tài)會(huì)依次經(jīng)歷連續(xù)流區(qū)(Kn<10-2)、滑移流區(qū)(10-2
(a)不同特征長(zhǎng)度下氣體流態(tài)劃分
(b)不同馬赫數(shù)和雷諾數(shù)下氣體流態(tài)劃分圖1 運(yùn)載火箭穿越大氣層時(shí)氣體流態(tài)隨大氣高度的分布
滑移流區(qū)是指氣體與氣體之間的碰撞頻率,雖然其比氣體與物體表面間的碰撞頻率高相當(dāng)多,但稀薄效應(yīng)已不能再忽略,在靠近壁面處的非平衡效應(yīng)開始占主導(dǎo)作用,此時(shí)氣體流動(dòng)與傳熱會(huì)呈現(xiàn)出兩個(gè)典型的特征:①速度滑移;②溫度跳躍。為了充分反映出該流域的典型特征,許多研究者采用壁面邊界條件修正的CFD技術(shù)進(jìn)行流動(dòng)與傳熱研究[5-10]。Kishore等通過(guò)速度滑移和定壁溫邊界條件修正對(duì)牛頓流體外掠球形顆粒的傳熱現(xiàn)象進(jìn)行了數(shù)值仿真[5-6]。Choi等開發(fā)了嵌入Maxwell滑移邊界條件的非結(jié)構(gòu)化網(wǎng)格N-S方程求解模型,對(duì)圓柱體表面的速度滑移和溫度跳躍進(jìn)行了研究[7]。索曉娜等采用連續(xù)介質(zhì)流動(dòng)方程與壁面速度滑移和溫度跳躍邊界條件相結(jié)合描述微通道內(nèi)稀薄氣體二維可壓縮滑移流動(dòng)[10]。Xie等對(duì)滑移流區(qū)稀薄空氣外掠圓柱體流動(dòng)與傳熱特性展開了研究,其雷諾數(shù)范圍為0.001~20,克努森數(shù)范圍為0.01~0.1[11]。雖然許多學(xué)者針對(duì)滑移流區(qū)采用N-S方程進(jìn)行了一定的研究,然而都沒(méi)有考慮壁面曲率對(duì)速度滑移和溫度跳躍的影響以及滑移邊界方程與CFD求解器之間坐標(biāo)系的轉(zhuǎn)換。Lockerby指出在目前稀薄氣體流動(dòng)中,Maxwell最初滑移邊界條件被廣泛地誤用,其認(rèn)為假如采用通用滑移邊界來(lái)模擬帶曲率的壁面問(wèn)題時(shí),一些關(guān)鍵的物理現(xiàn)象就會(huì)丟失:①通用滑移邊界條件計(jì)算出的切向滑移速度明顯低于Maxwell最初滑移條件得出的;②對(duì)于帶曲率的彎頭通道,采用Maxwell最初滑移邊界可以預(yù)測(cè)速度逆流現(xiàn)象,而通用滑移邊界條件則不能預(yù)測(cè)該流動(dòng)行為;③對(duì)于外掠一個(gè)理想光滑球體,假定所有的入射分子都發(fā)生鏡面反射,即α=0,則采用Maxwell最初滑移邊界推導(dǎo)的Basset阻力方程預(yù)測(cè)表面摩擦阻力系數(shù)應(yīng)該為0,但采用通用滑移邊界條件推導(dǎo)的方程預(yù)測(cè)結(jié)果為負(fù)值,出現(xiàn)了非物理解[12]。
課題組前期搭建了一套低密度風(fēng)洞試驗(yàn)裝置,開展了大量的實(shí)驗(yàn)研究,獲得了一些有價(jià)值的研究結(jié)果[13-14]。然而,這些實(shí)驗(yàn)結(jié)果并不能形象地展示出流場(chǎng)和溫度場(chǎng)的變化規(guī)律。當(dāng)克努森數(shù)增加時(shí),也不能獲得壓力系數(shù)和傳熱系數(shù)沿圓柱體壁面變化的規(guī)律,更不能揭示流體掠過(guò)圓柱體流動(dòng)與傳熱的機(jī)理。同時(shí),由于實(shí)驗(yàn)條件的限制,對(duì)總結(jié)出的經(jīng)驗(yàn)關(guān)聯(lián)式,其雷諾數(shù)使用范圍仍然比較窄。因此,有必要開發(fā)新的數(shù)值模型來(lái)深入研究滑移流區(qū)流體流動(dòng)與傳熱問(wèn)題。
為了精確預(yù)測(cè)滑移流區(qū)稀薄氣體流動(dòng)與傳熱問(wèn)題,本文擬基于有限容積法的ANSYS FLUENT 18.0平臺(tái),開發(fā)一套可用于該流域氣體外掠圓柱體流動(dòng)與傳熱仿真的數(shù)值模型,其考慮壁面曲率的影響,解決了CFD求解器中一些關(guān)鍵性技術(shù)問(wèn)題,如坐標(biāo)變換、UDF開發(fā)、UDF植入、求解發(fā)散等。在此基礎(chǔ)上,與課題組之前實(shí)驗(yàn)數(shù)據(jù)進(jìn)行詳細(xì)的對(duì)比驗(yàn)證,充分驗(yàn)證數(shù)值模型的可靠性和精確度,為計(jì)算流體動(dòng)力學(xué)實(shí)現(xiàn)滑移流區(qū)氣體流動(dòng)與傳熱問(wèn)題研究構(gòu)建一座橋梁。
本文所研究的問(wèn)題為稀薄氣體外掠圓柱體,直徑為D,流態(tài)區(qū)域尺寸及坐標(biāo)如圖2所示。圓柱體壁面溫度定義為Tw,空氣與圓柱體之間的溫度跳躍和速度滑移定義為Ts和us。
圖2 計(jì)算區(qū)域與邊界條件示意圖
對(duì)于滑移流區(qū),氣體橫掠圓柱體在主流區(qū)的流動(dòng)和傳熱現(xiàn)象可由N-S方程和能量方程來(lái)描述,其中二維、層流、不可壓縮流體的無(wú)量綱表達(dá)式如下
(1)
(2)
(3)
(4)
式中:u*、v*分別為x、y方向的無(wú)量綱速度分量;p*為無(wú)量綱壓力;T*為無(wú)量綱溫度。
無(wú)量綱變量定義為
式中:u、v分別為x、y方向速度分量;u∞為自由流速度;p為空氣壓力;ρ∞為氣體密度;T為氣體溫度;T∞為自由流溫度;Cp為壓力系數(shù);p∞為自由流靜壓力;Cf為表面摩擦系數(shù);τw為壁面剪切力。
滑移流區(qū)氣體橫掠圓柱體流動(dòng)與傳熱問(wèn)題的邊界條件如圖2所示,詳細(xì)設(shè)置如下。
在計(jì)算區(qū)域左側(cè),均勻的速度入口邊界條件和氣體溫度設(shè)置為
u*=1;v*=0;T*=0
在計(jì)算區(qū)域右側(cè),默認(rèn)氣體在出口處為充分發(fā)展?fàn)顟B(tài),壓力出口邊界條件和空氣溫度設(shè)置為
沿著計(jì)算區(qū)域頂端和底端,壁面設(shè)置為移動(dòng)壁面,速度等于入口速度,即
u*=1;v*=0;T*=0
在圓柱體壁面,速度滑移表達(dá)式以速度分量的形式植入到移動(dòng)壁面中,其具體的表達(dá)式為[15]
(5)
式中:un為壁面法向的氣體速度;n為壁面法向的坐標(biāo);s為壁面切向的坐標(biāo);σm為動(dòng)量協(xié)調(diào)系數(shù),具體取值見1.3小節(jié);Lm為氣體分子平均自由程。方程(5)右側(cè)第1項(xiàng)是一階等溫速度滑移,由圓柱體壁面處的切向應(yīng)力推導(dǎo)而來(lái),而第2項(xiàng)是由于切向溫度梯度引起的一個(gè)熱蠕動(dòng)。需要說(shuō)明的是,在帶曲率的物體中應(yīng)將導(dǎo)數(shù)項(xiàng)?un/?s考慮進(jìn)去,然而在大多數(shù)已公布的工作中,該項(xiàng)要么被忽略,要么被導(dǎo)數(shù)項(xiàng)?us/?s錯(cuò)誤地取代。Lockerby研究發(fā)現(xiàn),將導(dǎo)數(shù)項(xiàng)?un/?s代入滑移邊界方程之中,對(duì)于越稀薄的氣體,計(jì)算越準(zhǔn)確[12]。
溫度跳躍的表達(dá)式以壁面溫度修正的形式植入到壁面中,其具體的表達(dá)式為
(6)
式中:a為熱協(xié)調(diào)系數(shù),具體取值見1.3小節(jié)。
速度滑移邊界條件以一個(gè)移動(dòng)壁面速度的形式加入,而溫度跳躍的邊界條件則以一個(gè)修正的定壁溫形式添加。該UDF采用C語(yǔ)言來(lái)編寫,以編譯的形式在Fluent求解器中執(zhí)行,其中所涉及到的速度梯度和溫度梯度分別從求解器中抽取,而速度滑移和溫度跳躍根據(jù)方程(5)(6)計(jì)算得到。
由于圓柱體表面存在曲率,所以自然坐標(biāo)系與求解器中默認(rèn)的笛卡爾坐標(biāo)系是不統(tǒng)一的,正如圖3所示。為了表達(dá)自然坐標(biāo)系(s,n)中速度向量和所有的導(dǎo)數(shù)項(xiàng),必須進(jìn)行坐標(biāo)變換。笛卡爾坐標(biāo)系(x,y)中速度向量和所有的導(dǎo)數(shù)項(xiàng)可以從Fluent求解器中抽取。在坐標(biāo)變換之前,首先必須確定在壁面每個(gè)單元上兩個(gè)坐標(biāo)系之間的夾角θ。根據(jù)空間向量的性質(zhì),通過(guò)余弦函數(shù)則可以給出θ的表達(dá)式
(7)
式中:A=-An為邊界單元面積的法向向量,由Fluent?求解器提供;i、j為沿著x、y方向的單位向量。
圖3 在壁面處笛卡爾坐標(biāo)系與自然坐標(biāo)系及速度分量
使用鏈?zhǔn)椒▌t,可將方程(5)(6)的導(dǎo)數(shù)項(xiàng)表示為
(8)
(9)
(10)
(11)
自然坐標(biāo)系與笛卡爾坐標(biāo)系關(guān)系如下
(12)
在s、n方向?qū)、y分別求偏導(dǎo)數(shù)的表達(dá)式為
(13)
在笛卡爾坐標(biāo)系和自然坐標(biāo)系中,向量u可分別表達(dá)為
uxi+uyj
(14)
usk+ung
(15)
式中:k和g為s、n方向上的單位向量。
根據(jù)幾何關(guān)系,如4所示,可獲得向量(i,j)和(k,g)之間的關(guān)系如下
(16)
(a)自然坐標(biāo)系對(duì)i的表達(dá) (b)自然坐標(biāo)系對(duì)j的表達(dá)圖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ù)項(xiàng)即可整理為
(19)
(20)
(21)
(22)
通過(guò)訪問(wèn)ANSYS Fluent求解器中的數(shù)據(jù)內(nèi)存,可獲得速度和溫度的全部偏導(dǎo)數(shù)項(xiàng)(?ux/?x,?ux/?y,?uy/?x,?uy/?y,?T/?x,?T/?y)。對(duì)于角度偏導(dǎo)數(shù)(?θ/?x和?θ/?y),則需要通過(guò)極坐標(biāo)系來(lái)確定。
如圖5所示,笛卡爾坐標(biāo)系中(x,y)與極坐標(biāo)系中(r,θ)之間的關(guān)系可依據(jù)幾何關(guān)系來(lái)表示
(23)
式中r為圓的半徑。于是,角度偏導(dǎo)數(shù)可表示為
(24)
圖5 笛卡爾坐標(biāo)系與極坐標(biāo)系之間的變換
需要說(shuō)明的是,動(dòng)量協(xié)調(diào)系數(shù)σm和熱協(xié)調(diào)系數(shù)a一般取決于許多因素,如氣體與物體表面材料的種類、壁面光滑度、表面的清潔程度(氧化程度、是否有污染物覆蓋、是否吸附其他氣體)、溫度、壓力等以及一些尚未被人們知曉的因素,且兩個(gè)系數(shù)都在0(鏡面反射)到1(漫反射)之間變化。根據(jù)文獻(xiàn)[16-17]報(bào)道,對(duì)于大多數(shù)工程應(yīng)用來(lái)說(shuō),這兩個(gè)系數(shù)都非常接近于1,故在目前的工作中這兩個(gè)系數(shù)都取為1。
采用商業(yè)計(jì)算流體動(dòng)力學(xué)代碼——ANYSYS Fluent展開數(shù)值仿真,即聯(lián)合邊界條件求解控制方程(1)~方程(4),最終獲得了稀薄氣體外掠圓柱體的壓力場(chǎng)、速度場(chǎng)和溫度場(chǎng)??臻g離散梯度方案選擇基于單元最小二乘法梯度,壓力離散方案選擇二階格式,動(dòng)量方程和能量方程離散方案選擇二階迎風(fēng)格式。在Fluent中,默認(rèn)收斂標(biāo)準(zhǔn)是足以用來(lái)判斷收斂的,即連續(xù)方程和動(dòng)量方程的全局殘差均降到了10-11,而能量方程的全局殘差降到了10-15,則認(rèn)為迭代收斂。
圖6展示了連續(xù)流區(qū)低壓狀態(tài)下無(wú)滑移邊界數(shù)值預(yù)測(cè)值與滑移邊界數(shù)值預(yù)測(cè)值的比較結(jié)果,其中低壓狀態(tài)的壓力為50 kPa,所對(duì)應(yīng)的克努森數(shù)為0.006。從圖中對(duì)比可知,只要?dú)怏w流態(tài)處于連續(xù)流區(qū)時(shí),無(wú)論是常壓狀態(tài)還是低壓狀態(tài),無(wú)滑移邊界與滑移邊界的計(jì)算結(jié)果都非常相近,兩者相對(duì)誤差在1.5%以內(nèi)。
圖6 連續(xù)流區(qū)數(shù)值計(jì)算結(jié)果
圖7呈現(xiàn)了在滑移流區(qū)狀態(tài)下無(wú)滑移邊界數(shù)值預(yù)測(cè)值、滑移邊界數(shù)值預(yù)測(cè)值和實(shí)驗(yàn)結(jié)果的對(duì)比。滑移流區(qū)低壓狀態(tài)分別是指氣壓壓力為6 kPa和3 kPa,所對(duì)應(yīng)的克努森數(shù)分別為0.054和0.108。從圖中可看出,當(dāng)氣體流態(tài)由連續(xù)流區(qū)進(jìn)入到滑移流區(qū)后,采用無(wú)滑移邊界條件的計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果的誤差已達(dá)到±15%,尤其是壓力越低,克努森數(shù)越大,計(jì)算誤差越大。即使在同一壓力下,采用無(wú)滑移邊界條件的計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果的誤差也會(huì)隨著雷諾數(shù)增大而增大。采用滑移邊界條件的計(jì)算結(jié)果可將誤差控制在±2.5%以內(nèi),所以在稀薄氣體滑移區(qū)的流動(dòng)換熱問(wèn)題中,本文修正的滑移流區(qū)數(shù)值模型可保障數(shù)值計(jì)算的準(zhǔn)確性。
(a)6 kPa
(b)3 kPa圖7 滑移流區(qū)數(shù)值計(jì)算結(jié)果
圖8給出了過(guò)渡流區(qū)狀態(tài)下無(wú)滑移邊界數(shù)值預(yù)測(cè)值、滑移邊界數(shù)值預(yù)測(cè)值和實(shí)驗(yàn)結(jié)果的對(duì)比。過(guò)渡流區(qū)低壓狀態(tài)是指氣體壓力為1 kPa,所對(duì)應(yīng)的克努森數(shù)為0.323。從圖中可看出,當(dāng)?shù)蛪簹怏w的流態(tài)處于過(guò)渡流區(qū)時(shí),采用無(wú)滑移邊界條件的計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果的相對(duì)誤差已達(dá)到±60%,即使采用滑移邊界條件修正,也只能將相對(duì)誤差縮小到±20%以內(nèi)。
圖8 過(guò)渡流區(qū)數(shù)值計(jì)算結(jié)果
為了進(jìn)一步驗(yàn)證該數(shù)值模型的可靠性和精確度,維持氣體速度(3 m/s)不變,通過(guò)改變氣體壓力,使壓力在100 Pa~100 kPa之間變化,獲得了數(shù)值預(yù)測(cè)值、實(shí)驗(yàn)值,并與傳統(tǒng)經(jīng)驗(yàn)關(guān)聯(lián)式(Fand[18],Hilpert[19],McAdams[20],Collis[21])預(yù)測(cè)值進(jìn)行了比較,如圖9所示。從圖中可清楚地看到,在連續(xù)流區(qū),數(shù)值計(jì)算結(jié)果與本文實(shí)驗(yàn)值、Collis的經(jīng)驗(yàn)關(guān)聯(lián)式預(yù)測(cè)值和Fand的經(jīng)驗(yàn)關(guān)聯(lián)式預(yù)測(cè)值都吻合良好,說(shuō)明數(shù)值計(jì)算過(guò)程都是合理的,包括邊界條件設(shè)置、求解方法選擇、物性設(shè)置、區(qū)域無(wú)關(guān)性研究和網(wǎng)格無(wú)關(guān)性研究等。
圖9 跨流域變克努森數(shù)結(jié)果比較
在滑移流區(qū),采用滑移邊界條件修正的數(shù)值計(jì)算結(jié)果與本文實(shí)驗(yàn)值吻合良好,其最大相對(duì)偏差處于±2.5%以內(nèi),而采用無(wú)滑移邊界條件的數(shù)值計(jì)算結(jié)果逐漸偏離實(shí)驗(yàn)值,最終達(dá)到10%的相對(duì)偏差。也就是說(shuō),通過(guò)滑移邊界條件修正,數(shù)值計(jì)算結(jié)果的精度有了顯著提升。
在過(guò)渡流區(qū),采用滑移邊界條件修正的數(shù)值計(jì)算結(jié)果與本文實(shí)驗(yàn)結(jié)果相比,計(jì)算偏差逐漸增加,而采用無(wú)滑移邊界條件的數(shù)值計(jì)算結(jié)果則完全偏離了本文實(shí)驗(yàn)結(jié)果,因此在該流態(tài)下,數(shù)值計(jì)算模型需要進(jìn)一步修正。然而,本文側(cè)重于滑移流區(qū)低壓空氣外掠圓柱體流動(dòng)與傳熱特性研究,所以在此不考慮過(guò)渡流區(qū)。
圖10展現(xiàn)了在滑移流區(qū)當(dāng)雷諾數(shù)為1時(shí)無(wú)量綱滑移速度沿圓柱體壁面的分布規(guī)律。從圖中可看出,當(dāng)克努森數(shù)保持不變時(shí),無(wú)量綱滑移速度在前止滯點(diǎn)處為0,隨著角度α的增加,其逐漸增加到某一最大值,然后開始減小,直到在后止滯點(diǎn)處變?yōu)?。因?yàn)樵谇爸箿c(diǎn)處,空氣動(dòng)壓轉(zhuǎn)變?yōu)殪o壓,使速度變?yōu)?。隨后,靜壓會(huì)逐漸再次轉(zhuǎn)變?yōu)閯?dòng)壓,使壁面的切向速度最終達(dá)到一個(gè)最大值。由于稀薄效應(yīng)的影響,最大無(wú)量綱滑移速度會(huì)隨著克努森數(shù)的增加而增大,當(dāng)克努森數(shù)達(dá)到0.1時(shí),最大無(wú)量綱滑移速度可達(dá)到自由來(lái)流速度的16%。
圖10 無(wú)量綱滑移速度沿圓柱體壁面的分布
圖11描述了在滑移流區(qū)雷諾數(shù)為1時(shí)壓力系數(shù)沿圓柱體壁面的分布規(guī)律。當(dāng)克努森數(shù)保持不變時(shí),壓力系數(shù)最大值位于圓柱體前止滯點(diǎn)處,且沿著圓柱體壁面逐漸減小。然而,由于滑移效應(yīng)的影響,圓柱體前止滯點(diǎn)處的壓力系數(shù)會(huì)隨著克努森數(shù)的增加而逐漸減小,不過(guò)克努森數(shù)對(duì)壓力系數(shù)的影響其實(shí)并不大。
圖12 表面摩擦系數(shù)沿圓柱體壁面的分布
圖12描述了在滑移流區(qū)雷諾數(shù)為1時(shí)表面摩擦系數(shù)沿圓柱體壁面的分布規(guī)律。從圖中可知,當(dāng)克努森數(shù)一定時(shí),由于雷諾數(shù)較低,黏性力起主導(dǎo)作用,從而抑制邊界層分離,使圓柱體前止滯點(diǎn)與后止滯點(diǎn)速度為0,故表面摩擦系數(shù)沿著圓柱體壁面先從0增加至最大值,然后逐漸減小至0。由于稀薄效應(yīng)影響,最大表面摩擦系數(shù)隨著克努森數(shù)增加而逐漸減小,當(dāng)克努森數(shù)達(dá)到0.1時(shí),最大表面摩擦系數(shù)比連續(xù)流區(qū)最大表面摩擦系數(shù)減小了38%。
針對(duì)滑移流區(qū)氣體外掠圓柱體問(wèn)題,構(gòu)建了一套預(yù)測(cè)低壓氣體流動(dòng)與傳熱的數(shù)值模型,考慮了軸向溫度梯度引起的熱蠕動(dòng)影響和壁面曲率的影響,通過(guò)與實(shí)驗(yàn)數(shù)據(jù)對(duì)比,得出以下結(jié)論:
(1)當(dāng)氣體流態(tài)處于連續(xù)流區(qū)時(shí),滑移邊界修正模型的計(jì)算結(jié)果與無(wú)滑移邊界的計(jì)算結(jié)果相近;
(2)當(dāng)氣體流態(tài)處于滑移流區(qū)時(shí),預(yù)測(cè)外掠圓柱體傳熱的傳統(tǒng)經(jīng)驗(yàn)關(guān)聯(lián)式已不再適用,此時(shí)采用滑移邊界修正的數(shù)值模型,可實(shí)現(xiàn)對(duì)稀薄氣體滑移流區(qū)流動(dòng)與換熱的精確預(yù)測(cè),且數(shù)值計(jì)算結(jié)果明顯優(yōu)于傳統(tǒng)經(jīng)驗(yàn)關(guān)聯(lián)式預(yù)測(cè)結(jié)果;
(3)在滑移流區(qū),稀薄效應(yīng)對(duì)低壓氣體外掠圓柱體的流體動(dòng)力學(xué)特性影響顯著,隨著克努森數(shù)增加,最大無(wú)量綱滑移速度呈線性增加,而最大表面摩擦系數(shù)呈現(xiàn)線性減小。