張俊韜，侯中喜，柳兆偉，麻震宇

臨近空間飛艇內(nèi)部自然對(duì)流的流場(chǎng)特征仿真*

張俊韜，侯中喜，柳兆偉，麻震宇

(國防科技大學(xué) 航天科學(xué)與工程學(xué)院， 湖南 長(zhǎng)沙 410073)

針對(duì)臨近空間飛艇，對(duì)其內(nèi)部氦氣的自然對(duì)流特征進(jìn)行了數(shù)值模擬。采用計(jì)算流體力學(xué)方法，并利用自編的用戶自定義函數(shù)將外部的溫度邊界條件導(dǎo)入壁面網(wǎng)格。在不同的穩(wěn)態(tài)條件下，通過對(duì)內(nèi)部氦氣壓力、溫度、速度的分布等流場(chǎng)特征參數(shù)的研究，分析了臨近空間飛艇內(nèi)部氣體自然對(duì)流的運(yùn)動(dòng)特性及其影響規(guī)律，并對(duì)自然對(duì)流非穩(wěn)態(tài)變化過程進(jìn)行了初步的探索。仿真結(jié)果表明，在臨近空間環(huán)境下，飛艇內(nèi)部氦氣的自然對(duì)流，對(duì)于內(nèi)部氦氣自身的熱交換具有一定程度的促進(jìn)作用，而對(duì)蒙皮受力和結(jié)構(gòu)安全性影響很小。

自然對(duì)流；臨近空間飛艇；數(shù)值仿真；流場(chǎng)特征

(College of Aerospace Science and Engineering, National University of Defense Technology, Changsha 410073, China)

自20世紀(jì)70年代以來，以美國為首的西方國家開展了臨近空間飛艇的相關(guān)研究，分析了長(zhǎng)期駐空等難題[1]。由于臨近空間熱環(huán)境中的諸多熱源和影響因素，飛艇表面蒙皮的溫度呈高低溫區(qū)域復(fù)雜分布且隨時(shí)間變化，同時(shí)飛艇內(nèi)部的浮升氣體(氦氣)在受熱不均的條件下發(fā)生密度的變化，進(jìn)而由于重力作用而產(chǎn)生氣體的流動(dòng)，即自然對(duì)流的現(xiàn)象。而自然對(duì)流形成的流場(chǎng)特征和內(nèi)部氣體溫度、壓力等參數(shù)的分布變化反之也會(huì)影響表面蒙皮的溫度、壓力的分布，因此這是一個(gè)耦合影響的過程。

在之前眾多有代表性的相關(guān)研究中，Lambert和Smith等[2-3]指出，對(duì)臨近空間飛艇進(jìn)行熱分析、研究熱模型和仿真，并進(jìn)行熱特性試驗(yàn)是十分有意義的；Harada等[4]對(duì)低空飛艇模型進(jìn)行了相關(guān)的試驗(yàn)研究，測(cè)得飛艇駐空期間蒙皮周向的溫度分布；Lee等[5]對(duì)飛艇由于內(nèi)部氣體運(yùn)動(dòng)而導(dǎo)致的浮力和質(zhì)心變化進(jìn)行了低空試驗(yàn)，計(jì)算了蒙皮的三維穩(wěn)態(tài)溫度分布。對(duì)自然對(duì)流而言，研究多集中于較為規(guī)則的簡(jiǎn)單外形上[6-7]，如圓球或者圓柱。Marc等[8]采用旋轉(zhuǎn)加熱的方法對(duì)圓柱形飛艇內(nèi)自然對(duì)流進(jìn)行試驗(yàn)研究和仿真分析。對(duì)于實(shí)際飛艇外形，其內(nèi)部自然對(duì)流的研究還不是很透徹。國內(nèi)方面，方賢德等[9]用兩節(jié)點(diǎn)穩(wěn)態(tài)模型計(jì)算了平流層飛艇平衡溫度，得到表面蒙皮和內(nèi)部氣體的平均溫度。然而，由于飛艇表面積很大，用平均溫度代替表面蒙皮溫度會(huì)造成較大誤差，因此Li等[10]進(jìn)一步分析了半硬式飛艇蒙皮的溫度分布，建立了相應(yīng)的分析模型。而本文的著眼點(diǎn)在于，將上述方法得到的溫度分布作為初始邊界條件，并采用計(jì)算流體力學(xué)(Computational Fluid Dynamics, CFD)的方法，直接對(duì)整個(gè)內(nèi)部流場(chǎng)進(jìn)行數(shù)值計(jì)算。

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

