核反應(yīng)堆物理課設(shè)

1、 程設(shè)計(jì)（綜合實(shí)驗(yàn)）報(bào)告 課程設(shè)計(jì)(綜合實(shí)驗(yàn))報(bào)告名 稱： 核反應(yīng)堆物理分析 題 目：利用雙群理論求解堆芯參數(shù)院 系： 1111111111111 班 級： 111111111111111 學(xué) 號： 111111111111 學(xué)生姓名： 11111111111 指導(dǎo)教師： 111111111 設(shè)計(jì)周數(shù)： 11111111 成 績： 1、 課程設(shè)計(jì)(綜合實(shí)驗(yàn))的目的與要求1 課程設(shè)計(jì)的要求 課程設(shè)計(jì)是重要的實(shí)踐教學(xué)環(huán)節(jié)。它是根據(jù)教學(xué)計(jì)劃的要求，在教師指導(dǎo)下對學(xué)生進(jìn)行的階段基礎(chǔ)或?qū)I(yè)技術(shù)訓(xùn)練，該實(shí)踐環(huán)節(jié)著重培養(yǎng)學(xué)生綜合分析和解決實(shí)際問題的方法與能力，實(shí)現(xiàn)由知識向智能的初步轉(zhuǎn)化；是對前期理論與實(shí)踐教

2、學(xué)效果的檢驗(yàn)，也是對學(xué)生綜合分析能力與獨(dú)立工作能力的核反應(yīng)堆物理分析課程設(shè)計(jì)的目的是對理論課上學(xué)過的理論知識進(jìn)行踐應(yīng)用，進(jìn)而加深對前期理論知識的學(xué)習(xí)，是對學(xué)生綜合運(yùn)用核反應(yīng)堆物理分析知識和思想方法的綜合檢驗(yàn)過程。2 課程設(shè)計(jì)要求核反應(yīng)堆物理分析課程設(shè)計(jì)的要求有如下幾點(diǎn)：（1）學(xué)生必須修完課程設(shè)計(jì)的先修課程，才有資格做課程設(shè)計(jì)。（2）明確課程設(shè)計(jì)的目的和重要性，認(rèn)真領(lǐng)會課程設(shè)計(jì)的題目，讀懂課程設(shè)計(jì)指導(dǎo)書的要求，學(xué)會設(shè)計(jì)的基本方法與步驟，積極認(rèn)真地做好準(zhǔn)備工作。（3）在課程設(shè)計(jì)中，學(xué)會如何運(yùn)用前修知識與收集、歸納相關(guān)資料解決具體問題的方法。（4）學(xué)生必須在指導(dǎo)教師指導(dǎo)下獨(dú)立完成設(shè)計(jì)任務(wù)，嚴(yán)禁抄襲、

3、找人代做等，一經(jīng)發(fā)現(xiàn)成績記零分，按考試作弊處理。（5）課程設(shè)計(jì)報(bào)告學(xué)校有統(tǒng)一格式，學(xué)生必須按照此格式填寫，說明書、計(jì)算書要求簡潔、通順、計(jì)算正確，圖紙表達(dá)內(nèi)容完整、清楚、規(guī)范。二、設(shè)計(jì)（實(shí)驗(yàn)）正文1.課程設(shè)計(jì)題目：二維XY幾何模型，中間是堆芯（裂變區(qū)）：-50cm<x<+50cm, -55cm<y<+55cm,四周是15cm厚反射層，堆芯及反射層的群常數(shù)見下表：表1 堆芯各種參數(shù)群常數(shù)堆芯反射層第1群第2群第1群第2群1.00.01.00.00.00750.16510.00.00.01210.1210.00040.0200.0241-0.0493-1.260.3541.

