西安交通大学——温度场数值模拟(matlab)

通过在室内的某些位置布置适当的节点，采集回来室内的温湿度以及空气质量等实际参数。

首先对室内空间建模，用一个无限细化的三维矩阵来模拟出室内的温度分布情况，针对采集回来的数据，采用插值法和适当次数的拟合函数的拟合，得出三维矩阵的实际值的分布，最后结合matlab软件绘制出计算出的温度场的三维图像。

一．数据的采集与处理因为影响人的舒适感的温度层只是室内的某一高度范围内的温度，而温度传感器虽然是布置在一个平面内，但是采用插值法和拟合函数法是可以大致再现出影响人的舒适感的温度层的温度变化的。

同时，在构建出的三维模型中，用第三维表示传感器层面的温度。

在传感器层面，传感器分布矩阵如下：X=【7.5 36.5 65.5】（模型内单位为cm）Y=【5.5 32.5 59.5】Z=【z1 z2 z3；z4 z5 z6；z7 z8 z9；】（传感器采集到的实时参数）采用meshgrid（xi，yi，zi，…）产生网格矩阵；首先按照人的最小温度分辨值，将室内的分布矩阵按照同样的比例细化，均分，使取值点在坐标一定程度上也是接近于连续变化的，从而才能最大程度上使处理数据得来的分布值按最小分辨值连续变化！根据人体散热量计算公式：C=hc（tb-Ta）其中hc为对流交换系数；结合Gagge教授提出的TSENS热感觉指标可以计算出不同环境下人的对环境温度变化时人体温度感知分辨率，作为插值法的一个参考量，能使绘制出的温度场更加的符合人体的温度变化模式。

例如按照10cm的均差产生网格矩阵（实际上人对温度的分辨率是远远10cm大于这个值的，但是那样产生的网格矩阵也是异常庞大的，例如以0.5cm为例，那么就可以获得116*108=12528个元素，为方便说明现已10cm为例）：[xi yi]=meshgrid(7.5:10:65.5,5.5:10:59.5)xi =7.5000 17.5000 27.5000 37.5000 47.5000 57.50007.5000 17.5000 27.5000 37.5000 47.5000 57.50007.5000 17.5000 27.5000 37.5000 47.5000 57.50007.5000 17.5000 27.5000 37.5000 47.5000 57.50007.5000 17.5000 27.5000 37.5000 47.5000 57.50007.5000 17.5000 27.5000 37.5000 47.5000 57.5000yi =5.5000 5.5000 5.5000 5.5000 5.5000 5.500015.5000 15.5000 15.5000 15.5000 15.5000 15.500025.5000 25.5000 25.5000 25.5000 25.5000 25.500035.5000 35.5000 35.5000 35.5000 35.5000 35.500045.5000 45.5000 45.5000 45.5000 45.5000 45.500055.5000 55.5000 55.5000 55.5000 55.5000 55.5000产生网格矩阵之后，就可以在测得的实时数据的基础上，通过相关的温度场的专业的估算函数，以及相关的数值处理函数来估计整个分布面（有最小的分辨率）上的温度了。

温度场模拟matlab 代码：clear,clc,clfL1=8;L2=8;N=9;M=9;% 边长为8cm 的正方形划分为8*8 的格子T0=500;Tw=100; % 初始和稳态温度a=0.05; % 导温系数tmax=600;dt=0.2; % 时间限10min 和时间步长0.2s dx=L1/(M-1);dy=L2/(N-1); M仁a*dt/(dx A2);M2=a*dt/(dy A2);T=T0*ones(M,N); T1=T0*ones(M,N);t=0;l=0;k=0;Tc=zeros(1,600);% 中心点温度，每一秒采集一个点for i=1:9 for j=1:9if(i==1|i==9|j==1|j==9)T(i,j)=Tw;% 边界点温度为100CelseT(i,j)=T0;endendendif(2*M1+2*M2<=1) % 判断是否满足稳定性条件while(t<tmax+dt)t=t+dt;k=k+1;for i=2:8for j=2:8 T1(i,j)=M1*(T(i-1,j)+T(i+1,j))+M2*(T(i,j-1)+T(i,j+1))+(1-2*M1-2*M2)*T(i,j);endendfor i=2:8for j=2:8T(i,j)=T1(i,j);endendif(k==5)l=l+1;Tc(l)=T(5,5);k=0;endendi=1:9;j=1:9;[x,y]=meshgrid(i);figure(1);subplot(1,2,1);mesh(x,y,T(i,j))% 画出10min 后的温度场axis tight;xlabel('x','FontSize',14);ylabel('y','FontSize',14);zlabel('T/ C ','FontSize',14)title('1min 后二维温度场模拟图','FontSize',18)subplot(1,2,2);[C,H]=co ntour(x,y,T(i,j));clabel(C,H);axis square;xlabel('x','Fo ntSize',14);ylabel('y','Fo ntSize',14);title('1min 后模拟等温线图','FontSize',18)figure(2);xx=1:600;plot(xx,Tc,'k-','li newidth',2)xlabel('时间/s','FontSize',14);ylabel('温度/ C ','FontSize',14);title('中心点的冷却曲线','Fo ntSize',18)else disp('Error!') % 如果不满足稳定性条件，显示" Error!”end实验结果:中心点的冷却曲线1min 后二维温度场模拟图222X100.06 100.05 100.04 100.03 100.02 100.01 100 9 r'5min 后二维温度场模拟图1061051041min 后模拟等温线图1OCOT-00110310min 后二维温度场模拟图100.06 _100.05 .100.04 .100.03 .y1102881876y 54321 123456789X06a1 o 1031Oo 5 100.'100. 04100.0103^25min 后模拟等温线图9876y 54321123456789X—101031-I10510104屛102101 —10min 后模拟等温线图9876y 54321 123456789X100.0 1I n ^n v u l -I I 40Q U 15U .O 0O- 2 1O 0403O0010-2U.D 0 10a100.0 2 1 00 .01 - -10min后二维温度场模拟图（不满足稳定性条件）中心、点的冷却曲线（不满足稳定性条件）。