1.1 物理問題描述

為簡(jiǎn)化研究對(duì)象，假設(shè)飛艇只有一個(gè)主氣囊，內(nèi)部充滿一定質(zhì)量的氦氣(不考慮副氣囊以及閥門放氣等因素)，同時(shí)在任何時(shí)刻飛艇均處在最大體積狀態(tài)下(浮力不變從而不會(huì)發(fā)生上下運(yùn)動(dòng)，保持平飛巡航)。則飛艇在臨近空間巡航駐留的狀態(tài)下，外部的熱條件隨時(shí)間不斷變化，導(dǎo)致內(nèi)部氦氣產(chǎn)生自然對(duì)流現(xiàn)象。在這個(gè)過程中，自然對(duì)流是一個(gè)不斷變化發(fā)展的非穩(wěn)態(tài)過程，然而由于一天之中在相對(duì)較短的時(shí)間尺度范圍內(nèi)(1到2個(gè)小時(shí)之內(nèi))，熱環(huán)境的變化還是相對(duì)較慢的，因此可以把一整天的自然對(duì)流變化過程離散成若干近似穩(wěn)態(tài)的過程來研究，并通過分析這些不同時(shí)刻的自然對(duì)流狀態(tài)來分析其變化規(guī)律。

1.2 數(shù)值方法

.2.1 控制方程

適用于該問題的控制方程為三維定?？蓧篘-S方程：

? 連續(xù)方程：

div(ρU)=0

(1)

? 動(dòng)量方程：

(2)

? 能量方程：

ρ(U·)2T

(3)

其中：U是流體速度矢量；ρ為密度；div代表散度運(yùn)算符；代表Hamilton算子；F為徹體力矢量；μ為動(dòng)力粘度(為常數(shù))；P為壓力；grad代表梯度算子；T為溫度；k為熱導(dǎo)率；cp為定壓比熱。

以上方程組為數(shù)值計(jì)算的控制方程，通過商業(yè)軟件FLUENT的求解器進(jìn)行求解。

.2.2 邊界條件

建立蒙皮表面的網(wǎng)格化單元，如圖1所示；通過文獻(xiàn)[11]中的方法建立蒙皮單元的熱模型，并進(jìn)行迭代計(jì)算即可得到任意時(shí)刻各網(wǎng)格單元的溫度。

圖1 表面網(wǎng)格示意圖Fig.1 Surface mesh sketch map

將上述每個(gè)網(wǎng)格單元所獲得的溫度值，通過FLUENT的用戶自定義函數(shù)(User-Defined Function, UDF)進(jìn)行加載，從而得到了這一時(shí)刻的表面蒙皮的溫度分布(以下所有溫度的單位均為K)，以此作為CFD方法的邊界條件，其效果如圖2所示。

圖2 壁面溫度分布Fig.2 Wall temperature distribution

.2.3 確定流態(tài)

在純粹自然對(duì)流中，浮力所引致的流動(dòng)強(qiáng)度可用瑞利數(shù)Ra來判定：

(4)

式中：ΔT為溫差；L為特征長(zhǎng)度；β為熱膨脹系數(shù)。

(5)

α為熱擴(kuò)散率：

(6)

根據(jù)Ra的值所在的范圍，即可以判定流動(dòng)是層流還是湍流。其中當(dāng)Ra<108，自然對(duì)流基本為層流；當(dāng)108

表1 飛艇參數(shù)

2 計(jì)算結(jié)果與分析

2.1 流動(dòng)特征及時(shí)間變化的影響

