工程热力学课程设计(论文)模板_V2.0.docx

《工程热力学》课程设计（论文）

题 目

学号

姓名

院系

专业

完成日期

授课教师

得分

要求：
1. 论文内容（任选其中一类）：
a）新能源/可再生能源开发与利用现状综述（论文类）；
b）国防技术（航空、航天、航海）领域中的新兴动力技术（论文类）；
c）热力学研究进展（论文类）；
d）节能减排装置/系统设计（设计类）。
2. 不能从网上下载，抄袭；
3. 13周周五（本学期最后一次课）交论文，需同时提交打印版（需装订）和电子版（Word文档，Word2000以上版本）各1份。提交电子版时，文件题目格式：学号_姓名，如“001_张三”，由各班负责人收齐后提交。不按要求提交论文按无效提交处理。

设计类参考格式
1．总体要求
全文控制在8页A4纸以内，设计作品类应包括下列内容：作品背景、设计制作中解决的关键技术问题的描述、作品或模型的图片、创新特色、预计应用前景等、参考文献等。采用word 2000及以上版本编排。
2．页面要求
A4页面。页边距：上25mm，下25mm，左、右各20mm。正文采用小四号字体，标准字间距，单倍行间距。不要设置页眉，页码位于页面底部居中。
3．图表要求
插图按序编号，并加图名（位于图下方），采用嵌入型版式。图中文字用小五号宋体，符号用小五号Times New Roman（矢量、矩阵用黑斜体）；坐标图的横纵坐标应标注对应量的名称和符号/单位。
表格按序编号，并加表题（位于表上方）。采用三线表，必要时可加辅助线。
4．字号、字体要求（仅作参考）
家用电器节能系统设计说明书
学号：×××；班号：×××；设计者：×××
（空一行）
作品内容简介

通过实验设计了一套家用电器节能系统……（400—600字以内）
（空一行）
1 研制背景及意义

2 设计方案
2.1 电器控制
电器部分采用自动控制比较容易实现，考虑到电器元件易发热等问题……
……
2.2 机械部分
机械部分设计如图1所示，……
设计时考虑的主要问题：
……
3 理论设计计算
……
4 工作原理及性能分析
……
完成制作后，作品实物外形照片见图9。
5 创新点及应用
1）适用于不同类型家庭电器。
2）操作和控制简便，容易地使用它。
3）……。

在全国大中城市，家用电器普及数量很多，所有电器都有待在节能措施实施改进，因此应用前景很广。

图1 家用节能机构原理图
……
正文中表示物理量的符号，表示点、线、面的字母均用Times New Roman斜体；
表示法定计量单位、词头的符号、函数等，化学元素符号均用Times New Roman正体。

（空一行）
参考文献

[1] xxx，xxx．家用电器节能现状和发展．节能机械，2001，23(3)：275-279
[2] xxx．节能技术基础．xxxxxx出版社，1996：15-47
[3] xxxxx，xxxx．xxx，xxx译．机器人操作的数学导论．xxxx出版社，1998：11-67
[4] Lee H Y, Reinholtz C F. Inverse kinematics of serial-chain manipulators[J]. ASME Journal of Mechanical Design. 1996, 118(3): 396-404

论文类参考格式
1．总体要求
不少于3000字，全文控制在8页A4纸以内。论文应包括下列内容：摘要、关键词、引言、综述、结论、参考文献等。采用word 2000及以上版本编排。
2．页面要求
A4页面。页边距：上25mm，下25mm，左、右各20mm。正文采用小四号字体，标准字间距，单倍行间距。不要设置页眉，页码位于页面底部居中。
3．图表要求
插图按序编号，并加图名（位于图下方），采用嵌入型版式。图中文字用小五号宋体，符号用小五号Times New Roman（矢量、矩阵用黑斜体）；坐标图的横纵坐标应标注对应量的名称和符号/单位。
表格按序编号，并加表题（位于表上方）。采用三线表，必要时可加辅助线。
4．论文结构（仅作参考）