热处理过程中温度场的数值模拟及分析热处理是一种常用的金属加工工艺，通过控制金属材料的加热与冷却过程，可以改变金属材料的组织结构和性能。

温度场是热处理过程中重要的参数之一，直接影响着金属材料的组织和性能的形成与变化。

因此，准确地模拟和分析热处理过程中的温度场对于优化工艺、改善产品质量具有重要意义。

数值模拟是研究温度场的有效方法之一。

它基于数学模型和计算方法，通过计算机的数值计算来获得温度场的分布情况。

在热处理过程中，温度场的分布受到多个因素的影响，如加热功率、材料热导率、热辐射、对流散热等。

数值模拟通过建立数学模型，考虑这些因素，并进行相应的计算，可以得到较为准确的温度场分布。

首先，进行数值模拟需要选择适当的数学模型。

在热处理过程中，常用的模型有热传导方程、能量方程等。

热传导方程是研究物体内部温度分布的基本方程，它考虑了热传导过程中的温度梯度对热流的影响。

能量方程则是考虑了热源与物体之间的热交换过程，可以更全面地描述温度场的变化。

其次，进行数值模拟需要确定边界条件。

边界条件是指在模拟过程中与外界接触的部分，它对于温度场的分布起着重要的影响。

常见的边界条件有热流、热辐射和对流散热等。

热流边界条件是指物体表面受到的外部热量输入或输出，热辐射边界条件是指物体表面受到的辐射热量，而对流散热边界条件则是指物体与周围介质间的热交换。

然后，进行数值模拟需要进行网格剖分。

网格剖分是将模拟区域分成小的单元，用于离散方程和计算。

在温度场的数值模拟中，常用的网格剖分方法有结构化网格和非结构化网格。

结构化网格是指将模拟区域划分为规则的矩形或立方体单元，易于计算和分析。

非结构化网格则是将模拟区域划分为任意形状的单元，适用于复杂几何形状和不均匀材料性质的模拟。

最后，进行数值模拟需要选择合适的求解方法。

在热处理过程中，常用的求解方法有有限差分法、有限元法和边界元法等。

有限差分法是基于差分逼近的一种方法，将参与方程离散化成代数方程，并通过迭代计算得到数值解。

二维导热物体温度场的数值模拟一、物理问题有一个用砖砌成的长方形截面的冷空气通道，其截面尺寸如下图1-1所示，假设在垂直于纸面方向上用冷空气及砖墙的温度变化很小，可以近似地予以忽略。

在下列两种情况下试计算： 砖墙横截面上的温度分布；垂直于纸面方向的每米长度上通过砖墙的导热量。

第一种情况：内外壁分别均匀维持在0℃及30℃； 第二种情况：内外壁均为第三类边界条件，且已知：Km K m W h C t Km W h C t ∙=∙=︒=∙=︒=∞∞/35.0/93.3,10/35.10,30222211λ砖墙导热系数二、数学描写由对称的界面必是绝热面，可取左上方的四分之一墙角为研究对象，该问题为二维、稳态、无内热源的导热问题。

控制方程：02222=∂∂+∂∂y tx t边界条件： 第一种情况：由对称性知边界1绝热： 0=w q ； 边界2为等温边界，满足第一类边界条件： C t w ︒=0；1-1图2-1图边界3为等温边界，满足第一类边界条件： C t w ︒=30。

第一种情况：由对称性知边界1绝热： 0=w q ；边界2为对流边界，满足第三类边界条件： )()(2f w w w t t h n tq -=∂∂-=λ； 边界3为对流边界，满足第三类边界条件： )()(2f w w w t t h ntq -=∂∂-=λ。

三、方程离散用一系列与坐标轴平行的间隔0.1m 的二维网格线将温度区域划分为若干子区域，如图1-3所示。

采用热平衡法，利用傅里叶导热定律和能量守恒定律，按照以导入元体（m,n ）方向的热流量为正，列写每个节点代表的元体的代数方程，第一种情况： 边界点：边界1（绝热边界）：5~2)2(411,11,12,1,m =++=+-m t t t t m m m ， 11~8)2(411,161,16,15,16=++=+-n t t t t n n n n，3-1图边界2（等温内边界）： 7,16~7;7~1,6,0,=====n m n m t nm边界3（等温外边界）： 12,16~2;12~1,1,30,=====n m n m t n m内节点：11~8,15~6;11~2,5~2)(411,1,,1,1,====+++=-+-+n m n m t t t t t n m n m n m n m n m第二种情况 边界点：边界1（绝热边界）： 5~2)2(411,11,12,1,m =++=+-m t t t t m m m ， 11~8)2(411,161,16,15,16=++=+-n t t t t n n n n ，边界2（内对流边界）：6~1)2(222111,61,6,5,6=++++=∆∆-+n Bi t Bi t t t t n n n n ，16~7)2(2221117,17,18,7,=++++=∆∆-+m Bi t Bi t t t t m m m m ，边界3（外对流边界）：11~1)2(2222221,11,1,2,1=++++=∆∆-+n Bi t Bi t t t t n n n n，16~2)2(22222212,112,111,12,=++++=∆∆-+m Bi t Bi t t t t m m m m ，内角点： )3(22)(21116,67,78,67,57,6+++++=∆∆Bi t Bi t t t t t外角点：)1(222211,112,212,1+++=∆∆Bi t Bi t t t内节点：11~8,15~6;11~2,5~2);(411,1,,1,1,====+++=-+-+n m n m t t t t t n m n m n m n m n m（10,22121==∆=∞∆t t xh Bi λ；30,21212==∆=∞∆t t xh Bi λ）四、编程思路及流程图编程思路为设定两个二维数组t(i,j)、ta(i,j)分别表示本次迭代和上次迭代各节点的温度值，iter （实际编程时并未按照此名称来命名迭代步长）表示迭代进行的次数, 1Q 、2Q 分别表示外边界、内边界的散热量。

大学物理仿真实验---热敏电阻温度特性曲线实验实验名称：热敏电阻温度特性曲线实验一．实验简介：热敏电阻是由对温度非常敏感的半导体陶瓷质工作体构成的元件。

与一般常用的金属电阻相比，它有大得多的电阻温度系数值。

热敏电阻作为温度传感器具有用料省、成本低、体积小等优点，可以简便灵敏地测量微小温度的变化，在很多科学研究领域都有广泛的应用。

二．实验目的：了解热敏电阻的电阻—温度特性及测温原理，学习惠斯通电桥的原理及使用方法，学习坐标变换、曲线改直的技巧。

三．实验原理：半导体热敏电阻的电阻—温度特性热敏电阻的电阻值与温度的关系为：A，B是与半导体材料有关的常数，T为绝对温度，根据定义，电阻温度系数为：R t是在温度为t时的电阻值。

惠斯通电桥的工作原理如图所示：四个电阻R0，R1，R2，Rx组成一个四边形，即电桥的四个臂，其中Rx就是待测电阻。

在四边形的一对对角A和C之间连接电源，而在另一对对角B和D之间接入检流计G。

当B和D两点电位相等时，G中无电流通过，电桥便达到了平衡。

平衡时必有Rx = (R1/R2)·R0，(R1/R2)和R0都已知，Rx即可求出。

电桥灵敏度的定义为：式中ΔRx指的是在电桥平衡后Rx的微小改变量，Δｎ越大，说明电桥灵敏度越高。

实验仪器四．实验装置：直流单臂电桥、检流计、待测热敏电阻和温度计、调压器。

五．实验内容：从室温开始，每隔5°C测量一次Rt，直到85°C。

撤去电炉，使水慢慢冷却，测量降温过程中，各对应温度点的Rt。

求升温和降温时的各Rt的平均值，然后绘制出热敏电阻的Rt-t特性曲线。

求出t＝50°C点的电阻温度系数。

作ln Rt ~ (1 / T)曲线，时)。

