盧 鵬，徐昌貴

（西南交通大學(xué)數(shù)學(xué)學(xué)院，四川成都610031）

人們?cè)诟邷丨h(huán)境下工作時(shí)，需要穿著專用服裝以免灼傷.專用服裝由三層織物材料構(gòu)成，即I、II、III層，其中I 層與外界環(huán)境接觸，III 層與皮膚之間還存在空隙，記為IV 層.在設(shè)計(jì)專用服裝時(shí)，為了降低研發(fā)成本、縮短研發(fā)周期，將體內(nèi)溫度控制在37°C 的假人放置在實(shí)驗(yàn)室的高溫環(huán)境中，測(cè)量假人皮膚外側(cè)的溫度. 專用服裝材料的一些參數(shù)值如表1 所示.對(duì)環(huán)境溫度為75 °C，工作時(shí)間為90 分鐘的情形開展實(shí)驗(yàn)，測(cè)量得到假人皮膚外側(cè)的溫度（見表2，完成數(shù)據(jù)可查文獻(xiàn)[1]）. 本文首先對(duì)環(huán)境―人體―服裝傳熱的物理過程進(jìn)行推理和分析，綜合傳熱學(xué)知識(shí)[2-3]將服裝導(dǎo)熱過程簡(jiǎn)化為一維平板熱傳導(dǎo)模式，接著用微元法[4]推導(dǎo)出一維熱傳導(dǎo)偏微分方程，然后從材料的均勻性得到了不同材料在分界面處溫度連續(xù)、熱量連續(xù)等條件，并根據(jù)實(shí)驗(yàn)要求確定初始條件和邊值條件，由此建立了在穩(wěn)定環(huán)境中專用服裝熱傳導(dǎo)的偏微分方程組模型. 通過對(duì)模型進(jìn)行求解[5-7]，得出每層傳播的規(guī)律. 最后對(duì)結(jié)果進(jìn)行了驗(yàn)證，以期可用于指導(dǎo)服裝的生產(chǎn)與研制.

表1 專用服裝材料的參數(shù)值

表2 假人皮膚外側(cè)的測(cè)量溫度

1 問題分析

1.1 一維化思想的分析

由于熱量在傳遞過程中，最終熱量的對(duì)外表現(xiàn)形式為溫度場(chǎng)函數(shù)T(x,y,z,t)，而在本文中，高溫防護(hù)服裝材料的各處均在進(jìn)行熱量的傳導(dǎo)，對(duì)任意一個(gè)熱量微元，由于材料各向分布均勻，即各個(gè)熱元傳熱的速度是互相協(xié)調(diào)的，因此微元在垂直向的熱量變化是相同的，即在圖1 中Y方向上溫度處處相同，由此可以認(rèn)為熱量是沿一個(gè)方向傳遞，熱元的側(cè)面之間沒有進(jìn)行熱量交換，即二維溫度分布可進(jìn)一步表示為T(x,t)，即溫度為時(shí)間和單維度位置變量的函數(shù).

圖1 熱量微元的選取示意圖

1.2 一維熱傳導(dǎo)方程的推導(dǎo)

（1）熱量流動(dòng)速率的確定

在對(duì)服裝材料上任意位置取一個(gè)某種材料的圓柱體微元，該材料微元的橫截面均勻且材質(zhì)相同，熱量只能沿某特定方向進(jìn)行流動(dòng)，此處方向約定為x軸方向. 根據(jù)經(jīng)驗(yàn)原則，在某時(shí)刻t，通過點(diǎn)x橫截面上的單位面積熱流動(dòng)速率為[8]：

其中：λ為該材料的熱傳導(dǎo)率，v(x,t)為熱流速率.

圖2是熱量流入材料微元的示意圖.

圖2 熱量流入材料微元的示意圖

由圖2可知，通過兩個(gè)端面的熱量流動(dòng)速率為：

將v(x,t)代入得：

（2）單種材料微元吸收的熱量[9-10]

把該材料圓柱微元溫度從0 加熱到T所需的熱量為：

對(duì)上式兩邊進(jìn)行微分并應(yīng)用積分中值定理有：

其中為x與x+Δx之間中的某個(gè)值.

（3）一維熱傳導(dǎo)方程的確定

Q(t)表示把圓柱微元對(duì)應(yīng)于x與x+Δx的一小段溫度從0 度升到已知溫度T(x,t)所需的熱量，又因?yàn)闊崃鲀H從端面進(jìn)出，由此根據(jù)V=dQ/dt和熱量流動(dòng)V表達(dá)式得：

