仿真报告修改后Word格式.docx
《仿真报告修改后Word格式.docx》由会员分享，可在线阅读，更多相关《仿真报告修改后Word格式.docx（10页珍藏版）》请在冰豆网上搜索。
（
=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
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
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<
d0
图7d=0.025m整体温度场
显然可以看出此时保温层外壁的温度变高了，达到1400K以上。
与原图像对比，总体变化趋势不变，但温度上升变快了，当t=2000s时保温层内表面温度达1455.9K,比d=0.0050m时高了200K左右。
保温层内侧温度提高，其余基本一致。
2）d=0.060m>
对比之前两张图，当保温层变厚，可以看出保温效果变好，例如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;
forj=1:
150
A2（j,j+1）=-1;
A1（1,1）=500+3/dr;
A2（1,2）=-3/dr;
A3=zeros（151,151）;
A3（j+1,j）=dr/r（j+1）-1;
forj=52:
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:
b（q,1,n）=dr^2/（alpha1*dt）*x（q,j,n-1）;
b（51,1,n）=0;
forq=52:
b（q,1,n）=dr^2/（alpha2*dt）*x（q,j,n-1）;
b（151,1,n）=0;
149%设计循环，高斯-赛德尔迭代求解,符合要求时跳出循环,得到迭代次数j
ifj==1
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（:
else
T（:
n）=x（:
j,n）;
break
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））
r/m'
T/K'
升级会员 