有限差分法的一維熱傳導方程應用

陳金雄，張敏，沈丹梅，楊玲玲，羅翔文

（1.武夷學院 數學與計算機學院，福建 武夷山 354300;2.武夷學院 土木工程與建筑學院，福建 武夷山 354300;3.福州超德中學，福建 福州 350000;4.莆田私立實驗中學，福建 莆田 351100;5.福州靠譜云科技有限公司，福建 福州 350000）

服裝阻止人體皮膚與空氣環(huán)境直接接觸，相當于保護人體的第二層皮膚。當人們從事不同的生產活動時，在不同環(huán)境下人們需要穿著不同性能的服裝。在高溫環(huán)境工作時，人們一般要穿著高溫防護服，而高溫防護服通常由三層織物材料構成，記為I、II、III層，其中I層與外界環(huán)境接觸，III層與皮膚之間還存在空隙，將此空隙記為IV層。高溫防護服的作用原理就是減緩和阻止熱量傳遞，使得到達人體皮膚的熱量盡可能少的聚集，以起到保護皮膚不被燒傷或灼傷的作用[1]。因此，高溫防護服應該能夠起到保護人體的作用，同時也起到促使人體熱量散發(fā)、減小高溫空氣的傳遞速度、防止高溫中暑等作用[1]。而在高溫防護服的設計上要求高溫防護服要有較好的阻燃性和較強的隔熱性，為評價高溫防護服專用材料的防護性能的好壞需要進行大量的熱防護性能的測試，但是，在高溫環(huán)境中進行大量實驗是不現(xiàn)實的，其將需要耗費大量資金并且將造成資源浪費。

研究對多層高溫防護服的熱傳導問題，首先運用熱傳導的基本定律-傅里葉定律[2]建立一維熱傳導方程模型，并利用有限差分法對時間和空間節(jié)點作離散化處理，得到溫度分布的時間空間函數，進而得到最終的溫度分布數值解。接著，建立基于有限差分法的一維熱傳導方程模型并運用二分法得出最優(yōu)近似解數值，即高溫防護服第II層最佳厚度，為高溫防護服的合理設計提供依據。

1 研究問題及其主要思路

研究問題來源于2018年全國大學生數學建競賽A題，問題如下：為了研究高溫防護服的厚度設計問題，現(xiàn)將37℃的假人放置在高溫環(huán)境中并測得假人皮膚外側的溫度，下面我們利用借鑒數學模型[1-7]來確定假人皮膚外側的溫度變化狀況，并解決兩方面的問題：

(1)已知高溫防護服專用測量材料的某些參數值(見表1)，當環(huán)境溫度為T(x,t)=75℃時，各層專用材料厚度分別為LⅠ=0.6 mm、LⅡ=6 mm、LⅢ=3.6 mm、LⅣ=5 mm，作業(yè)時間為0≤t≤90min時，假人皮膚外側溫度T（x=xⅣ，t），建立數學模型，計算溫度分布。

表1 專用服裝材料的參數值Tab.1 Parameters of special clothing materials

(2)當環(huán)境溫度為T(x=0，t)=65℃、織物材料第四層厚度為LⅣ=5.5mm時，確定第II層的最優(yōu)厚度，確保當工作60min時，假人皮膚外側溫度不超過47℃，且超過44℃的時間不超過5分鐘。

