非穩態自然對流換熱系數計算方法及其在防護服隔熱預報中的運用

丁 寧， 林 潔

(1. 中國民航大學 航空工程學院， 天津 300300； 2. 中國民航大學 中歐航空工程師學院， 天津 300300)

在船舶制造行業，一線人員經常需要在雙層底、雙舷側以及其他密閉艙室內實施焊接、檢測等作業。這類作業環境內的溫度往往較高，傳熱方式以自然對流為主，此外，在建筑行業等其他領域，一線人員也同樣需要在這類環境中工作。為保護人體免遭損害，人員在進入這類環境時往往需要穿著防護服，因此防護服的隔熱性能至關重要。采用數值方法研究防護服的隔熱性能，因其成本和靈活性方面的優勢，越來越受到重視。

因隔熱性能與熱源密切相關，學者們根據不同的熱源特點，經適當假設后建立了多種數學模型，并相應采用了多種數值方法進行求解。這些熱源主要是石油化工行業的強熱流閃火環境和消防行業的中等熱流輻射環境[1]，傳熱方式主要是對流和輻射。其中，自然對流產生的換熱占比較小忽略不計，Torvi[2]對這類熱源條件進行了研究，將輻射和對流換熱合并為單一系數并視為常量給出了簡潔的表述形式，后續較多研究也均采納了這一表述[3-5]。為考慮自然對流時的傳熱情況，Mell等[6]采用Howard模型給出的經驗公式[7]，在常熱流條件下，利用Nusselt數計算了自然對流換熱系數進而獲得了織物溫度分布的數值解，并與解析解進行了對比。然而，Howard模型通常適用于常壁溫或常熱流的對流邊界條件，對于變壁溫、變熱流的對流邊界條件是否適用，并未有明確結論，也尚未見有其他相關研究成果，而在高溫環境中，情況完全不同，其傳熱方式以自然對流為主，在達到穩態前，自然對流換熱系數一直在變化，對流邊界上的壁溫和熱流也一直在變化，采用Howard模型難以計算。此外，Su[8]采用了Howard模型計算自然對流邊界條件下的冷卻問題，從與實驗結果的對比上看差異較大。從上述分析可以看出，盡管熱防護服的邊界條件是影響模型精度的重要因素[2]，但對非穩態條件下自然對流換熱系數的研究并未引起足夠的重視。如陳揚等[9]在自然對流傳熱條件下對織物的非穩態傳熱進行了模擬，并與實驗值進行了對比，其自然對流換熱系數則直接指定為10 W/(m2·K)，并未給出相關確定依據。Piotr等[10]研究了輻射傳熱占優的熱防護服非穩態傳熱問題，其自然對流換熱系數也未加說明而直接指定。目前，在防護服隔熱性能非穩態數值預報方面，強迫對流傳熱和輻射傳熱均有較為成熟的公式可以利用，而針對自然對流傳熱不可忽略的熱源環境，尚未見有計算非穩態自然對流傳熱系數的相關成果。

為此，本文針對置于高溫艙室中的熱防護服材料，建立了服裝隔熱性能預報的數學模型，提出了一種計算非穩態自然對流換熱系數的方法，利用有限差分法計算了模擬皮膚外側的溫度時程曲線，并與實驗測量值進行了對比，結果顯示文中自然對流換熱系數計算方法有效，可明顯改善防護服隔熱性能的數值預報精度，在壓縮防護服研發周期和實驗成本方面意義明顯。

1 傳熱模型

1.1 模型假設

圖1示出了傳熱系統的組成，為便于后續推導，構成熱防護服的3層材料由外到內依次編號為1、2、3，空氣層編號為4，模擬皮膚層編號為5。

對于圖1所示傳熱系統，假設：1)熱量沿垂直于模擬皮膚方向傳遞，為一維傳熱問題；2)能量從外界環境向熱防護服傳播時，考慮對流和輻射；3)各材料的熱屬性穩定，與溫度無關；4)材料4的厚度不大于6.4 mm，不計熱對流的影響[2]；5)恒溫內核溫度為37 ℃。

1.2 傳熱系統微分方程

基于以上假設，對上述瞬態傳熱系統的數學描述如下。各材料內部的熱量傳遞滿足：

(1)

式中：Ti(xi,t)為不同材料內部溫度分布，分別與相應標號的材料對應，K；xi為水平坐標，采用局部坐標系，即各材料的坐標原點在該材料的左端壁上，x軸正向向右，m；t為時間，s；λi為不同材料的熱傳導率，W/(m·K)；ρi為不同材料的密度，kg/m3；ci為不同材料的比熱容，J/(kg·K)。