4、1300.166利用數(shù)值求解下面問題，計(jì)算其快群和熱群中子通量密度空間分布及有效增殖因子。提示：（1） 采用內(nèi)外迭代方法求解；（2） 給出源迭代過程；（3） 給出內(nèi)迭代過程并給出堆芯雙群方程、反射層雙群方程以及邊界條件，詳見課本第五章；（4） 建立數(shù)學(xué)模型，由于問題具有對稱性，因此，只計(jì)算1/4堆芯即可。（5） 可用賽德爾迭代法求解；（詳見課本第五章）；（6） 給出計(jì)算流程（7） 編寫程序，給出結(jié)果，并用OREGIN畫圖把結(jié)果表現(xiàn)出來；（8） 結(jié)果的分析討論；（9） 寫出報(bào)告。2.1差分方程的建立為簡化計(jì)算，進(jìn)行以下假設(shè)：（1）二維XY幾何模型是一個整體，即堆芯和反射層交界處的中子通量密度和中

5、子流連續(xù)；（2）認(rèn)為該問題中反射層厚度已包含外推邊界；（3）認(rèn)為該問題中堆芯和反射層對稱分布。在帶反射層的多區(qū)雙群擴(kuò)散方程的數(shù)值解法中，通常是通過外-內(nèi)迭代過程來求解。外迭代是一個通過迭代求特征值的迭代過程，所以又成為源迭代。內(nèi)迭代則是對源迭代中出現(xiàn)的擴(kuò)散方程進(jìn)行具體的數(shù)值求解。計(jì)算中若不考慮中子自低能群的向上散射，多群擴(kuò)散方程可以寫為： （1）略取上下標(biāo)，該方程可以寫成以下標(biāo)準(zhǔn)形式： （2）本設(shè)計(jì)討二維（x，y）情況下的差分方程。對所取平面首先用 x = ， y = ，直線族把平面分成許多矩形網(wǎng)格（如圖一）。交點(diǎn)（，）稱為節(jié)點(diǎn)或網(wǎng)點(diǎn)，如圖一中共有（N+1）(M+1)個節(jié)點(diǎn)， 稱為網(wǎng)距。 圖1

6、 網(wǎng)點(diǎn)的劃分圖2 （i，j）節(jié)點(diǎn)示意圖按照數(shù)值解法的相關(guān)要求和定義，討論如圖二所示的（i，j）節(jié)點(diǎn)，根據(jù)二維幾何情況下擴(kuò)散方程導(dǎo)出（i，j）節(jié)點(diǎn)的差分方程為 （3）式中系數(shù)為 （4）對于邊界上的點(diǎn)需要根據(jù)邊界條件來確定，本設(shè)計(jì)由于反射層已經(jīng)包括了外推距離，所以外圍上的點(diǎn)的中子通量密度均設(shè)為0。上述（3）和（4）式便是求解單群中子擴(kuò)散的差分方程組，其系數(shù)，等均可以事先計(jì)算求得，對每次外迭代，其右端源項(xiàng)也是已知的。因此，它是一個含有的代數(shù)方程組，對于二維問題，其系數(shù)矩陣是一個五隊(duì)角的矩陣，可以很方便的應(yīng)用通常的線性方程組的求解方法進(jìn)行求解。2.2源迭代過程根據(jù)連續(xù)性的假設(shè)條件，堆芯的雙群方程與放射

7、層得雙群方程可以歸結(jié)為下面兩個方程：快群中子擴(kuò)散方程： （5）熱群中子擴(kuò)散方程： （6） 設(shè) （7）式（5）是一個關(guān)于的齊次方程組。式（5）除了中子通量密度外，和也是未知參數(shù)，并且這兩個參數(shù)又與中子通量密度有關(guān)。因此我們假定一個初始的裂變源分布，認(rèn)為在整個堆芯都等于1或某個給定常數(shù)，并猜測一個初始的值，同時把作為初始迭代源項(xiàng)代入到方程（5）的右端中。這樣，方程（1）的右端便是一個已知項(xiàng)，方程（1）便成為一個普通的非齊次方程組，它可以用通常的數(shù)值方法求解。假定由它求出的中子通量密度分布為，把代入方程（6），由方程（6）數(shù)值求解求出，將求得的、帶入（7）式便可求得第二次迭代裂變中子源分布，并由它可

