仿真报告修改后.docx

1、仿真报告修改后基于matlab对管式反应器内传热过程的研究1.热传导的基本原理 热传导是由物质内部分子、原子和自由电子等微观粒子的热运动而产生的热量传递现象。热传导的机理非常复杂，简而言之，非金属固体内部的热传导是通过相邻分子在碰撞时传递振动能实现的；金属固体的导热主要通过自由电子的迁移传递热量；在流体特别是气体中，热传导则是由于分子不规则的热运动引起的。1.1温度场和等温面任一瞬间物体或系统内各点温度分布的空间，称为温度场。在同一瞬间，具有相同温度的各点组成的面称为等温面。因为空间内任一点不可能同时具有一个以上的不同温度，所以温度不同的等温面不能相交。1.2温度梯度 从任一点开始，沿等温面移

2、动，如图1所示，因为在等温面上无温度变化，所以无热量传递；而沿和等温面相交的任何方向移动，都有温度变化，在与等温面垂直的方向上温度变化率最大。将相邻两等温面之间的温度差t与两等温面之间的垂直距离n之比的极限称为温度梯度，其数学定义式为：1.3傅里叶定律 导热的机理相当复杂，但其宏观规律可以用傅里叶定律来描述，其数学表达式为： 图1 温度梯度与傅里叶定律 傅里叶定律表明：在热传导时，其传热速率与温度梯度及传热面积成正比。2.模型描述 如图所示，一高温管式反应器由内径 5cm，壁厚 2.5cm 的陶瓷材料制成，为实现高温条件下的反应，需要对反应器进行保温，保温层厚度为 5cm，陶瓷和保温层的物性为

3、： 陶瓷：1=2600 kg/m3 c1=1150 J/(kgK ) k1=3.0 W/(mK) 保温层：2=600 kg/m3 c2=200 J/(kgK ) k2=0.2 W/(mK) 反应器初始温度是 300 K，内表面突然被暴露在 1500K 的反应物流中，管壁和反应物流间的对流换热系数为 500 W/(m2K)。请给出保温层内表面温度随时间的变化规律，以及保温层内温度分布，讨论保温层厚度的影响。（保温层外表面视为绝热） 图2 管式反应器的结构 3.构建数学模型 该问题是柱坐标系下的一维非稳态导热问题， 比较特殊之处在于该体系分陶瓷层、 保温层两部分，共有内、中、外三个边界条件。 建立

4、柱坐标系并传热微分方程进行化简，不难得到如下关系式： 其中，无量纲温度（=1500K,=300K） 根据题目给出的参数可校验发现该圆柱是厚壁物体， 则必须解以上偏微分方程， 考虑用差分方程的方法给出数值解。 为保证结果的稳定性不受步长影响，选择用时间向后、空间一阶向后、空间二阶中心差分方程。 以上上标n表示时间,下标i表示半径。t,r分别为时间、半径的步长。 由此对于给定时间点， 可以得到一个大规模方程组， 自变量为各个 r 位置的无量纲温度。 以上即可解决三个非边界点以外的所有差分点的列式，下面再给出边界点的列式。 (1) r=0.0025m 处，陶瓷材料与热流界面处在该位置可通过热平衡法得

5、到其表达式: 进而做差分近似处理得 对于无内热源的边界微元层，我们可以认为在足够小的微元层内无热量积累。则有： 进而做差分近似处理得 这样，方程求解会变得简单一些。 (2) r=0.0050m 处，陶瓷材料与保温层界面处，同理的可以由热量衡算给出其微分表达式。 进而做差分近似处理得 此处左边表达式用向前差分，右边为向后差分。 (3) r=0.010m 处，保温层外界面处，壁面处绝热，故近似认为，由此经整理可得到完整的差分方程组。 不妨选取r=0.0005m，t=1s,并用x代替，则有： 方程(1)： 方程(2)方程(50): 方程(51): 方程(52)方程(150): 方程(151): 1,

6、2分别为陶瓷层，保温层的导热系数。 借助 Matlab 对以上大矩阵做高斯-赛达尔迭代进行求解。4.结果分析4.1 整体图像 图3 整体温度场变化图像 取时间 t=1:2000s ，空间半径 r=0.0025:0.0005:0.0100，将无量纲温度转换为对温度 T。先从整体上观察温度场随时间、空间的变化规律。从图上可得到信息： 1）随着时间增加，同一空间点（同一薄圆柱壳层）的温度逐渐增加， 最初为 300K，在 t=2000s时各空间点均已经上升至 1200K 以上。 2）同一时间，随着半径的增大（远离热反应物流） ，温度逐渐降低。以 t=2000s 为例，陶瓷层内壁 1476.6K，保温层

7、外壁 1205.1K。 3）还可看出在两层导热介质交界处存在温度梯度突变；余位置温度梯度连续，曲线、曲面光滑。 4.2 保温层内表面温度 Tc 分析 图4 Tc-t关系图像 先做出 Tc-t 图像，再截取部分 Tc-t 观察其变化规律。 表1 Tc-t 变化数据 t/sTc/Kt/sTc/K1300.000051308.91472300.005052309.57833300.009753310.26824300.014354310.98415300.019055311.72596300.023656312.49337300.028357313.28618300.033058314.103993

8、00.037859314.946410300.042760315.813311300.047961316.704312300.053562317.618913300.059763318.556814300.066964319.517615300.075465320.500916300.085866321.506317300.098567321.506318300.114468322.533419300.134269324.651020300.158970325.740621300.189371326.850322300.226672327.979423300.271873329.1278243

