董文晉，朱家明

（1.暨南大學(xué) 經(jīng)濟(jì)學(xué)院，廣東 廣州 510632；2.安徽財(cái)經(jīng)大學(xué) 數(shù)量經(jīng)濟(jì)研究所，安徽 蚌埠 233030）

對(duì)于從事高溫作業(yè)的工作人員來(lái)說(shuō)，熱負(fù)荷是其生命安全的最大挑戰(zhàn)。熱防護(hù)服的厚度既要滿足能夠隔熱（安全的前提下）的功能，也要是其范圍內(nèi)的最薄厚度，方便作業(yè)人員工作。因此，對(duì)防護(hù)服的最優(yōu)厚度的研究具有著重要的意義。不僅關(guān)乎高溫作業(yè)人員的生命安全，也對(duì)防護(hù)服產(chǎn)業(yè)的發(fā)展具有借鑒意義。國(guó)內(nèi)有不少對(duì)這類(lèi)問(wèn)題的探討。付明等利用“暖體假人”實(shí)驗(yàn)進(jìn)行研究，得出高溫下防護(hù)服的熱阻并聯(lián)比串聯(lián)的要小的結(jié)論[1]。李紫含等建立了熱防護(hù)服降溫效果的指標(biāo)評(píng)價(jià)體系，對(duì)防護(hù)服降溫效果進(jìn)行了綜合評(píng)價(jià)[2]。汪明珮等運(yùn)用方程離散化以及元胞自動(dòng)機(jī)算法計(jì)算出各層織物材料的溫度分布[3]。丁寧等建立了模擬皮膚傳熱系統(tǒng)的非穩(wěn)態(tài)導(dǎo)熱方程，并結(jié)合實(shí)驗(yàn)結(jié)果，確定了模型中的參數(shù)，最后得出防護(hù)服的導(dǎo)熱特性[4]。本研究則是基于2018年全國(guó)大學(xué)生數(shù)學(xué)建模賽題，利用賽委會(huì)提供的相關(guān)數(shù)據(jù)進(jìn)行研究。本研究的創(chuàng)新之處在于對(duì)溫度分布函數(shù)進(jìn)行靈敏度分析，證明了函數(shù)計(jì)算結(jié)果具有穩(wěn)定性；同時(shí)還對(duì)方程求解的假人皮膚溫度擬合的結(jié)果與實(shí)際數(shù)據(jù)進(jìn)行差異性檢驗(yàn)，兩者關(guān)聯(lián)程度較高，證實(shí)了傅里葉熱傳導(dǎo)模型求解結(jié)果的可靠性。

1 數(shù)據(jù)來(lái)源與假設(shè)

本研究數(shù)據(jù)來(lái)源于2018年全國(guó)大學(xué)生數(shù)學(xué)建模比賽A題，為了便于解決問(wèn)題，提出以下幾條假設(shè)：①假定在實(shí)驗(yàn)室的各個(gè)區(qū)域內(nèi)，溫度不變；②假定暖體假人是一個(gè)幾何圓柱體，熱傳遞只在垂直于皮膚的表面進(jìn)行，沒(méi)有進(jìn)行外延散熱；③熱量在織物材料中的傳遞是穩(wěn)定均勻的，忽略各層織物材料之間不同點(diǎn)的溫度變化，只考慮不同層之間的溫度分布，且認(rèn)為同一層織物材料之間的溫度變化是連續(xù)的；④假設(shè)人的體溫恒定為37 ℃，將零時(shí)刻IV層的溫度視為人的體表溫度37 ℃；⑤假定熱流的速度足夠小，熱傳遞過(guò)程中只考慮熱傳遞過(guò)程，忽略熱輻射以及熱對(duì)流現(xiàn)象；⑥假設(shè)溫度傳導(dǎo)過(guò)程是一維的，忽略面料和空氣層厚度的變化；⑦假定防護(hù)服的每一層面料質(zhì)地均勻，密度處處相等。

2 基于傅里葉熱傳導(dǎo)定律對(duì)多層織物材料溫度分布的研究

2.1 研究思路

高溫作業(yè)的專業(yè)服裝由不同屬性的三層織物材料構(gòu)成，本研究將空氣層作為Ⅳ層。不同層次構(gòu)成的功能分工和熱傳導(dǎo)現(xiàn)象在不同介質(zhì)中的傳遞機(jī)制：最外層（Ⅰ）作為外殼，Ⅱ?qū)拥淖饔弥饕欠浪?，Ⅲ層的主要作用是隔熱，空氣層（Ⅳ）是在人體表皮與隔熱層之間起到緩沖作用。由于熱傳導(dǎo)是介質(zhì)內(nèi)無(wú)宏觀運(yùn)動(dòng)時(shí)的傳熱現(xiàn)象，其在固體、液體和氣體中均可發(fā)生，但在流體中靜止時(shí)也會(huì)由溫度梯度而形成的密度差，從而產(chǎn)生自然對(duì)流，所以本研究不需考慮對(duì)流，只考慮熱傳導(dǎo)。根據(jù)熱能守恒定理和傅里葉熱傳導(dǎo)定律[5]推斷出一維層次熱傳導(dǎo)方程，在給定的限制條件下對(duì)方程進(jìn)行求解，得到溫度隨時(shí)間變化的多層溫度分布模型，最后將不同分層的物質(zhì)屬性系數(shù)代入方程，得到每一層的溫度隨時(shí)間變化情況。

