《蒙特卡洛应用讲义 完整版》由会员分享，可在线阅读，更多相关《蒙特卡洛应用讲义 完整版（400页珍藏版）》请在金锄头文库上搜索。
1、蒙特卡罗方法蒙特卡罗方法在核技术中的应用在核技术中的应用林谦目目 录录第一章第一章 蒙特卡罗方法概述蒙特卡罗方法概述第二章第二章 随机数随机数第三章第三章 由已知分布的随机抽样由已知分布的随机抽样第四章第四章 蒙特卡罗方法解粒子输运问题蒙特卡罗方法解粒子输运问题教材教材l蒙特卡罗方法蒙特卡罗方法在在实验核物理中的应用实验核物理中的应用许淑艳编著原子能出版社l蒙特卡罗方法蒙特卡罗方法清华大学参考书参考书l蒙特卡罗方法蒙特卡罗方法及其在粒子输运问题中的应用及其在粒子输运问题中的应用裴鹿成张孝泽编著 科学出版社l蒙特卡罗方法蒙特卡罗方法徐钟济编著上海科学技术出版社联系方式联系方式l电话电话83918
2、l电子邮件电子邮件第一章第一章 蒙特卡罗方法概述蒙特卡罗方法概述1.蒙特卡罗方法的基本思想蒙特卡罗方法的基本思想2.蒙特卡罗方法的收敛性，误差蒙特卡罗方法的收敛性，误差3.蒙特卡罗方法的特点蒙特卡罗方法的特点4.蒙特卡罗方法的主要应用范围蒙特卡罗方法的主要应用范围作作 业业第一章第一章 蒙特卡罗方法概述蒙特卡罗方法概述蒙特卡罗方法又称随机抽样技巧或统计试验方法。半个多世纪以来，由于科学技术的发展和电子计算机的发明，这种方法作为一种独立的方法被提出来，并首先在核武器的试验与研制中得到了应用。蒙特卡罗方法是一种计算方法，但与一般数值计算方法有很大区别。它是以概率统计理论为基础的一种方法。由于蒙特卡
3、罗方法能够比较逼真地描述事物的特点及物理实验过程，解决一些数值方法难以解决的问题，因而该方法的应用领域日趋广泛。1.蒙特卡罗方法的基本思想蒙特卡罗方法的基本思想二十世纪四十年代中期，由于科学技术的发展和电子计算机的发明，蒙特卡罗方法作为一种独立的方法被提出来，并首先在核武器的试验与研制中得到了应用。但其基本思想并非新颖，人们在生产实践和科学试验中就已发现，并加以利用。两个例子例1.蒲丰氏问题例2.射击问题（打靶游戏）基本思想计算机模拟试验过程例1. 蒲丰氏问题为了求得圆周率值，在十九世纪后期，有很多人作了这样的试验：将长为2l的一根针任意投到地面上，用针与一组相间距离为2a（ la）的平行线相
4、交的频率代替概率P，再利用准确的关系式：求出值其中为投计次数，n为针与平行线相交次数。这就是古典概率论中著名的蒲丰氏问题。一些人进行了实验，其结果列于下表：实验者年份投计次数的实验值沃尔弗(Wolf)185050003.1596斯密思(Smith)185532043.1553福克斯(Fox)189411203.1419拉查里尼(Lazzarini)190134083.1415929例2. 射击问题（打靶游戏） 设r表示射击运动员的弹着点到靶心的距离，(r)表示击中r处相应的得分数（环数），f(r)为该运动员的弹着点的分布密度函数，它反映运动员的射击水平。该运动员的射击成绩为用概率语言来说，是随
5、机变量(r)的数学期望，即现假设该运动员进行了次射击，每次射击的弹着点依次为r1，r2，rN，则次得分g(r1)，g(r2)，g(rN)的算术平均值代表了该运动员的成绩。换言之，为积分的估计值，或近似值。在该例中，用次试验所得成绩的算术平均值作为数学期望的估计值（积分近似值）。基本思想 由以上两个例子可以看出，当所求问题的解是某个事件的概率，或者是某个随机变量的数学期望，或者是与概率、数学期望有关的量时，通过某种试验的方法，得出该事件发生的频率，或者该随机变量若干个具体观察值的算术平均值，通过它得到问题的解。这就是蒙特卡罗方法的基本思想。当随机变量的取值仅为1或0时，它的数学期望就是某个事件的
6、概率。或者说，某种事件的概率也是随机变量（仅取值为1或0）的数学期望。因此，可以通俗地说，蒙特卡罗方法是用随机试验的方法计算积分，即将所要计算的积分看作服从某种分布密度函数f(r)的随机变量(r)的数学期望通过某种试验，得到个观察值r1，r2，rN（用概率语言来说，从分布密度函数f(r)中抽取个子样r1，r2，rN，），将相应的个随机变量的值g(r1)，g(r2)，g(rN)的算术平均值作为积分的估计值（近似值）。为了得到具有一定精确度的近似解，所需试验的次数是很多的，通过人工方法作大量的试验相当困难，甚至是不可能的。因此，蒙特卡罗方法的基本思想虽然早已被人们提出，却很少被使用。本世纪四十年代
7、以来，由于电子计算机的出现，使得人们可以通过电子计算机来模拟随机试验过程，把巨大数目的随机试验交由计算机完成，使得蒙特卡罗方法得以广泛地应用，在现代化的科学技术中发挥应有的作用。计算机模拟试验过程 计算机模拟试验过程，就是将试验过程（如投针，射击）化为数学问题，在计算机上实现。以上述两个问题为例，分别加以说明。例1.蒲丰氏问题例2.射击问题（打靶游戏）由上面两个例题看出，蒙特卡罗方法常以一个“概率模型”为基础，按照它所描述的过程，使用由已知分布抽样的方法，得到部分试验结果的观察值，求得问题的近似解。例蒲丰氏问题设针投到地面上的位置可以用一组参数（x,）来描述，x为针中心的坐标，为针与平行线的夹
8、角，如图所示。任意投针，就是意味着x与都是任意取的，但x的范围限于0，a，夹角的范围限于0，。在此情况下，针与平行线相交的数学条件是针在平行线间的位置如何产生任意的（x,）？x在0，a上任意取值，表示x在0，a上是均匀分布的，其分布密度函数为：类似地，的分布密度函数为：因此，产生任意的（x,）的过程就变成了由f1(x)抽样x及由f2()抽样的过程了。由此得到：其中1，2均为（0,1）上均匀分布的随机变量。每次投针试验，实际上变成在计算机上从两个均匀分布的随机变量中抽样得到（x,），然后定义描述针与平行线相交状况的随机变量s(x,)，为如果投针次，则是针与平行线相交概率的估计值。事实上，于是有例
9、射击问题 设射击运动员的弹着点分布为用计算机作随机试验（射击）的方法为，选取一个随机数，按右边所列方法判断得到成绩。这样，就进行了一次随机试验（射击），得到了一次成绩(r)，作次试验后，得到该运动员射击成绩的近似值环数78910概率0.10.10.30.52.蒙特卡罗方法的收敛性，误差蒙特卡罗方法的收敛性，误差 蒙特卡罗方法作为一种计算方法，其收敛性与误差是普遍关心的一个重要问题。收敛性误差减小方差的各种技巧效率收敛性 由前面介绍可知，蒙特卡罗方法是由随机变量X的简单子样X1，X2，XN的算术平均值：作为所求解的近似值。由大数定律可知，如X1，X2，XN独立同分布，且具有有限期望值（E(X)）
10、，则即随机变量X的简单子样的算术平均值，当子样数充分大时，以概率1收敛于它的期望值E(X)。误差蒙特卡罗方法的近似值与真值的误差问题，概率论的中心极限定理给出了答案。该定理指出，如果随机变量序列X1，X2，XN独立同分布，且具有有限非零的方差2，即 f(X)是X的分布密度函数。则当N充分大时，有如下的近似式其中称为置信度，1称为置信水平。这表明，不等式近似地以概率1成立，且误差收敛速度的阶为。通常，蒙特卡罗方法的误差定义为上式中与置信度是一一对应的，根据问题的要求确定出置信水平后，查标准正态分布表，就可以确定出。下面给出几个常用的与的数值：关于蒙特卡罗方法的误差需说明两点：第一，蒙特卡罗方法的
11、误差为概率误差，这与其他数值计算方法是有区别的。第二，误差中的均方差是未知的，必须使用其估计值来代替，在计算所求量的同时，可计算出。0.50.050.0030.67451.963减小方差的各种技巧 显然，当给定置信度后，误差由和N决定。要减小，或者是增大N，或者是减小方差2。在固定的情况下，要把精度提高一个数量级，试验次数N需增加两个数量级。因此，单纯增大N不是一个有效的办法。另一方面，如能减小估计的均方差，比如降低一半，那误差就减小一半，这相当于N增大四倍的效果。因此降低方差的各种技巧，引起了人们的普遍注意。后面课程将会介绍一些降低方差的技巧。效率一般来说，降低方差的技巧，往往会使观察一个子
12、样的时间增加。在固定时间内，使观察的样本数减少。所以，一种方法的优劣，需要由方差和观察一个子样的费用（使用计算机的时间）两者来衡量。这就是蒙特卡罗方法中效率的概念。它定义为，其中c 是观察一个子样的平均费用。显然越小，方法越有效。3.蒙特卡罗方法的特点蒙特卡罗方法的特点优点1)能够比较逼真地描述具有随机性质的事物的特点及物理实验过程。2)受几何条件限制小。3)收敛速度与问题的维数无关。4)具有同时计算多个方案与多个未知量的能力。5)误差容易确定。6)程序结构简单，易于实现。缺点1)收敛速度慢。2)误差具有概率性。3)在粒子输运问题中，计算结果与系统大小有关。1)能够比较逼真地描述具有随机性质的
13、事物的特点及物理实验过程从这个意义上讲，蒙特卡罗方法可以部分代替物理实验，甚至可以得到物理实验难以得到的结果。用蒙特卡罗方法解决实际问题，可以直接从实际问题本身出发，而不从方程或数学表达式出发。它有直观、形象的特点。2)受几何条件限制小在计算s维空间中的任一区域Ds上的积分时，无论区域Ds的形状多么特殊，只要能给出描述Ds的几何特征的条件，就可以从Ds中均匀产生N个点，得到积分的近似值。其中Ds为区域Ds的体积。这是数值方法难以作到的。另外，在具有随机性质的问题中，如考虑的系统形状很复杂，难以用一般数值方法求解，而使用蒙特卡罗方法，不会有原则上的困难。3)收敛速度与问题的维数无关由误差定义可知
14、，在给定置信水平情况下，蒙特卡罗方法的收敛速度为，与问题本身的维数无关。维数的变化，只引起抽样时间及估计量计算时间的变化，不影响误差。也就是说，使用蒙特卡罗方法时，抽取的子样总数N与维数s无关。维数的增加，除了增加相应的计算量外，不影响问题的误差。这一特点，决定了蒙特卡罗方法对多维问题的适应性。而一般数值方法，比如计算定积分时，计算时间随维数的幂次方而增加，而且，由于分点数与维数的幂次方成正比，需占用相当数量的计算机内存，这些都是一般数值方法计算高维积分时难以克服的问题。4)具有同时计算多个方案与多个未知量的能力对于那些需要计算多个方案的问题，使用蒙特卡罗方法有时不需要像常规方法那样逐个计算，
15、而可以同时计算所有的方案，其全部计算量几乎与计算一个方案的计算量相当。例如，对于屏蔽层为均匀介质的平板几何，要计算若干种厚度的穿透概率时，只需计算最厚的一种情况，其他厚度的穿透概率在计算最厚一种情况时稍加处理便可同时得到。另外，使用蒙特卡罗方法还可以同时得到若干个所求量。例如，在模拟粒子过程中，可以同时得到不同区域的通量、能谱、角分布等，而不像常规方法那样，需要逐一计算所求量。5)误差容易确定对于一般计算方法，要给出计算结果与真值的误差并不是一件容易的事情，而蒙特卡罗方法则不然。根据蒙特卡罗方法的误差公式，可以在计算所求量的同时计算出误差。对干很复杂的蒙特卡罗方法计算问题，也是容易确定的。一般
16、计算方法常存在着有效位数损失问题，而要解决这一问题有时相当困难，蒙特卡罗方法则不存在这一问题。6)程序结构简单，易于实现在计算机上进行蒙特卡罗方法计算时，程序结构简单，分块性强，易于实现。1)收敛速度慢如前所述，蒙特卡罗方法的收敛速度为，一般不容易得到精确度较高的近似结果。对于维数少（三维以下）的问题，不如其他方法好。2)误差具有概率性由于蒙特卡罗方法的误差是在一定置信水平下估计的，所以它的误差具有概率性，而不是一般意义下的误差。3)在粒子输运问题中，计算结果与系统大小有关经验表明，只有当系统的大小与粒子的平均自由程可以相比较时（一般在十个平均自由程左右），蒙特卡罗方法计算的结果较为满意。但对
17、于大系统或小概率事件的计算问题，计算结果往往比真值偏低。而对于大系统，数值方法则是适用的。因此，在使用蒙特卡罗方法时，可以考虑把蒙特卡罗方法与解析（或数值）方法相结合，取长补短，既能解决解析（或数值）方法难以解决的问题，也可以解决单纯使用蒙特卡罗方法难以解决的问题。这样，可以发挥蒙特卡罗方法的特长，使其应用范围更加广泛。4.蒙特卡罗方法的主要应用范围蒙特卡罗方法的主要应用范围 蒙特卡罗方法所特有的优点，使得它的应用范围越来越广。它的主要应用范围包括：粒子输运问题，统计物理，典型数学问题，真空技术，激光技术以及医学，生物，探矿等方面。随着科学技术的发展，其应用范围将更加广泛。蒙特卡罗方法在粒子输
18、运问题中的应用范围主要包括：实验核物理，反应堆物理，高能物理等方面。蒙特卡罗方法在实验核物理中的应用范围主要包括：通量及反应率，中子探测效率，光子探测效率，光子能量沉积谱及响应函数，气体正比计数管反冲质子谱，多次散射与通量衰减修正等方面。作作 业业 1)用蒲丰投针法在计算机上计算值，取a=4、l=3。2)分别用理论计算和计算机模拟计算，求连续掷两颗骰子，点数之和大于6且第一次掷出的点数大于第二次掷出点数的概率。第二章第二章 随机数随机数1.随机数的定义及产生方法随机数的定义及产生方法 2.伪随机数伪随机数3.产生伪随机数的乘同余方法产生伪随机数的乘同余方法4.产生伪随机数的乘加同余方法产生伪随
19、机数的乘加同余方法5.产生伪随机数的其他方法产生伪随机数的其他方法6.伪随机数序列的均匀性和独立性伪随机数序列的均匀性和独立性作作 业业第二章第二章 随机数随机数由具有已知分布的总体中抽取简单子样，在蒙特卡罗方法中占有非常重要的地位。总体和子样的关系，属于一般和个别的关系，或者说属于共性和个性的关系。由具有已知分布的总体中产生简单子样，就是由简单子样中若干个性近似地反映总体的共性。随机数是实现由已知分布抽样的基本量，在由已知分布的抽样过程中，将随机数作为已知量，用适当的数学方法可以由它产生具有任意已知分布的简单子样。1.随机数的定义及产生方法随机数的定义及产生方法1)随机数的定义及性质2)随机
20、数表3)物理方法1)随机数的定义及性质 在连续型随机变量的分布中，最简单而且最基本的分布是单位均匀分布。由该分布抽取的简单子样称，随机数序列，其中每一个体称为随机数。单位均匀分布也称为0，1上的均匀分布，其分布密度函数为：分布函数为：由于随机数在蒙特卡罗方法中占有极其重要的位置，我们用专门的符号表示。由随机数序列的定义可知，1，2，是相互独立且具有相同单位均匀分布的随机数序列。也就是说，独立性独立性、均匀性均匀性是随机数必备的两个特点。随机数具有非常重要的性质：对于任意自然数s，由s个随机数组成的s维空间上的点（n+1，n+2，n+s）在s维空间的单位立方体Gs上均匀分布，即对任意的ai，如下
21、等式成立：其中P()表示事件发生的概率。反之，如果随机变量序列1, 2对于任意自然数s，由s个元素所组成的s维空间上的点（n+1，n+s）在Gs上均匀分布，则它们是随机数序列。由于随机数在蒙特卡罗方法中所处的特殊地位，它们虽然也属于由具有已知分布的总体中产生简单子样的问题，但就产生方法而言，却有着本质上的差别。2)随机数表为了产生随机数，可以使用随机数表。随机数表是由0，1，9十个数字组成，每个数字以0.1的等概率出现，数字之间相互独立。这些数字序列叫作随机数字序列。如果要得到n位有效数字的随机数，只需将表中每n个相邻的随机数字合并在一起，且在最高位的前边加上小数点即可。例如，某随机数表的第一
22、行数字为7634258910，要想得到三位有效数字的随机数依次为0.763，0.425，0.891。因为随机数表需在计算机中占有很大内存，而且也难以满足蒙特卡罗方法对随机数需要量非常大的要求，因此，该方法不适于在计算机上使用。3)物理方法用物理方法产生随机数的基本原理是：利用某些物理现象，在计算机上增加些特殊设备，可以在计算机上直接产生随机数。这些特殊设备称为随机数发生器。用来作为随机数发生器的物理源主要有两种：一种是根据放射性物质的放射性，另一种是利用计算机的固有噪声。一般情况下，任意一个随机数在计算机内总是用二进制的数表示的：其中i（i=1，2，m）或者为0，或者为1。因此，利用物理方法在
23、计算机上产生随机数，就是要产生只取0或1的随机数字序列，数字之间相互独立，每个数字取0或1的概率均为0.5。用物理方法产生的随机数序列无法重复实现，不能进行程序复算，给验证结果带来很大困难。而且，需要增加随机数发生器和电路联系等附加设备，费用昂贵。因此，该方法也不适合在计算机上使用。2.伪随机数伪随机数1)伪随机数2)伪随机数存在的两个问题3)伪随机数的周期和最大容量1)伪随机数在计算机上产生随机数最实用、最常见的方法是数学方法，即用如下递推公式：产生随机数序列。对于给定的初始值1,2，k，确定n+k，=1，2，。经常使用的是k=1的情况，其递推公式为：对于给定的初始值1，确定n+1，=,2)
24、伪随机数存在的两个问题用数学方法产生的随机数，存在两个问题：a)递推公式和初始值1,2，k确定后，整个随机数序列便被唯一确定。不满足随机数相互独立的要求。b)由于随机数序列是由递推公式确定的，而在计算机上所能表示的0，1上的数又是有限的，因此，这种方法产生的随机数序列就不可能不出现无限重复。一旦出现这样的n，n(n n)，使得下面等式成立：随机数序列便出现了周期性的循环现象。对于k=1的情况，只要有一个随机数重复，其后面的随机数全部重复，这与随机数的要求是不相符的。由于这两个问题的存在，常称用数学方法产生的随机数为伪随机数。对于以上存在的两个问题，作如下具体分析。关于第一个问题，不能从本质上加
25、以改变，但只要递推公式选得比较好，随机数间的相互独立性是可以近似满足的。至于第二个问题，则不是本质的。因为用蒙特卡罗方法解任何具体问题时，所使用的随机数的个数总是有限的，只要所用随机数的个数不超过伪随机数序列出现循环现象时的长度就可以了。用数学方法产生的伪随机数容易在计算机上得到，可以进行复算，而且不受计算机型号的限制。因此，这种方法虽然存在着一些问题，但仍然被广泛地在计算机上使用，是在计算机上产生伪随机数的主要方法。3)伪随机数的周期和最大容量 发生周期性循环现象的伪随机数的个数称为伪随机数的周期。对于前面介绍的情况，伪随机数的周期为nn。从伪随机数序列的初始值开始，到出现循环现象为止，所产
26、生的伪随机数的个数称为伪随机数的最大容量。前面的例子中，伪随机数的最大容量为n。3.产生伪随机数的乘同余方法产生伪随机数的乘同余方法乘同余方法是由Lehmer在1951年提出来的，它的一般形式是：对于任一初始值x1，伪随机数序列由下面递推公式确定：其中a为常数。1)乘同余方法的最大容量的上界 对于任意正整数M，根据数论中的标准分解定理，总可以分解成如下形式： 其中P0=2，P1， Pr表示不同的奇素数，0表示非负整数，1，r表示正整数。a无论取什么值，乘同余方法的最大容量的上界为： 的最小公倍数。其中：2)关于a与x1的取值 如果a与x1满足如下条件： 对于 ， x1与M互素，则乘同余方法产生
27、的伪随机数序列的最大容量达到最大可能值(M)。3)乘同余方法在计算机上的使用为了便于在计算机上使用，通常取：=2s其中s为计算机中二进制数的最大可能有效位数x1=奇数a=52k+1其中k为使52k+1在计算机上所能容纳的最大整数，即a为计算机上所能容纳的5的最大奇次幂。一般地，s=32时，a=513；s=48，a=515等。伪随机数序列的最大容量(M)=2s-2。 乘同余方法是使用的最多、最广的方法，在计算机上被广泛地使用。4.产生伪随机数的乘加同余方法产生伪随机数的乘加同余方法 产生伪随机数的乘加同余方法是由Rotenberg于1960年提出来的，由于这个方法有很多优点，已成为仅次于乘同余方
28、法产生伪随机数的另一主要方法。 乘加同余方法的一般形式是，对任意初始值x1，伪随机数序列由下面递推公式确定：其中a和c为常数。1)乘加同余方法的最大容量 关于乘加同余方法的最大容量问题，有如下结论：如果对于正整数M的所有素数因子P，下式均成立： 当M为4的倍数时，还有下式成立： c与M互素，则乘加同余方法所产生的伪随机数序列的最大容量达到最大可能值M。2)M，x1，a，c的取值 为了便于在计算机上使用，通常取M=2s 其中s为计算机中二进制数的最大可能有效位数。a=2b+1(b2)c =1 这样在计算中可以使用移位和指令加法，提高计算速度。 5.产生伪随机数的其他方法产生伪随机数的其他方法1)
29、取中方法2)加同余方法6.伪随机数序列的均匀性和独立性伪随机数序列的均匀性和独立性判断伪随机数序列是否满足均匀和相互独立的要求，要靠统计检验的方法实现。对于伪随机数的统计检验，一般包括两大类：均匀性检验和独立性检验。六十年代初，人们开始用定性的方法研究伪随机数序列的均匀性和独立性问题，简要叙述如下。1)伪随机数的均匀性 这里只考虑伪随机数序列1,2，n全体作为子样时的均匀性问题。其中n为伪随机数序列的最大容量。对 于 任 意 的 0x1， 令 Nn(x)表 示 伪 随 机 数 序 列1,2，n中适合不等式i0。对该分布的直接抽样方法如下：例例3. 掷骰子点数的抽样掷骰子点数的抽样 掷骰子点数X
30、=n的概率为：选取随机数，如则在等概率的情况下，可使用如下更简单的方法：其中表示取整数。例例4. 碰撞核种类的确定碰撞核种类的确定 中子或光子在介质中发生碰撞时，如介质是由多种元素组成，需要确定碰撞核的种类。假定介质中每种核的宏观总截面分别为1，2，n，则中子或光子与每种核碰撞的概率分别为：其中t12n。碰撞核种类的确定方法为：产生一个随机数，如果则中子或光子与第I种核发生碰撞。例例5. 中子与核的反应类型的确定中子与核的反应类型的确定 假设中子与核的反应类型有如下几种:弹性散射，非弹性散射，裂变，吸收，相应的反应截面分别为el，in，f，a。则发生每一种反应类型的概率依次为：其中反应总截面t
31、elinfa。反应类型的确定方法为：产生一个随机数2)连续型分布的直接抽样方法连续型分布的直接抽样方法 对于连续型分布，如果分布函数F(x)的反函数 F1(x)存在，则直接抽样方法是：例例6. 在在a，b上均匀分布的抽样上均匀分布的抽样 在a，b上均匀分布的分布函数为：则例例7. 分布分布为连续型分布，作为它的一个特例是：其分布函数为：则例例8. 指数分布指数分布为连续型分布，其一般形式如下：其分布函数为：则因为1也是随机数，可将上式简化为连续性分布函数的直接抽样方法对于分布函数的反函数存在且容易实现的情况，使用起来是很方便的。但是对于以下几种情况，直接抽样法是不合适的。1)分布函数无法用解析
32、形式给出，因而其反函数也无法给出。2)分布函数可以给出其解析形式，但是反函数给不出来。3)分布函数即使能够给出反函数，但运算量很大。下面叙述的挑选抽样方法是克服这些困难的比较好的方法。3.挑选抽样方法挑选抽样方法 为了实现从己知分布密度函数f(x)抽样，选取与f(x)取值范围相同的分布密度函数h(x)，如果则挑选抽样方法为：即从h(x)中抽样xh，以的概率接受它。下面证明xf 服从分布密度函数f(x)。证明：对于任意x使用挑选抽样方法时，要注意以下两点：选取h(x)时要使得h(x)容易抽样且M的值要尽量小。因为M小能提高抽样效率。抽样效率是指在挑选抽样方法中进行挑选时被选中的概率。按此定义，该
33、方法的抽样效率E为：所以，M越小，抽样效率越高。当 f(x)在0，1上定义时，取 h(x)=1，Xh=， 此时挑选抽样方法为例例9. 圆内均匀分布抽样令圆半径为R0，点到圆心的距离为r，则r的分布密度函数为分布函数为容易知道，该分布的直接抽样方法是由于开方运算在计算机上很费时间，该方法不是好方法。下面使用挑选抽样方法：取则抽样框图为显 然 ， 没 有 必 要 舍 弃 1 2的 情 况 ， 此 时 ， 只 需 取 就可以了，亦即另一方面，也可证明与具有相同的分布。4.复合抽样方法复合抽样方法 在实际问题中，经常有这样的随机变量，它服从的分布与一个参数有关，而该参数也是一个服从确定分布的随机变量，
34、称这样的随机变量服从复合分布。例如，分布密度函数是一个复合分布。其中Pn0，n=1，2，且 fn(x)为与参数n有关的分布密度函数，n=1，2，参数n服从如下分布复合分布的一般形式为：其中f2(x/y)表示与参数y有关的条件分布密度函数，F1(y)表示分布函数。复合分布的抽样方法为:首先由分布函数F1(y)或分布密度函数f1(y)中抽样YF1或Yf1，然后再由分布密度函数f2(x/YF1)中抽样确定Xf2(x/YF)证明：所以，Xf所服从的分布为f(x)。例例10. 指数函数分布的抽样指数函数分布的一般形式为：引入如下两个分布密度函数：则使用复合抽样方法，首先从f1(y)中抽取y再由f2(x/
35、YF1)中抽取x5.复合挑选抽样方法复合挑选抽样方法 考虑另一种形式的复合分布如下：其中0H(x,y)M，f2(x/y)表示与参数y有关的条件分布密度函数，F1(y)表示分布函数。抽样方法如下：证明：抽样效率为：E=1/M为了实现某个复杂的随机变量y 的抽样，将其表示成若干个简单的随机变量x1，x2，xn的函数得到x1，x2，xn 的抽样后，即可确定y 的抽样，这种方法叫作替换法抽样。即6.替换抽样方法例例11. 散射方位角余弦分布的抽样散射方位角在0,2上均匀分布，则其正弦和余弦sin和cos服从如下分布：直接抽样方法为：令=2，则在0,上均匀分布，作变换其中01，0，则(x,y)表示上半个
36、单位圆内的点。如果(x,y)在上半个单位圆内均匀分布，则在0,上均匀分布，由于因此抽样sin和cos的问题就变成在上半个单位圆内均匀抽样(x,y)的问题。为获得上半个单位圆内的均匀点，采用挑选法，在上半个单位圆的外切矩形内均匀投点（如图）。舍弃圆外的点，余下的就是所要求的点。抽样方法为：抽样效率E=/40.785为实现散射方位角余弦分布抽样，最重要的是在上半个单位圆内产生均匀分布点。下面这种方法，首先在单位圆的半个外切正六边形内产生均匀分布点，如图所示。于是便有了抽样效率更高的抽样方法：抽样效率例例12. 正态分布的抽样标准正态分布密度函数为：引入一个与标准正态随机变量X独立同分布的随机变量Y
37、，则（X,Y）的联合分布密度为：作变换则（,）的联合分布密度函数为：由此可知，与相互独立，其分布密度函数分别为分别抽取，：从而得到一对服从标准正态分布的随机变量X和Y：对于一般的正态分布密度函数N(,2)的抽样，其抽样结果为：例例13. 分布的抽样分布密度函数的一般形式为：其中n，k为整数。为了实现分布的抽样，将其看作一组简单的相互独立随机变量的函数，通过这些简单随机变量的抽样，实现分布的抽样。设x1，x2，xn为一组相互独立、具有相同分布F(x)的随机变量，k为x1，x2，xn按大小顺序排列后的第k个，记为：则k的分布函数为：当F(x)=x时，不难验证，k的分布密度函数为分布。因此，分布的抽
38、样可用如下方法实现：选取n个随机数，按大小顺序排列后取第k个，即7.随机抽样的一般方法随机抽样的一般方法 1)加抽样方法2)减抽样方法3)乘抽样方法4)乘加抽样方法5)乘减抽样方法6)对称抽样方法7)积分抽样方法1)加抽样方法加抽样方法是对如下加分布给出的一种抽样方法：其中Pn0， ，且 fn(x)为与参数n有关的分布密度函数，n=1，2，。由复合分布抽样方法可知，加分布的抽样方法为:首先抽样确定n，然后由fn(x)中抽样x，即：例例14. 多项式分布抽样多项式分布密度函数的一般形式为：将f(x)改写成如下形式：则该分布的抽样方法为：例例15. 球壳内均匀分布抽样设球壳内半径为R0，外半径为R
39、1，点到球心的距离为r，则r的分布密度函数为分布函数为该分布的直接抽样方法是为避免开立方根运算，作变换：则x0,1，其分布密度函数为：其中则x及r的抽样方法为：2)减抽样方法减抽样方法是对如下形式的分布密度所给出的一种抽样方法：其中A1、A2为非负实数，f1(x) 、f2(x)均为分布密度函数。减抽样方法分为以下两种形式：以上两种形式的抽样方法，究竟选择哪种好，要看f1(x) 、f2(x)哪一个容易抽样，如相差不多，选用第一种方法抽样效率高。（1）将f (x)表示为令m表示f2(x)f1(x)的下界，使用挑选法，从f1(x)中抽取Xf1抽样效率为：（2）将f (x)表示为使用挑选法，从f2(x
40、)中抽取Xf2抽样效率为：例例16. 分布抽样分布的一个特例：取A12，A21，f1(x)1，f2(x)2x，此时m0，则根据第一种形式的减抽样方法，有或由于11可用1代替，该抽样方法可简化为：对于21的情况，可取Xf1，因此与分布的推论相同。如下形式的分布称为乘分布：其中H(x)为非负函数， f1(x)为任意分布密度函数。令M为H(x)的上界，乘抽样方法如下：抽样效率为：3)乘抽样方法例例17. 倒数分布抽样倒数分布密度函数为：其直接抽样方法为：下面采用乘抽样方法，考虑如下分布族：其中i =1，2，该分布的直接抽样方法为：利用这一分布族，将倒数分布f(x)表示成:其中，乘法分布的抽样方法如下
41、:该分布的抽样效率为：例例18. 麦克斯韦(Maxwell)分布抽样麦克斯韦分布密度函数的一般形式为：使用乘抽样方法，令该分布的直接抽样方法为：此时则麦克斯韦分布的抽样方法为:该分布的抽样效率为：在实际问题中，经常会遇到如下形式的分布：其中Hn(x)为非负函数，fn(x)为任意分布密度函数，n=1，2，。不失一般性，只考虑n=2的情况：将f(x)改写成如下的加分布形式：4)乘加抽样方法其中乘加抽样方法为：该方法的抽样效率为：这种方法需要知道P1的值(P2=1P1），这对有些分布是很困难的。下面的方法可以不用计算P1：对于任意小于1的正数P1，令P2=1P1；则采用复合挑选抽样方法，有：当取时，
42、抽样效率最高这时，乘加抽样方法为：由于可知第一种方法比第二种方法的抽样效率高。例例19. 光子散射后能量分布的抽样令光子散射前后的能量分别为 和（以m0c2为单位，m0为电子静止质量，c为光速），则x的分布密度函数为：该分布即为光子散射能量分布，它是由著名的KlinNishina公式确定的。其中K()为归一因子：把光子散射能量分布改写成如下形式：在1,1+2上定义如下函数:则有使用乘加抽样方法：光子散射能量分布的抽样方法为：该方法的抽样效率为：乘减分布的形式为：其中H1(x)、H2(x)为非负函数，f1(x)、f2(x)为任意分布密度函数。与减抽样方法类似，乘减分布的抽样方法也分为两种。5)乘
43、减抽样方法（1）将f (x)表示为令H1(x)的上界为M1，的下界为m，使用乘抽样方法得到如下乘减抽样方法：（2）将f (x)表示为令H2(x)的上界为M2，使用乘抽样方法，得到另一种乘减抽样方法：例例20. 裂变中子谱分布抽样裂变中子谱分布的一般形式为：其中A，B，C，Emin，Emax均为与元素有关的量。令其中为归一因子，为任意参数。相应的H1(E)，H2(E)为：于是裂变中子谱分布可以表示成乘减分布形式:容易确定H1(E)的上界为：为提高抽样效率，应取使得M1达到最小，此时取m0，令则裂变中子谱分布的抽样方法为:抽样效率对称分布的一般形式为：其中f1(x)为任意分布密度函数，满足偶函数对
44、称条件，H(x)为任意奇函数，即对任意x满足：对称分布的抽样方法如下：取=216)对称抽样方法证明：因为=21，x相当于，因此例例21. 质心系各向同性散射角余弦分布抽样在质心系各向同性散射的假设下，为得到实验室系散射角余弦，需首先抽样确定质心条散射角余弦：再利用下面转换公式:得到实验室系散射角余弦L。其中A为碰撞核质量，C、L分别为质心系和实验室系散射角。为避免开方运算，可以使用对称分布抽样。根据转换公式可得：依照质心系散射各向同性的假定，可得到实验室系散射角余弦L的分布如下：该密度函数中的第一项为偶函数，第二项为奇函数，因而是对称分布。其中从f1(L)的抽样可使用挑选法然后再以的概率决定接
45、受或取负值。上述公式涉及开方运算，需要进一步简化。注意以下事实：对于任意0a1令则上述挑选抽样中的挑选条件简化为：另一方面，在即的条件下，2/a在1,1上均匀分布，故可令2/a，则最终决定取正负值的条件简化为：于是，得到质心系各向同性散射角余弦分布的抽样方法为：如下形式的分布密度函数称为积分分布密度函数，其中f0(x，y)为任意二维分布密度函数，H(x)为任意函数。该分布密度函数的抽样方法为：7)积分抽样方法证明：对于任意x例例22. 各向同性散射方向的抽样为了确定各向同性散射方向，根据公式：对于各向同性散射，cos在1,1上均匀分布，在0,2上均匀分布。由于直接抽样需要计算三角函数和开方。定
46、义两个随机变量：可以证明，当时，随机变量x和y服从如下分布:定义区域为：则wcos 的分布可以用上述分布表示成积分分布的形式：令，则属于上述积分限内的y一定满足条件。各向同性散射方向的抽样方法为：抽样效率为：8.随机抽样的其它方法随机抽样的其它方法 1)偏倚抽样方法2)近似抽样方法3)近似-修正抽样方法4)多维分布抽样方法5)指数分布的抽样使用蒙特卡罗方法计算积分时，可考虑将积分I改写为其中f *(x)为一个与f (x)有相同定义域的新的分布密度函数。于是可以这样计算积分I：这里Xi是从f *(x)中抽取的第i个子样。1)偏移抽样方法由此可以看出，原来由f (x)抽样，现改为由另一个分布密度函
47、数f *(x)抽样，并附带一个权重纠偏因子这种方法称为偏倚抽样方法。从f (x)中抽取的Xf，满足而对于偏倚抽样，有一般情况下，Xf是具有分布f (x)总体的简单子样的个体，只代表一个。Xf*是具有分布f *(x)总体的简单子样的个体，但不代表一个，而是代表W(Xf*)个，这时Xf*是带权W(Xf*)服从分布f (x)。在实际问题中，分布密度函数的形式有时是非常复杂的，有些甚至不能用解析形式给出，只能用数据或曲线形式给出。如中子散射角余弦分布多数是以曲线形式给出的。对于这样的分布，需要用近似分布密度函数代替原来的分布密度函数，用近似分布密度函数的抽样代替原分布密度函数的抽样，这种方法称为近似抽
48、样方法。2)近似抽样方法设fa(x)f (x)，即fa(x)是f (x)的一个近似分布密度函数。对于阶梯近似，有其中，x0，x1，xn为任意分点。在此情况下，近似抽样方法为：或a)阶梯近似对于梯形近似，有其中，c为归一因子，fif (xi)，x0，x1，xn为任意分点。根据对称抽样方法，梯形近似抽样方法为：b)梯形近似除了上述这种近似外，近似抽样方法还包括对直接抽样方法中分布函数反函数的近似处理，以及用具有近似分布的随机变量代替原分布的随机变量。例例23. 正态分布的近似抽样我们知道，随机数的期望值为1/2，方差为1/12，则随机变量渐近正态分布，因此，当n足够大时便可用Xn作为正态分布的近似
49、抽样。特别是n12时，有对于任意分布密度函数f (x)，设fa(x)是f (x)的一个近似分布密度函数，它的特点是抽样简单，运算量小。令则分布密度函数f(x)可以表示为乘加分布形式：其中H1(x)为非负函数，f1(x)为一分布密度函数。对f(x)而言，fa(x)是它的近似分布密度函数，而H1(x)f1(x)正好是这种近似的修正。3)近似-修正抽样方法近似-修正抽样方法如下：抽样效率由上述近似-修正抽样方法可以看出，如果近似分布密度函数fa(x)选得好，m接近1，这时有很大可能直接从fa(x)中抽取Xfa，而只有很少的情况需要计算与f(x)有关的函数H1(Xf1)。在乘抽样方法中，每一次都要计算
50、H(Xfa)f(Xfa)fa(Xfa)。因此，当f(x)比较复杂时，近似-修正抽样方法有很大好处。例例24. 裂变中子谱分布的近似-修正抽样裂变中子谱分布的一般形式为：其中A，B，C，Emin，Emax均为与元素有关的量。对于铀-235，A=0.965，B=2.29，C=0.453，Emin=0，Emax=。若采用乘减抽样方法，其抽样效率约为0.5。令相应的则从fa(x)的抽样为从f1(x)的抽样为参数的确定，使1A0，且使H1(E)的上界M1最小。裂变中子谱的近似修正抽样方法为对于铀-235，m0.8746，M0.2678，0.5543，抽样效率E0.9333。而且近似修正抽样方法有0.87
51、46的概率直接用近似分布抽样，只计算一次对数。因此，较之乘减抽样方法大大节省了计算时间，提高了抽样效率。为方便起见，这里仅讨论二维分布的情况，对于更高维数的分布，可用类似的方法处理。对于任意二维分布密度函数，总可以用其边缘分布密度函数和条件分布密度函数的乘积表示:其中fl(x)，f2(y|x)分别为分布f(x,y)的边缘分布密度函数和条件分布密度函数，即4)多维分布抽样方法二维分布密度函数的抽样方法是:首先由fl(x)中抽取Xf1，再由f2(y|Xf1)中抽样确定Yf2。对于多维分布密度函数，也可直接采用类似于一维分布密度函数的抽样方法。例如，对如下形式的二维分布密度函数：其中H(x,y)为非
52、负函数，f1(x,y)为任意二维分布密度函数。设M为H(x,y)的上界，则有二维分布的乘抽样方法如下：例例25. 下面二维分布密度函数的抽样将f(x,y)写为其中用直接抽样方法分别从fl(x)和f2(y|Xf1)中抽样，得到前面已经介绍了，指数分布的直接抽样为：这不仅需要计算对数，而且由于要使用伪随机数，受精度的限制，该抽样值在小概率处即数值较大处呈现明显得离散性。下面介绍两种抽样方法可以避免这些问题。5)指数分布的抽样所用随机数的平均个数Ne2/(e1)4.3方法一NY方法二NY作作 业业 1)光子散射后能量分布的抽样把光子散射能量分布改写成如下形式进行抽样：在1,1+2上定义如下函数:第四
53、章第四章 蒙特卡罗方法解粒子输运问题蒙特卡罗方法解粒子输运问题1.屏蔽问题模型屏蔽问题模型2.直接模拟方法直接模拟方法3.简单加权法简单加权法4.统计估计统计估计法法5.指数变换法指数变换法6.蒙特卡罗方法的效率蒙特卡罗方法的效率作作 业业第四章第四章 蒙特卡罗方法解辐射屏蔽问题蒙特卡罗方法解辐射屏蔽问题辐射（光子和中子）屏蔽问题是蒙特卡罗方法最早广泛应用的领域之一。本章主要从物理直观出发，说明蒙特卡罗方法解决这类粒子输运问题的基本方法和技巧。而这些方法和技巧对于诸如辐射传播、多次散射和通量计算等一般粒子输运问题都是适用的。1.屏蔽问题模型屏蔽问题模型在反应堆工程和辐射的测量与应用中，常常要用
54、一些吸收材料做成屏蔽物挡住光子或中子。我们所关心的是经过屏蔽后射线的强度及其能量分布，这就是屏蔽问题。当屏蔽物的形状复杂，散射各向异性，材料介质不均匀,核反应截面与能量、位置有关时，难以用数值方法求解，用蒙特卡罗方法能够得到满意的结果。粒子的输运问题带有明显的随机性质，粒子的输运过程是一个随机过程。粒子的运动规律是根据大量粒子的运动状况总结出来的，是一种统计规律。蒙特卡罗模拟，实际上就是模拟相当数量的粒子在介质中运动的状况，使粒子运动的统计规律得以重现。不过，这种模拟不是用实验方法，而是利用数值方法和技巧，即利用随机数来实现的。为方便起见，选用平板屏蔽模型，在厚度为a，长、宽无限的平板左侧放置
55、一个强度已知，具有已知能量、方向分布的辐射源S。求粒子穿透屏蔽概率（穿透率）及其能量、方向分布。穿透率就是由源发出的平均一个粒子穿透屏蔽的数目。同时，假定粒子在两次碰撞之间按直线运动,且粒子之间的相互作用可以忽略。2.直接模拟方法直接模拟方法直接模拟方法就是直接从物理问题出发，模拟粒子的真实物理过程。1)状态参数与状态序列2)模拟运动过程3)记录结果粒子在介质中的运动的状态，可用一组参数来描述，称之为状态参数状态参数。它通常包括：粒子的空间位置r， 能量E和运动方向，以S(r,E,)表示。有时还需要其他的参数，如粒子的时间t和附带的权重W，这时状态参数状态参数为S(r,E,t,W)。状态参数状
56、态参数通常要根据所求问题的类型和所用的方法来确定。对于无限平板几何，取S(z ,E,cos)其中z为粒子的位置坐标，为粒子的运动方向与Z轴的夹角。对于球对称几何,取S(r ,E,cos)其中r表示粒子所在位置到球心的距离，为粒子的运动方向与其所在位置的径向夹角。1)状态参数与状态序列粒子第m次碰撞后的状态参数为或它表示一个由源发出的粒子，在介质中经过m次碰撞后的状态，其中rm：粒子在第m次碰撞点的位置Em：粒子第m次碰撞后的能量m：粒子第m次碰撞后的运动方向tm：粒子到第m次碰撞时所经历的时间Wm：粒子第m次碰撞后的权重有时，也可选为粒子进入第m次碰撞时的状态参数。一个由源发出的粒子在介质中运
57、动，经过若干次碰撞后，直到其运动历史结束（如逃出系统或被吸收等）。假定粒子在两次碰撞之间按直线运动，其运动方向与能量均不改变，则粒子在介质中的运动过程可用以下碰撞点的状态序列状态序列描述：S0，S1，SM-1，SM或者更详细些,用来描述。这里S0为粒子由源出发的状态，称为初态，SM为粒子的终止状态。M称为粒子运动的链长。这样的序列称为粒子随机运动的历史，模拟一个粒子的运动过程，就变成确定状态序列的问题。为简单起见，这里以中子穿透均匀平板的模型来说明，这时状态参数状态参数取S(z ,E,cos)。模拟的步骤如下：(1)确定初始状态S0：确定粒子的初始状态，实际上就是要从中子源的空间位置、能量和方
58、向分布中抽样。设源分布为则分别从各自的分布中抽样确定初始状态。对于平板情况，抽样得到z00。2)模拟运动过程(2)确定下一个碰撞点：已知状态Sm-1，要确定状态Sm，首先要确定下一个碰撞点的位置 zm。在相邻两次碰撞之间，中子的输运长度 l 服从如下分布：对于平板模型，l 服从分布：其中，t 为介质的中子宏观总截面，积分 称为粒子输运的自由程数，系统的大小通常就是用系统的自由程数表示的。 显然，粒子输运的自由程数服从指数分布，因此从 f(l)中抽样确定l，就是要从积分方程中解出l。对于单一介质则下一个碰撞点的位置如果 zma，则中子穿透屏蔽，若 zm0, 则中子被反射出屏蔽。这两种情况，均视为
59、中子历史终止。(3)确定被碰撞的原子核：通常介质由几种原子核组成，中子与核碰撞时，要确定与哪一种核碰撞。设介质由A、B、C 三种原子核组成，其核密度分别为NA、NB、NC，则介质的宏观总截面为：其中 分别为核A、B、C 的宏观总截面。其定义如下： 分别表示()核的宏观总截面、核密度和微观总截面。由于中子截面表示中子与核碰撞可能性的大小，因此，很自然地，中子与A、B、C 核发生碰撞的几率分别为：利用离散型随机变量的抽样方法，确定碰撞核种类：(4)确定碰撞类型：确定了碰撞的核(比如B核)后，就要进一步确定碰撞类型。中子与核的反应类型有弹性散射、非弹性散射、(n,2n)反应，裂变和俘获等，它们的微观
60、截面分别为则有各种反应发生的几率分别为利用离散型随机变量的抽样方法，确定反应类型。在屏蔽问题中，中子与核反应常只有弹性散射和吸收两种类型，吸收截面为：这时，总截面为：发生弹性散射的几率为：若 ，则为弹性散射；否则为吸收，发生吸收反应意味着中子的历史终止。(5)确定碰撞后的能量与运动方向：如果中子被碰撞核吸收，则其输运历史结束。如果发生弹性散射，需要确定散射后中子的能量和运动方向。中子能量 Em 为：A是碰撞核的质量与中子质量之比，一般就取元素的原子量；C为质心系中中子散射前后方向间的夹角，即偏转角。 可从质心系中弹性散射角分布fC(C)中抽样产生。实验室系散射角L的余弦L为：如果给出实验室系散
61、射角余弦分布 fL(L)，可直接从 fL(L)中抽取L，此时能量Em与L的关系式为：确定了实验室系散射角L后，再使用球面三角公式确定cosm ：其中为在0,2上均匀分布的方位角。至此，由Sm-1完全可以确定Sm。 因此，当中子由源出发后，即S0确定后，重复步骤 (2)(5)，直到中子游动历史终止。于是得到了一个中子的随机游动历史 S0，S1，SM-1，SM，即也就是模拟了一个由源发出的中子的运动过程。以上模拟过程可分为两大步：第一步确定粒子的初始状态S0，第二步由状态Sm-1来确定状态Sm。这第二步又分为两个过程：第一个过程是确定碰撞点位置zm ，称为输运过程；第二个过程是确定碰撞后粒子的能量
62、及运动方向，称为碰撞过程。对于中子而言，碰撞过程是先确定散射角，进而确定能量和运动方向；而对于光子，碰撞过程是先确定能量，再确定散射角以及运动方向。重复这两个过程，直至粒子的历史终止。这种模拟过程，是解任何类型的粒子输运问题所共有的，它是蒙特卡罗方法解题的基本手段。在获得中子的随机游动历史后，我们要对所要计算的物理量进行估计。对于屏蔽问题，我们要计算中子的穿透率。考察每个中子的随机游动历史，它可能穿透屏蔽（zMa），可能被屏蔽发射回来（zM0），或者被吸收。设第n个中子对穿透的贡献为n，则如果我们共跟踪了N个中子，则穿透屏蔽的中子数为：3)记录结果则穿透屏蔽概率的近似值为：它是穿透率的一个无偏
63、估计。我们称这种直观地模拟过程和估计方法为直接模拟方法。在置信水平 10.95时， 的误差为：其中 为n的均方差，由于n是一个服从二项分布的随机变量，所以或为得到中子穿透屏蔽的能量、角分布，将能量、角度范围分成若干个间隔：其中Emax，Emin分别表示能量的上、下限，对于穿透屏蔽的中子按其能量、方向分间隔记录。设一穿透屏蔽的中子能量为EM，其运动方向与Z轴夹角为M，若能量EM属于第 i 个能量间隔Ei，角度M属于第 j 个角度间隔j，则分别在第 i 个能量计数器及第 j 个角度计数器中加 1。跟踪 N 个中子后，则分别为穿透中子的能量分布和角分布。其中N1,i 和 N2,i 分别为第 i 个能
64、量和第 j 个角度间隔的穿透中 子数。归一后分别为：3.简单加权法简单加权法 从模拟物理过程来说，直接模拟法是最简单、也是最基本的方法。但是，在直接模拟法中，不管中子在屏蔽中经过多少次碰撞，只要在介质中被吸收，对穿透的贡献就为零；因此在所跟踪的粒子中绝大部分都对穿透没有贡献。而在许多屏蔽问题中，穿透率的数量级在10-6到10-8。进一步，如果我们要求穿透率达相对误差小于1，即那么，N要大到惊人的数量级1010到1012。显然，这时用直接模拟法计算不是很有效。屏蔽物一般是由吸收强的介质组成，因此在每次碰撞时，粒子很有可能被吸收而停止跟踪。现在改变模拟方法，在判断碰撞类型时，可以认为粒子的部分是弹
65、性散射，而其余部分被吸收，即人为地把中子分成两部分，一部分弹性散射，一部分吸收。弹性散射这部分继续跟踪；吸收部分则停止跟踪。也就是说，我们利用中子权重的变化来反应继续弹性散射的部分。这就是简单加权法的基本思想。1)简单加权法显然，在加权法中中子的权重W已成为中子状态参数的组成部分。这时，中子历史成为：对源中子，取W0=1。经过碰撞中子权重的变化为：因子称为尚存因子。这时，第n个中子对穿透的贡献为：如果我们共跟踪了N个中子，则穿透率P的无偏估计为：类似地，可以得到穿透中子的能量分布和角分布。只不过在对各计数器进行的加 1 操作改为加WM。简单加权法的方差估计为：与直接模拟法相比，有注意到n1，有
66、这表明简单加权法的方差小于直接模拟法的方差。这是因为加权法比直接模拟法减少了一次随机抽样。2)简单加权法的方差加权法的思想在蒙特卡罗方法中用途很广泛。例如，对于具有中子增殖反应，如裂变，(n,2n)，(n,3n)反应的中子输运问题，一个中子与核发生碰撞后，根据反应的类型会产生不同数量的次级中子，每个次级中子又会产生新的次级中子，这样链锁反应下去，使得用直接模拟法模拟每一个中子是非常困难的。这种情况可以利用加权法来处理。3)权重方法的其它应用中子与核发生碰撞后，产生的次级中子平均数为：这里f为裂变次级中子数。于是，碰撞后的权重为：而决定碰撞类型的几率分别为：其中加权法的思想，还可以应用到连续分布
67、情况和偏倚抽样的问题4.统计估计法统计估计法 加权法虽然改进了直接模拟法，但它同样只关心中子是否穿透屏蔽这一信息，因此对每一个中子历史的信息利用得很不充分。统计估计法能够较多地利用中子的历史信息，因而能得到更好的结果。一个中子，可能在介质内不发生碰撞而直接穿透屏蔽，也可能在介质内发生一次碰撞后再穿透屏蔽，或经过二次碰撞穿透屏蔽，等等，这些事件是互不相容的，因此穿透概率P可表示为：其中Pm是中子恰好经过m次碰撞而穿透屏蔽的概率。这表明，可以用求Pm(m=0,1,)的方法得到P。这样，中子对穿透概率的贡献就不只限于末次碰撞了。01S0S1SmmP1P0PmZa0设中子的历史为：根据该中子的历史，我
68、们可以估计出中子恰好经过m次碰撞后，穿透屏蔽的部分显然，具有初态S0(0,E0,cos0,W0)的中子，未经碰撞直接穿透的部分是：类似地，在经过了第m次碰撞后的中子具有状态Sm(zm,Em,cosm,Wm)，其可能穿透的部分，正好是一个中子恰好经过m次碰撞穿透的部分：这里的这种估计技巧，由于是对每次碰撞后的状态，求其后未经碰撞直接穿透的贡献，因此该方法也称为最后自由飞行估计。于是得到该中子对穿透的贡献：如果我们共跟踪了N个中子，则穿透率P的估计为：其方差估计为：在直接模拟方法中，相对误差为其中为与置信水平1相应的量。如果构造一个新的概率模型，使得该模型的穿透率P*与原模型的穿透率P之间存在关系
69、：使用直接模拟方法,相对误差为5.指数变换法如果令*，即这意味着，达到同样的相对误差，跟踪粒子的数目缩小K倍，从而减少K倍的计算量。指数变换法就是构造一个新的概率模型的一个有效方法。构造如下伪过程：宏观总截面为散射截面仍为el(E)。其中Emin、Emax分别为能量的下限和上限，为粒子的运动方向与Z轴的夹角。可以证明这个伪过程的穿透概率P*与原过程的穿透概率P之间有如下关系:显然，。因伪过程与原过程的结果相差e指数，所以该方法称为指数变换法。分析一下伪过程的定义,可以明显看出P*增大的原因。当cos=1时，粒子运动方向与Z轴方向一致，其截面最小，粒子沿Z轴方向输运的距离较远；而当cos=1时，
70、粒子运动方向背向Z轴方向，这时其截面最大，粒子向后输运的距离较短。因此，截面变换的结果是加强了粒子向前运动的能力，因而使穿透概率增大。伪过程的构造与几何形状及所考虑的问题有关。比如，对球形几何，使用指数变换法求穿透概率时,所构造的宏观总截面与平板屏蔽的情况不同，粒子的模拟方法也较复杂。6.蒙特卡罗方法的效率蒙特卡罗方法的效率衡量一种蒙特卡罗技巧的好坏，除了看其方差大小外，还要看其所需费用（计算时间）多少，即从该技巧的效率Ef（方差与费用乘积的倒数）全面考虑：其中2为方差，T为所需费用。Ef大时，所用方法的效率高；否则，效率低。在一般情况下，有些方法虽然减小了方差，却增加了费用。例如，加权法、统
71、计估计法虽然较直接模拟方法减小了方差，却使每个粒子的运动链长增加，或记录贡献的计算时间增加。因此，不能认为方差小的方法一定好，要从方法的效率全面考虑。在有些情况下，直接模拟方法仍然是一个被广泛使用的方法。作作 业业 1)光子散射后能量分布的抽样把光子散射能量分布改写成如下形式进行抽样：在1,1+2上定义如下函数:作作 业业 1)给出分布密度函数的抽样方法。第五章第五章 蒙特卡罗方法在计算机上的实现蒙特卡罗方法在计算机上的实现1.源分布抽样过程源分布抽样过程2.空间、能量和运动方向的随机游动过程空间、能量和运动方向的随机游动过程3.记录贡献和分析结果过程记录贡献和分析结果过程4.核截面数据的引用
72、核截面数据的引用5.蒙特卡罗程序结构蒙特卡罗程序结构作作 业业第五章第五章 蒙特卡罗方法在计算机上的实现蒙特卡罗方法在计算机上的实现蒙特卡罗方法是随着计算机的出现和发展而逐步发展起来的。在计算机上能够产生符合要求的随机数，实现对已知分布的抽样，奠定了蒙特卡罗方法在计算机上得以实现的基础。在计算机上使用蒙特卡罗方法解粒子输运问题大致包括三个过程：源分布抽样过程，空间、能量和运动方向的随机游动过程以及记录、分析结果过程。1.源分布抽样过程源分布抽样过程源分布抽样的目的是产生粒子的初始状态。下面我们介绍一些常见的特定类型的源分布抽样方法。1)源粒子的位置常见分布的随机抽样(1)圆内均匀分布设圆半径为
73、R0，粒子在圆内均匀分布时，从发射点到中心的距离r的分布密度函数为：r的抽样方法为：(2)圆环内均匀分布设圆环的内半径为R0，外半径为R1，则粒子在该圆环内均匀分布时，从发射点到中心的距离r的分布密度函数为：r的抽样方法为：(3)球内均匀分布设球的半径为R，粒子在球内均匀分布时，从发射点到中心的距离r的分布密度函数为：r的抽样方法为：在直角坐标系下，抽样方法为：(4)球壳内均匀分布设球壳的内半径为R0，外半径为R1，在均匀分布时，从发射点到中心的距离r的分布密度函数为：r的抽样方法为：在直角坐标系下，球壳内点的坐标为：其中，r由前面的抽样方法确定，、服从各向同性分布，其抽样方法为：(5)圆柱内
74、均匀分布圆柱内均匀分布是指粒子发射点均匀地分布在底半径为R，高为2H的圆柱内。若固定圆柱的中心为原点，圆柱的轴向为z轴，则分布密度函数为：抽样方法为：(6)点源分布点源分布是指粒子由一固定点发射，其分布密度函数为：其中，为狄拉克函数，源粒子的抽样方法为：在球坐标系中，粒子发射点到球心的距离r的分布密度函数为：其中，为点源到球心的距离。源粒子的位置抽样为：(7)球外平行束源分布球外平行束源分布是指粒子平行入射到半径为R的球面上，或球外点源距离球很远，可以近似地看作平行束源。设r为粒子发射点到球心的距离,其分布密度函数为：r的抽样方法为：在直角坐标系中，抽样方法为：2)源粒子的能量常见分布的随机抽
75、样(1)单能源分布单能源分布是指粒子的发射能量为一固定值E0，其分布密度函数为：源粒子的能量为：(2)裂变中子谱分布裂变中子谱分布的一般形式为：其中A，B，C，Emin，Emax均为与元素有关的量。对于铀-235，A=0.965，B=2.29，C=0.453，Emin=0，Emax=。采用近似修正抽样，抽样方法为：其中，m0.8746，M10.2678，0.5543。此外，裂变谱分布也有以数值曲线形式给出的，此时，用数值曲线抽样方法抽取E。(3)麦克斯韦(Maxwell)谱分布麦克斯韦谱分布的一般形式为：该分布的抽样方法为3)源粒子运动方向常见分布的随机抽样(1)各向同性分布各向同性分布密度函
76、数为：其中，cos，为运动方向与z轴的夹角，为方位角。在直角坐标系下，各方向余弦u，v，w为：其抽样方法为：(2)半面各向同性分布不妨设在x0的半面方向上各向同性发射粒子，则在前述各向同性分布的抽样方法中，用2代替2就能得到所需分布的抽样。对于其它方向的情况，可用类似的方法处理。(3)球外平行束源分布令cos，为粒子运动方向的径向夹角，则分布密度函数为：的抽样方法为：(4)球外各向同性点源分布设球外点源S到球心的距离为D0。点源S到球的最大张角为*，则球外各向同性点源分布的抽样方法是：先抽样确定，再转换成。在直角坐标系下，取OS为z轴，抽样方法为：4)次级粒子的源分布在有关次级粒子（如裂变中子
77、，中子生成光子，光子生成中子）的输运过程中，次级粒子源分布的抽样方法，主要可分为以下两种：(1)直接生成法可将生成的次级粒子的位置、能量、方向、权重等参数直接作为源分布的抽样结果。也就是直接对生成的次级粒子进行跟踪。这种方法比较简单、直观。(2)离散分布法将生成的次级粒子的权重，按空间位置、能量、方向分别记录，得到次级粒子的空间、能量、运动方向的离散的近似分布。再根据该分布，利用各种抽样技巧，得到源分布的抽样，对抽样的源粒子进行跟踪、记录。当一个问题需要用两个以上的蒙卡程序处理时，可采用这种方法。2.空间、能量和运动方向的随机游动过程空间、能量和运动方向的随机游动过程粒子由状态Sm到状态Sm+
78、1时，需要确定粒子的空间位置rm+1，能量Em+1和运动方向m+1。1)碰撞点位置的计算公式设rm为粒子第m次碰撞点的位置，m为碰撞后的运动方向，则粒子第m+1次碰撞点的位置rm+1为：即其中为的方向余弦，L为两次碰撞点间的距离。L的分布密度函数为：由f(L)抽样确定L的方法通常有三种：(1)直接抽样方法确定L的直接抽样方法是：首先由自由程分布中抽取再由下列关系式解出L。对于均匀介质，有对于多层介质，如果则其中，为粒子由rm出发，沿m方向在顺序经过的第i个介质区域内走过的距离，为第i个介质区域的宏观总截面(i=1，2，Imax)。当时，意味着粒子穿出系统。(2)最大截面法对于多层介质，或其他介
79、质密度与位置有关的问题，在求(i=1，2，Imax)时，如果系统形状复杂，计算是非常烦杂的。在这种情况下，使用最大截面法更方便。最大截面抽样方法为：其中(3)限制抽样法当介质区域很小时，如使用直接抽样法抽取输运长度，粒子很容易穿出介质，此时使用限制抽样法确定自由程个数较好，的分布密度函数为：其中Dm为粒子由rm出发，沿m方向到达区域边界的自由程个数。的抽样方法是：然后用直接抽样法中根据计算L 的方法计算输运长度L。此时，粒子的权重需乘以纠偏因子。2)碰撞后能量Em+1的随机抽样粒子在介质中发生碰撞后，首先要确定与哪种原子核发生何种反应。粒子发生碰撞后（吸收除外）的能量Em+1一般只与其碰撞前后
80、运动方向的夹角（散射角）有关。粒子碰撞后常见的能量分布有下面几种情况。(1)裂变中子谱中子引起原子核裂变反应时，裂变中子的能量服从裂变谱分布。其抽样方法可参考以前的介绍。(2)中子弹性散射后能量的确定中子弹性散射后，能量与质心系散射角C的关系是：能量与实验室系散射角L的关系是：其中，A为碰撞核的质量，。或确定后，即可求出Em+1。(3)中子非弹性散射后能量的确定中子非弹性散射后，能量与质心系散射角C的关系是：其中，为第K个能级的阈能，为第K个能级的激发态能量。如果确定了实验室系散射角L，则根据下式确定后，再计算出Em+1。(4)光子康普顿（Compton）散射后能量的确定光子发生康普顿散射后，
81、其能量分布密度函数为：其中，K()为归一因子。，和分别为光子散射前后的能量，以m0c2为单位，m0为电子静止质量，c为光速。光子康普顿散射能量分布的抽样方法为：x的抽样确定后，散射后的能量为：3)碰撞后散射角的随机抽样粒子碰撞后运动方向m+1的确定，一般与散射角有关。由已知分布抽样确定散射角后，再确定m+1。常见的散射角分布有如下几种：(1)质心系各向同性分布散射角在质心系服从各向同性分布时，其抽样方法为。质心系散射角C抽样确定后，需转换成实验室系散射角L。在中子弹性散射情况下，转换公式为：其中A为碰撞核质量，。在中子非弹性散射情况下，转换公式为：其中，为第K个能级的阈能。(2)中子弹性散射勒
82、让德(Legendre)多项式分布中子弹性散射角分布常以勒让德多项式的展开形式给定。散射角余弦xcos的分布密度函数为：其中Pl(x)为l阶勒让德多项式。该分布即为n阶勒让德近似展开。勒让德多项式由以下递推公式确定：考虑新的分布：当选取x0，x1，xn为Pn+1(x)0的根，且时，fa(x)依照勒让德多项式展开的前n项与f(x)的展开形式相同。因此，可以用fa(x)作为f(x)的近似分布。在实际问题中，由于勒让德多项式展开项数不够，可能出现某个为负值的现象。此时可以采用如下近似分布：其中：对于该近似分布，可用加抽样方法进行抽样：此时，由于偏倚抽样而引起的纠偏因子为wK，也就是说，粒子的权重要乘
83、上wK。(3)光子康普顿散射角分布光子的康普顿散射角与其散射前后的能量有关,它的分布密度函数为：抽样方法为：4)碰撞后运动方向m+1的确定实验室系散射角L确定后，依据不同的坐标系的表现形式，有不同的确定方法。(1)确定方向余弦um+1，vm+1，wm+1其中，方位角 在0,2上均匀分布。当时，不能使用上述公式，可用下面的简单公式：(2)确定m+1的球坐标(m+1，m+1)设m的球坐标分别为(m,m)，其中，为粒子运动方向与z轴的夹角，为粒子运动方向在x y平面上投影的方位角。则m+1的球坐标(m+1,m+1)分别由下式确定：5)球形几何的随机游动公式一般几何的随机游动公式可以应用到球形几何，而
84、对球对称问题，使用特殊形式更为方便。(1)下次碰撞点的径向位置rm+1的确定两次碰撞点间的距离L确定之后，下次碰撞点的径向位置rm+1的计算公式为：设系统的外半径为R，如rm+1R，则粒子逃出系统。(2)粒子碰撞后瞬时运动方向的确定在球对称系统中，粒子运动方向用其与径向夹角余弦来描述。使用球面三角公式，粒子碰撞后瞬时运动方向与径向夹角余弦cosm+1的计算公式为：其中，为在0,2上均匀分布的方位角，为在rm+1点进入碰撞前瞬时运动方向与rm+1径向之间的夹角。6)点到给定边界面的距离在抽样确定输运距离、判断粒子是否穿透系统时，常遇到求由rm出发，沿m方向到达某个区域表面的距离问题。在记录对结果
85、的贡献时，也常使用类似的量。区域表面通常是平面或二次曲面。求到达区域表面的距离问题，实际上是求直线（或半射线）与平面或二次曲面的交点问题。这是蒙特卡罗方法解粒子输运的各种实际问题时,所遇到的基本几何问题。(1)点到平面的距离点沿方向的直线方程为：该直线到达方程为的平面的距离为：当与平面平行时，即直线与平面无交点。如果l为负值，直线与平面也无交点。这时，粒子的运动方向是背离平面的。(2)点到球面的距离在三维直角坐标系中，设球心为rc(xc,yc,zc)，球半径为R，则球面方程为：将直线方程代入该球面方程，得到点r0沿0方向到达球面的距离l：其中当R0R时，即r0点在球内，2，l只有一个正根：当R
86、0R时，即r0点在球外，分以下三种情况：a)若0，l无正实根，直线与球面无交点。b)若0，0，l无实根，直线与球面无交点。c)若0，0，l有两个正实根，直线与球面有两个交点。在球坐标系中，不失一般性，设球心为rc0，则球面方程为rR。当r0R时，即r0点在球内，有一个交点：其中0为0与r0的径向夹角。当r0R时，即r0点在球外，令当cos00时，直线与球面无交点。当cos00时，若dR，则直线与球面无交点。若dR，则有两个交点：(3)点到圆柱面的距离设圆柱面的方程为：其中(xc,yc,0)为圆柱的中心，R为圆柱底半径。点r0沿0方向到达圆柱面的距离l为：其中当R0R时，r0点在圆柱内，如果，则
87、l有一个正根：如果，即0平行于圆柱的对称轴，直线与圆柱面无交点。当R0R时，r0点在圆柱外，分以下三种情况：a)若0，l无正实根，直线与圆柱面无交点。b)若0，0，l无实根，直线与圆柱面无交点。c)若0，0且，l有两个正实根，直线与圆柱面有两个交点。d)在的情况下，直线与圆柱面不相交。(4)点到圆锥面的距离设圆锥顶点在原点，以z轴为对称轴，则圆锥面的方程为：点r0沿0方向到达圆锥面的距离l为：其中如果0与锥面某一母线平行，即，则(5)空腔处理在粒子输运问题中，所考虑的系统常有空腔存在，如中空的球壳,平板间的空隙等。粒子输运时，可有两种处理空腔的方法：a)将空腔作为宏观总截面t0的区域,按通常的
88、方法输运。b)设分别为由rm出发，沿m方向到空腔区域的近端和远端的交点，当粒子超过时，以为新的起点，重新开始输运。显然，这两种方法在统计上是等价的。7)等效的边界条件(1)全反射边界在反应堆活性区中，元件盒常常按正方形或六角形排列。假定元件盒足够多，每个盒结构相同，那么活性区中每个盒所占的栅胞的物理情况，可以代表整个活性区中的状况。进一步假定，元件盒是圆对称的，那么每个栅胞中情况，可以用更小的单位（栅元）来反映。比如对六角形栅胞可取其1/12的OAB来做代表；正方形栅胞可用其1/8的OAB来做代表。这样一来问题就大大简化了。现在的问题是怎样计算直角三角形栅元的物理量（如通量）。用蒙特卡罗方法如
89、何模拟中子在栅元内的运动，反映出整个活性区对它的影响。我们可把OA、OB、AB 作为全反射边界来处理。所谓全反射边界，就是当中子打到该边界上时，按镜面反射的方式，从边界上全部反射回来，中子的能量与权重均不改变。在这种边界上的反射条件，称之为全反射条件，就是通常的镜面反射条件。在全反射边界条件下，一条通过活性区若干个区域的中子径迹，可以用栅元OAB中的一条折线轨道来反映出来。反过来，在直角三角形栅元OAB中任一条反射成的折线轨道，都代表了中子在活性区内一条直线轨道的作用。由于系统的对称性，在活性区内，凡是与栅元内位置相当的地方，都有相同的物理情况，因此栅元内各处的情况，当然代表了整个活性区的情况
90、。(2)一般曲面全反射条件对于一般曲面的全反射，设入射方向为，入射点的内法线方向为n，则反射方向为：其中设则(3)平面全反射条件设三角形栅元的横截面OAB在X-Y平面上，OAB。则边界OA、OB、AB 上的反射都是平面全反射。在任一与X-Y平面垂直且与X轴成角的平面上，全反射条件为：由此就可得到OA、OB和AB边上的全反射条件，对于OB边，=；对于OA边，=0；对于AB边，=/2。(4)反射层边界条件对于具有大反射层的系统，如存放，运输和生产裂变物质的仓库、车厢和车间等，当中子从里面打到四周墙上或反射层时，还要继续对它进行跟踪。这种跟踪常常要花费很大的计算量，并且在结果中引起的方差也比较大。如
91、果在计算这种系统的不同方案中，反射层条件不变，那么这种大量重复的计算是很不经济的。中子射入反射层后，一部分被介质吸收，只有一部分返回，由于中子的散射慢化，损失一部分能量，因此反射回来的中子有一个能量方向分布。显然，对这种反射层，不能应用全反射条件。不过，我们仍然可以把它当做边界，在边界上按反射层的物理作用来处理。比如，如果反射层是一种平板几何，我们可以用数值方法或蒙特卡罗方法，预先算好在各种不同入射能量E下的反照率(E)，反射中子的能量分布RE(EE)。于是代替在反射层中眼踪中子，我们可在反射层边界上作如下处理：一旦中子打入反射层，立即返回，反射后权重为其中，E为射入反射层中子的能量，W为中子
92、的权重。反射后的能量E由反射能谱RE(EE)中抽样产生。反射后的方向由半平面各向同性分布或余弦分布中抽样。反射后的中子位置为入射时的位置。计算表明，对于大尺寸的反射层来说，这样的近似，引起的结果上的误差是可以忽略的，却能带来计算量的大量节省。3.记录贡献与分析结果过程记录贡献与分析结果过程在粒子输运问题中，除了得到某些量的积分结果外，还需要得到这些量的方差、协方差、以及这些量的空间、能量、方向和时间的分布。这些量可以利用分类记录手续同时得到。1)记录与结果为了得到所求量的估计值，在粒子输运过程中需进行记录，即求每个粒子对所求量的贡献。设模拟了N个粒子，所求量的估计值为：其中gi为第i个粒子的总
93、贡献。记录的贡献由所求量决定。对于同一个所求量，又随所用的蒙特卡罗技巧的不同而不同。例如，所求量是粒子穿透屏蔽概率，使用直接模拟法时，如粒子穿透屏蔽，在叠加记录单元加“1”(初始值为零)，否则没有贡献。使用加权法时，如粒子穿透屏蔽，在叠加记录单元加粒子的权重，否则没有贡献。使用统计估计法时，粒子每发生一次碰撞(包括零次碰撞)，都要记录贡献，等等。2)方差和协方差的估计估计量g和g的方差和协方差为：它们可以用下式估计：因此，要得到和的估计，只要对每一个历史记录结果的和进行记录，并加以累加即可。方差估计值确定后，可得到误差其中为置信限，它随置信水平而定。在通常情况下，取。3)位置、能量、方向、时间
94、分布在前面已经提到，用蒙特卡罗方法求某种量的空间、能量、方向和时间分布，实质上是得到这种分布的阶梯函数近似的估计值。而求这种估计值是很方便的，只要将跟踪过程中所得到的感兴趣量，按其状态的空间、能量、方向、时间特征，分别记录其权重，最后将这些记录结果适当处理即可。事先，将问题的空间、能量、方向（常按相对于某个方向的夹角余弦）、时间范围，各分为如下不同间隔：再用一批存贮单元A记录相应间隔上阶梯函数近似的累计值。4.核截面数据的引用核截面数据的引用用蒙特卡罗方法解粒子输运问题，需要介质所包含的各种原子核的核数据。以中子核数据为例，需要各种涉及到的核的微观总截面、弹性散射、非弹性散射、n-2n反应、裂
95、变、俘获等截面；也需要这些反应的相应能量、角度分布、次级粒子数，以及其它关心的粒子数及其能量、方向分布。从输运方程中可以知道，有了这些数据，问题就完全确定了；反映到蒙特卡罗模拟中，有了这些数据，就能决定宏观总截面，决定碰撞核的具体形式，就能实现抽样和跟踪。在蒙特卡罗计算中，引用的核数据有点截面、分段常数截面和群截面三种形式。1)点截面形式在跟踪粒子时，对粒子的每一种能量，先从截面库中取出需要的核数据，再用插值（或其它方式），求出相应能量的各种截面数据。这种方法是直接的，也是比较精确的。不少通用程序就是这样做的。这样做的条件是要有完备的核截面数据库，计算机有大而快的存贮能力。2)分段常数形式将问
96、题的能量范围（Emin，Emax）分成许多间隔，截面数据在每个间隔上看作与能量无关，即截面取分段常数形式。这种引用截面数据的好处是，数据量相对地少。问题是它要根据问题的物理特点来划分能量间隔。而为了保证误差较小，所取的间隔数一般是比较多的。3)群截面形式解中子输运问题，常常采用多群近似方法。在多群近似下，中子能量E用群指标i代替。为了实现蒙特卡罗跟踪，只需要以下群截面：显然，这种跟踪过程数据量小，程序简单。5.蒙特卡罗程序结构蒙特卡罗程序结构在电子计算机上，蒙特卡罗方法解粒子输运问题的程序，一般都可分为：源抽样，空间输运过程，碰撞过程，记录过程和结果的处理与输出等部分。开始数据预处理,各记录单
97、元清零取一个粒子历史源分布抽样输运过程碰撞过程历史终止否？统计处理做完给定历史数否？结果的处理与输出终止记录过程记录过程记录过程记录过程至于粒子历史终止条件，根据问题的几何条件、物理假定，处理方法，可归纳为以下几种：1)粒子从系统逃脱；2)粒子经碰撞被吸收；3)经俄国轮盘赌后，历史被终止；4)粒子能量低于给定能量（阈能）;5)粒子位置越过某一界面；6)粒子飞行时间超过给定时间；7)粒子权重小于某个小量。作作 业业 1)第六章第六章 蒙特卡罗方法在通量计算中的应用蒙特卡罗方法在通量计算中的应用1.通量的定义通量的定义2.通量的能谱和角分布通量的能谱和角分布3.计算体通量的模拟方法计算体通量的模拟
98、方法4.计算面通量的模拟方法计算面通量的模拟方法5.计算点通量的模拟方法计算点通量的模拟方法6.与与通量有关的物理量的计算通量有关的物理量的计算作作 业业第六章第六章 蒙特卡罗方法在通量计算中的应用蒙特卡罗方法在通量计算中的应用通量计算在粒子输运问题中占有非常重要的地位。很多问题，如碰撞率、反应率以及系统逃脱几率等都可以通过通量来计算。通量计算问题，包括点通量、面通量和体通量的计算问题。相对来说，点通量的计算要困难一些。1.通量的定义通量的定义设分别表示粒子的位置、能量和运动方向。则通量的定义为：在r点的体积元dV内，能量E和运动方向属于dE d的粒子平均径迹长度。1)点通量的定义给定点r0的
99、点通量为：点通量的含义为：在r0点的体积元dV内，粒子的平均径迹长度。2)面通量的定义给定曲面A0上的面通量为：面通量的含义为：沿曲面A0的法线方向增加厚度ds所组成的体积元的体积元A0ds中，粒子的平均径迹长度。3)体通量的定义给定体V0内的体通量为：体通量的含义为：在体V0内，粒子的平均径迹长度。4)粒子各次散射对通量的贡献通量可用粒子各次散射对通量的贡献和表示：其中为粒子n次散射后对通量的贡献，其含义为：粒子在第n次散射到第n1次散射之间，在r点的体积元dV内，能量E和运动方向属于dE d的粒子平均径迹长度。2.通量的能谱与角分布通量的能谱与角分布用蒙特卡罗方法计算通量的能谱与角分布，所
100、采用的手段与计算其它物理量一样，即把能量和方向分成若干个区间，分别按粒子状态所处的区间累积记录各自的贡献。现将能量分成I区：E1，E2，EI；方向分成J区：1，2，I。则有：3.计算体通量的模拟方法计算体通量的模拟方法在实际问题中，经常遇到要计算某一区域V0的体通量。在通量的定义部分已经介绍过，通量可以表示为粒子各次散射对通量的贡献和。因此，下面要介绍的各种估计方法，只叙述各次散射后的通量计算方法。计算体通量的方法主要有以下几种。1)解析(统计)估计方法粒子n次散射（n0时为源粒子）后的通量贡献为：其中，s1和s2分别为粒子由点rn出发，沿n方向到达区 域V0的近端和远端的交点的距离。如果点r
101、n在V0内，则s10。如果粒子沿n方向与V0有多段相交，则为每段相交线段的通量贡献之和。如果粒子沿n方向与V0不相交，则。解析估计方法就是把体通量的贡献表达式直接计算出来。当系统为均匀介质时，如果只是V0为均匀介质，则如果V0由多层介质组成，则需分段计算积分。在解析估计方法中，粒子每发生一次碰撞（包括零次散射），都要记录通量的贡献值。2)径迹长度方法设粒子从第n次散射到第n1次散射之间走过的径迹长度为s ，则 n次散射的通量贡献为：径迹长度方法就是把粒子在V0内走过的径迹长度记录下来。下面证明，径迹长度估计是无偏的。3)碰撞密度方法设粒子从第n次散射到第n1次散射之间走过的径迹长度为s ，则
102、n次散射的通量贡献为：碰撞密度方法就是把粒子在V0内发生的碰撞记录下来。下面证明，碰撞密度估计是无偏的。4)均匀径迹长度方法确定一个定义在s1,s2上的概率密度函数fn(s)，从fn(s)中抽样s*，则 n次散射通量贡献的估计为：fn(s)的最简单形式是均匀分布这时5)点通量代替方法设为在V0上定义的任一概率密度函数，则体通量可表示为：体通量的估计为：其中，r*为从中抽取的一个样本值。6)几种方法的比较(1)解析估计方法：直接计算体通量的贡献表达式，因此该方法的方差小，但计算时间长，需要计算指数函数的积分。(2)径迹长度方法：记录贡献方法简单，可与输运过程同时进行，只要粒子穿过记录区域就有贡献
103、。但该方法方差大些，对于较小的系统（如自由程个数小于2），该方法较好。(3)碰撞密度方法：由于只在记录区域内发生碰撞才有贡献，因此方差较大，尤其在记录区域较小时更是如此。但该方法省时间，适用于大的记录区域。(4)均匀径迹长度方法：在记录区域为多层介质时，较解析估计方法容易实现。但在记录贡献时仍需计算指数函数，也费时间。(5)点通量代替方法：可以较好地解决小区域的体通量计算问题。尤其是记录区域与粒子的输运区域分开时，更是如此。4.计算面通量的模拟方法计算面通量的模拟方法计算面通量的方法主要有以下几种。1)解析估计方法设经过n次散射的粒子，由点rn出发，沿n方向到达曲面域A0的距离为s1，与曲面相
104、交处曲面的法线方向为n，则n次散射粒子对该曲面的通量贡献为：如果粒子沿n方向与A0有多个交点，则为每个交点处的通量贡献之和。如果粒子沿n方向与A0没有交点，则。解析估计方法就是把面通量的贡献表达式直接计算出来。粒子每发生一次碰撞（包括零次散射），都要记录通量的贡献值。2)加权（径迹长度）方法设粒子从第n次散射到第n1次散射之间走过的径迹长度为s ，则 n次散射的通量贡献为：加权方法只有在粒子穿过曲面A0时，才对该曲面有通量贡献。3)点通量代替方法设为在A0上定义的任一概率密度函数，则面通量可表示为：面通量的估计为：其中，r*为从中抽取的一个样本值。4)体通量代替方法沿曲面A0的法线方向均匀地增
105、加一个厚度s，由此构成的体积为。的体通量为：A0的面通量为：因此，如取得足够小，有如下近似：5.计算点通量的模拟方法计算点通量的模拟方法与体通量、面通量的计算相比，点通量的计算最困难。这是因为，在大量的模拟粒子中，只能有很少的粒子穿过该点所包含的一个小区域，因此无法使用通常的通量计算方法。1)指向概率方法设n次散射后粒子的状态为，进入n次碰撞的粒子的状态为，表示粒子的碰撞核，其定义为：一个粒子在点r发生碰撞后，能量由E变为E的dE内，方向由变为的d内的粒子平均数。则n次散射的粒子对点r*的通量贡献为：其中当n0时，用源分布密度函数代替碰撞核。(1)光子问题的指向概率方法光子问题的碰撞核为：其中
106、光子能量E以电子静止能量mec20.511MeV为单位；K(EE/r)为KleinNishina公式，由下式确定N(r)表示在r处单位立方体内的原子数，z (r)表示在r处元素的原子序数，r0表示电子的经典半径。其中(2)中子问题的指向概率方法中子问题的碰撞核为：其中下标A和i分别表示不同的原子核和不同的反应；和分别表示能量为E的中子与第A种原子核发生第i种反应后产生的平均次级中子数和微观截面；NA(r)表示在r处第A种原子核的核密度；表示能量为E和方向为的中子与第A种原子核发生第i种反应后的能量E和方向的分布。则有中子的通量贡献为：2)关于估计量无界问题当r*点附近不含散射物质时（如真空），
107、也就是说，粒子的输运区域与记录点分开时，指向概率方法的估计量是有界的，因此是一种比较好的计算点通量的方法。不含散射物质的区域越大，指向概率方法的优点越明显。然而，当r*点附近含有散射物质时，由于在指向概率方法的估计量中含有无界因子因此，指向概率方法的估计量一般来说是无界的。6.与通量有关的物理量的计算与通量有关的物理量的计算一些物理量可以通过通量来计算。1)系统逃脱概率状态为的粒子的系统逃脱概率为：系统逃脱概率为：2)各种反应率(1)碰撞密度(2)裂变中子密度(3)吸收率(4)反应率作作 业业 1)第七章第七章 蒙特卡罗方法在积分计算中的应用蒙特卡罗方法在积分计算中的应用1.蒙特卡罗方法求积分
108、蒙特卡罗方法求积分2.重要抽样重要抽样3.俄国轮盘赌和分裂俄国轮盘赌和分裂4.半半解析方法解析方法5.系统抽样系统抽样6.分层抽样分层抽样第七章第七章 蒙特卡罗方法在积分计算中的应用蒙特卡罗方法在积分计算中的应用计算多重积分是蒙特卡罗方法的重要应用领域之一。本章着重介绍计算定积分的蒙特卡罗方法的各种基本技巧，而这些技巧在粒子输运问题中也是适用的。1.蒙特卡罗方法求积分蒙特卡罗方法求积分蒙特卡罗方法求积分的一般规则如下：任何一个积分，都可看作某个随机变量的期望值，因此，可以用这个随机变量的平均值来近似它。设欲求积分其中，PP(x1，x2，xs)表示s维空间的点，Vs表示积分区域。取Vs上任一联合
109、概率密度函数f(P)，令则即是随机变量g(P)的数学期望，P的分布密度函数为f(P)。 现从f(P)中抽取随机向量P 的N 个样本：Pi，i1，2，N，则就是的近似估计。2.重要抽样重要抽样1)偏倚抽样和权重因子取Vs上任一联合概率密度函数f1(P)，令则有现从f1(P)中抽样N 个点：Pi，i1，2，N，则就是的又一个无偏估计。2)重要抽样和零方差技巧要使最小，就是使泛函If1极小。利用变分原理，可以得到最优的f1(P)为特别地，当 g(P)0时，有这时即g1的方差为零。实际上，这时有不管那种情况，我们称从最优分布fl(P)的抽样为重要抽样，称函数|g(P)|为重要函数。3.俄国轮盘赌和分裂
110、俄国轮盘赌和分裂1)分裂设整数n1，令则于是计算的问题，可化为计算n 个i 的和来得到，而每个gi(P)为原来的估计g(P)的1/n，这就是分裂技巧。2)俄国轮盘赌令0q1，则于是变为一个两点分布的随机变量的期望值，的特性为：这样就可以通过模拟这个概率模型来得到，这就是俄国轮盘赌。3)重要区域和不重要区域我们往往称对积分贡献大的积分区域为重要区域，或感兴趣的区域；称对积分贡献小的区域为不重要区域，或不感兴趣的区域。考虑二重积分令R是V2上x的积分区域，表为RR1+R2，其中R1是重要区域，R2是不重要区域，两者互不相交。又命Q为V2上相应于y的积分区域。则通常蒙特卡罗方法，由f (x,y)抽样
111、(x,y)的步骤是：从fl(x)中抽取xi，再由f2(y|xi)中抽样确定yi，然后用作为的一个无偏估计。现在，改变抽样方案如下：(1)当xR1时，定义一个整数n(xi)1，对一个xi，抽取(2)n(xi)个yij，j1，2，n(xi)。以平均值(3)代替上述估计式中的 g(yi,xi)。(2)当xR2时，定义一个函数q(xi)，0q(xi)1，(3)以抽样值代替上述估计式中的 g(yi,xi)。这里是随机数。显然，这种抽样估计技巧，就是对xR1时，利用分裂技巧，而对xR2时，利用俄国轮盘赌，而使估计的期望值不变。由于对重要区域多抽样，对不重要区域少观察，因此能使估计的有效性增高。4.半解析（
112、数值）方法半解析（数值）方法考虑二重积分令则x为的无偏估计。x的方差为而由f (x,y)抽样(x,y)，用g (x,y)作为的估计，其方差为5.系统抽样系统抽样我们知道，由f (x,y)抽样(x,y)的步骤是：从fl(x)中抽取xi，再由f2(y|xi)中抽样确定yi，现在改变xi 的抽样方法如下：yi的抽样方法不变。其方差为与通常蒙特卡罗方法相比，方差减少了约6.分层抽样分层抽样考虑积分在（0，1）间插入J1个点001J-1J1令则有现在，用蒙特卡罗方法计算j ，对每个j 利用 fj(x)中的nj 个样本xij ，那么有第八章第八章 蒙蒙特卡罗方法应用程序介绍特卡罗方法应用程序介绍1.蒙特卡
113、罗方法应用软件的特点蒙特卡罗方法应用软件的特点2.常用的通用蒙特卡罗程序简介常用的通用蒙特卡罗程序简介3.MCNP程序输入的描述程序输入的描述第八章第八章 蒙蒙特卡罗方法应用程序介绍特卡罗方法应用程序介绍建立完善的通用蒙特卡罗程序可以避免大量的重复性工作，并且可以在程序的基础上，开展对于蒙特卡罗方法技巧的研究以及对于计算结果的改进和修正的研究，而这些研究成果反过来又可以进一步完善蒙特卡罗程序。1.蒙特卡罗方法应用软件的特点蒙特卡罗方法应用软件的特点通用蒙特卡罗程序通常具有以下特点：1)具有灵活的几何处理能力具有灵活的几何处理能力2)参数通用化，使用方便参数通用化，使用方便3)元素和介质材料数据
114、齐全元素和介质材料数据齐全4)能量范围广，功能强，输出量灵活全面能量范围广，功能强，输出量灵活全面5)含有简单可靠又能普遍适用的抽样技巧含有简单可靠又能普遍适用的抽样技巧6)具有较强的绘图功能具有较强的绘图功能2.常用的通用蒙特卡罗程序简介常用的通用蒙特卡罗程序简介1)MORSE程序程序较早开发的通用蒙特卡罗程序，可以解决中子、光子、中子光子的联合输运问题。采用组合几何结构，使用群截面数据，程序中包括了几种重要抽样技巧，如俄国轮盘赌和分裂技巧，指数变换技巧，统计估计技巧和能量偏移抽样等。程序提供用户程序，用户可根据需要编写源分布以及记录程序。2)EGS程序程序EGS是Electron-Gamm
115、aShower的缩写，它是一个用蒙特卡罗方法模拟在任意几何中，能量从几个KeV到几个TeV的电子-光子簇射过程的通用程序包。由美国Stanford LinearAcceleratorCenter提供。EGS于1979年第一次公开发表，提供使用。EGS4是1986年发表的EGS程序的最新版本。3)MCNP程序程序MCNP是美国LosAlamos国家实验室开发的大型多功能通用蒙特卡罗程序，可以计算中子、光子和电子的联合输运问题以及临界问题，中子能量范围从10-11MeV至20MeV，光子和电子的能量范围从1KeV至1000MeV。程序采用独特的曲面组合几何结构，使用点截面数据，程序通用性较强，与其
116、它程序相比，MCNP程序中的减方差技巧是比较多而全的。3.MCNP程序输入的描述程序输入的描述 MCNP的输入包括几个文件，但主要的一个是由用户编写的INP文件，该文件包括描述问题所必须的全部输入信息。文件采用卡片结构，每行代表一张卡片，文件由一系列卡片组成，对于任一特定的问题，只需用到INP全部输入卡片的一小部分。MCNP输入文件中物理量的单位输入文件中物理量的单位长度厘米能量MeV时间10-8秒温度MeV（kT）原子密度1024个原子/厘米3质量密度克/厘米3截面10-24厘米2原子量中子质量的1.008664967倍阿伏加德罗常数 6.02310231)输入文件的基本形式输入文件的基本形
117、式(1)信息块信息块(2)信息块的卡片放在INP文件中标题卡之前。信息块给出了MCNP的一些运行信息，信息块上各部分的意思和运行行信息是一样的，当运行行信息与信息块中所指定的信息相矛盾时，则忽略信息块中相应的信息，而以运行行信息为准。(3)信息块是可选的，信息块的第一张卡片，必须在第18 列写上“MESSAGE:”，从第一张卡片的第980 列到后续卡片的第180 列都可填写运行信息。在标题卡之前用一个空行分隔符结束信息块。(2)初始运行的输入文件初始运行的输入文件信息块空行分隔符选择项标题卡仅一行，占用第180 列。作为输出标题。栅元卡空行分隔符定义构成整个系统的各个基本介质单元以及相应的物理
118、信息。曲面卡空行分隔符定义组成栅元的曲面信息。数据卡空行分隔符其它数据，包括问题类型、源描述、材料描述、计数描述，问题截断条件等。其它选择项(3)接续运行的输入文件接续运行的输入文件接续运行必须在运行行信息或信息块中给出C项选择，即Cm，表示从RUNTPE文件中读出第m次转储的内容接着运算，如果m未指定，则读最后一次转储的数据。如果不需要改变内容，则不需要接续输入文件，仅需运行RUNTPE以及在运行行加上C选择。信息块空行分隔符选择项CONTINUE写在第18 列数据卡空行分隔符只允许部分数据卡。（FQ，DD，NPS，CTME，IDUM，RDUM，PRDMP，LOST，DBCN，PRINT，K
119、CODE，MPLOT，ZA，ZB，和ZC）其它选择项(4)卡片格式卡片格式INP输入文件的每一行（称之为一张卡片）都限于使用第180 列并构成卡片映象。大部分输入卡片按行填写；然而，对数据卡允许按列填写。$符号为它所在那行数据的结束符，在$符号后面的内容作为注释，它可从$符号后面的任一列开始。标题卡只占一行，整行都可填入用户需要的信息，也可以是空行。但要注意在其它地方使用空行是作为结束符或者分隔符。输入文件中，在标题卡之后及最后的空行结束卡之前的任何地方都可插入注释卡。注释卡必须是字母“C”写在15列中的任意位置，且至少用一个空格隔开后面的注释内容。a)行输入格式行输入格式栅元卡、曲面卡和数据
120、卡的书写格式是相同的。必须从15 列开始填写这些卡片相应的名字(或编号)和粒子标识符，后面填写用空格分隔的数据项。如果15 列为空，则表示它是前一张卡片的继续卡。如果在一行的末尾有一个用空格隔开的符号“&”，则表示下一行是该行的继续卡，数据可填写在 180 列。一个数据项必须在一张卡片上写完，不得跨到下一张卡片上。完全空白的一行则为两组卡片的分隔符。对任何给定的带有粒子标识符的类型卡只能有一张。需要整数的数据项必须填写整数，其它数据可填写为整数或浮点数以及MCNP能读的数据。为书写方便，可以使用四项书写功能：A.nR功能，表示将它前面的数据重复n次。例如：24R等同于22222B.nI功能，表
121、示在与其前后相邻的两个数之间，插入n个线性插值点。对于XnI Y的结构，如果X和Y是整数，且XY刚好是n+1的整倍数，则产生标准的整数插值，否则产生实数插值，但Y值直接存储。例如：1.52I3.01.52.02.532.0可能不精确而14I6123456都是精确定整数C.XM功能，它表示的数值为前面的数据乘上X。例如：112M2M4M2M11241632D.nJ功能，表示其后n个数据项使用缺省值。例如：DD.1（缺省值）1000DDJ1000如果nR、nI、及nJ项中缺省n，则假设n1。b)列输入格式列输入格式列输入块的格式： Si必须是MCNP卡片名字，它们必须全部是栅元参数、或者全部是曲面
122、参数、或者全部是其它参数。15 列672 列S1S2SmK1D11D12D1mK2D21D22D2mKnDn1Dn2Dnm(5)粒子标识符粒子标识符几个输入卡片都需要粒子标识符以区别中子、光子和电子的输入数据。这些卡片是：IMP、EXT、FCL、WWN、WWE、WWP、WWGE、DXT、DXC、F、F5X、F5Y、F5Z、PHYS、ELPT、ESPLT、CUT和PERT。粒子标识符由上述卡片名字后面的冒号、字母N、P或E组成。例如：中子重要性卡为IMP:N光子重要性卡为IMP:P(6)缺省值缺省值MCNP的许多输入参数都有缺省值，因此用户不需要每次都给出各个输入参量的值。当缺省值符合用户要求时
123、，便可不在输入文件中指定。当省略某张输入卡时，则该卡上的全部参数均使用缺省值。如果只想改变一张卡上的某一个特定参量时，则它前面的参量仍需指明，或者用nJ方式跳过前面那些使用缺省值的参量。例如：光子截断卡CUT:P3J-.10表示前3个参量使用缺省值，只改变第四项参量的值。(7)输入错误信息输入错误信息MCNP对输入文件出现的错误作广泛的检查，如果用户违反了输入说明的规定，将在终端上以及输出文件中打印致命错误信息，MCNP不再进行粒子输运计算，作业中断。第一个出现的致命错误是真的，而后面的错误可能不一定是真的，这取决于前面出现的致命错误的情况。若在MCNP运行行上指定FATAL项，则MCNP忽略
124、致命错误，照常运行。对于MCNP的警告信息，用户不应忽视，应搞清楚它们的含义。(8)检查几何错误检查几何错误MCNP在处理输入文件的数据时，不能检查一种非常重要的输入错误。即MCNP无法查出各栅元之间的重叠和空隙，只有当粒子丢失时，才会发现几何错误。即使如此，可能仍然无法准确判断错误性质。2)栅元描述卡栅元描述卡格式： jm d geom params或：jLIKE n BUT listj栅元号，1j99999，写在第15 列上。m栅元材料号，与材料卡（Mm）中的序号对应。m0为真空栅元。d栅元材料密度。正值为原子密度，负值为质量密度。对于真空栅元，该项缺省，不填写。geom栅元的几何说明。由
125、一系列带符号的曲面号经过布尔运算组成。params任选的栅元参数说明。n另一个栅元的名字（编号）。list描述栅元j和栅元n之间差别的栅元参数。在栅元的几何说明中，关于曲面的指向是一个很重要的概念。假定曲面S的曲面方程为f (x,y,z)0，则对于f (x,y,z)0的区域对于曲面S具有正的指向；而对于f (x,y,z)0的区域对于曲面S具有负的指向。正指向的区域用+S表示，“+”号可不写；负指向的区域用-S表示。栅元用各相关曲面的布尔运算表示，布尔算符包括交（用空格表示）、并（用冒号:表示）和非（用#表示）。缺省的运算顺序是先非，其次是交，最后是并，使用括号可控制布尔运算的次序。非运算有两种
126、形式：(1)#n，n是某个栅元号，#n表示一个由不在栅元n内的点组成的空间区域。(2)# ( -)，括号内是对某一个栅元进行描述的曲面栅元关系组，这一形式定义的几何区域由不属于括号内描述区域的点组成的空间。例如：30-12-4$定义栅元3#3$与下行相同#(-12-4)在栅元卡上可定义栅元参数以代替在输入文件中数据卡部分定义的栅元参数。格式为：关键词值。这儿允许的关键词是：带有粒子标识符的IMP、VOL、PWT、EXT、FCL、WWN、DXC、NONU、PD和TMP，以及关于重复结构的4个栅元参数卡：U卡、TRCL卡、LAT卡和FILL卡。在LIKE n BUT格式中，还有两个关键词MAT和R
127、HO，分别表示栅元的介质号和密度。例如：1016-4.21-23IMP:N=4IMP:P=8表示栅元10由曲面1的正面、曲面2的负面和曲面3的正面的交集组成，填充质量密度为4.2克/厘米3的16号材料。该栅元的中子重要性为4，光子重要性为8。例如：23-3.7-1IMP:N=2IMP:P=43LIKE2BUTTRCL=1IMP:N=103)曲面描述卡曲面描述卡(1)由方程定义曲面由方程定义曲面格式： jn a listj曲面号，1j99999，写在第15 列上。如果曲面号前有*号，则该曲面为反射面。n对应坐标变换卡TRn，表示该曲面是在辅助坐标系下描述的，而该辅助坐标系与基本坐标系之间的关系由
128、TRn卡给出。如果没有坐标变换，即曲面是在基本坐标系下描述的，则该项缺省。a曲面助记符。list曲面方程参数，110项，取决于曲面类型。参见MCNP手册，表3.1。(2)用点定义轴对称曲面用点定义轴对称曲面类型为X、Y或Z的曲面卡是用坐标点描述曲面而不是用方程系数描述。用这些卡描述的曲面必须是分别关于X、Y或Z轴对称的，并且如果该曲面是由多叶组成的，则指定的坐标点必须全都在同一个叶上。格式： jn a listj曲面号，1j99999，写在第15 列上。nTRn卡的号，如果没有坐标变换，则该项缺省。a字母X、Y或Z。list13对点的坐标。每一对坐标点定义这个曲面上的一个点。例如在一张Y卡上可
129、以给出：jYy1r1y2r2其中，()是第i点的坐标。给出的坐标点对数的不同，描述的曲面类型也不同。a)给出一对坐标，则定义一个平面（PX、PY或PZ）。b)给出二对坐标，则定义的是线性曲面（PX、PY、PZ、CX、CY、CZ、KX、KY或KZ）。c)给出三对坐标，则定义的是二次曲面（PX、PY、PZ、SO、SX、SY、SZ、CX、CY、CZ、KX、KY、KZ或SQ）。当用两点定义一个锥面时，只生成一个单叶锥面。曲面的指向与方程指定曲面（SQ除外）是一样的。对于SQ，远离对称轴的点具有正指向。而方程定义的SQ可以自由选取指向。(3)由三个点定义一般平面由三个点定义一般平面MCNP对用户指定的P
130、型曲面，将检查所给的数据个数，若是4项，则作一般斜置平面方程的系数理解，若多于4项时，便作为三维空间点的坐标值理解。每三个数定义空间一个点，MCNP将把它们转换成所需要的曲面系数以产生平面：AX+BY+GZD0格式： jn P X1 Y1 Z1 X2 Y2 Z2 X3 Y3 Z3j曲面号，1j99999，写在第15 列上。nTRn卡的号，如果没有坐标变换，则该项缺省。P该曲面卡的助记符。 （Xi,Yi,Zi）定义该平面的点坐标。3)数据卡数据卡在信息卡、栅元描述卡和曲面描述卡之后输入的是数据卡，数据卡可分为10类：(1)问题类型问题类型(2)几何卡几何卡(3)减方差减方差(4)源描述源描述(5
131、)计数描述计数描述(6)材料及截面描述材料及截面描述(7)能量及热处理能量及热处理(8)问题截断条件问题截断条件(9)用户数据数组用户数据数组(10)外围卡外围卡数据卡中，标识符必须从前5列开始填写。(1)问题类型（问题类型（MODE）卡卡如果不给出MODE卡，则缺省形式是MODEN，即缺省值是中子输运问题。格式： MODEx1 xi xi N，中子输运。P，光子输运。E，电子输运。(2)几何卡几何卡几何卡有以下几类：助记符卡片类型VOL栅元体积AREA曲面面积UUniverseTRCL栅元变换LAT格子FILL填充卡TR坐标变换坐标变换卡格式：TRnO1，O2，O3，B1，B2，B3，B4，
132、B5，B6，B7，B8，B9，Mn变换号，1n999，*TRn表示Bi是角度而非角度的余弦。O1，O2，O3坐标变换向量的位移。B1至B9坐标变换的坐标旋转矩阵。元素B1，B2，B3，B4，B5，B6，B7，B8，B9轴x,xy,xz,xx,yy,yz,yx,zy,zz,zM1，表示位移是辅助坐标系原点相对于基本坐标系的位移。1，表示位移是基本坐标系原点相对于辅助坐标系的位移。(3)减方差减方差MCNP运用以下卡片来减小方差：助记符卡片类型IMP栅元重要性ESPLT能量分裂和俄国轮盘赌PWT次级光子权重EXT指数变换VECT方向矢量定义FCL强迫碰撞助记符卡片类型WWE权重窗的能量或时间间隔W
133、WN权重窗的边界WWP权重窗的参数WWG权重窗生成器WWGE权重窗生成器的能量或时间间隔MESH重叠重要性网格权重生成器PD探测器贡献DXCDXTRAN贡献BBREM韧致辐射偏倚因子(4)源定义源定义助记符卡片类型SDEF通用源SIn源的信息SPn源的概率SBn源的偏倚DSn相关的源SCn源的注释SSW写曲面源SSR读曲面源KCODE临界源KSRC临界计算的源起始点ACODE特征值源通用源卡：格式(5)计数描述计数描述下列卡片用来记录计算结果：助记符卡片类型Fna计数类型FCn计数注释En计数能量间隔Tn计数时间间隔Cn计数方向余弦间隔FQn计数打印层次FMn计数乘子DEn/DFn剂量能量/剂
134、量函数EMn计数能量乘子助记符卡片类型TMn计数时间乘子CMn计数余弦乘子CFn计数栅元标志SFn计数曲面标志FSn计数片段划分SDn计数片段的体积/面积FUn子程序TALLYX输入TFn计数涨落打印DD探测器和DXTRAN诊断DXTDXTRAN球参数FTn计数特殊处理计数类型卡Fna格式：助记符类型说明Fn单位*Fn单位F1:(N、P、E)面流粒子MeVF2:(N、P、E)面通量粒子/cm2MeV/cm2F4:(N、P、E)体通量粒子/cm2MeV/cm2F5a:(N、P)点或环探测器通量粒子/cm2MeV/cm2F6:(N、P、N,P)平均沉积能量MeV/克109J/克F7:N平均裂变沉积
135、能量MeV/克109J/克F8:(P、E、P,E)脉冲MeV+F8:E沉积电荷电荷无(6)材料描述材料描述这组卡片用于指定在栅元中所使用的材料成分和使用那些截面数据。助记符卡片类型Mm材料成分DRXS离散反应截面TOTNU总裂变NONU裂变截断AWTAB原子量XSn截面文件VOID否定材料PIKMT次级光子产生偏倚MGOPT多群特征描述材料成分卡Mm格式：MmZAID1fr1ZAID2fr2keyword=valueZAIDi材料中第i种成份的截面数据，ZZZAAA.nnX或ZZZAAAZZZ是元素的原子序号，AAA是原子量，nn截面库标识号，X是数据分类。fri 材料中第i种成份的原子的分量
136、（负值表示重量比例）。AAA000表示自然元素。(7)能量和热处理方式指定能量和热处理方式指定这组卡片用于控制MCNP的能量以及其它物理状况。助记符卡片类型PHYS能量物理截断TMP自由气体模型热温度THTME热时间卡MTmS(,)材料卡(8)问题截断卡问题截断卡这组卡片在初始运行或接续运行的输入文件中均可使用，用于终止粒子的历史或中断计算。助记符卡片类型CUT截断粒子历史ELPT栅元能量截断NPS历史数截断CTIME计算时间截断(9)用户数据数组卡用户数据数组卡MCNP在其COMMON变量中定义了两个数组IDUM（整数）和RDUM（浮点数）供用户使用，每个数组可存放50个数据。这组卡片为这两
137、个用户数组提供输入数据。a)IDUM，整型数组卡格式：IDUMI1，I2，In（1n50）b)RDUM，实型数组卡格式：RDUMR1，R2，Rn（1n50）(10)外围卡外围卡这组卡片为用户提供方便，不影响MCNP的计算。助记符卡片类型PRDMP打印及转储周期LOST丢失粒子DBCN调试信息FILES建立用户文件PRINT打印控制MPLOT运行期间打印计数PTRAC粒子径迹输出PERT微扰卡蒙特卡罗中心服务器信息IP地址：166.111.32.63166.111.32.74166.111.32.79目前安装了MCNP4C，以后将陆续安装Egs4、Geant4、Fluka等蒙特卡罗程序。中心电话：62784552联系人：范佳锦、武祯
