輸運(yùn)問(wèn)題的蒙特卡羅模擬課件

北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所

JMCT團(tuán)隊(duì)

iapcm.ac鄧力粒子輸運(yùn)MC方法及軟件發(fā)展現(xiàn)狀科大，2019年10月31日北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所

JMCT團(tuán)隊(duì)

iapcm.報(bào)告內(nèi)容5.國(guó)內(nèi)現(xiàn)狀4.MC方法的應(yīng)用3.發(fā)展與挑戰(zhàn)2.粒子輸運(yùn)計(jì)算1.MC方法回顧

6.總結(jié)與展望2022/10/152報(bào)告內(nèi)容5.國(guó)內(nèi)現(xiàn)狀4.MC方法的應(yīng)用3.發(fā)展與一.MC方法回顧(1)

S.UlamiscreditedastheinventorofMonteCarlomethodin1940s,whichsolvesmathematicalproblemsusingstatisticalsampling.ThealgorithmbyMetropolis(andARosenbluth,MRosenbluth,ATellerandETeller,1953)hasbeencitedasamongthetop10algorithmshavingthe"greatestinfluenceonthedevelopmentandpracticeofscienceandengineeringinthe20thcentury."S.Ulam

Metropolis

2022/10/153一.MC方法回顧(1)S.Ulamiscredite一.MC方法回顧(2)

JohnvonNeumann(1903.12.28-1957.2.8)N.Kolmogorov

(1903.4.12-1987.10.20)

2022/10/154一.MC方法回顧(2)JohnvonNeumannN一.MC方法回顧(3)

Monte-Carlo,Monaco1944年VonNeumann等人把他們研制的第一個(gè)隨機(jī)模擬中子鏈?zhǔn)椒磻?yīng)的程序用摩納哥著名賭城“MonteCarlo”命名。

2022/10/155一.MC方法回顧(3)Monte-Carlo,Monac研究背景MonteCarlo（MC）粒子輸運(yùn)計(jì)算優(yōu)勢(shì)：適應(yīng)實(shí)際幾何形狀精確的結(jié)果面臨的問(wèn)題：復(fù)雜幾何處理計(jì)算量巨大精細(xì)模擬，存儲(chǔ)挑戰(zhàn)2022/10/156研究背景MonteCarlo（MC）粒子輸運(yùn)計(jì)算優(yōu)勢(shì)：20一.MC方法回顧(4)

核裝置具有系統(tǒng)小、結(jié)構(gòu)復(fù)雜、含輕重介質(zhì)、尺度小、系統(tǒng)反應(yīng)時(shí)間短(s量級(jí))、高溫、高壓、沖擊波等特點(diǎn)。數(shù)值模擬過(guò)程涉及流體力學(xué)與中子(光子)輸運(yùn)耦合計(jì)算，其中中子計(jì)算量占總計(jì)算量的98%以上，輸運(yùn)計(jì)算精度對(duì)整個(gè)核系統(tǒng)放能起著決定作用。在核科學(xué)領(lǐng)域，蒙特卡羅方法仍然是最主要的方法。美國(guó)能源部宣稱，自從高性能計(jì)算機(jī)應(yīng)用以來(lái)，蒙特卡羅模擬花費(fèi)了一半以上的機(jī)時(shí)。在我所情況也是如此，定態(tài)、動(dòng)態(tài)MC程序在武器物理研究中發(fā)揮積極而重要的作用。過(guò)去許多無(wú)法看到的物理圖象，如今在大規(guī)模并行計(jì)算條件下，可以看得比較清楚。核武器的發(fā)展帶動(dòng)了計(jì)算機(jī)的發(fā)展，也帶動(dòng)了我所計(jì)算方法及軟件的發(fā)展。2022/10/157一.MC方法回顧(4)核裝置具有系統(tǒng)小、結(jié)構(gòu)一.MC方法回顧(5)

MonteCarle方法模擬能力

精細(xì)幾何描述；使用連續(xù)點(diǎn)截面，基于連續(xù)能量的粒子物理；精細(xì)的角分布；大規(guī)模并行計(jì)算

反應(yīng)堆物理分析

1960年代－keff計(jì)算，理論為主；1970年代－keff計(jì)算，考慮組件功率計(jì)算；1980年代－keff計(jì)算，二維堆芯模擬；1990年代－keff計(jì)算，三維全堆芯模擬；2000年代－keff計(jì)算，精細(xì)三維全堆芯模擬、損耗及參數(shù)設(shè)計(jì)、連續(xù)點(diǎn)截面MC中子輸運(yùn)與燃耗耦合計(jì)算程序模擬全堆。

2022/10/158一.MC方法回顧(5)MonteCarle方法模擬能力反應(yīng)堆數(shù)值模擬發(fā)展歷程2022/10/159反應(yīng)堆數(shù)值模擬發(fā)展歷程2022/10/119一.MC方法回顧(6)