考虑温度反馈的点堆中子动力学求解
学号：×××；班号：×××；作者：×××
摘要：本文分别用采用绝热模型的幂级数法和Simulink动态仿真的方法，对考虑简单温度反馈的点堆中子动力学方程的进行求解。得到了考虑简单温度反馈时，反应性、反应堆功率、缓发中子先驱核浓度，反应堆温度随时间的变化情况，并通过对引入小阶跃反应性情况下计算结果的分析得出了绝热模型的适用范围。
关键词：温度反馈、幂级数法、Simulink仿真、绝热模型
1 引 言
在核反应堆实际的运行过程中，反应性是不断变化的，影响其变化的因素有很多，在中子动力学分析中注重的是反应堆在启动、停堆和功率调节过程中短期内中子密度随着时间的变化关系，因此可以不考虑燃料的燃耗和裂变产物的积累对反应性的影响，只需要考虑温度对反应性的影响。考虑温度反馈可以对反应堆功率的瞬态行为和其它反应堆堆芯系统变量进行估计，这些变量都是紧密耦合的，考虑温度反馈的点堆中子动力学方程为非线性方程组，求解比较困难。以前尝试的求解方法多数只适用于瞬发临界。目前对于考虑温度反馈且引入大阶跃反应性（反应堆处于超瞬发临界）的情况，主要采用诺德黑姆-福赫斯模型[1]进行处理，可得到较好的结果，但其只针对考虑单组缓发中子效应的问题。比较常见的反应性反馈的数学模型还有准稳态模型和绝热模型[2]，其中绝热模型因其形式简单，能大大降低方程求解的难度，得到了广泛地应用。但由于其忽略了反应堆的热量损失，在一些情况下并不适用。
本文在求解考虑温度反馈的点堆中子动力学方程时就采用了绝热模型，并应用求解点堆中子动力学方程综合性能较好的幂级数法[3]进行编程求解，并对所得结果进行分析。此外，还在不应用简化模型的情况下，对考虑温度反馈的点堆中子动力学方程中的各变化量进行Simulink动态仿真，得到了反应性、反应堆功率、缓发中子先驱核浓度，反应堆温度随时间的变化情况。
2 采用绝热模型的幂级数法
考虑温度对反应性的反馈作用时，反应性ρ不但是时间的函数，还是中子通量密度或中子密度n(t)的函数，因此考虑6组缓发中子效应的点堆中子动力学方程变为非线性方程组。不考虑外加中子源时，方程组具有如下形式：

(1)

(2)

(3)

(4)
式中，T(t)为核反应堆温度，K；T(0)为T(t)的初始值，K；Tc为堆芯平均冷却剂温度（不随时间变化），K；1/γ可被理解为热能由燃料传递给冷却剂的平均时间，s；ρ0为初始时引入的阶跃反应性；α为反应性温度系数α<0，K-1。
采用绝热模型对方程组进行简化，认为热量传递的时间常数1γ与功率偏离的时间相比很大，可忽略热量损失，即γ=0。则有

(5)
由(4)和(5)式可得：

(6)
将方程(1)、(2)、(6)改写为如下矩阵形式：

(7)
式中，y(t)为8维列向量，M(t)为(8)×(8)矩阵：

利用幂级数作为基函数对(7)式在一个时间步长内展开为：

(8)
式中，A0、A1、A2⋯均为8维列向量。
由初始条件，根据(8)式可得

(9)
另由(39)式可得

(10)
式中，

t=0时，由方程(9)和(10)可得

(11)
因此，要想求得Ai，只需求解由(11)式分别取n=1,2,⋯,k组成一个线性方程组即可（n为整数，是幂级数的阶数），再将Ai 代入(4-22)式即可求得该步长内点堆中子动力学方程组的解，然后将所得末点值作为下一个时间步长的初始值重新计算，如此迭代循环可以得到整个时间段内的解。
当迭代计算进行第一步时，均为零矩阵，则(11)式可化简为

(12)
进行后续迭代计算时，则有[4]：

且(i≥2)均为零矩阵，则(11)式可化简为：

(13)
其中，。

