蒙特卡罗方法及应用实验讲义2016资料

蒙特卡罗方法及应用

实验讲义

东华理工大学核工系

2016.8

实验一 蒙特卡罗方法基本思想

一、实验目的

1、了解蒙特卡罗方法方法的基本思想；

2、掌握蒙特卡罗方法计算面积、体积的方法；

3、掌握由已知分布的随机抽样方法。

二、实验原理

Monte Carlo 方法，又称统计模拟方法或计算机随机模拟方法，是一种基于“随机数”进行数值模拟的方法，一种采用统计抽样理论近似求解物理或数学问题的方法。

如待求量可以表述成某些特征量的期望值、某些事件出现的概率或两者的函数形式，那么可采用蒙特卡罗方法求解。在求解某些特征量的期望值或某些事件出现的概率时，必须构建合符实际的数学模型。例如采用蒙特卡罗方法计算某函数所围面积时，构建的数学模型是构造一已知面积的可均匀抽样区域，在该区域投点，由伯努利定理大数定理可知，进入待求区域投点的频率依概率1收敛于该事件出现的概率（面积之比）。

由已知分布的随机抽样方法指的是由已知分布的总体中抽取简单子样。具体方法很多，详见教材第三章。 三、实验内容

1、安装所需计算工具（MATLAB 、fortran 、C++等）；

2、学习使用rand(m,n)、unifrnd(a,b,m,n)函数

3、求解下列问题：

3.0、蒲丰氏投针求圆周率。

3.1、给定曲线y =2 – x 2 和曲线y 3 = x 2，曲线的交点为：P 1( – 1，1 )、P 2( 1，1 )。曲线围成平面有限区域，用蒙特卡罗方法计算区域面积；

3.2

、计算1z z ?≥??≤??所围体积

其中{(,,)|11,11,02}x y z x y z Ω=-≤≤-≤≤≤≤。

4、对以下已知分布进行随机抽样：

4.1、()()[]2

3

321,0,12

f x x x x =+

-∈； 4.2、()()

()[]11

,1,21E f x f x x E k E =

?∈+ 其中()()()()()2

123

221111211411ln 212221E x f x E x x x x E k E E E E E ?+-??=+-+? ?????

?+???=-?+++-

???+???

。

四、实验报告编写

1、给出各题的抽样程序并解释语句的含义；

2、给出3.1和3.2抽样结果误差随抽样次数的关系图，并解释原因；

表1 实验记录表

3、给出4.1和4.2的抽样框图、试验累积频率与理论累积频率关系图，并给出抽样次数（>106）与抽样时间。

实验二 由已知分布的随机抽样方法

一、实验目的

1、掌握由已知分布的随机抽样方法。

2、用编程语言实现某具体随机抽样方法。 二、实验原理

由已知分布的随机抽样方法指的是由已知分布的总体中抽取简单子样。具体方法很多，本实验综合直接抽样方法、挑选抽样方法和替换抽样方法，以散射方位角余弦分布的抽样为例。实验原理详见教材对应章节。 1.连续型分布的直接抽样方法

对于连续型分布，如果分布函数F(x) 的反函数F －1(x)存在，则直接抽样方法是：

1()F X F ξ-=

2.挑选抽样方法

为了实现从己知分布密度函数f(x)抽样，选取与f(x)取值范围相同的分布密度函数h(x)，如果

()

sup

()

x f x M h x -∞<<∞

=<∞ 则挑选抽样方法为：

3.替换法抽样方法

为了实现某个复杂的随机变量 y 的抽样，将其表示成若干个简单的随机变量x 1，x 2，…，x n 的函数

12(,,

,)n y g x x x =

得到 x 1，x 2，…，x n 的抽样后，即可确定 y 的抽样，这种方法叫作替换法

抽样。

蒲丰氏问题的算法 如何产生任意的（x,θ）？

x 在［0，a ］上任意取值，表示x 在［0，a ］上是均匀分布的，其分布密度函数为：

11/,0()0,

a x a

f x ≤≤?=?

?其他 类似地，θ的分布密度函数为：

21/,0()0,f πθπ

θ≤≤?=?

?

其他 因此，产生任意的（x,θ）的过程就变成了由f 1(x)抽样x 及由f 2(θ)抽样θ的过程了。由此得到：

1

2

x a ξθπξ==

其中ξ1，ξ2均为（0,1）上均匀分布的随机变量。

每次投针试验，实际上变成在计算机上从两个均匀分布的随机变量中抽样得到（x,θ），然后定义描述针与平行线相交状况的随机变量s(x,θ)，为

1,sin (,)0x l s x θ

θ≤??=?

?

,其他 如果投针Ｎ次，则