MCNP程序是由美國(guó)LosAlamos國(guó)家實(shí)驗(yàn)室研制的一個(gè)大型、多功能、多粒子輸運(yùn)蒙特卡羅程序，可計(jì)算任意三維復(fù)雜幾何系統(tǒng)中子、光子、電子或中子-光子-電子耦合輸運(yùn)問(wèn)題，還可計(jì)算臨界系統(tǒng)的本征值問(wèn)題。MCNP程序具有超強(qiáng)的幾何處理能力，幾何系統(tǒng)由幾何塊(cell)組成，而幾何塊的界面由平面、二次曲面及特殊的四次橢圓環(huán)曲面組成。幾何塊中的材料由包括同位素在內(nèi)的多種核組成，使用精確點(diǎn)截面參數(shù)，對(duì)特定的評(píng)價(jià)庫(kù)(ENDF/B-IV、V、VI庫(kù)或ENDL851庫(kù))，考慮了該庫(kù)給出的所有中子反應(yīng)類(lèi)型。在截面數(shù)據(jù)文件中收集了多種評(píng)價(jià)庫(kù)的數(shù)據(jù)。對(duì)熱中子還配備了相應(yīng)的截面數(shù)據(jù)，可按自由氣體?；騍(,)模處理。對(duì)光子考慮了相干和非相干散射，處理了光電吸收后可能的熒光發(fā)射和電子對(duì)產(chǎn)生。[1]MCNP–AGeneralMonteCarloN-ParticleTransportCode,Version5,LA-UR-03-1987,LANL-X-5MonteCarloTeam2022/10/1510一.MC方法回顧(6)MCNP程序是由美國(guó)一.MC方法回顧(7)2022/10/1511一.MC方法回顧(7)2022/10/1111

考慮溫度效應(yīng)的動(dòng)力學(xué)模擬一.MC方法回顧(8)2022/10/1512考慮溫度效應(yīng)的動(dòng)力學(xué)模擬一.MC方法回顧(8)2022/MCNP存在的不足計(jì)算方法陳舊，計(jì)算效率不高，可擴(kuò)展性差模擬的粒子數(shù)少<232容許的最大幾何塊105、幾何面106次臨界本征值計(jì)算不完善多群計(jì)算存在諸多問(wèn)題并行效率和可擴(kuò)展性不佳MCNP6S(,)有較多改進(jìn)，核素?cái)?shù)增加到20個(gè)本征值計(jì)算有所改進(jìn)，考慮了緩發(fā)幾何塊擴(kuò)充到109隨機(jī)數(shù)周期擴(kuò)到264敏感性分析[2]ForrestBrownetal.,AdvancedMonteCarloforReactorPhysicsCoreAnalysis,LA-UR-12-20397,PHYSOR20192022/10/1513MCNP存在的不足計(jì)算方法陳舊，計(jì)算效率不高，可擴(kuò)展性差MC一.MC方法回顧(9)

計(jì)算方法

1970年代－三維組合幾何、反應(yīng)堆堆芯重復(fù)結(jié)構(gòu)描述，粒子徑跡計(jì)算等；[3]Carter&Cashwell,“Particle-TransportSimulationwiththeMonteCarloMethod”,ERDACriticalReviewSeries,1975,TID-26607.

MC基本原理及誤差理論；各種抽樣方法；降低方差技巧；通量估計(jì)方法。[4]格拉斯頓著，黃祖恰、千里譯，”反應(yīng)堆理論”，原子能出版社，1974.2022/10/1514一.MC方法回顧(9)計(jì)算方法1970年代－三維組合一.MC方法回顧(10)

應(yīng)用領(lǐng)域

各種粒子(中子、光子、電子、粒子、質(zhì)子)輸運(yùn)問(wèn)題模擬；

反應(yīng)堆臨界安全分析、燃料循環(huán)；反應(yīng)堆動(dòng)力學(xué)；輻射輸運(yùn)；核探測(cè)；核醫(yī)學(xué)；

其它領(lǐng)域

