恒温箱数学建模

全国第X届研究生数学建模竞赛

题目恒温箱温度变化的数学建模

摘要

恒温箱是航空、汽车、家电、科研等领域必备的测试设备，用于测试和确定电工、电子及其他产品及材料进行高温试验的温度环境变化后的参数及性能。。因此,对温度进行测量和控制也是科学实验和工业生产中经常需要解决的重要问题。随着计算机技术和控制理论的发展,以及对产品质量要求的提高,人们对高精度的温度测量和控制的要求也越来越高。电加热设备温度特性复杂,其温度的测量和控制亦显得尤为重要和复杂。因此，本问题具有很强的研究背景和实际应用价值。

为描述保温箱温度的变化规律，文章首先对环境温度进行三次样条插值，将环境温度的采样间隔变成和恒温箱一样，也就是一分钟。然后根据能量平衡原理和热力学知识列出恒温箱和隔热层从室温开始加热这一阶段的热量平衡方程组，方程组一共有6个未知量，取78个分钟的数据带入方程组中得到156个方程式，用matlab求解这个超定方程组，求得6个未知量的值，也就是恒温箱和隔热层的各项参数值。随后将剩余两个分钟的数据带入求得的方程组中进行检验，误差在接受范围内，说明建模合理。

接下来，将热量平衡代数方程组变为功率平衡微分方程组，在考虑恒温箱控制条件的情况下用龙格库塔法进行求解，得到接下来恒温箱和隔热层的温度变换规律，画出温度变换曲线。最后对恒温箱的性能进行了评价。

关键词：恒温箱隔热层热学平衡方程式超定方程组龙格库塔法

参赛队号

队员姓名陈文哲管俊孙鹏伟

参赛密码 （由组委会填写） XX大学承办

恒温箱温度变化的数学建模

一、问题重述

恒温箱是航空、汽车、家电、科研等领域必备的测试设备，用于测试和确定电工、电子及其他产品及材料进行高温试验的温度环境变化后的参数及性能。。因此,对温度进行测量和控制也是科学实验和工业生产中经常需要解决的重要问题。随着计算机技术和控制理论的发展,以及对产品质量要求的提高,人们对高精度的温度测量和控制的要求也越来越高。电加热设备温度特性复杂,其温度的测量和控制亦显得尤为重要和复杂。因此，本问题具有很强的研究背景和实际应用价值。

模型要求热源即加热器以恒定的速率加热恒温箱。恒温箱体内装有一个自动的温度控制器，在温度高于68.2华氏度会自动关闭加热系统, 在温度低于67.8华氏度时会自动打开加热系统，但温度控制器读取温度有时间间隔，设定温度控制器读取温度有时间间隔为每1分钟1次。

根据上述分析，本文拟解决如下问题：

（1）、根据已有的资料数据，利用所学知识及相关参考资料，通过研究温度传播规律,确定加热方案,并建立恒温箱系统的热数学模型。； （2）、建立恒温箱体的温度变化模型。并根据小型试验数据及相关数学、物理方法，验证模型准确度和适用性；

（3）、对恒温箱产品的质量优劣做出评价；

（4）、对建立的模型进行误差分析，通过模型参数调节优化模型效果，并根据修改模型预测最有效数学模型。根据误差原因分析等对未来需要进行的试验和研究工作提出了一些建议；

（5）、根据相似准则设计相关实验，对已建立模型进行推广和扩充。以达到理论-实验-工程三者的互相验证。

二、基本假设

（1）、环境温度不受热源、恒温箱和隔热层的影响； （2）、忽略空气与恒温箱及隔热层之间的对流换热、热辐射； （3）、所有的导热过程仅考虑稳态导热； （4）、所有物体的物性不随温度变化； （5）、环境温度始终低于67.8华氏度；

（6）、恒温箱体与箱顶隔热层温度始终不低于外界环境温度；

三、符号约定

c1恒温箱比热容 m1恒温箱质量

?Tbox?t时间内恒温箱温度增量

?1恒温箱箱体热导率 Tbox恒温箱某一时刻温度 Tair环境温度

A1恒温箱与外界环境的接触表面积

?1恒温箱体的厚度 Tgere隔热层某时刻温度

