1、全国第 X 届研究生数学建模竞赛题 目 恒温箱温度变化的数学建模摘 要恒温箱是航空、汽车、家电、科研等领域必备的测试设备，用于测试和确定电工、电子及其他产品及材料进行高温试验的温度环境变化后的参数及性能。 。因此,对温度进行测量和控制也是科学实验和工业生产中经常需要解决的重要问题。随着计算机技术和控制理论的发展,以及对产品质量要求的提高,人们对高精度的温度测量和控制的要求也越来越高。电加热设备温度特性复杂,其温度的测量和控制亦显得尤为重要和复杂。因此，本问题具有很强的研究背景和实际应用价值。 为描述保温箱温度的变化规律，文章首先对环境温度进行三次样条插值，将环境温度的采样间隔变成和恒温箱一样，
2、也就是一分钟。然后根据能量平衡原理和热力学知识列出恒温箱和隔热层从室温开始加热这一阶段的热量平衡方程组，方程组一共有6 个未知量，取 78 个分钟的数据带入方程组中得到 156 个方程式，用 matlab 求解这个超定方程组，求得 6 个未知量的值，也就是恒温箱和隔热层的各项参数值。随后将剩余两个分钟的数据带入求得的方程组中进行检验，误差在接受范围内，说明建模合理。接下来，将热量平衡代数方程组变为功率平衡微分方程组，在考虑恒温箱控制条件的情况下用龙格库塔法进行求解，得到接下来恒温箱和隔热层的温度变换规律，画出温度变换曲线。最后对恒温箱的性能进行了评价。关键词：恒温箱 隔热层 热学平衡方程式 超
3、定方程组 龙格库塔法参赛队号 队员姓名 陈文哲 管俊 孙鹏伟 参赛密码 （由组委会填写）XX 大学承办 恒温箱温度变化的数学建模一、问题重述恒温箱是航空、汽车、家电、科研等领域必备的测试设备，用于测试和确定电工、电子及其他产品及材料进行高温试验的温度环境变化后的参数及性能。 。因此,对温度进行测量和控制也是科学实验和工业生产中经常需要解决的重要问题。随着计算机技术和控制理论的发展,以及对产品质量要求的提高,人们对高精度的温度测量和控制的要求也越来越高。电加热设备温度特性复杂,其温度的测量和控制亦显得尤为重要和复杂。因此，本问题具有很强的研究背景和实际应用价值。 模型要求热源即加热器以恒定的速率
4、加热恒温箱。恒温箱体内装有一个自动的温度控制器，在温度高于 68.2 华氏度会自动关闭加热系统, 在温度低于 67.8 华氏度时会自动打开加热系统，但温度控制器读取温度有时间间隔，设定温度控制器读取温度有时间间隔为每 1 分钟 1 次。 根据上述分析，本文拟解决如下问题：（1） 、根据已有的资料数据，利用所学知识及相关参考资料，通过研究温度传播规律,确定加热方案,并建立恒温箱系统的热数学模型。 ；（2） 、建立恒温箱体的温度变化模型。并根据小型试验数据及相关数学、物理方法，验证模型准确度和适用性；（3） 、对恒温箱产品的质量优劣做出评价；（4） 、对建立的模型进行误差分析，通过模型参数调节优化
5、模型效果，并根据修改模型预测最有效数学模型。根据误差原因分析等对未来需要进行的试验和研究工作提出了一些建议；（5） 、根据相似准则设计相关实验，对已建立模型进行推广和扩充。以达到理论-实验- 工程三者的互相验证。 二、基本假设（1） 、环境温度不受热源、恒温箱和隔热层的影响；（2） 、忽略空气与恒温箱及隔热层之间的对流换热、热辐射；（3） 、所有的导热过程仅考虑稳态导热；（4） 、所有物体的物性不随温度变化；（5） 、环境温度始终低于 67.8 华氏度；（6） 、恒温箱体与箱顶隔热层温度始终不低于外界环境温度；三、符号约定1c恒温箱比热容m恒温箱质量boxTt时间内恒温箱温度增量1恒温箱箱体热
6、导率box恒温箱某一时刻温度airT环境温度1A恒温箱与外界环境的接触表面积恒温箱体的厚度gerT隔热层某时刻温度2A恒温箱体和隔热层的接触表面积q热源 t时间内所释放的热量2c隔热层的比热容m隔热层的质量gerTt时间隔热层温度的增量2隔热层的热导率隔热层的厚度6四、模型的建立与求解基本理论（1）热传导只要有温差,热量就会自发地从高温物体传到低温物体,所以传热是一种最常见的物理过程。实际上,热是一种扩散效应,物体将能量通过热的方式释放出来,因此热其实是能量的传递。热的传递方式主要包括以下三种:热传导、热对流和热辐射。在对恒温箱的数学建模中，我们忽略热对流和热辐射，仅考虑热传导的主要作用。热传
7、导是指当物体内具有温差或者不同温度的物体相互接触时,在各部分之间没有发生相对宏观位移情况下,通过物质微观粒子(分子、原子或自由电子)的热运动而进行的热量传递。导热是固体内唯一的热量传递方式;在气体和液体中 ,尽管也同样存在着导热现象。图 2 恒温箱壁面导热示意图如图所示恒温箱壁，一块厚度为 ，横截面积为 的平板，两侧表面分别维持着均匀A温度 和 。实验表明：单位时间内从表面 1 传递到表面 2 的热量 与壁面间的温1tw2 Q差 以及导热面积 成正比，而与平板的厚度 成反比，即AtQ=方程中的比例系数 反映了材料导热能力的大小，叫做导热系数或者热导率，单位是 。导热率是一个物性参数，不同材料的
8、导热率不同，即使是同一种材料，W(mK)在不同温度下也具有不同的热导率。为了确保恒温箱内的光谱仪能够安全可靠的工作,我们需要将恒温箱内多余的热量排放到外面，从而保证恒温箱内的光谱仪能够获得适宜的工作温度，这是系统级热设计。总的来说,恒温箱的热设计可以分为 3 个主要方面。第一：组件级热设计，即对元器件级的热设计。第二：封装级设计，即对于电子控制模块、散热设备、PCB 板的热设计。第三：系统级设计，即对于恒温箱内胆及外壳的热设计。系统级的热设计主要研究电子设备所处环境的温度对其产生的影响，环境温度是系统级热分析的重要边界条件，其热设计是采取措施控制环境温度，使电子设备在适宜的温度环境下进行工作。
9、随着电子技术的发展，尤其是微电子技术的发展，电子元器件和设备的尺寸正迅速缩小，而功率却一直在增大，使得单位体积容纳的热量越来越多，特别是在恶劣环境下电子产品以及国防电子仪器要求必须达到良好的电磁兼容性能及三防性能等。（2）三次样条差值将环境温度每隔 15min 采集到的数据利用三次样条差值成每隔 1min 的数据。算法原理：设在区间 上给定 个节点 ，在节点 处的函数值,ab1n01()i nxaxb ix为 。若函数 满足如下三条：()0,)iiyfx )S（1） 在每个子区间 上， 是三次多项式；1,(,2ixn ()Sx（2） ；()iiSy（3） 在区间 上， 的二阶导数 连续。,ab
10、(Sx()x则称 为函数 在区间 上的三次样条插值函数。()x)yf,ab子区间 上的 的表达式为：1,2i n ()x33211112()() )666(),i i iii iiiiii iixxhxSMyMhhyx关于参数 的方程组（三弯矩方程组）：i11 11112,2,6()6,ii iiiii iii ii ii iii dnhhhxyydf 牛顿插值多项式： 0100101(),(),(), .nk kn nNxffxfxxf 构造牛顿插值多项式首先列出差商表，进而由差商表写出牛顿插值多项式。阶差商为：n0101101,nnnfxfxfx 零阶差商和一阶差商： 11(),iiiii
11、ifxffxyfxf（3）龙格库塔法算法原理：考虑用函数 f(x,y)在若干点上的函数值的线性组合来构造近似公式，构造时要求近似公式在(xi,yi)处的 Taylor 展开式与解 y(x)在 xi 处的 Taylor 展开式的前面几项重合，从而使近似公式达到所需要的阶数。既避免求高阶导数,又提高了计算方法精度的阶数。或者说,在xi,xi+1这一步内多计算几个点的斜率值，然后将其进行加权平均作为平均斜率，则可构造出更高精度的计算格式，这就是龙格库塔（Runge-Kutta ）法的基本思想。用类似上述的处理方法，只需在区间xi,xi+1上用四个点处的斜率加权平均作为平均斜率 K*的近似值，构成一系
12、列四阶龙格库塔公式。具有四阶精度，即局部截断误差是 O(h5)。对于高阶常微分方程的初值问题： ()(1) (1)000, ,()mmmyfxyaxbay将其转化为一阶微分方程组求解： 1231 12 (1)1000(,)mmmyyfxyaya 标准的四级四阶龙格库塔法的向量形式： 11234213243()6,)(),()iiiiiiiiiyKhfxyfxKhK其分量形式： ,1,123411212 223143233()6;,)(,),;,( )jiijjjjjjjiini njjiii ijjiiinijjiiiiyKhfxyKyyjfxKhyKy n（4）超定方程的最小二乘解小二乘法广
13、泛地应用于工程计算中，用最小二乘法消除（平滑）误差，用最小二乘法从有噪声的数据中提取信号，从海量数据中找出数据变化的趋势，。甚至利用简单函数计算复杂函数的近似值，我们并不期望它的近似值多么精确（事实上很多时候也不用很精确） ，尽管如此还是希望计算出的近似数据与原始数据之间有相似之处。如果从线性代数角度来理解最小二乘法，实际上是将一个高维空间的向量投影到低维子空间所涉及的工作。超定方程组的最小二乘解当方程组 GX=b 的方程数多于未知数个数时，对应的系数矩阵 G 的行数大于列数，此时方程组被称为是超定方程组。设 G=(giu)mn，当 mn 时即所谓的高矩阵，绝大多数情况下，超定方程组没有古典意
14、义下的解。超定方程组的最小二乘解是一种广义解，是指使残差 r = b GX 的 2-范数达取极小值的解，即 22*|i| XbGXbmR超定方程组的最小二乘解是指正规方程组GTGX=GTb的解。如果系数矩阵(G TG)可逆，则正规方程组有唯一解。此时，最小二乘解可以形式地写为如下形式X=(GTG)-1GTb两种常用的方法如下1用对称矩阵的三角分解法解正规方程组 GTGX = GTb；记 A=GTG，则 A 是对称矩阵，由三角分解 A = L D LT，其中 L 是下三角矩阵，D 是对角矩阵。将这一算法写过三个过程：解下三角方程组：LY 1 = GTb；解对角方程组：DY 2 = Y1 ；解上三
15、角方程组：L TY3 = Y22用矩阵的 QR 分解直接求解超定方程组由 QR 分解（正交三角分解）G=QR ，其中 Q 是正交矩阵，R 是上三角矩阵。将 QR 分解代入最小二乘解表达式中，得X=( R T QTQR)-1(QR)Tb = R-1QTb由此可知，用分解方法求超定方程组的最小二乘解只需求解上三角方程组RX=QTb模型 1 热源给恒温箱体加热，恒温箱体给箱顶隔热层加热；恒温箱及箱顶隔热层向外界空气散热。此时，恒温箱未达到 68.2 华氏度，恒温箱体温度高于箱顶隔热层，如图恒温箱体箱顶隔热层热源图 3 模型 1 恒温箱热量传播方向图模型 2恒温箱达到 68.2 华氏度后，箱顶隔热层未
16、达到 67.8 华氏度，热源停止加热；恒温箱体给箱顶隔热层加热，恒温箱及箱顶隔热层向外界空气散热，如图恒温箱体箱顶隔热层热源图 4 模型 2 恒温箱热量传播方向图模型 3恒温箱及箱顶隔热层温度均高于 67.8 华氏度，热源停止加热；恒温箱温度低于箱顶隔热层，箱顶隔热层给恒温箱加热；恒温箱及箱顶隔热层向外界空气散热，如图恒温箱体箱顶隔热层热源图 5 模型 3 恒温箱热量传播方向图模型 4恒温箱体温度低于 67.8 华氏度后，此时箱顶隔热层温度高于 67.8 华氏度，热源加热，热源与箱顶隔热层给恒温箱体加热，直到箱顶隔热层与恒温箱体温度相等。恒温箱体箱顶隔热层热源图 6 模型 4 恒温箱热量传播方
17、向图模型建模根据已测的环境温度、恒温箱体和箱顶隔热层温度数据，在对每隔 15min 测量一次的环境温度数据进行三次样条差值后得到每隔 1min 一个数据，如图 7 所示。初始时刻，恒温箱体与隔热层温度等于环境温度。此时热源给恒温箱体加热，恒温箱体给箱顶隔热层加热；恒温箱及箱顶隔热层向外界空气散热。这段过程，恒温箱未达到68.2 华氏度，恒温箱体温度高于箱顶隔热层。根据能量守恒，恒温箱体增加的热量与恒温箱对环境空气和箱顶隔热层的散热之和等于热源提供的热量。 111121()/()/boxboxairboxgercmTAtTAtqt（1）恒温箱体对箱顶隔热层的散热等于箱顶隔热层增加的热量与箱顶隔热
18、层对环境空气的散热之和。212122()/()/gerboxgergeraict t （2）式中各变量含义分别为 c是恒温箱比热容， m是恒温箱质量， boxT是 t时间内恒温箱温度增量， boxT是 t时间内恒温箱温度增量， 1是恒温箱箱体热导率， box是恒温箱某一时刻温度， air是环境温度， 1A是恒温箱与外界环境的接触表面积， 1是恒温箱体的厚度， ger是隔热层某时刻温度， 2是恒温箱体和隔热层的接触表面积，q是热源 t时间内所释放的热量， 2c是隔热层的比热容， 2m是隔热层的质量，gerT是 时间隔热层温度的增量， 是隔热层的热导率， 是隔热层的厚度。将上式进一步简化为()()
19、boxboxairboxgerAtBTCtTt （3）gergeri rDE （4）式中 12qA， 12cm， 12A， 21cmD， 21E； t定为 1min取差值后的 80 组数据中的前 78 组数据，求解超定方程中的 A、B 、C、D 、E 五个参数。A=22.4495B=63.9245C=0.0366D=122.5048E=0.4669将能量守恒方程变换为功率守恒方程：()()boxboxgerboxairTBATCTt （5）gerboxgergeraiDEt （6）代入参数值可得：63.9245.4951.0360.3box boxgerairdTTTt （7）12.50481
20、.4690.6gerboxgerairdTTt （8）将第 79 和 80 组数据带入上面的方程式，绝对误差分别为 0.5624，0.3112 和0.602，0.331，在接受范围内，验证了方程式的正确性。图 7 环境温度实际测量与三次采样根据式（7）和（8）反复迭代计算恒温箱与箱顶隔热层的温度，计算前 80min 的数据，可以得到恒温箱体与箱顶隔热层温度变化曲线，如图 8 和图 9 所示。图 8 前 80min 恒温箱温度拟合曲线图 9 前 80min 箱顶隔热层温度拟合曲线由图可知，由建模得到的恒温箱体与箱顶隔热层温度变化，与温度实际测量曲线拟合得很好。再根据模型 2、3、4 讨论的情况采
21、用比较恒温箱、箱顶隔热层温度与临界温度的高低，再选择模型，具体程序见附录。由此可以得到恒温箱体和箱顶隔热层温度的变化曲线。仿真结果见图 10 和图 11。图 10 恒温箱与箱顶隔热层温度变化曲线图 11 恒温箱温度变化曲线由建模结果可知，箱顶隔热层始终低于 68 华氏度，箱顶隔热层与环境温度变化趋势相近。恒温箱体在刚开始持续受热，约在 260 分钟加热到 68 华氏度，然后在 67.8华氏度与 68.2 华氏度附近来回震荡。对比图 10 和图 11 可知，随着环境温度的升高，恒温箱体温度震荡频率降低。环境温度决定了恒温箱温度的震荡频率，即热源开关的频率。恒温箱箱体温度趋近于稳定后（300 分钟
22、1440 分钟） ，加热时间 Theat 约为 317 分钟，不加热时间约为 814 分钟,加热时间与不加热时间比值为 317/814=0.4017。恒温箱箱体温度趋近于稳定后（300 分钟1440 分钟） ，最高温度 68.4457 华氏度，最低温度 67.6212 华氏度。模型的改进 恒温箱温度控制算法设计建模前提是热源提供的速率一致，可以采用 PID 控制调节热源的加热速率。使恒温箱加热频率降低，响应速率变快。PID 控制算法简单、鲁棒性能好、可靠性高，被广泛应用于工业过程并取得了良好的控制效果。下图是 PID 控制的原理框图。将给定值 与实际输出值 相减，得到控制偏差：rtctertc
23、将偏差经过比例、积分、微分作用后线性相加，得到 PID 控制器输出的控制量，输入到被控对象进行控制： 01tp DIdetutKetT其传递函数形式为： 1cppIGssT式中， 是比例系数， 是积分时间常数， 是微分时间常数。pKITD比例环节：即时成比例地反应控制系统的偏差信号 ，偏差一旦产生，控制器et立即动作以减小误差。比例作用大，可以加快调节，减少误差。 （比例也不能过大，过大的比例会使系统的稳定性下降，甚至造成系统的不稳定。 ）积分环节：主要作用是消除静差，提高系统的控制精度。积分作用的强弱取决于积分时间常数 ，其值越大，积分作用越弱；反之越强，引入积分调节可使系统稳定IT性下降，
24、动态响应变慢。微分环节：能反应偏差信号的变化趋势（变化速率） ，并能在偏差信号变得太大之前引入一个修正信号，从而减小系统的超调，加快系统的过渡过程，减小调节时间。但微分作用对噪声干扰有放大作用。 参考文献1 李嘉林. 基于 TEC 制冷加热管制热的恒温箱设计D.复旦大学,2012. 2 陈国初,张琳 ,郝宁眉等. 恒温恒液位系统的动态机理建模与仿真J.青岛大学学报(工程技术版)，2003,18(2):46-50.3 叶庆银,刘斌,臧润清等.小型恒温箱恒温性能的理论及实验研究J.制冷与空调,2007,7(4):48-50,81.DOI:10.3969/j.issn.1009-8402.2007.
25、04.012.4 傅秦生. 热工基础与应用M.北京：机械工业出版社，2007.5 李乃成，梅立泉 . 数值分析M.北京：科学出版社，2011.附录三次样条插值代码：hwx,grc=textread(hwxwd.txt,%f%f,headerlines,1)wmwd=textread(wmwd.txt,%f)figureplot(t1,wmwd)T_interp=0:1:1440;wmwd1=interp1(t1,wmwd,T_interp,linear);wmwd2=interp1(t1,wmwd,T_interp,spline);plot(t1(1:7),wmwd(1:7),o,T_inte
26、rp(1:90),wmwd1(1:90),r+,T_interp(1:90),wmwd2(1:90),g*)温度曲线拟合代码：调用函数：Tair=wmwd2;wdqx(0,1000,34.6256,34.6256,Tair,1,1000)龙格库塔法数值积分函数：function hwxwd,grcwd = wdqx( x0,xn,Tbox0,Tgere0,Tair,h,n )%定义初值%h:步长；n:节点个数x=x0:h:xn;y(1,1)=Tbox0;%赋恒温箱、隔热层温度初值y(2,1)=Tgere0;i=1;while i67.8)K1(:,1)=h*func2(y(1,i),y(2,i
27、),Tair(i),h*func3(y(1,i),y(2,i),Tair(i);K2(:,1)=h*func2(y(1,i)+K1(1,1)/2,y(2,i)+K1(2,1)/2,Tair(i),h*func3(y(1,i)+K1(1,1)/2,y(2,i)+K1(2,1)/2,Tair(i);K3(:,1)=h*func2(y(1,i)+K2(1,1)/2,y(2,i)+K2(2,1)/2,Tair(i),h*func3(y(1,i)+K2(1,1)/2,y(2,i)+K2(2,1)/2,Tair(i);K4(:,1)=h*func2(y(1,i)+K3(1,1),y(2,i)+K3(2,1
28、),Tair(i),h*func3(y(1,i)+K3(1,1),y(2,i)+K3(2,1),Tair(i); y(:,i+1)=y(:,i)+1/6*(K1+2*K2+2*K3+K4);%不加热循环i=i+1;end endhwxwd=y(1,:);grcwd=y(2,:);end微分方程 1：function delta_hwx= func1( Tbox,Tgere,Tair )%加热方程时的恒温箱温度微分方程delta_hwx=1/63.9245*(22.4495-1.0366*Tbox+Tgere+0.0366*Tair);end微分方程 2：function delta_hwx1= func2( Tbox,Tgere,Tair )%不加热方程时的恒温箱温度微分方程delta_hwx1=1/63.9245*(-1.0366*Tbox+Tgere+0.0366*Tair);end微分方程 3：function delta_grc= func3( Tbox,Tgere,Tair )%隔热层温度微分方程delta_grc=1/122.5048*(Tbox-1.4669*Tgere+0.4669*Tair);end