統(tǒng)計(jì)物理、分子動(dòng)力學(xué)、量子力學(xué)、生物工程、金融、信息和航天等2022/10/1515一.MC方法回顧(10)應(yīng)用領(lǐng)域各種粒子(中子、光子、()10.682720.954530.997340.9999

MC方法誤差置信度估計(jì)量置信度應(yīng)用xi=F-1(ξi),i=1,2,,N二.粒子輸運(yùn)計(jì)算(1)2022/10/1516()10.682720.954530.997340.9Namema(multiplier)cperiodANSIC231110351524512345231Park-Miller231-1168070231-2drand48()2482521490391711248Hayes64-bit26463641362238467930051264二.粒子輸運(yùn)計(jì)算(3)Oneoftheearliestandfastestalgorithm:

xn+1=(a

xn+c)modm

n+1=xn+1/m

where0≤x<m,misthemodulus,aismultiplier,cisincrement.Allofthemareintegers.Choiceofa,c,mmustbedonewithspecialcare.2022/10/1517Namema(multiplier)cperiodANSI二.粒子輸運(yùn)計(jì)算(4)Boltzmann方程求解：微分-積分形式其中：(2.1)(2.2)2022/10/1518二.粒子輸運(yùn)計(jì)算(4)Boltzmann方程求解：其中：(二.粒子輸運(yùn)計(jì)算(5)中子速度與能量E滿足動(dòng)能方程(不考慮相對(duì)論效應(yīng))(2.3)2022/10/1519二.粒子輸運(yùn)計(jì)算(5)中子速度與能量E滿足動(dòng)能方程(不考慮定態(tài)：定常非定常動(dòng)態(tài)：(時(shí)間行為來(lái)自源)輸運(yùn)問(wèn)題分類(lèi)：系統(tǒng)狀態(tài)隨時(shí)間變化！系統(tǒng)狀態(tài)不隨時(shí)間變化！二.粒子輸運(yùn)計(jì)算(6)2022/10/1520定態(tài)：定常非定常動(dòng)態(tài)：(時(shí)間行為來(lái)自源)輸運(yùn)問(wèn)題分類(lèi)：系統(tǒng)狀二.粒子輸運(yùn)計(jì)算(7)Boltzmann方程等價(jià)的積分方程（MC）(2.4)定義輸運(yùn)核(或稱遷移核)和碰撞核：(2.5)(2.6)其中：(2.7)2022/10/1521二.粒子輸運(yùn)計(jì)算(7)(2.4)定義輸運(yùn)核(或稱遷移核)和二.粒子輸運(yùn)計(jì)算(8)此為特征線方程。特征線為：(2.8)(2.9)角通量的全空間形式2022/10/1522二.粒子輸運(yùn)計(jì)算(8)(2.8)(2.9)2022/10/

(1)從MC隨機(jī)模擬過(guò)程抽初始狀態(tài)：，初始權(quán)(2)從確定下一個(gè)碰撞點(diǎn)及到達(dá)時(shí)間，則中子穿出系統(tǒng)，歷史結(jié)束；若，則中子歷史結(jié)束；若二.粒子輸運(yùn)計(jì)算(9)2022/10/1523

(1)從MC隨機(jī)模擬過(guò)程抽初始狀態(tài)：，初始權(quán)(2)從

(4)確定新的中子權(quán)(3)確定碰撞核及反應(yīng)類(lèi)型

從確定碰撞核；從確定反應(yīng)類(lèi)型；

從確定新的中子能量及方向。若或，則中子歷史結(jié)束；否則回到(2).二.粒子輸運(yùn)計(jì)算(10)2022/10/1524

(4)確定新的中子權(quán)(3)確定碰撞核及反應(yīng)類(lèi)型從(6)給出統(tǒng)計(jì)平均結(jié)果(5)記錄貢獻(xiàn)

發(fā)射密度：通量：誤差：

發(fā)射密度：通量：反應(yīng)率：二.粒子輸運(yùn)計(jì)算(11)2022/10/1525(6)給出統(tǒng)計(jì)平均結(jié)果(5)記錄貢獻(xiàn)發(fā)射密度：發(fā)射密確定論方法從微分-積分方程出發(fā)，通過(guò)變量離散，把單個(gè)方程變?yōu)榫€性代數(shù)方程組。由于存儲(chǔ)每個(gè)變量的數(shù)組需要較大的內(nèi)存，早期受計(jì)算機(jī)內(nèi)存限制，對(duì)角度采用擴(kuò)散近似，對(duì)能量作多群近似，對(duì)空間作粗網(wǎng)近似，對(duì)輕介質(zhì)系統(tǒng)模擬精度較差。隨著計(jì)算機(jī)速度和存儲(chǔ)的提高，方向數(shù)和能群數(shù)進(jìn)一步增加，空間從一維發(fā)展到三維，網(wǎng)格從粗網(wǎng)格發(fā)展到細(xì)網(wǎng)格，計(jì)算方法從擴(kuò)散發(fā)展到輸運(yùn)。2.1確定論方法

二.粒子輸運(yùn)計(jì)算(12)2022/10/1526確定論方法從微分-積分方程出發(fā)，通過(guò)變量離散，把單個(gè)二.粒子輸運(yùn)計(jì)算(13)

確定論方法對(duì)一維、二維問(wèn)題是成功的，但對(duì)三維問(wèn)題，只有簡(jiǎn)單幾何問(wèn)題能夠模擬，對(duì)大變形問(wèn)題或網(wǎng)格長(zhǎng)、寬比超過(guò)3的模型，求解仍然存在困難；確定論方法計(jì)算速度快，計(jì)算精度高，深穿透問(wèn)題容易計(jì)算，適合多群伴隨輸運(yùn)方程的求解。確定論方法的近似主要來(lái)自對(duì)變量的離散，主要不足是離散過(guò)程引入了一定近似，特別當(dāng)靶核運(yùn)動(dòng)和流體運(yùn)動(dòng)對(duì)核運(yùn)動(dòng)有直接影響時(shí)，確定論方法只能近似考慮這些影響；高能部分存在射線效應(yīng)，負(fù)通量存在；無(wú)法在算出通量

的同時(shí)算出經(jīng)通量響應(yīng)量。確定論方法的優(yōu)缺點(diǎn)2022/10/1527二.粒子輸運(yùn)計(jì)算(13)確定論方法對(duì)一維、二維問(wèn)題是成功二.粒子輸運(yùn)計(jì)算(14)

MC方法從Bolzmann積分方程出發(fā)，采用源迭代求解：

2.2MC方法(2.10)其中

(2.11)每迭代一次中子就經(jīng)歷了一次輸運(yùn)和碰撞。共有4種估計(jì)：(1)碰撞估計(jì)；(2)吸收估計(jì)；(3)徑跡長(zhǎng)度估計(jì)；(4)點(diǎn)估計(jì)。

2022/10/1528二.粒子輸運(yùn)計(jì)算(14)MC方法從Bolzmann積分方程二.粒子輸運(yùn)計(jì)算(15)

包括點(diǎn)、面、體通量、各種反應(yīng)率、劑量率、反照率、能量沉積率等

缺點(diǎn)：收斂慢和誤差的概率性；深穿透問(wèn)題算不好，計(jì)算結(jié)果往往低于實(shí)驗(yàn)結(jié)果。MC方法可以在算出通量的同時(shí)算出各種通量響應(yīng)量：(2.12)優(yōu)點(diǎn)：[5]杜書(shū)華等，“輸運(yùn)問(wèn)題的計(jì)算機(jī)模擬”，湖南科學(xué)技術(shù)出版社，1989；[6]謝仲生，鄧力，“中子輸運(yùn)理論數(shù)值計(jì)算方法”，西工大出版社，2019.

2022/10/1529二.粒子輸運(yùn)計(jì)算(15)包括點(diǎn)、面、體通量、各種反應(yīng)率、劑二.粒子輸運(yùn)計(jì)算(16)

2.3MC方法與SN方法耦合計(jì)算SN

(MC)方法解伴隨方程為MC計(jì)算產(chǎn)生價(jià)值函數(shù)，指導(dǎo)MC的賭、分裂微擾、敏感性分析、eff……都需要伴隨通量(2.13)(2.14)(2.15)(2.16)2022/10/1530二.粒子輸運(yùn)計(jì)算(16)2.3MC方法與SN方法耦合計(jì)算二.粒子輸運(yùn)計(jì)算(17)

MORSE-DOMINO-DOT程序[7-9]，計(jì)算一個(gè)石油測(cè)井問(wèn)題，先用二維離散縱標(biāo)DOT程序算出中子通量在地層中的分布，以此作為MC程序MORSE的價(jià)值函數(shù)，通過(guò)分裂增加朝探測(cè)器方向飛行的粒子軌跡，計(jì)算結(jié)果表明后者較直接模擬效率提高一個(gè)量級(jí)。[7]EmmettMB.MORSE-CGA,AMonteCarloRadiationTransportCodewithArrayGeometryCapability.ORNL-6174,1985.[8]EmmetMB,BurgartCE,HoffmallTJ.DOMINO,AGeneralPurposeCodeforCouplingDiscreteOrdinatesandMonteCarloRadiationTransportCalculations.ORNL-4853,1973.[9]RoadesWA&MynattFR.TheDOTIIITwo-DimensionalDiscreteOrdinatesTransportCode.ORNL-TM-4280,1973.IiIj2022/10/1531二.粒子輸運(yùn)計(jì)算(17)MORSE-DOMINO-DOT程二.粒子輸運(yùn)計(jì)算(18)

MC與SN耦合計(jì)算把問(wèn)題區(qū)域劃分為MC適合區(qū)和SN適合區(qū)，交界面作為源區(qū)，兩種方法交替提供源項(xiàng)，為保證結(jié)果無(wú)偏，需在交界面兩邊建立適當(dāng)?shù)倪^(guò)渡區(qū)劉連焰博士發(fā)展了MC價(jià)值函數(shù)的自動(dòng)產(chǎn)生方法，在MCNP程序中應(yīng)用，取得良好的效果[10-11]MCSNS過(guò)渡區(qū)[10]LianyanLiu,Self-OptimizingMonteCarloMethodforNuclearWellLoggingSimulation[Ph.D.Thesis],NorthCarolinaStateUniv.,USA,2019.[11]LianyanLiu,R.P.Gardner,AGemetry-independentFine-Mesh-BasedMonteCarloImportanceGenerator,Nucl.Sci.Eng.125,188-195,2019.2022/10/1532二.粒子輸運(yùn)計(jì)算(18)MC與SN耦合計(jì)算MCSNS過(guò)渡二.粒子輸運(yùn)計(jì)算(19)

該方法在“堆外探測(cè)室響應(yīng)”計(jì)算中得到應(yīng)用[7]；[7]竹生東、鄧力、李樹(shù)等，堆外核儀表系統(tǒng)(RPN)的預(yù)設(shè)效驗(yàn)系數(shù)理論計(jì)算，核動(dòng)力工程，25(2):152-155,2019.本所的蔡少輝研究員和中科院等離子體所FDS研究組也開(kāi)展過(guò)這方面的研究。SNMC應(yīng)用2022/10/1533二.粒子輸運(yùn)計(jì)算(19)該方法在“堆外探測(cè)室響應(yīng)”計(jì)算中得二.粒子輸運(yùn)計(jì)算(20)

粒子標(biāo)識(shí)(SN/MC)利用SN方法計(jì)算全系統(tǒng)網(wǎng)格通量精度較高的優(yōu)勢(shì)，把SN計(jì)算的網(wǎng)格通量作為MC輸運(yùn)計(jì)算的源，算出感興趣的計(jì)數(shù)量及各種標(biāo)識(shí)量，如粒子出生地、產(chǎn)地、反應(yīng)類(lèi)型等，補(bǔ)充MCNP程序沒(méi)有的計(jì)數(shù)功能[13-14]。

[13]上官丹驊，許海燕，鄧力，基于MonteCarlo方法的中子標(biāo)識(shí)輸運(yùn)計(jì)算，清華大學(xué)學(xué)報(bào)，47(S1)，p1062-1064，2019.[14]鄧力，謝仲生，碳氧比能譜測(cè)井的蒙特卡羅模擬，地球物理學(xué)報(bào)，44，p252-264，2019.*********粒子身份卡2022/10/1534二.粒子輸運(yùn)計(jì)算(20)粒子標(biāo)識(shí)(SN/MC)[13]四.MC方法的應(yīng)用(5)

2022/10/1535四.MC方法的應(yīng)用(5)2022/10/1135二.粒子輸運(yùn)計(jì)算(21)

2.4(t)特征值計(jì)算

含時(shí)無(wú)源齊次中子輸運(yùn)方程

(2.17)假設(shè)可以關(guān)于時(shí)間變量t分離變量

(2.18)代(2.18)式入(2.17)式得(2.19)等式左端為t的函數(shù)，右端為(r,E,)的函數(shù)，必為一常數(shù)，于是有2022/10/1536二.粒子輸運(yùn)計(jì)算(21)2.4(t)特征值計(jì)算含時(shí)無(wú)二.粒子輸運(yùn)計(jì)算(22)

(2.20)和

(2.21)解方程(2.19)、(2.20)得一系列特征值{i}和特征函數(shù){i}，則(2.22)其中1為最大本征值，1>0為超臨界；1<0為次臨界；1=0為臨界。為求本征值1，方程(2.21)可以寫(xiě)為(2.23)2022/10/1537二.粒子輸運(yùn)計(jì)算(22)(2.20)和(2.21)解方二.粒子輸運(yùn)計(jì)算(23)

當(dāng)>0時(shí)，被視為時(shí)間吸收，MC模擬不存在任何困難。然而當(dāng)t+/v<0，時(shí)，MC估計(jì)方法無(wú)效。此時(shí)引入(0<<1)，方程改寫(xiě)為(2.24)顯然方程(2.17)與方程(2.16)等價(jià)，的選取要保證。與原方程比較，相當(dāng)于方程右端散射部分增加了一項(xiàng)-散射。當(dāng)相對(duì)散射項(xiàng)較小時(shí)，方程解是存在唯一的。2022/10/1538二.粒子輸運(yùn)計(jì)算(23)當(dāng)>0時(shí)，被視為時(shí)間吸收，MC模二.粒子輸運(yùn)計(jì)算(24)

(t)本征值計(jì)算仍然存在較大不確定性，特別對(duì)含氫材料的次臨界問(wèn)題實(shí)驗(yàn)與數(shù)值模擬結(jié)果存在分叉現(xiàn)象。由于MCNP程序計(jì)算本征值是通過(guò)迭代keff本征值得到的((k-1)/)，而兩者之間只有臨界附近存在轉(zhuǎn)換關(guān)系，對(duì)深次臨界或超高臨界，顯然這種迭代是不成立的。[13]BrockwayD,SoranP&WhalenP,MonteCarloEigenvalueCalculation,LectureNotesinPhysics,Monte-CarloMethodsandApplicationsinNeutronics,Photonics,andStatisticalPhysics,Proceedings,Cadarache,Castle,France,1985,Springer-Verlag,Berlin,ISBN3-540-16070,1985.[14]邱有恒，鄧力，次臨界本征值計(jì)算中的時(shí)間選取方法，核科學(xué)與工程，27(2)，p182-186，2019.[15]邱有恒，應(yīng)陽(yáng)均，鄧力，次臨界系統(tǒng)中本征值與k本征值的關(guān)系,清華大學(xué)學(xué)報(bào)，47(S1)，p967-971,2019.問(wèn)題：對(duì)深次臨界問(wèn)題，除最大本征值外，是否應(yīng)該考慮次本征值？2022/10/1539二.粒子輸運(yùn)計(jì)算(24)(t)本征值計(jì)算仍然存在較大不確二.粒子輸運(yùn)計(jì)算(25)

本征值迭代算法：

[t0,,t1,t2,,ts](2.25)對(duì)深次臨界問(wèn)題仍然存在不確定性！改進(jìn)(本所)實(shí)踐證明這種處理基本有效，但(t)存在一定的波動(dòng)，時(shí)間步的選取很重要。2022/10/1540二.粒子輸運(yùn)計(jì)算(25)本征值迭代算法：[t0,,t二.粒子輸運(yùn)計(jì)算(26)

2.5占優(yōu)比(DominanceRatio)[16]針對(duì)穩(wěn)態(tài)本征值計(jì)算，占優(yōu)比決定了迭代收斂的速度，是目前臨界計(jì)算重點(diǎn)關(guān)心的問(wèn)題。設(shè)k0為確切的本征值，k1為最大本征值，則占優(yōu)比定義為[16]BrainNease,ForrestBrown,TaroUeki,DominanceratiocalculationswithMCNP,InternationalConferenceonthePhysicsofReactors,Interlaken,Switzerland,September14-19,2019.(2.26)越小，keff本征值迭代收斂越快。的方差為(2.27)2022/10/1541二.粒子輸運(yùn)計(jì)算(26)2.5占優(yōu)比(Dominan二.粒子輸運(yùn)計(jì)算(27)

keff滿足本征方程(2.28)其中L表示漏失，T表示碰撞，S表示散射，F(xiàn)表示裂變。

(2.29)方程(2.28)兩端同時(shí)減去F/k0，構(gòu)造迭代格式定義參數(shù)ke=k0+,>0。

Wielandt加速迭代格式定義為

(2.30)中子游動(dòng)期間，在裂變區(qū)的發(fā)生的每次碰撞，將產(chǎn)生2022/10/1542二.粒子輸運(yùn)計(jì)算(27)keff滿足本征方程(2.28)其二.粒子輸運(yùn)計(jì)算(28)

(2.31)個(gè)中子，這里wgt為中子權(quán)重。存入(2.32)個(gè)中子作為下一迭代的源。

Wielandt迭代得到的Wielandt與源迭代得到的之間滿足

(2.33)因?yàn)閣ielandt<，因此Wielandt迭代較源迭代收斂快[17]。2022/10/1543二.粒子輸運(yùn)計(jì)算(28)(2.31)個(gè)中子，這里wgt為中二.粒子輸運(yùn)計(jì)算(29)

[17]BrockwayD,SoranP&WhalenP,MonteCarloEigenvalueCalculation,LectureNotesinPhysics,Monte-CarloMethodsandApplicationsinNeutronics,Photonics,andStatisticalPhysics,Proceedings,Cadarache,Castle,France,1985,Springer-Verlag,Berlin,ISBN3-540-16070,1985.問(wèn)題：的選?。?022/10/1544二.粒子輸運(yùn)計(jì)算(29)[17]BrockwayD,S二.粒子輸運(yùn)計(jì)算(30)

2.6確定型MC方法

(a)直接抽樣：假定有n個(gè)離散隨機(jī)事件E1、E2、、En；每個(gè)事件發(fā)生的概率分別為p1,p2,…,pn。抽隨機(jī)數(shù)，求出滿足不等式(2.34)的j，則事件Ej以權(quán)重w發(fā)生。

(b)偏倚抽樣:

對(duì)上面提到的n個(gè)事件，構(gòu)造一組偏倚概率，。對(duì)同樣的隨機(jī)數(shù)，求出滿足不等式(2.35)2022/10/1545二.粒子輸運(yùn)計(jì)算(30)2.6確定型MC方法(a)直接的k，則事件Ek以權(quán)重發(fā)生。

顯然Ek可以不同于Ej。偏倚抽樣就是加大感興趣事件發(fā)生的概率，減少對(duì)計(jì)數(shù)貢獻(xiàn)小事件的概率，通過(guò)無(wú)偏修正保證計(jì)算結(jié)果無(wú)偏。但偏倚概率的選取需要有一定的理論指導(dǎo)。(c)確定型MC方法：對(duì)相同的n個(gè)事件，讓每個(gè)事件以各自的概率和權(quán)重發(fā)生:這種處理不會(huì)漏掉任何事件，適合共振區(qū)的處理。但計(jì)算時(shí)間和存儲(chǔ)量要相應(yīng)增加，適合并行處理。二.粒子輸運(yùn)計(jì)算(31)

2022/10/1546的k，則事件Ek以權(quán)重發(fā)生。顯然Ek可以不四.MC方法的應(yīng)用(3)特征的計(jì)算2022/10/1547四.MC方法的應(yīng)用(3)特征的計(jì)算2022/10/11472.7輸運(yùn)問(wèn)題反問(wèn)題（正算＋反演）中子輸運(yùn)的反問(wèn)題，已知照射量I(r),反演介質(zhì)密度(r)問(wèn)題：變量多，解不唯一，代數(shù)方程組奇異，……(2.36)?源靶成像二.粒子輸運(yùn)計(jì)算(32)

2022/10/15482.7輸運(yùn)問(wèn)題反問(wèn)題（正算＋反演）中子輸運(yùn)的反問(wèn)題，已知照2.8輻射輸運(yùn)

在ICF中涉及流體力學(xué)與光子輸運(yùn)耦合計(jì)算。由于系統(tǒng)隨時(shí)間變化，對(duì)時(shí)間變量離散后，在每一個(gè)時(shí)間步內(nèi)需要解流體力學(xué)與光子輸運(yùn)耦合方程組。由于溫度效應(yīng)的影響，需要解非線性輻射輸運(yùn)方程，采用迭代求解，計(jì)算量具大。早期求解多采用確定論方法，隨著并行計(jì)算機(jī)的發(fā)展，用隱式MC方法解輻射輸運(yùn)問(wèn)題成為可能[18]。[18]ToddUrbatsch,StrategiesforParallelImplicitMonteCarlo(U),LA-NM-87545,2019.二.粒子輸運(yùn)計(jì)算(33)

2022/10/15492.8輻射輸運(yùn)在ICF中涉及流體力學(xué)與光子輸運(yùn)耦合2.9探測(cè)器響應(yīng)計(jì)算

二.粒子輸運(yùn)計(jì)算(34)

圖13”3”NaI閃爍探測(cè)器響應(yīng)函數(shù)2022/10/15502.9探測(cè)器響應(yīng)計(jì)算二.粒子輸運(yùn)計(jì)算(34)圖13”2.10并行計(jì)算

MC方法具有數(shù)據(jù)獨(dú)立、通訊少、計(jì)算粒度大等特點(diǎn)，很適合并行計(jì)算，利用高性能并行計(jì)算機(jī)，可以縮短計(jì)算的周期；設(shè)計(jì)合理的隨機(jī)數(shù)發(fā)生器，可以確保并行和串行計(jì)算結(jié)果一致；長(zhǎng)周期隨機(jī)數(shù)發(fā)生器設(shè)計(jì)，以滿足現(xiàn)代問(wèn)題的需要；

MCNP-4C、MCNP-5MPI并行計(jì)算[19-20]；

GPU并行[21]。二.粒子輸運(yùn)計(jì)算(35)

[19]DengLi,etal.,ParallelizationofMCNPMonteCarloNeutronandPhotonTransportCodeinPVMandMPI,JNucl.Sci.Tech.,36(7):626-629,2000.[20]鄧力，張文勇，徐涵，蒙特卡羅程序MCNP-II與MCNP-5并行效率比較，計(jì)算機(jī)工程與科學(xué)，31(A1):p185-1875，2009.[21]ChunyeGong,DengLi,etal.,AcceleratingMCNP-bastedMonteCarloSimulationsforNeutronTransportonGPU,JournalofRadiationOncologyBiologyPhysics,Vol.81,Issue2,p8157-815,2019.

2022/10/15512.10并行計(jì)算MC方法具有數(shù)據(jù)獨(dú)立、通訊少、計(jì)算粒度大三.發(fā)展與挑戰(zhàn)(1)

隨著計(jì)算機(jī)的快速發(fā)展，MC方法應(yīng)用于反應(yīng)堆物理分析獲得巨大增長(zhǎng)，過(guò)去MC方法主要用于臨界系統(tǒng)的keff計(jì)算，最近在反應(yīng)堆功率分布，各種參數(shù)的優(yōu)化設(shè)計(jì)(如eff、、eff、反應(yīng)性系數(shù)、多普勒效應(yīng)、占優(yōu)比等)；三維功率分布收斂性的確定；消除置信區(qū)間偏倚；裂變概率的多次迭代；多普勒溫度效應(yīng)反饋；燃耗和裂變產(chǎn)物累計(jì)的模擬等。MC方法變得通用靈活，但給研制者帶來(lái)新的挑戰(zhàn)[22]。[22]F.B.Brown,W.R.Martin,R.D.Mosteller,MonteCarlo–AdvanceandChallenges,PHYSOR’08,MontecarloWorkshop,Interlaken,Switzerland,Sept14-19,2019.2022/10/1552三.發(fā)展與挑戰(zhàn)(1)隨著計(jì)算機(jī)的快速發(fā)展，MC[23]ToddUrbatsch,StrategiesforParallelImplicitMonteCarlo(U),LA-NM-87545,2019.[24]2019年美國(guó)LosAlamos實(shí)驗(yàn)室在10萬(wàn)億次超級(jí)并行機(jī)，采用隱式MC方法和區(qū)域分解，使用2000個(gè)處理器，進(jìn)行三維全物理動(dòng)態(tài)數(shù)值實(shí)驗(yàn)，用1個(gè)月時(shí)間完成了模擬計(jì)算.[25]HowdoesASCIactuallycompletemulti-month1000-processormilestonesimulations?Pitfall:failedjobstarts;badCPUS;badinterconnect;badI/O;baddumpfiles.

Restart!Restart!Restart!

高性能計(jì)算與武器物理，第3卷第8期，2019CTBT后，數(shù)值模擬需要由實(shí)驗(yàn)上升到理論，對(duì)計(jì)算方法提出了新的挑戰(zhàn)！三.發(fā)展與挑戰(zhàn)(2)

2022/10/1553[23]ToddUrbatsch,Strategies三.發(fā)展與挑戰(zhàn)(3)

由于上世紀(jì)六、七十年代建成的許多反應(yīng)堆臨近壽期，需要延壽或重建，近年MC方法在建模、損耗及參數(shù)優(yōu)化設(shè)計(jì)，特別是針對(duì)新堆設(shè)計(jì)上發(fā)揮重要作用。

燃耗計(jì)算中的裂變產(chǎn)物；挑戰(zhàn)程序規(guī)模極限；多物理、多尺度、時(shí)空動(dòng)力學(xué)模擬；PinbyPin問(wèn)題；核數(shù)據(jù)不確定性；CTBT帶來(lái)的挑戰(zhàn)；2022/10/1554三.發(fā)展與挑戰(zhàn)(3)由于上世紀(jì)六、七十年代建成的許多三.發(fā)展與挑戰(zhàn)(5)挑戰(zhàn)計(jì)算能力的問(wèn)題

2022/10/1555三.發(fā)展與挑戰(zhàn)(5)挑戰(zhàn)計(jì)算能力的問(wèn)題2022/10/115三.發(fā)展與挑戰(zhàn)(6)

復(fù)雜物理量的分量計(jì)算2022/10/1556三.發(fā)展與挑戰(zhàn)(6)復(fù)雜物理量的分量計(jì)算2022/10/1三.發(fā)展與挑戰(zhàn)(7)

取消二次中子源項(xiàng)目的計(jì)算

堆芯徑向探測(cè)器堆芯軸向組建結(jié)構(gòu)2022/10/1557三.發(fā)展與挑戰(zhàn)(7)堆芯徑向探測(cè)器堆芯軸向組建結(jié)構(gòu)2022三.發(fā)展與挑戰(zhàn)(8)

其它

穩(wěn)定性診斷連續(xù)材料及計(jì)數(shù)自動(dòng)降方差技巧MC損耗源微擾顆粒問(wèn)題的模擬連續(xù)能量的微擾計(jì)算

………2022/10/1558三.發(fā)展與挑戰(zhàn)(8)其它穩(wěn)定性診斷2022/10/1四.MC方法的應(yīng)用(1)核探測(cè)

煤質(zhì)成分分析；毒品和隱藏爆炸物探測(cè)；集裝箱探測(cè)；利用中子(或光子)源照射被探測(cè)物，通過(guò)探測(cè)器響應(yīng)判斷被探測(cè)物的元素成份[26－28]。

[26]R.P.Gardner,etal.,ANewDirectionBiasingApproachforMonteCarloSimulation,Nucl.Sci.Eng.Vol.98,51-63,1988.[27]R.P.Gardner,etal.,AMonteCarloDirectionBiasingApproachintheLaboratorySystemforIsotropicNeutronCenter-of-MassScatteringIncludingHydrogen,Nucl.Sci.Eng.Vol.108,240-246,1991.[28]DengLi,etal.,AMonteCarloModelforGamma-RayKlein-NishinaScatteringProbabilitiestoFiniteDetectors,JNucl.Sci.Tech.,33(9):736-740,2019.2022/10/1559四.MC方法的應(yīng)用(1)核探測(cè)煤質(zhì)成分分析；[26]R.

特征計(jì)算直穿計(jì)數(shù)分線光和群光(即連續(xù)光)兩部分：

期望估計(jì)產(chǎn)光指向概率估計(jì)四.MC方法的應(yīng)用(2)2022/10/1560特征計(jì)算直穿計(jì)數(shù)分線光和群光(即連續(xù)光)兩部分：期Ex1.

隱藏爆炸物探測(cè)四.MC方法的應(yīng)用(3)2022/10/1561Ex1.隱藏爆炸物探測(cè)四.MC方法的應(yīng)用(3)2022/1四.MC方法的應(yīng)用(3)特征的計(jì)算2022/10/1562四.MC方法的應(yīng)用(3)特征的計(jì)算2022/10/1162Ex2.

C/O碳氧比能譜測(cè)井圖

C/O能譜測(cè)井儀器及地層結(jié)構(gòu)圖

四.MC方法的應(yīng)用(4)2022/10/1563Ex2.C/O碳氧比能譜測(cè)井圖C/O能譜測(cè)井儀器及地層四.MC方法的應(yīng)用(5)

2022/10/1564四.MC方法的應(yīng)用(5)2022/10/1164醫(yī)學(xué)劑量計(jì)算

核技術(shù)在醫(yī)學(xué)上已有廣泛應(yīng)用。如硼中子俘獲治療(BNCT)，需要從醫(yī)學(xué)CT、MRI提取人體組織的密度和邊界，在臨床規(guī)定的時(shí)間和精度要求下，完成腫瘤劑量的計(jì)算，并確定照射部位和照射時(shí)間，實(shí)施治療。由于計(jì)數(shù)目標(biāo)量多，過(guò)程復(fù)雜，傳統(tǒng)MC程序(如MCNP)很難滿足臨床對(duì)計(jì)算精度和計(jì)算時(shí)間的要求。構(gòu)建特殊網(wǎng)格模型，發(fā)展特殊網(wǎng)格下的快速粒子徑跡算法和計(jì)數(shù)方法是MC研究的另一新領(lǐng)域[29-30]。[29]M.W.Frandsen,RapidGeometryInterrogationforUniformVolumeElement-BasedBNCTMonteCarloParticleTransportSimulation,ProceedingsoftheEighthInternationalSymposiumonNeutronCaptureTherapy,PlenumPress,NewYork,2019.[30]DengLi,LiGang,etal.,MCDBMonteCarloCodewithFastTrackTechniqueandMeshTallyMatrixforBNCT.

JNucl.Sci.Tech.44,1518–1525,2019.四.MC方法的應(yīng)用(6)

2022/10/1565醫(yī)學(xué)劑量計(jì)算核技術(shù)在醫(yī)學(xué)上已有廣泛應(yīng)用。如硼中子俘獲治療(BPA+BSH四.MC方法的應(yīng)用(7)

2022/10/1566BPA+BSH四.MC方法的應(yīng)用(7)2022/10/11CTMRI臨床治療MonteCarloSimulationVoxelmodelIrradiationlocationHistogramofdose3dreconstruction快速?gòu)桔E算法MCDB劑量計(jì)算四.MC方法的應(yīng)用(8)

2022/10/1567CTMRI臨床治療MonteCarloSimulatioDoseDistributionDICOMdataMCinputCT、MRIimagesdefcardMCrunCross-sectionoutputKERMAfactor

MCDB_TPSPro-processorPost-processorNeutronbeamAccelerator2022/10/1568DoseDistributionDICOMdataMC快速?gòu)桔E算法neighborrelationsamongthemeshes

(x0,y0)(i,j)i+1j+1i+2j+3dxdyx=|1/u|y=|1/v|i+3j+2(u,v)CollisionPoint??dj=dls2-dls12022/10/1569快速?gòu)桔E算法neighborrelationsamongVoxelmodelAnalyticmodel2022/10/1570VoxelmodelAnalyticmodel2022/1