A2恒温箱体和隔热层的接触表面积

q热源?t时间内所释放的热量

c2隔热层的比热容 m2隔热层的质量

?Tgere?t时间隔热层温度的增量

?2隔热层的热导率 ?2隔热层的厚度

四、模型的建立与求解

基本理论

（1）热传导

只要有温差,热量就会自发地从高温物体传到低温物体,所以传热是一种最常见的物理过程。实际上,热是一种扩散效应,物体将能量通过热的方式释放出来,因此热其实是能量的传递。热的传递方式主要包括以下三种:热传导、热对流和热辐射。在对恒温箱的数学建模中，我们忽略热对流和热辐射，仅考虑热传导的主要作用。

热传导是指当物体内具有温差或者不同温度的物体相互接触时,在各部分之间没有发生相对宏观位移情况下,通过物质微观粒子(分子、原子或自由电子)的热运动而进行的热量传递。导热是固体内唯一的热量传递方式;在气体和液体中,尽管也同样存在着导热现象。

图2 恒温箱壁面导热示意图

如图所示恒温箱壁，一块厚度为?，横截面积为A的平板，两侧表面分别维持着均匀温度tw1和tw2。实验表明：单位时间内从表面1 传递到表面2的热量Q与壁面间的温差

?t?tw1?tw2以及导热面积A成正比，而与平板的厚度?成反比，即

?方程中的比例系数?反映了材料导热能力的大小，叫做导热系数或者热导率，单位是W(m?K)。导热率是一个物性参数，不同材料的导热率不同，即使是同一种材料，在不同温度下也具有不同的热导率。

为了确保恒温箱内的光谱仪能够安全可靠的工作,我们需要将恒温箱内多余的热量排放到外面，从而保证恒温箱内的光谱仪能够获得适宜的工作温度，这是系统级热设计。

总的来说,恒温箱的热设计可以分为3个主要方面。第一：组件级热设计，即对元器件级的热设计。第二：封装级设计，即对于电子控制模块、散热设备、PCB板的热设计。第三：系统级设计，即对于恒温箱内胆及外壳的热设计。系统级的热设计主要研究电子设备所处环境的温度对其产生的影响，环境温度是系统级热分析的重要边界条件，

6

Q=?A?t

其热设计是采取措施控制环境温度，使电子设备在适宜的温度环境下进行工作。随着电子技术的发展，尤其是微电子技术的发展，电子元器件和设备的尺寸正迅速缩小，而功率却一直在增大，使得单位体积容纳的热量越来越多，特别是在恶劣环境下电子产品以及国防电子仪器要求必须达到良好的电磁兼容性能及三防性能等。

（2）三次样条差值

将环境温度每隔15min采集到的数据利用三次样条差值成每隔1min的数据。 算法原理：

设在区间[a,b]上给定n?1个节点xi(a?x0?x1???xn?b)，在节点xi处的函数值 为yi?f(xi)(i?0,1,?,n)。若函数S(x)满足如下三条：

（1） 在每个子区间[xi?1,xi](i?1,2,?,n)上，S(x)是三次多项式； （2） S(xi)?yi(i?1,2,?,n)；

（3） 在区间[a,b]上，S(x)的二阶导数S''(x)连续。

则称S(x)为函数y?f(x)在区间[a,b]上的三次样条插值函数。 子区间[xi?1,xi](i?1,2,?,n)上的S(x)的表达式为：

S(x)?(x3i?x)M(x?xi?1)3yh2ixi?x6hi?1?Mi?(i?1?6Mi?1)i6hihi2

?(yhix?xi?1i?6Mi)h,xi?1?x?xii关于参数Mi的方程组（三弯矩方程组）：

?iMi?1?2Mi??iMi?1?di,i?1,2,?,n?1??hihi?1ih,?i?hi?xi?xi?1i?hi?1h?1??i,i?hi?1d6i?h(yi?1?yi?yi?yi?1)?6f[xi?1,xi,xi?1]i?hi?1hi?1hi

牛顿插值多项式：

Nn(x)?f[x0]?f[x0,x1](x?x0)?f[x0,x1](x?x0)(x?x1)???f[x0,x1,?,xk](x?x0)(x?x1)?(x?xk?1)??

