蒙特卡洛方法在中子输运中的应用

1 中子输运理论与数值方法中子输运理论与数值方法 课程作业课程作业 蒙特卡洛方法蒙特卡洛方法 2 目录目录 1 前言 3 2 蒙特卡洛方法概述 3 2 1 蒙特卡洛方法的基本思想 4 2 2 蒙特卡洛方法的收敛性 误差 4 2 2 1 蒙特卡洛方法的收敛性 4 2 2 2 蒙特卡洛方法的误差 5 2 3 蒙特卡洛方法的特点 6 2 4 蒙特卡洛方法的主要应用范围 7 3 随机数 7 3 1 线性乘同余方法 9 3 2 伪随机数序列的均匀性和独立性 9 3 2 1 伪随机数的均匀性 9 3 2 2 伪随机数的独立性 10 4 蒙特卡洛方法在粒子输运上的应用 10 4 1 屏蔽问题模型 10 4 2 直接模拟方法 11 4 2 1 状态参数与状态序列 11 4 2 2 模拟运动过程 12 4 2 3 记录结果 15 4 3 蒙特卡洛方法的效率 16 5 蒙特卡洛方法应用程序 MCNP 17 5 1 MCNP 简述 17 5 2 MCNP 误差的估计 18 5 3 MCNP 效率因素 19 6 结论 19 参考文献 20 3 1 1 前言前言 半个多世纪以来 由于科学技术的发展和电子计算机的发明 蒙特卡洛 Monte Carlo 方法作为一种独立的方法被提出来 并首先在核武器的试验与研 制中得到了应用 蒙特卡洛方法是一种计算方法 但与一般数值计算方法有很 大区别 它是以概率统计理论为基础的一种方法 由于蒙特卡洛方法能够比较 逼真地描述事物的特点及物理实验过程 解决一些数值方法难以解决的问题 因而该方法的应用领域日趋广泛 蒙特卡洛模拟计算是解决中子在介质中输运 较为成熟 有效的方法 对于原子能 辐射防护 剂量学和辐射生物物理学等 研究领域实际问题的计算 都可以利用蒙特卡洛方法予以实现 粒子输运过程可以用玻耳兹曼方程加以描述 然而 以此基础上发展起来 的近似数值方法如扩散近似法 离散坐标方法在处理截面与能量相关以及散射 各向异性介质 复杂几何条件问题时碰到了较大困难 而蒙特卡洛方法在处理 这类问题时得心应手 有很强的解题能力 并且近似较少 接近于真实情况 粒子辐射问题计算通常有输运方程法 蒙特卡洛法 MC 法 实验测量法以 及经验法等几种方法 蒙特卡洛计算法又称随机抽样法或统计试验法 是基于 计算机模拟的思想 抓住物理过程的数量和几何特征 进行数字模拟试验 该 方法是求解辐射输运问题的一种相当成熟和有效的方法 而且它对于各种复杂 问题 具有良好的通用性 实用性相当广泛 几乎涉及核科学的各个领域 本 文主要介绍蒙特卡洛的概念 原理和应用及研究现状 2 2 蒙特卡洛方法概述蒙特卡洛方法概述 蒙特卡洛方法又称随机抽样技巧或统计试验方法 半个多世纪以来 由于 科学技术的发展和电子计算机的发明 这种方法作为一种独立的方法被提出来 并首先在核武器的试验与研制中得到了应用 蒙特卡洛方法是一种计算方法 但与一般数值计算方法有很大区别 它是以概率统计理论为基础的一种方法 由于蒙特卡洛方法能够比较逼真地描述事物的特点及物理实验过程 解决一些 数值方法难以解决的问题 因而该方法的应用领域日趋广泛 蒙特卡洛方法的主要组成部分有 1 概率密度函数 pdf 必须给出描述一个物理系统的一组概率密度函数 4 2 随机数产生器 能够产生在区间 0 1 上均匀分布的随机数 3 抽样规则 如何从在区间 0 1 上均匀分布的随机数出发 随机抽取服从给定 的pdf的随机变量 4 模拟结果记录 记录一些感兴趣的量的模拟结果 5 误差估计 必须确定统计误差 或方差 随模拟次数以及其它一些量的变 化 6 减少方差的技术 利用该技术可减少模拟过程中计算的次数 7 并行和矢量化 可以在先进的并行计算机上运行的有效算法 2 12 1 蒙特卡洛方法的基本思想蒙特卡洛方法的基本思想 可以通俗地说 蒙特卡洛方法是用随机试验的方法计算积分 即将所要计 算的积分看作服从某种分布密度函数 f r 的随机变量 r 的数学期望 MERGEFORMAT 2 1 0 gg r f r dr 通过某种试验 得到 N 个观察值 r1 r2 rN 用概率语言来说 从分布密 度函数 f r 中抽取 N 个子样 r1 r2 rN 将相应的 N 个随机变量的值 g r1 g r2 g rN 的算术平均值 作为积分的估计值 近似值 1 1 N Ni i gg r N 为了得到具有一定精确度的近似解 所需试验的次数是很多的 通过人工 方法作大量的试验相当困难 甚至是不可能的 因此 蒙特卡洛方法的基本思 想虽然早已被人们提出 却很少被使用 本世纪四十年代以来 由于电子计算 机的出现 使得人们可以通过电子计算机来模拟随机试验过程 把巨大数目的 随机试验交由计算机完成 使得蒙特卡洛方法得以广泛地应用 在现代化的科 学技术中发挥应有的作用 2 22 2 蒙特卡洛方法的收敛性 误差蒙特卡洛方法的收敛性 误差 蒙特卡洛方法作为一种计算方法 其收敛性与误差是普遍关心的一个重要 问题 2 2 1 蒙特卡洛方法的收敛性 由前面介绍可知 蒙特卡洛方法是由随机变量 X 的简单子样 X1 X2 XN的算术平均值 作为所求解的近似值 由大数定 1 1 N Ni i XX N 5 律可知 如 X1 X2 XN独立同分布 且具有有限期望值 则 即随机变量 X 的简单子样的算术平均值 当子样 1 lim N N PXE X N X 数 N 充分大时 以概率 1 收敛于它的期望值 E X 2 2 2 蒙特卡洛方法的误差 蒙特卡洛方法的近似值与真值的误差问题 概率论的中心极限定理给出了 答案 该定理指出 如果随机变量序列 X1 X2 XN独立同分布 且具有 有限非零的方差 2 即 f X 是 X 的分布密度 22 0 xE Xf x dx 函数 则 2 21 2 lim x t N x N N PXE Xxedt MERGEFORMAT 2 2 当 N 充分大时 有如下的近似式 2 2 0 2 1 2 t N PXE Xedt N MERGEFORMAT 2 3 其中 称为置信度 1 称为置信水平 这表明 不等式 近似地以概率 1 成立 且误差收敛速度的阶为 N XE X N 1 2 O N 通常 蒙特卡洛方法的误差 定义为 MERGEFORMAT 2 4 N 上式中与置信度 是一一对应的 根据问题的要求确定出置信水平后 查标准正态分布表 就可以确定出 常用的 与的对应关系为 0 5 0 6745 0 05 0 96 0 003 3 蒙特卡洛方法的误差为概率误 差 这与其他数值计算方法是有区别的 误差中的均方差 是未知的 必须使 用其估计值 6 MERGEFORMAT 2 5 22 11 11 NN ii ii XX NN 来代替以求出均方差 由式 MERGEFORMAT 2 4 可知当给定置信度 后 误差 由 和 N 决定 要减小 或者是增大 N 或者是减小方差 2 在 固 定的情况下 要把精度提高一个数量级 试验次数 N 需增加两个数量级 因此 单纯增大 N 不是一个有效的办法 另一方面 如能减小估计的均方差 比如 降低一半 那误差就减小一半 这相当于 N 增大四倍的效果 因此降低方差的 各种技巧 引起了人们的普遍注意 2 32 3 蒙特卡洛方法的特点蒙特卡洛方法的特点 作为一种统计试验方法 蒙特卡洛方法因其优点在诸多领域内有着广泛 但同时存在一些缺点 蒙特卡洛的主要优点有 1 能够比较逼真地描述具有随机性质的事物的特点及物理实验过程 蒙特 卡洛方法可以部分代替物理实验 甚至可以得到物理实验难以得到的结果 用 蒙特卡洛方法解决实际问题 可以直接从实际问题本身出发 而不从方程或数 学表达式出发 它有直观 形象的特点 2 受几何条件限制小 在计算 s 维空间中的任一区域 Ds上的积分 时 无论区域 Ds的形状多么特殊 只要能给 1212 s ss D gg x xx dx dxdx 出描述 Ds的几何特征的条件 就可以从 Ds中均匀产生 N 个点 12 iii s xxx 得到积分的近似值 其中 Ds为区域 Ds的体积 这 12 1 N iii s Ns i D gg xxx N 是数值方法难以作到的 3 收敛速度与问题的维数无关 由误差定义可知 在给定置信水平情况下 蒙特卡洛方法的收敛速度为 与问题本身的维数无关 维数的变化 1 2 O N 只引起抽样时间及估计量计算时间的变化 不影响误差 也就是说 使用蒙特 卡洛方法时 抽取的子样总数 N 与维数 s 无关 维数的增加 除了增加相应的 计算量外 不影响问题的误差 这一特点 决定了蒙特卡洛方法对多维问题的 适应性 4 具有同时计算多个方案与多个未知量的能力 对于那些需要计算多个方 7 案的问题 使用蒙特卡洛方法有时不需要像常规方法那样逐个计算 而可以同 时计算所有的方案 其全部计算量几乎与计算一个方案的计算量相当 例如 对于屏蔽层为均匀介质的平板几何 要计算若干种厚度的穿透概率时 只需计 算最厚的一种情况 其他厚度的穿透概率在计算最厚一种情况时稍加处理便可 同时得到 5 误差容易确定 对于一般计算方法 要给出计算结果与真值的误差并不 是一件容易的事情 而蒙特卡洛方法则不然 根据蒙特卡洛方法的误差公式 可以在计算所求量的同时计算出误差 对干很复杂的蒙特卡洛方法计算问题 也是容易确定的 6 程序结构简单 易于实现 在计算机上进行蒙特卡洛方法计算时 程序 结构简单 分块性强 易于实现 蒙特卡洛的主要缺点有 1 收敛速度慢 如前所述 蒙特卡洛方法的收敛速度为 一般不 1 2 O N 容易得到精确度较高的近似结果 对于维数少 三维以下 的问题 不如其他 方法好 2 误差具有概率性 由于蒙特卡洛方法的误差是在一定置信水平下估计的 所以它的误差具有概率性 而不是一般意义下的误差 3 在粒子输运问题中 计算结果与系统大小有关 经验表明 只有当系统 的大小与粒子的平均自由程可以相比较时 一般在十个平均自由程左右 蒙特 卡洛方法计算的结果较为满意 但对于大系统或小概率事件的计算问题 计算 结果往往比真值偏低 而对于大系统 数值方法则是适用的 因此 在使用蒙 特卡洛方法时 可以考虑把蒙特卡洛方法与解析 或数值 方法相结合 取长 补短 2 42 4 蒙特卡洛方法的主要应用范围蒙特卡洛方法的主要应用范围 蒙特卡洛方法所特有的优点 使得它的应用范围越来越广 它的主要应用 范围包括 粒子输运问题 统计物理 典型数学问题 真空技术 激光技术以 及医学 生物 探矿等方面 随着科学技术的发展 其应用范围将更加广泛 蒙特卡洛方法在粒子输运问题中的应用范围主要包括 实验核物理 反应 堆物理 高能物理等方面 蒙特卡洛方法在实验核物理中的应用范围主要包括 8 通量及反应率 中子探测效率 光子探测效率 光子能量沉积谱及响应函数 气体正比计数管反冲质子谱 多次散射与通量衰减修正等方面 3 3 随机数随机数 随机数是蒙特卡洛方法的主要组成部分之一 随机数是指一个数列 其中 的每一个体称为随机数 其值与数列中的其它数无关 在一个均匀分布的随机 数中 每一个体出现的概率是均等的 物理中的很多过程需要随机数确定 比 如出射粒子的能量 方向等属性 粒子与介质的相互作用等等 所模拟的物理 过程要求随机数应具有下列特性 1 随机数序列应是独立的 互不相关的 uncorrelated 即序列中的任一子 序列应与其它的子序列无关 2 长的周期 long period 实际应用中 随机数都是用数学方法计算出来 的 这些算法具有周期性 即当序列达到一定长度后会重复 3 均匀分布的随机数应满足均匀性 Uniformity 随机数序列应是均匀的 无偏的 即 如果两个子区间的 面积 相等 则落于这两个子区间内的 随机数的个数应相等 4 有效性 Efficiency 模拟结果可靠 随机数的产生必须快速 有效 最 好能够进行并行计算 为了产生随机数 可以使用随机数表 随机数表是由0 1 9十个数字 组成 每个数字以0 1的等概率出现 数字之间相互独立 这些数字序列叫作随 机数字序列 如果要得到n位有效数字的随机数 只需将表中每n个相邻的随机 数字合并在一起 且在最高位的前边加上小数点即可 例如 某随机数表的第 一行数字为7634258910 要想得到三位有效数字的随机数依次为 0 763 0 425 0 891 可以使用物理方法产生随机数 用来作为随机数发生器的物理源主要有两 种 一种是根据放射性物质的放射性 另一种是利用计算机的固有噪声 但在 计算机上产生随机数最实用 最常见的方法是数学方法 即用如下递推公式 3 1 11 1 2 n knnn k Tn 产生随机数序列 对于给定的初始值 1 2 k 确定 n k 1 2 经常 使用的是k 1的情况 9 a 用数学方法产生的随机数有两个特点 即 递推公式和初始值 1 2 k 确定后 整个随机数序列便被唯一确定 不满足随机数相互独立的要求 b 由于随机数序列是由递推公式确定的 而在计算机上所能表示的 0 1 上 的数又是有限的 因此 这种方法产生的随机数序列就不可能不出现无 限重复 一旦出现这样的n n n n 使得 成立随机数序列便出现了周期性的循环现象 对1 2 nini ik 于k 1的情况 只要有一个随机数重复 其后面的随机数全部重复 这与 随机数的要求是不相符的 由于这两个问题的存在 常称用数学方法产 生的随机数为伪随机数 3 13 1 线性乘同余方法线性乘同余方法 线性乘同余方法是由Lehmer在1951年提出来的 是一种最常用的产生伪 随机数的方法 乘同余方法中采用的递推公式为 3 2 1 mod nn IaIcm 其中 I0为初始值 a为乘法器 c为增值 m为模数 mod为取模运算 除以m后的余数 a c m皆为整数 实型随机序列 n aIc 3 3 0 1 n nn I rIm float m 3 4 0 1 1 1 n nn I rIm float m 上式中 独立性和均匀性取决于参数a和c的选择 m应尽可能地大 因为序列 的周期不可能大于m 通常将m取为计算机所能表示的最大的整型量 在32位机 上 1961年 M Greenberger证明 用线性乘同余方法产生的 319 22 10m 随机数序列具有周期m的条件是 1 c和m互为质数 2 a 1是质数p的倍数 其 中p是a 1和m的公约数 3 如果m是4的倍数 a 1也是4的倍数 3 23 2 伪随机数序列的均匀性和独立性伪随机数序列的均匀性和独立性 3 2 1 伪随机数的均匀性 10 这里只考虑伪随机数序列 1 2 n全体作为子样时的均匀性问题 其中n 为伪随机数序列的最大容量 对于任意的0 x 1 令Nn x 表示伪随机数序列 1 2 n中适合不等式 i x i 1 2 n的个数 则 3 5 01 sup n x Nx nx n 将伪随机数序列 1 2 n从小至大重新排列 令 12n 01 0 1 n 则由 n 的定义 容易证明 很明显 对于固定 1 01 maxii x ii n nn 的 n 的值越小越好 它是描述伪随机数序列均匀程度的基本量 对于任 意随机数序列 均有不等式成立 当成立时 所对应的伪 1 2 n n 1 2 n n 随机数序列为最佳分布 3 2 2 伪随机数的独立性 对于任意 令表示 1 2 2 3 n n 1 中适合不等式0 1x y n Nx y 的个数 根据随机变量间相互独立的定义和频率近似概率的方 1 ii xy 法 令 3 6 0 1 sup nnn x y Nx yNx Ny n nnn 则 n 标志伪随机数序列 1 2 n的独立程度 简称为独立偏度 对于固定的 n n 的值越接近于零 伪随机数序列的独立性越好 4 4 蒙特卡洛方法在粒子输运上的应用蒙特卡洛方法在粒子输运上的应用 辐射 光子和中子 屏蔽问题是蒙特卡洛方法最早广泛应用的领域之一 现主要从物理直观出发 说明蒙特卡洛方法解决这类粒子输运问题的基本方法 和技巧 解决屏蔽问题时可采取多种方法 如直接模拟方法 简单加权法 统 计估计法 指数变换法等 这里只对直接模拟方法做介绍 4 14 1 屏蔽问题模型屏蔽问题模型 在反应堆工程和辐射的测量与应用中 常常要用一些吸收材料做成屏蔽物 挡住光子或中子 我们所关心的是经过屏蔽后射线的强度及其能量分布 这就 11 是屏蔽问题 当屏蔽物的形状复杂 散射各向异性 材料介质不均匀 核反应 截面与能量 位置有关时 难以用数值方法求解 用蒙特卡洛方法能够得到满 意的结果 粒子的输运问题带有明显的随机性质 粒子的输运过程是一个随机 过程 粒子的运动规律是根据大量粒子的运动状况总结出来的 是一种统计规 律 蒙特卡洛模拟 实际上就是模拟相当数量的粒子在介质中运动的状况 使 粒子运动的统计规律得以重现 不过 这种模拟不是用实验方法 而是利用数 值方法和技巧 即利用随机数来实现的 为方便起见 选用平板屏蔽模型 在厚度为 a 长 宽无限的平板左侧放置 一个强度已知 具有已知能量 方向分布的辐射源 S 见图 4 1 求粒子穿透屏 蔽概率 穿透率 及其能量 方向分布 穿透率就是由源发出的平均一个粒子 穿透屏蔽的数目 同时 假定粒子在两次碰撞之间按直线运动 且粒子之间的 相互作用可以忽略 图4 1 屏蔽问题模型 4 24 2 直接模拟方法直接模拟方法 直接模拟方法就是直接从物理问题出发 模拟粒子的真实物理过程 4 2 1 状态参数与状态序列 粒子在介质中的运动的状态 可用一组参数来描述 称之为状态参数 它 通常包括 粒子的空间位置r 能量E和运动方向 以S r E 表示 有时 还需要其他的参数 如粒子的时间t和附带的权重W 这时状态参数为 S r E t W 状态参数通常要根据所求问题的类型和所用的方法来确定 对于无 限平板几何 取S z E cos 其中z为粒子的位置坐标 为粒子的运动方向 与Z轴的夹角 对于球对称几何 取 S r E cos 其中r表示粒子所在位置 到球心的距离 为粒子的运动方向与其所在位置的径向夹角 粒子第m次碰撞后的状态参数为 或 mmmm E Sr 12 它表示一个由源发出的粒子 在介质中经过 m 次碰 mmmmmm EtW Sr 撞后的状态 其中 rm 粒子在第 m 次碰撞点的位置 Em 粒子第 m 次碰撞后的能量 m 粒子第 m 次碰撞后的运动方向 tm 粒子到第 m 次碰撞时所经历的时间 Wm 粒子第 m 次碰撞后的权重 一个由源发出的粒子在介质中运动 经过若干次碰撞后 直到其运动历史 结束 如逃出系统或被吸收等 假定粒子在两次碰撞之间按直线运动 其运动 方向与能量均不改变 则粒子在介质中的运动过程可用以下碰撞点的状态序列 描述 即S0 S1 SM 1 SM 或 011 011 011 MM MM MM EEEE rrrr 来描述 这里S0为粒子由源出发的状态 称为初态 SM 为粒子的终止状态 M 称为粒子运动的链长 4 2 2 模拟运动过程 这里以中子穿透均匀平板的模型来说明 这时状态参数取S z E cos 模 拟的步骤如下 1 确定初始状态S0 确定粒子的初始状态 实际上就是要从中子源的空间位置 能量和方向分 布中抽样 设源分布为 则分别从各自 000102030 cos cos f z Ef zfEf 的分布中抽样确定初始状态 2 确定下一个碰撞点 已知状态Sm 1 要确定状态Sm 首先要确定下一个碰撞点的位置zm 在相邻 两次碰撞之间 中子的输运长度l服从如下分布 4 1 111111 0 exp l tmmmtmmm f llElEdl r r 对于平板模型 l服从分布 13 4 2 111111 0 cos exp cos l tmmmtmmm f lzlEzlEdl 其中 t为介质的中子宏观总截面 积分称为粒子输 111 0 l tmmm lEdl r 运的自由程数 显然 粒子输运的自由程数服从指数分布 因此从f l 中抽样确 定l 就是要从积分方程中解出l 对于单一介 111 0 ln l tmmm lEdl r 质 则下一个碰撞点的位置为 1 ln tm l E 4 3 1111 1 ln coscos mmmmm tm zzlz E 如果zm a 则中子穿透屏蔽 若zm 0 则中子被反射出屏蔽 这两种情况 均视 为中子历史终止 3 确定被碰撞的原子核 通常介质由几种原子核组成 中子与核碰撞时 要确定与哪一种核碰撞 设介质由A B C 三种原子核组成 其核密度分别为NA NB NC 则介质的 宏观总截面为 4 4 1111 ABC tmtmtmtm EEEE 其中分别为核A B C的宏观总截面 其定义如下 ABC ttt 分别表示 核的宏观总截面 1 1 tmtm ENE 1 1 tmtm ENE 核密度和微观总截面 由于中子截面表示中子与核碰撞可能性的大小 因此 很自然地 中子与 A B C核发生碰撞的几率分别为 4 5 111 111 ABC tmtmtm ABC tmtmtm EEE PPP EEE 若 则中子与A核碰撞 若 则中子与B核碰撞 若 A P AB PP AB PP 则与C核碰撞 4 确定碰撞类型 确定了碰撞的核 比如B核 后 就要进一步确定碰撞类型 中子与核的反应 类型有弹性散射 非弹性散射 n 2n 反应 裂变和俘获等 它们的微观截面分 14 别为 则有 11 2 111 BBBBB elminmnnmfmcm EEEEE 和 4 6 111 2 111 BBBBBB tmelminmnnmfmcm EEEEEE 各种反应发生的几率分别为 4 7 11 11 2 2 11 11 11 BB elelmtm BB ininmtm BB nnnnmtm BB ffmtm BB ccmtm PEE PEE PEE PEE PEE 利用离散型随机变量的抽样方法 确定反应类型 在屏蔽问题中 中子与核反 应常只有弹性散射和吸收两种类型 吸收截面为 111 BBB amfmcm EEE 这时 总截面为 4 8 111 BBB tmelmam EEE 发生弹性散射的几率为 若 则为弹性散射 否则为吸收 1 1 B elm el B tm E P E el P 发生吸收反应意味着中子的历史终止 5 确定碰撞后的能量与运动方向 如果中子被碰撞核吸收 则其输运历史结束 如果发生弹性散射 需要确 定散射后中子的能量和运动方向 中子能量Em为 4 9 1 11cos 2 m mC E Err 其中 A是碰撞核的质量与中子质量之比 一般就取元素的原子量 21 1 A r A C为质心系中中子散射前后方向间的夹角 即偏转角 可从质心系cos CC 中弹性散射角分布fC C 中抽样产生 实验室系散射角 L的余弦 L为 2 1 12 C L C A AA 如果给出实验室系散射角余弦分布fL L 可直接从fL L 中抽取 L 此时能 量Em与 L的关系式为 确定了实验室系散射 2 22 1 2 1 1 m mLL E EA A 15 角 L后 再使用球面三角公式确定cos m 11 coscoscossinsincos mmLmL 各角度关系如图4 2所示 图4 2 角度关系示意图 至此 由Sm 1完全可以确定Sm 因此 当中子由源出发后 即S0确定后 重复步 骤 2 5 直到中子游动历史终止 于是得到了一个中子的随机游动历史 S0 S1 SM 1 SM 即 011 011 011 cos cos cos cos MM MM MM zzzz EEEE 也就是模拟了一个由源发出的中子的运动过程 4 2 3 记录结果 在获得中子的随机游动历史后 我们要对所要计算的物理量进行估计 对 于屏蔽问题 我们要计算中子的穿透率 考察每个中子的随机游动历史 它可 能穿透屏蔽 zM a 可能被屏蔽发射回来 zM 0 或者被吸收 设第n个中子 对穿透的贡献为 n 则 1 0 0 M n M za z 当 当 如果我们共跟踪了N 个中子 则穿透屏蔽的中子数为 则穿透屏 1 1 N n n N 蔽概率的近似值为 4 10 1 1 1 1 N Nn n N P NN 我们称这种直观地模拟过程和估计方法为直接模拟方法 在置信水平 1 0 95时 的误差为 1 N P 16 4 11 1 1 96 N PP NN 其中为 n的均方差 由于 n是一个服从二项分布的随机变量 所以 4 12 2 1 PP 或 4 13 2 1 1 1 NN PP 为得到中子穿透屏蔽的能量 角分布 将能量 角度范围分成若干个间隔 其中Emax Emin分别 min10maxI EEEEE 01 02 J 表示能量的上 下限 对于穿透屏蔽的中子按其能量 方向分间隔记录 设一 穿透屏蔽的中子能量为EM 其运动方向与Z轴夹角为 M 若能量EM属于第 i 个 能量间隔 Ei 角度 M属于第 j 个角度间隔 j 则分别在第 i 个能量计数器及 第 j 个角度计数器中加1 跟踪N个中子后 则 4 14 1 1 1 i i N i N P NE 1 2 iI 4 15 2 2 1 j j N j N P N 1 2 jJ 分别为穿透中子的能量分布和角分布 其中N1 i和N2 i分别为第 i 个能量和第 j 个角度间隔的穿透中子数 归一后分别为 4 16 1 1 1 1 1 1 1 i i N i N i N P N P NEP 1 2 iI 4 17 2 2 1 2 1 1 1 j j N j N j N P N P NP 1 2 jJ 4 34 3 蒙特卡洛方法的效率蒙特卡洛方法的效率 衡量蒙特卡洛技巧的好坏 除了看其方差大小外 还要看其所需费用 计 算时间 多少 即从该技巧的效率Ef 方差与费用乘积的倒数 全面考虑 4 18 2 1 f E T 17 其中 2为方差 T为所需费用 Ef大时 所用方法的效率高 否则 效率低 在一般情况下 直接模拟方法 简单加权法 统计估计法 指数变换法等方法 中有些方法虽然减小了方差 却增加了费用 例如 加权法 统计估计法虽然 较直接模拟方法减小了方差 却使每个粒子的运动链长增加 或记录贡献的计 算时间增加 因此 不能认为方差小的方法一定好 要从方法的效率全面考虑 在有些情况下 直接模拟方法仍然是一个被广泛使用的方法 5 5 蒙特卡洛方法应用程序蒙特卡洛方法应用程序 MCNPMCNP 5 15 1 MCNPMCNP 简述简述 MCNP A General Monte Carlo Code for Neutron and Particle Transport 是一套 通用的 三维空间中连续能量中子 光子和带电粒子 离子 联合输运过程模 拟程序 在军事和工业领域有着广泛应用 是基于蒙特卡洛方法的用于计算三 维复杂几何结构中的中子 光子 电子或者耦合中子 光子 电子输运问题的 通用软件包 也具有计算核临界系统 包括次临界和超临界系统 本征值问题 的能力 该软件包通过FORTRAN语言编程实现 MCNP程序具有超强的几何处理能力 几何系统由几何空间单元 cell 组成 而几何空间单元的界面 surface 由平面 二次曲面及特殊的四次椭圆环曲面组成 几何空间单元中的材料由包括同位素在内的多种核素组成 使用精确的点截面 参数 对特定的评价库 ENDF B IV V V VI库或ENDL851库 考虑了该库给出 的所有中子反应类型 在截面数据文件中收集了多种评价库的数据 对热中子 还配备了相应的截面数据 可按自由气体模型或S模型处理 对光子考虑了相干 和非相干散射 并处理了光电吸收后可能的荧光发射或电子对产生 MCNP3版 1983 年 和3A版 1985年 发行后 这一软件就成为用蒙特卡洛方 法模拟核过程最流行的通用程序 程序在计算辐射能量沉积和辐射计量等方面 取得成功 88年出版的 MCNP3B程序具有重复构造和结构的能力 能够解决特 征 谱线的问题 可以很好地模拟中子和光子的联合输运问题 使用的主要核数 据库是ENDF B 4 91年MCNP4版问世 这时程序可以模拟中子 光子 带电 粒子 离子 的联合输运过程 可以模拟探测器的测量结果 MCNP4版使用了更 新的ENDF B 6评价核数据库 加入了脉冲中子源功能等 MCNP5版 2003年 提高了彩色描点能力 64种颜色 提高了处理中性粒子照相问题能力 为源增 18 加了新选项 并对广泛应用的windows系统有了更好的支持 该程序是目前国 际上在核技术领域中应用最广泛 效果较佳 具有通用性的蒙特卡洛模拟计算 程序 许多核反应蒙特卡洛专用程序都引用该程序的核心部分 MCNP程序涉及面如此之多 关键是通过读入一个经用户创建的称为INP的 输入文件来进行计算 该文件必须遵循按照栅元卡 card 的格式进行组织 指定 描述空间问题的信息 具体地有 1 空间几何体的描述说明 2 几何体的使用 材料描述和交叉区域的选择估计 3 中子 光子以及电子这3种粒子源的位置 和特性说明 4 必要的回答卡和标记卡的类型 5 任何必需的冗余量消除技术 以提高计算效率 目前 MCNP以其灵活 通用的特点以及强大的功能被广泛 应用于辐射防护与射线测定 辐射屏蔽设计优化 反应堆设计 次 临界装置 实验 医学以及检测器设计与分析等学科领域 并得到一致认可 5 25 2 MCNPMCNP 误差的估计误差的估计 蒙特卡洛方法的结果方法的结果代表被抽样的许多历史过程贡献的平均值 假定P x 是选择一个随机步的几率密度函数 x是这个随机步产生的被估计的记 录值 其平均值记为 5 1 E xxP x dx E x 近似期望值可以通过MC方法得到 其中N是粒子数目 xi是 1 1 N i i xx N 从P x 中第i个历史的值 从加强大数定理 x的方差是离散度的度量 定义为 lim 1 n PxE x 为标准差 MCNP方法可以估计这 2 222 xE xP x dxE xE x 个值