基于傅里葉熱傳導(dǎo)定律，建立溫度、位置和時(shí)間的方程，并對(duì)微分方程進(jìn)行求解，解出三者之間的偏導(dǎo)等式，再建立關(guān)于等式的初值條件和邊界條件，求出差分后式子中的參數(shù)，得出三者的具體函數(shù)，代入每一層的厚度和導(dǎo)熱率，運(yùn)用MATLAB軟件編程，求出每一層的溫度分布圖，最后在模型中代入皮膚溫度進(jìn)行擬合度的檢驗(yàn)。研究思路流程圖如圖1所示。

圖1 思路流程圖Figure 1 Ideas flow chart

2.2 研究方法

本研究先把防護(hù)服理想化成n個(gè)圓柱體，類(lèi)似于具有一定橫截面積M的桿，其方向?yàn)樗剑傞L(zhǎng)度設(shè)為L(zhǎng)，位置的變量為l，設(shè)單位面積的熱能量為未知變量，叫作熱能密度，記為E(l,t)。

觀察桿是介于l和l+ Δl之間的薄片，在Δl非常小時(shí)，E(l,t)在薄片內(nèi)可近似為常數(shù)，這樣得出薄片的體積為M?Δl，此時(shí)熱能記為E(l,t)?M?Δl。在l和l+ Δl之間的熱能隨時(shí)間的變化都是由薄片兩端的熱能引起的。由于假設(shè)側(cè)面是絕熱的，所以在側(cè)面上沒(méi)有熱能變化，基本的熱流過(guò)程[6]表達(dá)為：

熱能的瞬時(shí)變化率=單位時(shí)間流過(guò)邊界的熱能，

這種現(xiàn)象稱之為熱能守恒。對(duì)小薄片而言，記熱能Q的變化率為η，表達(dá)式為

熱能的流向可為左，也可為右，在單位時(shí)間內(nèi)熱能流向單位表面積右邊的熱量稱之為熱通量，記為φ(l,t)。根據(jù)熱能守恒，熱能變化率[7]是由流過(guò)邊界的熱能和內(nèi)部熱源造成的，表達(dá)式為

由于對(duì)小橫截面積，通常會(huì)被近似看成常數(shù)，方程⑴不是精確的，于是本研究采用極限思想，得到

其中，常數(shù)橫截面積被消去了，因此這個(gè)結(jié)果斷然消除了誤差，所以用等號(hào)代替方程（1）中的約等號(hào)。在極限過(guò)程中t是固定的，因此，由偏導(dǎo)數(shù)的定義得

根據(jù)實(shí)驗(yàn)，通常對(duì)于限制溫度的區(qū)間，比熱容與溫度無(wú)關(guān)，因此建立方程c=c(l)。由于單層面料由一種物質(zhì)均勻構(gòu)成，因此就把比熱容看作常數(shù)，單位質(zhì)量的熱能就是T(l,t)?c(l)，這樣就要引入密度ρ(l)，允許它隨l變化，在任意切片內(nèi)的熱能是Q（Q=T(l,t)?c(l)?ρ?M?Δl），得到式子

這樣就解釋了熱能跟溫度之間的關(guān)系，同時(shí)單位體積的熱能等于單位質(zhì)量單位的熱能乘以密度，消去熱能密度后，熱能守恒式（4）變?yōu)?/p>

2.3 結(jié)果分析

求解一維熱傳導(dǎo)模型，先要確定初值條件和邊界條件。假人皮膚的溫度趨勢(shì)如圖2所示。

通過(guò)STATA軟件擬合得到函數(shù)表達(dá)式

根據(jù)不同層的密度、比熱容和熱傳導(dǎo)率求得不同的導(dǎo)熱能力，具體數(shù)值見(jiàn)表1。

聯(lián)合式（6）、（7），可得到偏微分方程

并運(yùn)用MATLAB軟件解出該方程，得到每層不同時(shí)間段溫度分布的效果圖（見(jiàn)圖3）以及每層不同時(shí)間段溫度分布的部分?jǐn)?shù)據(jù)（見(jiàn)表2）。

圖2 假人皮膚溫度變化趨勢(shì)圖Figure 2 Changing trend of skin temperature of dummy

表1 各層導(dǎo)熱能力分布表Table 1 Thermal conductivity distribution table of each layer

圖3 各層不同時(shí)間段溫度分布圖Figure 3 Temperature distribution map of each layer in different time periods

表2 各層溫度分布表Table 2 Temperature distribution table of each layer

3 基于有限差分法及粒子群算法對(duì)單層厚度的優(yōu)化設(shè)計(jì)