基于上述兩種不同的熱量速度推導(dǎo)方法得出的V相等，于是有：

所以，可以得出

令Δx→0，即xˉ→x，可以得到熱量的一維傳導(dǎo)方程：

2 模型建立

2.1 坐標(biāo)系建立及熱量傳遞機(jī)理分析

在忽略沿假人皮膚平面方向?qū)岬臈l件下，可把皮膚內(nèi)的導(dǎo)熱問題認(rèn)為是無限大平板問題，即一維導(dǎo)熱問題. 由題目可知，在高溫下工作的人體皮膚溫度必須在一定的范圍內(nèi)，人體才能正常工作，而人體皮膚溫度升高是由于環(huán)境溫度向內(nèi)傳遞，才不斷導(dǎo)致溫度越來越高，而高溫防護(hù)服的特點(diǎn)就在于能夠有效地減少熱量向人體傳遞，控制服裝內(nèi)部的溫度，維持人體舒適溫度，而熱量向內(nèi)傳遞是一個(gè)逐步傳遞過程，通過一層層防護(hù)服材料的過程.

建立如圖3 所示坐標(biāo)系，根據(jù)前述分析，坐標(biāo)系的兩個(gè)維度為材料的位置坐標(biāo)和時(shí)間. 溫度分布是坐標(biāo)位置和時(shí)間的二維函數(shù)，坐標(biāo)軸的原點(diǎn)位置及熱量傳遞過程如圖3所示.

圖3 熱量傳遞過程示意圖

2.2 各織物材料層內(nèi)一維熱傳導(dǎo)偏微分方程的建立

熱量在特定的某種材料中傳遞符合一維熱傳導(dǎo)規(guī)律，由此，建立任意層中熱量傳遞的偏微分方程.第i層織物材料中的一維熱傳導(dǎo)偏微分方程為：

其中：Ti是第i層織物材料中的溫度場(chǎng)函數(shù)，為二元函數(shù)Ti(x,t)，x是織物材料在厚度方向的位置坐標(biāo).

2.3 各織物材料在分界面處傳導(dǎo)模式分析

能量在同種介質(zhì)中按照一維熱傳導(dǎo)偏微分方程的遞變規(guī)律進(jìn)行傳遞，根據(jù)熱量傳遞的連續(xù)性特征可知，溫度向下傳遞時(shí)，在分界面上二者的溫度應(yīng)該相等，則有如下等式：

同時(shí)在交界面處熱流密度也是連續(xù)變化的，從而有如下等式：

2.4 邊值條件的確定

由高溫防護(hù)服的工作環(huán)境可知，服裝最外層面料是直接暴露在空氣中的，最里層即第IV 層空隙與人體皮膚接觸，因此在溫度場(chǎng)函數(shù)中，當(dāng)位置坐標(biāo)確定后有如下邊界條件：

（?。┛椢锊牧蠈拥淖筮吔鐥l件：T1(0,t)=T0.

（ⅱ）織物材料層的右邊界條件：T4(l4,t)=T5.

其中：T0=75 ℃為外界環(huán)境的溫度；T5為人體皮膚外側(cè)溫度且數(shù)據(jù)由表2 給出；di表示第i層織物材料的厚度，li表示前i層織物材料的厚度和，即l4=d1+d2+d3+d4.

2.5 初值條件的確定

高溫防護(hù)服裝在進(jìn)入實(shí)驗(yàn)室之前，應(yīng)該長(zhǎng)期處于一個(gè)穩(wěn)定環(huán)境場(chǎng)中，同時(shí)由于假人體內(nèi)溫度一直控制為37 ℃，故有初值條件為：

2.6 偏微分方程組模型的建立

基于上述熱量傳遞的分析過程以及相關(guān)的定解條件，可以建立起整個(gè)熱量傳導(dǎo)的正向偏微分方程組模型：

3 模型求解

在偏微分方程領(lǐng)域，這是一個(gè)一維拋物型方程，考慮到定解條件中既包含開始時(shí)刻織物材料層溫度分布的初始條件，又包含織物材料層溫度分布的邊值條件. 因此，這是一個(gè)一維熱傳導(dǎo)混合問題[11]. 由于混合問題的求解過程往往是比較復(fù)雜的，有時(shí)很難有解析解. 即使該一維熱傳導(dǎo)混合問題的解析解可以求出，其解的形式也通常是一個(gè)無窮級(jí)數(shù)形式，不易于對(duì)外界環(huán)境―服裝―人體整個(gè)系統(tǒng)熱量傳導(dǎo)及溫度變化過程進(jìn)行分析. 因此，本文利用有限差分法對(duì)該偏微分方程組進(jìn)行求解.