在文獻(xiàn)[12]對(duì)單一時(shí)刻的自然對(duì)流狀態(tài)進(jìn)行分析的基礎(chǔ)上，進(jìn)行多時(shí)刻點(diǎn)的仿真計(jì)算和分析，具體計(jì)算參數(shù)依然如表1所示。此外，計(jì)算條件還包括：地點(diǎn)為北緯40°附近某地上空20km處，時(shí)間為2013年7月5日一整天的若干個(gè)典型時(shí)刻的狀態(tài)(2小時(shí)間隔，即4:00，6:00，……，18:00，20:00，其中6:00前與18:00后的狀態(tài)是一致的，因?yàn)檫@些時(shí)間段沒有太陽輻射)。針對(duì)以上情況的自然對(duì)流狀態(tài)進(jìn)行仿真計(jì)算，仿真過程中，設(shè)置的飛艇飛行朝向?yàn)檎戏较颉?/p>

2.1.1 蒙皮溫度極大值的變化

由圖3可以看到，蒙皮溫度的極大值曲線以正午12:00為中心呈近似對(duì)稱結(jié)構(gòu)分布。日出前和日落后的壁面最高溫度在224K附近，同時(shí)，在有日照的日間時(shí)間段，壁面溫度的極大值并不在正午12:00，而是在接近早晨8:00和傍晚18:00附近，約為254K。

圖3 蒙皮溫度極大值的變化曲線Fig.3 Curve of maximum temperature of skin

分析這種現(xiàn)象的主要原因，具體來講：對(duì)于某個(gè)蒙皮單元的熱平衡方程：

(7)

式中：δ為厚度，ρ為密度，Cp為比熱，τ為時(shí)間；hin是內(nèi)部自然對(duì)流換熱系數(shù)；hex為外部強(qiáng)迫對(duì)流換熱系數(shù)；Tg是內(nèi)部氦氣在這一瞬時(shí)的溫度；Ta是來流溫度(恒定值)；T是單元溫度；qt是單元外部吸收的總的輻射熱流，包括直射太陽輻射、反射和大氣散射以及地表紅外輻射；εex是材料的表面發(fā)射率。

由式(7)可知，一個(gè)最高溫度的局部單元需要獲得的總熱量最高。此單元在各個(gè)時(shí)刻所獲得的太陽輻射、散射、反射輻射以及紅外輻射等熱量基本相同(忽略不同時(shí)刻透射率不同的影響)，鑒于外界氣體溫度恒溫，則外部的對(duì)流換熱因素影響一致，因此影響溫度極大值變化的因素就主要集中在內(nèi)部的對(duì)流換熱項(xiàng)以及蒙皮自身向外的紅外輻射項(xiàng)。就對(duì)流換熱項(xiàng)而言，這個(gè)溫度最高的單元向內(nèi)部氣體進(jìn)行對(duì)流換熱量最少的時(shí)刻就是內(nèi)部氣體溫度最高的時(shí)刻，這時(shí)這個(gè)單元的溫度才會(huì)達(dá)到極大值。

圖4 內(nèi)部氦氣平均溫度的變化曲線Fig.4 Curve of average temperature of internal helium

同時(shí)，內(nèi)部氣體的平均溫度變化曲線如圖4所示，正午12:00時(shí)刻內(nèi)部氣體的平均溫度并不是極大值，而是一個(gè)低谷，這是因?yàn)閮?nèi)部氣體的平均溫度的高低主要取決于太陽高度角和方位角以及飛艇朝向的綜合效果而導(dǎo)致的太陽輻射在飛艇上的總面積，面積越大總熱量就越多，而前面已經(jīng)提到，算例中飛艇的方向朝南，則可以推出剛剛?cè)粘龊涂煲章鋾r(shí)的輻射面積最大，因而氦氣平均溫度的極大值就出現(xiàn)在早晨和傍晚附近。與此同時(shí)，蒙皮溫度越高，向外輻射的熱量也越多，反而會(huì)阻礙溫度的進(jìn)一步上升，因此蒙皮溫度極大值的分布是以上影響因素的綜合效果。

2.1.2 內(nèi)部氣體的溫度分布

如圖5和圖6所示，雖然蒙皮邊界上的溫度分布不同，但內(nèi)部氣體的溫度分布狀態(tài)在不同的時(shí)間點(diǎn)是類似的(比較10:00和12:00，其他時(shí)間狀態(tài)未列出，但也類似)，即由于自然對(duì)流的存在使得處于中間部分的氣體溫度分布出現(xiàn)一定程度的分層現(xiàn)象，靠上部的氣體更加接近高溫而靠下部的氣體更加接近較低的低溫，并且無論壁面的熱邊界條件呈何種角度，內(nèi)部氦氣的溫度分布基本都在Z軸方向且由上到下逐漸分層過渡，溫度由高到低。至此，說明蒙皮溫度的角度分布變化對(duì)內(nèi)部氣體的溫度分布規(guī)律沒有影響，只影響內(nèi)部氣體的溫度具體數(shù)值的高低。