8、以求得有效增殖因數(shù)的新的估計(jì)值和第二代迭代源項(xiàng)來，。依此類推，逐次的迭代下去，對于第n次迭代計(jì)算有 （8）式中： （9） 根據(jù)的物理意義（它等于新生一代的裂變總中子數(shù)與上一代源中子總數(shù)之比）可以由下式算出的估計(jì)值 （10）方程（3），（4）和（5）便是源迭代法的計(jì)算公式。從物理上看，上面所確定的迭代過程相當(dāng)于將反應(yīng)堆內(nèi)的中子數(shù)人為地區(qū)分為不同的“代”，每一“代”中子有裂變形成之時作為起始。每一次源迭代則相當(dāng)于老一代源中子通過擴(kuò)散、慢化、俘獲、裂變，最后又形成新一代源中子的過程。式（5）中的顯然反映了中子每經(jīng)過一代的倍增。因而可以預(yù)期，經(jīng)過了足夠的循環(huán)代數(shù)后，任意的初始中子源必將在堆內(nèi)形成一個穩(wěn)

9、定的最終中子源分布。多群擴(kuò)散方程的特征值問題在數(shù)學(xué)上也已被詳細(xì)研究過，并證明了上述迭代過程的收斂性。2.3內(nèi)迭代過程對于一次內(nèi)迭代過程，假定一個初始的裂變源分布及一個初始的數(shù)值，得如下計(jì)算步驟：（1）芯部中子通量密度的計(jì)算 芯部快群（第1群）芯部快群中子擴(kuò)散方程為：內(nèi)迭代所需群常數(shù)為： 內(nèi)迭代相應(yīng)的系數(shù)為： 由節(jié)點(diǎn)（i，j）處的差分方程為：最后得相應(yīng)點(diǎn)處的快中子通量密度： 邊界部分的點(diǎn)的中子通量密度認(rèn)為為0，不必計(jì)算。 芯部熱群（第2群）芯部熱群中子擴(kuò)散方程為：內(nèi)迭代所需群常數(shù)為： 內(nèi)迭代相應(yīng)的系數(shù)為： 由節(jié)點(diǎn)（i，j）處的差分方程為：最后得相應(yīng)點(diǎn)處的熱中子通量密度： 邊界部分的點(diǎn)的中子通量密

10、度認(rèn)為為0，不必計(jì)算。（2）反射層的中子通量密度反射層快群（第一群）和熱群（第二群）的擴(kuò)散方程為： 2.4有效增殖因子的計(jì)算將上述芯部計(jì)算過程得到的快中子和熱中子的通量密度和 同時代入下式根據(jù)的計(jì)算式得分別將本次中子源和中子增殖因子與上一次的中子源和中子增殖因子進(jìn)行比較，觀察是否滿足下列收斂準(zhǔn)則（1） 特征值收斂準(zhǔn)則: （2） 中子通量密度分布收斂準(zhǔn)則：在本設(shè)計(jì)中，和分別取為和。經(jīng)過計(jì)算，若上述兩個收斂準(zhǔn)則同時得到滿足，則此次計(jì)算得到的中子通量密度分布和有效增殖因子為所求的解；反之，則將此次計(jì)算得到的中子通量密度分布和有效增殖因子作為新的裂變源分布和新的數(shù)值再重復(fù)上述計(jì)算步驟。2.5建模由于問

11、題具有對稱性，因此，故計(jì)算1/4 堆芯即可（如下圖3所示） 圖3 示意圖2.6 計(jì)算結(jié)果分析通過計(jì)算機(jī)程序求得： 可以近似認(rèn)為處于臨界狀態(tài)。三、課程設(shè)計(jì)（綜合實(shí)驗(yàn)）總結(jié)或結(jié)論1 結(jié)論(1)帶有反射層反應(yīng)堆內(nèi)中子源的分布由中心向四周逐漸減小，在反射層內(nèi)為零。但是其在靠近反射層附近出現(xiàn)一個小的峰值，尤其是在邊角這個峰值最大，出現(xiàn)這種現(xiàn)象的原因可能是快中子進(jìn)入反射層后被慢化為慢中子回到堆芯，在靠近反射層的地方產(chǎn)生較大的反應(yīng)界面，產(chǎn)生更多的熱中子。 (2)如圖7、8、9、10所示，熱中子通量密度曲線在堆芯與反射層的交界處后有一個明顯的峰值，而快中子通量密度曲線有一個接近水平的平緩變化的段。主要由于反射