采用一定的網(wǎng)格劃分離散化溫度場(chǎng)，把實(shí)際連續(xù)的溫度場(chǎng)離散為有限數(shù)量的點(diǎn). 用這些點(diǎn)上的溫度值近似描述連續(xù)的溫度場(chǎng). 針對(duì)x-t平面，取Δx和Δt分別為x方向和t方向的步長(zhǎng)，分別做一組平行于x軸和平行于t軸的直線，將x坐標(biāo)等分為N份，將t坐標(biāo)等分為M份. 則坐標(biāo)可表示為(kΔx,jΔt)(k=1,2,…,N;j=1,2,…,M)的網(wǎng)格點(diǎn)，簡(jiǎn)記為(k,j)，網(wǎng)格上每個(gè)格點(diǎn)對(duì)應(yīng)一個(gè)溫度，則Ti(x,t)在此格點(diǎn)上的取值記為Ti(k,j)，如圖4所示.

圖4 網(wǎng)格劃分示意圖

3.1 顯式差分算法

用二階中心差分近似代替溫度分布函數(shù)對(duì)空間的偏微分，即

用向前差分近似代替溫度分布函數(shù)對(duì)時(shí)間的偏微分，即

根據(jù)上述建立的差分格式，可將該偏微分方程改寫為：

化簡(jiǎn)為：

算法的計(jì)算過程是由Ti(k-1,j)、Ti(k,j)、Ti(k+1,j)三個(gè)點(diǎn)去計(jì)算Ti(k,j+1)，具體如圖5所示：

圖5 顯式差分格式

這樣根據(jù)上面公式，每次算完一行以后往上迭代就可以算出網(wǎng)格中所有點(diǎn)的值. 但需要注意是總共有四層，每層的ai取值不同；顯式方法在時(shí)，為數(shù)值穩(wěn)定且收斂；所以在時(shí)間步長(zhǎng)與位移步長(zhǎng)取值時(shí)需要滿足條件；為了滿足收斂，每個(gè)部分的步長(zhǎng)有可能不同；材料分界面上的網(wǎng)格節(jié)點(diǎn)(k,j)使用下面方程進(jìn)行計(jì)算（溫度相等，熱量相等）：

離散化第二個(gè)方程，化簡(jiǎn)可得

其中νi=λi/Δxi.

若要用顯式方式求解，時(shí)間步長(zhǎng)要取為0.01 s，位移步長(zhǎng)I、II、III層取為0.1 mm，IV 層取為1 mm，再按照上述公式計(jì)算可得結(jié)果.

3.2 隱式差分算法

用二階中心差分近似代替溫度分布函數(shù)對(duì)空間的偏微分，即

用向后差分近似代替溫度分布函數(shù)對(duì)時(shí)間的偏微分，即

根據(jù)上述建立的差分格式，可將該偏微分方程改寫為：

化簡(jiǎn)為：

計(jì)算過程是由Ti(k-1,j+1)、Ti(k,j+1)、Ti(k+1,j+1)與Ti(k,j)建立一個(gè)線性方程組進(jìn)行計(jì)算，具體如圖6所示：

圖6 隱式差分格式

根據(jù)上面公式，每次求解一個(gè)聯(lián)立線性方程組后，再往上迭代就可以算出網(wǎng)格中所有點(diǎn)的值. 具體如下：

第一步：當(dāng)j=0 時(shí)求解如下線性方程組，得出Ti(k,1)(i,=1,2,3,4;k=1,2,…,N-1)：

其中：

若規(guī)模較大的線性方程組直接采用高斯消元法求解，對(duì)于微型計(jì)算機(jī)來說運(yùn)算量與存儲(chǔ)量難以承受，本文根據(jù)方程組系數(shù)矩陣的特點(diǎn)，采用追趕法（需要滿足一定條件，本題已滿足）求解，大大節(jié)省了計(jì)算時(shí)間.

第二步：當(dāng)j=1,2,…,M-1 時(shí)求解M-1 個(gè)線性方程組，得出所有Ti.

需要注意的是總共有四層，每層的ai取值不同，在求解線性方程組時(shí)注意邊界條件；隱式方法不論ai的大小，都數(shù)值穩(wěn)定且收斂，但計(jì)算量會(huì)較顯式方法要大，因?yàn)槊壳斑M(jìn)一個(gè)時(shí)間間隔，就需要求解一個(gè)聯(lián)立的線性方程組；方法本身收斂，每個(gè)部分在取步長(zhǎng)時(shí)可以相同，即Δx=Δxi；材料分界面上的網(wǎng)格節(jié)點(diǎn)(k,j)由下面方程給出：