圖5 在10:00的溫度分布Fig.5 Temperature distribution at 10:00

圖6 在12:00的溫度分布Fig.6 Temperature distribution at 12:00

2.1.3 壓力分布的不均勻性

壓力不均勻性是指在某一時(shí)刻，內(nèi)部氣體壓力的最小值與最大值的差值。雖然飛艇所經(jīng)歷的熱環(huán)境在時(shí)刻變化，然而，由于自然對(duì)流所引起的腔內(nèi)壓力不均勻性的量級(jí)非常小，只有帕(Pa)以下的量級(jí)，如圖7所示。

圖7 在10:00和12:00的壓力分布Fig.7 Pressure distribution at 10:00 and 12:00

與此同時(shí)，在各個(gè)時(shí)刻，如圖7所示，內(nèi)部氣體的壓力分布形式基本上保持一致(依然以10:00和12:00比較為例)，呈現(xiàn)出兩端高中間低的分布狀態(tài)，只是不同時(shí)刻對(duì)應(yīng)的具體壓力值略有不同。具體來講，無論壁面的熱邊界條件呈何種角度，壓力分布都在沿Z軸(重力所在的豎直方向)由中間向上下升高，但這種分布的不均勻性很小，因此可以說明自然對(duì)流所引起的壓力分布不均基本不會(huì)對(duì)蒙皮的受力以及結(jié)構(gòu)特性造成嚴(yán)重的影響。

2.1.4 速度最大值的變化

如圖8所示，從一晝夜速度最大值的變化曲線來看，整個(gè)過程自然對(duì)流引起的氦氣內(nèi)部流動(dòng)的最大速度不超過2m/s。

圖8 速度最大值的變化曲線Fig.8 Curve of maximum velocity

而且，如圖9所示，最大速度集中分布在壁面附近的周向和縱軸向的環(huán)流，而其余大部分內(nèi)部區(qū)域的流速小于1m/s。因此可以認(rèn)為，飛艇駐留段內(nèi)部自然對(duì)流所導(dǎo)致的內(nèi)部氦氣流動(dòng)從速度的量級(jí)上來看很小，流動(dòng)比較微弱。

圖9 在12:00的速度分布切片圖Fig.9 Slices of velocity distribution at 12:00

進(jìn)一步分析，與蒙皮溫度的極大值類似，氦氣的最大速度同樣出現(xiàn)在接近早晨8:00和傍晚18:00附近，最小速度出現(xiàn)在接近上午10:00和下午14:00附近，這種現(xiàn)象說明：在早晨8:00和傍晚18:00附近，氦氣的自然對(duì)流較為強(qiáng)烈；而在上午10:00和下午14:00附近，氦氣的自然對(duì)流相對(duì)微弱一些。

2.2 內(nèi)部自然對(duì)流的非穩(wěn)態(tài)運(yùn)動(dòng)過程研究

設(shè)定邊界條件從早上6:00時(shí)刻為初值，下一時(shí)刻變?yōu)槿粘龊蟮纳衔?:00的情況。溫度邊界條件的變化如圖10所示。

圖10 從早上6:00至8:00的壁面溫度變換Fig.10 Change of the wall temperature from 6:00 a.m. to 8:00 a.m.

以10s為時(shí)間間隔，仿真模擬了15min的內(nèi)部自然對(duì)流非穩(wěn)態(tài)流動(dòng)過程。這里需要說明的是，實(shí)際情況下，8:00的邊界條件是從6:00緩慢變化而來的，而這里為簡(jiǎn)化研究問題，直接將8：00的蒙皮溫度分布加載到在6:00已經(jīng)達(dá)到穩(wěn)定狀態(tài)的自然對(duì)流的邊界條件上去，即認(rèn)為邊界條件是由一個(gè)狀態(tài)瞬時(shí)變化到了下一個(gè)狀態(tài)。這樣，雖然所設(shè)置的條件與實(shí)際情況有所出入，但是最初和最終的狀態(tài)是一致的，即中間自然對(duì)流過程只是過渡的時(shí)間有所不同，因此不影響變化過程的規(guī)律研究。以下主要從流動(dòng)現(xiàn)象來分析溫度分布的變化規(guī)律。

