[Matlab]2018国赛 A题思路以及Matlab代码实现

一些关于CUMCM的思考

去年参加的CUMCM，觉得代码闲置在附录里边比较可惜就打算分享出来了。
国赛的题目想要求出结果并不是很难，使得论文的写作以及验证等环节愈发重要，语言表达是否精练以及模型的评价与改进是否合理有效成为了评价参赛作品的一项重要指标。由于去年作品要求通过客户端上传，一定要注意上传支撑材料的大小不要超过20M，注意参赛时间安排以及qq账号的清洁（！）。
本文给出的思路，旨在辅助对于题目的理解并帮助解释贴出的Matlab代码，不代表最优做法。本文思路及代码仅供参考，思路中涉及到的公式、方法以及详细步骤等请移步至参考文献。代码中变量命名比较乱，仅供参考。

CUMCM2018 A题

A题 高温作业专用服装设计
在高温环境下工作时，人们需要穿着专用服装以避免灼伤。专用服装通常由三层织物材料构成，记为I、II、III层，其中I层与外界环境接触，III层与皮肤之间还存在空隙，将此空隙记为IV层。
为设计专用服装，将体内温度控制在37ºC的假人放置在实验室的高温环境中，测量假人皮肤外侧的温度。为了降低研发成本、缩短研发周期，请你们利用数学模型来确定假人皮肤外侧的温度变化情况，并解决以下问题：
(1) 专用服装材料的某些参数值由附件1给出，对环境温度为75ºC、II层厚度为6 mm、IV层厚度为5 mm、工作时间为90分钟的情形开展实验，测量得到假人皮肤外侧的温度（见附件2）。建立数学模型，计算温度分布，并生成温度分布的Excel文件（文件名为problem1.xlsx）。
(2) 当环境温度为65ºC、IV层的厚度为5.5 mm时，确定II层的最优厚度，确保工作60分钟时，假人皮肤外侧温度不超过47ºC，且超过44ºC的时间不超过5分钟。
(3) 当环境温度为80 时，确定II层和IV层的最优厚度，确保工作30分钟时，假人皮肤外侧温度不超过47ºC，且超过44ºC的时间不超过5分钟。

附件1. 专用服装材料的参数值
附件2. 假人皮肤外侧的测量温度

对于题目的分析

对于该题目可以建立如下假设：
1.假设外界环境与假人体内的气温保持不变。
2.假设辐射传热的影响可以忽略。
3.假设人的体型对传热的影响可以忽略，人体可以等效为圆柱体。
4.假设热量在材料交接交界面处的传导方向为交界面法线方向。

这是那种第一问做出来之后第二问第三问都迎刃而解的题目，找准建立模型的物理理论依据，再确保了解决偏微分方程的求解手段就可以得到答案了。

第一问

1.1 总体思路

问题1要求在给定外界温度和高温工作服各层材料厚度的条件下求解工作服的温度分布情况。根据热学分析可以建立出关于各层材料传热情况的偏微分方程组模型。工作服的初始温度以及工作过程中材料各交界面的温度作为各层偏微分方程边界条件。由于此微分方程组模型较为复杂，运用有限差分法将微分方程化为差分方程，再将差分方程组以矩阵的形式写出，并用追赶法（Thomas算法）求解差分方程组模型可得出工作服的温度分布。
针对微分的边界条件确定，最内层材料与假人皮肤交界面的温度同最外层材料与外界空气交界面的温度依据牛顿冷却定律通过外界环境气温、假人内部温度计算，运用牛顿冷却定律时需要确定对流换热系数这一未知参数，此外工作服的初始温度也是未知量。对于这两个未知量，先大致估计这两个未知量的范围，再按照一定步长遍历搜索。对于每个遍历到的值，代入模型中运用有限差分法计算，得出相应的皮肤表面的温度，与附件二给出的实测值进行对比，取残差和最小时的值作为初始温度和对流换热系数的理想值。

1.2 模型求解

1.2.1 偏微分方程的建立
热传导在不同温度的物体相互接触或者同一物体不同部分具有不同温度时产生，符合本题中几种材料互相接触传递热量的情形。将假人抽象为多层互相包裹的同轴圆柱体后，可以对上述三维几何模型进行降维处理，建立一维几何模型。且该题可忽略对流传热与辐射传热的影响，物体内部的热传递由傅里叶定律确定，而物体表面和内部的热量流动又满足能量守恒定律，依据题目给定部分边界条件可建立一维热传导模型。
物体在边界上与其他传热介质接触时，物体表面与介质的温度往往不相同，此时运用牛顿冷却定律：单位时间从物体流入介质中的热量和两者的温度差成正比，可确定热传导的第三类边界条件。
1.2.2 有限差分法——化为差分方程
生成网格–>微分化差分–>差分方程矩阵化–>边界条件矩阵化–>合成矩阵差分方程组（这一步比较考验耐心）–>追赶法
具体做法可参考这篇文章：史策.热传导方程有限差分法的MATLAB实现[J].咸阳师范学院学报,2009,24(04):27-29+36.
1.2.3 追赶法——解差分方程
实际建立出方程组的系数矩阵为三对角矩阵，即非零元素都集中在主对角线和其相邻两条次对角线上，故容易分解成两个三角阵，因此可以用追赶法求解。
追赶法的理论可参考这篇博文：追赶法求解三对角方程组
1.2.4 输出结果至Excel文件。

1.3 Matlab实现

代码如下：

%%为应用牛顿冷却定律确定边值条件，需要在测试范围内遍历k21取值，计算边值条件与附件二中数据进行比对；
%%选取得到最小平均残差，并将最大残差控制在一定范围内的k21值（此程序得到k21=-0.0046）；
clc;
clear;
load('tX4.mat');%导入附件2假人皮肤外侧温度变化数据，事先已存储为变量，在此问作为参考值使用
%%常量及边界条件定义
d1=0.6; d2=6; d3=3.6; d4=5;%单位为mm
k1=0.082; k2=0.370; k3=0.045; k4=0.028;%热传导率
c1=1377; c2=2100; c3=1726; c4&#