计算时采用参考文献[5]中给出的反应堆特性参数及初始条件。平均中子代时间Λ=5×10-5s，缓发中子总份额β=0.0065，6组缓发中子份额及先驱核衰变常数如表1所示。α=-0.00306β，K=7.6886×10-5K/J，初始条件：ρ0=0.043β，n0=0.01W。
表1 六组缓发中子份额及先驱核衰变常数
组数i
1
2
3
4
5
6
份额βi(×10-3)
0.22
1.42
1.27
2.57
0.75
0.27
衰变常数λi(s-1)
0.0124
0.0305
0.111
0.301
1.14
3.01
将计算得到的100min内方程组的解随时间变化的部分结果列于表2中。参考文献[5]中对考虑简单温度反馈的点堆中子动力学方程求解时没有采用简化模型，其计算结果还是比较准确的，将其作为标准参考值与自己用幂级数法采用绝热模型编程求得的结果进行比较。参考文献[5]中给出的100min内方程组的解随时间变化的部分结果见表3。
由表2和表3可以看出，在取四位有效数字的情况下，前30min内自己用幂级数法编程求得的结果与参考文献[5]中给出的结果基本完全相同，而30min后的计算结果与文献[5]中的结果开始产生明显的差距。可见绝热模型在阶跃反应性引入的初始一段时间还是非常适用的。
为了更清楚地了解采用绝热模型的幂级数法的计算精度，求取其计算结果与参考文献[5]中给出结果的相对误差，得到相对误差随时间变化的曲线如图1至图3所示。
表2 幂级数法计算结果
t(min)
ρ($)
n(W)
C1
C2
C3
C4
C5
C6
0
0.043
0.01
3.548
9.311
2.288
1.708
0.1316
0.01794
10
0.04299
0.1033
28.23
85.78
22.88
17.43
1.355
0.1851
20
0.04294
0.9514
260.1
790.1
210.7
160.5
12.48
1.705
30
0.04245
8.652
2371
7196
1917
1460
113.5
15.50
40
0.03820
70.42
1.975×104
5.926×104
1.566×104
1.190×104
924.0
126.2
50
0.01426
278.8
8.930×104
2.503×105
6.322×104
4.745×104
3665
500.0
60
-0.024
194.3
7.876×104
1.919×105
4.521×104
3.339×104
2561
348.9
70
-0.0391
46.24
2.084×104
4.723×104
1.085×104
7968
609.8
83.02
80
-0.0423
8.505
3925
8744
1998
1466
112.2
15.27
90
-0.0429
1.493
692.0
1537
350.8
257.4
19.70
2.681
100
-0.043
0.2600
120.6
267.6
61.09
44.83
3.429
0.4668
表3 参考文献[5]中给出的计算结果
t(min)
ρ($)
n(W)
C1
C2
C3
C4
C5
C6
0
0.043
0.01
3.548
9.311
2.288
1.708
0.1316
0.01794
10
0.04299
0.1033
28.23
85.78
22.88
17.43
1.355
0.1851
20
0.04294
0.9515
260.1
790.2
210.7
160.5
12.48
1.705
30
0.04245
8.654
2372
7197
1917
1460
113.5
15.51
40
0.03825
70.51
1.976×104
5.933×104
1.568×104
1.191×104
925.2
126.4
50
0.01468
282.2
9.018×104
2.531×105
6.396×104
4.802×104
3709
506.0
60
-0.0224
205.2
8.243×104
2.019×105
4.768×104
3.524×104
2704
368.3
70
-0.0355
53.99
2.382×104
5.469×104
1.264×104
9295
711.9
96.93
80
-0.0356
12.09
5362
1.225×104
2830
2082
159.5
21.71
90
-0.0328
2.917
1272
2936
681.5
501.9
38.46
6.237
100
-0.0295
0.7945
339.2
793.7
185.2
136.6
10.47
1.426

图1 反应性的相对误差随时间变化

图2 功率的相对误差随时间变化

图3 缓发中子先驱核浓度的相对误差随时间变化
由图1至图3可以看出由采用绝热模型的幂级数法计算得到前40min之内方程的结果与参考文献[5]中的结果间的相对误差非常小，但之后随着时间的增长相对误差明显增大，100min时反应性的相对误差已达到约0.46，反应堆功率和缓发中子先驱核浓度的相对误差已达到约0.68。
分析误差产生的主要原因应该是计算时采用了绝热模型这一简化处理，由于初始时引入了一个较小的阶跃反应性，中子的增殖或功率的增长要在较长时间内完成，在开始的较短时间内可以忽略反应堆的热量损失，这时绝热模型是适用的，因此在采用绝热模型计算的前一段时间的结果还是非常准确的，但随着时间的增长，反应堆有足够的时间将热量传递出去，此时采用绝热模型就会使求得的温度随时间变化比实际情况下快，从而使反应性的变化也快于实际情况。
3 基于Simulink的动态响应仿真
创建动态系统Simulink模型的一般步骤为：首先建立理论数学模型，列出描写系统动力学的全部方程。然后，根据理论数学模型，在Simulink模块库中选择所需的模块并用鼠标勾画各模块间的信号连线，建成仿真系统。最后，根据所研究的具体问题设置各模块的参数和仿真时间并进行仿真和调试得到满足需要的Simulink模型。
由于参考文献[5]中给出的考虑温度反馈的方程并未采用任何化简模型，而是较充分地考虑了反应堆的热量传递等情况，因此，采用参考文献[5]中所列出的考虑温度反馈的点堆中子动力学方程组作为理论数学模型。忽略外加中子源项，并考虑引入阶跃反应性的情况，其形式如下：

(14)

(15)

(16)

(17)
其中

式中，H、D分别为反应堆的高度、直径，m；T为反应堆的温度，K；CP为反应堆的热容，J/K；ε为裂变能转化为热能的分数；Tc为冷却剂的温度，K；Kc为冷却剂传热系数，J∙m-1∙s-1∙K-1；ρc为冷却剂密度，kg/m3；μ为冷却剂的质量流速，kg∙m-1∙s-1；
为冷却剂的比热容，J∙kg-1∙K-1。
根据理论数学模型搭建完成的Simulink系统仿真模型如图4所示。

