考慮蒙皮透射率的飛艇熱力學模型及其熱特性

程 晨， 王曉亮

(上海交通大學 航空航天學院， 上海 200240)

平流層飛艇在工作過程中，外界熱環境復雜且具有時變性.平流層空氣稀薄，因此飛艇一般體積都很龐大，在駐空飛行過程中，其熱特性受自身材料熱特性、外部環境熱源及囊體內外換熱方式的影響極大.強烈的太陽輻射和微弱的對流換熱導致飛艇蒙皮及內部氣囊內溫度場晝夜溫差很大，從而影響飛艇的浮力和蒙皮的應力水平，改變其駐空高度、飛行姿態及軌跡，局部過熱和應力增大也可能會導致蒙皮的膨脹變形乃至破壞.因此建立準確的飛艇熱特性模型，掌握及預測飛艇蒙皮及內部填充氣體的溫度，對于飛艇設計及控制有著重要的意義.

20世紀80年代，國外許多學者就開始研究飛艇熱特性.近年來，國內學者對平流層飛艇熱特性模型的研究也在持續發展，取得了一定的成果[1].文獻[2]建立了瞬態運動方程和傳熱模型.文獻[3]用Fluent軟件對飛艇進行了仿真分析.文獻[4]研究發現相同條件下南瓜形氣球的氦氣溫度及蒙皮平均溫度均高于球形氣球.文獻[5]詳細研究了飛艇晴空時的太陽輻射、長波輻射以及強迫換熱模型.文獻[6]分析了駐空期間太陽電池等效面積熱阻、轉換效率及鋪裝面積對飛艇熱溫度晝夜變化規律的影響.

目前飛艇熱特性的主要研究方法可以分成兩大類：集總參數法和分布參數法.集總參數法是將飛艇看成一個或幾個溫度均勻分布的部分進行研究，飛艇蒙皮及內部填充氣體的熱特性僅僅隨著時刻的變化而變化，與具體的坐標位置無關.分布參數法則首先需要將蒙皮表面網格化，考慮熱源強度與蒙皮空間位置的相關性，再分別計算每一個蒙皮單元的熱特性.集總參數法的計算結果較為粗略，但是計算速度較快，適用于計算上升、下降時的飛艇狀態[7-9]，分布參數法相較于集總參數法計算量更大，但是計算結果更為準確，更適用于計算飛艇駐空狀態時的熱特性[10-12].

在現有的研究中，絕大多數學者在分析內部填充氣體的熱特性時，只計算了氣體與蒙皮內表面之間的自然對流換熱，未考慮實際情況下飛艇蒙皮材料的透射特性以及內部填充氣體對于輻射的吸收率與發射率，忽略了各類輻射對內部填充氣體熱特性的直接影響，因此傳統模型只適用于計算蒙皮為完全不透明材料的飛艇熱力學特性.本文在建立飛艇熱力學特性模型時，考慮了蒙皮的透射特性，給出了新的蒙皮單元遮擋系數的計算方法，同時也研究了外部輻射熱源對于飛艇內部填充氣體熱特性的影響，并補充計算了蒙皮與內部填充氣體之間的輻射換熱.本文采用有限拆分法來實現對模型的數值計算，即在分析內部填充氣體熱特性時采用集總參數法，在分析表面熱特性時采用分布參數法，在保證計算精度的前提下，提高了計算效率.

1 考慮蒙皮透射率的飛艇熱力學模型

影響平流層飛艇熱力學特性的熱源及換熱方式可分為太陽短波輻射、太陽長波輻射以及對流換熱.

1.1 飛艇熱特性模型耦合機制

目前，大部分學者在研究平流層飛艇熱特性模型時，對飛艇蒙皮材料以及內部填充氣體的特性做了以下簡化處理：

(1) 忽略了蒙皮對輻射的透射特性，外部輻射無法透射穿過蒙皮表面，只會作用在不被遮擋的蒙皮單元表面，且不會影響飛艇內部填充氣體的換熱.

(2) 將飛艇內部填充氣體看作完全透明氣體，不考慮內部氣體對輻射的吸收及發射.

未考慮蒙皮透射特性的傳統飛艇熱環境模型及換熱機制如圖1和2所示.

圖1 傳統平流層飛艇熱環境

圖2 傳統飛艇熱模型換熱機制

然而，在實際計算中，根據蒙皮材料特性的不同，外部熱源輻射會或多或少地透射穿過蒙皮，作用在內部填充氣體以及對面蒙皮上.考慮蒙皮透射特性時，外部輻射經過飛艇的路徑示意圖如圖3所示.