12、層對熱中子的吸收比堆芯弱，但對快中子的慢化比堆芯強(qiáng)，由堆芯泄泄漏的快中子被反射層慢化成熱中子的緣故。2 總結(jié)通過此次核反應(yīng)堆物理分析課程設(shè)計(jì)，我收獲了很多東西。首先，我對課本上很多當(dāng)初不懂的知識有了更深刻的認(rèn)識；其次，對于很多以前課本上的結(jié)論公式，知道了其來源和用處；然后，也對以前學(xué)過的MATLAB軟件的編程和數(shù)據(jù)處理有了更好的掌握，也會通過ORIGIN軟件畫圖，節(jié)省了很多的計(jì)算過程；并且通過課程設(shè)計(jì)了解了這門的處理問題的思想，同時也敬佩做研究人的那種處理抽象問題的能力。最后，感謝老師的辛苦教導(dǎo)。四、參考文獻(xiàn)1謝仲生 主編. 核反應(yīng)堆物理分析. 修訂本.西安交通大學(xué)出版社.原子能出版社.201

13、2.72劉衛(wèi)國 主編 .MATLAB程序設(shè)計(jì)與運(yùn)用. 第二版.高等教育出版社. 2006.73劉衛(wèi)國 主編 .MATLAB程序設(shè)計(jì)教程.第二版.中國水利水電出版社.2010.2附錄（設(shè)計(jì)流程圖、程序、表格、數(shù)據(jù)等）幾何或成分調(diào)整初始，n=0群常數(shù)n=n+1臨界計(jì)算g=0初始；g=g+1臨 界 搜 索差分方程的迭代求解外迭代內(nèi)迭代否g=G是計(jì)算否否是是停機(jī)圖4（二）三維圖像 圖4 1/4帶反射層的堆芯內(nèi)中子源的空間分布 圖5 1/4帶反射層的堆芯內(nèi)中子源沿X軸方向分布 圖6 1/4帶反射層的堆芯內(nèi)中子源沿Y軸方向分布 圖7 1/4帶反射層的堆芯內(nèi)中子通量密沿X軸方向分布 圖8 1/4帶反射層的堆

14、芯內(nèi)中子通量密沿Y軸方向分布圖9 1/4帶反射層的堆芯內(nèi)快中子通量密度的空間分布圖10 1/4帶反射層的堆芯內(nèi)熱中子通量密度的空間分布（三）程序%定義堆芯參數(shù)for y=1:55; for x=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)

15、.Ea2; endend%定義反射層參數(shù)for y=1:70; for x=1:65; if x>50|y>55 N(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; end end end %設(shè)置迭代初值Q1=0.45*ones(65,70);keff1=

16、12.8;F11=ones(65,70);F12=ones(65,70);F21=ones(65,70);F22=ones(65,70);%定義收斂準(zhǔn)則初值e1=1;e2=1;e3=1;e4=1;%設(shè)置邊界條件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) %內(nèi)迭代過程 while(e3>0.1) for y=2:69; for x=2:64; %快群 a(x,y)=-

17、(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,

18、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

19、+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); end end e3=sum(sum(abs(F12-F11)/(sum(sum(abs(F12); end F11=F12;%熱群 while(e4>0.1) for y=2:69; for x=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

20、; 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*

21、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); end end e4=sum(sum(abs(F22-F21)/sum(sum(abs(F22); end F21=F22; for y=1:70 for x=1:65 Q2(x,y)=(N(x,y).vEf1)*F11(x,y)+(N(x,y).vEf2)*F21(x,y); end end keff2=keff1*(sum(sum(Q2)/sum(sum(Q1); e1=abs(keff2-keff1)/keff2); e2=abs(sum(sum(Q2)-sum(sum(Q1)/sum(sum(Q2); Q1=Q2; keff1=keff2; e3=1; e4=1;endkeff1%畫圖程序figure(1);mesh(F11);figure(2);mesh(F21);figure(3);mesh(Q1);a=linspace(1,65,65);b=linspace