多層熱防護(hù)服一維瞬態(tài)導(dǎo)熱優(yōu)化模型

朱家明，左正東，杜佳軒，王浩歌

（1.安徽財(cái)經(jīng)大學(xué) 統(tǒng)計(jì)與應(yīng)用數(shù)學(xué)學(xué)院，蚌埠 233030；2.安徽財(cái)經(jīng)大學(xué) 金融學(xué)院，蚌埠 233030）

大量的科學(xué)統(tǒng)計(jì)研究表明，高溫極端環(huán)境會使暴露其中的工作人員發(fā)生中暑等癥狀，各地方政府推出了“高溫限工令”，即規(guī)定氣溫超過35℃的，要避免11：00-15：00室外露天和高處作業(yè)，此項(xiàng)舉措對保護(hù)工作人員的健康和權(quán)益方面有一定的積極作用［1］。體表溫度作為人體生理指標(biāo)之一［2］，在室內(nèi)熱環(huán)境學(xué)、生理學(xué)等領(lǐng)域有著重要的研究意義。通過探究多層熱防護(hù)服一維瞬態(tài)導(dǎo)熱規(guī)律，在生產(chǎn)設(shè)計(jì)高溫?zé)岱雷o(hù)服方面有一定的借鑒意義。

數(shù)據(jù)來源于2018年高教杯全國大學(xué)生數(shù)學(xué)建模比賽A題，為了便于解決問題，提出了如下假設(shè)：（1）無內(nèi)熱源且初始溫度均勻分布；（2）材料成本與厚度成正比；（3）熱量傳遞的方向視為一維垂直于皮膚表面輻射；（4）熱防護(hù)服裝的織物各項(xiàng)指標(biāo)性質(zhì)相同；（5）該模型僅考慮熱傳遞；（6）傳導(dǎo)和輻射在織物中熱傳遞是均勻的；（7）防護(hù)服最外層和人體皮膚傳遞過程中的輻射可以忽略；（8）空氣層的厚度值不超過6.4mm（忽略熱對流影響）；（9）織物與織物間，織物與外界空氣間，空氣與人體皮膚之間的溫度變化連續(xù)的且梯度跳躍。

目前關(guān)于熱防護(hù)服內(nèi)部傳熱模型根據(jù)熱防護(hù)服材料層數(shù)分為單層和多層模型。單層熱防服模型只有外殼，國內(nèi)外學(xué)者主要研究外界環(huán)境的輻射熱量、織物性質(zhì)、空氣層厚度對防護(hù)服熱性能的影響［3］。Gibson［4］首先提出高溫單層多孔介質(zhì)傳熱傳質(zhì)模型，Torvi［5］考慮不同輻射強(qiáng)度對熱防服內(nèi)部熱傳遞的影響。Ghazy［6］利用單層織物中熱傳導(dǎo)及比熱采用經(jīng)驗(yàn)公式，將常量轉(zhuǎn)為變量，更好地?cái)M合了高溫情況的熱傳遞過程。

許多學(xué)者在單層模型的基礎(chǔ)上改進(jìn)并探究了熱防護(hù)服熱濕傳遞多層模型。Mell［7］提出了多層面料模型層與層之間的傳熱模型。Mercer［8］和Ahmed Elgafy［9］考慮到由于相變材料對防護(hù)服的熱防護(hù)效果的影響，建立了內(nèi)含相變材料的多層動態(tài)熱傳遞模型。

1 基于傅里葉定律和COMSOL仿真模型對溫度分布的計(jì)算

1.1 研究思路

首先分析單層無限長圓筒溫度分布的一般情況。以此為基礎(chǔ)進(jìn)一步探討一維平板穩(wěn)態(tài)線性變溫導(dǎo)熱率的情況下的溫度分布，根據(jù)邊界條件引出四種組合，再建立非線性方程組，通過求解雅克比矩陣，并使用牛頓迭代法求解，給出溫度位置分布函數(shù)。考慮到不同層厚度差異較大，同層不同厚度的位置其溫度存在著差別，為了更加準(zhǔn)確地計(jì)算溫度分布狀況，將層數(shù)進(jìn)行進(jìn)一步細(xì)化，利用COMSOL軟件進(jìn)行仿真，得到各個厚度位置的每秒具體的溫度值，可以直觀地看出傳熱過程。織物厚度越厚，隔熱效果越好，皮膚表面溫度也越低。對第II層進(jìn)行掃描。改變第II層和第IV層的厚度時，從節(jié)約材料的角度考慮，應(yīng)該盡可能減少第II層的厚度，增大第IV層的厚度，以達(dá)到相同的效果。首先固定第II層的厚度，改變第IV層的厚度。最大的第IV層厚度的同時也須保證第II層厚度最小，此時為最優(yōu)值。

1.2 模型的建立與分析

1.2.1 溫度的時間空間分布