圖3 考慮蒙皮透射特性的外部輻射路徑示意圖

本文在傳統飛艇熱特性模型的基礎上，補充考慮了蒙皮的透射特性，研究了外部輻射熱源對于飛艇內部填充氣體熱特性的影響，并計算了蒙皮與內部填充氣體之間的輻射換熱.本文提出的考慮蒙皮透射特性的新飛艇換熱機制如圖4所示.

圖4 新飛艇熱模型傳熱機制

1.2 飛艇蒙皮熱力學方程

由于平流層空氣稀薄，飛艇體積龐大，且外部熱源在蒙皮上分布不均勻，因此同一時刻蒙皮不同處的溫度差別很大，在研究飛艇蒙皮的受熱情況時，不宜將蒙皮看作一個整體進行分析，本文采取分布參數法，將蒙皮表面劃分成若干個三角形單元進行熱特性的分析.為了在保證結果準確性的前提下提高計算效率，本文在計算飛艇傳熱過程時采用了如下的基本假設：

(1) 由于浮空器蒙皮只能吸收一部分熱輻射而反射其余部分，所以可將蒙皮近似為灰體, 平衡狀態下蒙皮材料的紅外發射率與其紅外吸收率相等.

(2) 由于對于一般的漫灰體，在工程上的溫度變化范圍內，其紅外吸收率以及紅外發射率與波長和溫度無關，所以可假設蒙皮表面紅外發射率不隨波長變化；由于短波波長范圍很小，所以也可以假設蒙皮表面短波輻射吸收率不隨波長變化.

(3) 由于蒙皮單元的大小相較于浮空器非常小，每個單元上的溫度梯度極小，所以蒙皮各單元內的溫度可近似看成均勻分布.

(4) 由于本文設計的浮空器為內部被浮升氣體填滿的橢球體，所以可假設蒙皮內表面所有單元均為非凹面.

(5) 由于蒙皮厚度極小且導熱率低，所以沿厚度方向的熱傳導以及不同蒙皮單元之間的熱傳導可以忽略不計.

考慮透射穿過蒙皮的輻射在飛艇內部反復被蒙皮及內部填充氣體反射、透射及吸收的過程，為了保證計算的準確性，本文使用內部填充氣體的等效物理參數來計算此換熱過程中的熱特性.根據文獻[13]給出的等效參數計算公式，可得內部填充氣體等效太陽輻射吸收率αgaseff和等效紅外發射率εgaseff，計算方法如下：

(1)

(2)

式中：αgas、εgas分別為內部填充氣體常態時的太陽輻射吸收率以及紅外發射率；τskinsol、τskinIR分別為蒙皮對太陽輻射以及紅外輻射的透射率；rskinsol、rskinIR分別為蒙皮對太陽輻射以及紅外輻射的反射率.

考慮蒙皮透射率的飛艇表面蒙皮單元i的瞬態溫度控制方程為

QIRSky,i+QIRin,i+QIRgas,i+Qci,i+Qco,i

(3)

式中：cskin,i、mskin,i及Tskin,i分別為蒙皮單元i的比熱容、質量及溫度；式右端為單位時間內蒙皮單元i上吸收的所有熱量的總和，Qs,i、Qsunsa,i及Qg,i分別為太陽短波輻射中直接太陽輻射、天空散射輻射以及地面及云層反射輻射在單位時間內被飛艇蒙皮單元i吸收的熱量；QIREarth,i、QIRSky,i、QIRin,i及QIRgas,i分別為長波輻射中蒙皮外表面吸收的地面長波輻射和天空長波輻射、飛艇蒙皮內表面單元之間輻射換熱以及飛艇蒙皮與內部填充氣體之間的輻射換熱在單位時間內被飛艇蒙皮單元i吸收的熱量；Qci,i、Qco,i分別為蒙皮內表面自然對流換熱以及外表面混合對流換熱過程中在單位時間內被飛艇蒙皮單元i吸收的熱量.

飛艇表面蒙皮單元i在單位時間內吸收的直接太陽輻射量計算式為

Qs,i=λ1αskinIsAskin,icosβ1

(4)

式中：λ1為直接太陽輻射的遮擋系數，考慮蒙皮的透射率，當直接太陽輻射直接作用在蒙皮單元i上時，λ1=1，反之λ1=τskinsol(1-αgaseff)；αskin為蒙皮對太陽直射輻射的有效吸收率；Is為直接太陽輻射強度[14]，與大氣透射率以及大氣層外邊界處的太陽輻射強度相關；Askin,i是蒙皮單元i的面積;β1為蒙皮單元法向與太陽光線的夾角.