首先假設人體的熱量不會向外傳遞，而假人皮膚最外層溫度最終穩(wěn)定在48.08℃，再提出假設:熱量傳遞形式只有熱傳導，并且熱傳遞是沿著垂直于皮膚方向進行的，故視為一維的[3]。其一，針對“環(huán)境-高溫防護服-人體皮膚”系統(tǒng)，利用熱傳導的基本定律-傅里葉定律，建立一維熱傳導方程模型，然后利用有限差分對溫度的空間和時間節(jié)點進行離散化處理。利用MATLAB軟件對問題模型求解得到的溫度分布數據表和已測量的數據進行對比，找到擬合程度最好的一組數據，驗證該厚度結合下的高溫防護服的實際有效性。其二，在一定條件下確定第Ⅱ層的最優(yōu)厚度。假設人體體溫恒定為37℃，邊界條件如下：環(huán)境溫度T(x=0,t)=65℃，右側邊界溫度T(x=xd，t)=37℃，假人的各參數信息是已知的，則可以在問題(1)的模型基礎上利用有限差分法求得假人在不同時刻的體表溫度分布。但實際上題目并沒有給出假人的各項參數，即假人皮膚厚度Ld，熱傳導系數Rd，密度ρd，比熱cd等信息。該問題上可以假設假人皮膚厚度，利用熱阻理論可以求出假人的熱傳導系數;而假人密度不妨設與真人一致，利用二分法原理求出假人皮膚比熱。再在問題(1)的模型基礎上利用有限差分法，確定出第Ⅱ層的最優(yōu)厚度。

2 模型建立與求解

2.1 問題一的模型建立與求解

2.1.1 模型分析

熱傳導過程中溫度隨時間或空間坐標變化而變化，熱傳導過程分為Ⅰ層、Ⅱ層、Ⅲ層、Ⅳ層（空氣層）和皮膚層。因此考慮建立一維熱傳導方程模型，用有限差分法對發(fā)生熱傳導現(xiàn)象的空間和時間節(jié)點進行離散化處理，并建立關于時間和空間節(jié)點的溫度分布函數。

2.1.2 模型建立

將高溫防護服的三層材料平壁，假設熱量穿過層間接觸面時，不會有熱阻作用，溫度不下降。建立“環(huán)境-高溫防護服-人體皮膚”系統(tǒng)。取最表層織物材料上的任一點為原點O，建立平面直角坐標系，過這點豎直方向為y軸正方向、垂直于織物材料和皮膚方向為x軸正方向。圖1為該導熱系統(tǒng)的截面示意圖。

圖1 “環(huán)境-高溫防護服-人體皮膚”系統(tǒng)截面示意圖Fig.1 “Environment-high temperature protective clothing-human skin”system section diagram

熱傳導遵循傅里葉定律[4，5]

式中:q為熱流量，x表示一維空間坐標，t為時間，k為熱傳導系數，T為溫度。

由傅里葉定律知，要計算物體的熱流量，除了要知道導熱系數還要知道物體溫度分布[6]，因此必須建立熱傳導微分方程。

在一維非穩(wěn)態(tài)無熱源的熱傳導過程中，一維無熱源熱傳導方程表示為

式中:ρ是織物材料密度，c是織物材料比熱，k是熱傳導系數，T為溫度，t為時間。一維空間中的ρ、c、k對t而言都是常數，則方程（2）變?yōu)闊醾鲗Х匠?/p>

現(xiàn)應用有限差分法[6-7]對溫度T（x,t）分布的空間和時間節(jié)點進行離散化處理，再求離散點對應的溫度，圖2為離散后的空間和時間節(jié)點示意圖。

圖2 離散后的空間和時間節(jié)點示意圖Fig.2 Discrete sketch of space and time nodes

對連續(xù)方程進行離散化處理，任意取四個整數N，a，b，c，時間步長為h=，三層防護服專用材料所對應的空間步長分別為，則時間節(jié)點和空間節(jié)點分別為

假設T（xi，tj）的數值解為,初始條件=（xi，0）,離散后可得到織物材料層內溫度的有限差分格式如下

根據表1的數據和公式(4)、(5)定義的空間和時間節(jié)點，應用有限差分法求解熱傳導方程(3),通過MATLAB[8]軟件獲得不同時刻的溫度分布圖(如圖3所示),截取了0、1、2、4、8、16、45 min等七個時間的溫度分布圖,如下圖3中（a）、（b）、（c）、（d）、（e）、（f）、（g）所示；最終得到穩(wěn)態(tài)時的溫度分布數值解，如下圖3中（h）所示。此外，根據不同時刻溫度分布值畫出不同時刻的空間溫度分布圖,如下圖4示。

