有限差分法的一维热传导方程应用

2022-05-21陈金雄张敏沈丹梅杨玲玲罗翔文

武夷学院学报 2022年3期

陈金雄，张敏，沈丹梅，杨玲玲，罗翔文

（1.武夷学院数学与计算机学院，福建武夷山 354300;2.武夷学院土木工程与建筑学院，福建武夷山 354300;3.福州超德中学，福建福州 350000;4.莆田私立实验中学，福建莆田 351100;5.福州靠谱云科技有限公司，福建福州 350000）

服装阻止人体皮肤与空气环境直接接触，相当于保护人体的第二层皮肤。当人们从事不同的生产活动时，在不同环境下人们需要穿着不同性能的服装。在高温环境工作时，人们一般要穿着高温防护服，而高温防护服通常由三层织物材料构成，记为I、II、III层，其中I层与外界环境接触，III层与皮肤之间还存在空隙，将此空隙记为IV层。高温防护服的作用原理就是减缓和阻止热量传递，使得到达人体皮肤的热量尽可能少的聚集，以起到保护皮肤不被烧伤或灼伤的作用[1]。因此，高温防护服应该能够起到保护人体的作用，同时也起到促使人体热量散发、减小高温空气的传递速度、防止高温中暑等作用[1]。而在高温防护服的设计上要求高温防护服要有较好的阻燃性和较强的隔热性，为评价高温防护服专用材料的防护性能的好坏需要进行大量的热防护性能的测试，但是，在高温环境中进行大量实验是不现实的，其将需要耗费大量资金并且将造成资源浪费。

研究对多层高温防护服的热传导问题，首先运用热传导的基本定律-傅里叶定律[2]建立一维热传导方程模型，并利用有限差分法对时间和空间节点作离散化处理，得到温度分布的时间空间函数，进而得到最终的温度分布数值解。接着，建立基于有限差分法的一维热传导方程模型并运用二分法得出最优近似解数值，即高温防护服第II层最佳厚度，为高温防护服的合理设计提供依据。

1 研究问题及其主要思路

研究问题来源于2018年全国大学生数学建竞赛A题，问题如下：为了研究高温防护服的厚度设计问题，现将37℃的假人放置在高温环境中并测得假人皮肤外侧的温度，下面我们利用借鉴数学模型[1-7]来确定假人皮肤外侧的温度变化状况，并解决两方面的问题：

(1)已知高温防护服专用测量材料的某些参数值(见表1)，当环境温度为T(x,t)=75℃时，各层专用材料厚度分别为LⅠ=0.6 mm、LⅡ=6 mm、LⅢ=3.6 mm、LⅣ=5 mm，作业时间为0≤t≤90min时，假人皮肤外侧温度T（x=xⅣ，t），建立数学模型，计算温度分布。

表1 专用服装材料的参数值Tab.1 Parameters of special clothing materials

(2)当环境温度为T(x=0，t)=65℃、织物材料第四层厚度为LⅣ=5.5mm时，确定第II层的最优厚度，确保当工作60min时，假人皮肤外侧温度不超过47℃，且超过44℃的时间不超过5分钟。

首先假设人体的热量不会向外传递，而假人皮肤最外层温度最终稳定在48.08℃，再提出假设:热量传递形式只有热传导，并且热传递是沿着垂直于皮肤方向进行的，故视为一维的[3]。其一，针对“环境-高温防护服-人体皮肤”系统，利用热传导的基本定律-傅里叶定律，建立一维热传导方程模型，然后利用有限差分对温度的空间和时间节点进行离散化处理。利用MATLAB软件对问题模型求解得到的温度分布数据表和已测量的数据进行对比，找到拟合程度最好的一组数据，验证该厚度结合下的高温防护服的实际有效性。其二，在一定条件下确定第Ⅱ层的最优厚度。假设人体体温恒定为37℃，边界条件如下：环境温度T(x=0,t)=65℃，右侧边界温度T(x=xd，t)=37℃，假人的各参数信息是已知的，则可以在问题(1)的模型基础上利用有限差分法求得假人在不同时刻的体表温度分布。但实际上题目并没有给出假人的各项参数，即假人皮肤厚度Ld，热传导系数Rd，密度ρd，比热cd等信息。该问题上可以假设假人皮肤厚度，利用热阻理论可以求出假人的热传导系数;而假人密度不妨设与真人一致，利用二分法原理求出假人皮肤比热。再在问题(1)的模型基础上利用有限差分法，确定出第Ⅱ层的最优厚度。

2 模型建立与求解

2.1 问题一的模型建立与求解

2.1.1 模型分析

