本發(fā)明涉及核反應堆堆芯設計和反應堆物理計算領域，具體涉及一種改進反應堆堆芯三維中子通量密度精細分布的方法。
背景技術：
：反應堆堆芯設計和運行需要準確快速地計算出反應堆及相關的設備內三維中子通量密度分布的情況。目前廣泛采用的傳統(tǒng)反應堆物理分析計算方法，隨著反應堆堆芯設計越來越復雜，安全要求越來越高，逐漸不能滿足工程計算要求。所謂傳統(tǒng)反應堆物理分析計算方法又稱為“兩步法”，第一步是在全反射邊界條件下對各種非均勻組件進行多群中子輸運計算，得到組件的少群等效均勻化群常數(shù)及不連續(xù)因子等物理量；第二步，是根據(jù)前一步生成的均勻化參數(shù)，采用粗網(wǎng)節(jié)塊方法對堆芯進行少群中子擴散或輸運計算，附加功率重構等計算，便可實現(xiàn)對三維堆芯功率分布的模擬?！皟刹椒ā庇捎谥豢紤]有限孤立組件，導致真實堆芯布置中組件之間的相互影響無法充分考慮，以及組件的燃耗歷史效應等重要因素，其理論已經(jīng)限制了對計算精度的追求?！耙徊健比研痉蔷鶆蛑凶虞斶\計算可以從根本上解決上述問題。目前的計算條件下，直接對大型核反應堆堆芯進行三維中子輸運計算還不現(xiàn)實。目前的二維/一維耦合方法被提出，核心思想是耦合二維、一維計算來替代三維計算。但是目前該方法中關于泄漏項的處理有些簡單，需要做一些改進。技術實現(xiàn)要素：為了改進反應堆三維中子通量密度精細分布獲取精度，本發(fā)明的目的在于提供一種改進反應堆堆芯三維中子通量密度精細分布的方法，本發(fā)明方法將二維/一維耦合計算中徑向泄漏項的平近似改為正交多項式展開，軸向泄漏項用通量計算出來的分布來代替平近似，能夠很好的改善反應堆堆芯三維中子通量密度的精細分布的計算精度，給反應堆堆芯設計和安全提供可靠的信息。為了實現(xiàn)上述目的，本發(fā)明采用如下技術方案：一種改進反應堆堆芯三維中子通量密度精細分布的方法，步驟如下：步驟1：對于所涉及的反應堆堆芯進行幾何建模，劃分計算區(qū)域，離散角度空間，生成特征線信息，指定各計算區(qū)域的材料，獲取材料的宏觀截面參數(shù)，為計算區(qū)域的中子通量密度通量、反應堆邊界條件、特征值設置初始值；a)先將反應堆堆芯按照軸向劃分成若干層，每一層劃分成若干計算區(qū)域；b)根據(jù)計算需求對每一層計算區(qū)域的三維角度空間進行離散；c)在每一層內生成特征線信息；d)根據(jù)所涉及反應堆堆芯材料信息，讀取各層中每一個計算細區(qū)的宏觀截面參數(shù)；e)對于二維、一維計算區(qū)域的中子通量密度，反應堆邊界條件，特征值賦初值；步驟2：更新徑向泄漏項，根據(jù)步驟1中所得幾何尺寸、材料信息計算一維離散縱標差分方程，求得每個柵元的每一群中子通量密度，求得軸向泄漏項；在二維/一維耦合算法中，徑向泄漏項是由二維求解器計算提供的，因此在求解徑向一維方程時，認為該項是已知源項；一般計算中，每一層的寬度為10cm-20cm，二維提供的徑向泄漏項在每一層內是常數(shù)，引入這樣一個平泄漏近似對于計算精度是有影響的；所以給一維計算提供的徑向泄漏項是對二維提供的平泄漏做了勒讓德多項式多項式展開，得到徑向泄漏沿軸向的分布，對于最后的計算結果有明顯改善；步驟3：更新軸向泄漏項，根據(jù)步驟1中所得幾何尺寸、材料信息、特征線信息對軸向的每一層進行二維輸運計算，迭代求解每一層的二維中子通量密度，計算柵元均勻化截面，計算徑向泄漏項；在計算過程中，利用平源區(qū)通量計算得到一個柵元內軸向泄漏項的空間分布；在二維/一維耦合算法中，軸向泄漏項是由一維軸向求解計算提供的，因此在求解徑向二維方程時，認為該項是已知源項；步驟2中一維計算的區(qū)域選擇的是一個柵元，所以上述軸向泄漏項在一個柵元中是一個平分布；真實情況是軸向泄露在一個柵元中也存在分布；利用上一個迭代中二維計算得到的通量分布計算得到柵元內軸向泄漏項的空間分布，對于整個二維一維耦合計算有精度上的改善；用平源區(qū)上下表面的標通量之差計算得到軸向泄漏項在柵元內的分布；平源區(qū)的上下表面的標通量則是通過將通量沿軸向展開得到的；同樣采用二階的勒讓德多項式，計算得到上下表面的標通量；最終，柵元內每個平源區(qū)的差值組成的分布乘以柵元的軸向泄漏項得到柵元內軸向泄漏項的差值；從平分布近似到有分布的精確描述，提高了計算的精度；步驟4：判斷特征值和三維中子通量密度是否收斂，如果不收斂，則轉到步驟2繼續(xù)迭代，并且計算中采用上一次迭代中求得的中子通量密度和特征值，直至中子通量密度和特征值收斂，就能夠得到三維中子通量密度精細分布。與現(xiàn)有技術相比，本發(fā)明有如下突出優(yōu)點：1.一維計算中對徑向泄漏項進行多項式展開，使用帶分布的徑向泄漏項，比平泄露近似精度高。2.二維計算中求得軸向泄漏項的二維分布，使得計算更接近真實情況。附圖說明圖1為二維/一維耦合計算流程圖。圖2為二維/一維堆芯劃分示意圖。圖3為反應堆堆芯幾何布置示意圖。圖4為計算區(qū)域角度空間離散示意圖。圖5為計算區(qū)域特定方向特征線示意圖。圖6為徑向泄漏項多項式展開示意圖。圖7為軸向泄漏項多項式展開示意圖。圖8為三維非均勻壓水堆堆芯相對棒功率分布圖。具體實施方式本發(fā)明將二維/一維耦合計算中徑向泄漏項的平近似改為正交多項式展開，軸向泄漏項用通量計算出來的分布來代替平近似，它能夠很好的改善反應堆堆芯三維中子通量密度的精細分布的計算精度，給反應堆堆芯設計和安全提供可靠的信息，具體實施方式如下所示。圖1所示為二維/一維耦合方法的總體流程圖。步驟1：對于所涉及的反應堆堆芯進行幾何建模，劃分計算區(qū)域，離散角度空間，生成特征線信息，指定各計算區(qū)域的材料，獲取材料的宏觀截面參數(shù)，為計算區(qū)域的中子通量密度通量、反應堆邊界條件、特征值設置初始值，圖2為二維/一維堆芯劃分示意圖；f)先將反應堆堆芯按照軸向劃分成若干層，每一層劃分成若干計算區(qū)域，圖3為壓水堆堆芯幾何布置示意圖；g)根據(jù)計算需求對每一層計算區(qū)域的三維角度空間進行離散，如圖4所示的角度空間離散，圖中用了9個離散方向等效八分之一角度空間，θ代表極角，代表輻角；h)在每一層內生成特征線信息，如圖5所示的某區(qū)域中某一個方向的所有特征線，圖中實線為特征線，虛線為兩條相鄰特征線的間距中心線，是出射中子角通量，是入射中子角通量密度，si,k是特征線的長度，δAk是特征線線寬；i)根據(jù)所涉及反應堆堆芯材料信息，讀取各層中每一個計算細區(qū)的宏觀截面參數(shù)；j)對于二維、一維計算區(qū)域的中子通量密度，反應堆邊界條件，特征值賦初值；步驟2：更新徑向泄漏項，根據(jù)步驟1中所得幾何尺寸、材料信息計算一維離散縱標差分方程，求得每個柵元的每一群中子通量密度，求得軸向泄漏項。在二維/一維耦合算法中，徑向泄漏項是由二維求解器計算提供的，因此在求解徑向一維方程時，認為該項是已知源項。一般計算中，每一層的寬度在10cm-20cm左右，二維提供的徑向泄漏項在每一層內是常數(shù)，引入這樣一個平泄漏近似對于計算精度是有影響的。所以本發(fā)明中在給一維計算提供的徑向泄漏項是對二維提供的平泄漏做了勒讓德多項式多項式展開，得到徑向泄漏沿軸向的分布，對于最后的計算結果有明顯改善；步驟3：更新軸向泄漏項，根據(jù)步驟1中所得幾何尺寸、材料信息、特征線信息對軸向的每一層進行二維輸運計算，迭代求解每一層的二維中子通量密度，計算柵元均勻化截面，計算徑向泄漏項；在計算過程中，利用平源區(qū)通量計算得到一個柵元內軸向泄漏項的空間分布；在二維/一維耦合算法中，軸向泄漏項是由一維軸向求解計算提供的，因此在求解徑向二維方程時，認為該項是已知源項。步驟2中一維計算的區(qū)域選擇的是一個柵元，所以上述軸向泄漏項在一個柵元中是一個平分布；真實情況是軸向泄露在一個柵元中也存在分布。本發(fā)明利用上一個迭代中二維計算得到的通量分布計算得到柵元內軸向泄漏項的空間分布，對于整個二維一維耦合計算有精度上的改善。本發(fā)明用平源區(qū)上下表面的標通量之差計算得到軸向泄漏項在柵元內的分布。平源區(qū)的上下表面的標通量則是通過將通量沿軸向展開得到的。同樣采用二階的勒讓德多項式，計算得到上下表面的標通量。最終，柵元內每個平源區(qū)的差值組成的分布乘以柵元的軸向泄漏項得到柵元內軸向泄漏項的差值。從平分布近似到有分布的精確描述，提高了計算的精度；步驟4：判斷特征值和三維中子通量密度是否收斂，如果不收斂，則轉到步驟2繼續(xù)迭代，并且計算中采用上一次迭代中求得的中子通量密度和特征值，直至中子通量密度和特征值收斂，就可以得到三維中子通量密度精細分布。步驟1中二維/一維耦合算法中的幾何描述、計算區(qū)域劃分包括軸向分層、徑向劃分子區(qū)、子區(qū)內劃分計算區(qū)域、所有計算區(qū)域的角度空間離散、子區(qū)內特征線生成均需要根據(jù)不同的反應堆堆芯和不同的計算條件選擇合適的方案。一般壓水堆堆芯軸向分為10到15層，徑向每一層分為9到25個子區(qū)，每個子區(qū)中幾百個計算區(qū)域不等，通常使用48到80個角度方向去離散角度空間，角度空間的離散方案確定了特征線的方向，特征線線寬一般采用0.001cm到0.05cm。步驟2中求解離散縱標差分方程，徑向泄漏項的展開，推導過程如下所示。三維直角坐標系中，穩(wěn)態(tài)多群中子輸運方程如公式(1)所示，Ω·▽ψg(r,Ω)+Σt,g(r)ψg(r,Ω)＝Qg(r,Ω),g＝1,...,G公式(1)式中：ψg——第g能群中子角通量密度g——能群標識；G——能群總數(shù)；Q——中子輸運方程源項。以特征值問題為例，第g能群源項的具體形式為：式中：φg——第g能群中子通量密度——第g能群的裂變源項；——第g能群的散射源項。經(jīng)過角度離散后，在第g群、第m離散方向上的形式如公式(2)所示。根據(jù)反應堆堆芯進行的幾何劃分，對每個柵元進行積分，并除以柵元面積，最終整理得：ξmΔx[ψg,m,x+p(z)-ψg,m,x-p(z)]+ηmΔy[ψg,m,y+p(z)-ψg,m,y-p(z)]+μmdψg,mp(z)dz+Σt,g,p(z)ψg,mp(z)=Qgp(z)]]>在該方程中，徑向泄漏項定義為：TLg,m,pRadial(z)=ξmΔx[ψg,m,x+p(z)-ψg,m,x-p(z)]+ηmΔy[ψg,m,y+p(z)-ψg,m,y-p(z)]]]>其物理意義是第k層沿軸向在離散方向m上的泄漏。在二維/一維耦合算法中，徑向泄漏項是由二維求解器計算提供的，因此在求解徑向一維方程時，認為該項是已知源項。一般計算中，每一層的寬度在10cm-20cm左右，二維提供的徑向泄漏項在每一層內是常數(shù)，引入這樣一個平泄漏近似對于計算精度是有影響的。所以本發(fā)明中在給一維計算提供的徑向泄漏項是對二維提供的平泄露做了多項式展開，得到徑向泄漏沿軸向的分布，對于最后的計算結果有明顯改善。由于徑向泄漏項在每一層是平分布，為了得到徑向泄漏項精細的軸向分布，本發(fā)明將徑向泄漏項沿軸向展開，采用勒讓德多項式，表達式如下所示。Ln(ξ)=Σi=02linPi(ξ)]]>whereP0＝1,P1＝ξ,P2＝3ξ2-1/4l0k=TLg,kRadial]]>l1k=Δzkdzk-1(TLg,k+1Radial-TLg,kRadial)+zk+1(TLg,kRadial-TLg,k-1Radial)]]>l2k=(Δzk)2d(Δzk-1+Δzk)(TLg,k+1Radial-TLg,kRadial)-(Δzk+1+Δzk)(TLg,kRadial-TLg,k-1Radial)]]>wherezi＝(Δzk+Δzi)(Δzk+2Δzi)d＝(Δzk-1+Δzk)(Δzk+Δzk+1)(Δzk-1+Δzk+Δzk+1)其中，L表示某個方向的泄漏項，經(jīng)驗表明采取勒讓德二階近似就可以取得很好的結果；0階、1階、2階的系數(shù)分別如上述所示。得到這三個系數(shù)后，就可以輕松計算得到每個細網(wǎng)上的徑向泄漏值。圖6為徑向泄漏項展開前后示意圖，展開之前在每層內徑向泄漏項是平分布，展開之后可以計算得到每個平源區(qū)網(wǎng)格上的徑向泄漏項。一維Sn差分方程求解步驟：1)用勒讓德多項式沿軸向展開徑向泄漏項，得到每個差分網(wǎng)格的徑向泄漏項；2)根據(jù)上次計算得到的Keff和通量更新右端源項；3)根據(jù)差分格式一次計算每個網(wǎng)格的、每一個方向、每一個能群的角通量；4)判斷步驟3)計算得到的通量是否收斂，不收斂則繼續(xù)進行2)-3)的迭代計算，若通量收斂則結束計算。收斂準則如下：||φ‾gk-φ‾gk-1φ‾gk||≤ϵ1]]>ε1為很小的正數(shù)，也就是收斂條件。步驟3中根據(jù)步驟1中反應堆堆芯進行的幾何劃分，對軸向的每一層進行積分，并處以網(wǎng)格寬度，最終整理得：在該方程中，定義軸向泄漏項為：其物理意義是第k層沿軸向在離散方向m上的泄漏。在二維/一維耦合算法中，軸向泄漏項是由一維軸向求解計算提供的，因此在求解徑向二維方程時，認為該項是已知源項。步驟2中一維計算的區(qū)域選擇的是一個柵元，所以上述軸向泄漏項在一個柵元中是一個平分布。真實情況是軸向泄露在一個柵元中也存在分布。本發(fā)明利用上一個迭代中二維計算得到的通量分布計算得到柵元內軸向泄漏項的空間分布，對于整個二維一維耦合計算有精度上的改善。TLg,m,kAxial(x,y)=μmΔzk[ψg,m,k+1/2(x,y)-ψg,m,k-1/2(x,y)]]]>TL_shapekAxial(x,y)=μmΔzk[φk+1/2(x,y)-φk-1/2(x,y)]]]>φ(x,y,ξ)=Σi=02linPi(ξ)]]>其中,P0＝1,P1＝ξ,P2＝3ξ2-1/4φ0k=φk]]>φ1k=Δzkdzk-1(φk+1-φk)+zk+1(φk-φk-1)]]>φ2k=(Δzk)2d(Δzk-1+Δzk)(φk+1-φk)-(Δzk+1+Δzk)(φk-φk-1)]]>其中,zi＝(Δzk+Δzi)(Δzk+2Δzi)d＝(Δzk-1+Δzk)(Δzk+Δzk+1)(Δzk-1+Δzk+Δzk+1)φk+1/2(x,y)＝φ(x,y,1/2)φk-1/2(x,y)＝φ(x,y,-1/2)上述表達式中，本發(fā)明用平源區(qū)上下表面的標通量之差計算得到軸向泄漏項在柵元內的分布。平源區(qū)的上下表面的標通量則是通過將通量沿軸向展開得到的。同樣采用二階的勒讓德多項式，計算得到上下表面的標通量。最終，柵元內每個平源區(qū)的差值組成的分布乘以柵元的軸向泄漏項得到柵元內軸向泄漏項的差值。從平分布近似到有分布的精確描述，提高了計算的精度。圖7為軸向泄漏項展開前后示意圖，重構之前在每個柵元內軸向泄漏項是平分布，展開之后可以計算得到每個平源區(qū)網(wǎng)格上的軸向泄漏項。MOC方法是將公式(2)轉換成沿著某條特征線的微分方程，如下所示：認為Σt,g,Qg在某個區(qū)域內是常數(shù)，則公式(6)是一個常微分方程，可以有解析解。步驟1中二維的計算區(qū)域在合適的范圍內，MOC方法可以得到非常精確的輸運解。沿著某一段特征線的平均中子角通量密度表達式如下所示：對同一方向所有特征線的平均中子角通量密度進行體積加權平均得到該計算區(qū)域這個方向的平均中子角通量密度，對該計算區(qū)域所有方向的中子角通量密度進行加權平均就可以得到該計算區(qū)域的平均中子通量密度。由于公式(7)的右端含有源項，該源項需要通過中子通量密度進行求解，所以一般采用迭代計算求解MOC方程。二維MOC方程求解流程如下：1)更新軸向泄漏項，該泄漏項由一維計算提供，然后根據(jù)上述重構方式，得到軸向泄漏項在一個柵元內的空間分布；2)根據(jù)上一步計算得到的Keff、入射角通量和標通量更新右端源項；3)根據(jù)MOC方程解的表達形式迭代計算每個平源區(qū)的標通量；4)判斷平源區(qū)標通量是否收斂，若不收斂，繼續(xù)迭代求解，若收斂更新徑向泄漏項和柵元均勻化截面；判斷所有能群是否都已經(jīng)計算，對各個能群進行迭代，直到所有區(qū)域通量收斂。步驟4判斷特征值和三維通量是否收斂，如果不收斂，則轉到步驟3繼續(xù)迭代，并且計算中采用上一次迭代中步驟2-步驟3求得的中子通量密度和特征值，反復迭代直至中子通量密度和特征值收斂，就可以得到三維中子通量精細分布。所述收斂條件為：||φgi,(l)-φgi,(l-1)φgi,(l)||≤ϵ1,||keff(l)-keff(l-1)keff(l)||≤ϵ2]]>其中：為第l次迭代第i區(qū)第g群通量為第l次迭代特征值ε1，ε2為兩個很小的正數(shù)，也就是收斂條件。大量數(shù)值驗證結果顯示，本發(fā)明具有可靠的精度、很好的效率、對于平泄漏近似的改善非常有效，適應工程實際中的計算要求。圖8是典型的壓水堆堆芯三維功率分布，其中外圍黑色部分是反射層沒有功率，內部的控制棒區(qū)域也沒有功率。當前第1頁1 2 3