43piecesCT，thedistanceis3mmbetweentwopiecesofCT，512×512pixel/eachCT.

CT8voxelmodel:8×8pixels/meshmeshes=64×64×43=176128voxelsize=0.3703×0.3703×0.3(mm3)

CT4voxelmodel:4×4pixels/meshmeshes=128×128×43=704512voxelsize=0.1852×0.1852×0.3(mm3)Ex.Patientmodel四.MC方法的應(yīng)用(8)

2022/10/157143piecesCT，thedistancei四.MC方法的應(yīng)用(9)

2022/10/1572四.MC方法的應(yīng)用(9)2022/10/1172

CT4voxelmodel3-Dreconstruction四.MC方法的應(yīng)用(10)

2022/10/1573CT4voxelmodel3-DreconstrucMCDBMCNPT.N.MCNPMCDBMCNPMCDBF.N.P.四.MC方法的應(yīng)用(11)

2022/10/1574MCDBT.N.MCNPMCDBMCNPMCDBF.N.P.Speedup=3.47四.MC方法的應(yīng)用(12)

2022/10/1575Speedup=3.47四.MC方法的應(yīng)用(12)23DDoseSurfacePhotondosesurfaceNeutrondosesurface2022/10/15763DDoseSurfacePhotondosesur五.國(guó)內(nèi)現(xiàn)狀(1)