热传导过程中温度随时间或空间坐标变化而变化，热传导过程分为Ⅰ层、Ⅱ层、Ⅲ层、Ⅳ层（空气层）和皮肤层。因此考虑建立一维热传导方程模型，用有限差分法对发生热传导现象的空间和时间节点进行离散化处理，并建立关于时间和空间节点的温度分布函数。

2.1.2 模型建立

将高温防护服的三层材料平壁，假设热量穿过层间接触面时，不会有热阻作用，温度不下降。建立“环境-高温防护服-人体皮肤”系统。取最表层织物材料上的任一点为原点O，建立平面直角坐标系，过这点竖直方向为y轴正方向、垂直于织物材料和皮肤方向为x轴正方向。图1为该导热系统的截面示意图。

图1 “环境-高温防护服-人体皮肤”系统截面示意图Fig.1 “Environment-high temperature protective clothing-human skin”system section diagram

热传导遵循傅里叶定律[4，5]

式中:q为热流量，x表示一维空间坐标，t为时间，k为热传导系数，T为温度。

由傅里叶定律知，要计算物体的热流量，除了要知道导热系数还要知道物体温度分布[6]，因此必须建立热传导微分方程。

在一维非稳态无热源的热传导过程中，一维无热源热传导方程表示为

式中:ρ是织物材料密度，c是织物材料比热，k是热传导系数，T为温度，t为时间。一维空间中的ρ、c、k对t而言都是常数，则方程（2）变为热传导方程

现应用有限差分法[6-7]对温度T（x,t）分布的空间和时间节点进行离散化处理，再求离散点对应的温度，图2为离散后的空间和时间节点示意图。

图2 离散后的空间和时间节点示意图Fig.2 Discrete sketch of space and time nodes

对连续方程进行离散化处理，任意取四个整数N，a，b，c，时间步长为h=，三层防护服专用材料所对应的空间步长分别为，则时间节点和空间节点分别为

假设T（xi，tj）的数值解为,初始条件=（xi，0）,离散后可得到织物材料层内温度的有限差分格式如下

根据表1的数据和公式(4)、(5)定义的空间和时间节点，应用有限差分法求解热传导方程(3),通过MATLAB[8]软件获得不同时刻的温度分布图(如图3所示),截取了0、1、2、4、8、16、45 min等七个时间的温度分布图,如下图3中（a）、（b）、（c）、（d）、（e）、（f）、（g）所示；最终得到稳态时的温度分布数值解，如下图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等信息。该问题上可以假设假人皮肤厚度，利用热阻理论可以求出假人的热传导系数;而假人密度不妨设与真人一致，利用二分法原理求出假人皮肤比热。再在问题一的基础上利用有限差分法确定第Ⅱ层的最优厚度。

2.2.2 模型求解

2.2.2.1 热阻

采用类比电流、电压和电阻间的关系可快速计算最终温度分布的稳态情况。热阻系数亦称为热阻，是指在有温度差的情形下，导热材料的热阻系数增大，则导热材料对热传导的阻碍能力就增强。因此在一维无内热源的情况下，对于两端温差为ΔT（类比于电压），热流密度为qd=k×ΔT/Ld（类比于电流），厚为Ld的导热层，其热阻Rd（类比于电阻）定义如下

在问题一中用的是一维热传导方程模型，可以看成每一层的热阻是一样的，即

2.2.2.2 二分法原理[9]

(1)f（x）的有根开区间为[a,b]，在两端点处的函数值为f（a）,f（b）；

(4)重复进行②③步，直到区间[a,b]的长度小于允许的误差值，此时中点即为所求根。

2.2.2.3 模型求解

求假人皮肤表面温度T（xⅣ，t），但必须确定第Ⅱ层的最优厚度。假设人体体温恒定37℃，边界条件如下：环境温度T（x=0,t）=65℃，右侧边界温度T（x=xd，t）=37℃，假如假人的各参数信息是已知的，则可以在第一问的基础上利用有限差分法确定第Ⅱ层的最优厚度。但实际上题目并没有给出假人的各项参数即假人皮肤厚度Ld，热传导系数kd，密度ρd，比热cd等信息。该问题上可以假设假人皮肤厚度，利用热阻理论可以求出假人的热传导系数;而假人密度设与真人一致，利用二分法原理求出假人皮肤比热。再在基于第一问的基础上利用有限差分法确定第Ⅱ层的最优厚度。

高温防护服专用材料部分参数已知如表1所示。利用MATLAB软件进行求解得到皮肤热阻Rd=0.116 1，进而求得假人皮肤比热cd=1.0×103，最终确定得到第Ⅱ层的最优厚度LⅡ=19.090 6mm。