基于蒙特卡羅算法的高溫作業防護服優化設計

王 麗 燕， 李 慶 杰， 李 雁 宙， 馮 丹 亭， 嚴 瀟， 蔡 方 赟

( 1.大連大學 信息工程學院, 遼寧 大連 116622；2.大連大學 環境與化學工程學院, 遼寧 大連 116622 )

0 引 言

隨著社會的發展，越來越多的工作需要在高溫環境下完成，而高溫作業帶有極高的危險性．在高溫環境下工作時，人們需要穿著專用的高溫作業防護服，才能避免高溫灼熱帶來的傷害．盧琳珍[1]基于非線性變系數拋物型方程，建立了多層熱防護服-空氣層-皮膚系統的熱傳遞模型，利用有限差分法對所建模型進行了收斂性分析，給出了各層接觸面之間的關系．但該系統僅停留在理論研究層面，沒有進一步運用到服裝設計上．陳柔羲[2]從人體熱濕生理和運動需求出發，對PCM(相變材料)在高溫防護服內的分布情況進行了優化設計，并依據熱濕舒適性對模型進行了改進．但所選取的樣本數比較少，試驗結果具有偶然性．Ghazy等[3]通過建立一個瞬態傳熱有限體積模型，研究了防護服系統的瞬態溫度分布和能量含量，但忽略了防護服中空氣間隙對防護服性能的影響．毛瑤瑤等[4]在溫度和時間的約束條件下，利用遺傳算法求出了防護服第2層的最優厚度，但遺傳算法不能及時利用網絡的反饋信息，計算效率較低．本文圍繞2018年全國大學生數學建模競賽A題“高溫作業專用服裝設計”(http://www.mcm.edu.cn/html_cn/node/7cec7725b9a0ea07b4dfd175e8042c33.html)，在降低研發成本的條件下，對高溫作業專用服裝的厚度進行優化求解．所研究的高溫作業專用服裝由3層織物材料構成，其中第1層與外界環境接觸，第3層與皮膚之間還存在空隙，將此空隙記為第4層．利用蒙特卡羅算法求解防護服第2層與第4層的最優厚度．在求解過程中，不僅考慮到空氣間隙對防護服性能的影響，而且在全局隨機取樣，避免陷入局部最優解．

1 防護服溫度分布的求解

由文獻[5]可知，當皮膚外側的溫度不再改變時，整個系統達到穩態，即物體各點的溫度不再隨時間變動，而之前時間段皮膚外側傳熱過程均為非穩態導熱過程．因此，在求解關系式的過程中，要考慮至少兩個狀態下的導熱情況．以后的所有分析都基于這一前提．本章通過求解熱傳導偏微分方程[5]，建立離散型熱傳導模型，對高溫作業防護服仿真設計，求解每一層的溫度分布．

1.1 熱傳導方程離散化

將高溫作業防護服進行一維切片化處理后得到高溫作業防護服熱傳導示意圖，如圖1所示．

圖1 高溫作業防護服熱傳導三維等溫示意圖

Fig.1 Three-dimensional isothermal diagram of heat conduction of high-temperature operation protective clothing

從圖1可以直觀看出溫度的傳導方向和傳導過程．熱傳導方程式[6]如下：

(1)

其中ρ、Cρ均為常數，由于每一層的熱傳導率不同且是離散的，式(1)中k值不是常數．對熱傳導方程(1)進行離散化處理，離散化處理過程如圖2所示．

圖2 離散化處理示意圖

圖中的Δx和Δt分別表示厚度的步長和時間的步長．i表示離散點的厚度軸坐標，j表示離散點的時間軸坐標．則由式(1)可得

(2)

(3)

將式(2)、(3)同時代入熱傳導方程(1)中，得到離散化的熱傳導方程為

(4)

在將偏微分方程(1)離散化處理之前，應先設置初值，即各層防護服在t=0時均為一常數值，初值設立如下：

T1(x，0)=348

(5)

T2(x，0)=310

(6)

T3(x，0)=310

(7)

T4(x，0)=310

(8)

在求解過程中，經過離散化熱傳導方程的不斷迭代，可得到每一層的、每一時間點上的各個溫度值．

1.2 防護服各層溫度分布的求解

分別將密度、比熱容、熱傳導率、環境溫度和假人外表皮溫度代入式(4)中．假設橫軸達到0.6 mm 時作為第1層，橫軸為0.6～6.6 mm時作為第2層，橫軸為6.6～10.2 mm時作為第3層，橫軸為10.2～15.2 mm時作為第4層．用MATLAB求解式(4)，可得到每一層材料各處的溫度值．

在此基礎上，基于得到的每一層材料各處、各時間點對應的溫度值，找出相鄰兩層之間不同時間點對應的溫度值并作圖，即得到每兩層間溫度隨時間變化的二維平面圖，如圖3所示．

圖3 每一層溫度隨時間變化的分布示意圖

從圖3中可以看出：溫度從第1層開始到第4層，先后逐漸趨于穩態，并且在傳遞熱量的過程中，能量被損耗，最終使得熱量的傳遞值遞減，每層的穩態溫度都低于其上一層．為方便看出各層處于穩態時的溫度，做出5 400 s時的各層溫度情況，如圖4所示．

圖4 各層達到穩態時的溫度分布情況圖

由圖4可以看出，第4層空氣層溫度下降最多，說明空氣層對人體隔熱防護更為重要．

1.3 誤差分析

將各層溫度數據代入均方根誤差[7]公式：

(9)

對溫度分布結果進行分析．把用MATLAB計算得出的溫度分布值與測量得到的真實值進行對比，得到表1．

表1 誤差分析表

由此得出均方差σ=1.021 0，誤差較小．因此，本文所建立的溫度分布模型合理、可行．

2 基于蒙特卡羅算法求解防護服最優厚度

求解最優厚度的方法有很多，例如窮舉法、模擬退火、遺傳算法等．窮舉法雖然方法簡單，但是實施起來過于繁瑣復雜；啟發式算法雖然常用，但是容易陷入局部最優解．對于本文所研究的最優厚度求解來說，蒙特卡羅算法不僅方法簡單，而且在求解過程中，通過隨機投點，可以避免陷入局部最優解，使結論更具有說服力．因此，本文采用蒙特卡羅算法，求解防護服的最優厚度．

假設邊界空氣層溫度T0=353 K，假人皮膚外側溫度約束條件有

(10)

在滿足約束條件的情況下，為了降低研發成本和縮短研發周期，服裝的總厚度值應該最小，即

min(d)=min(d1+d2+d3+d4)

(11)

根據蒙特卡羅算法，將厚度進行隨機投點，則由d2和d4生成n個隨機數(其中d2的范圍為0.6～25 mm，d4的范圍為0.6～6.4 mm)，代入式(4)并求出其溫度分布．在滿足約束條件的情況下，隨著每次d2和d4數值的改變，會生成不同的目標函數，同時記錄當前最優解．當投點次數n取100，101，102，103，104，…時得到材料厚度收斂圖，如圖5所示．

(a) 第2層

(b) 第4層

圖5 厚度收斂圖像

Fig.5 Thickness convergence image

當t=1 800 s時，溫度不斷向320 K逼近；當t=1 500 s時，溫度不斷向317 K逼近．第2層和第4層厚度隨時間變化不斷向最優厚度收斂，此時可以得到第2層的最優厚度值收斂于18.68 mm左右，第4層的最優厚度值收斂于6.4 mm左右．

由于蒙特卡羅算法生成的隨機數帶有一定的隨機性，每次的結果都具有一定的誤差，本文選擇多次測量取平均值的方法，確定最終的最優厚度．而上文所提到的投點次數n的選取會嚴重影響代碼的計算時間，在保證準確性并且縮短計算時間的前提下，經過多次試驗，當選取投點次數n為1×106時，各提取一組樣本(共7個)如表2、3所示．

表2 迭代次數為1×106時第2層的厚度

表3 迭代次數為1×106時第4層的厚度

將樣本中的最優厚度求算術平均值，計算得出第2層的最優厚度為18.68 mm，第4層的最優厚度為6.38 mm．此時滿足當邊界空氣層溫度為353 K 時，確保假人皮膚外側溫度在工作30 min時，不超過320 K，且超過317 K的時間不超過5 min．

3 結 語

在利用傅里葉熱傳導定律求解防護服各層的溫度分布時，由于每一層的熱傳導率不同且是離散的，本文通過穩態熱傳導方程、傳熱學原理等構建了一維穩態熱傳導模型．利用有限差分法求解出人體皮膚外側每一層織物材料的溫度分布情況，對所求得的結果進行了誤差分析，驗證了模型的可靠性．并基于以上模型，對防護服第2層與第4層的最優厚度進行了求解．借助MATLAB軟件，運用蒙特卡羅算法，先對厚度值進行不同次數的隨機模擬取點，在滿足約束條件的情況下，求解出不同次數下對應的厚度值，然后對厚度值進行收斂分析，發現當投點次數為1×106時，收斂曲線趨于平滑，此時計算結果精確度較高，計算時間較短．由于蒙特卡羅算法獲得的結果具有一定的隨機性，在投點次數為1×106的條件下獲得7組樣本數據，然后通過樣本均值排除誤差，得到防護服第2層與第4層的最優厚度．此時的最優厚度為防護服滿足約束條件下的最小值，降低了研發成本，縮短了研發周期．但是為了方便求解，沒有考慮人體皮膚表面的汗液、水汽等因素的影響，因此，模型還有較大的優化空間．但可為相關高溫熱防護服設計及其衍生產業發展提供一定的參考．