基于分层差分方程的热防护服热量分布规律的模拟 Simulation of Heat Distribution Law of Thermal Protective Clothing Based on Lay-ered Difference Equation

doi:10.12677/MOS.2019.81002

Modeling and Simulation
Vol. 08 No. 01 ( 2019 ), Article ID: 28404 , 8 pages
10.12677/MOS.2019.81002

Simulation of Heat Distribution Law of Thermal Protective Clothing Based on Layered Difference Equation

Bixuan Li¹, Kexin Zhu¹, Xuanlin Chen¹, Fang Wang^1,2*

●How to Cite this Article

¹College of Information Science and Technology, Hunan Agricultural University, Changsha Hunan

²Agricultural Mathematical Model and Data Process Center, Hunan Agricultural University, Changsha Hunan

Received: Dec. 17^th, 2018; accepted: Jan. 2^nd, 2019; published: Jan. 10^th, 2019

ABSTRACT

In order to study the heat distribution of heat-resistant garments in high-temperature environments and further to reduce the production cost of heat-resistant garments under high-temperature environment, and shorten the development cycle, in this paper, we first establish a model of high-temperature protective clothing-air-skin for different material layers. The heat conduction equation is used to study the heat distribution of different layers, where the thermal diffusivity in the heat transfer equation depends on the parameter values of the dielectric material. We propose two models to solve the boundary condition problem in partial differential equations. The concept of thermal resistance is introduced in Model I. According to the Crank-Nicolson implicit scheme for solving partial differential equations in numerical analysis, the temperature-space-time correlation results in different material layers can be obtained. In Model II, four different material layers are considered as a whole. Thus, the classical heat differential equations are used to solve the partial differentials of the heat conduction equations for the four layers. Numerical solution of the equation shows that the methods are workable.

Keywords:Crank-Nicolson Implicit Format Method, Classical Explicit Difference Method, Hierarchical Difference Equation

基于分层差分方程的热防护服热量分布规律的模拟

李碧璇¹，朱可馨¹，陈宣霖¹，王访^1,2*

¹湖南农业大学信息科学技术学院，湖南长沙

²湖南农业大学农业数学建模与数据处理中心，湖南长沙

收稿日期：2018年12月17日；录用日期：2019年1月2日；发布日期：2019年1月10日

摘要

为研究耐热服装在高温环境中的热量分布情况，降低在高温环境工作下耐热服装的制作成本，且缩短其研发周期，本文首先建立对于不同材质层的“高温作业防护服–空气–皮肤”热传导方程来研究不同层的热量分布情况，其中热传导方程中热扩散率依据介质材料的参数值而定。我们利用两种方法来解决偏微分方程的边界条件问题。在方法一中引入热阻的概念，根据数值分析中求解偏微分方程的Crank-Nicolson隐式格式，即可求出不同材质层中温度–空间–时间的相关结果。在模型二中，考虑将四种不同材质层视为一个整体，分别对四种不同材质层中热传导方程运用古典显式差分法迭代进行求解偏微分方程数值解。结果表明我们的方法是可行的。

关键词 :Crank-Nicolson隐式格式法，古典显式差分法，分层差分方程

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

1. 引言

随着我国工业生产的繁荣发展，消防、化工、冶金等行业的高温作业问题也日趋严峻，长期在高温和强烈的辐射环境下工作有许多潜在的安全隐患，如热中暑、烧伤和灼伤等危害直接影响着人体的健康安全。高温工作防护服装作为人体与高温环境之间的一层屏障，具有一定程度的隔热和避免灼烧伤功能，是高温作业环境中一项重要的防护措施。它的防护原理是降低热转移速度，让外界高温度的热量缓慢而少量的转移至皮肤。在现今的热防护服的研究领域中，亟待解决的问题是：增强热防护性能与减少人体新陈代谢热负荷总是相互矛盾的 [1] 。基于高温工作防护服装的重要性以及其研究现状，本文通过建立分层热传导方程以及求解其温度–空间–时间分布来分析高温环境下热防护服的不同材质层中的温度变化规律，以辅助于高温工作防护服装的设计问题。模型测试结果显示，与外界相接触的材质层的降温幅度最小，与皮肤相接触的空气层的降温幅度最大。

2. 分层热传导模型

为简化热传递的过程，本文考虑热传递方向是一维的，即只在单位面积皮肤的内法线方向进行热量传递，且在防护服的四层结构中都只存在热传导，而不考虑热辐射和热对流的效应。假定在防护服四层材料之间的连接是无缝隙的，即上一层的下边界与下一层的上边界温度分布相同(从上至下方向为从I层至IV层的方向)。