确定式(1)中常数A和B,再由（2）式求α (50°Ｃ六．实验所测数据：?不同T所对应的Rt 值????R t均值，1 / T，及ln R t的值七．数据处理：1.热敏电阻的R t-t特性曲线数据点连线作图在图上找到T=50所对应的点做切线，可以求得切线的斜率：K=（500-0）/(0-85)=5.88 由由此计算出：α＝-0.031二次拟合的曲线：在图上找到T=50所对应的点做切线，可以求得切线的斜率：K=（495-0）/(0-84)=5.89由由此计算出：α＝--0.030.2.ln R t -- (1 / T)曲线仿真实验画出图线如下图所示但计算机仿真实验画出的曲线图中A的值计算有误，正确的A=0.0153.将图修正后如下：A=0.0153，B=3047.5383由此写出R t= 0.0153由此当T=50时，α=-0.030八、思考题1. 如何提高电桥的灵敏度？答：电桥的灵敏度和电源电压，检流计的灵敏度成正比，因此提高电源电压，检流计的灵敏度能提高电桥灵敏度。

实验报告（二）完成人：L.W.Yohann注：本次实验主要学习了用MATLAB绘制二维、三维图形的基本命令、图形的标识与修饰以及用符号函数绘图，在学习完成后小组对52页的上机练习题进行了程序编辑和运行。

1.绘制数列变化趋势图.解：在编辑窗口输入：n=1:100;an=(1+1./n).^n;plot(n,an,'r*')grid并保存，命名为lab1；在命令窗口中输入lab1,得：2.绘制数列变化趋势图.解：在编辑窗口输入：n=1:0.1:50;an=n.^(1./n);plot(n,an,'r*')grid并保存，命名为lab2；在命令窗口中输入lab2,得：3.绘制函数在无定义点处的变化趋势.解：在编辑窗口输入：x=-10:0.05:10;y=sin(x)./x;plot(x,y,'r*')grid并保存，命名为lab3；在命令窗口中输入lab3,得：4.在同一坐标系中画出函数及其Taylor多项式的图像解：y=sinx在编辑窗口输入：syms xf=sin(x);T6=taylor(f,x);T8=taylor(f,x,'Order',8);T10=taylor(f,x,'Order',10);T12=taylor(f,x,'Order',12);fplot([T6 T8 T10 T12 f])xlim([-8 8])grid onlegend('approximation of sin(x) up to O(x^6)',...'approximation of sin(x) up to O(x^8)',...'approximation of sin(x) up to O(x^{10})',...'approximation of sin(x) up to O(x^{12})',...'sin(x)','Location','Best')title('Taylor Series Expansion')并保存，命名为lab4sin；在命令窗口中输入lab4sin,得：y=exp(x)在编辑窗口输入：syms xf=exp(x);T6=taylor(f,x);T8=taylor(f,x,'Order',8);T10=taylor(f,x,'Order',10);T12=taylor(f,x,'Order',12);fplot([T6 T8 T10 T12 f])xlim([-8 8])grid onlegend('approximation of exp(x) up to o(x^6)',...'approximation of exp(x) up to o(x^8)',...'approximation of exp(x) up to o(x^{10})',...'approximation of exp(x) up to o(x^{12})',...'exp(x)','Location','Best')title('Taylor Series Expansion')并保存，命名为lab4exp；在命令窗口中输入lab4exp,得：5.符号函数绘图.注：在matlab r2010b 和matlab r2019b中对绘制函数图像的输入方法有不同的要求，故此类题分两个版本来求解。

通过在室内的某些位置布置适当的节点，采集回来室内的温湿度以及空气质量等实际参数。

首先对室内空间建模，用一个无限细化的三维矩阵来模拟出室内的温度分布情况，针对采集回来的数据，采用插值法和适当次数的拟合函数的拟合，得出三维矩阵的实际值的分布，最后结合matlab软件绘制出计算出的温度场的三维图像。

一．数据的采集与处理因为影响人的舒适感的温度层只是室内的某一高度范围内的温度，而温度传感器虽然是布置在一个平面内，但是采用插值法和拟合函数法是可以大致再现出影响人的舒适感的温度层的温度变化的。

同时，在构建出的三维模型中，用第三维表示传感器层面的温度。

在传感器层面，传感器分布矩阵如下：X=【7.5 36.5 65.5】（模型内单位为cm）Y=【5.5 32.5 59.5】Z=【z1 z2 z3；z4 z5 z6；z7 z8 z9；】（传感器采集到的实时参数）采用meshgrid（xi，yi，zi，…）产生网格矩阵；首先按照人的最小温度分辨值，将室内的分布矩阵按照同样的比例细化，均分，使取值点在坐标一定程度上也是接近于连续变化的，从而才能最大程度上使处理数据得来的分布值按最小分辨值连续变化！根据人体散热量计算公式：C=hc（tb-Ta）其中hc为对流交换系数；结合Gagge教授提出的TSENS热感觉指标可以计算出不同环境下人的对环境温度变化时人体温度感知分辨率，作为插值法的一个参考量，能使绘制出的温度场更加的符合人体的温度变化模式。

例如按照10cm的均差产生网格矩阵（实际上人对温度的分辨率是远远10cm大于这个值的，但是那样产生的网格矩阵也是异常庞大的，例如以0.5cm为例，那么就可以获得116*108=12528个元素，为方便说明现已10cm为例）：[xi yi]=meshgrid(7.5:10:65.5,5.5:10:59.5)xi =7.5000 17.5000 27.5000 37.5000 47.5000 57.50007.5000 17.5000 27.5000 37.5000 47.5000 57.50007.5000 17.5000 27.5000 37.5000 47.5000 57.50007.5000 17.5000 27.5000 37.5000 47.5000 57.50007.5000 17.5000 27.5000 37.5000 47.5000 57.50007.5000 17.5000 27.5000 37.5000 47.5000 57.5000yi =5.5000 5.5000 5.5000 5.5000 5.5000 5.500015.5000 15.5000 15.5000 15.5000 15.5000 15.500025.5000 25.5000 25.5000 25.5000 25.5000 25.500035.5000 35.5000 35.5000 35.5000 35.5000 35.500045.5000 45.5000 45.5000 45.5000 45.5000 45.500055.5000 55.5000 55.5000 55.5000 55.5000 55.5000产生网格矩阵之后，就可以在测得的实时数据的基础上，通过相关的温度场的专业的估算函数，以及相关的数值处理函数来估计整个分布面（有最小的分辨率）上的温度了。

通过在室内的某些位置布置适当的节点，采集回来室内的温湿度以及空气质量等实际参数。

首先对室内空间建模，用一个无限细化的三维矩阵来模拟出室内的温度分布情况，针对采集回来的数据，采用插值法和适当次数的拟合函数的拟合，得出三维矩阵的实际值的分布，最后结合matlab软件绘制出计算出的温度场的三维图像。

一．数据的采集与处理因为影响人的舒适感的温度层只是室内的某一高度范围内的温度，而温度传感器虽然是布置在一个平面内，但是采用插值法和拟合函数法是可以大致再现出影响人的舒适感的温度层的温度变化的。

同时，在构建出的三维模型中，用第三维表示传感器层面的温度。

在传感器层面，传感器分布矩阵如下：X=【7.5 36.5 65.5】（模型内单位为cm）Y=【5.5 32.5 59.5】Z=【z1 z2 z3；z4 z5 z6；z7 z8 z9；】（传感器采集到的实时参数）采用meshgrid（xi，yi，zi，…）产生网格矩阵；首先按照人的最小温度分辨值，将室内的分布矩阵按照同样的比例细化，均分，使取值点在坐标一定程度上也是接近于连续变化的，从而才能最大程度上使处理数据得来的分布值按最小分辨值连续变化！根据人体散热量计算公式：C=hc（tb-Ta）其中hc为对流交换系数；结合Gagge教授提出的TSENS热感觉指标可以计算出不同环境下人的对环境温度变化时人体温度感知分辨率，作为插值法的一个参考量，能使绘制出的温度场更加的符合人体的温度变化模式。

例如按照10cm的均差产生网格矩阵（实际上人对温度的分辨率是远远10cm大于这个值的，但是那样产生的网格矩阵也是异常庞大的，例如以0.5cm为例，那么就可以获得116*108=12528个元素，为方便说明现已10cm为例）：[xi yi]=meshgrid(7.5:10:65.5,5.5:10:59.5)xi =7.5000 17.5000 27.5000 37.5000 47.5000 57.50007.5000 17.5000 27.5000 37.5000 47.5000 57.50007.5000 17.5000 27.5000 37.5000 47.5000 57.50007.5000 17.5000 27.5000 37.5000 47.5000 57.50007.5000 17.5000 27.5000 37.5000 47.5000 57.50007.5000 17.5000 27.5000 37.5000 47.5000 57.5000yi =5.5000 5.5000 5.5000 5.5000 5.5000 5.500015.5000 15.5000 15.5000 15.5000 15.5000 15.500025.5000 25.5000 25.5000 25.5000 25.5000 25.500035.5000 35.5000 35.5000 35.5000 35.5000 35.500045.5000 45.5000 45.5000 45.5000 45.5000 45.500055.5000 55.5000 55.5000 55.5000 55.5000 55.5000产生网格矩阵之后，就可以在测得的实时数据的基础上，通过相关的温度场的专业的估算函数，以及相关的数值处理函数来估计整个分布面（有最小的分辨率）上的温度了。

西安交⼤传热学上机实验报告传热学上机实验报告⼆维导热物体温度场的数值模拟学院：化⼯学院姓名：沈佳磊学号：2110307016班级：装备11⼀、物理问题有⼀个⽤砖砌成的长⽅形截⾯的冷空⽓空道，其截⾯尺⼨如下图所⽰，假设在垂直于纸⾯⽅向上冷空⽓及砖墙的温度变化很⼩，可以近似地予以忽略。

在下列两种情况下试计算：（1）砖墙横截⾯上的温度分布；（2）垂直于纸⾯⽅向的每⽶长度上通过砖墙的导热量。

外矩形长为3.0m，宽为2.2m；内矩形长为2.0m，宽为1.2m。

第⼀种情况：内外壁分别均匀地维持在0℃及30℃；第⼆种情况：内外表⾯均为第三类边界条件，且已知：外壁：30℃，h1=10W/m2·℃,内壁：10℃，h2= 4 W/m2·℃砖墙的导热系数λ=0.53 W/m·℃由于对称性，仅研究1/4部分即可。

⼆、数学描写对于⼆维稳态导热问题，描写物体温度分布的微分⽅程为拉普拉斯⽅程22220t t x x ??+=??这是描写实验情景的控制⽅程。

三、⽅程离散⽤⼀系列与坐标轴平⾏的⽹格线把求解区域划分成许多⼦区域，以⽹格线的交点作为确定温度值的空间位置，即节点。

每⼀个节点都可以看成是以它为中⼼的⼀个⼩区域的代表。

由于对称性，仅研究1/4部分即可。

依照实验时得点划分⽹格。

建⽴节点物理量的代数⽅程对于内部节点，由?x=?y ，有,1,1,,1,11()4m n m n m n m n m n t t t t t +-+-=+++由于本实验为恒壁温，不涉及对流，故内⾓点，边界点代数⽅程与该式相同。

设⽴迭代初场，求解代数⽅程组图中，除边界上各节点温度为已知且不变外，其余各节点均需建⽴类似3中的离散⽅程，构成⼀个封闭的代数⽅程组。

以t ? =0°C 为场的初始温度，代⼊⽅程组迭代，直⾄相邻两次内外传热值之差⼩于0.01，认为已达到迭代收敛。

四、编程及结果program mainimplicit nonereal ,dimension(1:16,1:12)::treal ,dimension(1:16,1:12)::t1real q,q1,q2,q3,q4,q5,q6,q7,q8,q9,q10,q11,a integer m,n,z logical::converged=.false.z=1t=0a=0.53do n=1,12t(1,n)=30end dodo m=2,16t(m,12)=30end dodo n=1,7t(6,n)=0end dodo m=7,16t(m,7)=0end dodo while(.not.converged.and.z<10000)t1=tdo m=2,5do n=1,11if( n==1 )thent(m,n)=0.25*(t(m-1,n)+t(m+1,n)+2*t(m,n+1))elset(m,n)=0.25*(t(m-1,n)+t(m+1,n)+t(m,n-1)+t(m,n+1)) end if end doend dodo n=8,11do m=6,16if (m==16) thent(m,n)=0.25*(t(m,n-1)+t(m,n+1)+2*t(m-1,n)) elset(m,n)=0.25*(t(m-1,n)+t(m+1,n)+t(m,n-1)+t(m,n+1)) end if end doend doz=z+1do m=1,16do n=1,12if(abs(t(m,n)-t1(m,n))>0.000001) thenconverged=.false.exitelseconverged=.true.end ifend doend doend dowrite(*,'(16f5.1)',advance='no')((t(m,n),m=1,16),n=12,7,-1) write(*,*) write(*,'(6f5.1)',advance='no')((t(m,n),m=1,6),n=6,1,-1)do n=2,11q1=(t(1,n)-t(2,n))*a+q1end dodo m=2,15q2=(t(m,12)-t(m,11))*a+q2end doq3=(t(1,1)-t(2,1))*a*0.5q4=(t(16,12)-t(16,11))*a*0.5q10=q1+q2+q3+q4write(*,*)do n=2,6q5=(t(5,n)-t(6,n))*a+q5end dodo m=7,15q6=(t(m,8)-t(m,7))*a+q6end doq7=(t(5,1)-t(6,1))*a*0.5q8=(t(16,8)-t(16,7))*a*0.5q9=(t(5,7)-t(6,7))*a*2q11=q5+q6+q7+q8+q9q=(q10+q11)*0.5*4print*,"外表⾯导量=",q10,"内表⾯导热量",q11,"每⽶⾼砖墙导热量",q end结果截图：将以上结果⽤matlab画图⼯具绘制出如下图像：。

实用文档传热大作业二维导热物体温度场的数值模拟：璇班级：能动A02学号：10031096一．物理问题有一个用砖砌成的长方形截面的冷空气通道，其截面尺寸如下图所示，假设在垂直于纸面方向上用冷空气及砖墙的温度变化很小，可以近似地予以忽略。

在下列两种情况下试计算：（1）砖墙横截面上的温度分布；（2）垂直于纸面方向的每米长度上通过砖墙的导热量。

第一种情况：外壁分别均与地维持在0℃及30℃；第二种情况：外壁均为第三类边界条件，且已知：t ∞1=30℃,ℎ1=10wm2∙℃t ∞2=10℃,ℎ2=4wm2∙℃砖墙的导热系数λ=0.53 Wm∙℃二．数学描写由对称的界面必是绝热面，可取左上方的四分之一墙角为研究对象，该问题为二维、稳态、无热源的导热问题，其控制方程和边界条件如下：ðt2ðx2+ðt2ðy2=0边界条件（情况一） t(x,0)=30 0≤x≤1.5t(0,y)=30 0≤y≤1.1t(0.5,y)=0 0.5≤y≤1.1t(x,0.5)=0 0.5≤x≤1.5ðt(1.5,y)=0 0≤y≤0.5ðy∂t(x,1.1)=0 0≤x≤0.5ℎ(t−t f1) x=0,0≤y≤1.11=ℎ(t−t f2) x=0.5,0.5≤y≤1.12ℎ(t−t f1) y=0,0≤x≤1.51ℎ(t−t f2) y=0,0.5≤x≤21.50 0≤y≤0.5=0 0≤x≤0.5∂x三．网格划分网格划分与传热学实验指导书中“二维导热物体温度场的电模拟实验”一致，如下图所示：四．方程离散对于节点，离散方程t[i][j]=0.25*(t[i+1][j]+t[i-1][j]+t[i][j+1]+t[i][j-1])对于边界节点，则应对一、二两种情况分开讨论：情况一:绝热平直边界点：t[15][j]=0.25*(2*t[14][j]+t[15][j-1]+t[15][j+1])1≤j≤4t[i][11]=0.25*(2*t[i][10]+t[i-1][11]+t[i+1][11]) 1≤i≤4外等温边界点：t[i][j]=30等温边界点：t[i][j]=0情况二：（Bi1,Bi2为网格Bi数，Bi1=ℎ1∆xλ Bi2=ℎ2∆xλ）绝热平直边界点：t[15][j]=0.25*(2*t[14][j]+t[15][j-1]+t[15][j+1])1≤j≤4t[i][11]=0.25*(2*t[i][10]+t[i-1][11]+t[i+1][11]) 1≤i≤4外侧对流平直边界：t[i][0]=(2*t[i][1]+t[i+1][0]+t[i-1][0]+2*Bi1*tf1)/(2*Bi1+4) 1≤i≤14t[0][j]=(2*t[1][j]+t[0][j+1]+t[0][j-1]+2*Bi1*tf1)/(2*Bi1+4) 1≤j≤10侧对流平直边界：t[i][5]=(2*t[i][4]+t[i+1][5]+t[i-1][5]+2*Bi2*tf2)/(2*Bi2+4) 6≤i≤14t[5][j]=(2*t[4][j]+t[5][j+1]+t[5][j-1]+2*Bi2*tf2)/(2*Bi2+4) 6≤j≤10特殊点：a点t[15][0]=(t[14][0]+t[15][1]+tf1*Bi1)/(Bi1+2)b点t[15][5]=(t[14][5]+t[15][4]+tf2*Bi2)/(Bi2+2)c点t[5][5]=(2*t[4][5]+2*t[5][4]+t[5][6]+t[6][5]+3*Bi2*tf2)/(2*Bi2+6) d点t[5][11]=(t[5][10]+t[4][11]+tf2*Bi2)/(Bi2+2)e点t[0][11]=(t[0][10]+t[1][11]+tf1*Bi1)/(Bi1+2)f点t[0][0]=(t[0][1]+t[1][0]+tf1*Bi1*2)/(2*Bi1+2)五．编程思路及流程图编程思路为设定两个二维数组t[i][j]、ta[i][j]分别表示本次迭代和上次迭代各节点的温度值，iter表示迭代进行的次数,daore_in、daore_out分别表示外边界的散热量。

通过在室内的某些位置布置适当的节点，采集回来室内的温湿度以及空气质量等实际参数。

首先对室内空间建模，用一个无限细化的三维矩阵来模拟出室内的温度分布情况，针对采集回来的数据，采用插值法和适当次数的拟合函数的拟合，得出三维矩阵的实际值的分布，最后结合matlab软件绘制出计算出的温度场的三维图像。

一.数据的采集与处理因为影响人的舒适感的温度层只是室内的某一高度范围内的温度，而温度传感器虽然是布置在一个平面内，但是采用插值法和拟合函数法是可以大致再现出影响人的舒适感的温度层的温度变化的。

同时，在构建出的三维模型中，用第三维表示传感器层面的温度。

在传感器层面，传感器分布矩阵如下：X=【7.5 36.5 65.5】(模型内单位为cm)Y=【5.5 32.5 59.5 】Z=【z1 z2 z3;z4 z5 z6;z7 z8 z9;】(传感器采集到的实时参数)采用meshgrid (xi, yi, zi,…)产生网格矩阵；首先按照人的最小温度分辨值，将室内的分布矩阵按照同样的比例细化，均分，使取值点在坐标一定程度上也是接近于连续变化的，从而才能最大程度上使处理数据得来的分布值按最小分辨值连续变化！沃於厘rr見公玖】C =根据人体散热量计算公式：C=hc (tb-Ta) 「‘ i L-；,*其中hc为对流交换系数；结合Gagge教授提出的TSENS热感觉指标可以计算出不同环境下人的对环境温度变化时人体温度感知分辨率，作为插值法的一个参考量，能使绘制出的温度场更加的符合人体的温度变化模式。

例如按照10cm的均差产生网格矩阵(实际上人对温度的分辨率是远远10cm大于这个值的，但是那样产生的网格矩阵也是异常庞大的，例如以0.5cm为例，那么就可以获得116*108=12528个元素，为方便说明现已10cm为例)：[xi yi]=meshgrid(7.5:10:65.5,5.5:10:59.5)xi =7.5000 17.5000 27.5000 37.5000 47.5000 57.50007.5000 17.5000 27.5000 37.5000 47.5000 57.50007.5000 17.5000 27.5000 37.5000 47.5000 57.50007.5000 17.5000 27.5000 37.5000 47.5000 57.50007.5000 17.5000 27.5000 37.5000 47.5000 57.50007.5000 17.5000 27.5000 37.5000 47.5000 57.5000yi =5.5000 5.5000 5.5000 5.5000 5.5000 5.500015.5000 15.5000 15.5000 15.5000 15.5000 15.500025.5000 25.5000 25.5000 25.5000 25.5000 25.500035.5000 35.5000 35.5000 35.5000 35.5000 35.500045.5000 45.5000 45.5000 45.5000 45.5000 45.500055.5000 55.5000 55.5000 55.5000 55.5000 55.5000产生网格矩阵之后，就可以在测得的实时数据的基础上，通过相关的温度场的专业的估算函数，以及相关的数值处理函数来估计整个分布面（有最小的分辨率）上的温度了。

基于Matlab导热问题的数值模拟徐凯;石利娜;吴东垠【摘要】基于Matlab软件,采用有限差分法、pdepe函数法和pdetool(工具箱)法对3种不同形式的导热问题进行数值模拟,具体包括二维稳态导热、一维非稳态导热和二维非稳态导热,得出温度分布的数值解和图形解.通过求解过程可以看出,数值模拟方法对分析导热物体的温度分布非常直观和方便.【期刊名称】《上海工程技术大学学报》【年(卷),期】2016(030)004【总页数】6页(P353-358)【关键词】导热;Matlab软件;温度;数值模拟【作者】徐凯;石利娜;吴东垠【作者单位】西安交通大学能源与动力工程学院,西安710049;榆林学院管理学院,榆林719000;西安交通大学能源与动力工程学院,西安710049【正文语种】中文【中图分类】TK11导热是传热学三大热量传递的途径之一,在工程实际中有着广泛的应用,许多工程问题都要求解物体内部温度的分布和非稳态传热中温度随时间的变化[1-2].导热问题大多为多维的偏微分方程,在已知的定解条件(初始条件和边界条件)下求解析解.对于定解条件和几何尺寸相对简单的导热问题,可以求其解析解,但是对于比较复杂的导热问题,求解析解就显得非常麻烦,有时候甚至难以求解[3-4].伴随计算机的发展,数值模拟方法为求解导热问题提供了新的思路,这种方法由于能够处理解析解不能处理的问题,并且准确度高,因而越来越受到广大科技工作者的重视,在处理导热问题中起着不可替代的作用[5-6].本文以Matlab软件为载体,运用不同的数值模拟方法,包括有限差分法、pdepe函数法和pdetool(工具箱)法对不同的导热过程(二维稳态导热、一维非稳态导热和二维非稳态导热)进行数值模拟[7],得出其温度分布的数值解和温度分布图.对于任意一个导热问题,不管求解析解还是数值解,数学描写都是不可或缺的.为了得到导热物体温度场的数学描写,根据能量守恒定律和傅里叶定律建立导热微分方程,得到笛卡尔坐标系下的三维非稳态导热微分方程为式中:ρ为密度,kg/m3；c为比热容，J/(kg·K)；t为温度，℃;τ为时间,s;λ为导热系数为内热源，W/m3.根据不同的情形简化式(1)可以得到不同的导热过程.为了求解式(1),需要添加定解条件,包括初始条件和边界条件.导热微分方程和定解条件构成导热问题完整的数学模型.2.1 有限差分法有限差分法,即将求解区域划分为差分网格,用有限个网格节点代替连续的求解区域.再将偏微分方程的偏导数用差商代替,推导含有离散点上有限个未知数的差分方程.求解差分方程,就是微分方程定解问题的近似数值解,如图1所示.一个无内热源的二维稳态导热物体,其上凹面、下表面分别维持在t1=0 ℃,t2=30 ℃,其余表面绝热,求其温度分布.物理方程和定解条件为t1=0,t2=30对几何结构进行离散化后,再对偏微分方程离散,写出差分方程,如图2和式(2)所示. 为方便起见,取步长Δx=Δy=1.用Matlab编程,分析图2.由于其不是矩形区域,所以对于网格要分区域进行处理,以m=1∶10和n=1∶21为例,其设计的计算流程图如图3所示.其余区域与此类似,不作重复.在Matlab软件下运行编制的程序,经过数次迭代后,产生数值解(由于篇幅限制,仅给出部分节点数值解)见表1.温度分布图如图4所示.图4(a)为等温边界温度分布图,图4(b)为等温边界温度分布云图.2.2 pdepe函数法Matlab语言提供了pdepe函数,可直接调用其标准形式.设一细棒由各项均匀材料组成,长度为2 m,热扩散率为0.005 m2/s,细棒两端绝热,初始时间细棒温度分布为T0=30+10[1-cos(πx)],求非稳态下细棒的温度分布.物理方程和定解条件为pdepe函数的标准形式为式中,m=0,1,2分别为平面、圆柱和球.调用形式为边界条件为p(x,t,u)+q(x,t,u).调用形式为[pl,ql,pr,qr]=bcpde(xl,ul,xr,xr,t)初始条件为调用形式为经过分析得至此,可以编程来求解此问题,具体程序在此省略.由于pdepe函数编程过程中可以调用已有的函数,所以主程序的条数相对于有限差分法简单了许多.在编程过程中,借调sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t)求解.在Matlab中运行程序.图5为非稳态下细棒的温度分布.调用pdeval(m,x,u,xout)画出不同时间的数值解图像,图6为不同时间下细棒的温度分布.从图6可以看出,初始时间温度分布呈正弦曲线,这是由初始时间给定的温度函数决定的.随着时间的推移,各节点温度大小分布的差异逐渐减小,当t>50 s时,温度分布将趋于稳定,在最后时间的温度分布基本不变,保持为40 ℃.2.3 pdetool法对于二维稳态和非稳态导热问题,偏微分方程,其解析解往往非常复杂,有时候甚至无法得到解析解,只能采用数值模拟的方法.但是编制程序时却很复杂.Matlab中提供的pdetool法不用编程就可以解出导热问题的数值解,它是基于有限元法解偏微分方程.但是pdetool法只能解决二维模型.对于一维的要扩成二维,三维的缩成二维,时间维不计算在内.pdetool法可以解出4种偏微分方程式(4)～(7)分别为椭圆形方程、抛物型方程、双曲型方程和特征值方程.式中：u为域Ω上的求解变量;λ为特征值;d、c、a、f为常数或变量；t为时间变量. pdetool提供了两类边界条件,即式中：n为垂直于边界的单位矢量；g、q、h、r为定义在边界上的函数.式(8)为Dirichlet条件,式(9)为Neumann条件,至此,可以应用pdetool法来处理具体问题.一直经为0.15 m、高为0.05 m的平板玻璃圆盘,送入退火炉中消除应力,其初始温度为30 ℃,炉中温度为450 ℃,设该玻璃盘在炉内时，各表面均可受到加热,表面传热系数为9.5 W/(m2·K),按工艺要求,须加热到盘内各处温度均为400 ℃以上,估计所需的加热时间.已知该盘的导热系数λ=0.78 W/(m·K),质量密度ρ=2 700 kg/m3,比定压热容cp=835 J/(kg·K).物理方程为这是一个三维问题,在pdetool中,需要将三维缩为二维.因此,取柱坐标(r,θ,z).由于其关于轴的对称性,故与θ无关,从而简化成仅关于(r,z)的二维方程,即与式(5)的抛物型方程对比后得出d=ρcr=2 700×835×0.075=169 087.5c=λr=0.78×0.075=0.058 5a=0,f=0对边界条件有此条件为Neumann条件,将其转化为标准形式,得到由c=λr,代入(12),有得出Neumann条件中q=hr=9.5×0.075=0.712 5g=450×hr=450×0.712 5=320.625至此,就可以应用pdetool法求解此问题,前期工作为二维图形的绘制、边界条件的输入、方程形式的选取,在绘制二维图时,选用以角点方式画矩形,边界条件选择Neumann条件,方程形式选择抛物型,输入上述所求的对应参数,并对求解区域画出网格,图7为计算网格.在solve选项中输入初始温度和不同的时间,得出不同时间的温度图,图8为玻璃圆盘在不同时间的温度分布图.由图8可以看出,大约经过9 500 s,即可达到盘内各处温度均为400 ℃以上,符合题设要求.当然,还可以解出不同时间各个节点处的温度数值解,由于篇幅有限,这里不再列出.本文基于Matlab软件,采用有限差分法、pdepe函数法和pdetool法,对3种不同形式的导热问题进行了数值模拟.研究表明，对于导热问题的数值解法,不仅求解过程简单,而且结果更加直观明显,随着计算机技术的发展,数值模拟方法还将会得到越来越广泛的应用.对于同一种问题,可以有不同的模拟方法,比如本文用有限差分法处理的二维稳态问题,同样也可以用pdetool法处理；用pdepe函数法处理的一维非稳态问题也可以用有限差分法来编程完成,只是编程过程略显繁琐.对于具体问题,应根据实际情况选出最合适的方法来处理.【相关文献】[1] 陶文铨.数值传热学[M].2版.西安:西安交通大学出版社,2001.[2] 安德森约翰D.计算流体力学基础及其应用[M].吴颂平，刘赵淼，译.北京:机械工业出版社,2007.[3] 王福军.计算流体动力学分析:CFD软件原理与分析[M].北京:清华大学出版社,2011.[4] 田禾.关于二维非稳态导热的可视化研究[D].天津:天津师范大学,2003.[5] 杨世铭,陶文铨.传热学[M].4版.北京:高等教育出版社,2006.[6] 李明.偏微分方程的MATLAB解法[J].湖南农机(学术版),2010,37(3):89-91.[7] 彭东玲,张义,方慧,等.日光温室墙体一维导热的MATLAB模拟与热流分析[J].中国农业大学学报,2014,19(5):174-179.。

