仿真报告修改后.docx
《仿真报告修改后.docx》由会员分享，可在线阅读，更多相关《仿真报告修改后.docx（10页珍藏版）》请在冰豆网上搜索。
仿真报告修改后
基于matlab对管式反应器内传热过程的研究
1.热传导的基本原理
热传导是由物质内部分子、原子和自由电子等微观粒子的热运动而产生的热量传递现象。
热传导的机理非常复杂，简而言之，非金属固体内部的热传导是通过相邻分子在碰撞时传递振动能实现的；金属固体的导热主要通过自由电子的迁移传递热量；在流体特别是气体中，热传导则是由于分子不规则的热运动引起的。
1.1温度场和等温面
任一瞬间物体或系统内各点温度分布的空间，称为温度场。
在同一瞬间，具有相同温度的各点组成的面称为等温面。
因为空间内任一点不可能同时具有一个以上的不同温度，所以温度不同的等温面不能相交。
1.2温度梯度
从任一点开始，沿等温面移动，如图1所示，因为在等温面上无温度变化，所以无热量传递；而沿和等温面相交的任何方向移动，都有温度变化，在与等温面垂直的方向上温度变化率最大。
将相邻两等温面之间的温度差△t与两等温面之间的垂直距离△n之比的极限称为温度梯度，其数学定义式为：
1.3傅里叶定律
导热的机理相当复杂，但其宏观规律可以用傅里叶定律来描述，其数学表达式为：
图1温度梯度与傅里叶定律
傅里叶定律表明：
在热传导时，其传热速率与温度梯度及传热面积成正比。
2.模型描述
如图所示，一高温管式反应器由内径5cm，壁厚2.5cm的陶瓷材料制成，为实现高温条件下的反应，需要对反应器进行保温，保温层厚度为5cm，陶瓷和保温层的物性为：
陶瓷：
ρ1=2600kg/m3c1=1150J/（kg·K）k1=3.0W/（m·K）
保温层：
ρ2=600kg/m3c2=200J/（kg·K）k2=0.2W/（m·K）
反应器初始温度是300K，内表面突然被暴露在1500K的反应物流中，管壁和反应物流间的对流换热系数为500W/（m2·K）。
请给出保温层内表面温度随时间的变化规律，以及保温层内温度分布，讨论保温层厚度的影响。
（保温层外表面视为绝热）
图2管式反应器的结构
3.构建数学模型
该问题是柱坐标系下的一维非稳态导热问题，比较特殊之处在于该体系分陶瓷层、保温层两部分，共有内、中、外三个边界条件。
建立柱坐标系并传热微分方程进行化简，不难得到如下关系式：
其中，无量纲温度
（
=1500K,
=300K）
根据题目给出的参数可校验发现该圆柱是厚壁物体，则必须解以上偏微分方程，考虑用差分方程的方法给出数值解。
为保证结果的稳定性不受步长影响，选择用时间向后、空间一阶向后、空间二阶中心差分方程。
以上上标n表示时间,下标i表示半径。
Δt,Δr分别为时间、半径的步长。
由此对于给定时间点，可以得到一个大规模方程组，自变量为各个r位置的无量纲温度。
以上即可解决三个非边界点以外的所有差分点的列式，下面再给出边界点的列式。
（1）r=0.0025m处，陶瓷材料与热流界面处在该位置可通过热平衡法得到其表达式:
进而做差分近似处理得
对于无内热源的边界微元层，我们可以认为在足够小的微元层内无热量积累。
则有：
进而做差分近似处理得
这样，方程求解会变得简单一些。
（2）r=0.0050m处，陶瓷材料与保温层界面处，同理的可以由热量衡算给出其微分表达式。
进而做差分近似处理得
此处左边表达式用向前差分，右边为向后差分。
（3）r=0.010m处，保温层外界面处，壁面处绝热，故近似认为
，由此经整理可得到完整的差分方程组。
不妨选取Δr=0.0005m，Δt=1s,并用x代替Θ，则有：
方程
（1）：
方程
（2）~方程（50）:
方程（51）:
方程（52）~方程（150）:
方程（151）:
α1,α2分别为陶瓷层，保温层的导热系数。
借助Matlab对以上大矩阵做高斯-赛达尔迭代进行求解。
4.结果分析
4.1整体图像
取时间t=1:
2000s，空间半径r=0.0025:
0.0005:
0.0100，将无量纲温度转换为对温度T。
先从整体上观察温度场随时间、空间的变化规律。
从图上可得到信息：
1）随着时间增加，同一空间点（同一薄圆柱壳层）的温度逐渐增加，最初为300K，在t=2000s时各空间点均已经上升至1200K以上。
2）同一时间，随着半径的增大（远离热反应物流），温度逐渐降低。
以t=2000s为例，陶瓷层内壁1476.6K，保温层外壁1205.1K。
3）还可看出在两层导热介质交界处存在温度梯度突变；余位置温度梯度连续，曲线、曲面光滑。
4.2保温层内表面温度Tc分析
先做出Tc-t图像，再截取部分Tc-t观察其变化规律。
表1Tc-t变化数据
t/s
Tc/K
t/s
Tc/K
1
300.0000
51
308.9147
2
300.0050
52
309.5783
3
300.0097
53
310.2682
4
300.0143
54
310.9841
5
300.0190
55
311.7259
6
300.0236
56
312.4933
7
300.0283
57
313.2861
8
300.0330
58
314.1039
9
300.0378
59
314.9464
10
300.0427
60
315.8133
11
300.0479
61
316.7043
12
300.0535
62
317.6189
13
300.0597
63
318.5568
14
300.0669
64
319.5176
15
300.0754
65
320.5009
16
300.0858
66
321.5063
17
300.0985
67
321.5063
18
300.1144
68
322.5334
19
300.1342
69
324.6510
20
300.1589
70
325.7406
21
300.1893
71
326.8503
22
300.2266
72
327.9794
23
300.2718
73
329.1278
24
300.3261
74
330.2948
25
300.3908
75
331.4801
26
300.4669
76
332.6832
27
300.5558
77
333.9038
28
300.6586
78
335.1414
29
300.7765
79
336.3955
30
300.9107
80
337.6658
31
301.0623
81
338.9519
32
301.2324
82
340.2534
33
301.4219
83
341.5698
34
301.6318
84
342.9008
35
301.8630
85
344.2459
36
302.1163
86
345.6048
37
302.3924
87
346.9771
38
302.6921
88
348.3625
39
303.0159
89
349.7605
40
303.3644
90
351.1708
41
303.7380
91
352.5930
42
304.1371
92
354.0269
43
304.5621
93
355.4720
44
305.0133
94
356.9280
45
305.4908
95
358.3945
46
305.9949
96
359.8714
47
306.5255
97
361.3582
48
307.0829
98
362.8546
49
307.6669
99
364.3603
50
308.2775
100
365.8750
截取了前100s的数据进行分析，发现保温层内壁处温度变化还是比较缓慢的。
经过100s的时间温度仅从300K升至365K左右。
当t=2000s时该位置的温度上升至1423.9K,已经比较接近热流温度。
之所以说其温度变化缓慢，是因为陶瓷层内温度梯度很陡，以t=100s为例给出温度场图.看以下图5可发现,保温层内温度变化相比陶瓷层十分缓慢，明显这是因为后者有更小的导热系数k。
4.3保温层内温度场分析
保温层内温度随位置、时间的变化规律与之前整体图像规律一致。
但是与前面对比易知，在保温层内温度无论是随时间还是随空间的变化都比陶瓷层内平缓。
4.4保温层厚度的影响
定性的分析我们很容易知道：
1.随保温层厚度增加，保温层外表面温度会降低，单位表面积散热量会减小
2.但是保温层变厚后，保温层外表面积增大，故总的散热量是否减小则不一定
3.但基于该问题，题目认为保温层是绝热的，那么保温层的厚度如何改变都不会影响散热。
下面试取边界层厚度d=0.025m,d=0.060m（原先d=0.050m）
观察整体温度场图像以及保温层内表面的温度变化曲线、t=100s时整体温度场。
将结果与之前比较。
1）d=0.025m
图7d=0.025m整体温度场
显然可以看出此时保温层外壁的温度变高了，达到1400K以上。
与原图像对比，总体变化趋势不变，但温度上升变快了，当t=2000s时保温层内表面温度达1455.9K,比d=0.0050m时高了200K左右。
保温层内侧温度提高，其余基本一致。
2）d=0.060m>d0
对比之前两张图，当保温层变厚，可以看出保温效果变好，例如t=2000s时保温层外壁温度低至1100K以下。
与前两个d的情形相比，变化趋势不变；与d=0.05m相比，陶瓷层内温度降低不明显，例如保温层内壁温度此时为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;
M=2600*1150*（dr）^2/dt;
N=600*200*（dr）^2/dt;
fori=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）;
forj=2:
50
A1（j,j）=（dr^2）/（alpha1*dt）-dr/r（j）+2;
end
forj=1:
150
A2（j,j+1）=-1;
end
A1（1,1）=500+3/dr;
A2（1,2）=-3/dr;
A3=zeros（151,151）;
forj=1:
150
A3（j+1,j）=dr/r（j+1）-1;
end
forj=52:
150
A1（j,j）=（dr^2）/（alpha2*dt）-dr/r（j）+2;
end
A1（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=1;
x（:
:
1）=1;
forn=2:
2000
b（1,1,n）=0;%由前一时间循环结果给下一时间的方程组b赋值
forq=2:
50
b（q,1,n）=dr^2/（alpha1*dt）*x（q,j,n-1）;
end
b（51,1,n）=0;
forq=52:
150
b（q,1,n）=dr^2/（alpha2*dt）*x（q,j,n-1）;
end
b（151,1,n）=0;
forj=1:
149%设计循环，高斯-赛德尔迭代求解,符合要求时跳出循环,得到迭代次数j
ifj==1
x（:
j+1,n）=inv（D-L）*U*x（:
j,n）+inv（D-L）*b（:
:
n）;
elseifnorm（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'）
升级会员 