2018A：高温作业专用服装设计(1)【公众号：数模加油站】.pdf

高温作业专用服装设计 摘要 本文主要研究高温作业专用服装设计，以 Fourier 定律和能量守恒定律为理论依据，建立了基于热传导方程的温度分布模型，借助追赶法求解。在问题一温度分布模型求解中。首先，基于 Fourier 定律和能量守恒定律，建立的基于热传导方程的温度分布模型。基于牛顿冷却定律，借助枚举法，确定空气与皮肤表面的转化系数，并给出初值温度 37 度、左边界 Dirichlet 边值条件，右边界 Robin 边值条件，及基于临界面热流量密度和温度相等的耦合条件。其次，将连续定解区域作网格剖分，用隐式向后差分格式对原微分方程组离散化，得到三对角线性方程组，借助追赶法求解，得到时间与空间维度下的温度分布，见problem1.xlsx。最后，对模型进行误差分析，定义偏差指数 f 并求得其值为 0.4593，最大误差为 1.99。在问题二求 II 介质最优厚度问题中，建立单目标优化模型。首先，基于对服装成本和舒适度的考虑，制定 II 介质厚度的“最优”准则最小厚度为最优，进而确定优化目标；其次，确定约束条件：初值温度 37 度、左边界 Dirichlet 边值条件、右边界 Robin 边值条件、基于临界面热流量密度和温度相等的耦合条件及题目对于温度的限制条件；然后，用循环遍历的枚举法，借助 matlab 搜索出II 介质的最优厚度为 19.3mm。最后，对单目标优化模型作灵敏性分析，最优厚度与温度呈现线性关系。在问题三求 II、IV 介质厚度的问题中，建立多目标优化模型。首先，从成本与穿着舒适度两方面，制定“最优”准则，并确定两个不同的优化目标；然后，借助 matlab 采用双重 for 循环枚举遍历，搜索出 II 介质的最优厚度为 21.7mm，IV 层介质的最优厚度为 6.4mm。关键词：Fourier定律 热传导方程 追赶法 枚举法 向后差分 获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：5444576572 一、问题重述 在高温环境下工作时，人们需要穿着专用服装以避免灼伤。专用服装通常由三层织物材料构成，记为I、II、III层，其中I层与外界环境接触，III层与皮肤之间还存在空隙，将此空隙记为IV层。为设计专用服装，将体内温度控制在37C的假人放置在实验室的高温环境中，测量假人皮肤外侧的温度。为了降低研发成本、缩短研发周期，请你们利用数学模型来确定假人皮肤外侧的温度变化情况，并解决以下问题：(1)专用服装材料的某些参数值由附件1给出，对环境温度为75C、II层厚度为6 mm、IV层厚度为5 mm、工作时间为90分钟的情形开展实验，测量得到假人皮肤外侧的温度（见附件2）。建立数学模型，计算温度分布，并生成温度分布的Excel文件（文件名为problem1.xlsx）。(2)当环境温度为65C、IV层的厚度为5.5 mm时，确定II层的最优厚度，确保工作60分钟时，假人皮肤外侧温度不超过47C，且超过44C的时间不超过5分钟。(3)当环境温度为80时，确定II层和IV层的最优厚度，确保工作30分钟时，假人皮肤外侧温度不超过47C，且超过44C的时间不超过5分钟。二、问题分析 2.1 问题一分析 问题一，建立基于热传导方程的温度分布模型，确定在一维空间中介质在不同时刻，不同厚度下的温度。在模型建立时本文首先借助导热基本定律傅里叶定律和能量守恒定律推导热传导方程。其次简化问题，将四个介质层视为两个新的介质层介质 a、介质 b，在只有 a、b 介质层的情况下，基于临界处温度与热流量密度相同进行二层耦合。最后，将二层耦合推广到四层耦合，建立基于热传导方程的的温度分布模型。模型的求解采用隐式向后有限差分近似对方程进行离散化处理，给出方程的差分格式并整理得到代数方程组，采用追赶法求解方程组，得到时间与空间维度下的温度分布情况。2.2 问题二分析 问题二，求解 II 层介质最优厚度是一个最优化问题，首先从服装成本与穿着舒适度两个方面讨论“最优”标准的制定，确定优化问题的目标为 II 层介质厚度最小。其次，考虑问题二关于假人皮肤外侧温度的两个要求，同时结合问题一建立的基于热传导方程的温度分布模型，确定最优化问题的约束条件，从而建立 II层最优厚度的单目标优化模型。问题二模型的求解利用循环遍历的变步长枚举法，对 II 层介质的所有可能厚度进行遍历，求出满足约束条件的最小厚度。2.3 问题三分析 问题三，求解 II，IV 两层的最优厚度是一个多目标的优化问题。首先，从服装成本与穿着舒适度两个方面考虑，分别制定出不同方面下的“最优”厚度标准，确定多目标优化问题的两个不同的目标。获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：5444576573 其次，基于问题一建立的基于热传导方程的温度分布模型，考虑问题三提出的两个要求，给出最优化问题的约束条件，分别建立目标是服装成本最低和穿着舒适度最高的两个 II，IV 层厚度优化模型。问题三模型的求解采用双重循环遍历的枚举法，借助 matlab 对 II 介质与 IV介质厚度同时进行双重循环遍历，搜寻服装成本最低和穿着舒适度最高这两个目标下 II，IV 层介质的最优厚度。三、模型假设 1假设四层介质均匀，且保持各项同性。2假设每层介质的热传导率在各个方向相同。3假设在第四层介质中，不考虑空气对流。4假设外界无辐射。四、符号说明与名词解释 4.1符号说明 符号 说明 iL 第 i 层介质的厚度，i=1,2,3,4(,)iu x t 第 i 层介质在 t 时刻厚度 x 下的温度，i=1,2,3,4(,)ijux t 第 i 层介质在第 j 时间层中 t 时刻下的厚度 x 时的温度，i=1,2,3,4 q 热流量密度 ic 第 i 层介质的比热，i=1,2,3,4 i 第 i 层介质的密度，i=1,2,3,4 i 第 i 层介质的热传导率，i=1,2,3,4 Q 单位时间通过截面的热量 h 空气与皮肤的转化系数 i 第 i 层与第 i+1 层介质的临界面，i=1,2,3 4.2名词解释（1）比热：是指没有相变化和化学变化时，一定量均相物质温度升高 1K 所需获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：5444576574 的热量，单位为：oJ kgC。（2）热传导率：是物性参数，在数值上等于单位温度梯度作用下单位时间内单位面积的热量，单位为：oW mC。它与物质的结构与状态密切相关，与几何形状无关。（3）冷却系数：是空气与皮肤的转化系数，空气自然对流下大致范围是5,25。五、模型建立与求解 实际问题的数学化处理：在高温作业过程1中，实际为三维立体的空间温度扩散。由于不计空气对流，不考虑人的体积，因此忽略高度与宽度这两个维数，将服装的三层材料与空气层，简化为只与厚度 L 有关的一维空间，如图 1 下：图 1 一维简化图 其中，I、II、III 为服装的三层织物材料，IV 为 III 层与皮肤间的空气层，且I 层与外界环境接触。并分别将 I、II、III、IV 层记为四层介质。5.1 问题一：确定温度分布情况 针对问题一，需要建立数学模型，计算温度分布。由一维简化图可知：只需要确定在一维空间中介质在不同时刻，不同厚度下的温度。首先，借助导热基本定律Fourier 定律和能量守恒定律推导热传导方程。其次，简化问题，将 I、II、III、IV 四个介质层视为两个新的介质层介质 a、介质 b，在只有 a、b 介质层的情况下，进行二层耦合。最后，将二层耦合推广到四层耦合，建立基于热传导方程的温度分布模型。模型的求解采用隐式向后有限差分近似对方程进行离散化处理，给出方程的差分格式并整理得到代数方程组，采用追赶法求解方程组，得到时间与空间维度下的温度分布情况 5.1.1两大定律Fourier定律和牛顿冷却定律 在求解温度分布过程中，由于热量随时间进行扩散，因此本文考虑导热现象和冷却现象。（1）Fourier 定律 热量与热流密度 获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：5444576575 在导热现象中，单位时间内通过截面面积为 S 的截面所传递的热量 Q，正比例于垂直于该截面方向上的温度变化率，但热量传递的方向与温度升高的方向相反，如下图 2 所示：图 2 Fourier 定律示意图 傅里叶定律表达式为：=QuuQSSxx (1)用热流密度表示为：uqx=(2)其中：负号表示热量传递的方向与温度升高的方向相反；x 表示从外界到皮肤的厚度，为空间坐标；u=u(x,t)表示关于厚度 x 和时间 t 的函数；ux 表示温度沿 x 轴方向的变化率；表示热传导率，均匀介质中为一固定数值。一维空间中热流密度矢量 本文考虑在一维空间下温度分布，即只有一个坐标 x（表示厚度），具体见下图 3：图 3 单介质示意图 则热流密度矢量的形式为：获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：5444576576 uqgradunn=(3)其中：gradu 是一维空间某厚度下的温度梯度；n 是在临界面的外法向向量。（2）牛顿冷却定律 牛顿冷却定律是温度高于周围环境的物体向周围媒质传递热量逐渐冷却时所遵循的规律4。具体的表述为：当物体表面与周围存在温差时，单位时间从单位面积散失的热量与温度差成正比，比例系数为热传导率。计算公式为：(37)uh un=空。5.1.2热传导方程的推导（1）基于 Fourier 热传导定律的热量的微元算式 在上述 Fourier 定律推导过程3中，得到热量的微元算式如下：udQdsdtu dSdtn=(4)（2）流入热量与吸收热量 从12tt 时间段内，流入介质内部的热量为：211=ttSuQds dtn (5)并将（5）式借助高斯公式化简为：22111=tttStuuQds dtdx dtnx (6)介质内温度升高吸收的热量为：221=(,)(,)Qcu x tu x tdx (7)其中，c 为比热，为密度。并将（6）式化简为：获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：5444576577 2121221=(,)(,)(,)=(,)=ttttQcu x tu x tdxu x tcdt dxtu x tcdx dtt (8)（3）基于能量守恒定律建立能量关系 由能量守恒定律，有12QQ=，即：2211(,)=ttttuu x tdx dtcdx dtxt (9)因此，有22uuctx=，即222=uuatx，且2=ac。其中，表示为介质的热传导率；c 表示为介质的比热，表示为介质的密度。最终推导的热传导方程为：222=uuatx，2=ac 5.1.3 二层耦合介质温度分布模型的建立 首先，在只有介质层 a、b 的情况下，建立二层耦合的温度分布模型。先确定在不同介质中的热传导方程，再根据在同一临界面具有相同的热流量密度和温度进行二层耦合。其次，将二层耦合扩展为四层耦合，确定完整的温度分布模型。详细图解如下图 4：图 4 两层耦合图 Step1：确定介质：确定介质 a、介质、介质 b 的热传导方程。的热传导方程。获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：5444576578 由 5.1.2 的热传导方程推导过程可知：2=ac，则带入相应的热传导率，比热 c，密度，分别确定介质 a、b 的热传导方程，见表 1。表 1 介质 a、介质 b 的热传导方程 热传导方程 方程参数 空间范围 介质 a 221112uuatx=21111=ac 10,xd 介质 b 222222uuatx=22222=ac 12,xd d Step2：确定耦合条件。：确定耦合条件。由介质临界面的热流量密度相同，确定第一耦合条件。根据 Fourier 热传导定律，在导热过程中，两相邻介质临界面处热流量密度相同，得到：11121212dduunn=(10)其中，1n,2n分别表示介质 a、介质 b 在临界面上的外法线方向；1，2分别表示介质 a、介质 b 上的热传导率。由介质临界面的温度相同，确定第二耦合条件。根据两介质临界面处温度相同，得到任意时刻的临界面处温度等价关系：1112dduu=(11)完整耦合条件的确立。由（9）式和（10）式可得，两层耦合状态下的耦合条件为：121112121212dddduunnuu=(12)Step3：确定两层耦合介质的初边值条件。：确定两层耦合介质的初边值条件。确定初值条件。获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：5444576579 在 t=0 时刻，介质 a、介质 b 的温度均与假人的恒定皮肤外侧温度相同，即均为 37oC，由此确定初值条件为：11212(,0)37 0,(,0)37 ,u xxduxxd d=(13)确定边值条件。第一，左边界条件的确定。介质 a 左侧与外界接触，因而其温度与恒定的外界温度相同，故方程的左边界 Dirichlet 边值条件为：1(0,)75 0,uttT=(14)第二，右边界条件的确定。介质 a 接受了温度为 75oC的外界传递的热量，并将热量继续传递给介质 b，并由附件 2 提供的数据可知：介质 b 右边界的温度始终高于假人皮肤的温度 37oC。进而，介质 b 右边界与假人皮肤发生热量交换，相当于对介质 b 产生冷却作用。此时，假人相当于冷却源，这个热量传递过程遵循牛顿冷却定律。依据牛顿冷却定律，给出 Robin 右边界条件：454()huh uun=(15)其中，hu 表示为介质 b 在贴近假人皮肤处的温度值；5u表示为假人内部恒定温度，即 37oC；h 为介质与皮肤之间的转化系数，即冷却系数。由（13）式和（14）式，确定的热传导方程的完整的边值条件为：1454(0,)75 0,()huttTuh uun=(16)Step4：确定转化系数：确定转化系数 h 的值。的值。Step5：确定二层耦合介质温度分布模型。：确定二层耦合介质温度分布模型。根据以上三个步骤，分别得到了热传导方程以及初值条件、边值条件和耦合获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：54445765710 条件，则在只有两介质的情况下，确立基于热传导方程的二层耦合温度分布模型：1211221111222222122211112222212121212112121 0,=(,0)37 0,(,0)37 ,(0,)dddduuaxdtxuuaxd dtxacacuunnuuu xxduxxd dut=45475 0,()btTuh uun=(17)5.1.4 基于热传导方程的温度分布模型的建立 在二层耦合介质的基础上，用相同的方法，扩展为四层介质。四层耦合介质图如下所示：图 5 四层耦合介质图 Step1：确定四层介质的热传导方程。：确定四层介质的热传导方程。24210212 ,0 1,2,3,4iiiiiiiiiiuuaxLLLtxaic=(18)Step2：确定耦合条件、初值条件及边值条件。：确定耦合条件、初值条件及边值条件。（1）相邻两介质的临界面共有 3 个，即：1，2，3。在每一个临界面获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：54445765711 的热流量密度和温度相同，可得到两个耦合条件。因此在三个临界面上共确定 6个耦合条件，详见温度分布模型；（2）由 t=0 时刻，四层介质的温度均与假人的恒定皮肤外侧温度相同，确定初值条件：(,0)37 1,2,3,4iu xi=(19)（3）由介质 I 左侧与外界温度相同，确定左边界 Dirichlet 边值条件：1(0,)75ut=(20)由介质 IV 右侧与假人进行热量传递，依据牛顿冷却定律，确定右边界 Robin边值条件：445()huh uux=(21)Step3：确定转化系数：确定转化系数 h 的值。的值。首先，查阅文献资料知：转化系数 h 的范围是5,25。其次，通过附件 2 给出的不同时刻下假人皮肤表面的温度值，借助变步长多次枚举法确定最佳转化系数的值。具体过程如下：粗略估计转化系数 h 的值。将 h=8，h=9，与实际数据用 matlab 进行数据处理绘制如下对比图：图 6 在不同转化系数h下的皮肤温度变化图 由上图可知：红色实线为实际数据温度变化情况；h=8 即黄色实线，在实际数据的上方；h=9 即蓝色实线，在实际数据的下方。因此得到最佳吻合的转化系数 h 值介于区间8,9之中。获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：54445765712 设置的步长为 0.01，在8,9范围内，经 matlab 枚举遍历，确定 h 的值在8.62 左右。为达到 0.0001 的精度，进一步缩小 h 的步长为 0.0001，8.1,8.3范围内，再次经 matlab 枚举遍历，确定最佳转化系数 h 的精确值为 8.6125。故，空气与皮肤的转化系数 h=8.6125。Step4：得到最终的基于热传导方程的温度分布模型：得到最终的基于热传导方程的温度分布模型：1111222421021212121232232 ,0 1,2,3,4 iiiiiiiiiiuuaxLLLtxaicuuxxuuuuxxu=22333333434341445 (,0)37 1,2,3,4 (0,)75()ihuuuxxuuu xiutuh uux=(22)说明：222iiiuuatx=（2iiiiac=），表示在第 i 个介质中的热传导方程；hu 表示为 IV 层介质在贴近假人皮肤处的温度值；5u 表示为假人内部恒定温度，即 37oC；h 表示为空气与皮肤的转化系数，即空气冷却系数，其值为 8.6124。1111121212 uuuuxx=，表示介质 I 与介质 II 的两个耦合条件，获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：54445765713 2222322323uuuuxx=，表示介质 II 与介质 III 的两个耦合条件；3333343434 uuuuxx=，表示介质 III 与介质 IV 的两个耦合条件；(,0)37 1,2,3,4 iu xi=，表示热传导方程的初值条件；1(0,)75ut=，表示热传导方程的左边界 Dirichlet 边值条件；445()huh uux=，表示热传导方程的右边界 Robin 边值条件。5.1.5 基于热传导方程的温度分布模型的求解（1）求解方法分析）求解方法分析 在 5.1.4 中建立的温度分布模型（20）属于抛物型方程，因边值条件条件复杂难以求得解析解。本文采用有限差分法：将连续的定解区域用有限个离散点构成的网格来代替；把定解区域上的连续变量的函数用网格上定义的离散变量的函数来近似；把原方程和定解条件中的微商用差商来近似。最终，把原微分方程和定解条件用代数方程组来代替，即有限差分方程组。解此方程组就可以得到原问题在离散点上的近似值。常用的差分格式有：向前差分格式、向后差分格式和C-N差分格式（即Crank-Nicolson差分格式）。在向前差分格式中，要求时间步长与空间步长的比0.5r ，即时间步长要比空间步长小得多。若不满足此条件，极易发生解的爆破。在本文温度分布模型中，时间步长与空间步长达到 50，远远大于 0.5，因此不适合用向前差分格式。此外，向前差分格式为显式格式，难以进行多层介质的耦合。在 C-N 差分格式中，稳定性好精度高，计算在交叉点处的函数值，具体如下图 7 所示：图 7 C-N 差分格式图 在向后差分格式中，具有无条件稳定的特点，选取第 j+1 时间层相邻的三个结点进行 u 对 x 二阶偏导的近似，选取第 j+1 时间层与第 j 时间层的相邻两个结点进行对 u 对 t 一阶偏导的近似，具体如下图 8 所示：获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：54445765714 图 8 向后差分格式图（2）温度分布模型的具体求解过程：）温度分布模型的具体求解过程：经过上述三种差分格式的比较，本文采用向后差分格式求解，并借助迭代法求解。求解流程如下图 9 所示：图 9 温度分布模型求解流程图 Step1：对求解区域进行网格剖分。：对求解区域进行网格剖分。求解区域记为=(,)0,0 x txLtT，其中41=iiLL=且iL为第 i 层介质的厚度。对此求解区域进行网格剖分，具体如下：在空间维度上，将区间10,L做1M等分，将区间12,L L做2M等分，将区间23,L L做3M等分，将区间34,L L做4M等分。记1 1,2,3,4iillmMi=，。记iiiLxM=，ix 表示为第 i 层介质的空间步长。获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：54445765715 在时间维度上，将区间0,T作 n 等分，并记Ttn=，t为时间步长。用两组平行直线簇,1,1ijxximttjn=将分割成矩形网格，见下图 10：图 10 网格剖分图 Step2：建立隐式向后差分格式。：建立隐式向后差分格式。定义上的网格函数，411,11ijuuimjn=+，其中(,)ijijuu x t=，且有411,11imjn+。考虑温度分布模型（20）中的热传导方程中的右端项，略去小量项，用二阶中心差商代替 u 对 x 的偏导数，得到：211221(,)(,)2(,)(,)()ijijijijux tu xtu x tu xtxx+=+(23)考虑温度分布模型（20）中的热传导方程中的左端项，用一阶向前差商代替u 对 t 的一阶偏导数，得到：11(,)(,)(,)ijijijux tu x tu x ttt+=(24)建立差分格式的具体步骤如下：在结点处考虑不同介质下的热传导方程，得到介质内部的差分格式。第 i 层介质的热传导方程为：222,1,2,3,4kkkuuaktx=，则相应的差分格式为：1111,1,12220()iiiiikjkjk jkjk jkuuuuuatx+=(25)获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：54445765716 其中，ikju 表示第 k 层介质在(,)ijx t结点处的温度，1,2,3,4;1,2,kjm=。依据热流量密度的耦合条件，得到介质边界处的差分格式。在临界面123、处，关于热流量密度的耦合条件为：11,1,2,3kkkkkkuukxx+=则得到相应的差分格式为：111,1,1,21,111,1,2,3kkkkiiiik mk mkmkmkkkkuuuukxx+=(26)依据第 IV 层介质右侧 Robin 边值条件，得到边界4 上的差分格式。第 IV 层介质边界4处，满足温度分布模型（20）的条件445()huh uux=，可得到相应的差分格式为：444111144(37)iimmimuuh ux+=(27)综合，得到温度分布模型（20）的隐式向后有限差分近似如下：4441111,1,122111,1,1,21,1111111441120 1,2,3,4,1,2,()1,2,3(37)75 kkkkiiiiikjkjk jkjk jkiiiik mk mkmkmkkkkiimmimiuuuuuakjmtxuuuukxxuuh uxu+=1,11,1 1,2,37 1,2,3,4,1,2,1 kkkjjjk mkminukjmuu+=+=1,2,3,1,2,1 kjn=+(28)其中，（26）式中第一式、二式及三式在步骤、步骤及步骤中已做详细说明；（26）式第四式，为基于热传导方程的温度分布模型（20）中左边界Dirichlet 边值条件；（26）式第五式，为温度分布模型（20）中初值条件；（26）式第六式，为温度分布模型（20）中基于温度的耦合条件条件。Step3：将差分格式整理为代数方程组。：将差分格式整理为代数方程组。差分格式（26）中，第一式可整理为：获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：54445765717 ,11,1(12)k ijjjjkkik ikir ur uruu+=(29)其中，=1,2,3,4 ,2 ,2kimjm 。2()kkkatrx=，表示为第 k 层介质剖分的步长比。差分格式（26）式中，第二式可整理为：11111,11,211()0kkkiiikkkkk mk mkmkkkkuuuxxxx+=(30)差分格式（26）式中，第三式可整理为：4411444,14,44()37iimmuh uhxx+=(31)由可将差分格式（26）写成如下矩阵形式：1111211113111111111121111112133343,334344444444441+21+21+21 21 2iiimimmurrurrrrrruuxxxxuxxxxrrrrrrhxx+1333344441121 1111311113,11114,14,2114,14,1114,4,0iiiimiimiimmiimmiimmuruuuuuuuuuu+=Step4：利用追赶法解三对角线性方程组：利用追赶法解三对角线性方程组 在 step3 中，得到的代数矩阵为三对角线性方程组 AU=b，一般可用 matlab求逆命令求解 U，但因本题中系数矩阵规模较大，且除主对角线和两个次对角线之外，其余元素均为 0，因而求逆解方程组会导致浪费内存、程序运行缓慢。因此本文利用追赶法求解三对角线性方程组。获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：54445765718 算法步骤如下图 11 所示：图 11 追赶法求解三对角线性方程组 Step5：每一时间层上的温度分布情况：每一时间层上的温度分布情况 根据追赶法，求得不同时刻不同厚度下的温度分布情况。具体 5 个临界面处的温度分布情况见表 problem1.xlsx，如下为部分温分布情况。表 2 部分温度分布表 时间/s 临界面 I 临界面 II 临界面 III 临界面 IV 0 37 37 37 37 360 72.77 69.34 62.11 46.89 720 74.18 72.48 64.88 47.98 1080 74.29 72.73 65.10 48.07 1440 74.3 72.75 65.12 48.08 1800 74.3 72.75 65.12 48.08 2160 74.3 72.75 65.12 48.08 2520 74.3 72.75 65.12 48.08 2880 74.3 72.75 65.12 48.08 3240 74.3 72.75 65.12 48.08 表 2 中的数据，是四层介质的四个临界面相同时间间隔360ts=下的温度值。绘制成温度关于时间的二维曲线，如下图 12 所示：获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：54445765719 图 12 在不同时刻下临界面处的温度变化情况 由表格 2 和图 12 分析可知：I、II 介质在临界面处温度随时间缓慢降低，当时间 t=1440s 时，I 介质温度达到稳定值 74.3oC；当时间 t=1440s 时，II 介质温度达到稳定值 72.75oC。II 介质在临界面处的温度随时间降低，但比 I、II 介质降得快,当时间 t=1440s 时，温度达到稳定值 65,12oC；空气层 IV 介质温度降低且速度最快，当时间 t=1440s 时，温度达到稳定值 48.08oC。在转化系数 h=8.6125 的条件下，不同时刻不同厚度的温度变化情况如下图所示：图 13 不同时刻不同厚度的温度变化图 根据图 13，将表格 problem1.xlsx 中的离散数据在 matlab 中绘制成三维曲面图。其中，x 轴表示对时间维度作剖分 05400s；y 轴表示对空间维度作剖分，将I 介质的厚度1L做 6 等分，将 II 介质的厚度2L做60等分，将 III 介质的厚度3L做获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：54445765720 36 等分，将区间4L做 50 等分；z 轴表示温度的变化情况。曲面图形在不同时刻临界面处的温度与表格反映的内容一致。5.1.6 问题一模型的检验 问题一中，通过求解四层耦合介质的温度分布模型，可得到任意厚度的介质在每一时刻的温度。将 IV 层介质右侧与皮肤直接接触的临界面的温度与题目中附件 2 假人皮肤外侧的测量温度进行对比，可检验求解结果的正确性。假设时间层it下通过模型解得的温度为iu，附件 2 所给的温度为iv，本文定义求解结果与题目所给数据之间的偏差指数 f，用于表示求解得到的温度与附件2 温度的差值平方和的平均值，540121()5401iiiuvf=(32)对 c 的数值开根号进行标准化处理，得到标准化后的偏差 540121()5401iiiuvf=(33)通过利用 matlab 编程计算，得到标准化后的偏差为0.459=3=2 A(i,i-1)=-r1;end end A(m1,m1)=(lam_1/derta_x1+lam_2/derta_x2);A(m1,m1-1)=-lam_1/derta_x1;附录二：在已知空气系数范围下求解问题一中空气系数的精确值（matlab 程序）clear;%清除工作区变量 clc;%清屏 close all;%关闭所有图形窗口 z=;for h=8.61:0.0001:8.63%确定空气交换系数%材料参数输入 m1=6;m2=60;m3=36;m4=50;%分别对四种介质分割 m=m1+m2+m3+m4;%介质分割和 n=5400;%对时间分割 t=5400;%总时长 l1=0.6/1000;l2=6/1000;l3=3.6/1000;l4=5/1000;%四种材料厚度 lam_1=0.082;lam_2=0.37;lam_3=0.045;lam_4=0.028;%四种材料的热传导率 de_1=300;de_2=862;de_3=74.2;de_4=1.18;%四种材料的密度 c1=1377;c2=2100;c3=1726;c4=1005;%四种材料的比热容%计算热扩散率 a1=lam_1/(c1*de_1);%I 层材料的热扩散率 a2=lam_2/(c2*de_2);%II 层材料的热扩散率 a3=lam_3/(c3*de_3);%III 层材料的热扩散率 a4=lam_4/(c4*de_4);%IV 层材料的热扩散率%材料长度分割和时间步长分割求解 derta_x1=l1/m1;%I 层材料的分割长度 derta_x2=l2/m2;%II 层材料的分割长度 derta_x3=l3/m3;%III 层材料的分割长度 derta_x4=l4/m4;%IV 层材料的分割长度 derta_t=t/n;%时间步长分割%计算各层介质剖分的步长比 获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：54445765731 r1=derta_t/derta_x12*a1;%第 I 层介质剖分的步长比 r2=derta_t/derta_x22*a2;%第 II 层介质剖分的步长比 r3=derta_t/derta_x32*a3;%第 III 层介质剖分的步长比 r4=derta_t/derta_x42*a4;%第 IV 层介质剖分的步长比 u=zeros(m+1,n+1);%定义四层耦合介质温度分布矩阵%初始条件和边界条件 u(:,1)=37;%初始条件 u(1,:)=75;%边界条件%差分格式的系数矩阵的构造 A=zeros(m,m);for i=1:m1-1 A(i,i)=1+2*r1;A(i,i+1)=-r1;if i=2 A(i,i-1)=-r1;end end A(m1,m1)=(lam_1/derta_x1+lam_2/derta_x2);A(m1,m1-1)=-lam_1/derta_x1;A(m1,m1+1)=-lam_2/derta_x2;for i=m1+1:m1+m2-1 A(i,i)=1+2*r2;A(i,i+1)=-r2;A(i,i-1)=-r2;end A(m1+m2,m1+m2)=(lam_2/derta_x2+lam_3/derta_x3);A(m1+m2,m1+m2-1)=-lam_2/derta_x2;A(m1+m2,m1+m2+1)=-lam_3/derta_x3;for i=m1+m2+1:m1+m2+m3-1 A(i,i)=1+2*r3;A(i,i+1)=-r3;A(i,i-1)=-r3;end A(m1+m2+m3,m1+m2+m3)=(lam_3/derta_x3+lam_4/derta_x4);A(m1+m2+m3,m1+m2+m3-1)=-lam_3/derta_x3;A(m1+m2+m3,m1+m2+m3+1)=-lam_4/derta_x4;for i=m1+m2+m3+1:m1+m2+m3+m4-1 获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：54445765732 A(i,i)=1+2*r4;A(i,i-1)=-r4;A(i,i+1)=-r4;end A(m,m)=h+lam_4/derta_x4;A(m,m-1)=-lam_4/derta_x4;%构造右端项 for k=2:n+1 b=zeros(m,1);for i=2:m-1 b(i,1)=u(i+1,k-1);end b(1,1)=u(2,k-1)+r1*u(1,k);b(m1,1)=0;b(m1+m2,1)=0;b(m1+m2+m3,1)=0;b(m,1)=37*h;%追赶法求解 bb=diag(A);aa=0,diag(A,-1);c=diag(A,1);N=length(bb);L=zeros(N);uu0=0;y0=0;aa(1)=0;L(1)=bb(1)-aa(1)*uu0;y(1)=(b(1)-y0*aa(1)/L(1);uu(1)=c(1)/L(1);for i=2:(N-1)L(i)=bb