9、00.326174330.294825300.390875331.480126300.466976332.683227300.555877333.903828300.658678335.141429300.776579336.395530300.910780337.665831301.062381338.951932301.232482340.253433301.421983341.569834301.631884342.900835301.863085344.245936302.116386345.604837302.392487346.977138302.692188348.3625393

10、03.015989349.760540303.364490351.170841303.738091352.593042304.137192354.026943304.562193355.472044305.013394356.928045305.490895358.394546305.994996359.871447306.525597361.358248307.082998362.854649307.666999364.360350308.2775100365.8750 截取了前100s的数据进行分析，发现保温层内壁处温度变化还是比较缓慢的。经过100s的时间温度仅从300K升至365K左右

11、。 当t=2000s时 该位置的温度上升至1423.9K,已经比较接近热流温度。 之所以说其温度变化缓慢，是因为陶瓷层内温度梯度很陡，以t=100s为例给出温度场图.看以下图5可发现,保温层内温度变化相比陶瓷层十分缓慢，明显这是因为后者有更小的导热系数k。 图5 t=100s时整体温度场4.3 保温层内温度场分析 图6 保温层内温度场 保温层内温度随位置、时间的变化规律与之前整体图像规律一致。但是与前面对比易知，在保温层内温度无论是随时间还是随空间的变化都比陶瓷层内平缓。4.4 保温层厚度的影响 定性的分析我们很容易知道： 1.随保温层厚度增加，保温层外表面温度会降低，单位表面积散热量会减小

12、2.但是保温层变厚后，保温层外表面积增大，故总的散热量是否减小则不一定 3.但基于该问题， 题目认为保温层是绝热的， 那么保温层的厚度如何改变都不会影响散热。 下面试取边界层厚度 d=0.025 m, d=0.060 m（原先 d=0.050 m） 观察整体温度场图像 以及保温层内表面的温度变化曲线、 t=100s 时整体温度场。 将结果与之前比较。 1）d=0.025 m d0 对比之前两张图，当保温层变厚，可以看出保温效果变好，例如 t=2000s 时保温层外壁温度低至1100K以下。 图12 d=0.060m 保温层内壁温度变化 与前两个d的情形相比，变化趋势不变；与d=0.05m相比，

13、陶瓷层内温度降低不明显，例如保温层内壁温度此时为365.875K，与d=0.05m的数据持平。但是因为其保温层变厚，其外壁温度更低，在t=2000s时外壁仅1035K。 对比发现 这条曲线与d=0.05m时的曲线基本一致，t=2000s时保温层内壁温度也与之基本持平。 可见此时增加保温层厚度对陶瓷层的温度场影响不大，但是使保温层的温度场“延长”了，以致外壁温度更低。故可得出结论： 1. 原题情形下，t=2000s 时 保温层内壁 1476.6K，外壁 1205.1K 2. 随着保温层厚度增加，保温层外壁温度降低，但保温层内壁温度未必受明显影响 5.matlab程序dr=0.0005; dt=1

14、; M=2600*1150*( dr)2/ dt; N=600*200*( dr)2/ dt; for i=1:151 r(i)=0.025+(i-1)* dr; end T=ones(151,2500); x=zeros(151,150,2500); A1=eye(151,151); A2=zeros(151,151); alpha1=3/(2600*1150); alpha2=0.2/(600*200); for j=2:50 A1(j,j)=(dr2)/(alpha1* dt)- dr/r(j)+2; end for j=1:150 A2(j,j+1)=-1; end A1(1,1)=5

15、00+3/dr; A2(1,2)=-3/dr; A3=zeros(151,151); for j=1:150 A3(j+1,j)= dr /r(j+1)-1; end for j=52:150 A1(j,j)=( dr 2)/(alpha2* dt)- dr /r(j)+2; endA1(51,51)=-(3+0.2); A2(51,52)=0.2; A3(51,50)=3; A=A1+A2+A3;A(151,151)=1; A(151,150)=-1;b=zeros(151,1,2500);L=-1*tril(A,-1); U=-1*triu(A,1); D=diag(diag(A); j=

16、1; x(:,:,1)=1; for n=2:2000 b(1,1,n)=0; %由前一时间循环结果给下一时间的方程组b赋值 for q=2:50 b(q,1,n)= dr 2/(alpha1* dt) *x(q,j,n-1); end b(51,1,n)=0; for q=52:150 b(q,1,n)= dr2/(alpha2* dt) *x(q,j,n-1); end b(151,1,n)=0;for j=1:149 %设计循环，高斯-赛德尔迭代求解,符合要求时跳出循环,得到迭代次数 j if j=1 x(:,j+1,n)=inv(D-L)*U*x(:,j,n)+inv(D-L)*b(:,:,n); elseif norm(x(:,j,n)-x(:,j-1,n),inf) 10-5 x(:,j+1,n)=inv(D-L)*U*x(:,j,n)+inv(D-L)*b(:,:,n); else T(:,n)=x(:,j,n); break end end end T2=1500-1200*T; t=1:2000; X,Y=meshgrid(t,r); T2=T2(:,1:2000); mesh(X,Y,T2) xlabel(t/s) ylabel(r/m) zlabel(T/K)plot(r,T2(:,100) xlabel(r/m ) ylabel(T/K )