MCNP3B说明书下.doc

MCNP-3B/PC 程序技术说明书 西安交通大学 一九九八年五月 前言 MCNP是美国 Los Alamos实验室应用理论物理部(X部)的Monte Carlo小组（X-6组） 研制的用于计算复杂三维几何结构中的粒子输运的大型多功能蒙特卡罗程序。 MCNP可用于计算中子、光子以及中子-光子耦合的输运问题，也可计算临界系统（包括次临界及超临界）的本征值问题。 MCNP可以处理任意三维材料结构的问题，栅元的界面可以是平面、二次曲面及某些特殊的四阶曲面（如椭圆环面）。栅元中的材料可由任意多种同位素组合而成。 MCNP使用精细的点截面数据。考虑了ENDF/B-V库给出的所有中子反应类型。对于热中子，可选用自由气体及S(,)两种模型处理。对于光子，考虑了相干与非相干散射，并处理了光电吸收之后可能有的荧光发射及电子对产生后的就地轫致辐射光子。 MCNP的通用性很强，使用也较容易。它已为用户配置了多种标准形式的源，同时留有接口，允许用户定义自己的源。为方便用户对几何输入卡的检查，配备了几何绘图程序。MCNP的记数部分是精心设计的，除有标准类型的记数外，也为用户准备了接口，用户想要的任何物理量几乎都能够计算。MCNP中汇集了非常丰富的降低方差技巧，对截面数据也进行了广泛的收集。 MCNP在世界范围内得到广泛应用。MCNP程序已用于分析物理试验和反应堆设计；设计核保护无损检测系统、辐射屏蔽和核仪器；计算材料活化和磁聚变中子学；进行临界分析；从事物理保健问题的研究等等。它可很好地用于跟踪计算、决定辐射剂量、物理实验、宇宙辐射模拟及辐射探伤研究等。 MCNP-3B完全是按照美国国家标准学会(ANSI)的标准FORTRAN77书写的，以前程序文本中依附于 Los Alamos软件及操作系统的内容大都于予以替换。广泛使用可调数组；也免除了很多原有的限制（如栅元及曲面数的限制等）。因此，MCNP-3B的可移植性很强，可以很方便地安装在CRAY、IBM、CDC、VAX、PRIME、SIEMENS、RIDGE、APOLLO和 SUN等各类计算机上运行。本程序为我们移植到微机上的版本称为MCNP-3B/PC。 MCNP程序的理论基础及主要参考文献是Carter及Cashwell的书Particle-Transport Simulation with the Monte Carlo Method2。MCNP程序的最新版本是MCNP-4A6。 i目 录 第一章 几何描述与处理.1 1.l 栅元(Cell).l 1.1.1 用“交”运算定义栅元2 1.1.2 用“联”运算定义栅元2 1.1.3 用“余”运算定义栅元.3 1.2 曲面(Surface)5 1.2.1 虚设曲面.6 l.2.2 反射曲面.6 1.2.3 曲面的描述方法7 1.3 粒子飞行轨线的计算.7 1.4 体积与面积的计算.8 1.4.1 对称体积与面积的计算8 1.4.2 不可计算的体积与面积.10 1.4.3 非对称的体积、面积之随机估计10 第二章 核反应的物理处理.12 2.l 截面数据12 2.1.l 中子截面数据.12 2.1.2 光子截面数据.13 2.1.3 中子热S(,)截面数据.14 2.1.4 中子剂量截面数据.14 2.2 核反应的模拟.14 2.2.1 粒子权重.14 2.2.2 粒子轨迹.14 2.2.3 源的描述及处理.15 2.2.4 中子反应.15 2.2.5 光子反应.22 第三章 基本物理量的记数.26 3.1 界面的粒子通量密度与流密度.27 3.2 平均通量密度记数.28 3.3 栅元中的能量沉积记数.28 3.4 探测器通量.29 3.4.1 点探测器.30 3.4.2 OMCFE探测器.32 3.4.3 环探测器.33 第四章 MCNP精度的估计.37 第五章 降低方差技巧37 ii5.1 几何分裂与轮盘赌.39 5.2 粒子的截断处理.40 5.2.1 能量截断.40 5.2.2 时间截断.40 5.2.3 权截断.41 5.3 指数变换.41 5.4 强迫碰撞.42 5.5 权窗.43 5.6 俘获的隐式处理.43 5.7 在探测点上的通量.44 5.8 DXTRAN44 5.9 能量分裂与轮盘赌.47 5.9.1 分裂.47 5.9.2 轮盘赌.47 5.10 源的偏倚抽样.47 5.10.1 能量偏倚47 5.10.2 方向偏倚47 5.11 相关抽样.49 第六章 临界计算.50 6.1 初始处理.50 6.2 输运计算阶段.50 6.3 增值因子.53 参考资料.54 iii第一章 几何描述与处理 MCNP能够处理任意的三维几何结构问题，各栅元可用笛卡尔坐标系下的一阶、二阶曲面及某些特殊的四阶曲面（如椭圆环曲面）所界定。笛卡尔坐标系由用户随意定义，但一般都采用右手系。而且，对于轴对称系统，常取Z轴作对称轴（并非必须）。 1.1 栅元(Cells) MCNP对问题的几何描述采用了组合几何的思想。但与 MORSE等程序的习惯作法不同，它并不限于对有限几种预定形状几何体的组合，而是给用户以极大的灵活性。用于组合的几何体完全由用户通过输入卡自行定义。凡是以MCNP规定的曲面所界定的几何区域均可参与组合，组合的方法采用对空间点集的布尔运算：联（或）、交（与）和余（非）。 在定义栅元时，栅元对一个曲面的坐向是重要的概念。假如S=f(x,y,z)=0是所解问题中的一个界面方程。对几何空间中任意一个点(x,y,z)，若有f(x,y,z)=0，则称该点位于曲面s上。对于不在曲面s上的点，若有f(x,y,z)0，则称该点对于曲面s具有正的坐向；反之，f(x,y,z)y的一切点对该平面均具有正的坐向。 MCNP描述栅元的方法是：对每一栅元由用户给出一张栅元卡，在其上依次列出该栅元用户定义的序号、栅元中材料的序号、材料的密度，接着再列出若干曲面-栅元关系组；每组给出一个界定该栅元的曲面号。显然，该栅元由几个不同的曲面界定，便应有同样数目的关系组。栅元卡上列出的每一曲面，都将整个几何空间划分为两个区域，分别对这一曲面具有正的和负的坐向。而所列曲面号的正/负，恰好规定了该栅元选用的是正/负坐向区域。几何卡所描述的栅元就是由卡上选定的这样一些区域的交、联及/或余所组成。 几何空间进行分栅元的原则，不只是应使每个栅元中只含一种均匀材料（或真空），还可能需要从一些其它因素来考虑。如：要便于抽样及使用降低方差技巧（实施几何分裂及轮盘赌），要实现对题目的几何系统的确切描述，要满足对记数量的空间分布描述的需要，等等。但当对栅元描述使用了“联”运算，对栅元及曲面的记数使用了分段记数方式后，上这些额外的分栅元，既使不能消除也可大为减少。 图 l-l 对题目的几何系统做出几何描述时，务必慎重。不仅要避免错误，还应注意到节省机器计算费用。 如图l-l a 它的几何结构十分简单，只是由很多具有相同顶部与底部平面的平 1行圆柱体组成。若将阴影部分定义成一个栅元，它便有14个界面。按照MCNP的算法，对 于进入该栅元的粒子，要求出到达下一个碰撞点的距离，首先要逐个计算粒子当前飞行线与这14个界面的交点，这显然是很费机时的。对这一题目，若改用图1-1b的分栅元设计，将原来的阴影区剖分成几个小栅元，计算效率将提高近4倍！ 1.1.1 用“交”运算定义栅元 “交”运算符是隐式的，栅元卡上，在两个曲面一栅元关系组之间置以空格表示。 当一个栅元仅用“交”运算定义时，则对界定该栅元的每一曲面而言，此栅元中的所有点都必须对其具有相同的坐向。直观地讲，单纯用“交”运算定义的栅元，不允许存在凹角。如图l-2，点1与点2关于曲面3或曲面4都具有不同的坐向，因此仅用“交”运算无法使实线内的空间定义成一个栅元。一种权宜办法是在凹角顶画出一个平面（虚线），用它将实线内的空间人为剖分成两个栅元。按照栅元卡的书写规则，列出、栅元的栅元卡如下（三栅元都是真空的，且按二维考虑。）： l 0 l -2 -3 6 2 0 1 -6 -4 5 3 0 -1 2 3 4 -5 图l-2 图1-3 栅元表示实线以外的全部空间，它是一个外部栅元。这里对栅元的定义也只使用了“交运算。因栅元也有凹角，这里的定义显然是不合法的，要确切地定义它，必须使用“联”或“余”运算。如果外部栅元的重要性指定为0，以此保证粒子一旦走入此栅元必定被杀死，则也是可以允许的。 如图1-3，若仅用“交”运算来定义有凹角的外部栅元，也是比较难处理的一个问题。最容易的方法是使用“联”、“涂”运算，对所有栅元都给以正确的描述。 1.1.2 用“联”运算定义栅元 在栅元卡上，用“：”指明“联”运算。使用“联”运算，可以对含有凹角的栅元产出正确的定义，甚至可以把一些互不连通的区域定义成一个栅元。 “联”与“交”都是二进制布尔运算，可以在栅元卡上联合作用，其运算次序是遵循先“交”后“联”的规则。同时也允许使用括号来改变规定的次序，因为括号内的项是先行处理的。 “联”、“交”运算的概念与通常点集中的“联”、“交”概念完全相同，下面给出几个例子说明之。 例1-1 图1-4中给出两个相交曲面1及2，则按栅元卡书写规则，整个阴影区应表示为-1:-2，即曲面1以内的区域与曲面2以内区域之和。这里，-1:-2描述了一个带有凹角的区域。 2 例1-2 图1-5中用两个平面1及2将几何空间分成了、两栅元，则栅元卡应写成： 1 0 1 -2 2 0 -l:2 这里假定两栅元都是真空区，若将的栅元卡写成： 2 0 -1 2 它定义的实际是图l-5中的双斜线区，因此未能正确地描述栅元。 图l4 图l5 例1-3 看图1-2，若定义实线以内的区域为栅元，实线以外的所有空间为栅元，那么，栅元卡则应写成： 1 0 1 -2 (-3:-4) 5 2 0 -5 : -1 : 2 : 3 4 这里，栅元与栅元都是含有凹角的栅元，由于使用了“交”与“联”，它们便得到了确切的描述。还应注意，在栅元的栅元卡上隐含着“先交后联”的原则，而栅元的栅元卡上由使用括号实现了特定的先“联”后“交”，而且此时是必须的。 从例1-3可以看出，两个栅元，若其中之一被另个所包围，则它们的栅元卡之间应有一个简单的自恰关系：一个卡上的“交”，在另一个卡上全部变成“联”，反之亦然。对同一曲面的坐向，二卡恰好相反。这种自恰性正是两个互余栅元的正确描述。 例1-3还告诉我们如下的普遍规律： i）凡是凸角的描述，均用“交”实现；凡是凹角的描述，均由“联”实现。使用“联”运算可以减少栅元个数，从而简化了对问题的栅元描述。 ii）若把“交”看作“乘”，把“联”看作“加”，并考虑括号的作用，则栅元卡上的运算层次便完全遵循着算术四则运算的同样法则，因此用户是容易掌握的。 1.1.3 用“余”运算定义栅元 “余”算将用于几何描述，并没有什么超越“交”及“联”的功能，它只是为某些栅元的几何描述提供了一种简写手段，而实际上它仍是隐含着若干“ 交 ”与“ 联 ”的运算。MCNP中以符号“”表示“余”算将。在直观上，可用“不在其内”的说法来理解“余”的含意。 算符“”有两种基本用法： i）#n，n是某个栅元号。n表示一个由不在栅元n内的点组成空间区域。 ii）#（），括号内是对某一几何区域进行描述的曲面-栅元关系组。这一形式定义的几何区域由不属于括号内描述区域的点组成空间。 例l-4 考虑一个立方体，其六个平面顺序编号为l，6。立方体内称为栅元，整个空间的其余部分定义为栅元。使用“联”、“交”算符进行描述，栅元卡可写为： 1 0 1 -2 -3 4 -5 6 3 2 0 -1:2:3:-4:5:6 这里假定两栅元都是真空区，若使用算符“”的第一种形式，便可将栅元表述成： 2 01 若使用“”的第二种形式，则应写成： 2 0(l -2 -3 4 -5 6) 显然使用第二种形式不如第一种形式简短。这里应注意，n与(n)具有完全不同的意义。若用 # (1)来定义这里的栅元则是不正确的。因为它描述的是对平面1具有负坐向的半空间。 MCNP对于第一种形式的“余”运算，要顺序做五项操作：i）移去符号“”； ii）接收栅元n栅元卡上的界面表，并用一对括号括上；iii）将其中原有的“交”全部改成“联”；iv）将其中原有的“：”全部换成一对背括号“)(”（相当于改成“交”）；v)将原有的曲面坐向全部变号。对第二种形式的“余”运算，除去上述第二项外，MCNP也将完成同样的步骤对于MCNP的这一处理方法，用户在使用“余”算符号时应作恰当的考虑。 图1-6 例1-5考虑一个立方体，其中有两个同心的球面（如图1-6）。立方体六个平面顺序编号为1，6。两个球面编号分别为7与8。 一组正确的栅元卡是： 1 0 -7 2 0 7 -8 3 0 1 -2 -3 4 -5 6 4 0#3#2#1 虽然外部栅元并不接触栅元、，但若定义其为： 4 0#3 则是错误的，因为这样就将图中的三栅元、及定义成一个栅元。 上述这组栅元卡虽然是正确的。但都存在很大的的缺点。因为按照 MCNP对“余”算符的处理，栅元的上述定义将被展开成： 4 0(-8:-1:2:3:-4:5:-6)(-7:8)(7) 界面表不仅过长，而且包含不必要的界面。然而在输运计算过程中，如果粒子进入栅元或在其中发生和碰撞后，要计算粒子飞行线与列出的每一界面的交点，显然其中增加了很多额外的计算量。因此不要像上述那样随意滥用“余”运算，将栅元的栅元卡改为： 4 0#(1 -2 -3 4 -5 6) 或直接写成： 4 0 -1:2:3:-4:5:-6 为好。 另外，在有些情况下，使用“余”算将会使MCNP计算栅元体积及曲面面积的必要条 4件遭到破坏（见有关体积、面积的计算），因此“余”运算必须使用得当。 有一种对几何做快速检查的方法，即去掉几何体中各栅元的材料，把各栅元暂时作成空腔，随后对这样的假系统做粒子输运，如有几何描述错误，很快便可暴露出来。 1.2 曲面（Surface） MCNP中共处理了29种曲面，包括一般的一阶、二阶曲面，如平面、球面、柱面、锥面及双曲面等，同时也包括了某些特殊的四阶曲面，如椭形圆环面等。 图1-7 如同双曲面，锥面方程描述的曲面也是分为两支的，一支的斜率为正，另一支为负（如 图1-7），以锥面为栅元界面时，通常只用到锥面中的一支。为区分锥面的两个分支，常用的方法是在锥点虚设一个平面（如图1-7上二虚线）以对于虚设平面的不同坐向来区分锥面的正、负两支。MCNP对此作了简化处理。当一个锥面只有其中的一支有用时，则可以定义该锥面的描述卡上添加一项+1或-1，用以指明该卡定义的是正斜锥或负斜锥。 MCNP计算中处理四阶曲面是困难的，因为求粒子飞行线与四阶曲面的交点，及计算体积与面积时求四阶曲面与其它曲面的交点，都需要求解四次方程的根。方程形式归结为： x4+Bx3+Cx2+Dx+E=0 由于舍人误差的问题，在计算机上求四次方程的根有一定难度。多年来，MCNP处理这一问题都要求30位(+进制数)的精度。即使这样，其存在的舍入误差也仍然需要用Newton-Raphsem迭代加以修正。若对四次方程使用单精度的解算公式，再跟之以牛顿迭代予以修正；是不可行的。因为，若提供给牛顿迭代的初值太不精确；而方程的根又靠得很近的话，迭代常常是发散的。 MCNP依照不同情况，采用不同方法求解四次方程。在计算体积与面积时。因为求解四次方程并不频繁，计算速度是不重要的，因此使用了全迭代的单精度方案。这一方案很精确，但也很慢。MCNP的输运计算部分，解四次方程的速度非常重要，如果仍使用上述方案，计算速度可能比双精度方案慢10100倍。因此在输运部分采用了另一种单精度（15位十进制数）方案。 MCNP输运部分求解四次方程的计算方案基本上是遵循Cashwell-Everett的方法3。当方程的根是散开的情况，这种改进的Newton-Raphson迭代很快就收敛了。但这一方法的关键是将重根或靠得很紧的根简单地扔掉了，因为重报意味着粒子轨线与环面（Torus）相切，可以近似假定粒子没有碰上这个环面。在极端少有的情况下，这一假定可能是不好的，会使粒子丢失，此时就用一种后备的较慢的全迭代方案。为了节省机时，在解四次方程的子程序中附加了一些精细的处理，如：慎重选取求根的有限区间；只要粒子是从环面出发时，不再解四次方程，而化简为解三次方程；对方程系数作全面检查，以便先确定是否存在正实根。这样就使得这一单精度方案比双精度方案的计算速度明显地快，而且有时还更为精确。 5由于方程解法的改进及其它一些修改，使得MCNP对环曲面的处理能力提高了，它可以处理任意具有对称轴（限于平行坐标轴）的椭形环回，而且可以计算环面切割的各类记数的体积与面积。 1.2.l 虚设曲面 描述一个栅元最基本的要求，就是对它表示哪栅元空间区域不能有任何含糊不清，即对于进入一个栅元的粒子，应当从各界面的坐向能够唯一确定它进入的是哪个栅元。如图1-8，如不指定一个虚设曲面，就少掉了一个栅元。假定图中几何体以y轴为旋转对称轴。当一个粒子由球内区域进人栅元时，也可以认为它是进入了栅元，因为按坐标坐向检验，老满足栅元的描述，则也会满足栅元的描述。为了消除这种不确定性，可以引入一个虚设曲面，如平面y=0。此时对栅元、栅元的描述，除要列出它们原有的界面表外，都必须多列出一项对曲面的坐向指定。栅元与栅元中的点对曲面而言，显然具有不同的坐向，因此进入栅元的粒子不可能再被判断为进入栅元。 虚设曲面不一定要是栅元的真正界面，但它也常常确实是个几何界面；而且可能是问题中别的栅元的真正界面。对一个特殊栅元的定义，还可能会需要不止一个虚设曲面。 图1-9将有助于进一步理解虚设界面的意义。图中，曲面1、3为同心球面，曲面2是环绕y轴的柱面。原意是要描述图1-9a的几何，但如果不加虚设曲面，则可能会被描述成图1-9b的几何，因为乍一看，很容易把栅元的定义错误地写成-1:2:3。倘若以平面y=0为虚设曲面，将其编号为4，则栅元的正确定义应是-1:(4:-2):3。 图l-8 图1-9 从用户角度看，虚拟曲面与其它曲面一样，它必须在曲面描述卡上给出定义，而且在某些栅元卡上以一个带数符的界面号被列出。唯一的限制是，虚拟曲面不得作为反射曲面使用。但对程序中的处理而言，若一个虚设曲面只在一张栅元卡上出现时，MCNP才把它纯粹当作虚拟曲面用，如果它出现在不止一张几何卡上时，它虽然仍起虚拟曲面的作用，但在TRACK子程序中，将要计算飞行线与它的交点，从而会多用一些机时。 1.2.2 反射曲面 任何曲面（虚设曲面除外），只要在其描述卡上的曲面序号前加一个“*”号，便被当作反射曲面使用。任何粒子打在反射曲面上便按镜面反射处理。使用反射界面的价值在于可以简化问题的几何描述及跟踪粒子的过程。但使用不当会带来麻烦，可能导致得不到正确的结果。例如，考虑一个由碳元素组成的立方体，边长为 10cm。在立方体底部放一个体积均匀的 5 MeV中子源，源栅元厚1cm，正好铺满立方体的底部。以 MCNP的面通量记数（F2）得到任一侧面的总通量为0.150(0.5)，而对侧面中心的点通量记数（F5），所得值为1.55E-3(1) ，这两个结果都是归一化到单位源中子的。若将该立方体竖着切掉一半，以 6切面作为反射手面，同时源栅元也切去一半，而其余条件均不变。此时计算的面通量为0.302(0.5)，是不用反射平面时的2倍；而侧面中心的点通量为2.58E-3(1)，与不用反射平面相比，少于2倍。 对同一个问题，使用或不使用反射平面，同一个记数量应当有相同的结果。上例中的面通量记数是正确的，因为在有反射平面的情况下，源体积被减去了一半，对记数做归一化处理时，源粒子的出发权本该是无反射平面时的一半。但探测器的记数在有反射平面的情况下，肯定是错了，而且总是算低了，其原因将在以后指出。因此，在有探测器记数或DXTRAN球的问题中，木能使用反射界面。 对一般的问题，使用反射曲面时，归一化因子应如何处理是与题目本身有关的，要仔细地推敲，保证计算结果的正确性。因此告诫用户，在没有弄得很清楚的情况下，决不要使用反射曲面。 1.2.3 曲面的描述方法 MCNP对曲面的描述有两种方法：i）按照表述曲面方程的需要，提供恰当的方程系数；ii）对以坐标轴为对称轴的曲面，指定其上已知几何点的坐标。 1指定方程系数 采用这一方法需先从曲面类型表中选定一种助记符号，然后算出定义这一曲面的方程系数。例如，对于中心位于点（4，l，3）半径为3.62的球面（其助记符号为S），其描述卡很容易给出为： S 4 1 -3 3.62 但有些曲面要指定它的方程系数就较为复杂。例如考虑一个椭球面，而它的对称轴不与坐标轴平行。此时只能将它作为一般H次曲面（助记符号为GQ）来描述，共需给出十项系数。直接算出这些系数并不容易，但可以把问题加以简化。先定义一个辅助坐标系，以棉球中心为原点，以椭球轴为坐标轴，该椭球在这一坐标系下便很容易给出方程系数。然后，通过使用坐标变换卡TRn，可以实现坐标系的平移及旋转，由计算机给出该棉球面在基本坐标系下的方程系数。 描述一个曲面，采用SQ（特殊H次曲面）型或GQ型，取决于曲面对称轴的定向。如果能以简单形式定义的曲面，决不应当用复杂形式去定义，否则会不必要地增加机时。 2指定几何点 指定几何点的方法有时对用户尤为便利。例如，用户从结构设计图上知道了曲面上的一些点或某些曲面交点的坐标，则通过指定这些点来定义曲面就很方便。尤其当几个曲面共同交于一点时，若用这一方法指定曲面，并使用了这一公共点，就能保证这些曲面交点的准确性。倘若用给定系数的方法定义曲面，就难于保证几个复杂曲面能准确地交于一点，有时这就会导致粒子的丢失。 用指定几何点来定义曲面有一些特殊的限制，而对用系数定义曲面的方法来说，则并不受这些限制。规定的限制条件是，用点描述的曲面必须以某一坐标轴为对称轴，而且它必须是唯一的真正的连续曲面。例如，在一个双曲面的两叶上都指定了几何点是不允许的，因为这样定义的曲面是不连续的。然而，若指定的点都在双曲面的同一叶上，则定义是有效的。 1.3 粒子飞行轨线的计算 对问题的几何结构给出了曲面及栅元的描述后，MCNP便可在这一几何系统下跟踪粒子的游动轨迹。 7 在一个碰撞点（或源粒子出生点）上，粒子当前所在栅元ICL的各界面是知道的，假设共有m个界面，其曲面号分别为J1, Jm。依据这些曲面的类型，通过求一次、二次或三、四次方程的最小正实根，首先要算出粒子沿其当前飞行线首次交于各界面的距离，设为L1， ，Lm。其中对某个n（lnm），若Ln=则表示粒子飞行线不可能穿过曲面Jn。取DLS=min(L1,Lm)。如果抽样得到粒子到达下一个碰撞点的飞行距离是PMF，且DLSPMF，则粒子便离开当前栅元（假定不考虑DXTRAN球），穿过曲面Jm(若有DLS=Lm， 1mm)进入邻近栅元。粒子究竟进入哪个栅元，要经过验证判定。首先要将粒子的几何坐标(x,y,z)沿飞行方向(u,v,w)向前移动距离DLS+，使新的坐标点（x,y,z）位于曲面J m的另一侧。假定栅元ICL的描述卡上列在曲面号J m后面的栅元号为I1，Ik。依据所列的次序，将对这k个栅元逐个进行验证。例如对栅元Ik（lkk），要对它的栅元描述卡上列出的所有界面（Jm除外），进行对点(x,y,z)的坐向检查，同时依据这些界面之间的“交”、“联”关系，确定点(x,y,z)是否位于栅元Ik中。一旦找到了这样的栅元，便确认粒子穿越曲面J m进入了栅元Ik。当在k个栅元均被否定时，当前粒子就被丢掉了。但当在J m后只列出一个相邻栅元号（即k=l），而且该栅元的重要性为“0”时，即使验证被否定，也认定粒子是进入了这个“0”重要性的栅元。对于未列出界面另侧相邻栅元号的情况，MCNP将对几何系统中凡有曲面号Jm为界面的所有栅元（ICL栅元除外）进行验证，判定点(x,y,z)究竟属于哪个栅元。 当确定了粒子越过曲面J m所进人的栅元后，粒子轨线的计算便可重复以上步骤，直到求出粒子的下一个碰撞点，或因进入了重要性为“0”的栅元而终止粒子的历史。 1.4 体积与面积的计算 对于蒙特卡罗输运问题，体积与面积的计算是很重要的，因为在很多情况下都需要知道问题中各几何区域的体积和表面积。例如，用单位体积的径迹长度估计粒子通量就需要知道几何区域的体积；用穿过曲面的单位面积粒子数估计粒子通量需要用到曲面的面积；计算栅元中单位体积或单位质量物质的受热量，需要给出栅元的材料密度和质量，从而也要用到栅元的体积等等。特别地，如果用户知道实际问题的材料装量，再通过核对程序算出来的装量，便是检查几何描述正确性的一个较好手段。这样也要求程序能正确地计算几何体积。 在近代蒙特卡罗输运程序中，体积与面积的计算并不是一件寻常的事情。特别是MCNP，它处理的几何界面包括了二阶甚至四阶曲面，而由这些曲面界定的各式各样几何区域，以“交”和/或“联”进行任意组合构成用户希望的各个栅元。定义的曲面，其方位可以互不相同，而且按记数要求可能还需分成若干段。而定义的栅元甚至可能是由几个互不相连的小栅元构成的。几何描述上的这种一般性，大大提高了MCNP应用的灵活性，但就体积与面积的计算而言，显然要求其计算程序设计得格外精细。 对于具有旋转对称轴（即使轴是斜置的）的体积与面积，MCNP是自动予以计算的。它也计算多面体栅元的体积与表面积，但侧平面必须精确相交。例如立方体的三个侧面本应只有一个公共的交点，不能因描述上的误差，使其没有严格的公共交点、对于其它情况，体积与面积的计算，MCNP只能根据用户的指示以随机方法实现。 1.4.1 对称体积与面积的计算 对于一大类广泛使用的几何，MCNP都建立了体积与曲面面积的算法，其性能如下： 计算曲面的面积。 计算体积的栅元，或者旋转体，其对称轴是任意的，甚至可以斜置；或者是多面体。 体积与多数面积的积分计算是使用公式，木采用数值积分方法。 8 可以计算由n个不相连的子栅元“联”成的栅元的体积。 为了记数的需要算出栅元及曲面（满足对称性要求）的分段体积与面积。 计算步骤大体如下： l）识别给定栅元的所有界面，二次曲面均按MCNP的笛卡尔坐标系统一表述为一般形式： Ax2+By2+Cz2+Dxy+Eyz+Fxz+Gx+Hy+Jz+k=0 （1-1） 对环形曲面(Torus)作为特殊情况另作处理。 2）对若存在辅助坐标系（x，y，z），使得给定栅元在该坐标系中是旋转对称的，则应确定出这一坐标系。 3）将给定栅元的所有界面经过旋转与平移变换到（x，y，z）坐标系，使曲面方程简化为二维柱坐标形式： ar2+br+cs2+ds+e=0 （1-2） 其中，r2=x2+y2;s=z。或者简单表示为：r=f(s)。 4）求出所有界面彼此的交点，只保留其中构成栅元角顶的点。识别交点是否是一个角顶的方法，如图1-10所示，这是两个相交曲面在（r,s）坐标系中的示意图。图中两个曲面将空间分成了四部分。可以想象得到，一个栅元总超不出这四个区域任意组合的范围。那么，两个曲面的交点定义成栅元的一个角顶，而且仅当 （1-3） =41)1(2*iiif不能被3整除。这里， （1-4） =10i区中。当栅元在区中；当栅元不在ii例如，若一个栅元出现在l、3区中，而不在2、4区中，那么1=3=1,2=4=0，从而有f=5，不能被3整除，因此曲面的交点便认为是该栅元的一个角顶。 图 1-10 图 1-11 5）按标准的积分公式，在角顶之间对曲面求积分： V （求体积） （l-5） =dsri2 +=dssrrAi2)(12 （求曲面面积） （l-6） 对环形曲面，Ai将必须用数值积分计算。 96）对所得积分值进行适当的加、减，便可得到每个栅元的体积及每个曲面的总有效面积。面积积分实际上是算了两遍，曲面的每侧各算了一次。由此可知，由旋转对称曲面界定的非对称栅元仍是可以考虑的。 1.4.2 不可计算的体积与面积 概括前节的内容，计算体积的栅元只能是两类，一类是完全由一些旋转对称曲面界定的栅元，而且这些曲面必须有一个公共对称轴，轴向任定；另一类是全部由平面界定的多面体。而要计算面积的曲面，至少有一侧都是一些可计算体积的栅元。 对于体积与面积之算法，MCNP本着宁可不算也不要算错的原则，设计得很保守。这样一来，在某些场合下，即使初看起来是符合上述要求的，但体积与面积仍不予计算。有的是因为舍入误差的问题，这种情况较少。有的是由于对栅元的描述使用了“余”算符，导致在曲面两侧都放上了一个不可计算体积的栅元。不论什么情况，对不予计算的体积与面积，MCNP都会打印出原因（若使用了PRINT卡）。此外，如果栅元是互不相连的，或者有个洞，MCNP会告诉你这个栅元共有几段，有助于用户鉴别是否是几何描述错误。 例l-6 图1-11给出一个立方体，其中置一球。球面编号为1，立方体六个面编号依次为2，7。如果将栅元卡写成： 1 0 -l 2 0 1 2 3 -4 -5 -6 7 3 0#1#2 则出现一个意外情况，此时曲面1的面积木予计算，原因是曲面1的两边都放有不可计算体积的栅元。曲面1外侧的栅元是不能算体积的，因为它既不是单纯由平面围成的多面体，也不是全部由旋转对称曲面所界定。曲面1的内侧是栅元，它的体积的确可以计算，但由“余”算符的间接意义，栅元也认作曲面1内的栅元，而它是无限几何体，体积无法计算。MCNP规定，两侧均放有不可计算体积的栅元的曲面，其面积便不予计算。倘若的栅元卡上不用#1，而直接写成-l，曲面1的面积便可计算了。 对于多面体，由几个平面形成的角顶必须是严格的公共的交点，才符合计算多面体体积与表面积的条件。在某些情况下，由于定义平面的方程系数有效位给得不够多，或者由于计算机的舍入误差，可能会导致MCNP拒绝对该多面体体积与面积的计算。 1.4.3 非对称的体积、面积之随机估计 对于没有旋转对称性质的栅元及曲面，MCNP可用轨线方法对体积与面积做出随机估算。一种做法是使用球面源。将要计算体积与面积的区域用一个适当半径的球面围住，在它上面置一个单位强度的源（中子或光子），粒子能量随意给。所有栅元都改成真空。在源球 以内的栅元，重要性均置为“l”；球外各栅元，重要性都置为“0”。对计算体积的栅元指定F4型记数（经迹长度估计），对计算面积的曲面指定F2型记数（面通量估计）。随后便可运行，经过“输运”计算，分别从F4型及F2型记数的输出结果便得到了体积与面积的估算值。源球的半径要给得适当，给得过大，统计误差也大。当然，用球面划出所关心的几何区域可能是不方便的限制。 另一种做法是使用平面源。对于某些情况，它可能更为方便。此时要设置一个平面源，它从一个圆形窗口发射源粒子，方向垂直于圆窗平面，源粒子初始权为R2。圆窗半径R应选择恰当，既应大到足以使平行线束源能映射到整个所关心的区域，又不可过大而不必要地增大了统计误差。在使用平面源计算体积与面积时，应使源束映射范围内的栅元重要性均置为“l”，而以外部分则置为“0”。其余处理与使用球面源时相同。 如果输入文件中使用了VOID卡，MCNP将把所有栅元自动改换成真空栅元，并将所 10有非“0”重要性的栅元，重置重要性为“l”。因此临时使用它，可减少用户对原有输入文件的改动。 以随机方式算出来的体积与面积，用户可列入VOL、AREA卡，供正式输运计算使用。 应注意到，对于MCNP能够解析计算体积与面积的栅元及曲面，此时打印出来的F4 及F2型记数量将近似为“l”。 11第二章 核反应的物理处理 2.l 截面数据 MCNP使用的截面数据，以广泛的原始数据库为基础，经LOS Alamos的T-2小组加工处理而成。MCNP程序包含有五种核截面数据：（l）连续能量中子截面数据；（2）离散中子截面数据；（3）光子截面数据；（4）中子S(,)热截面数据；（5）中子剂量截面数据。1的附录F相当详细地给出了MCNP使用的截面数据文件的格式说明。 MCNP读取的截面文件可以有三种形式，一般采用以下两种形式： i）标准的有格式文件，它是按80列的卡片映像BCD格式给出。它是便携的，但读取时很慢，而且文件很庞大。也可使用辅助程序MAKXSF进行加。 ii）标准的无格式文件，它是一个二进制数据文件，是由用户使用辅助程序MAKXSF将第一种形式的文件加工而成。它显然是不便携的，但读取很快。 2.1.l 中子截面数据 中子截面数据中子截面库主要有四个来源： ENDF/B库、ENDL库、AWRE（Atomic Weapons Research Establishment）库以及特殊的Los Alamos子库。 MCNP使用的截面数据，包括整个共振区在内，都是点截面形式。MCNP考虑了特定库（例如ENDF/B）中所有的各种