國(guó)內(nèi)使用MC程序的單位及個(gè)人很多，但多數(shù)程序是引進(jìn)開(kāi)發(fā)的。國(guó)內(nèi)在壓水堆、快堆、聚變堆、高溫堆方面的發(fā)展帶動(dòng)方法和軟件的發(fā)展；在堆物理方面，西安交大、清華大學(xué)、中科大、上海交大、哈工程、華電等建立了核科學(xué)工程學(xué)院，培養(yǎng)了大量人才；核動(dòng)力院、原子能院、中科院等離子所、728、585、西北核技術(shù)研究所、中物院等都在發(fā)展MC方法和程序，解決工程中的問(wèn)題；北師大、吉大等在統(tǒng)計(jì)物理方面應(yīng)用MC解決問(wèn)題的優(yōu)勢(shì)明顯。

我所在方法研究及應(yīng)用程序研制上有近50年的歷史。MC方法和軟件在產(chǎn)品設(shè)計(jì)、中子物理學(xué)、照相、測(cè)試、輻射輸運(yùn)、ICF、聚變堆包層設(shè)計(jì)、ITER等領(lǐng)域都有應(yīng)用，有多個(gè)自研和開(kāi)發(fā)的MC程序。2022/10/1577五.國(guó)內(nèi)現(xiàn)狀(1)國(guó)內(nèi)使用MC程序的單位及1987:MORSE-CGT程序(在ORNL-RSIC入庫(kù))*2019:MCCO高分辨率-探測(cè)MC程序

(JNucl.Sci.Tech.,

33(9),2019)2019:實(shí)現(xiàn)MCNP程序的PVM和MPI并行化**

(JNucl.Sci.Tech.,

36(7),2000)2019:MCMG多群P3中子輸運(yùn)MC程序(JNucl.Sci.Tech.

37(7),2000)2000:MCMG多群P3中子輸運(yùn)-燃耗耦合計(jì)算MC程序

(AnnNucl.Energy,27(10),201