3.1 研究思路

基于多層織物材料溫度分布模型，結(jié)合Henriques積分作為皮膚燒傷程度的評(píng)判依據(jù)，為熱防護(hù)服設(shè)計(jì)提出反問(wèn)題，以安全為前提，以減少傷害程度為目的，最優(yōu)決定II層厚度。為此，我們運(yùn)用有限差分法建立了單層厚度優(yōu)化模型，并且同時(shí)使用粒子群算法求解。

3.2 研究方法

（1）熱傷害量化模型

用Henriques建立的熱傷害量化模型為

式中，Ω為熱損傷指數(shù)，t為時(shí)間（s），P為環(huán)境溫度和體溫有關(guān)的常數(shù)或者叫頻率因子（s-1），ΔE為激活能（J/mol），R為理想氣體常數(shù)（8.314 J/（mol·k）），T為絕對(duì)溫度。

為了改進(jìn)和提高熱防服系統(tǒng)的整體性能，確定熱防服的第II層最優(yōu)厚度，本研究運(yùn)用Henriques積分作為人體皮膚安全的評(píng)測(cè)標(biāo)準(zhǔn)，因此，當(dāng)人的表皮絕對(duì)溫度達(dá)到44 ℃時(shí)，就認(rèn)定其皮膚受到了熱損傷，通過(guò)將皮膚表層與接觸面的溫度帶入公式

因此，本研究在建立模型二過(guò)程中所提出的反問(wèn)題是：在外界環(huán)境溫度為65 ℃、IV層厚度為5.5 mm的前提下，求出滿足安全性條件Φ'（工作60 min 時(shí)，人皮膚外側(cè)溫度不超過(guò)47 ℃，且超過(guò)44 ℃的時(shí)間小于5 min）下的II層最優(yōu)決定厚度，因此根據(jù)約束條件和已知條件本研究可以建立基于溫度分布的第II層厚度的最優(yōu)解模型[11]，即：

3.3 結(jié)果分析

將上述式子改寫(xiě)成離散形式可得

將上述公式進(jìn)行MATLAB運(yùn)算后得出結(jié)果為8.8 mm，見(jiàn)圖4。

圖4 運(yùn)行結(jié)果圖Figure 4 Running result diagram

4 基于遺傳算法對(duì)高溫防護(hù)服厚度的優(yōu)化設(shè)計(jì)

4.1 研究思路

運(yùn)用遺傳算法來(lái)求解約束優(yōu)化問(wèn)題時(shí)，可通過(guò)懲罰函數(shù)方法轉(zhuǎn)換到目標(biāo)函數(shù)中，這樣可將初始問(wèn)題轉(zhuǎn)化成無(wú)約束最優(yōu)化問(wèn)題，再利用遺傳算法進(jìn)行求解。對(duì)于具體的優(yōu)化問(wèn)題，可以用圖5的流程進(jìn)行求解[12]。

圖5 思路流程圖Figure 5 Flowchart of thinking

4.2 研究方法

首先對(duì)公式（13）進(jìn)行正則化處理，再使用遺傳算法進(jìn)行求解，會(huì)做出如下目標(biāo)函數(shù)[13]：

其中參數(shù)α稱為正則化參數(shù)，該目標(biāo)函數(shù)的正則化解正是滿足安全性指標(biāo)和厚度限制條件的最優(yōu)解[14]。由此，可以厘清最優(yōu)解的求解過(guò)程，為接下來(lái)模型的求解做鋪墊。

4.3 結(jié)果分析

選取最小二乘解L作為GA搜索空間，跟一般的正則化問(wèn)題很類(lèi)似。在求解的過(guò)程中，對(duì)α處理直接影響結(jié)果的誤差。本研究將運(yùn)行50次的最小值作為最優(yōu)解，由此可以得到L=（11，4.8）為最優(yōu)解，即Ⅱ?qū)拥淖顑?yōu)厚度為11 mm，Ⅳ層的最優(yōu)厚度為4.8 mm。

5 結(jié)論

根據(jù)熱能守恒定理和傅里葉熱傳導(dǎo)定律，構(gòu)建出多層溫度分布模型，運(yùn)用MATLAB 軟件對(duì)模型方程進(jìn)行求解，并得出溫度隨時(shí)間的變化情況；然后結(jié)合Henriques積分提出反問(wèn)題，運(yùn)用有限差分法建立單層厚度優(yōu)化模型，并通過(guò)粒子群算法求解出第Ⅱ?qū)拥淖顑?yōu)厚度為8.8 mm；最后，再運(yùn)用動(dòng)態(tài)規(guī)劃算法（DP）建立基于遺傳算法的優(yōu)化模型[15]。在外界環(huán)境溫度為80 ℃，且Ⅰ層和Ⅲ層厚度固定的前提下對(duì)模型進(jìn)行數(shù)值積分。用MATLAB軟件求解出Ⅱ?qū)拥淖顑?yōu)厚度為11 mm，IV層的最優(yōu)厚度為4.8 mm。