為了簡化問題，假設(shè)人為半徑220mm，高1700mm的圓柱，當(dāng)人穿著防護(hù)服時，在柱坐標(biāo)系下類似于多層圓筒壁（見圖1），一維非穩(wěn)態(tài)導(dǎo)熱實(shí)際上是二維問題，需要考慮時間（單向）和空間坐標(biāo)。

圖1 假人高溫作業(yè)專用服裝穿著仿真圖

在圖1中給出了由服裝——皮膚及空氣層組成的系統(tǒng)，其中，防護(hù)服通常由三層織物材料構(gòu)成，記為I、II、III層，其中I層與外界環(huán)境接觸，III層與皮膚之間還存在空隙，將此空隙記為IV層。

從外到內(nèi)，織物距離皮膚的距離以次為d1、d2、d3，材料的熱傳導(dǎo)率以次為k1、k2、k3，表面溫度值以次為T1、T2、T3、T4。如圖2所示。

圖2 多層結(jié)構(gòu)溫度分布邊界結(jié)構(gòu)簡圖

其中T1=75℃，T5=37℃ ，d2-d3=6mm ，d4-d5=5mm。

導(dǎo)熱系數(shù)的計(jì)算式是由傅里葉導(dǎo)熱定律的數(shù)學(xué)表達(dá)式導(dǎo)而出，即：

導(dǎo)熱系數(shù)是表示物質(zhì)導(dǎo)熱能力的大小的一項(xiàng)基本宏觀物理性質(zhì)。本題中導(dǎo)熱系數(shù)被當(dāng)作常數(shù)對待。大部分情況下，導(dǎo)熱系數(shù)與溫度呈非線性關(guān)系的。通常在實(shí)際計(jì)算中，若溫度變化范圍不大，其數(shù)值可近似采用線性關(guān)系，即

式中，T為溫度，b代表該溫度范圍內(nèi)導(dǎo)熱系數(shù)隨溫度的變化，其值為常量，b＞0時熱導(dǎo)率與溫度成正比，b＜0時熱導(dǎo)率與溫度成反比，b=0時熱導(dǎo)率是常數(shù)［10］。

下式中的系數(shù)Ai、Bi等值為由邊界條件確定的系數(shù)。

由題可知，兩層織物的熱流密度相等，即

由式（3）、（4）、（5）、（6）整理得：

方程（9）是一個關(guān)于中間界面溫度T2、T3的非線性方程組，可用牛頓法解出：

雅克比矩陣為：

代入（11）迭代，可求得T2、T3、T4，然后計(jì)算溫度分布等。

第IV層溫度分布：

考慮到不同層厚度差異較大，同層不同厚度的位置其溫度存在著差別，在（13）-（16）式的基礎(chǔ)之上，為了更加準(zhǔn)確地計(jì)算溫度分布狀況，將層數(shù)進(jìn)行了進(jìn)一步細(xì)化，選取了52個節(jié)點(diǎn)，將0～15.2mm的厚度劃分為了51個細(xì)微厚度層，在厚度層內(nèi)，溫度基本無差異。而多層圓筒壁導(dǎo)熱的一維問題，在垂直導(dǎo)熱的假設(shè)前提下，可以進(jìn)一步簡化為圖3，圖3從左到右衡量分別為I、II、III、IV層，且在51個區(qū)間內(nèi)每秒溫度存在著瞬時變化，每個區(qū)間的垂直圓筒壁的溫度在每一秒下是相等的。

圖3 一維簡化仿真圖

利用COMSOL軟件進(jìn)行仿真，可以得到各個厚度位置的每秒具體的溫度值，將其繪制為溫度隨厚度與時間變化的三維立體圖（圖4），可以更直觀的看出傳熱過程，導(dǎo)熱進(jìn)行50秒后，初始條件37℃對系統(tǒng)中無量綱分布的影響趨近于0，溫度分布此時主要取決于第三邊界條件的影響。其在1100mm時基本達(dá)到穩(wěn)態(tài)，溫度基本保持穩(wěn)定。溫度時間分布誤差分析如圖5所示。

圖4 溫度厚度時間三維坐標(biāo)系圖

圖5 溫度時間分布誤差分析圖

將所得到的假人皮膚外側(cè)的實(shí)時溫度值（即厚度為15.2mm處）與數(shù)據(jù)中給的測定值進(jìn)行比較（見圖6），可以看出本模型的熱傳導(dǎo)預(yù)測效果較好，基本吻合實(shí)測值，只在1000秒到1500秒之間略有偏差，因此本模型準(zhǔn)確度較高。

圖6 預(yù)測溫度與實(shí)際溫度差值相關(guān)圖

再利用MATLAB繪制三維等溫線圖（見圖7），在三維空間內(nèi)的同一等溫線上溫度相等，到達(dá)穩(wěn)態(tài)的實(shí)際根據(jù)厚度的不同而變化。