圖3 不同時刻溫度分布圖Fig.3 Temperature distribution at different time

圖4 不同時刻的空間溫度分布Fig.4 Spatial temperature distribution at different time

通過Matlab畫出上面求解得到的溫度與對應時間的關系圖、賽題附件二的假人皮膚外側溫度與時間的關系圖，如圖5所示，其中粗曲線是賽題附件二的溫度曲線，細曲線是求解得到的溫度曲線。

圖5 溫度對比圖Fig.5 Temperature chart

由圖5可以看出得到的假人皮膚外側溫度分布曲線和賽題的附件2中假人皮膚外側溫度曲線擬合良好，并通過Matlab運用插值擬合得到溫度與時間的關系式為

2.2 問題二的模型建立與求解

2.2.1 模型分析

問題一已經得到了假人外側皮膚在不同時刻的溫度分布圖。但實際上題目并沒有給出假人的皮膚厚度Ld、熱傳導系數Rd、密度ρd、比熱cd等信息。該問題上可以假設假人皮膚厚度，利用熱阻理論可以求出假人的熱傳導系數;而假人密度不妨設與真人一致，利用二分法原理求出假人皮膚比熱。再在問題一的基礎上利用有限差分法確定第Ⅱ層的最優(yōu)厚度。

2.2.2 模型求解

2.2.2.1 熱阻

采用類比電流、電壓和電阻間的關系可快速計算最終溫度分布的穩(wěn)態(tài)情況。熱阻系數亦稱為熱阻，是指在有溫度差的情形下，導熱材料的熱阻系數增大，則導熱材料對熱傳導的阻礙能力就增強。因此在一維無內熱源的情況下，對于兩端溫差為ΔT（類比于電壓），熱流密度為qd=k×ΔT/Ld（類比于電流），厚為Ld的導熱層，其熱阻Rd（類比于電阻）定義如下

在問題一中用的是一維熱傳導方程模型，可以看成每一層的熱阻是一樣的，即

2.2.2.2 二分法原理[9]

(1)f（x）的有根開區(qū)間為[a,b]，在兩端點處的函數值為f（a）,f（b）；

(4)重復進行②③步，直到區(qū)間[a,b]的長度小于允許的誤差值，此時中點即為所求根。

2.2.2.3 模型求解

求假人皮膚表面溫度T（xⅣ，t），但必須確定第Ⅱ層的最優(yōu)厚度。假設人體體溫恒定37℃，邊界條件如下：環(huán)境溫度T（x=0,t）=65℃，右側邊界溫度T（x=xd，t）=37℃，假如假人的各參數信息是已知的，則可以在第一問的基礎上利用有限差分法確定第Ⅱ層的最優(yōu)厚度。但實際上題目并沒有給出假人的各項參數即假人皮膚厚度Ld，熱傳導系數kd，密度ρd，比熱cd等信息。該問題上可以假設假人皮膚厚度，利用熱阻理論可以求出假人的熱傳導系數;而假人密度設與真人一致，利用二分法原理求出假人皮膚比熱。再在基于第一問的基礎上利用有限差分法確定第Ⅱ層的最優(yōu)厚度。

高溫防護服專用材料部分參數已知如表1所示。利用MATLAB軟件進行求解得到皮膚熱阻Rd=0.116 1，進而求得假人皮膚比熱cd=1.0×103，最終確定得到第Ⅱ層的最優(yōu)厚度LⅡ=19.090 6mm。

3 結論

在合理假定的前提下，應用傅里葉定律、有限差分法、二分法等理論將高溫防護服的厚度設計問題轉化為高溫防護服熱傳遞模型，利用MATLAB軟件編程解決了該問題，為高溫防護服設計者提供參考。此外，可將該模型應用到皮膚燒傷度預測等領域。