利用热传导方程研究高温作业专用服设计的数学模型

Study on Mathematical Model of High Temperature Apparel Design by Using Heat Conduction Equation

摘要: 本文利用热力学知识建立数学模型，研究在实验室环境下，三层隔热材料组成的专用作业服装的隔热效果，计算假人皮肤外侧的温度随时间变化情况。首先根据Fourier试验定律，建立关于温度函数u(x,t)的一维分段热传导方程，使用有限差分法化简该热传导方程，得到有限差分方程组，近似计算原热传导方程的数值解。其次，当第II层厚度不确定时，结合建立的热传导方程，使用二分法计算出的最优厚度。最后进行灵敏度分析对模型的稳定性进行验证，同时确定出当第II和第IV层厚度都不确定时两层的最优厚度。

Abstract: In this paper, we use the thermodynamics knowledge to establish the mathematical model in order to study the thermal insulation effect of the special operation clothing composed of three layers’ insulation materials in the laboratory environment, and calculate the temperature change outside the mannequins’ skin with time. Firstly, according to the Fourier test law, a one-dimensional segmental heat conduction equation for the temperature function u(x,t) is established. Simpli-fied the heat conduction equation by the finite difference method, the finite difference equations are obtained to approximate the numerical solution of the original heat conduction equation. Secondly, when the thickness of the second layer is uncertain, combined with the established heat conduction equation, the optimal thickness is calculated by the dichotomy method. Finally, the sensitivity analysis is performed to verify the stability of the model, and the optimal thickness of the two layers is determined when the thicknesses of the second and fourth layers are both uncertain.

1. 引言

1.1. 问题重述

高温防护服通常由三层物料组成，第I层直接与外界环境接触，第II层为防止人体免受高温侵害起到重要作用，第III层靠近人体皮肤表面，但不直接接触，其与皮肤间的空气层为第IV层 [1] [2] 。现解决以下问题

1) 首先建立热传导方程模型，根据附件，需要计算温度分布，进而生成温度分布的Excel文件。

2) 当第II层厚度不确定时，依据边值条件、约束条件和建立的热传导方程模型等其他因素，明确第II层厚度满足什么条件才是最优的，进而求解出第II层的最优厚度。

3) 当第II和第IV层厚度都不确定时，需要探究两个分层之间的杠杆关系，明确第II和第VI层的厚度满足什么条件时才是最优的，然后由建立的热传导方程模型，依据边值条件和约束条件，求出最优的第II和第IV层的厚度。

1.2. 问题分析

针对问题一：因为热量传导方向共线，而且外界温度较低、第IV层厚度较小，所以忽略热辐射和热对流的影响 [3] [4] 。依据Fourier实验定律 [5] [6] [7] ，建立一维分段热传导方程。然后根据附件2中的数据求解热传导方程模型中未知的参数。原热传导方程组是无法求出解析解的，所以考虑有限差分法或者有限元法进行求解。最终依据边界条件，计算出温度分布，生成温度分布的Excel文件。

针对问题二：考虑到成本等因素，第II层的最优厚度就是在满足约束条件的情况下，第II层能达到的最小厚度。然后结合边界条件、约束条件和一维分段热传导方程，使用二分法求解出第II层的最优厚度。

针对问题三：结合成本、人体工程学等因素，推导出第II层和第IV层之间的联系，使用逻辑推导来简化问题，并最终可以转化为问题二，然后改变第IV层的厚度。

2. 模型假设

1) 只考虑热传导效应，不考虑热对流、热辐射效应。热力传递的基本方式有：热对流，热传导，热辐射三种，而生产生活中所遇到的热传递现象常是这三种基本方式的不同主次组合。因为外界环境较低，所以忽略热辐射效应。热对流只发生在液体和流体中，而第IV层的最大厚度为6.4 mm，所以可以忽略热对流效应。

2) 实验开始前，外界环境温度、人体各处温度均为37℃，且在实验开始后瞬间外界环境温度发生突变。

3. 模型的建立与求解

3.1. 一维分段热传导方程的建立

设 $u\left(x,t\right)$ 温度关于位置x和时间t的二元函数。推导一维分段热传导方程需要依据Fourier实验定律：

在 $\Delta t$ 时间内，通过面积元s流入小体积元的热量Q与沿面积元外法线方向的温度变化率 $\frac{\partial u}{\partial n}$ 成正比。

因为热量传导的方向共线，所以在本文中只考虑一维情况，所以在 $\Delta t$ 时间内，通过面积元s流入小

体积元的热量Q与沿长度微元法线方向 $\frac{\partial u}{\partial x}$ 成正比。