1

1(,)N

N i i i s s x N θ==∑

是针与平行线相交概率Ｐ的估计值。事实上，

12sin 0

0(,)()()d d d d 2ππl P s x f x f x x l

a a

πθθθθθ===

???

?

于是有

22πN

l l

aP as =

≈

1、给出源程序程序并解释语句的含义；

2、作出抽样框图、试验累积频率与理论累积频率关系图，并给出抽样次数（>106）与抽样时间。

实验三MCNP方法在实验核物理中的应用

一、实验目的

1、了解MCNP程序运行流程；

2、掌握MCNP输入文件编写规范；

3、理解模拟内容、并能编写输入文件、运行，并获得计算结果；

二、实验原理

MCNP是一种常见的粒子输运模拟软件，软件的安装、运行和输入文件编写方法详见相关参考资料。

MCNP输入文件编写完成后，先确认输入模型是否正确，在DOS环境下进行，打开运行DOS环境，进行以下操作：

DOS命令操作命令含义

Mcnp i=name.inp o=name.o 打开画图框

PX vx 输出模型在x=vx面上的切面

PY vy 输出模型在y=vy面上的切面

PZ vz 输出模型在z=vz面上的切面FACTOR m 将输出图放大1/m倍

Extent a b 切面沿两坐标轴方向分别放大

ORIGIN X Y Z 定义画图中心位置（X,Y,Z）

三、实验内容

1、学习MCNP程序常见各种运行方法；

2、编写以下问题的输入文件；

2.1对课堂讲解的实例，模拟溴化镧探测器对点源的能谱，实验做一遍。

2.2有一HPGe探测器，结构如图1所示。分别给出位于探测器轴心、距离探测器晶体中心25cm处的137Cs源、60Co、131I源对应特征γ射线的探测效率（计算时相应特征射线的源粒子至少为107个），并给出三者混合源（活度比为1:1:3）的能谱图（源发射总粒子数大于3×108个）。

Cu 的尺寸

1.75cm

0.09cm

Al 壳0.05cm

图1 HPGe 探测器结构图

四、实验报告

1、给出2.1和2.2的MCNP 输入文件并解释每一行的含义；

2、分别运行实例，给出实验结果，并对结果进行分析。

实验四 MCNP 模拟计算γ射线造成的剂量

一、实验目的

1、掌握应用软件MCNP 、应用范围以及在辐射剂量计算和防护中的作用；

2、进一步掌握MCNP 程序基本用法；

3、利用MCNP 解决一个简单的求解γ射线在空气、组织等效材料(肌肉)中造成的剂量沉积的计算问题，并进行结果分析，得出结论；

4、利用MCNP 程序解决实际工作中碰到的实际问题； 二、实验内容

1、学习MCNP 程序的基本组成、操作方法以及问题描述文件的写法；

2、利用MCNP 程序计算简单的γ射线源在空气、肌肉模型中的剂量沉积分布，并对计算结果进行分析并绘图，得出结论，调整数据重新计算，并与理论计算结果进行比较； 三、内容简介

1、MCNP 程序的计算流程如下图1所示：

2、MCNP 输入文件

通过这个文件描述并建立一个蒙特卡罗计算问题，

对问题的几何结构、材料、记数要求等给以描述，如果需要，便可直接运行。该文件的格式如下： 栅元卡1 栅元卡2 。。。

栅元卡n

空行分隔符

曲面卡1

曲面卡2

。。。

曲面卡n

空行分隔符

数据卡1

数据卡2

。。。

数据卡n

空行分隔符(optional)

其它选择项(optional)

其中栅元卡用来描述由不同的封闭曲面分割的立体空间区域，并用独有的数字ID号加以标示，同时在各个栅元卡中说明包围该区域的曲面类型(曲面卡)、填充该区域的材料类型(材料卡)以及对应的材料密度等；

曲面卡是用来描述不同类型曲面的，并用独有的数字ID号加以标示，最终曲面卡被应用在栅元卡中，并利用交(与)、联(或)、补(非)这些逻辑运算符号联合不同曲面组成所需要的复杂的栅元。在mcnp中支持的常见曲面类型见参考文献[3,4]。

数据卡类型很多，主要有粒子类型标识卡mode、重要性卡imp、通用源卡sdef、粒子计数器卡Fn、材料描述卡Mn以及粒子截断卡(nps或ctme)等，数据卡的类型涉及到了方方面面，类型很多，具体请见参考文献[3,4]。

下面利用一个简单的例子来配合说明mcnp中的输入卡(inp)的编写格式。

3、一个简单的说明例子

为说明如何填写INP文件，这里例举一个简

单问题。如图3所示，在一个边长10cm的石墨立