離散化第二個(gè)方程，化簡(jiǎn)可得

其中νi=λi/Δxi.

3.3 結(jié)果說明

為了計(jì)算方便及穩(wěn)定，采用隱式差分算法進(jìn)行計(jì)算，時(shí)間步長(zhǎng)取為1 s，四層的位移步長(zhǎng)都取為0.1 mm. MATLAB 軟件[12]編程計(jì)算偏微分方程組即可得到溫度分布數(shù)據(jù)表格（表3）.

表3 部分溫度分布數(shù)據(jù)

溫度分布函數(shù)Ti(x,t)是時(shí)間參數(shù)t和空間參數(shù)x的二元函數(shù)，繪制溫度分布三維圖（圖7）.由圖7可見，初始時(shí)刻服裝材料各處溫度均為37 ℃. 假人進(jìn)入實(shí)驗(yàn)室高溫環(huán)境后，服裝材料溫度迅速上升，其中I層材料溫度變化幅度最大，因?yàn)槠渑c外界環(huán)境直接接觸. 針對(duì)同一時(shí)間，隨著坐標(biāo)位置的增加，節(jié)點(diǎn)處溫度上升幅度逐漸減小；針對(duì)同一節(jié)點(diǎn)，隨著時(shí)間增加，該節(jié)點(diǎn)處溫度逐漸上升，但溫度變化率逐漸降低，最后所有的節(jié)點(diǎn)溫度趨于某一穩(wěn)定溫度值.

圖7 溫度分布三維圖

圖8 反映了不同時(shí)刻下織物材料內(nèi)部溫度的變化. 由圖8 可以看出，隨著時(shí)間的增加，各個(gè)位置處網(wǎng)格節(jié)點(diǎn)的溫度在逐漸趨于穩(wěn)定，即熱傳遞達(dá)到了穩(wěn)定狀態(tài)；且在穩(wěn)定狀態(tài)下，各個(gè)網(wǎng)格節(jié)點(diǎn)的溫度隨位置坐標(biāo)呈線性變化.

圖8 不同時(shí)刻下織物材料內(nèi)部溫度的變化

圖9 的分界面處織物材料溫度隨時(shí)間的變化表明，位于分界處的網(wǎng)格節(jié)點(diǎn)，隨著時(shí)間的增加，這些節(jié)點(diǎn)上的溫度均有一個(gè)快速增長(zhǎng)的過程，但在該過程中溫度變化的速率逐漸變緩，最后，節(jié)點(diǎn)上的溫度都達(dá)到了穩(wěn)定狀態(tài).

圖9 分界面處織物材料溫度隨時(shí)間的變化

4 模型驗(yàn)證

由表2數(shù)據(jù)和圖9可以發(fā)現(xiàn)，當(dāng)假人放入穩(wěn)定的環(huán)境中，從1 645 s以后溫度的變化將趨于穩(wěn)定，這個(gè)時(shí)候溫度和時(shí)間就沒有關(guān)系了，只和位移有關(guān)系.故可以建立穩(wěn)定狀態(tài)下的數(shù)學(xué)模型，并用來檢驗(yàn)非穩(wěn)態(tài)模型所算出的結(jié)果.

穩(wěn)定狀態(tài)下的數(shù)學(xué)模型如下：

求解可得：

由原模型計(jì)算可知，當(dāng)1 645 s后溫度趨于穩(wěn)定，表4 給出了穩(wěn)定狀態(tài)和原模型1 645 s后不同位置上的溫度，同時(shí)也給出了兩種結(jié)果的散點(diǎn)圖（圖10）.從表4 和圖10 可以看出，兩種方式的計(jì)算結(jié)果幾乎沒有差別，所以本文所建立的偏微分方程組模型以及求解的結(jié)果是正確的.

表4 兩種模型溫度結(jié)果對(duì)比

圖10 兩種模型的溫度變化曲線

5 結(jié)語(yǔ)

本文建立了在穩(wěn)定環(huán)境中專用服裝熱傳導(dǎo)的偏微分方程組模型，由此可以掌握溫度變化的規(guī)律，同時(shí)改變材料中的某些參數(shù)值，運(yùn)用模型和算法可得更多溫度變化的結(jié)果，對(duì)這些結(jié)果進(jìn)行研究與分析，可以為生產(chǎn)專用服裝的科研人員提供參考，以使其降低研發(fā)成本，提高制作工藝.