基于组合几何与自定义网格的粒子输运模拟方法

张玲玉，李 瑞，李 刚，付元光，邓 力

(1. 北京应用物理与计算数学研究所；2. 中国工程物理研究院 高性能数值模拟软件中心： 北京 100088)

蒙特卡罗粒子输运方法能较逼真地描述真实的物理实验过程[1]。在一定意义上，它可部分替代物理实验。因此，该方法应用的领域日趋广泛，成为解决实验核物理、反应堆物理、高能物理等实际问题非常有效的工具。传统的蒙特卡罗粒子输运方法通过组合几何布尔运算解决3维复杂几何的描述和建模难题，通过几何实体中粒子的定位和求交运算实现粒子追踪模拟[2]，优点是几何仿真能力强、物理建模完善和模拟近似程度少。目前，常用的蒙特卡罗粒子输运软件大多采用组合几何方法，如美国洛斯阿拉莫斯国家实验室研制的MCNP程序[3]，其几何块由1阶、2阶及包括某些特殊的4阶表面所包围，每个几何单元可通过包围表面的交、和及余来定义。北京应用物理与计算数学研究所与中国工程物理研究院高性能数值模拟软件中心研制的JMCT程序[4]，具备CAD可视化输入输出界面[5]，通过球、长方体和圆柱体等基本体素的布尔运算得到复杂的装置模型，并可通过编辑、显示、分析和定制等功能对模型进行各种操作。

但在某些复杂几何复杂输运、大尺度场景长程输运和高精度高分辨率统计计数应用中，仅使用组合几何的建模及粒子输运方式不能满足需求。如，在复杂内爆压缩客体辐射照相全过程模拟中，内爆压缩客体模型十分复杂[6]，采用组合几何方法建模不仅繁琐，易出错，且频繁地求交和定位运算，使计算效率较低；在大气场景辐射输运模拟中[7]，为描述大气密度随高度的渐近变化，需要将大气划分为多个几何体，并设置相应的材料密度，采用组合几何的建模方式描述大气几何体将不够灵活；在辐射线缆系统电磁脉冲光电输运模拟中[8]，采用组合几何建模方法，只能得到交界面上的电子流分布，为得到线缆上精细的电荷分布，需将线缆划分为具有更小体积的实体，此时传统的组合几何建模方法会影响粒子真实的游走过程，计算结果会产生较大偏差。

针对上述应用特点，如采用组合几何与自定义网格相结合的方法进行建模及输运模拟，对模型简单且相对固定的部分，采用组合几何建模及输运模拟，具有简单、便捷和高效的优势；对于模型几何结构、材料、计数复杂的部分，采用自定义网格的方法建模及输运模拟，具有建模灵活、计算精度高的优势。因此，本文基于自主研发的JMCT软件，开展了基于组合几何与自定义网格的粒子输运模拟研究。

1 基于自定义网格的粒子输运算法

1.1 基于自定义网格的自动建模方法

复杂几何一般具有结构复杂及材料繁多的特点，使用几何实体建模易出错，效率低。另外，有些几何模型来自流体动力学计算的结构网格或非结构网格[9]，如蒙特卡罗输运模拟采用与流体力学相同的网格，则计算效率非常低。本文采用基准网格与细化网格相结合的建模方法，如图1所示。

基准网格：首先建立一个包围几何体的网格盒子，再在盒子内进行网格划分，网格划分可均匀划分，也可局部加密划分，进而产生自定义基准网格，如图1蓝色网格部分所示。在每个小网格内，采用中心点网格材料模型[10]确定网格的材料和密度，即用网格中心点的材料和密度作为整个小网格的材料和密度，有效避免材料混合引起的材料数目繁多的缺点。中心点网格材料模型是一种简化处理方式，能减少程序前处理的时间和模拟计算需要的内存，但必须进行质量守恒检查，保证模型的合理性。

细化网格：在模型重要区域或复杂区域将网格进行细化，每个细化网格与基准网格构造方法相同，如图1红色网格部分所示。细化网格可与基准网格或实体几何重叠，后续输运中根据几何层叠的优先级触发输运事件。细化网格与基准网格构成自定义网格。