现考虑长度微元 $\Delta x$ ，如果温度 $u\left(x+\Delta x,t\right)>u\left(x,t\right)$ ，那么 $\Delta Q$ 将是负值。规定热量从x处流向 $x+\Delta x$ 处的方向为正向，如果 $x+\Delta x$ 处的温度高于x处的温度，则热量必然会从 $x+\Delta x$ 处流向x处，让 $\Delta x\to \text{0}$ ，

差商项 $\frac{u\left(x+\Delta x\right)-u\left(x\right)}{\Delta x}$ 将接近 $\frac{\partial u}{\partial x}$ 。

定义 $x=0$ 处为第I层与外界接触的位置。

长度微元吸收的热量：

$Q_{吸}=-k\frac{\partial u}{\partial x}\cdot s\cdot \Delta t$ (1)

长度微元放出的热量：

$Q_{放}=-k{\frac{\partial u}{\partial x}|}_{x+\Delta x}\cdot s\cdot \Delta t$ (2)

由(1)、(2)得出长度微元在 $\Delta t$ 时间内获得的热量为：

$Q_{吸}-Q_{放}=k\left({\frac{\partial u}{\partial x}|}_{x+\Delta x}-\frac{\partial u}{\partial x}\right)\cdot s\cdot \Delta t$ (3)

又由热量计算公式 $Q=c\cdot m\cdot \Delta t$ ，得：

$Q_{吸}-Q_{放}=c\cdot \Delta m\left(u\left(x,t+\Delta t\right)-u\left(x,t\right)\right)=k\frac{{\partial }^{2}u}{\partial x^2}\cdot \Delta x\cdot s$ (4)

其中： $\Delta m=\rho \cdot s\cdot \Delta x$

联立(3)、(4)解得

$\frac{\partial u}{\partial t}={\alpha }_i^2\frac{{\partial }^{2}u}{\partial x^2}\left(i=1,2,3,4\right)$ (5)

其中