天空散射輻射、地面及云層反射輻射、周圍大氣長波輻射以及地面長波輻射這4個熱源的輻射方向不隨時刻改變(天空散射輻射以及周圍大氣長波輻射的方向始終豎直向下，而地面及云層反射輻射以及地面長波輻射的方向始終豎直向上)，因此在計算蒙皮單元吸收的熱量時，蒙皮遮擋系數有一定的共性.可首先根據蒙皮單元i的法向與豎直向下方向的夾角β2判斷出半透明蒙皮材料飛艇表面蒙皮單元i所在位置，然后按照對應表面的熱量計算式計算蒙皮單元吸收的熱量.

(1) 若飛艇表面蒙皮單元在飛艇上表面，則蒙皮單元i在單位時間內吸收的天空散射輻射、周圍大氣長波輻射、地面及云層反射輻射以及地面長波輻射計算式如下[2,5,15-16]：

(5)

(6)

(7)

QIREarth,i=τskinsol(1-εgaseff)εskin_ex(IIREarth-

(8)

式中：Isunsa、Ig分別為天空散射輻射強度以及地面及云層反射輻射強度；σ為玻爾茲曼常數；εskin_ex為蒙皮外表面紅外發射率；IIREarth和IIRSky分別為等效地球長波輻射強度以及等效大氣長波輻射強度.

(2) 若飛艇表面蒙皮單元在飛艇下表面，則計算式分別為

Qsunsa,i=τskinsol(1-αgaseff)αskinIsunsa×

(9)

QIRSky,i=τskinsol(1-εgaseff)εskin_ex(IIRSky-

(10)

(11)

(12)

飛艇蒙皮內表面間的輻射換熱的計算式為

(13)

式中：εskin_in為蒙皮內表面的紅外發射率; 第i個單元發射的熱輻射Ji可以表示為

(14)

i=1,2,…,N

Xi,j是i單元相對于j單元的角系數[17].

蒙皮內表面單元i單位時間內吸收的內部填充氣體熱輻射量計算式為

QIRgas,i=εskin_inIIRgasAskin,i

(15)

式中：IIRgas為內部填充氣體發射的紅外輻射強度.在計算內表面蒙皮與內部填充氣體之間的輻射換熱時，將氣體及蒙皮單元等效成灰體，根據灰體輻射的四次方定律，可以得到內部填充氣體發射的輻射強度為

(16)

式中:Tgas為內部填充氣體的溫度.

單位時間內，蒙皮單元i內表面吸收的自然對流換熱量以及外表面吸收的混合對流換熱量計算式分別為

Qci,i=hin(Tgas-Tskin,i)Askin,i

(17)

Qco,i=hex(Tatm-Tskin,i)Askin,i

(18)

式中：Tatm為外部大氣的環境溫度；hin和hex分別為飛艇內部自然對流換熱系數及外部混合對流換熱系數[18]，計算式如下

(19)

(20)

式中：λg為內部填充氣體的熱導率；D為蒙皮單元所在飛艇截面的直徑；參數Cl=0.515；Ra為局部瑞利數；n為與流動特征有關的變量.由于飛艇外表面與大氣既存在一定的相對運動，又有溫度上的差別，因此飛艇外表面與空氣之間的換熱為混合對流換熱，要綜合考慮自然對流和強迫對流兩種換熱，hfree_ex及hforced_ex分別為外部自然對流及外部強迫對流換熱系數，當自然對流與強迫對流同向或橫向時取正值，逆向時取負值.一般情況下n取值為3，但在涉及水平平板和圓柱的橫向流動時，n分別取3.5和4更為合適.本文在對飛艇熱特性仿真時采用正號，n取4.

由于蒙皮外表面與大氣存在一定的相對運動，需要考慮外部強迫對流換熱，外表面局部強迫對流換熱系數計算公式為[19-20]：

hforced_ex=

在計算飛艇外部自然對流換熱時，可將飛艇按截面分成許多個微元柱進行分析，根據流動特性的不同可將換熱分為層流和湍流.自然對流初期，邊界層厚度遠遠小于飛艇截面直徑，換熱近似為層流傳熱.隨著流動距離的增加，邊界層厚度越來越大，與此相對應，表面傳熱系數越來越小.當邊界層厚度達到臨界值時，層流傳熱將轉變為湍流傳熱，邊界層厚度保持恒定，不再與流動距離相關.層流湍流的分界點對于外部自然對流換熱的計算至關重要.本文流體介質為空氣，其普朗特數約為 0.733 6.

圖5所示為為圓柱體層流湍流轉化分界點與瑞利數關系，圖中：RaD為直徑D處的局部瑞利數；θc為臨界瑞利角度；θ為表面單元傾斜角；g為重力加速度.根據圖5可得在Pr=0.72時層流湍流臨界角度θc與瑞利數之間的關系，從而判斷飛艇表面單元的自然對流換熱是層流還是湍流.

圖5 圓柱體層流湍流轉化分界點與瑞利數關系[21]

對直徑為D的圓柱，壁面溫度和來流溫度均為常數，x為沿流線的距離.圖6所示為飛艇外部自然對流截面參數示意圖，圖中：r為計算點與圓柱軸線的距離；Δl為計算點與圓柱壁面的距離.

圖6 飛艇外部自然對流截面參數示意圖

外部自然層流換熱與湍流換熱的對流換熱系數以及努塞爾數的計算式如表1所示，表中：NuD為直徑D處的努塞爾數；Nux為直徑x處的努塞爾數；Rax為x處的局部瑞利數；kair為空氣的導熱系數；Cc=0.49[Pr/(0.861+Pr)]1/4;Ct=0.14Pr0.084；A(θ)是與θ有關的函數：

表1 層流及湍流對流換熱系數及努塞爾數計算式

(21)

1.3 內部填充氣體熱力學方程

飛艇內部填充氣體受太陽輻射、長波輻射以及與內表面蒙皮之間的自然對流換熱的作用，因此熱力學方程式可以寫成：

(22)

式中：cgas、mgas及Tgas分別為內部填充氣體的比熱容、質量及溫度；∑Qci,i為單位時間內填充氣體與所有蒙皮單元自然對流換熱量總量；內部填充氣體在單位時間內吸收的太陽輻射總量以及長波輻射總量分別為

(23)

(24)

(25)

2 計算方法與模型驗證

首先通過CFD軟件建立了所設計飛艇的3D幾何模型，將飛艇表面蒙皮劃分成若干個三角形單元，然后采用Tecplot對單元數據進行轉換，形成了MATLAB仿真程序需要的飛艇模型數據.基于上述的飛艇熱力學特性的數學模型，使用MATLAB編程，實現了整個熱模型的仿真.最后采用Tecplot以及Origin軟件處理了計算結果，并對結果進行了一系列的分析.

2.1 時間步長及網格對計算結果的影響分析

本文數值仿真的計算程序中，有蒙皮單元以及時間兩個離散量.為在保證計算結果的準確性的前提下提高計算的效率，分析了蒙皮網格數與時間步長對飛艇熱特性計算的影響.以兩個晝夜(48 h)內的飛艇熱特性變化為1個算例，初始計算參數如表2所示，填充氣體為He.

表2 算例基本輸入參數

分別計算了不同時間步長dt以及不同表面網格數Nm下的飛艇蒙皮及內部填充氣體的熱特性.圖7所示為不同表面網格數的直觀網格剖分圖.內部填充氦氣溫度以及蒙皮最高超熱溫度隨時刻的變化分別如圖8和9所示.圖中：Tmax為蒙皮最高溫度；Tmin為蒙皮最低溫度；t為實驗總時間.

圖7 不同表面網格數的直觀網格剖分圖

圖8 不同時間步長下飛艇熱特性48 h內的變化

由此可知，當時間步長在200 s以內時，時間步長對飛艇表面蒙皮以及內部填充氣體熱特性沒有影響，但當時間步長達到225 s時， 內部填充He的溫度與其余時間步長下的計算值有了細微偏差，位于He最高溫度右側，且蒙皮最高超熱溫度在白天也開始出現了振蕩， 可以判斷225 s的時間步長對于計算模型來說過大，會導致計算結果不準確.而不同網格數對填充氣體的熱特性基本無影響，對蒙皮溫度的影響也很小，最大溫度偏差為0.7 K.

圖9 不同網格數下飛艇熱特性48 h內的變化

因此，本文在之后的計算中采用200 s作為計算的時間步長，采用 4 666 作為計算的網格數，既能保證計算結果的準確性，也能提高計算效率.

2.2 模型驗證與結果分析

為了驗證本文飛艇熱特性模型、計算程序及方法的準確性，針對文獻[16]中的實驗進行了仿真，并將數值計算結果與文獻中的實驗結果進行對比分析.文獻[16]模型實驗通過太陽模擬器模擬太陽輻射，通過環繞在飛艇表面及內部的18個銅-康銅熱電偶采集蒙皮和內部氣體的溫度，并用風機以及引氣道控制飛艇的外部對流環境.具體實驗參數見表3.

表3 文獻[16]實驗基本幾何參數

從圖10中可知，數值計算出的飛艇內部空氣溫度的平均值在17，18號傳感器測量的空氣溫度值之間，與測量值吻合較好且變化趨勢基本相同.

圖10 內部填充氣體平均溫度的計算值與實驗值對比

來流速度為0，計算時間為300 s時，由實驗測得的上部填充氣體溫度(18號傳感器30.7 ℃)及下部填充氣體溫度(17號傳感器17.0 ℃)，可得填充氣體平均溫度實驗值為23.9 ℃，而填充氣體平均溫度計算值為25.7 ℃，比填充氣體平均溫度的實驗值高1.8 ℃，相差7.5%.計算時間為 600 s 時，由實驗測量得到的上部及下部填充氣體溫度(18號傳感器 38.7 ℃，17號傳感器19.6 ℃)，可得填充氣體平均溫度實驗值為29.2 ℃，而填充氣體平均溫度計算值為31.6 ℃，比填充氣體平均溫度實驗值高 2.4 ℃，相差8.2%.

當來流速度為7 m/s，計算時間為600 s時，實驗測量得到的上部及下部填充氣體溫度分別為21.3 ℃和16.0 ℃， 因此填充氣體平均溫度實驗值為18.7 ℃，而填充氣體平均溫度計算值為18.6 ℃，比填充氣體平均溫度實驗值低0.1 ℃，相差0.5%.計算時間為 1 200 s時，實驗測量得到的上部及下部填充氣體溫度分別為21.3 ℃及15.9 ℃，填充氣體平均溫度實驗值為18.6 ℃，填充氣體平均溫度計算值為19.3 ℃，比填充氣體平均溫度的實驗值高出 0.7 ℃，相差3.8%.

考慮到文獻[16]模擬實驗中，太陽模擬器的誤差約為±5%，風速儀的誤差約為±2%，熱電偶的測量誤差約為±1%，再綜合考慮實驗中存在的其他環境與人員操作誤差，可認為數值計算的結果與實驗結果較好地吻合.

3 飛艇熱力學特性分析

本文研究了基于相關參考文獻確定了一些可靠的初始輸入參數[13,16].蒙皮材料的熱輻射特性如表4所示.其余具體初始輸入參數如表2所示.

表4 蒙皮材料熱輻射特性

飛艇囊體的超冷超熱特性是蒙皮自身以及內部氣體熱輻射特性綜合作用的結果，通過圖11和12給出的3種典型材料下飛艇內部He和外蒙皮溫度變化可知，透明材料因其具有低的吸收率和高透射率，可使飛艇囊體內部He和蒙皮的超冷超熱量均低于半透明材料及不透明材料，是更加理想的飛艇蒙皮材料. 考慮材料的透射特性以及內部填充氣體對輻射的吸收與發射特性，能夠更加準確地獲得飛艇內部氣體及蒙皮的溫度范圍.

圖11 不同蒙皮材料的飛艇內部填充氣體溫度48 h內的變化

圖12 不同蒙皮材料的飛艇蒙皮最高溫度48 h內的變化

4 結語

平流層飛艇在駐空期間的熱特性會影響其浮力、蒙皮應力及整艇變形，浮力的變化和囊體變形引起的整艇阻力的顯著變化會對飛艇的飛行性能產生顯著的影響，故準確的飛艇熱模型及其熱特性分析成為其關鍵技術之一.

本文在已有研究的基礎上，根據常用飛艇蒙皮材料以及內部填充氣體的特性，將蒙皮的透射特性和內部填充氣體對輻射吸收與發射特性這兩方面的因素進行考慮，推導出飛艇表面蒙皮及內部填充氣體的熱力學方程.基于有限拆分法建立了考慮蒙皮透射率的飛艇熱力學仿真模型，并對其求解精度進行了驗證.

最終基于典型蒙皮材料的熱輻射特性，給出飛艇填充氣體及外蒙皮的熱特性及其變化規律，可為飛艇設計以及駐空期間的飛行性能評估提供參考和指導.