将组合几何模型与自定义网格模型结合，形成粒子输运可计算模型。其中，组合几何模型主要描述规则结构几何体，如四面体和柱体等；自定义网格主要描述几何结构及材料复杂的几何体。该模型的优点是可精确建立复杂几何模型，并可与其他程序建立的复杂模型对接。

目前JMCT支持3维笛卡尔坐标系和3维柱坐标系的自定义网格。2种坐标系下，要求网格是矩形，大小可不均匀；且JMCT支持多套层叠自定义网格；为检测自定义网格模型的正确性，JMCT还支持目标区域网格材料的可视化输出，用于检验及协助建模。

1.2 基于自定义网格的粒子输运方法

粒子输运包括3个关键技术，分别是粒子定位、粒子求交和粒子追踪。

1.2.1 粒子定位

用(x,y,z,u,v,w,E,W)表示粒子当前的状态，x,y,z,u,v,w,E,W分别表示粒子的x方向坐标、y方向坐标、z方向坐标、x方向速度、y方向速度、z方向速度，粒子能量和权重。单个结构网格上的粒子定位就是确定粒子x,y,z在3维网格的索引(i1,i2,i3)，由于3个维度网格坐标均为单调递增形式，在网格内部采用二分查找即可实现。

目前由于采用了光线追踪粒子输运算法[11]，在每个网格边界上会强制触发事件，粒子会在边界上停留。为避免粒子定位在网格边界出现定位死锁的问题，需保证粒子飞行时正常跨越网格边界。因此，在网格定位时，根据粒子的飞行方向，对粒子坐标沿飞行方向进行了微小的调整，可表示为

(1)

其中：l′为微小的位移，取1×10-10cm，x′,y′,z′分别为边界处事件发生时粒子的位置。

由于蒙特卡罗方法的随机性，出现粒子飞行方向与网格相切情况的概率极小。此时，仅利用位置平移无法解决定位死锁问题。因此，在连续平移数次之后，对粒子飞行方向也进行了微小偏转调整，进而保证粒子飞行时正常跨越网格边界，旋转后粒子的坐标可表示为

(2)

其中，θ为微小的旋转角度，取0.001 rad。

1.2.2 粒子求交

粒子求交的目的是计算粒子到边界的距离。当粒子处于自定义网格内部时，根据粒子飞行方向，计算粒子特征线与所在网格边界面的交点和距离，进而得到粒子到网格边界的距离。计算点到面距离的本质是联立求解粒子特征线方程与面方程。对于直角坐标系，根据粒子位置(x,y,z)与飞行方向网格位置(xboud,yboud,zboud)可得到3个方向的点面距离lx,ly,lz，可表示为

(3)

其中，lx,ly,lz的最小值便是粒子到当前网格边界的距离。

对于柱坐标，求解特征线与r=Rmin,r=Rmax曲面、θ=θmin,θ=θmax平面及z=zmin,z=zmax平面的点面距离，可表示为

(4)

其中，lr，min,lr，max,lθ，min,lθ，max,lz，min,lz，max的最小值便是粒子到当前网格边界的距离。

1.2.3 粒子追踪

基于组合几何与自定义网格模型的粒子追踪流程如图2所示。

基于组合几何与自定义网格模型的粒子追踪，着重要判断粒子是否在网格内，解决网格与网格之间、网格与组合几何之间的重叠问题。本文在每个自定义网格外都有一个包围盒子，采用安全距离提高输运效率。根据几何层叠的优先级，设计输运事件响应机制，避免几何重叠并提高输运效率，目前JMCT中自定义网格模型优先级高于组合几何模型。

1.3 基于自定义网格的统计计数方法

粒子输运模拟中，需将输运过程中关键的物理信息进行记录，并得到统计平均值，用于后续分析。常见的统计平均计数有基于几何实体的面通量、体通量、点通量及能量沉积等[3]。这些计数的特点是基于几何实体中的几何面、几何点或几何体进行计数。如果要得到局部或全局精细的计数，必须增加几何实体或探测点的数目，提高计数分辨率。然而高分辨率的几何实体建模和储存是比较困难的，且会影响粒子真实的游走过程，使计算效率显著降低。网格计数能很好地解决以上问题，目前常用的蒙特卡罗程序支持基于网格的体通量计数[12]，并被用于反应堆屏蔽设计和探测器响应计算等实际问题。但受限于实体几何和面几何的建模分辨率，对能量沉积计数和电荷统计计数的计算效率和计算精度较低。因此，本文也研制了基于自定义网格的统计计数方法。