本文建立的“高温作业防护服–空气–皮肤”的导热模型 [1] 的示意图如图1所示：

Figure 1. Schematic diagram of thermal conduction model

图1. 热传导模型示意图

如图上箭头所示，其方向为热量传递方向。设 $x$ 为层内导热方向上每一个位置点与靠近热源的边界之间的距离。在初始状态即 $t = 0$ 时，认为防护服四层每一个位置节点的温度，都经过假人恒定的体温传导为37℃，外界环境温度固定为一恒定值 $T (0, t) = T_{0}$ ℃。设温度为 $T (℃)$ ，时间为 $t (s)$ ；引入热扩散率 $a = λ / (c ρ)$ 。建立每层热传导二阶偏微分方程如下：

(1)

其中 $x_{j 0}$ ， $t_{j 0}$ 表示第 $j$ 层中的初始空间位置和时间， $x_{j - 1 e n d}$ ， $t_{j - 1 e n d}$ 表示第 $j - 1$ 层中的终点空间位置和时间。

3. 模型的求解方法

在偏微分方程组的求解过程中，每一个偏微分方程的边界条件是一个很重要的已知条件。若已知不同材质层的边界温度，则该模型可以直接通过偏微分的数值解法直接求得其温度-空间-时间分布。在一般环境中，由于无法通过测量得知每一层的边界温度值，仅可测量出外界环境温度和人体皮肤表层温度，因此，本文分别通过引入热阻的概念来确定不同材质层的边界温度以运用Crank-Nicolson隐式格式求解偏微分方程以及将服装整体化运用古典显式差分法迭代进行求解偏微分方程数值解的方法，来分析高温环境下热防护服的不同材质层中的温度变化规律。

3.1. 基于热阻法求解边界条件以及Crank-Nicolson隐式格式求解偏微分方程

本文考虑类似于电阻在所连支路中限流的作用，引入热阻的概念，将防护服的每一层材料视作具有不同的热量阻碍能力的热阻 [2] 。导热热阻的计算公式为 $R = \frac{L}{K S}$ ，其中， $L$ 为某层的厚度， $S$ 为垂直于热传递方向的横截面积， $K$ 为材料的热导率。由相关文献得到：热阻与温度差成正比： $R = \frac{Δ T}{P}$ ，由不同层之间的热阻比例之间的关系得到热传递过程中每一层的温度差，有： $R_{i} : R_{j} \propto Δ T_{i j}, 1 \leq i < j \leq 4$ 。

假定外界环境温度固定为一恒定值，经过一段时间后，假人皮肤外侧温度达到热平衡，可求得I层与IV层的温度差，通过热阻计算方法中的比例关系即可求解出不同材质层中的温度改变量，由此可得到每一层的边界条件 [3] 。

利用Sigmoid函数 $f (x) = \frac{c}{a + e^{- x / b}}$ 对不同材质边界处的数据进行函数拟合 [4] ，将所得到的函数作为边界条件，再运用Crank-Nicolson隐式格式差分法 [5] 来求解分层热传导偏微分方程的数值解。

在利用上述差分法求出微分方程解的近似值之前，先对网格的建立进行说明。如图2所示，横向 $x$ 轴和纵向时间 $t$ 轴的步长分别为 $h$ 和 $k$ ，定义网格比 $r = k / h^{2}$ ；绿色边和两条橙色边分别为一个初始条件和两个边界条件。记蓝色边框内的待求解节点矩阵为

$T_{i} (x, t) = (T_{i} (1, 1), \dots, T_{i} (A x - 2, 1), T_{i} (1, 2), \dots T_{i} (A x - 2, 2), \dots T_{i} (1, B t - 1), \dots, T_{i} (A x - 2, B t - 1))$ (2)

其中表示时刻，与靠热源较近的边界距离为 $x$ 的网格节点的值。

Figure 2. Implicit schematic diagram of Crank-Nicolson

图2. Crank-Nicolson隐式格式示意图

由Crank-Nicolson隐式格式得到的差分形式 [5] 如下：

$\begin{array}{l} - \frac{1}{2} r T_{i} (x - 1, t + 1) + (1 + r) T_{i} (x, t + 1) - \frac{1}{2} r T_{i} (x + 1, t + 1) \\ = \frac{1}{2} r T_{i} (x - 1, t) + (1 - r) T_{i} (x, t) + \frac{1}{2} r T_{i} (x + 1, t) \end{array}$ (3)

