核反应堆物理课设

课程设计(综合实验)报告名称：核反响堆物理分析题目：利用双群理论求解堆芯参数院系：班级：学号：学生姓名：指导教师：111111111设计周数：11111111成绩：课程设计(综合实验)的目的与要求1课程设计的要求课程设计是重要的实践教学环节。它是根据教学方案的要求，在教师指导下对学生进行的阶段根底或专业技术训练，该实践环节着重培养学生综合分析和解决实际问题的方法与能力，实现由知识向智能的初步转化；是对前期理论与实践教学效果的检验，也是对学生综合分析能力与独立工作能力的核反响堆物理分析课程设计的目的是对理论课上学过的理论知识进行践应用，进而加深对前期理论知识的学习，是对学生综合运用核反响堆物理分析知识和思想方法的综合检验过程。2课程设计要求核反响堆物理分析课程设计的要求有如下几点：〔1〕学生必须修完课程设计的先修课程，才有资格做课程设计。〔2〕明确课程设计的目的和重要性，认真领会课程设计的题目，读懂课程设计指导书的要求，学会设计的根本方法与步骤，积极认真地做好准备工作。〔3〕在课程设计中，学会如何运用前修知识与收集、归纳相关资料解决具体问题的方法。〔4〕学生必须在指导教师指导下独立完成设计任务，严禁抄袭、找人代做等，一经发现成绩记零分，按考试作弊处理。〔5〕课程设计报告学校有统一格式，学生必须按照此格式填写，说明书、计算书要求简洁、通顺、计算正确，图纸表达内容完整、清楚、标准。二、设计〔实验〕正文1.课程设计题目：二维XY几何模型，中间是堆芯〔裂变区〕：-50cm<x<+50cm,-55cm<y<+55cm,四周是15cm厚反射层，堆芯及反射层的群常数见下表：表1堆芯各种参数群常数堆芯反射层第1群第2群第1群第2群1.00.01.00.00.00750.16510.00.00.01210.1210.00040.0200.0241-0.0493-1.260.3541.1300.166利用数值求解下面问题，计算其快群和热群中子通量密度空间分布及有效增殖因子。提示：采用内外迭代方法求解；给出源迭代过程；给出内迭代过程并给出堆芯双群方程、反射层双群方程以及边界条件，详见课本第五章；建立数学模型，由于问题具有对称性，因此，只计算1/4堆芯即可。可用赛德尔迭代法求解；〔详见课本第五章〕；给出计算流程编写程序，给出结果，并用OREGIN画图把结果表现出来；结果的分析讨论；写出报告。2.1差分方程的建立为简化计算，进行以下假设：〔1〕二维XY几何模型是一个整体，即堆芯和反射层交界处的中子通量密度和中子流连续；〔2〕认为该问题中反射层厚度已包含外推边界；〔3〕认为该问题中堆芯和反射层对称分布。在带反射层的多区双群扩散方程的数值解法中，通常是通过外-内迭代过程来求解。外迭代是一个通过迭代求特征值的迭代过程，所以又成为源迭代。内迭代那么是对源迭代中出现的扩散方程进行具体的数值求解。计算中假设不考虑中子自低能群的向上散射，多群扩散方程可以写为：〔1〕略取上下标，该方程可以写成以下标准形式：〔2〕本设计讨二维〔x，y〕情况下的差分方程。对所取平面首先用x=，，…，，…，y=，，…，，…，直线族把平面分成许多矩形网格〔如图一〕。交点〔，〕称为节点或网点，如图一中共有〔N+1〕(M+1)个节点，称为网距。图1网点的划分图2〔i，j〕节点示意图按照数值解法的相关要求和定义，讨论如图二所示的〔i，j〕节点，根据二维几何情况下扩散方程导出〔i，j〕节点的差分方程为〔3〕式中系数为〔4〕对于边界上的点需要根据边界条件来确定，本设计由于反射层已经包括了外推距离，所以外围上的点的中子通量密度均设为0。上述〔3〕和〔4〕式便是求解单群中子扩散的差分方程组，其系数，，…，等均可以事先计算求得，对每次外迭代，其右端源项也是的。因此，它是一个含有的代数方程组，对于二维问题，其系数矩阵是一个五队角的矩阵，可以很方便的应用通常的线性方程组的求解方法进行求解。2.2源迭代过程根据连续性的假设条件，堆芯的双群方程与放射层得双群方程可以归结为下面两个方程：快群中子扩散方程：〔5〕热群中子扩散方程：〔6〕设〔7〕式〔5〕是一个关于的齐次方程组。式〔5〕除了中子通量密度外，和也是未知参数，并且这两个参数又与中子通量密度有关。因此我们假定一个初始的裂变源分布，认为在整个堆芯都等于1或某个给定常数，并猜想一个初始的值，同时把作为初始迭代源项代入到方程〔5〕的右端中。这样，方程〔1〕的右端便是一个项，方程〔1〕便成为一个普通的非齐次方程组，它可以用通常的数值方法求解。假定由它求出的中子通量密度分布为，把代入方程〔6〕，由方程〔6〕数值求解求出，将求得的、带入〔7〕式便可求得第二次迭代裂变中子源分布，并由它可以求得有效增殖因数的新的估计值和第二代迭代源项来，……。依此类推，逐次的迭代下去，对于第n次迭代计算有〔8〕式中：〔9〕根据的物理意义〔它等于新生一代的裂变总中子数与上一代源中子总数之比〕可以由下式算出的估计值〔10〕方程〔3〕，〔4〕和〔5〕便是源迭代法的计算公式。从物理上看，上面所确定的迭代过程相当于将反响堆内的中子数人为地区分为不同的“代”，每一“代”中子有裂变形成之时作为起始。每一次源迭代那么相当于老一代源中子通过扩散、慢化、俘获、裂变，最后又形成新一代源中子的过程。式〔5〕中的显然反映了中子每经过一代的倍增。因而可以预期，经过了足够的循环代数后，任意的初始中子源必将在堆内形成一个稳定的最终中子源分布。多群扩散方程的特征值问题在数学上也已被详细研究过，并证明了上述迭代过程的收敛性。2.3内迭代过程对于一次内迭代过程，假定一个初始的裂变源分布及一个初始的数值，得如下计算步骤：〔1〕芯部中子通量密度的计算芯部快群〔第1群〕芯部快群中子扩散方程为：内迭代所需群常数为：内迭代相应的系数为：由节点〔i，j〕处的差分方程为：最后得相应点处的快中子通量密度：边界局部的点的中子通量密度认为为0，不必计算。芯部热群〔第2群〕芯部热群中子扩散方程为：内迭代所需群常数为：内迭代相应的系数为：由节点〔i，j〕处的差分方程为：最后得相应点处的热中子通量密度：边界局部的点的中子通量密度认为为0，不必计算。〔2〕反射层的中子通量密度反射层快群〔第一群〕和热群〔第二群〕的扩散方程为：2.4有效增殖因子的计算将上述芯部计算过程得到的快中子和热中子的通量密度和同时代入下式根据的计算式得分别将本次中子源和中子增殖因子与上一次的中子源和中子增殖因子进行比拟，观察是否满足以下收敛准那么特征值收敛准那么:中子通量密度分布收敛准那么：在本设计中，和分别取为和。经过计算，假设上述两个收敛准那么同时得到满足，那么此次计算得到的中子通量密度分布和有效增殖因子为所求的解；反之，那么将此次计算得到的中子通量密度分布和有效增殖因子作为新的裂变源分布和新的数值再重复上述计算步骤。2.5建模由于问题具有对称性，因此，故计算1/4堆芯即可〔如以下图3所示〕图3示意图2.6计算结果分析通过计算机程序求得：可以近似认为处于临界状态。三、课程设计〔综合实验〕总结或结论1结论(1)带有反射层反响堆内中子源的分布由中心向四周逐渐减小，在反射层内为零。但是其在靠近反射层附近出现一个小的峰值，尤其是在边角这个峰值最大，出现这种现象的原因可能是快中子进入反射层后被慢化为慢中子回到堆芯，在靠近反射层的地方产生较大的反响界面，产生更多的热中子。(2)如图7、8、9、10所示，热中子通量密度曲线在堆芯与反射层的交界处后有一个明显的峰值，而快中子通量密度曲线有一个接近水平的平缓变化的段。主要由于反射层对热中子的吸收比堆芯弱，但对快中子的慢化比堆芯强，由堆芯泄泄漏的快中子被反射层慢化成热中子的缘故。2总结通过此次核反响堆物理分析课程设计，我收获了很多东西。首先，我对课本上很多当初不懂的知识有了更深刻的认识；其次，对于很多以前课本上的结论公式，知道了其来源和用处；然后，也对以前学过的MATLAB软件的编程和数据处理有了更好的掌握，也会通过ORIGIN软件画图，节省了很多的计算过程；并且通过课程设计了解了这门的处理问题的思想，同时也敬佩做研究人的那种处理抽象问题的能力。最后，感谢老师的辛苦教导。四、参考文献[1]谢仲生主编.《核反响堆物理分析》.修订本.西安交通大学出版社.原子能出版社.2012.7[2]刘卫国主编.《MATLAB程序设计与运用》.第二版.高等教育出版社.2006.7[3]刘卫国主编.《MATLAB程序设计教程》.第二版.中国水利水电出版社.2010.2附录〔设计流程图、程序、表格、数据等〕几何或成分调整初始，n=0几何或成分调整初始，n=0群常数群常数n=n+1n=n+1临界计算临界计算g=0g=0初始；g=g+1初始；g=g+1临界搜索临界搜索差分方程的迭代求解外迭代差分方程的迭代求解外迭代内迭代内迭代否g=G否g=G是是计算计算否否否是否是是是停机停机图4〔二〕三维图像图41/4带反射层的堆芯内中子源的空间分布图51/4带反射层的堆芯内中子源沿X轴方向分布图61/4带反射层的堆芯内中子源沿Y轴方向分布图71/4带反射层的堆芯内中子通量密沿X轴方向分布图81/4带反射层的堆芯内中子通量密沿Y轴方向分布图91/4带反射层的堆芯内快中子通量密度的空间分布图101/4带反射层的堆芯内热中子通量密度的空间分布〔三〕程序%定义堆芯参数fory=1:55;forx=1:50;N(x,y).X1=1;N(x,y).vEf1=0.0075;N(x,y).Ea1=0.0121;N(x,y).Es1=0.0241;N(x,y).D1=1.26;N(x,y).Er1=N(x,y).Ea1+N(x,y).Es1;N(x,y).X2=0;N(x,y).vEf2=0.1651;N(x,y).Ea2=0.121;N(x,y).D2=0.354;N(x,y).Er2=N(x,y).Ea2;endend%定义反射层参数fory=1:70;forx=1:65;ifx>50|y>55N(x,y).X1=1;N(x,y).vEf1=0;N(x,y).Ea1=0.0004;N(x,y).Es1=0.0493;N(x,y).D1=1.13;N(x,y).Er1=N(x,y).Ea1+N(x,y).Es1;N(x,y).X2=0;N(x,y).vEf2=0;N(x,y).Ea2=0.02;N(x,y).D2=0.166;N(x,y).Er2=N(x,y).Ea2;endendend%设置迭代初值Q1=0.45*ones(65,70);keff1=12.8;F11=ones(65,70);F12=ones(65,70);F21=ones(65,70);F22=ones(65,70);%定义收敛准那么初值e1=1;e2=1;e3=1;e4=1;%设置边界条件F11(65,:)=0;F12(65,:)=0;F11(:,70)=0;F12(:,70)=0;F21(65,:)=0;F22(65,:)=0;F21(:,70)=0;F22(:,70)=0;%迭代过程while(e1>1e-5&e2>1e-4)%内迭代过程while(e3>0.1)fory=2:69;forx=2:64;%快群a(x,y)=-(N(x-1,y-1).D1+N(x+1,y-1).D1)/2;b(x,y)=-(N(x-1,y-1).D1+N(x-1,y+1).D1)/2;c(x,y)=-(N(x+1,y-1).D1+N(x+1,y+1).D1)/2;d(x,y)=-(N(x+1,y+1).D1+N(x-1,y+1).D1)/2;e(x,y)=1/4*(N(x-1,y-1).Er1+N(x+1,y-1).Er1+N(x+1,y+1).Er1+N(x-1,y+1).Er1)-a(x,y)-b(x,y)-c(x,y)-d(x,y);s1=N(x-1,y-1).X1*(N(x-1,y-1).vEf1*F11(x-1,y-1)+N(x-1,y-1).vEf2*F21(x-1,y-1))/keff1;s2=N(x+1,y-1).X1*(N(x+1,y-1).vEf1*F11(x+1,y-1)+N(x+1,y-1).vEf2*F21(x+1,y-1))/keff1;s3=N(x+1,y+1).X1*(N(x+1,y+1).vEf1*F11(x+1,y+1)+N(x+1,y+1).vEf2*F21(x+1,y+1))/keff1;s4=N(x-1,y+1).X1*(N(x-1,y+1).vEf1*F11(x-1,y+1)+N(x-1,y+1).vEf2*F21(x-1,y+1))/keff1;s(x,y)=(s1+s2+s3+s4)/4;F12(x,y)=-(a(x,y)*F11(x,y-1)+b(x,y)*F11(x-1,y)+c(x,y)*F11(x+1,y)+d(x,y)*F11(x,y+1)-s(x,y))/e(x,y);endende3=sum(sum(abs(F12-F11)))/(sum(sum(abs(F12))));endF11=F12;%热群while(e4>0.1)fory=2:69;forx=2:64;a(x,y)=-(N(x-1,y-1).D2+N(x+1,y-1).D2)/2;b(x,y)=-(N(x-1,y-1).D2+N(x-1,y+1).D2)/2;c(x,y)=-(N(x+1,y-1).D2+N(x+1,y+1).D2)/2;d(x,y)=-(N(x+1,y+1).D2+N(x-1,y+1).D2)/2;e(x,y)=1/4*(N(x-1,y-1).Er2+N(x+1,y-1).Er2+N(x+1,y+1).Er2+N(x-1,y+1).Er2)-a(x,y)-b(x,y)-c(x,y)-d(x,y);s1=N(x-1,y-1).Es1*F11(x-1,y-1);s2=N(x+1,y-1).Es1*F11(x+1,y-1);s3=N(x+1,y+1).Es1*F11(x+1,y+1);s4=N(x-1,y+1).Es1*F11(x-1,y+1);s(x,y)=(s1+s2+s3+s4)/4;F22(x,y)=-(a(x,y)*F21(x,y-1)+b(x,y)*F21(x-1,y)+c(x,y)*F21(x+1,y)+d(x,y)*F21(x,y+1)-s(x,y))/e(x,y);endende4=sum(sum(abs(F22-F21)))/sum(sum(abs(F22)));endF21=F22;fory=1:70forx=1:65Q2(x,y)=(N(x,y).vEf1)*F11(x,y)+(N(x,y).vEf2)*F21(x,y);endendkeff2=keff1*(sum(sum(Q2))/sum(sum(Q1)));e1=abs((keff2