实用MCNP3B教程

使用MCNP3B教程 大型通用中子-光子联合 输运蒙特卡罗模拟程序 郑华 编 大庆生产测井研究所 1998年10月 前言 本教程用于培训那些曾经、正在或将要从事核测井研究工作的人员，使他们掌握一种在计算机上模拟核测井过程的技术。 目前研究核测井问题的计算机模拟方法有两种，一是扩散理论方法，二是蒙特卡罗方法。扩散理论法对核测井过程建立玻尔兹曼扩散方程，然后求解。由于测井问题比较复杂，这个含有多个变量的微积分方程一般不能得到解析解，还要用数值方法进行求解，在求解时对粒子运动方向的角坐标可用级数展开或取离散角度处理，对空间和时间变量多采用差分处理。蒙特卡罗方法则对给定的问题建立相应的随机抽样模型，并用以一系列随机数跟踪大量粒子历程的方法完成对粒子输运的模拟。与确定论相比，蒙特卡罗方法能更好地适应复杂的几何条件。 蒙特卡罗方法模拟辐射输运的思想在40年代由美国Los Alamos实验室的科学家提出，76年开发出了通用程序MCNP。从80年代开始，美国一些核武器研究实验室和大学受雇于各大测井公司进行核测井数值模拟研究，编写出一批核测井专用蒙特卡罗模拟程序。专用程序的计算时效高，但针对性强，后来随着计算机速度的提高，一般人们还多使用通用程序进行计算。目前最常用的两个通用程序是MCNP和英国原子能署的MCBEND。 MCNP具有较强的通用性，可处理任意三维几何问题，提供了多种源分布和记数方式，使用精细的点截面数据。MCNP3版(83年)和3A版(85年)发行后，就成为用蒙特卡罗方法模拟核测井的最流行的通用程序，它解决了特征伽马谱线的问题，可以较好地模拟中子和光子联合输运，使用的主要核数据库是ENDF/B-4。3B是89年发行的版本，对3A版做了许多改进，并增加了几何重复定义功能。91年问世的4版可以联合模拟中子、光子、带电粒子(离子)的输运，可模拟探测器的测量结果。4版与以前版本相比在性能上又有较大改善，程序代码时效有显著提高，使用了更新的ENDF/B-6评价核数据库。 MCNP程序3B版传入中国已多年，可从中国原子能科学研究院计算机应用研究所得到拷贝。金文绵、李素梅等人将MCNP3B IBM 3033版移植到微机上，并翻译了原说明书第三章，对中国用户提供了很大帮助。 本教程主要来源于MCNP3B说明书的前三章，并加上部分本人应用心得。涉及数学方法和物理基础的部分仅作简要介绍，其中多数定义并不严格，目的是让读者易于接受。有关数学基础和蒙特卡罗方法在粒子输运问题中应用方面的详细知识可参阅裴鹿成和张孝泽的专著，中子、光子与物质相互作用的知识可从核物理教材中了解。 鉴于本教程有特定的应用范围，其中并不全面覆盖使用MCNP3B程序的所有知识。有关计算裂变核反应(如临界问题)、如何修改源程序和用户编写输入子程序(如定义用户源、用户记数)的内容本教程中一概不予介绍，还有一些我们不常用的MCNP功能(如几何描述重复定义、坐标变换等)也不介绍。对这些内容感兴趣的读者可以阅读MCNP3A和3B使用说明书。 由于本人才疏学浅，使用MCNP程序的经验又不多，编写过程也比较匆忙，本教程中不免存在错误和疏漏。如读者发现其中有差错之处或觉得部分叙述晦涩难懂，请与编者联系，电话是5970055(家)、5592052(办)。 郑华 1998/10/5 目录 页码 第一章 入门 1 §1.1 用蒙特卡罗方法模拟粒子输运 §1.2 MCNP简介 A. 程序的发展过程和应用领域 B. 程序的特点 §1.3 MCNP输入文件 2 A. 输入文件的基本形式 B. 一个简单的例子 C. 接续运行的输入文件 D. 卡片格式 3 E. 输入错误信息 §1.4 如何运行MCNP3B微机版 A. MCNP3B程序包的主要文件 B. 安装和运行MCNP3B的操作步骤 4 第二章 几何 §2.1 基础知识 A.​ 概述 B.​ 栅元定义中的一些概念 C.​ 有关曲面的一些知识 §2.2 几何描述卡 5 A.​ 栅元卡 B.​ 曲面卡 §2.3 有效地构建几何 6 A.​ 定义栅元的原则 B.​ 检查几何错误 第三章 数理基础 §3.1 物理 A.​ 粒子权重 B.​ 粒子径迹 C.​ 中子与物质作用 D.​ 光子与物质作用 8 §3.2 记数 9 A.​ 面流量记数 B.​ 通量记数 C.​ 栅元能量沉积记数 D.​ 探测器记数 E.​ 记数精度 10 §3.3 减小方差技巧 11 A.​ 统筹考虑 B.​ 能量截断 C.​ 时间截断 D.​ 几何分裂和轮盘赌 E.​ 能量分裂和轮盘赌 F.​ 暗含俘获和权重截断 G.​ 强迫碰撞 12 H.​ DXTRAN I.​ 源变量偏倚 J.​ 权重窗口 13 K.​ 指数变换 L.​ 相关抽样 M.​ 点探测器 第四章 数据卡 §4.1 问题类型卡 14 §4.2 栅元参数和曲面参数卡 A. IMP 栅元重要性卡 B. VOL 栅元体积卡 C. AREA 曲面面积卡 D. PWT 光子产生权重卡 E. EXT 指数变换卡 F. VECT 矢量输入卡 G. FCL 强迫碰撞卡 15 H. WWE 权窗能量卡 I. WWN 权窗边界卡 J. WWP 权窗参数卡 K. WWG 权窗产生器卡 L. WWGE 权窗产生器能量卡 M. PDn 探测器贡献卡 N. DXC DXTRAN贡献卡 16 §4.3 源的描述 A. SDEF 通用源卡 B. SIn 源信息卡 17 C. SPn 源概率卡 D. SBn 源偏倚卡 E. DSn 相关源分布卡 18 F. SCn 源注释卡 §4.4 记数方式的指定 19 A. Fna 记数类型卡 B. FCn 记数注释卡 20 C. En 记数能量卡 D. Tn 记数时间卡 E. Cn 记数余弦卡 F. EMn 记数能量乘子卡 G. TMn 记数时间乘子卡 H. CMn 记数余弦乘子卡 I. DEn/DFn 剂量能量/剂量函数卡 J. CFn 记数栅元标志卡 K. SFn 记数界面标志卡 L. FSn 记数片段划分卡 M. SDn 记数片段的体积/面积卡 N. FQn 记数打印层次卡 21 O. TFn 记数涨落打印卡 P. DD 探测器和DXTRAN诊断指定卡 Q. DXT DXTRAN的参量卡 R. FTn 记数特殊处理卡 §4.5 的指定 22 A. Mm 材料成份卡 B. DRXS 离散反应截面卡 23 C. AWTAB 原子量卡 D. VOID 否定材料卡 §4.6 能量和热处理方式的指定 A. PHYS 能量物理截断卡 B. ESPLT 能量分裂和轮盘赌卡 C. TMP 自由气体热处理卡 D. THTME 热时间卡 E. MTm S(,)材料卡 24 §4.7 问题截断卡 A. CUT 截断卡 B. NPS 历史截断卡 C. CTME 计算时间截断卡 §4.8 外围卡 A. PRDMP 打印及转储周期卡 B. LOST 丢失粒子卡 C. DBCN 调试信息卡 D. PRINT 打印控制卡 25 §4.9 MCNP输入文件综述 A. 输入卡 B. 存储限制 第五章 经验 26 §5.1 一般应用步骤 §5.2 需注意的问题 附录 连续能量中子截面库ENDL851数据目录 中子热截面库BMCC1数据目录 离散中子截面库D91数据目录 27 光子截面库MCPLIB1数据目录 特殊材料S(,)热截面库TMCC1数据目录 参考文献 第一章 入门 MCNP3B版是一套模拟中子和光子联合输运的通用蒙特卡罗程序，它具有连续能量、三维几何和与时相关的处理能力。本章仅简要介绍该程序的使用方法，详细内容将在后续各章中叙述。 §1.1 用蒙特卡罗方法模拟粒子输运 图1.1 单个中子随机历程示意 中子和光子在物质中输运的宏观表现是大量粒子与原子核微观作用的平均结果，蒙特卡罗方法通过逐一模拟和单个粒子的历程来求解输运问题。要得到比较合理的平均结果需要跟踪大量的粒子，至于单个粒子在其生命中的某一阶段如何度过，可以在已知统计分布规律的前提下通过抽取随机数来决定。这就象掷骰子赌博一样，因而得名蒙特卡罗方法。 图1.1显示了模拟中一个中子射入物质后的随机历程。首先根据中子与物质作用的物理规律(分布函数)，选取一个随机数决定中子在何处与原子核碰撞，本例中在1点碰撞；然后再用抽取随机数的方法决定中子与原子核发生了哪种反应，这里抽出的是非弹性散射反应；散射中子的能量和向哪个方向飞行也是用抽取随机数的方法从已知分布函数中决定的；碰撞过程中是否产生光子以及光子的能量、飞行方向等参数还是要通过抽取随机数从已知分布中决定，这里产生了一个光子。跟踪光子，确定它在7点与原子核碰撞并被吸收。散射后的中子在2点与原子核发生(n,2n)反应，其中一个出射中子射向探测器，另一个中子在3点被吸收。在2点的碰撞还产生了一个光子，它在5点又与原子核发生了一次散射反应，并离开物质。这一入射中子的历史过程结束了，有一个中子到达了探测器，感兴趣的结果被记录下来。跟踪越来越多的入射粒子历程后，平均结果就能反映出宏观效果。通过以上描述，读者不难领略蒙特卡罗方法如何通过跟踪粒子历程的方法计算问题，也了解了随机数在蒙特卡罗计算中的独特作用。 §1.2 MCNP简介 A. 程序的发展过程和应用领域 40年代美国Los Alamos实验室的Fermi、von Neumann和Ulam等人提出用蒙特卡罗方法模拟辐射输运的思想。47年Fermi发明了第一台用蒙特卡罗方法计算中子链式反应的机器。从50年代开始，von Neumann领导一个小组研究输运问题的蒙特卡罗处理方法，编写出模拟中子输运的程序MCS。63年蒙特卡罗方法描述语言化。65年完成的中子输运程序MCN有了很大改进，使用了标准的截面库，并且具有复杂几何描述功能。后来，Los Alamos实验室又开发了模拟光子输运的程序MCG(高能)和MCP(能量低至1keV)。73年MCN和MCG合并成MCNG，为MCNP的雏形。MCNP于76年开发成功，77年6月发行。3B版之前还发行过1A、1B、2、2A、2B、2C、2D、3和3A版。 MCNP软件包(a general Monte Carlo code for Neutron and Particle transport)是一套通用的、模拟三维空间中连续能量的中子、粒子联合输运的程序，其名字早先来源于the analog Monte Carlo method for Neutron and Proton's transport的缩写。MCNP3版(83年)和3A版(85年)发行后，这一软件在核测井领域逐渐成为最流行的通用程序，程序模拟结果和模型井实验结果较好地吻合，此时程序使用的主要核数据库是ENDF/B-4。88年发行的3B版程序增加了几何重构功能。91年4版问世，加入了模拟带电粒子(离子)输运部分，可以模拟探测器的测量结果，使用了新的ENDF/B-6评价核数据库。 MCNP程序的应用范围十分广阔，主要包括：反应堆设计、核临界安全、辐射屏蔽和核防护、探测器的设计与分析、核测井、个人剂量与物理保健、加速器靶的设计、医学物理与放射性治疗、国家防御、废物处理、射线探伤等。 B. 程序的特点 MCNP3B版程序能处理的中子能量范围是10-11MeV至20MeV，光子能量范围1keV至100MeV。程序中物理量的单位规定为： 物理量 单位 长度 cm 能量 MeV 时间 shake, =10-8sec 温度 MeV (kT) 原子密度 atoms/barn-cm 质量密度 g/cm3 截面 barns (10-24cm2) 加热量 MeV/collision 此外，原子质量按中子质量为1.0计算，这种单位下阿佛伽德罗常数是0.59703109；程序运行时间以分钟为单位。 MCNP程序的源代码是用FORTRAN语言编写的，3B版主程序有6万多行。程序包中携带了大量的核反应数据库文件，我们得到的库文件约占80Mb磁盘空间。 MCNP3B版具有很强的通用性，主要体现在： 1.​ 可以处理任意三维几何结构的问题。在输入文件INP中，空间被曲面(surface)分割成相互邻接的区域，称为栅元(cell)，可以给栅元填充各种物质。栅元的界面可以是平面、二阶曲面或某些四阶曲面(如椭圆环状面)。有关MCNP几何描述的知识将在第二章介绍。 2.​ 可以模拟中子输运、光子输运和二者联合输运。 3.​ 用户可以非常方便地在任何位置指定体源、面源、线源或点源，设置源粒子位置、能量、时间、飞行方向等参数的分布。源的指定方法将在第四章中介绍。 4.​ 程序提供多种记录模拟结果方法，包括通过某一界面的粒子流量或通量、进入某一栅元的通量、沉积能量和点通量。模拟结果在MCNP中称为记数(tally)，可以按位置、能量、时间、粒子来向和粒子种类记数。有关记数的指定方法将在第四章中介绍。 5.​ 程序包携带了大量核反应数据库文件，包括连续和离散的中子截面库、光子点截面库、热中子点截面库等，几乎可对所有天然物质进行计算。程序能比较精细地模拟中子和光子输运过程，并对一些特定的物理过程允许用户选择使用哪种方式进行处理，如对热中子处理可选用自由气体模型或S(,)模型，对低能光子处理可以考虑或忽略相干散射等。 6.​ 为了提高计算时效，给用户提供了许多可选用的减小方差(variance)技巧，主要包括：重要抽样、权重截断和轮盘赌、时间和能量截断、模拟俘获、指数变换、强迫碰撞、能量分裂和轮盘赌、源的偏倚、点探测器记数、确定论输运、权窗等。有关减小方差技巧的理论基础在第三章中阐述，应用方法在第四章中介绍。 7.​ 用户可通过设置源粒子数或运行时间来通知程序何时终止运行，还可以在原有计算结果的基础上接续运行程序。一些结果不会因计算的意外中断而丢失。 8.​ 在输出文件OUTP中给用户提供丰富的信息，包括输入列表、使用的截面表、粒子生成和丢失表、栅元中的粒子活动情况、中子诱发光子表、记数和记数涨落表等，还可以根据用户要求给出其它信息。 9.​ 提供了简单的问题调试工具。 §1.3 MCNP输入文件 若用户不修改源程序，MCNP的输入文件包括截面数据库文件、XSDIR文件、INP文件等。INP文件是用户要填写的主要输入文件，一般把该文件特指为输入文件。OUTP是MCNP的主要输出文件，其它输出文件还有转储文件RUNTPE、运行信息文件OUTPUT等。 A. 输入文件的基本形式 INP文件有初始运行和接续运行两种形式，用户须在INP中描述问题的几何、材料、源、记数和其它要求。INP由一些被空行分隔的输入块组成，主要的输入块是信息块、标题和栅元块、曲面块和数据块等。输入块又由一些被称为卡的输入行组成。初始运行输入文件的格式如下： 信息卡 } 可选 空行分隔符 标题卡 栅元卡 . . . 空行分隔符 曲面卡 . . . 空行分隔符 数据卡 . . . 空行结束符 其它 信息卡的1~8列应填写MESSAGE:，后面跟着用空格分隔的参数项。可用A=B参数项更改输出文件名，如OUTP = MYOUT。信息块是可选的。 在信息块之后的第一行是问题的标题卡，它仅限于一行，占用1~80列，可以是任何信息，将作为OUTP文件中各个输出表的标题被复制。 用户在栅元块和曲面块中描述问题的几何。栅元由栅元卡描述。空间必须由彼此相邻的栅元填满，栅元之间不能重叠，也不能出现无栅元的空区，否则会出现错误。构建栅元的曲面由曲面卡定义，曲面卡在曲面块中给出。曲面卡和栅元卡的填写方法见第二章。 曲面块之后是数据块，在数据块中用户描述源、记数方式、材料等。数据卡在第四章详细介绍。 无论数据块后是否有空行结束符，MCNP都能运行。用户可以把希望保留的一些附加信息写在数据块的空行之后，MCNP会将它们复制到OUTP文件末尾。 B. 一个简单的例子 为说明如何填写INP文件，这里例举一个简单问题。如图1.2所示，在一个边长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 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 N IMP:N 1 1 1 0 SDEF POS=0 –4 –2.5 ERG=14 F2 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 图1.2 例子的几何示意图 NPS 10000 空行 本例中没有信息块，第一行是标题卡，之后至空格前为栅元块。栅元卡上依次填写栅元号、材料号、密度和构成栅元界面的曲面号(带正负号)，这里定义了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卡 运行源粒子数10000 以上例子仅用于说明INP文件格式，有关各输入卡的详细内容，读者在此不必深究。 C. 接续运行的输入文件 为了继续一个早先被终止的计算，或想对早先终止的计算进行重新编辑打印(指向文件输出，即PRINT语句，不是向打印机输出)，可以接续运行MCNP。 接续运行INP文件中必须有信息块，信息卡中必须有C或C m参数，其中m表示从第m次转储开始接续运行，m缺省表示从最后一次开始。接续运行输入文件的格式为： MESSAGE: C m 空行分隔符 CONTINUE 数据卡 . . 空行结束符 其它 此时INP中不能有标题卡、栅元块和曲面块。数据块第一行必须是CONTINUE，允许使用的数据卡只是NPS、CTME、PRINT、FQ、DD、PRDMP、LOST和DBCN卡。NPS卡参数为累积运行的总源粒子数，CTME卡参数相对于本次运行的开始时间。如果给NPS置负参数，表示仅从转储文件中提取数据重新打印。 接续运行的另一个必要条件是上次运行的转储文件RUNTPE还没被删除。因每次运行都生成新的转储文件，若想对一次运行结果多次重新打印，则应把RUNTPE文件做备份。各次转储时运行的源粒子数和时间可从OUTPUT文件中查到。 D. 卡片格式 INP文件的每一行都限于使用1~80列。卡片都可以按行填写，数据卡也允许按列填写。 a) 行输入格式 通常卡片的1~5列用于填写栅元号、曲面号或数据卡的助记名，6~72列填写卡片参数，73~80列为注释，$符号之后也为注释。序号或卡片助记名可以写在1~5列的任何地方。带有粒子标识符助记名可能需要5列以上，但冒号必须写在6列以内。如果1~5列空白，表示本行为上一行的接续行。6列之后可以写数据项，多个数据项之间用空格分隔。一个数据项必须在一行上写完，不得接续到下一行。相同编号的卡片只能有一张。 填写卡片参数时可以使用以下书写功能： 1.​ nR功能，表示它前面的数据再重复n遍。 2.​ nI功能，表示与其前后相临两个数之间插进n个线性插值点。 3.​ xM功能，表示数值等于它前面数据的x倍。 4.​ nJ功能，表示从本项开始的n个数据使用缺省值。 b) 列输入格式 列输入只能用于数据块中，对栅元参数和源的描述比较有用。按行输入的栅元重要性、体积、权窗等数据项可读性较差，而且增加或删除栅元时要在行输入卡上仔细寻找相应项。列输入的可读性有很大提高，删除或增加与某一栅元相对应的数据项时也比较方便。 列输入格式的第一行以#开始，#可以放在1~5列的任意位置，卡片助记名逐个放在该行6列以后，在这些助记名之下按列给出数据项。同一个列输入格式块中的卡片必须是同一类卡片，比如都是栅元参数卡、都是曲面参数卡或都是源参数卡等，在#号下面的1~5列放置栅元号、曲面号或源分布号。 c) 粒子标识符 粒子标识符是卡片助记名的一部分，:N表示中子，:P表示光子，有时也能遇到表示中子-光子联合输运的:N P。下面一些数据卡需要粒子标识符：IMP、EXT、FCL、WWN、WWE、WWP、WWCE、DXT、DXC、Fna、PHYS、ESPLT和CUT。 d) 缺省值 MCNP许多输入卡的参数项有缺省值，用户不必每次都给出这些参数，如果卡片输入项有固定顺序，可以使用nJ功能跳过n个输入项。如果卡片上所有数据项都想缺省，只给出卡片助记名即可。有些卡片不给出也有缺省值，如MODE N卡就可以省略。 E. 输入错误信息 MCNP对输入文件作广泛的检查，如果用户违反输入规定，将在输出文件OUTP中打印致命错误信息(FATAL ERROR)，程序停止运行。第一个致命错误是真实的，后面的致命错误不一定是真实的。 MCNP还给出许多警告性的提示(WARNING)，对这些提示用户也不可忽视。 如果遇到任何代码损坏，如零作除数，MCNP将给出BAD TROUBLE信息，并终止运行。 §1.4 如何运行MCNP3B微机版 通常编辑INP文件后，键入RUN386 MCNP并回车，就可运行MCNP程序。3B微机版源程序用NDP FORTRAN 386编译器编译，要求386以上的微机和DOS3.3之后的操作系统，不能在Windows系统或保护方式下运行。 A. MCNP3B程序包的主要文件 1.​ MCNP.DAT：MCNP3B源程序MCNP.F的缩写型文件。文件中所有公用说明语句只列出一次，便于修改。 2.​ PRPR.F：辅助程序，把MCNP.DAT展成MCNP.F。 3.​ MAKXSF.F：该程序用于对NCMP所需截面数据加工处理，从原有的截面数据文件中挑选若干常用的截面数据，生成新的较小文件，可以节省磁盘空间和缩短程序运行时的数据读取速度。 4.​ SPECS：是MAKXS.F程序的输入文件，用来指定加工截面数据文件的任务清单，包括加工格式、输出文件名和核素。 例如，某SPECS文件的内容是： XSDIR1 XSDIR2 MCPLIB2 2 2048 512 0 1000.01P 2000.01P 1600.01P BMCCS2 2 2048 512 0 1003.03C 2004.03C 92235.11C 表示原截面库目录文件名是XSDIR1，新截面库目录文件是XSDIR2(为二进制的直接存取文件)，MCPLIB2中包含H、He和S的光子截面数据，BMCCS2中包含3H、4He和235U的连续能量中子截面数据，它们必须在XSDIR1目录规定的数据库中存在。 5.​ XSDIR1：它也是MAKXS.F程序的输入文件，是现有全部截面数据库目录。 6.​ BMCCS1：连续能量中子截面数据。 7.​ D91：离散能量中子截面数据。 8.​ MCPLIB1：光子截面数据。 9.​ TMCCS1：热处理中子的数据。 10.​ ENDL851：连续能量中子截面数据。 11.​ 531DOS1：中子剂量数据。 12.​ 532DOS1：中子剂量数据。 13.​ LLLDOS1：中子剂量数据。 14.​ PROB1~PROB18：例题。 15.​ OUTP1~OUTP18：例题结果。 B. 安装和运行MCNP3B的操作步骤 1.​ 拷贝原盘上的文件PRPR.F、MAKXSF.F、MCNP.DAT、XSDIR1及所需截面数据库文件。 2.​ 安装NDP FORTRAN 386编译器。 3.​ 建立MCNP输入文件INP。 4.​ 建立MAKXSF的输入文件SPECS。 5.​ 编译PRPR.F和MAKXSF.F： F77 PRPR.F F77 MAKXSF.F 6.​ 运行MAKXSF.EXP： RUN386 MAKXSF 生成新的截面库目录文件XSDIR2和新截面数据库文件。 7.​ 展开MCNP.DAT： RUN386 PRPR 生成MCNP.F文件。 8.​ 编译运行MCNP.F： COPY XSDIR2 XSDIR F77 MCNP.F RUN386 MCNP 运行结束生成结果文件OUTP、转储文件RUNTPE、运行信息文件OUTPUT等。 第二章 几何 本章介绍MCNP几何构建方法，INP文件中栅元卡和曲面卡格式，以及用户定义栅元时应注意的事项。 §2.1 基础知识 A. 概述 图2.1 迪卡尔坐标系 用户在三维正交坐标系(见图2.1)下定义曲面，用曲面围成栅元，在栅元中填充物质。栅元可由一些曲面包围区域的交集和并集构成。曲面可通过设置解析方程参数或指定曲面上一些已知点定义。 B. 栅元定义中的一些概念 a)​ 坐向 (sense) 定义栅元时，一个重要的概念是坐标点相对于曲面的坐向。假设s=f(x,y,z)=0是一曲面S的方程，对于任意一点(x’,y’,z’)，代入s=f(x,y,z)后，若s=0则该点在曲面S上。若该点不在曲面S上，有两种情况：s>0时称该点相对于曲面S具有正的坐向，反之称该点相对于S具有负的坐向。例如，球面以外的点相对于球面具有正坐向，沿坐标轴正向离开垂直于坐标轴平面的点相对该平面具有正坐向。 任何曲面都把空间分成两部分，一部分相对于该曲面具有正坐向，另一部分具有负坐向。填写栅元卡时，在曲面号之前用正负号表示栅元中的点相对于该曲面的坐向，正号可省略。 b)​ 交和并 (intersection & union) 图2.2 交和并 交集是两个集合的公共部分，如图2.2a所示。并集是两个集合的合集，见图2.2b。MCNP中交算符是隐含的，在栅元卡上两个曲面号之间只有空格就表示交运算。并算符用冒号 : 表示。 图2.3 例子 若想仅用交算符定义栅元，栅元内所有点对特定的界面必须有相同的坐向。如图2.3所示，1、2、3、4、5面围成一个实体，用密度为1.0g/cm3的材料1填充，外部是真空。实体内P1和P2点相对于4面有负坐向，P3点相对于4面有正坐向；P1点相对于3面有正坐向，P2和P3点相对于3面有负坐向。这时最好借助另一平面6把实体分成两部分，从而顺利定义栅元： 1 1 -1.0 1 –2 –3 6 2 1 -1.0 1 –6 –4 5 然而，只用交算符还不便定义实体之外的区域，这时用并算符比较方便，实体外可由3、4面以上区域的交集和其它几个面以外区域的并集构成一个栅元： 3 0 -5 : -1 : 2 : 3 4 当然实体之内的区域也可以借助并算符定义，把1、2栅元合并(这样定义会降低计算速度)。由于交运算优先，所以要用括号把并运算括起来： 1 1 -1.0 1 -2 (-3 : -4) 5 c) 补 (complement) 补算符#的用法有两种： 1. #n 表示当前定义的栅元是栅元n其余的部分。 2. #(...) 表示括号内曲面号定义区域之外的部分。 引入坐向概念后，补运算就不是交运算和并运算之外的新概念，而是交和并的另一种表示方法。例如，对上面第一种情况的处理过程是，去掉#号，把n用括号括上，n中的交算符变成并算符，n中的并算符变成交算符而且加上背靠背的括号) (，曲面之前的坐向取相反符号。数学中常用表示交运算、用+表示并运算、ˉ表示补运算，有 ， ，此处 和a则用坐向指定。 使用补算符时要格外小心，用户常因对补算符理解不透而弄出几何逻辑错误，因此建议不使用补算符。另外，#(n)是不合规定的。 C. 曲面定义中的一些概念 a) 模糊面 (ambiguity surface) 图2.4 模糊面的作用 有时需借助栅元界面之外的曲面来描述栅元。如图2.4所示，假设图中实体绕y轴旋转对称，栅元1和2相对于球面都有正坐向、相对于平行于y轴的柱面都有负坐向。这时，可用垂直于y轴的平面a(y=0)来区分栅元1和2，平面a被称为模糊面。模糊面是为了使栅元描述不模糊而设置的曲面，它不是本栅元界面，但通常是其它栅元的界面，有时需要借助不止一个模糊面。 b) 叶 (sheet) 图2.5 二叶锥面 锥面和椭环面分多叶。通常用方程定义的锥面含有两个叶，相对于锥面旋转对称轴所平行的坐标轴而言，一个具有正的展开方向，另一个具有负向斜率。一个使用二叶锥面描述的栅元常要借助模糊面来区分栅元在锥面的哪一个叶。为了减少使用模糊面，可以在曲面卡上只定义锥面的一个叶，用+1表示正斜率的叶，-1表示负斜率的叶。这种描述方法仅适用于旋转对称轴平行于坐标轴的锥面。 D. 思考题 1. MCNP用什么方法描述几何。 2. 如何用曲面描述栅元几何。 3. 模糊面有什么作用。 §2.2 几何描述卡 在输入文件INP中，栅元卡用于栅元块，曲面卡用于曲面块。材料卡出现在数据块中，将在第四章中介绍它。 A. 栅元卡 (cell card) 格式: j m d geom param j: 在1~5列上填写的栅元标号，可以不连续。MCNP按照读入的顺序对栅元另行编号，称为栅元的程序编号，从1开始单调上升。为避免程序编号和栅元标号混淆，建议按顺序给出栅元标号 m: 该栅元的材料号，是材料卡Mm中相应的序号。真空栅元m=0 d: 栅元材料的密度。填正值表示原子密度，单位atoms/barn-cm。负值表示质量密度，g/cm3。对于真空栅元，该项不填写 geom: 栅元几何说明，列出界定该栅元的曲面的号(有符号，表坐向)，和描述这些曲面包围区域之间关系的布尔算符。布尔算符包括交算符(缺省)、并算符( : )和补算符(#)。缺省运算顺序是，先补，再交，最后并。可用括号改变运算次序 param: 可选的栅元参数说明，形式为：关键字=某一值 在栅元卡上可以定义栅元参数，以替代在输入块中定义栅元参数。这里允许的关键字是：带有粒子标识符的IMP、VOL、PWT、EXT、ECL、WWN、DXC和TMP。栅元卡参数项上的等号可用空格代替。若在栅元卡上指定了某种栅元参数，则在数据块中不能重复指定。 B. 曲面卡 (surface card) a) 由方程定义曲面 格式: j n list j: 在1~5列上填写的曲面标号，可以不连续 n: 方程助记名 list: 按规定次序给出的方程数据项 下表列出了MCNP识别的曲面的类型、方程、助记名和数据项次序： 助记名 类型 说明 方程 数据项 P 平面 一般c) Ax+By+Cz-D=0 A B C D PX 垂直于x轴 x-D=0 D PY 垂直于y轴 y-D=0 D PZ 垂直于z轴 z-D=0 D SO 球 球心在原点 x2+y2+z2-R2=0 R S 一般 (x-x0)2+(y-y0)2+(z-z0)2-R2=0 x0 y0 z0 R SX 球心在x轴 (x-x0)2+y2+z2-R2=0 x0 R SY 球心在y轴 x2+(y-y0)2+z2-R2=0 y0 R SZ 球心在z轴 x2+y2+(z-z0)2-R2=0 z0 R C/X 圆柱 平行于x轴 (y-y0)2+(z-z0)2-R2=0 y0 z0 R C/Y 平行于y轴 (x-x0)2+(z-z0)2-R2=0 x0 z0 R C/Z 平行于z轴 (x-x0)2+(y-y0)2-R2=0 x0 y0 R CX 轴心为x轴 y2+z2-R2=0 R CY 轴心为y轴 x2+z2-R2=0 R CZ 轴心为z轴 x2+y2-R2=0 R K/X 锥 平行于x轴 [(y-y0)2+(z-z0)2]1/2-t(x-x0)=0 x0 y0 z0 t2 1 K/Y 平行于y轴 [(x-x0)2+(z-z0)2]1/2-t(y-y0)=0 x0 y0 z0 t2 1 K/Z 平行于z轴 [(x-x0)2+(y-y0)2]1/2-t(z-z0)=0 x0 y0 z0 t2 1 KX 轴心为x轴 [y2+z2]1/2-t(x-x0)=0 x0 t2 1 KY 轴心为y轴 [x2+z2]1/2-t(y-y0)=0 y0 t2 1 KZ 轴心为z轴 [x2+y2]1/2-t(z-z0)=0 z0 t2 1 SQ 椭圆面 双曲面 抛物面 轴心平行于x、y或z轴 A(x-x0)2+B(y-y0)2+C(z-z0)2 +2D(x-x0)+2E(y-y0)+2F(z-z0) +G=0 A B C D E F G x0 y0 z0 GQ 柱面、锥面椭圆面、双曲面或抛物面， 轴心不平行于坐标轴 Ax2+By2+Cz2 +Dxy+Eyz+Fzx +Gx+Hy+Jz+K=0 A B C D E F G H J K TX 平行于坐标轴的椭圆或环形 (x-x0)2/B2+{[(y-y0)2+(z-z0)2]1/2-A}2/C2-1=0 x0 y0 z0 A B C TY (y-y0) 2/B2+{[(x-x0)2+(z-z0)2]1/2-A}2/C2-1=0 x0 y0 z0 A B C TZ (z-z0) 2/B2+{[(x-x0)2+(y-y0)2]1/2-A}2/C2-1=0 x0 y0 z0 A B C b) 用点定义对称曲面 格式: j n list j: 在1~5列上填写的曲面标号，可以不连续 n: =X、Y或Z， list: 一至三个点的坐标，每个点用一对值表示 类型为X、Y或Z的曲面卡用坐标点而不是用方程来分别描述关于X、Y或Z轴对称旋转的曲面。如果这个曲面由一叶以上组成，指定的点必须都在一叶上。每对坐标值定义曲面上的一个点，例如，在Y卡上可以给出 j Y y1 r1 y2 r2 其中(yi, ri=(xi2+zi2)1/2)是i点的坐标。给出点数不同，曲面类型也不同：i) 一对坐标定义一个平面(PX、PY、PZ)；ii) 两对坐标定义一个线性曲面(PX、PY、PZ、CX、CY、CZ、KX、KY、KZ)；iii) 三对坐标定义一个二次曲面(除以上类型还有SO、SX、SY、SZ和SQ)。 c) 用点定义一般平面 格式: j P x1 y1 z1 x2 y2 z2 x3 y3 z3 j: 在1~5列上填写的曲面标号，可以不连续 x1 y1 z1: 定义该平面的坐标点 MCNP对用户指定的P型曲面，将检查数据项个数，若是4项，则作一般平面方程Ax+By+Cz-D=0的系数理解；若是9项，则作为平面上三个点的坐标值理解，所产生的平面方程系数遵循以下原则：1) 应使坐标原点相对于该平面具有负坐向；2) 当平面通过原点(D=0)，应使(0,0, )相对该平面具有正坐向；3) 若以上两项无法满足(D=C=0)，应使(0, ,0)相对该平面具有正坐向；4) 若以上三项无法满足(D=B=C=0)，应使( ,0,0)相对该平面具有正坐向；5)若以上四项都无法满足，那么用户指定的三个点一定处于一条直线上，则打印警告信息。 C. 思考题 1. 栅元卡的格式如何。为什么最好按顺序指定栅元标号。卡上材料号对应什么，密度数据项填负值表示什么。 2. 如何用曲面卡定义曲面。 §2.3 有效地构建几何 即使是十分有经验且细心的用户，在构建比较复杂的几何时，也经常会犯一些错误。有些错误程序可以检查出来，有些则不能，以致造成错误结果。几何构建还直接影响程序运行效率，不合适的几何构建将浪费机时。 A.​ 定义栅元的原则 原则1：栅元界面要少 图2.6 栅元界面示意 尽管可用大量的并算符定义一个几何上比较复杂的栅元，但这样做并不明智。问题在于，每次计算栅元内碰撞之间的径迹时，都计算径迹和栅元所有界面是否相交以判断粒子是否进入其它栅元，太多的栅元界面会浪费机时。如，图2.6a中的几何只是一些平行的柱面，较容易定义，然而两个大柱面之间的栅元有14个界面(包括前后界面)；图2.6b则是比较有效的栅元定义，复杂界面的栅元被分隔成一些小栅元，每个栅元的界面有所减少。 原则2：栅元大小要合适 构建物理空间时，除按不同材料区域分割栅元外，还必须考虑效率因素。为了节省运算时间，栅元不能定义得太大或太小。在比较重要的区域，建议调整栅元大小，使粒子进入相邻栅元后数量约减少一半，以利于使用几何分裂和轮盘赌技巧。 原则3：尽量用简单曲面定义栅元 如前所述，MCNP频繁计算碰撞点到界面距离。点到二阶以上曲面的距离一般不能用解析法求解，用牛顿下山法求解时，多次叠代很浪费时间。如果问题条件允许，尽量使用二阶以下曲面分段替代二阶以上曲面。 原则4：避免曲面和曲面相切 即使使用双精度变量，也不能计算出点到曲面的确切距离，因而在曲面和曲面相切处有时会因计算精度不够导致栅元重叠或空缺，造成粒子丢失。虽然这种粒子丢失几率很小，一般不会影响记数结果，但计算大量粒子问题时可能会因为丢失粒子数量超过限制，使程序停止运行。 B.​ 检查几何错误 MCNP检查INP文件的数据时，不能检查出几何逻辑错误。而运行中粒子丢失时，才能检查出栅元的重叠或栅元之间的空缺。因此在正式运算前，应采取以下办法检查几何逻辑错误： 1)​ 增加一个VOID卡。这个卡将废弃这个问题的材料说明，将全部栅元设为真空，并把加热记数转化为通量记数。 2)​ 对这个问题增加一个大球面，这个球面包围需计算的系统，球面把系统之外的栅元分成两个栅元。球内面是现在问题的边界，内部所有栅元的重要性都置为1。 3)​ 源描述改为：SDEF SUR=m NRM=-1，其中m是刚才定义大球面的标号。 由于上述改动使粒子在计算系统中没有碰撞，一个短时间运行就会使大量粒子的轨迹穿过系统，可能在出现几何错误的地方引起粒子丢失。 粒子丢失时，OUTP文件上将产生丢失粒子打印表。表中列出丢失粒子跨越的所有曲面，并告诉你粒子在什么位置上向什么方向运行后丢失。你可以从这些信息中推断出粒子丢失的原因。 C. 习题 图2.7 井况模型 1. 请描述如图2.7所示井况模型，为高1m、半径70cm的圆柱体。地层半径为10cm至70cm，材料号为1，密度2.65g/cm3。水泥环厚30mm，材料号为2，密度1.95g/cm3。套管外径140mm，厚7.72mm，材料号为3，密度7.86g/cm3。井内介质材料号为4。为缩短计算时间，地层栅元的大小取高5cm、环距7.5cm，其它栅元高度也取5cm。本题中共有220个栅元。 第三章 数理基础 本章介绍粒子在物质中输运的处理，记数的设置，以及减小方差的技巧。 §3.1 物理 A. 粒子权重 (particle weight) 若严格模拟实际过程，则一个MCNP粒子仅代表一个实际粒子，且具有为1的权重。然而为了提高计算效率，MCNP允许用户使用许多并非实际过程的偏倚技巧。偏倚就是使一些感兴趣区域(包括位置、时间、能量等方面的区域)中粒子数增加，另一些区域中粒子数减少。为使计算结果不出现偏差，人为地改变粒子数后，相应地要调整每个粒子的权重。所谓粒子权重，就是每个模拟粒子携带的一个数w，它代表该粒子对最终结果有w份贡献，或者说该粒子相当于w个实际粒子。w不一定是整数，有时小于1。 不论模拟多少个源粒子，MCNP都给出相对于一个源粒子的记数。如果用户希望得到一定源强条件下的结果，可在SDEF卡上改变源粒子权重，指定WGT=w，这时记数结果是相对于发射w个源粒子。 B. 粒子径迹 (particle tracks) 径迹是源粒子和诱发粒子在其历程中的运行轨迹，粒子在两次与原子核碰撞之间的运行线路是直线，粒子自身之间的作用忽略不计。在任何给定材料的栅元中，MCNP计算沿径迹方向到下次碰撞点的距离而不是计算单位体积内碰撞几率。假设微观总截面是t，要抽样的径迹长度是l，可以证明，l和抽取出的[0,1)区间随机数之间的关系是 l=-ln/t (3.1) C. 中子与物质作用 中子不带电，它与物质中原子的相互作用主要表现为与原子核作用，与核外电子的作用可忽略。不同能量的中子与原子核的作用有不同特点：1)能量在0~1keV之间的慢中子与原子核主要发生弹性散射和俘获反应。热中子能量kT=0.025eV(k是玻尔兹曼常数，T是绝对温度，取293K)，它是和原子热运动达到平衡时的能量，其速度分布接近麦克斯韦分布。能量0.5eV的中子称为超热中子。能量在1eV至1keV之间的中子称为共振中子，这种中子被原子核强烈地共振吸收，吸收截面很大。2)能量在1keV~0.5MeV之间的中能中子与原子核的主要作用形式是弹性散射。3) 快中子的能量在0.5MeV~10MeV之间，与原子核主要发射弹性散射和非弹性散射反应。4)能量高于10MeV的特快中子与原子核发生弹性散射、非弹散射反应等。MCNP3B模拟中子的能量范围是10-11~20MeV。 当能量很高的中子与大块物质作用时，经过多次与原子核碰撞，能量不断减小，最后或者被俘获，或者变成热中子。低能中子和原子核碰撞时受热运动的影响，也受碰撞原子核周围原子的影响。MCNP一般用自由气体模型处理原子核热运动，这种模型下弹性散射截面在零温度截面基础上作修正；其它中子反应的截面不随温度变化。对于某些常用原子，MCNP还提供一些温度条件下的考虑了晶格和化学键影响的精细S(,)表。如果中子能量足够低，使用S(,)表后，总截面是常规截面表中俘获截面以及S(,)表中弹性截面、非弹性截面的总和。否则，总截面从常规截面表中得出，随温度变化作调节。 MCNP模拟中子输运时，要跟踪径迹反复考虑如下事项：a) 碰撞前中子的能量、飞行方向和径迹长度；b)与哪种核碰撞，被碰撞核的速度与方向(除非使用S(,)处理)；c) 碰撞过程是哪种反应，如弹性散射、(n,xn)、(n,p)、(n,)、俘获等；c) 碰撞过程选择地产生光子，模拟每个光子的初始能量、运动方向、权重；d) 碰撞后中子(或几个中子)的能量、飞行方向、权重，下一次在什么位置和哪种核碰撞。如此反复，直至中子能量足够低或权重足够小。 1) 选择碰撞核 设介质由n种核组成，是一个[0,1)上的随机数，当 (3.2) 成立时，第k种核被选择为碰撞核，ti是栅元中第i种核的微观总截面。 2) 靶核速度抽样 如果不使用S(,)处理，需要考虑靶核速度。对轻核，若中子能量高于400kT(10eV)，对重核高于40kT(1eV)，靶核的速度对碰撞截面的影响可忽略不计。其它情况下，要按麦克斯韦分布抽取靶核速度， (3.3) 其中=(AMn/2kT)1/2，A是靶核质量，Mn是以MeV-sh2/cm2为单位的中子质量。靶核运动方向从各向同性分布中抽取。然后得到动能为E的中子的等效散射截面 (3.4) 其中vrel是中子和靶核间的相对速度，s(vrel)为相对速度的散射截面(从截面数据库中查出，未作温度修正)，是中子入射和出射方向夹角的余弦。相对速度的计算是， (3.5) 其中vn是中子速度，V是靶核速度，T是二者速度夹角的余弦。 3) 弹性散射和非弹性散射 如果碰撞时中子未被俘获，它就参与弹性散射或非弹散射反应，广义的非弹散射可以是(n,n’)、(n,np)、(n,2n)等。弹性散射以 的概率被选择，非弹散射以余下的概率被选择，其中el是弹性散射截面, in是非弹散射截面。 出射粒子方向：对弹性和非弹散射，抽取出射粒子飞行方向的方法是一样的，入射和出射方向夹角的余弦从截面文件的角分布表中抽取。对每种特定的中子入射能量，角分布表都有32个分立的箱。如果入射能量是E，介于有角分布表的能量En和En+1之间，则以(E-En)/(En+1-En)的几率从n表中抽取角分布，以(En+1-E)/(En+1-En)的几率从表n+1中抽取角分布。抽取随机数，找到第i个余弦箱，使得i-132