圖7 空間直角坐標(biāo)材料三維等溫線分布圖

利用COMSOL Multiphysics軟件的一維瞬態(tài)無內(nèi)熱源導(dǎo)入公式（17）（18）：

以上是熱傳導(dǎo)方程，Ac是一維線的橫截面積，ρ是密度，Cp是比熱容，T是溫度，t是時間，u是速度，固體不流動，不用考慮。k是導(dǎo)熱系數(shù)。等式左邊第一項(xiàng)代表材料溫度改變，第三項(xiàng)代表熱傳遞貢獻(xiàn)，等式右邊是熱量產(chǎn)生，在本模型種并沒有熱量生成，不需要考慮。

1.2.2 基于溫度分布曲線對防熱服厚度的最優(yōu)設(shè)計(jì)

織物厚度越厚，隔熱效果越好，假人皮膚外側(cè)溫度也越低。

利用之前的函數(shù)：

利用局部差分法，對該方程進(jìn)行二次積分，求出溫度分布的函數(shù)：

若要求在確保工作60分鐘，假人皮膚外側(cè)溫度不超過47℃，且超過44℃的時間不超過5分鐘的情況下，確定第II層的最優(yōu)厚度。最優(yōu)厚度即為滿足上述約束條件下的該層最薄厚度，這樣節(jié)省材料費(fèi)用，達(dá)到最優(yōu)，即為求min（d2-d3）。約束條件為：

其中，t(T,d)為T(t,d)的相關(guān)函數(shù)。

1.3 模型結(jié)果及分析

在其他變量不變的情況下，在第II層厚度的取值區(qū)間內(nèi)，不斷變換第II層厚度，得到不同厚度下的假人皮膚外側(cè)實(shí)時溫度分布狀況繪制如圖8所示。

圖8 假人皮膚外側(cè)溫度分布曲線圖

結(jié)果發(fā)現(xiàn)當(dāng)?shù)贗I層厚度達(dá)到19mm時，在55分鐘時，溫度達(dá)到44度，并且在60分鐘時，溫度也只是稍微高于44度，達(dá)到臨界值，此時的溫度分布曲線如圖9所示。

圖9 假人皮膚外側(cè)溫度分布曲線圖

為考慮節(jié)約材料，最優(yōu)厚度應(yīng)該定位在稍厚于19mm的19.5mm。

表1是當(dāng)?shù)贗I層厚度為19.5mm時，部分時間的溫度（單位為開爾文）變化。

當(dāng)?shù)贗I層和第IV層的厚度時可以改變的，從節(jié)約材料的角度考慮，應(yīng)該盡可能減少第II層的厚度，增大第IV層（該層為空氣層）的厚度，以達(dá)到相同的效果。首先固定第II層的厚度，改變第IV層的厚度。假設(shè)固定第II層為10mm，逐漸增加第IV層的厚度，通過COMSOL Multiphysics軟件掃描發(fā)現(xiàn)當(dāng)?shù)贗V層厚度達(dá)到最大值6.4mm時依然不滿足確保工作30分鐘時，假人皮膚外側(cè)溫度不超過47℃，且超過44℃的時間不超過5分鐘的條件。具體溫度變化如圖10所示。

圖10 固定第II層厚度改變第IV層厚度的溫度變化線圖

表1 假人皮膚外側(cè)溫度隨時間變化一覽表

因此逐漸增加第II層厚度，使得第IV層厚度在最大值6.4mm時，滿足要求條件。最大的第IV層厚度也保證第II層厚度最小，從而達(dá)到節(jié)約材料的目的。具體溫度變化如圖11所示。

結(jié)果發(fā)現(xiàn)當(dāng)?shù)贗I層厚度為22mm，第IV層厚度為6.4mm時，滿足要求條件，為最優(yōu)值，溫度曲線如圖12所示。

圖11 固定第IV層厚度改變第II層厚度的溫度變化線圖

圖12 第IV層、第II層最優(yōu)厚度的溫度變化線圖

2 結(jié)論

基于傳熱學(xué)的理論，對不同假設(shè)下、不同初始條件、以及不同形式的存在耦合關(guān)系的多層圓筒體模型進(jìn)行探究，以界面溫度為解耦參數(shù)，將多層傳熱簡化成單層傳熱問題，從而推導(dǎo)出多層圓筒體在不同邊界條件下的溫度分布表達(dá)式。通過借助COMSOL 5.3a軟件，以有限元法為基礎(chǔ)，通過求解偏微分方程來實(shí)現(xiàn)真實(shí)導(dǎo)熱現(xiàn)象的仿真，將復(fù)雜的多層圓筒壁瞬態(tài)導(dǎo)熱模型進(jìn)行簡化，大大降低了計(jì)算難度。文中建立的模型與緊密結(jié)合實(shí)際情況，突出了傳熱系數(shù)的作用，充分考慮物理學(xué)原理，從而使模型更加通用、易懂。