${\alpha }_i=\{\begin{array}{l}\sqrt{k_1/\left({\rho }_1{c}_1\right)}i\in \left[0,n_1\right]\\ \sqrt{k_2/\left({\rho }_2{c}_2\right)}i\in \left[n_1,n_1+n_2\right]\\ \sqrt{k_3/\left({\rho }_3{c}_3\right)}i\in \left[n_1+n_2,n_1+n_2+n_3\right]\\ \sqrt{k_4/\left({\rho }_4{c}_4\right)}i\in \left[n_1+n_2+n_3,n\right]\end{array}$ (6)

由题目信息可以得出条件：

据Fourier实验定律得：

${\frac{\text{d}u}{\text{d}x}|}_{x=d}=-\lambda \left(u\left(d,t\right)-37\right)$ (7)

$\lambda$ 为第IV层向第V层(皮肤外侧到体内温度恒定为37℃的区域)传导热量的比例系数。

$t=0$ 时，第I层到第V层的温度都为37℃：

$u\left(x,0\right)=37,0\le x\le d$ (8)

$t>0$ 时，外界温度突变为75℃：

$u\left(0,t\right)=75,t>0$ (9)

综上所述，我们建立的一维分段热传导方程为：

$\{\begin{array}{l}\frac{\partial u}{\partial t}={\alpha }_i^2\frac{{\partial }^{2}u}{\partial x^2}，i=1,2,3,4\\ u\left(x,0\right)=37，0\le x\le d\\ u\left(0,t\right)=75，t\ge 0\\ {\frac{\text{d}u}{\text{d}x}|}_{x=d}=-\lambda \left(u\left(d,t\right)-37\right)\end{array}$ (10)

其中 ${\alpha }_i$ 如上文所述。

3.2. 有限差分法简化热传导方程

求解一维热传导方程的方法主要为有限差分法和有限元法 [8] [9] 。有限差分法的基本理论是使用有限个网格节点构成的细小网格来替代连续的定解区域，将热传导方程和定解条件使用差商来近似，于是方程和定解条件就可以使用有限差分方程组进行近似计算。有限差分法理论成熟、可以选择精度，而且易于编程，所以选择有限差分法来求解上文建立的数学模型。

首先将x离散化，取步长为 $\Delta x={10}^{-4}\text{m}$ 。

第I层被分为 $n_1$ 个小格 $n_1=\frac{x_1}{\Delta x}$ ；第II层被分为 $n_2$ 个小格 $n_2=\frac{x_2}{\Delta x}$ 。

第III层被分为 $n_3$ 个小格 $n_3=\frac{x_{\text{3}}}{\Delta x}$ ；第IV层被分为 $n_4$ 个小格 $n_4=\frac{x_{\text{4}}}{\Delta x}$ 。

定义辅助函数 $u_i\left(t\right)=u\left(i\left(\Delta x\right),t\right)$ ，i为位置x被划分的第i个小格。

由Lagrange中值定理得出当满足 $2\le i\le n-1$ 的微分方程：

$u_{i+1}\left(t\right)-u_i\left(t\right)=u\left(\left(i+1\right)\left(\Delta x\right),t\right)-u\left(i\left(\Delta x\right),t\right)\approx {\frac{\partial u}{\partial x}|}_{\left(i+1/2\right)\Delta x}\left(t\right)\Delta x$ (11)

$u_i\left(t\right)-u_{i-1}\left(t\right)=u\left(i\left(\Delta x\right),t\right)-u\left(\left(i-1\right)\left(\Delta x\right),t\right)\approx {\frac{\partial u}{\partial x}|}_{\left(i-1/2\right)\Delta x}\left(t\right)\Delta x$ (12)

$u_{i+1}\left(t\right)+u_{i-1}\left(t\right)-2u_i\left(t\right)\approx {\frac{{\partial }^{2}u}{\partial x^2}|}_{i\Delta x}\cdot {\left(\Delta x\right)}^2$ (13)

原方程为

$\frac{\partial u}{\partial t}={\alpha }^2\frac{{\partial }^{2}u}{\partial x^2}$ (14)

将原方程离散化，x替换为 $i\Delta x$ 得出：

$\frac{\partial u_i\left(t\right)}{\partial t}={\alpha }_i^2\frac{u_{i+1}\left(t\right)+u_{i-1}\left(t\right)-2u_i\left(t\right)}{{\left(\Delta x\right)}^2}$ (15)

又因为

$\frac{\partial u_i\left(t\right)}{\partial t}=\frac{u_i\left(t_{j+1}\right)-u_i\left(t_j\right)}{\Delta t}=\frac{{\alpha }_i^2}{{\left(\Delta x\right)}^2}\left[u_{i+1}\left(t_j\right)+u_{i-1}\left(t_j\right)-2u_i\left(t_j\right)\right]$ (16)

结合边值条件最终得到：

我们分三类进行讨论

i) 第1个小格满足的方程为：

$u_1\left(t_1\right)=u_1\left(t_0\right)+\frac{{\alpha }_1^2}{{\left(\Delta x\right)}^2}\left(u_2\left(t_0\right)+75-2u_1\left(t_0\right)\right)\Delta t$ (17)

其中： $u_1\left(t_0\right)=37$ ， $u_2\left(t_0\right)$ 为第I层第2个小格在 $t=0$ 的温度

ii) 第2到第 $n-1$ 个小格满足的方程为：

$u_i\left(t_{j+1}\right)=u_i\left(t_j\right)+\frac{{\alpha }_i^2}{{\left(\Delta x\right)}^2}\left[u_{i+1}\left(t_j\right)+u_{i-1}\left(t_j\right)-2u_i\left(t_j\right)\right]\cdot \Delta t$ (18)

iii) 第n个小格(也即最后一个小格)满足的方程为

$u_n\left(t_{j+1}\right)=u_n\left(t_j\right)+\frac{{\alpha }_4^2}{{\left(\Delta x\right)}^2}\left(u_{n-1}\left(t_j\right)-u_n\left(t_j\right)-\beta \left(u_n\left(t_j\right)-37\right)\right)\Delta t$ (19)

由题意给出如下定解条件：

${\frac{\text{d}u}{\text{d}x}|}_{x=d}=-\lambda \left(u\left(d,t\right)-37\right)$ (20)

$\lambda$ 为第IV层向第V层传导热量的比例系数

$t=0$ 时，第I层到第V层的温度恒为37℃

$u\left(x,0\right)=37,0\le x\le d$ (21)

环境温度为75℃

$u\left(0,t\right)=75,t>0$ (22)

通过上述递推关系，可建立如下映射关系，进而求得方程的解：

$\left(\begin{array}{c}37\\ 37\\ ⋮\\ 37\end{array}\right)\stackrel{u_i\left(t\right)}{\to }\left(\begin{array}{c}{u}_1\left(t_1\right)\\ u_2\left(t_1\right)\\ ⋮\\ u_d\left(t_1\right)\end{array}\right)\stackrel{u_i\left(t\right)}{\to }\cdots \stackrel{u_i\left(t\right)}{\to }\left(\begin{array}{c}{u}_1\left(t_n\right)\\ u_2\left(t_n\right)\\ ⋮\\ u_d\left(t_n\right)\end{array}\right)$ (23)