?f[x0,x1,?,xn](x?x0)(x?x1)?(x?xn?1).构造牛顿插值多项式首先列出差商表，进而由差商表写出牛顿插值多项式。n阶差商为：

f[xx0,x1,?,xn]?f[x0,x1,?,xn?1]0,x1,?,xn]?f[x

n?x0

零阶差商和一阶差商：

f[xi]?yi?f(xi)f[xi,xi?1]?f[xi?1]?f[xi]xi?1?xi

（3）龙格库塔法

算法原理：

考虑用函数f(x,y)在若干点上的函数值的线性组合来构造近似公式，构造时要求近似公式在(xi,yi)处的Taylor展开式与解y(x)在xi处的Taylor展开式的前面几项重合，从而使近似公式达到所需要的阶数。既避免求高阶导数,又提高了计算方法精度的阶数。或者说,在[xi,xi+1]这一步内多计算几个点的斜率值，然后将其进行加权平均作为平均斜率，则可构造出更高精度的计算格式，这就是龙格—库塔（Runge-Kutta）法的基本思想。用类似上述的处理方法，只需在区间[xi,xi+1]上用四个点处的斜率加权平均作为平均斜率K*的近似值，构成一系列四阶龙格—库塔公式。具有四阶精度，即局部截断误差是O(h5)。

对于高阶常微分方程的初值问题：

(m)'''(m?1)?),a?x?b?y?f(x,y,y,y,?,y ?''(m?1)(m?1)(a)?y0,??y(a)?y0,y(a)?y0,?,y将其转化为一阶微分方程组求解：

?y1'?y2?'?y2?y3???? ?'y?ym?m?1?y'?f(x,y,y,?,y)12m?m'(m?1)??y1(a)?y0,y2(a)?y0,?,ym(a)?y0标准的四级四阶龙格—库塔法的向量形式：

1?y?y?(K1?K2?K3?K4)i?i?16??K1?hf(xi,yi)?11? K?hf(x?h,y?K1)?2ii22?11?K?hf(x?h,y?K2)ii?322???K4?hf(xi?h,yi?K3)其分量形式：

1?y?y?(Kj1?Kj2?Kj3?Kj4)i,j?j,i?16??Kj1?hfj(xi;y1i,y2i,?,yni,)?Kn1K11K21h?K?hf(x?;y?,y?,?,y?)j?1,2,?,n ?j2ji1i2ini2222?Kn2K12K22h?K?hf(x?;y?,y?,?,y?)ji1i2ini?j32222?K?hfj(xi?h;y1i?K13,y2i?K23,?,yni?Kn3)??j4

（4）超定方程的最小二乘解

小二乘法广泛地应用于工程计算中，用最小二乘法消除（平滑）误差，用最小二乘法从有噪声的数据中提取信号，从海量数据中找出数据变化的趋势，……。甚至利用简单函数计算复杂函数的近似值，我们并不期望它的近似值多么精确（事实上很多时候也不用很精确），尽管如此还是希望计算出的近似数据与原始数据之间有相似之处。如果从线性代数角度来理解最小二乘法，实际上是将一个高维空间的向量投影到低维子空间所涉及的工作。

超定方程组的最小二乘解

当方程组GX=b的方程数多于未知数个数时，对应的系数矩阵G的行数大于列数，此时方程组被称为是超定方程组。设G=(giu)m×n，当m>n时即所谓的高矩阵，绝大多数情况下，超定方程组没有古典意义下的解。超定方程组的最小二乘解是一种广义解，是指使残差r = b – GX 的2-范数达取极小值的解，即

||b?GX*||2?min||b?GX||2 mX?R超定方程组的最小二乘解是指正规方程组

GTGX=GTb

的解。如果系数矩阵(GTG)可逆，则正规方程组有唯一解。此时，最小二乘解可以形式地写为如下形式

X=(GTG)-1GTb

两种常用的方法如下

1．用对称矩阵的三角分解法解正规方程组GTGX = GTb；记A=GTG，则A是对称矩阵，由三角分解A = L D LT，其中L是下三角矩阵，D是对角矩阵。将这一算法写过三个过程： ①解下三角方程组：LY1 = GTb； ②解对角方程组：DY2 = Y1； ③解上三角方程组：LTY3 = Y2

2．用矩阵的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华氏度后，箱顶隔热层未达到67.8华氏度，热源停止加热；恒温箱体给箱顶隔热层加热，恒温箱及箱顶隔热层向外界空气散热，如图

箱顶隔热层恒温箱体热源图4模型2恒温箱热量传播方向图

模型3

恒温箱及箱顶隔热层温度均高于67.8华氏度，热源停止加热；恒温箱温度低于箱顶隔热层，箱顶隔热层给恒温箱加热；恒温箱及箱顶隔热层向外界空气散热，如图

箱顶隔热层恒温箱体热源图5模型3恒温箱热量传播方向图

模型4

恒温箱体温度低于67.8华氏度后，此时箱顶隔热层温度高于67.8华氏度，热源加热，热源与箱顶隔热层给恒温箱体加热，直到箱顶隔热层与恒温箱体温度相等。

箱顶隔热层恒温箱体热源图6模型4恒温箱热量传播方向图

模型建模

根据已测的环境温度、恒温箱体和箱顶隔热层温度数据，在对每隔15min测量一次的环境温度数据进行三次样条差值后得到每隔1min一个数据，如图7所示。初始时刻，恒温箱体与隔热层温度等于环境温度。此时热源给恒温箱体加热，恒温箱体给箱顶隔热层加热；恒温箱及箱顶隔热层向外界空气散热。这段过程，恒温箱未达到68.2华氏度，恒温箱体温度高于箱顶隔热层。

根据能量守恒，恒温箱体增加的热量与恒温箱对环境空气和箱顶隔热层的散热之和

等于热源提供的热量。

c1m1?Tbox??1(Tbox?Tair)A1?t/?1??1(Tbox?Tgere)A2?t/?1?q?t（1）

恒温箱体对箱顶隔热层的散热等于箱顶隔热层增加的热量与箱顶隔热层对环境空气的散热之和。

c2m2?Tgere??1(Tbox?Tgere)A2?t/?1??2(Tgere?Tair)A2?t/?2（2）

式中各变量含义分别为c1是恒温箱比热容，m1是恒温箱质量，?Tbox是?t时间内恒温箱温度增量，?Tbox是?t时间内恒温箱温度增量，?1是恒温箱箱体热导率，Tbox是恒温箱某一时刻温度，Tair是环境温度，A1是恒温箱与外界环境的接触表面积，?1是恒温箱体的厚度，Tgere是隔热层某时刻温度，A2是恒温箱体和隔热层的接触表面积，q是热源?t时间内所释放的热量，c2是隔热层的比热容，m2是隔热层的质量，?Tgere是?t时间隔热层温度的增量，?2是隔热层的热导率，?2是隔热层的厚度。

将上式进一步简化为

A?t?B?Tbox?C(Tbox?Tair)?t?(Tbox?Tgere)?t（3） D?Tgere?E(Tgere?Tair)?t?(Tbox?Tgere)?t（4）

式中A?q?1cm?Acm???，B?111，C?1，D?221，E?21；?t定为1min ?1A2?1A2?1A2A2?1?2取差值后的80组数据中的前78组数据，求解超定方程中的A、B、C、D、E五个

参数。

A=22.4495 B=63.9245 C=0.0366 D=122.5048 E=0.4669

将能量守恒方程变换为功率守恒方程：

?TBbox?A?(Tbox?Tgere)?C(Tbox?Tair)（5）

?tD?Tgere?t?(Tbox?Tgere)?E(Tgere?Tair)（6）

代入参数值可得：

dT63.9245box?22.4495?1.0366Tbox?Tgere?0.0366Tair（7）

dt

122.5048dTgeredt?Tbox?1.4669Tgere?0.4669Tair（8）

将第79和80组数据带入上面的方程式，绝对误差分别为0.5624，0.3112和0.602，0.331，在接受范围内，验证了方程式的正确性。

图7环境温度实际测量与三次采样

根据式（7）和（8）反复迭代计算恒温箱与箱顶隔热层的温度，计算前80min的数据，可以得到恒温箱体与箱顶隔热层温度变化曲线，如图8和图9所示。

图8前80min恒温箱温度拟合曲线

图9前80min箱顶隔热层温度拟合曲线

由图可知，由建模得到的恒温箱体与箱顶隔热层温度变化，与温度实际测量曲线拟合得很好。

再根据模型2、3、4讨论的情况采用比较恒温箱、箱顶隔热层温度与临界温度的高低，再选择模型，具体程序见附录。由此可以得到恒温箱体和箱顶隔热层温度的变化曲线。仿真结果见图10和图11。

图10恒温箱与箱顶隔热层温度变化曲线

图11恒温箱温度变化曲线

由建模结果可知，箱顶隔热层始终低于68华氏度，箱顶隔热层与环境温度变化趋势相近。恒温箱体在刚开始持续受热，约在260分钟加热到68华氏度，然后在67.8华氏度与68.2华氏度附近来回震荡。对比图10和图11可知，随着环境温度的升高，恒温箱体温度震荡频率降低。环境温度决定了恒温箱温度的震荡频率，即热源开关的频率。

恒温箱箱体温度趋近于稳定后（300分钟~1440分钟），加热时间Theat约为317分钟，不加热时间约为814分钟,加热时间与不加热时间比值为317/814=0.4017。

恒温箱箱体温度趋近于稳定后（300分钟~1440分钟），最高温度68.4457华氏度，最低温度67.6212华氏度。

模型的改进恒温箱温度控制算法设计

建模前提是热源提供的速率一致，可以采用PID控制调节热源的加热速率。使恒温箱加热频率降低，响应速率变快。PID 控制算法简单、鲁棒性能好、可靠性高，被广泛应用于工业过程并取得了良好的控制效果。下图是PID 控制的原理框图。

将给定值r?t?与实际输出值c?t?相减，得到控制偏差：

e?t??r?t??c?t?

将偏差经过比例、积分、微分作用后线性相加，得到PID控制器输出的控制量，输入到被控对象进行控制：

?1u?t??Kp?e?t??TI?其传递函数形式为：

?t0e?t?dt?TDde?t??? dt???1Gc?s??Kp?1??Tps?

?TIs?式中，Kp是比例系数，TI是积分时间常数，TD是微分时间常数。

比例环节：即时成比例地反应控制系统的偏差信号e?t?，偏差一旦产生，控制器立即动作以减小误差。比例作用大，可以加快调节，减少误差。（比例也不能过大，过大

的比例会使系统的稳定性下降，甚至造成系统的不稳定。）

积分环节：主要作用是消除静差，提高系统的控制精度。积分作用的强弱取决于积分时间常数TI，其值越大，积分作用越弱；反之越强，引入积分调节可使系统稳定性下降，动态响应变慢。

微分环节：能反应偏差信号的变化趋势（变化速率），并能在偏差信号变得太大之前引入一个修正信号，从而减小系统的超调，加快系统的过渡过程，减小调节时间。但微分作用对噪声干扰有放大作用。

参考文献

[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.04.012. [4] 傅秦生. 热工基础与应用[M].北京：机械工业出版社，2007. [5] 李乃成，梅立泉. 数值分析[M].北京：科学出版社，2011.

附录

三次样条插值代码：

[hwx,grc]=textread('hwxwd.txt','%f%f','headerlines',1) wmwd=textread('wmwd.txt','%f') figure

plot(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_interp(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;

whilei<=n

while (i<=n)&&(y(1,i)<68.2)

K1(:,1)=[h*func1(y(1,i),y(2,i),Tair(i)),h*func3(y(1,i),y(2,i),Tair(i))]';

K2(:,1)=[h*func1(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,Ta

ir(i))]';

K3(:,1)=[h*func1(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*func1(y(1,i)+K3(1,1),y(2,i)+K3(2,1),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

while (i<=n)&&(y(1,i)>67.8)

K1(:,1)=[h*func2(y(1,i),y(2,i),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),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 end

hwxwd=y(1,:); grcwd=y(2,:); end

微分方程1：

functiondelta_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：

functiondelta_grc= func3( Tbox,Tgere,Tair ) %隔热层温度微分方程

delta_grc=1/122.5048*(Tbox-1.4669*Tgere+0.4669*Tair); end

本文来源：https://www.bwwdw.com/article/f5jw.html