由上述格式和边界条件，构造反映时间层 $t$ 与时间层 $t - 1$ 之间的关系方程组如下：

$\begin{array}{l} (\begin{matrix} 1 + 2 r & - r & \dots & 0 \\ - r & 1 + 2 r & ⋱ & ⋮ \\ 0 & ⋱ & ⋱ & - r \\ 0 & 0 & - r & 1 + 2 r \end{matrix}) \cdot (\begin{matrix} T_{i} (1, t + 1) \\ T_{i} (2, t + 1) \\ ⋮ \\ T_{i} (A x - 2, t + 1) \end{matrix}) \\ = (\begin{matrix} 1 + 2 r & - r & \dots & 0 \\ - r & 1 + 2 r & ⋱ & ⋮ \\ 0 & ⋱ & ⋱ & - r \\ 0 & 0 & - r & 1 + 2 r \end{matrix}) \cdot (\begin{matrix} T_{i} (1, t) \\ T_{i} (2, t) \\ ⋮ \\ T_{i} (A x - 2, t) \end{matrix}) \\ + (\begin{matrix} 0.5 r T_{i} (0, t) + 0.5 r T_{i} (0, t + 1) \\ 0 \\ ⋮ \\ 0.5 r T_{i} (A x - 2, t) + 0.5 r T_{i} (A x - 2, t + 1) \end{matrix}) \end{array}$ (4)

因此，利用第0时间层（初始条件）解出第1时间层的值，由第1时间层解出第2时间层直至第 $B t - 1$ 时间层，这种求解方法称为追赶法。由此可得温度-时间-空间的分布结果，利用2018年全国大学生数学建模竞赛A题 [6] 中的相关数据带入该方法中求解得到不同材质层的温度分布图如图3所示。

Figure 3. Temperature-time-space Distribution of Layer I-IV Based on the Implicit Method of Solving Partial Differential Equation of CN

图3. 基于CN隐式格式偏微分方程求解法的第I-第IV层温度-时间-空间的分布

从图3、图4可以看出，沿着热传递方向温度不断下降，表明四个不同的材质层对热量的传递都有一定的阻碍作用，且经过一段时间后各自都达到热平衡，即温度不再随时间变化。此外，可看出IV层外壁的温度下降幅度较大，与用热阻判断温度差幅度的结论相符合。

Figure 4. The change of temperature at the last moment (5400 s) with the position based on the method of solving partial differential equation with CN implicit scheme

图4. 基于CN隐式格式偏微分方程求解法的最后时刻(第5400 s)温度随位置的变化

3.2. 基于整体化法以及古典显式差分法求解偏微分方程

首先考虑将四种不同材质层视为一个整体，在已测量出外界环境温度以及人体皮肤表层温度随时间变化的数据的条件下，分别对四种不同材质层中热传导方程运用古典显式差分法 [7] 迭代进行求解偏微分方程数值解。

因古典显格式稳定收敛的充分必要条件是网格比 $r < 0.5$ ，若不满足此条件，则可采取扩大步长、扩充矩阵的方法来提高模型求解的稳定性。运用古典显格式进行差分，在节点 $T_{i} (x, t)$ 上，用向前差商法得到：

$\frac{T_{i} (x, t + 1) - T_{i} (x, t)}{k} - \frac{T_{i} (x + 1, t) - 2 T_{i} (x, t) + T_{i} (x - 1, t)}{h^{2}} = 0$

整理可得：

$T_{i} (x, t + 1) = r T_{i} (x - 1, t) + (1 - 2 r) T_{i} (x, t) + r T_{i} (x + 1, t)$ (5)

由上式可以看出，第 $t + 1$ 时间层的节点可以由第时间层的节点直接推出来，且第0时间层的值是所已知的初始条件，故可以利用这个格式层层递推求解。考虑将防护服内的温度分布作为一个整体，假设四层材料直接连接是无缝隙的，即相邻两层的内外侧温度分布相同，根据这个条件采用迭代的方法，将第 $i$ 层迭代的结果作为第 $i + 1$ 层的初始条件，求解得到每一层的温度分布，由 [6] 提供的数据求解得到其整体温度分布结果图如图5。

Figure 5. Temperature-time-space distribution of layer I-IV based on classical explicit partial differential equation

图5. 基于古典显式偏微分方程求解法的第I-第IV层温度随时间和空间位置的分布

由图5可以看出，在III层位置呈现“下凹”形态，是由于服装材料在恒温环境下不断吸收热量，服装温度随之升高，从每个升温的瞬间来看，服装相对于身体来说相当于一个热源，因此服装的热量在材质中间某一段会比皮肤表层的温度略小。