材料1外側(左側)熱流邊界條件：

(2)

材料1與材料2及材料2與材料3的熱流邊界條件：

(3)

式中：δi為不同材料厚度，m。

材料3與材料4(空氣層)間的熱流邊界條件：

(4)

式中：ε=1/(1/εin+1/εs-1)，εin為材料3發射率，εs為模擬皮膚發射率。

材料4(空氣層)與材料5間的熱流邊界條件：

(5)

材料間的接觸邊界溫度邊界條件：

Ti|xi=δi=Ti+1|xi+1=0(i=1，…，4)

(6)

各材料的初始條件為：

Ti|t=0=C0(i=1，…，5)

(7)

T5|x5=δ5=C0

(8)

式中：C0為恒溫內核溫度，K。

1.3 非穩態自然對流換熱系數

實際上，外部自然對流換熱系數與初始和穩態自然對流換熱系數關系十分密切，在自然對流狀態下，邊界層內的空氣從初始時的靜止到穩態時的運動，其速度變化較為平緩，相應自然對流換熱系數也應是平緩、連續的變化，因邊界層內的空氣流動是因空氣溫度降低而至空氣收縮和密度增加所致，而邊界層內的空氣溫度取決于環境溫度和防護服的表面溫度，在環境溫度恒定的情況下，防護服的表面溫度是引起邊界層內空氣流動的直接因素，因此，有理由認為非穩態條件下的自然對流換熱系數是關于防護服表面溫度的連續函數，該連續函數有2個值已知，分別對應于初始時和穩態值。

為此，首先確定穩態時的傳熱系統的自然對流換熱系數，其數學描述如下：

各材料內部的熱量傳遞滿足：

(9)

材料1外側(右側)熱流邊界條件：

hc(T1|x1=0)·(Tair-T1|x1=0)

(10)

式中：hc(T1|x=0)為穩態時外部環境自然對流換熱系數，是T1|x=0的函數[7]，W/(m2·K)，其他參數如前所述。

材料1與材料2及材料2與材料3的熱流邊界條件見式(3)；材料3與材料4(空氣層)間的熱流邊界條件見式(4)；材料4(空氣層)與材料5間的熱流邊界條件見式(5)；材料間的接觸邊界溫度邊界條件見式(6)。通過解方程可求出穩態時熱防護服表面溫度T1|x=0(令其為C1)以及穩態時的自然對流換熱系數hc(C1)。

在初始條件下，上述傳熱系統被置于高溫艙室中，邊界層內的空氣尚未開始運動，其熱傳遞形式主要是熱傳導，將空氣層離散，離散步長為Δxair，并設空氣的熱導率為λair，根據傅里葉定律有

(11)

整理后可得初始條件下的自然對流換熱系數為

(12)

(13)

將(13)代入式(3)，方程(1)～(8)可解。

2 傳熱系統數值算法

2.1 穩態傳熱系統數值算法

構造函數：f(T1,1)=q0(T1,1)-q′0(T1,1)，

T1,1∈(Tair,C0)

式中：Ti,j表示第i種材料第j個節點的溫度，將公式(9)～(10)、(3)～(6)寫成差分格式整理后可得上式中其他符號的含義，如下：

T4,N4=δ5q0/λ5+C0

穩態時應有f(T1,1)=0，可以證明，q′0為T1,1的單增函數，q0為T1,1的單減函數，采用迭代法可求解方程f(T1,1)=0，得到T1,1(即C1)，進而得到hc(C1)，由式(17)可得h′c(T1|x=0)。

2.2 瞬態傳熱系統數值算法

將式(1)寫成差分格式：

(i=1，…，4)

(14)

式中：Ti,j,k表示材料i的節點j在tk時刻的溫度。根據式(7)給出的初始條件，4種材料在t2時刻的內部節點溫度均可由式(13)確定，t2時刻材料邊界處的溫度則需單獨求解。

因T1,2,2為內部節點，則其值為已知，將式(2)左端寫成向前差分格式，整理后可得：

(15)

式(15)是關于T1,1,2的4次方程，可采用迭代法求解。

將式(3)寫成差分格式，并結合式(6)，整理后可得

(16)

(17)

將式(4)、(5)寫成差分格式，整理后可得：

(18)

T4,N4,2=-P1T3,N3,2+P2

(19)

式中：Ck3=λ3/Δx3、Ck41=λ4/Δx4、Ck4=λ4/Δx4、Ck5=λ5/δ5，

P1=(Ck3+Ck41)/(Ck4+Ck5)

