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

課程設(shè)計(jì)(綜合實(shí)驗(yàn))報(bào)告名稱：核反響堆物理分析題目：利用雙群理論求解堆芯參數(shù)院系：班級(jí)：學(xué)號(hào)：學(xué)生姓名：指導(dǎo)教師：111111111設(shè)計(jì)周數(shù)：11111111成績(jī)：課程設(shè)計(jì)(綜合實(shí)驗(yàn))的目的與要求1課程設(shè)計(jì)的要求課程設(shè)計(jì)是重要的實(shí)踐教學(xué)環(huán)節(jié)。它是根據(jù)教學(xué)方案的要求，在教師指導(dǎo)下對(duì)學(xué)生進(jìn)行的階段根底或?qū)I(yè)技術(shù)訓(xùn)練，該實(shí)踐環(huán)節(jié)著重培養(yǎng)學(xué)生綜合分析和解決實(shí)際問題的方法與能力，實(shí)現(xiàn)由知識(shí)向智能的初步轉(zhuǎn)化；是對(duì)前期理論與實(shí)踐教學(xué)效果的檢驗(yàn)，也是對(duì)學(xué)生綜合分析能力與獨(dú)立工作能力的核反響堆物理分析課程設(shè)計(jì)的目的是對(duì)理論課上學(xué)過的理論知識(shí)進(jìn)行踐應(yīng)用，進(jìn)而加深對(duì)前期理論知識(shí)的學(xué)習(xí)，是對(duì)學(xué)生綜合運(yùn)用核反響堆物理分析知識(shí)和思想方法的綜合檢驗(yàn)過程。2課程設(shè)計(jì)要求核反響堆物理分析課程設(shè)計(jì)的要求有如下幾點(diǎn)：〔1〕學(xué)生必須修完課程設(shè)計(jì)的先修課程，才有資格做課程設(shè)計(jì)?！?〕明確課程設(shè)計(jì)的目的和重要性，認(rèn)真領(lǐng)會(huì)課程設(shè)計(jì)的題目，讀懂課程設(shè)計(jì)指導(dǎo)書的要求，學(xué)會(huì)設(shè)計(jì)的根本方法與步驟，積極認(rèn)真地做好準(zhǔn)備工作。〔3〕在課程設(shè)計(jì)中，學(xué)會(huì)如何運(yùn)用前修知識(shí)與收集、歸納相關(guān)資料解決具體問題的方法。〔4〕學(xué)生必須在指導(dǎo)教師指導(dǎo)下獨(dú)立完成設(shè)計(jì)任務(wù)，嚴(yán)禁抄襲、找人代做等，一經(jīng)發(fā)現(xiàn)成績(jī)記零分，按考試作弊處理?！?〕課程設(shè)計(jì)報(bào)告學(xué)校有統(tǒng)一格式，學(xué)生必須按照此格式填寫，說明書、計(jì)算書要求簡(jiǎn)潔、通順、計(jì)算正確，圖紙表達(dá)內(nèi)容完整、清楚、標(biāo)準(zhǔn)。二、設(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.1300.166利用數(shù)值求解下面問題，計(jì)算其快群和熱群中子通量密度空間分布及有效增殖因子。提示：采用內(nèi)外迭代方法求解；給出源迭代過程；給出內(nèi)迭代過程并給出堆芯雙群方程、反射層雙群方程以及邊界條件，詳見課本第五章；建立數(shù)學(xué)模型，由于問題具有對(duì)稱性，因此，只計(jì)算1/4堆芯即可?？捎觅惖聽柕ㄇ蠼猓弧苍斠娬n本第五章〕；給出計(jì)算流程編寫程序，給出結(jié)果，并用OREGIN畫圖把結(jié)果表現(xiàn)出來；結(jié)果的分析討論；寫出報(bào)告。2.1差分方程的建立為簡(jiǎn)化計(jì)算，進(jìn)行以下假設(shè)：〔1〕二維XY幾何模型是一個(gè)整體，即堆芯和反射層交界處的中子通量密度和中子流連續(xù)；〔2〕認(rèn)為該問題中反射層厚度已包含外推邊界；〔3〕認(rèn)為該問題中堆芯和反射層對(duì)稱分布。在帶反射層的多區(qū)雙群擴(kuò)散方程的數(shù)值解法中，通常是通過外-內(nèi)迭代過程來求解。外迭代是一個(gè)通過迭代求特征值的迭代過程，所以又成為源迭代。內(nèi)迭代那么是對(duì)源迭代中出現(xiàn)的擴(kuò)散方程進(jìn)行具體的數(shù)值求解。計(jì)算中假設(shè)不考慮中子自低能群的向上散射，多群擴(kuò)散方程可以寫為：〔1〕略取上下標(biāo)，該方程可以寫成以下標(biāo)準(zhǔn)形式：〔2〕本設(shè)計(jì)討二維〔x，y〕情況下的差分方程。對(duì)所取平面首先用x=，，…，，…，y=，，…，，…，直線族把平面分成許多矩形網(wǎng)格〔如圖一〕。交點(diǎn)〔，〕稱為節(jié)點(diǎn)或網(wǎng)點(diǎn)，如圖一中共有〔N+1〕(M+1)個(gè)節(jié)點(diǎn)，稱為網(wǎng)距。圖1網(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〕對(duì)于邊界上的點(diǎn)需要根據(jù)邊界條件來確定，本設(shè)計(jì)由于反射層已經(jīng)包括了外推距離，所以外圍上的點(diǎn)的中子通量密度均設(shè)為0。上述〔3〕和〔4〕式便是求解單群中子擴(kuò)散的差分方程組，其系數(shù)，，…，等均可以事先計(jì)算求得，對(duì)每次外迭代，其右端源項(xiàng)也是的。因此，它是一個(gè)含有的代數(shù)方程組，對(duì)于二維問題，其系數(shù)矩陣是一個(gè)五隊(duì)角的矩陣，可以很方便的應(yīng)用通常的線性方程組的求解方法進(jìn)行求解。2.2源迭代過程根據(jù)連續(xù)性的假設(shè)條件，堆芯的雙群方程與放射層得雙群方程可以歸結(jié)為下面兩個(gè)方程：快群中子擴(kuò)散方程：〔5〕熱群中子擴(kuò)散方程：〔6〕設(shè)〔7〕式〔5〕是一個(gè)關(guān)于的齊次方程組。式〔5〕除了中子通量密度外，和也是未知參數(shù)，并且這兩個(gè)參數(shù)又與中子通量密度有關(guān)。因此我們假定一個(gè)初始的裂變?cè)捶植迹J(rèn)為在整個(gè)堆芯都等于1或某個(gè)給定常數(shù)，并猜想一個(gè)初始的值，同時(shí)把作為初始迭代源項(xiàng)代入到方程〔5〕的右端中。這樣，方程〔1〕的右端便是一個(gè)項(xiàng)，方程〔1〕便成為一個(gè)普通的非齊次方程組，它可以用通常的數(shù)值方法求解。假定由它求出的中子通量密度分布為，把代入方程〔6〕，由方程〔6〕數(shù)值求解求出，將求得的、帶入〔7〕式便可求得第二次迭代裂變中子源分布，并由它可以求得有效增殖因數(shù)的新的估計(jì)值和第二代迭代源項(xiàng)來，……。依此類推，逐次的迭代下去，對(duì)于第n次迭代計(jì)算有〔8〕式中：〔9〕根據(jù)的物理意義〔它等于新生一代的裂變總中子數(shù)與上一代源中子總數(shù)之比〕可以由下式算出的估計(jì)值〔10〕方程〔3〕，〔4〕和〔5〕便是源迭代法的計(jì)算公式。從物理上看，上面所確定的迭代過程相當(dāng)于將反響堆內(nèi)的中子數(shù)人為地區(qū)分為不同的“代”，每一“代”中子有裂變形成之時(shí)作為起始。每一次源迭代那么相當(dāng)于老一代源中子通過擴(kuò)散、慢化、俘獲、裂變，最后又形成新一代源中子的過程。式〔5〕中的顯然反映了中子每經(jīng)過一代的倍增。因而可以預(yù)期，經(jīng)過了足夠的循環(huán)代數(shù)后，任意的初始中子源必將在堆內(nèi)形成一個(gè)穩(wěn)定的最終中子源分布。多群擴(kuò)散方程的特征值問題在數(shù)學(xué)上也已被詳細(xì)研究過，并證明了上述迭代過程的收斂性。2.3內(nèi)迭代過程對(duì)于一次內(nèi)迭代過程，假定一個(gè)初始的裂變?cè)捶植技耙粋€(gè)初始的數(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)的中子通量密度認(rèn)為為0，不必計(jì)算?！?〕反射層的中子通量密度反射層快群〔第一群〕和熱群〔第二群〕的擴(kuò)散方程為：2.4有效增殖因子的計(jì)算將上述芯部計(jì)算過程得到的快中子和熱中子的通量密度和同時(shí)代入下式根據(jù)的計(jì)算式得分別將本次中子源和中子增殖因子與上一次的中子源和中子增殖因子進(jìn)行比擬，觀察是否滿足以下收斂準(zhǔn)那么特征值收斂準(zhǔn)那么:中子通量密度分布收斂準(zhǔn)那么：在本設(shè)計(jì)中，和分別取為和。經(jīng)過計(jì)算，假設(shè)上述兩個(gè)收斂準(zhǔn)那么同時(shí)得到滿足，那么此次計(jì)算得到的中子通量密度分布和有效增殖因子為所求的解；反之，那么將此次計(jì)算得到的中子通量密度分布和有效增殖因子作為新的裂變?cè)捶植己托碌臄?shù)值再重復(fù)上述計(jì)算步驟。2.5建模由于問題具有對(duì)稱性，因此，故計(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)帶有反射層反響堆內(nèi)中子源的分布由中心向四周逐漸減小，在反射層內(nèi)為零。但是其在靠近反射層附近出現(xiàn)一個(gè)小的峰值，尤其是在邊角這個(gè)峰值最大，出現(xiàn)這種現(xiàn)象的原因可能是快中子進(jìn)入反射層后被慢化為慢中子回到堆芯，在靠近反射層的地方產(chǎn)生較大的反響界面，產(chǎn)生更多的熱中子。(2)如圖7、8、9、10所示，熱中子通量密度曲線在堆芯與反射層的交界處后有一個(gè)明顯的峰值，而快中子通量密度曲線有一個(gè)接近水平的平緩變化的段。主要由于反射層對(duì)熱中子的吸收比堆芯弱，但對(duì)快中子的慢化比堆芯強(qiáng)，由堆芯泄泄漏的快中子被反射層慢化成熱中子的緣故。2總結(jié)通過此次核反響堆物理分析課程設(shè)計(jì)，我收獲了很多東西。首先，我對(duì)課本上很多當(dāng)初不懂的知識(shí)有了更深刻的認(rèn)識(shí)；其次，對(duì)于很多以前課本上的結(jié)論公式，知道了其來源和用處；然后，也對(duì)以前學(xué)過的MATLAB軟件的編程和數(shù)據(jù)處理有了更好的掌握，也會(huì)通過ORIGIN軟件畫圖，節(jié)省了很多的計(jì)算過程；并且通過課程設(shè)計(jì)了解了這門的處理問題的思想，同時(shí)也敬佩做研究人的那種處理抽象問題的能力。最后，感謝老師的辛苦教導(dǎo)。四、參考文獻(xiàn)[1]謝仲生主編.《核反響堆物理分析》.修訂本.西安交通大學(xué)出版社.原子能出版社.2012.7[2]劉衛(wèi)國主編.《MATLAB程序設(shè)計(jì)與運(yùn)用》.第二版.高等教育出版社.2006.7[3]劉衛(wèi)國主編.《MATLAB程序設(shè)計(jì)教程》.第二版.中國水利水電出版社.2010.2附錄〔設(shè)計(jì)流程圖、程序、表格、數(shù)據(jù)等〕幾何或成分調(diào)整初始，n=0幾何或成分調(diào)整初始，n=0群常數(shù)群常數(shù)n=n+1n=n+1臨界計(jì)算臨界計(jì)算g=0g=0初始；g=g+1初始；g=g+1臨界搜索臨界搜索差分方程的迭代求解外迭代差分方程的迭代求解外迭代內(nèi)迭代內(nèi)迭代否g=G否g=G是是計(jì)算計(jì)算否否否是否是是是停機(jī)停機(jī)圖4〔二〕三維圖像圖41/4帶反射層的堆芯內(nèi)中子源的空間分布圖51/4帶反射層的堆芯內(nèi)中子源沿X軸方向分布圖61/4帶反射層的堆芯內(nèi)中子源沿Y軸方向分布圖71/4帶反射層的堆芯內(nèi)中子通量密沿X軸方向分布圖81/4帶反射層的堆芯內(nèi)中子通量密沿Y軸方向分布圖91/4帶反射層的堆芯內(nèi)快中子通量密度的空間分布圖101/4帶反射層的堆芯內(nèi)熱中子通量密度的空間分布〔三〕程序%定義堆芯參數(shù)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%定義反射層參數(shù)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%設(shè)置迭代初值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);%定義收斂準(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)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