1.3.1 基于自定义网格的能量沉积计数

一个中子或光子从产生到消失，能量沉积可表示为

ΔE=E0-∑Es+Ecut

(5)

其中：E0为源粒子的能量；Es为次级粒子的能量；Ecut为中子或光子的截断能量。

电子受库仑力影响，一个电子从产生到消失，能量沉积可表示为

(6)

用自定义网格罩住计算模型中需统计能量沉积的区域，粒子仍在原来模型的几何体中输运。当粒子产生能量沉积时，判断粒子所在网格，并将能量沉积计入相应网格内。

每一次碰撞j更新一次统计量，对任意网格k，能量沉积均值及方差为

(7)

其中：N为模拟的粒子数；n为某一个粒子在第k个网格内的碰撞次数；ΔEj(i)为第i个粒子在第j次碰撞中沉积在第k个网格上的能量。网格的定位与粒子输运中的定位方法相同。

1.3.2 基于自定义网格的电荷统计计数

网格电荷统计计数主要是统计输出自定义网格上的电子数目、正电荷数目和净电荷数目。传统的蒙特卡罗模拟方法可通过电子面流量计数，推算出几何体的电荷计数。针对计算电荷计数精细分布的需求，需建立多个几何实体来获得分界面，这种计数方式计算量非常大，效率很低。

通过研究光电耦合输运过程发现，只有在电子产生及历史结束的位置才会引起电子及正电荷数目的变化。因此，网格电荷统计计数原理是：电子产生处，正电荷和电子数都增加1；电子死亡处，电子数加1，相应的电子产生处的电子数减1。统计方法不关心电子输运过程，只关心电子产生及死亡的位置和方式。任意网格k上的正电荷数和负电荷数可表示为

(8)

其中：Bk,i为第i个粒子在第k个网格内产生的电子数；Dl,i为第i个粒子在第l个网格内历史结束的电子数；nc为网格数；Dk,j是第j个粒子在第k个网格内历史结束的电子数。网格的定位与粒子输运中的定位方法相同。

2 实际应用

2.1 线缆系统电磁脉冲光电输运模拟

线缆系统电磁脉冲是电磁脉冲的一个分类，是指X或γ射线与线缆相互作用产生的电磁脉冲。飞行器表面材料及线缆受X射线照射后产生光电子，形成光电流及线缆电流，以干扰的形式进入到与线缆相连接的电路中，对电子系统产生干扰或毁伤。线缆材料受X射线照射产生的电荷分布，是后续系统电磁脉冲研究的基础。

本文针对典型线缆模型，计算了线缆上的精细电荷分布。图3为典型线缆模型示意图。模型由3层构成，从内到外分别为铜制芯线、聚乙烯介质层和铝制屏蔽层，外径依次为0.01，0.05，0.07 cm。电缆长度为1 cm。能量为10 keV的X射线沿线缆轴线方向均匀辐照圆柱侧面。

采用电子面流量计数及自定义网格电荷计数分别进行统计。表1给出了采用JMCT与MCNP计算芯线外边界和介质外边界的电子面电流计数的计算结果。由表1可知，JMCT的模拟结果与MCNP的模拟结果相对偏差小于5%，结果可信。

表1 电子面流量计数Tab.1 Electron surface flow tally

为验证自定义网格电荷统计的正确性，本文采用自定义柱坐标网格统计线缆上的电荷分布，网格与实体几何模型重合。由于电子面流量计数无法直接统计实体中的净电荷，本文采用电子流出与流入的差值进行计算，线缆上的静电荷计数如表2所列。

表2 线缆上的静电荷计数Tab.2 Net charge tally