圖11 Y-Z平面非穩(wěn)態(tài)溫度變化過程Fig.11 Process of unstable state temperature at Y-Z plane

如圖11(a)～(f)所示，顏色越淺說明溫度越高。取一個(gè)Y-Z平面，在初始狀態(tài)下，顯然由于熱平衡，內(nèi)部氣體仍然存在穩(wěn)定的自然對(duì)流狀態(tài)，同時(shí)氣體的整體溫度較低(見圖11(a))。當(dāng)出現(xiàn)太陽輻照后，外界熱邊界條件瞬間變化，如圖10所示，熱壁面在右側(cè)，冷壁面在左上側(cè)，導(dǎo)致內(nèi)部氦氣在最貼近這兩處的地方溫度首先發(fā)生變化，分別趨向于高溫和低溫(見圖11(b)～(c))，進(jìn)而高溫氣體由于密度較低則沿壁面上升，低溫氣體密度較大沿壁面下降，此即為自然對(duì)流的起始動(dòng)因，同時(shí)由于流動(dòng)的慣性，在較早的時(shí)間內(nèi)，高低溫氣體的分布與垂直方向存在一小的角度(見圖11(d))，經(jīng)過一段時(shí)間的過渡后，氣體溫度分布趨于均勻，同時(shí)氦氣的整體溫度也有明顯的提高(見圖11(e)～(f))。

圖12 X-Z平面非穩(wěn)態(tài)溫度變化過程Fig.12 Process of unstable state temperature at X-Z plane

如圖12的溫度分布變化過程所示，同樣顏色越淺代表溫度越高。取一個(gè)X-Z平面，情況與Y-Z平面基本類似，而且溫度變化的情況要更簡(jiǎn)單直觀。由于在此平面內(nèi)，溫度的邊界條件即為上高下低，則上下兩端的氣體正好位于這兩個(gè)區(qū)域，溫度分層隨時(shí)間逐漸明顯變化(由(a)～(d)的變化)，同時(shí)整體溫度升高。

3 結(jié)論

在某一時(shí)刻的穩(wěn)態(tài)自然對(duì)流條件下，內(nèi)部氦氣的壓力呈現(xiàn)出有規(guī)律的分層現(xiàn)象，同時(shí)，其不均勻性很小。這表明自然對(duì)流所導(dǎo)致的內(nèi)部壓力不均勻性可以忽略不計(jì)；溫度則呈現(xiàn)出從上到下由高溫到低溫逐漸過渡的情況；內(nèi)部氣體的流動(dòng)速度較低，最大值不超過2m/s，其中大部分區(qū)域小于1m/s，與飛艇百米以上的尺度來衡量，可以說是比較緩慢和微弱的，不會(huì)對(duì)蒙皮受力和結(jié)構(gòu)安全性造成影響。

對(duì)于不同時(shí)刻，盡管外蒙皮的溫度分布不同，但其內(nèi)部氣體由于自然對(duì)流所引起的溫度、壓力、速度的分布規(guī)律沒有明顯的變化，說明在一晝夜內(nèi)，自然對(duì)流的總體特性基本保持不變，即自然對(duì)

流使內(nèi)部氦氣的溫度呈分層狀態(tài)，而壓力的不均勻性和內(nèi)部流場(chǎng)的速度量級(jí)均較小。

對(duì)于非穩(wěn)態(tài)過程，通過對(duì)不同時(shí)刻流場(chǎng)狀態(tài)的分析，說明變化的溫度邊界條件是導(dǎo)致自然對(duì)流狀態(tài)發(fā)生改變的直接動(dòng)因。

References)

[1] Stefan K. Thermal effects on a high altitude airship [C]//Proceedings of 5th Lighter-Than-Air Systems Conference, 1983.

[2] Lambert C. Developing architectures for high altitude airship[C]//Proceedings of AIAA′s 3rd Annual Aviation Technology, Integration, and Operations (ATIO) Forum,AIAA-2003-6781, 2003.