图4 Simulink系统仿真模型

采用参考文献[5]中给出的各物理量的值对各模块的非结构参数进行设置。其中反应堆特性参数及初始条件已在上一节中列出，其它各量的值为：ε=1，A=17.52，H=0.23m，D=0.20m，Cp=1800πD2/4∙H∙1000，Tc=20℃，α=-0.00306β。
在用Simulink对点堆中子动力学方程组进行求解时可能会涉及到刚性问题，为保持求解的稳定性，选择专门用于处理刚性问题的变步长解法ode45[6,7]。设置仿真时间为6000s，仿真结果可随时由Scope模块显示出清晰图形以供分析。100min内求得的各物理量随时间变化如图5至8所示。

图5 反应性随时间的变化

图6 功率随时间的变化

图7 先驱核浓度随时间的变化

图8 温度随时间的变化
为了检验仿真结果的准确性，图5和图6中不但给出了仿真的结果，同时还给出了参考文献[5]给出的结果，从图中可以看出，仿真得到的反应性和功率的结果与参考文献[5]中采用Fortran语言编程求得的结果吻合的非常好，这说明基于Simulink的动态响应仿真具有很好的计算精度，而且可以避开繁琐的程序编写，具有快捷、高效的特点。
由图5至图8可以看到，在引入小的正阶跃反应性后反应堆功率（中子密度）和缓发中子先驱核浓度开始缓慢增长，因而反应堆堆芯温度也缓慢升高，由于存在温度反馈且温度系数为负数，温度升高导致了反应性减小。当反应性减小到零，中子密度和缓发中子先驱核浓度也就不再增加达到了最大值。此时，由于裂变反应释放的热量仍大于堆芯传出的热量，反应堆温度继续上升，反应性继续减小，其值将小于零中子密度和缓发中子先驱核浓度开始减小。直到由于裂变反应释放的热量仍小于堆芯传出的热量，反应堆温度开始下降，反应性开始增大。
4 结 论
在实际反应堆运行过程中，反应性一般是随时间变化的。在中子动力学分析中注重的是短时间内中子密度随着时间的变化关系，因此主要考虑温度对反应性的影响。考虑简单的温度反馈之后，点堆动力学方程变为非线性方程组，直接编程求解比较困难。本文采用比较简单且常用的绝热模型对方程进行化简，并应用幂级数法编程计算得到结果，分析该结果得到在引入较小的正阶跃反应性后，中子的增殖或功率的增长需要较长的时间，在开始的较短时间内可以忽略反应堆的热量损失，这时应用绝热模型可以得到非常准确的结果的。但随着时间的增长，反应堆有足够的时间将热量传递出去，此时采用绝热模型后求得各物理量随时间的变化速度明显快于实际情况，绝热模型不再适用。基于Simulink的动态响应仿真不但可以避开繁琐的程序编写，而且具有很好的计算精度。从仿真结果看出，当温度系数为负数，反应堆出现反应性扰动时，堆内各物理量能通过负反馈效应自行调整达到另一平衡态，反应堆具有内在的稳定性，这对反应堆的调节和运行安全都具有重大意义。
参考文献
[1] 蔡章生. 核动力反应堆中子动力学[M]. 北京：国防工业出版社, 2005.
[2] 黄祖洽. 核反应堆动力学基础[M]. 第二版. 北京：北京大学出版社，2007.
[3] Aboanber A E, Hamada Y M. Power series solution (pws) of nuclear reactor dynamics with newtonian temperature feedback[J]. Annals of Nuclear Energy, 2003, 30: 1111-1122.
[4] Abdallah A Nahla. An efficient technique for the point reactor kinetics equations with Newtonian temperature feedback effects[J]. Annals of Nuclear Energy, 2011, 38: 2810-2817.
[5] H.L.Dodds,Jr, R.M.Westfall. SKINATH—A Computer Program for Solving the Reactor Point Kinetics Equations with Simple Thermal-Hydraulic Feedback[R]. Printed in the United States of America, Available from National Technical Information Service. 1984.
[6] 王乔，黎浩峰，陈文振. 基于Simulink输入阶跃反应性时有温度和毒物反馈的反应堆动态响应仿真[J]. 原子能科学技术. 2009，43(9)：823-827.
[7] Hamieh S D, Saidinezhad M. Analytical solution of the point reactor kinetics equations with temperature feedback[J]. Annals of Nuclear Energy, 2012, 42: 148-152.