综上所述，得出用有限差分法简化的热传导方程

$\{\begin{array}{l}{u}_1\left(t_1\right)=u_1\left(t_0\right)+\frac{{\alpha }_1^2}{{\left(\Delta x\right)}^2}\left(u_2\left(t_0\right)-2u_1\left(t_0\right)+75\right)\Delta t\\ u_i\left(t_{j+1}\right)=u_i\left(t_j\right)+\frac{{\alpha }_i{}^2}{{\left(\Delta x\right)}^2}\left[u_{i+1}\left(t_j\right)+u_{i-1}\left(t_j\right)-2u_i\left(t_j\right)\right]\cdot \Delta t\\ u_n\left(t_{j+1}\right)=u_n\left(t_j\right)+\frac{{\alpha }_4^2}{{\left(\Delta x\right)}^2}\left(u_{n-1}\left(t_j\right)-u_n\left(t_j\right)-\beta \left(u_n\left(t_j\right)-37\right)\right)\Delta t\end{array}$ (24)

其对应的边界条件为

$\{\begin{array}{l}u\left(x,0\right)=37,0\le x\le d\\ u\left(0,t\right)=75,t>0\\ {\frac{\text{d}u}{\text{d}x}|}_{x=d}=-\lambda \left(u\left(d,t\right)-37\right)\end{array}$ (25)

3.3. 问题一的求解

使用(17)、(18)、(19)式和(20)、(21)、(22)初值条件，取定不同的 $\beta$ 值，得出第n层0-5400，步长为1( $\Delta t=1\text{s}$ )秒时各时刻所对应的温度值，这一温度值是理论所得的温度值；结合附件中实际测量的温度值，计算二者之间距离和的最小值来确定最佳的 $\beta$ ，距离最小值往往很难确定，根据最小二乘原理，将距离转化为二者差的平方，二者差的平方和最小时对应的 $\beta$ 值即为最优解，记为 ${\beta }_0$ 。

首先我们给定几个 $\beta$ 值，观察曲线的大致走向，大致找到误差平方和最小时 $\beta$ 的取值区间，再在此区间中进一步搜索，找到当 $\beta =0.01645$ 时，误差平方和最小。如图1所示，将此时的 $\beta$ 值记为 ${\beta }_0$ ，此后的计算 $\beta$ 均取为 ${\beta }_0$ .使用有限差分方程组(13)近似求出温度随位置、时间的分布。

Figure 1. The relationship between the sum squares of errors and β

图1. 误差平方和与β的关系

由图1求得当 $\beta =0.01645$ 时，误差平方和最小，将此时的 $\beta$ 值记为 ${\beta }_0$ ，此后的计算 $\beta$ 均取为 ${\beta }_0$ 。使用有限差分方程组(24)近似求出温度随位置、时间的分布，将计算生成的数据保存于文件名“problem1.xlsx”的Excel文件。

3.4. 问题二的求解

当环境温度为65℃，IV层厚度为5.5 mm时，为确定第II层的最优厚度，忽略高温作业专用服的物理性质，仅考虑成本，最优的厚度就是在达到高温防护的目的后，尽可能小的厚度。因而求第II层的最优厚度问题实质上为求一个尽可能小的厚度值。

在(1)所建立模型的基础上，对方程的边值条件进行修改，此处注意到，温度随时间的变化关系是一个增函数，工作60分钟时，假人皮肤外侧温度不超过47℃，只需令第n个小格在60分钟时温度小于47℃即可。此时，满足 $u_n\left(t_1\right)\le T_1$ ，其中： $t_1=3600\text{s}$ ， $T_1=47˚\text{C}$ ，为确保超过44℃的时间不超过5分钟，需满足： $u_n\left(t_2\right)\le T_2$ ，其中： $t_2=3300\text{s}$ ， $T_2=44˚\text{C}$ ，并且满足如下初始条件：

$\{\begin{array}{l}u\left(0,t\right)=65\\ u\left(x,0\right)=37\end{array}$

使用二分法求解得出第I层的最优厚度为：8.50 mm。

3.5. 问题三的求解

需要改变第II层和第IV层的厚度，满足一系列条件，进而确定第II层和第IV层的最优厚度。

必须满足的条件有：

1、在环境温度为80℃和工作为30分钟时，一直不能超过47℃。即第n个小分块在30分钟时温度不能超过47℃， $u_n\left(t_1\right)\le T_1$ ，其中 $t_1=1800\text{s}$ ， $T_1=47˚\text{C}$ .

2、在环境温度为80℃和工作为30分钟时，假人的皮肤外侧温度超过44℃的时间不能超过5分钟。即： $u_n\left(t_2\right)\le T_2$ ，其中 $t_2=1500\text{s}$ ， $T_2=44˚\text{C}$ .

在满足上述两个条件后，尽可能需要满足的条件有：

3、高温作业专用服的成本最低。因为第IV层是空气，所以成本只与第II层有关。在满足第一和第二个条件时，第II层越薄，高温作业专用服的成本越低。

4、符合人体工程学，既高温作业专用服重量最轻,要使得高温作业专用服尽可能轻，因为第IV层是空气，密度相对于第II层很低，所以高温作业的重量只与第II层有关。在满足第一和第二个条件时，第II层越薄，高温作业专用服的重量越小，所以不妨让第IV层的厚度取最大，即取 $x_4=6.4\text{mm}$ ，这样可以最大限度减小第II层的厚度。

综上所述，确定第II层和第 IV 层的最优厚度就转化为：在满足第一和第二个条件的前提下，设定第III层厚度取最大，然后求解第II层的最小厚度，这就将问题三转化为问题二。

并且将最优的厚度仍需满足下面的定解条件：

$\{\begin{array}{l}u\left(0,t\right)=80\\ u\left(x,0\right)=37\\ u_n\left(t_1\right)\le T_1\\ u_n\left(t_2\right)\le T_2\end{array}$ (26)

其中 $t_1=1800\text{s}$ ， $T_1=47˚\text{C}$ ，其中 $t_2=1500\text{s}$ ， $T_2=44˚\text{C}$ 。

类似于问题二中求解第II层最小厚度的问题，当第IV层的厚度取6.4 mm时，求得第II层的厚度为18.50 mm。

3.6. 灵敏度分析

此处的灵敏度分析用于对上述模型的稳定性进行检验，上文在第IV层厚度取最大值的条件下，求出了第II层的最小厚度，现在考虑逐渐减小第IV层厚度，计算第II层厚度的变化幅度，直到不满足约束条件，则取上一次计算的第II和第IV层厚度。

Figure 2. Contrast chart of temperature change (at this time, the thickness of Layer II changed and Layer IV remain unchanged)

图2. 温度变化对比图(此时第II层的厚度变化，第IV层厚度不变)

由图2知第II层关于厚度的变化非常敏感，所以第II层的最优厚度为18.50 mm，同样的第IV层关于厚度的变化非常敏感，所以第IV层的最优厚度为6.40 mm。

上文建立的模型非常敏感，主要原因是小分格取得不够细，但是将小分格进一步取细，计算量就会特别大，所以下文进行模型的改进。

4. 模型的改进

在上文中，将一维分段热传导方程离散两次，原偏微分方程使用有限差分方程组进行近似求救，虽然有限差分法的精度可以自己选择，将网格分的越细精度越高，但是随之而来的问题是程序的运行时间大大增加，并且使用近似的思想本身就存在一定的误差。因此考虑将一维分段热传导方程仅离散一次，原微分化为常微分方程。

$\{\begin{array}{l}{\frac{\text{d}u}{\text{d}x}|}_{x=d}=-\lambda \left(u\left(d,t\right)-37\right)\\ \frac{\text{d}{u}_1\left(t\right)}{\text{d}t}=\frac{{\alpha }_i^2}{{\left(\Delta x\right)}^2}\left(75+u_2\left(t\right)-2u_1\left(t\right)\right)\\ \frac{\text{d}{u}_i\left(t\right)}{\text{d}t}=\frac{{\alpha }_i^2}{{\left(\Delta x\right)}^2}\left(u_{i+1}\left(t\right)+u_{i-1}\left(t\right)-2u_i\left(t\right)\right)\\ \frac{\text{d}{u}_{n_1+n_2+n_3+n_4}\left(t\right)}{\text{d}t}=\frac{{\alpha }_4^2}{{\left(\Delta x\right)}^2}\left(u_{n-1}\left(t\right)-u_n\left(t\right)-\mu \left(u_n-37\right)\right)\end{array}$ (27)

致谢

首先感谢天津师范大学数学科学学院和教务处的诸位老师，老师们严谨治学、勤奋严谨，是我们学习的楷模。其次感谢我们的同学对我们在学习和生活上的帮助。最后感谢我们的父母，父母一直是我们坚强的后盾，感谢他们的支持。

基金项目

天津师范大学大学生创新创业训练计划项目资助(项目号：201810065353)。