由表2可知，自定义网格净电荷计数结果与电子面流量计数结果自洽，与文献[8]模拟结果相对偏差小于5%，表明自定义网格电荷统计计数方法可行，模拟结果可信。

本算例采用基于组合几何与自定义网格的建模方法、粒子输运方法和统计计数方法，对线缆结构采用组合几何方法建立可计算模型，采用自定义网格电荷分布统计方法得到线缆上的精细电荷分布。该方法充分发挥了组合几何建模简单和自定义网格计数灵活的优势。

2.2 高能闪光照相模拟

高能闪光照相是利用高能X射线对客体进行投射照相[13]，研究透射辐射形成的图像来推断客体性质(边界和密度)的技术。高能闪光照相是诊断流体动力学实验的最主要手段之一，可用于武器内爆压缩物理特性和几何结构的诊断及物理规律的研究。

闪光照相系统几何模型包括2大类。(1)辐射照相设备模型。包括轫致辐射靶、准直孔、网栅、转换屏和探测器等，特点是模型简单且相对固定。针对这类模型，可使用可视化组合几何建模方法，生成相对固定的可计算模型，模型建好后可长期使用。(2) 内爆压缩客体模型。这类模型大多来自流体动力学计算中结构网格或非结构网格产生的客体模型，特点是几何结构和材料复杂，需与流体动力学计算结果对接，变化较多。针对这类模型，使用自定义网格自动建模方法，生成适合蒙特卡罗辐射输运的计算模型。

本文针对典型闪光照相系统模型进行了基于组合几何与自定义网格输运功能的测试。辐射照相设备包含高密度准直孔、实验客体、Al板、胶片接收系统和将X光转换为可见光的转换屏。闪光照相客体是典型的法国实验客体(FTO)，由3层结构组成：半径为1 cm的真空层、半径为4.5 cm的铀层和半径为6.5 cm的铜层。X光源为4 MeV的单能X射线，位于圆点，方向垂直向上。图4为典型闪光照相系统照相模型。采用2种方式模拟：(1) 基于组合几何的输运模拟，照相设备和辐射客体均采用组合几何建模及输运模拟，如图 4(a)所示；(2) 基于组合几何与自定义网格的输运模拟，照相设备采用可视化组合几何建模，如图4(b)所示。辐射客体采用自定义网格建模，1 000×1 000的网格，如图 4(c)所示。

Al板前的直穿照射量随探测器在x轴的位置的变化关系，如图 5所示。图5中，黑色曲线为组合几何模型计算的直穿照射量，红色曲线为组合几何与自定义网格模型计算的直穿照射量。由图5可见：曲线分为3段，分别对应FTO的3层结构，模拟结果符合物理规律；基于组合几何与自定义网格模拟的FTO直穿照射量模拟结果与基于组合几何的模拟结果吻合较好，最大相对偏差小于5%。模拟结果表明，基于组合几何与自定义网格粒子输运方法可行，模拟结果可信。

本算例采用了基于组合几何与自定义网格的建模方法、粒子输运方法和统计计数方法，对辐射照相设备部分采用组合几何方法建立可计算模型，对于内爆压缩客体部分采用自定义网格自动建模方法建立可计算模型。该方法充分发挥了组合几何建模便携、计算效率高，自定义网格建模灵活和计算精度高的优势。

3 结论

在复杂几何复杂输运、大尺度场景长程输运、高精度高分辨率统计计数应用中，仅使用组合几何的建模及粒子输运方式，存在建模过程繁琐、易出错和计算效率低的问题。本文研究了组合几何与自定义网格相结合的建模、粒子输运及统计计数方法：对模型简单且相对固定的部分，采用组合几何建模及输运模拟，具有简单，便捷和高效的优势；对几何形状、材料和计数复杂的部分，采用自定义网格的方法建模及输运模拟，具有建模灵活和计算精度高的优势。本文将该方法应用到了线缆系统电磁脉冲光电输运模拟和高能闪光照相直穿照射量模拟中。模拟结果表明，本文方法计算结果与参考解相对偏差小于5%，模拟方法可行，结果可信。基于组合几何与自定义网格的粒子输运方法可有效提高复杂几何的计算精度和计数分辨率。