为了进一步研究边界处温度变化规律，其I-IV层边界温度随时间变化的数据进行曲线拟合，其结果图如图6所示。

Figure 6. On the basis of the classical explicit partial differential equation solution method, the temperature change of the boundary layer I-IV with time

图6. 基于古典显式偏微分方程求解法的第I-第IV层边界温度随时间的变化

从图6可以看出，四层边界处所求解得结果的拟合效果均较好。其中第II层外壁以及皮肤表面温度经过一段时间后达到稳态，第III层外壁以及第IV层外壁表面温度随时间推移而呈几近线性增长。

4. 模型检验

本文对基于CN隐式格式偏微分方程求解模型进行检验 [8] [9] ，由 [6] 提供的数据，将皮肤外侧温度的测量值进行曲线拟合，比较同一时刻下的计算值与拟合曲线上温度取值之间的差值，误差图如图7。

Figure 7. The calculation and measurement of skin temperature and its error based on the method of solving implicit partial differential equation of CN

图7. 基于CN隐式格式偏微分方程求解法的皮肤温度计算值与测量值及其误差

从误差曲线来看，在实验时间的前段，皮肤外侧的计算值与实际测量值的误差较大，但也都控制在0.04℃以内，说明模型的拟合效果较好。当时间达到大约四百秒之后，误差减小，且在达到热平衡的时间点后，两曲线几乎完全吻合。另一方面，本文假设热传递方向为一维且人体表皮温度为一恒定值，若将该模型扩展至多维方向上且考虑人体表皮温度有变化，只需将人体不同部位的温度数据代入模型中即可 [10] 。若在模型中考虑热辐射和热对流的效应，其结果将会更加精准。

5. 结语

本文建立分层热传导模型来分析热防护服的热量变化规律。利用热阻的概念和整体化两种方法分别拟定偏微分方程组的边界条件，运用Crank-Nicolson隐式格式以及古典显式差分格式分别对该模型进行求解。由于古典显格式稳定收敛的充分必要条件是网格比 $r < 0.5$ ，若不满足此条件，需采取扩大步长、扩充矩阵的方法来提高模型求解的稳定性。因此，Crank-Nicolson隐式格式求解法更加简便，且在模型检验中发现，隐式格式求解的数据在达到热平衡后几乎与实际测量值吻合。

基金项目

本文由湖南农业大学大学生创新性实验计划项目(XCX18007)资助。

文章引用

李碧璇,朱可馨,陈宣霖,王访. 基于分层差分方程的热防护服热量分布规律的模拟
Simulation of Heat Distribution Law of Thermal Protective Clothing Based on Lay-ered Difference Equation[J]. 建模与仿真, 2019, 08(01): 6-13. https://doi.org/10.12677/MOS.2019.81002

参考文献

1. Liu, J. (2002) Numerical Solution of Forward and Backward Problem for 2-D Heat Conduction Equation. Journal of Computational & Applied Mathematics, 145, 459-482. https://doi.org/10.1016/S0377-0427(01)00595-7

2. 张海峰, 葛新石, 叶宏. 预测复合材料导热系数的热阻网络法[J]. 功能材料, 2005, 36(5): 757-759.

3. 熊平, 艾红雷, 卢涛, 等. 一维非稳态导热反问题反演管道内壁面温度波动[J]. 核动力工程, 2018, 39(2): 96-100.

4. 赵兰萍, 徐烈. 固体界面间接触导热的机理和应用研究[J]. 低温工程, 2000(4): 29-29.

5. 孙鸿烈. 解一维热传导方程的高精度的隐式差分格式[J]. 数学的实践与认识, 1998(3): 197-200.

6. http://www.mcm.edu.cn/html_cn/node/7cec7725b9a0ea07b4dfd175e8042c33.html

7. 金承日. 解抛物型方程的高精度显式差分格式[J]. 计算数学, 1991, 13(1): 38-44.

8. 张小峰, 张红武. Crank-Nicolson格式精度的改进[J]. 水科学进展, 2001, 12(1): 33-38.

9. 马继涌, 谢鸣. 多层壁导热反问题的解析公式及其应用[J]. 哈尔滨建筑大学学报, 1995(6): 71-75.

10. 潘斌. 热防护服装热传递数学建模及参数决定反问题[D]: [硕士学位论文]. 杭州: 浙江理工大学, 2017.

NOTES

^*通讯作者。

期刊菜单