式(18)是關于T3,N3,2四次方程，可采用迭代法求解，由式(19)可確定T4,N3,2。

3 算例及分析

各材料的物理屬性[13]見表1，環境溫度Tair=75 ℃、恒溫內核溫度C0=37 ℃，其他材料初始溫度為37 ℃。

表1 材料的物理屬性Tab.1 Thermophysical/geometrical properties of different material

利用上述方法對圖1所示傳熱系統進行分析，服裝內表面發射率εin取0.8，模擬皮膚發射率εs取0.97，服裝外表面發射率εout取0.1，λ5/δ5=5.433 W/(m·K)，標準大氣壓下空氣的物理屬性由美國國家標準局氣體熱性能表No. 564根據空氣溫度插值確定。

圖2示出防護服外表面在不同邊界條件和不同時間步長下模擬皮膚外側溫度的計算結果。

虛線是防護服外表為恒溫75 ℃時的情況，計算結果均高于實驗值[13]，趨于穩態時最大差值為3.92 ℃，但趨于穩態的時間與實驗值吻合較好，該邊界條件高估了防護服外表面的熱傳遞，總體上與實驗結果偏差較大；點劃線是防護服外表面自然對流系數按Howard模型計算時的結果，計算結果均低于實驗值，非穩態階段最大差值為6.6 ℃，趨于穩態時溫度值與實驗結果吻合較好，但在時間上有較大的滯后，該邊界條件低估了防護服外表面的熱傳遞，總體上與實驗結果偏差更大，盡管如此，利用該邊界條件計算穩態溫度情況時卻可以獲得令人滿意的結果；實線是上述2種情況計算結果的加權平均，經該方式處理后的數據在一定程度上得到了改善，但與實驗結果仍有不小的差距，尤其是在穩態階段。從上述分析可以看出，防護服外表面的傳熱邊界條件對計算結果有著非常重要的影響，是影響防護服隔熱性能模擬精度的關鍵因素，從對2種典型邊界條件的計算可以看出，模擬結果與實驗結果差別較大，雖在一定程度上反映出了溫度隨時間的變化情況，但仍然不夠理想。

圖3示出時間步長為0.001 s時計算結果相對實驗測量結果的誤差。可以看出，隨著時間步長的加密(空間步長相應加密)計算結果迅速收斂，當時間步長為0.001 s時，計算結果與實驗結果高度一致，從圖3可以看出整體誤差小于0.1 ℃，計算結果令人滿意。即便是時間步長為0.1 s的情況，其計算結果也優于圖2中利用Howard模型所得的結果。

圖4示出時間步長為0.001 s時熱流密度隨時間的變化。

從圖4可以看出，熱量主要來源于外界環境的對流傳熱，輻射傳熱僅占較小的比例，起始時的總熱流密度最大，接近1.7 kW/m2，幾乎全部來自對流傳熱，隨著時間的變化，自然對流的熱流密度迅速下降，輻射部分的熱流密度變化不大，相應總熱流密度迅速下降，趨于穩態時，自然對流的熱流密度達到最小值，約為46 W/m2，輻射熱流密度約為14 W/m2，總熱流密度達到最小值，約為60 W/m2。熱流密度的變化符合實際情況，也證實了高溫環境條件下，熱量主要來源于自然對流傳熱，在自然對流傳熱不可忽略的熱源環境下，自然對流傳熱系數計算的合理性將直接影響模擬結果的準確性。

4 結束語

自然對流傳熱不可忽略甚至占優的高溫環境是熱防護服一種較為典型的應用場景，文中探討了此時瞬態自然對流換熱系數的一種計算方法，依此方法預報了防護服隔熱性能并與實驗測量值進行了對比，驗證了該自然對流換熱系數計算方法，解決了自然對流傳熱不可忽略甚至占優時，防護服非穩態隔熱性能預報中的一個問題。對于作業于非穩態傳熱階段的防護服，在進行防護服隔熱性能數值模擬時，若高估自然對流的影響，則溫度分布會提前進入穩態，且穩態溫度高于實際溫度，這會導致防護服設計的過于保守；若低估自然對流的影響，則溫度分布會延后進入穩態，此時的設計將不足以隔熱，嚴重時會導致人體受到損傷。

此外，文中假設材料的物理屬性不隨溫度變化，但當環境溫度很高時，外側材料的物理屬性首先受到溫度的影響，下一步應將材料物理屬性隨溫度變化時的情況納入研究，以適應更惡劣的環境溫度。

致謝感謝天津海太科技有限公司馬曉梅總工在本文修改過程中提出的寶貴意見。