温度场模拟matlab代码：clear,clc,clfL1=8;L2=8;N=9;M=9;% 边长为8cm的正方形划分为8*8的格子T0=500;Tw=100; % 初始和稳态温度a=0.05; % 导温系数tmax=600;dt=0.2; % 时间限10min和时间步长0.2sdx=L1/(M-1);dy=L2/(N-1);M1=a*dt/(dx^2);M2=a*dt/(dy^2);T=T0*ones(M,N);T1=T0*ones(M,N);t=0;l=0;k=0;Tc=zeros(1,600);% 中心点温度，每一秒采集一个点for i=1:9for j=1:9if(i==1|i==9|j==1|j==9)T(i,j)=Tw;% 边界点温度为100℃elseT(i,j)=T0;endendendif(2*M1+2*M2<=1) % 判断是否满足稳定性条件while(t<tmax+dt)t=t+dt;k=k+1;for i=2:8for j=2:8T1(i,j)=M1*(T(i-1,j)+T(i+1,j))+M2*(T(i,j-1)+T(i,j+1))+(1-2*M1-2*M2)*T(i,j);endendfor i=2:8for j=2:8T(i,j)=T1(i,j);endendif(k==5)l=l+1;Tc(l)=T(5,5);k=0;endendi=1:9;j=1:9;[x,y]=meshgrid(i); figure(1);subplot(1,2,1);mesh(x,y,T(i,j))% 画出10min 后的温度场 axis tight;xlabel('x','FontSize',14);ylabel('y','FontSize',14);zlabel('T/℃','FontSize',14) title('1min 后二维温度场模拟图','FontSize',18) subplot(1,2,2);[C,H]=contour(x,y,T(i,j)); clabel(C,H);axis square;xlabel('x','FontSize',14);ylabel('y','FontSize',14); title('1min 后模拟等温线图','FontSize',18) figure(2); xx=1:600;plot(xx,Tc,'k-','linewidth',2)xlabel('时间/s','FontSize',14);ylabel('温度/℃','FontSize',14);title('中心点的冷却曲线','FontSize',18)else disp('Error!') % 如果不满足稳定性条件，显示“Error ！” end实验结果：时间/s温度/℃中心点的冷却曲线x1min后二维温度场模拟图T /℃xy1min 后模拟等温线图x5min 后二维温度场模拟图T /℃xy5min 后模拟等温线图x10min后二维温度场模拟图T /℃xy10min 后模拟等温线图x10min 后二维温度场模拟图（不满足稳定性条件）yT /℃21时间/s温度/℃中心点的冷却曲线（不满足稳定性条件）。