方体3中有两个半径0.5cm的球形空间，球1中充

满氧气，球2是铁球。在球1中置一14MeV各向

同性中子点源，计算球2外表面与能量相关的中

子通量。建立的INP文件如下：

SAMPLE PROBLEM INPUT DECK

1 1 –0.0014 -7

2 2 –7.86 -8

3 3 –1.60 1 –2 –3

4 –

5

6

7 8

4 0 -1:2:3:-4:5:-6

空行

1 PZ -5

4

3

2

1

8

7

43

2

1

y 图3 例子的几何示意图

2 PZ 5

3 PY 5

4 PY -5

5 PX 5

6 PX -5

7 S 0 –4 –2.5 .5

8 S 0 4 4 .5

←空行

MODE p

IMP:p 1 1 1 0

SDEF POS=0 –4 –2.5 ERG=14

F2:n 8

E0 1E-5 1E-4 1E-3 .01 .1 1 2 3 4 5 6 7 8 9 10 11 12 13

M1 8016 1

M2 26000 1

M3 6000 1

NPS 100000

←空行

本例中没有信息块，第一行是标题卡，之后至空格前为栅元块。栅元卡上依次填写栅元号、材料号、密度和构成栅元界面的曲面号(带正负号)，这里定义了4个栅元：栅元1由球面7围成，里面填充材料1(16O2气体)，密度是0.0014g/cm3；栅元2由球面8围成，填充材料2(铁)，密度7.86g/cm3；栅元3由平面1、2、3、4、5、6围成，不包括球面7、8以内的空间，填充材料3(石墨)，密度1.6g/cm3；栅元4是栅元3以外的空间，为真空。

曲面卡上需要填写曲面号、曲面类型和曲面参数，本例中定义了8个曲面，前6个为与原点距离5cm垂直于各坐标轴的平面，后两个是半径0.5cm的球面，球心分别在(0,-4,-2.5)和(0,4,4)。

数据块中指定了问题类型、源、记数方式、材料和运行粒子数，各卡数据项的意义如下：

MODE卡问题类型是中子输运

IMP卡4个栅元的中子重要性分别是1 1 1 0

SDEF卡位于(0,-4,-2.5)、能量14MeV的各向同性点源

F2卡在曲面8上做中子通量记数

E0卡对记数能量分区，1~14MeV之间间隔为1MeV，1MeV~10-5MeV之间间隔为一个数量级

M1卡材料1是16O核素

M2卡材料2是Fe元素

M3卡材料3是C元素

NPS卡运行源粒子数100000

以上例子仅用于说明INP文件格式，有关各输入卡的详细内容，具体使用方法见参考文献[3,4]。

四、实验步骤

1、编写对应任务的输入文件，如实例Dose.inp。

2、运行实例Dose，得到计算结果文件Dose.o，对Dose.o文件进行分析，了解并熟悉计算过程，得到并分析计算结果；

3、理论计算空气比释动能值，并与上述计算结果计算得到的空气比释动能值进行比较，分析差异，给出分析结果。

4、在上述例子的基础上，对计算空间进行栅元细分，得到更为精细的计算结果，并按照前面顺序对计算结果进行剂量沉积分析。

五、问题求解

1、问题

1个1Ci的Cs-137源，计算距离其30cm处的空气比释动能K、肌肉材料的吸收剂量D。

2、设计要求：

（1）所设计的源为各向同性通用源：即用SDEF卡设置源；

（2）计算光子在肌肉和空气中所沉积的能量；

（3）所描述的几何有一定数量的栅元，几何形状可以按照自己的兴趣选择；

3、光子作用的材料

肌肉材料如下：

1001. -0.102

6000. -0.143

7014. -0.034

8016. -0.710

11023. -0.001

15031. -0.002

16032. -0.003

17000. -0.001

19000. -0.004

空气材料如下

6000. -0.000124

7014. -0.7553

8016. -0.2318

18000. -0.0128989

4、输入卡的编写

5、运行结果

六、思考问题

1、什么是栅元？MCNP用什么方法描述栅元几何？

2、初次编写inp文件，难免会出现各种错误，请试列举你在计算过程中遇到的问题，并给出你的解决方案。

3、把实例中的栅元3、

4、5分别修改位置、半径大小，运行获得相关栅元图形和各个栅元的吸收剂量。

七、实验报告

1、给出本实验的MCNP输入文件并解释每一句的含义；

2、运行实例，给出实验结果，并对结果进行分析。

附录

表1：人体组织材料中光子的注量到剂量率的转换因子

(ANSI/ANS结果(左)，ICRP21号报告(右))

表2：空气材料中光子的注量到比释动能的转换因子