[3] Smith M S, Rainwater E L. Applications of scientific ballooning technology to high altitude airships [C]// Proceedings of AIAA′s 3rd Annual Aviation Technology, Integration, and Operations (ATIO) Forum, AIAA-2003-6711, 2003.

[4] Harada K, Eguchi K, Sano M, et al. Experimental study of thermal modeling for stratospheric platform airship[C]//Proceedings of AIAA′s 3rd Annual Aviation Technology, Integration, and Operations (ATIO) Forum,AIAA-2003-6833, 2003.

[5] Lee Y G, Kim D M, Yeom C H. Development of Korean high altitude platform systems[J].International Journal of Wireless Information Networks, 2006, 13(1): 31-42.

[6] Rodger E F. Balloon ascent: 3-D simulation tool for the ascent and float of high-altitude balloons[C]//Proceedings of AIAA 5th ATIO and 16th Lighter-Than-Air Systems Technology and Balloon Systems Conferences, AIAA-2005-7412, 2005.[7] Das T, Mukherjee R, Cameron J. Optimal trajectory planning for hot-air balloons in linear wind fields [J]. Journal of Guidance, Control, and Dynamics, 2003, 26(3): 416-424.

[8] Marc H, Bernhard W, Von Wolfersdorf J. Natural convection inside airship [C]//Proceedings of 9th AIAA/ASME Joint Thermophysics and Heat Transfer Conference, AIAA-2006-3798, 2006.

[9] 方賢德, 王偉志, 李小建. 平流層飛艇熱仿真初步探討 [J]. 航天返回與遙感, 2007, 28(2): 5-9.

FANG Xiande, WANG Zhiwei, LI Xiaojian. A study of thermal simulation of stratospheric airships [J]. Spacecraft Recover & Remote Sensing, 2007, 28(2)：5-9. (in Chinese)

[10] Li X J, Fang X D, Dai Q M, et al. Modeling and analysis of floating performances of stratospheric semi-rigid airships [J]. Advances in Space Research, 2012, 50(7): 881-890.

[11] 夏新林, 李德富, 楊小川. 復(fù)合熱條件下橢球形封閉腔內(nèi)低壓氣體的自然對(duì)流[J]. 航空學(xué)報(bào), 2010, 31(3): 453-458.

XIA Xinlin, LI Defu, YANG Xiaochuan. Natural convection of low pressure gas in ellipsoidal enclosure induced by combined thermal conditions [J]. Acta Aeronautica Et Astronautica Sinica, 2010, 31(3): 453-458. (in Chinese)

[12] 張俊韜,侯中喜,柳兆偉.平流層飛艇內(nèi)部氦氣的自然對(duì)流研究 [C].第二屆高分辨率對(duì)地觀測(cè)學(xué)術(shù)年會(huì), 長(zhǎng)沙,2013.

ZHANG Juntao, HOU Zhongxi, LIU Zhaowei. Research on natural convection of helium inside stratospheric airship [C]. 2nd China High Resolution Earth Observation Conference, Changsha, 2013. (in Chinese)

Simulation on flow field characteristics of natural convection inside near space airship

ZHANG Juntao, HOU Zhongxi, LIU Zhaowei, MA Zhenyu

Aiming at the near space airship, the natural convection characteristics of the helium inside an airship were simulated numerically. On the basis of the computational fluid dynamics method and user-defined function, external temperature boundary condition was imported into the wall grids. Under conditions of different steady states, the motion characteristics and influence regulations of gas natural convection inside the airship were analyzed by researching the characteristic parameters of flow field, such as distribution of pressure, temperature, speed. Besides, the unstable change process of natural convection was also studied preliminarily. The simulation results show that natural convection of helium inside a near space airship has some active influence on heat exchange of inside helium, but the influence is very weak on stressed skin and structural safety.

natural convection; near space airship; numerical simulation; flow field characteristics

2015-04-10

航空科學(xué)基金資助項(xiàng)目(20145788006)

張俊韜(1986—)，男，天津人，博士研究生，E-mail：zzt136@163.com；侯中喜(通信作者)，男，教授，博士，博士生導(dǎo)師，E-mail：hzx@sina.com

10.11887/j.cn.201504002

http://journal.nudt.edu.cn

O351.2

A

1001-2486(2015)04-009-06