CN104317772B - 一种基于空间网格分割的蒙特卡罗粒子输运快速几何处理方法 - Google Patents

Abstract

本发明公开了一种基于空间网格分割的蒙特卡罗粒子输运快速几何处理方法，基于空间网格分割的思想以及结合蒙特卡罗粒子输运几何处理的特点，提出了一种基于代价函数的空间网格分割方法，根据最小代价原则，依次沿各坐标轴对三维几何模型进行空间网格分割，并确定代价最小时各网格包含的几何体，最终得到的空间网格模型即为最优模型。采用该模型进行蒙特卡罗粒子输运几何处理时，可根据粒子空间坐标快速进行空间网格定位，然后再遍历该网格所包含的几何体，可极大程度地减少查询的候选几何体的数目，快速确定粒子所在几何体，从而加速蒙特卡罗粒子输运模拟计算，提高蒙特卡罗粒子输运的几何处理效率。

Description

一种基于空间网格分割的蒙特卡罗粒子输运快速几何处理方法

技术领域

本发明涉及一种基于空间网格分割的蒙特卡罗粒子输运快速几何处理方法，属于核科学中蒙特卡罗粒子输运计算与数值模拟研究领域。

背景技术

蒙特卡罗粒子输运方法是以统计学理论为基础，通过对大量具有随机性质特点的粒子物理事件和物理过程进行随机模拟来求解粒子输运方程的一种非确定论方法，它广泛应用反应堆物理、医学物理、高能物理和核探测等领域。相比传统的确定论方法而言，蒙特卡罗粒子输运方法的最大优点在于，它能够非常逼真地描述具有随机性质的事物特点及过程，并且对几何模型和材料限制小，可以精确模拟复杂源、复杂材料和复杂几何下的粒子输运问题。然而，它的缺点也非常明显，其主要在于计算非常耗时，收敛速度较慢。因此，如何提高蒙特卡罗粒子输运方法的计算效率成为该领域近几年发展的重点之一。

据统计，蒙特卡罗粒子输运方法在进行粒子输运模拟计算时，30％～80％的计算时间主要耗费在几何处理上，因此，如何采用高效的几何处理方法是提高蒙特卡罗粒子输运模拟计算效率的重要研究内容之一。在进行蒙特卡罗粒子输运模拟时，粒子每随机游走一步，都需要对其进行空间几何体定位和几何步长计算。几何体定位指根据粒子的三维空间坐标数据，判断粒子在整个几何模型中处于哪一个具体的几何体中；几何步长计算指从当前粒子坐标所在处，沿某一给定方向向量计算粒子到它所处几何体边界的直线距离。因此，由于在蒙特卡罗粒子输运模拟时，频繁地进行粒子空间定位和几何步长计算是导致计算耗时的重要原因。

对传统的蒙特卡罗粒子输运几何处理方法而言，在进行粒子空间定位时，通常需要遍历整个模型中的所有几何体，其处理的时间复杂度为O(N)，因此对于蒙特卡罗粒子输运几何处理大规模复杂几何模型时，如核反应堆百万级Pin-by-Pin模型问题，其几何处理效率将会非常低下，从而严重制约蒙特卡罗粒子输运计算的效率，也成为蒙特卡罗粒子输运计算方法走向工程实践的瓶颈之一。

发明内容

本发明要解决的技术问题为：本发明的目的在于改善传统蒙特卡罗粒子输运计算方法在几何处理效率上的不足，提供了一种基于空间网格分割的蒙特卡罗粒子输运快速几何处理方法，从而提高蒙特卡罗粒子输运计算的几何处理效率。

本发明解决上述技术问题采用的技术方案为：一种基于空间网格分割的蒙特卡罗粒子输运快速几何处理方法，包括以下步骤：

步骤(1)、计算几何模型中各个几何体的轴向包围盒：

分析待处理的几何模型，根据计算机图形学中的轴向包围盒AABB技术，为几何模型中的各个几何体构建轴向包围盒，轴向包围盒为各边平行于笛卡尔直角坐标轴的平行六面体，每个包围盒有且完全包含一个几何体，因此描述一个几何体的轴向包围盒只需六个标量即可，即一个几何体的轴向包围盒可表示为：

AABB_J＝[(X_J,min,X_J,max)；(Y_J,min,Y_J,max)；(Z_J,min,Z_J,max)](J＝0～N)

其中J表示几何模型中几何体的编号，N表示整个模型中几何体的总数，AABB₀表示整个模型的轴向包围盒，AABB_J表示模型中第J个几何体的轴向包围盒，(X_J,min,X_J,max)表示几何体J沿X轴的尺寸范围从坐标值X_J,min到坐标值X_J,max，(Y_J,min,Y_J,max)表示几何体J沿Y轴的尺寸范围从坐标值Y_J,min到坐标值Y_J,max，(Z_J,min,Z_J,max)表示几何体J沿Z轴的尺寸范围从坐标值Z_J,min到坐标值Z_J,max；

步骤(2)、基于代价函数最小原则，对给定几何模型进行空间网格分割，其步骤如下：

第一，根据第(1)步得到的整个几何模型的轴向包围盒AABB₀，先假设对其沿X、Y、Z轴按照一定步长依次进行空间网格分割，其分割步长一般选取沿各轴方向的某一最小几何体的包围盒尺寸作为沿对应轴分割的步长，即由如下公式获得：

$X_{\mathrm{step}}=\underset{J=1~N}{\mathrm{Min}}\left\{{\mathrm{AABB}}_{J,X}\right\}=\underset{J=1~N}{\mathrm{Min}}\left\{\right|X_{J,\max }-X_{J,\min }\left|\right\},$

$Y_{\mathrm{step}}=\underset{J=1~N}{\mathrm{Min}}\left\{{\mathrm{AABB}}_{J,Y}\right\}=\underset{J=1~N}{\mathrm{Min}}\left\{\right|Y_{J,\max }-Y_{J,\min }\left|\right\},$

$Z_{\mathrm{step}}=\underset{J=1~N}{\mathrm{Min}}\left\{{\mathrm{AABB}}_{J,Z}\right\}=\underset{J=1~N}{\mathrm{Min}}\left\{\right|Z_{J,\max }-Z_{J,\min }\left|\right\},$

其中X_step、Y_step、Z_step为将整个几何模型沿X、Y、Z轴进行空间网格分割的步长。

第二，假设以X_step、Y_step、Z_step为步长依次对整个几何模型沿X、Y、Z轴进行空间网格分割，得到三种空间网格分割方案，这三种空间网格分割方案得到网格数目N_X、N_Y、N_Z分别为：

$N_X=\frac{|X_{0,\max }-X_{0,\min }|}{X_{\mathrm{step}}}$

$N_Y=\frac{|Y_{0,\max }-Y_{0,\min }|}{Y_{\mathrm{step}}}$

$N_Z=\frac{|Z_{0,\max }-Z_{0,\min }|}{Z_{\mathrm{step}}}$

其中，X_0,min和X_0,max为整个模型沿X轴的最小尺寸和最大尺寸坐标值，Y_0,min和Y_0,max为整个模型沿Y轴的最小尺寸和最大尺寸坐标值，Z_0,min和Z_0,max为整个模型沿Z轴的最小尺寸和最大尺寸坐标值；于是三种网格分割方案的代价函数F(X)、F(Y)、F(Z)分别为：

$F\left(X\right)=\frac{\underset{J=1~N_X}{\mathrm{\Σ}}{M}_{J,X}}{\underset{J=1~N_X}{\mathrm{\Σ}}{N}_{J,X}}$

$F\left(Y\right)=\frac{\underset{J=1~N_Y}{\mathrm{\Σ}}{M}_{J,Y}}{\underset{J=1~N_Y}{\mathrm{\Σ}}{N}_{J,Y}}$

$F\left(Z\right)=\frac{\underset{J=1~N_Z}{\mathrm{\Σ}}{M}_{J,Z}}{\underset{J=1~N_Z}{\mathrm{\Σ}}{N}_{J,Z}}$

其中，F(X)为沿X轴分割得到的网格模型的代价函数，F(Y)为沿Y轴分割得到的网格模型的代价函数，F(Z)为沿Z轴分割得到的网格模型的代价函数；M_J,X、M_J,Y、M_J,Z表示分别沿X、Y、Z轴进行空间网格分割得到的网格模型中第J个网格中包含的几何体数目，N_J,X、N_J,Y、N_J,Z表示分别沿X、Y、Z轴进行空间网格分割得到的网格模型中第J个网格是否为空网格，即该网格是否含有几何体，它们的关系如下：

$N_{J,X}=\left\{\begin{array}{cc}1& \mathrm{if}(M_{J,X}>0)\\ 0& \mathrm{if}(M_{J,X}=0)\end{array}\right.$

$N_{J,Y}=\left\{\begin{array}{cc}1& \mathrm{if}(M_{J,Y}>0)\\ 0& \mathrm{if}(M_{J,Y}=0)\end{array}\right.$

$N_{J,Z}=\left\{\begin{array}{cc}1& \mathrm{if}(M_{J,Z}>0)\\ 0& \mathrm{if}(M_{J,Z}=0)\end{array}\right.$

由此可知，代价函数的意义在于：衡量一个给定的空间网格分割模型，平均每个网格包含的几何体数目；

第三，选择代价函数最小的网格分割为最优的空间网格分割方案：

根据第二步得到的三种空间网格分割方案，选取代价函数最小的Min{F(X),F(Y),F(Z)}方案作为第(2)步最终的空间网格分割方案；

步骤(3)、针对第(2)步得到的空间网格分割得到的每个网格，再重复第(2)步骤，并依次递归下去，直到待分割的网格中包含的几何体数目少于一定数目时，停止分割；第(3)步和第(2)步最大的区别在于：第(2)步以整个几何模型为分割对象，而第(3)步则针对某一个特定网格继续进行空间细分，但其本质是一样的；

步骤(4)、基于空间网格分割得到的网格模型进行快速几何体查询：

第一步，根据粒子坐标，确定包含该点的网格编号；在进行蒙特卡罗粒子输运模拟时，在给定粒子三维空间坐标数据P(X_p,Y_p,Z_p)时，首先根据以下公式，快速进行网格定位，查找该点所在的网格编号：

其中C_i为点P(X_p,Y_p,Z_p)在坐标轴i方向上的分量，L_i为整个模型或某一网格沿坐标轴i的最小边界坐标值，|C_i-L_i|为给定点沿坐标轴i方向到模型边界的距离，S_i为整个模型或某一网格沿坐标轴i进行网格分割的步长，表示取大于的最小正整数，则[ID(X),ID(Y),ID(Z)]即为点P(X_p,Y_p,Z_p)所在的网格编号；

第二步，假设网格[ID(X),ID(Y),ID(Z)]包含的几何体为G₁,G₂,...,G_M，然后依次对这M个几何体进行点包容性计算，从而便可快速查找到粒子所在几何体。

本发明与现有技术相比的优点在于：在将整个几何模型进行空间网格分割时，提出了代价函数的概念，根据代价函数最小原则，可以使最终分割得到的网格模型为最优网格模型；在基于最优网格模型进行蒙特卡罗粒子输运几何处理时，可以将粒子几何体定位操作的复杂度从Ο(N)降低到Ο(F)，其中F为网格模型的平均代价函数：

$F=\frac{\underset{J=1}{\overset{N_T}{\mathrm{\Σ}}}{N}_J}{N_T}$

其中，N_T为最终划分得到的网格模型中总网格个数，N_J为第J个网格包含的几何体数目，N为整个模型中几何体总数，且满足F<<N条件，极大地降低了候选几何体的数目，从而提高了蒙特卡罗粒子输运几何处理效率。本发明通过基于空间网格分割得到的网格模型，在进行蒙特卡罗粒子输运几何处理时，为了确定粒子所在的几何体，可以极大地极少不必要的几何查找数目，从而快速确定粒子所在几何体，这样不仅可以提高蒙特卡罗粒子输运的几何处理效率，还提高了对大规模几何模型处理的适应性。

附图说明

图1为本发明的几何模型示意图；

图2为本发明几何体的轴向包围盒示意图；

图3为本发明的以X轴为主轴进行网格分割得到的网格模型示意图；

图4为本发明的X轴为主轴、Y轴为次轴进行网格分割得到的网格模型示意图；

图5为本发明基于网格模型进行蒙特卡罗粒子输运几何空间定位的流程示意图；

图6为本发明基于网格模型进行蒙特卡罗粒子输运几何步长计算的流程示意图；

图7为本发明计算实例几何模型示意图。

图中，1为待处理的几何模型，2为几何模型中分布的填充几何体，3为几何体的轴向包围盒，4为网格分割线，5为X型网格，6为Y型网格，7为组件模型中的连接棒，8为组件模型中的水棒，9为组件模型中的燃料棒。

具体实施方式

下面结合附图和具体实施例进一步说明本发明。

一种基于空间网格分割的蒙特卡罗粒子输运快速几何处理方法：

1.计算几何模型中各几何体的轴向包围盒：

包围盒方法在计算机图形学中已经是发展非常成熟的方法，因此本发明可以直接利用国际上发展成熟的三维实体造型软件ACIS来计算几何模型中各几何体的轴向包围盒。

第一步，创建三维实体几何模型。如图1所示，其中1为待处理的几何模型，2为几何模型中分布的填充几何体；根据用户给定的几何模型中各几何体的尺寸数据文件，在ACIS软件中依次创建出每一个几何体的实体模型，将几何模型的尺寸数据文件完整地映射到软件中，创建出一个完整的三维实体几何模型；

第二步，计算几何模型中各几何体的轴向包围盒。利用第一步创建的三维实体几何模型，采用轴向包围盒方法，直接调用ACIS软件的包围盒计算功能，依次为几何模型中的各几何体生成轴向包围盒；如图2所示，其中2为几何模型中分布的填充几何体，3为几何体的轴向包围盒；一个几何体的轴向包围盒是一个各边平行于坐标轴的平行六面体，其中(X_min,X_max)、(Y_min,Y_max)和(Z_min,Z_max)分别表示该几何体在X、Y、Z轴向的几何尺寸范围，最终将所有几何体的包围盒信息输出到外部的包围盒数据文件中。

2.准备网格分割关键参数的数据描述文件：

第一步，几何模型中各几何体的尺寸数据文件。该文件可由用户手工撰写或者采用自动建模工具进行模型转换而获得。

第二步，几何模型中各几何体的轴向包围盒数据文件。可以直接利用第1步中生成的包围盒尺寸数据文件。

第三步，网格分割关键参数文件。该文件主要由用户提供，一个X型网格可填充几何体的最大数目N_X,max，一个Y型网格可填充几何体的最大数目N_Y,max，一个Z型网格可填充几何体的最大数目N_Z,max。

3.确定网格分割的沿各轴的步长，并择次沿X、Y、Z轴进行网格分割：

网格分割基本思想：首先根据代价函数最小值原则，从X、Y、Z轴三个轴的网格分割方案中选择一个轴(如Y轴)为主轴进行网格分割，分割得到的网格为Y型网格；然后再针对每个Y型网格，同样根据代价函数最小值原则，从剩余的X、Z轴两个轴的网格分割方案中选择一个轴(如X轴)为次轴，分割得到的网格为X型网格；最后再针对每一个X型网格，同样根据代价函数最小值原则，对剩下的一个轴(如Z轴)进行网格分割，分割得到的为Z型网格，到此得到的网格模型即为最优网格模型。具体步骤如下：

第一步，从X、Y、Z轴三个坐标轴选择一个轴为主轴进行网格分割。首先，根据几何模型中各几何体的包围盒数据，计算沿X、Y、Z轴进行网格分割时的分割步长，可由如下计算公式获得：

$X_{\mathrm{step}}=\underset{J=1~N}{\mathrm{Min}}\left\{{\mathrm{AABB}}_{J,X}\right\}=\underset{J=1~N}{\mathrm{Min}}\left\{\right|X_{J,\max }-X_{J,\min }\left|\right\},$

$Y_{\mathrm{step}}=\underset{J=1~N}{\mathrm{Min}}\left\{{\mathrm{AABB}}_{J,Y}\right\}=\underset{J=1~N}{\mathrm{Min}}\left\{\right|Y_{J,\max }-Y_{J,\min }\left|\right\},$

$Z_{\mathrm{step}}=\underset{J=1~N}{\mathrm{Min}}\left\{{\mathrm{AABB}}_{J,Z}\right\}=\underset{J=1~N}{\mathrm{Min}}\left\{\right|Z_{J,\max }-Z_{J,\min }\left|\right\},$

其中AABB₀＝[(X_0,min,X_0,max)；(Y_0,min,Y_0,max)；(Z_0,min,Z_0,max)]为待分割几何模型的包围盒尺寸，X_step、Y_step、Z_step为将整个几何模型沿X、Y、Z轴进行网格分割的分割步长。

确定各轴网格分割的分割步长后，可依次按X_step、Y_step、Z_step沿X、Y、Z轴分别进行网格分割，可得到三种分割方案，则每种网格分割的网格总数目N_X、N_Y、N_Z分别为：

$N_X=\frac{|X_{0,\max }-X_{0,\min }|}{X_{\mathrm{step}}}$

$N_Y=\frac{|Y_{0,\max }-Y_{0,\min }|}{Y_{\mathrm{step}}}$

$N_Z=\frac{|Z_{0,\max }-Z_{0,\min }|}{Z_{\mathrm{step}}}$

三种网格分割方案的代价函数F(X)、F(Y)、F(Z)分别为：

$F\left(X\right)=\frac{\underset{J=1~N_X}{\mathrm{\&Sigma;}}{M}_{J,X}}{\underset{J=1~N_X}{\mathrm{\&Sigma;}}{N}_{J,X}}$

$F\left(Y\right)=\frac{\underset{J=1~N_Y}{\mathrm{\&Sigma;}}{M}_{J,Y}}{\underset{J=1~N_Y}{\mathrm{\&Sigma;}}{N}_{J,Y}}$

$F\left(Z\right)=\frac{\underset{J=1~N_Z}{\mathrm{\&Sigma;}}{M}_{J,Z}}{\underset{J=1~N_Z}{\mathrm{\&Sigma;}}{N}_{J,Z}}$

其中M_J,X、M_J,Y、M_J,Z分别表示第J个X、Y、Z型网格中填充几何体的数目，N_J,X、N_J,Y、N_J,Z分别表示第J个X、Y、Z型网格是否为空网格，它们的关系如下：

$N_{J,X}=\left\{\begin{array}{cc}1& \mathrm{if}(M_{J,X}>0)\\ 0& \mathrm{if}(M_{J,X}=0)\end{array}\right.$

$N_{J,Y}=\left\{\begin{array}{cc}1& \mathrm{if}(M_{J,Y}>0)\\ 0& \mathrm{if}(M_{J,Y}=0)\end{array}\right.$

$N_{J,Z}=\left\{\begin{array}{cc}1& \mathrm{if}(M_{J,Z}>0)\\ 0& \mathrm{if}(M_{J,Z}=0)\end{array}\right.$

如图3所示，其中4为网格分割线，5为X型网格。假设代价函数最小的为F(X)，则第一次网格分割便以X轴为主轴进行网格分割，得到的分割模型如图3所示，图3所示为将具有相同几何体的网格归并之后的网格分割模型，总共沿X轴方向分割了7个X型网格。

第二步，针对第一步按照主轴划分的所有X型网格，若其中某一X型网格中的填充几何体数目大于N_X,max，则从剩下的两个Y、Z轴中选一条坐标轴作为次轴，继续进行网格分割，其次轴的筛选思想同筛选主轴的思想一致，即按两个轴进行网格划分时，选取代价函数最小的轴为分割次轴。唯一不同的地方在于，此时计算分割步长时采用的是对应X型网格的尺寸作为边界尺寸进行计算。可以看到，当假设N_X,max＝3时，只有第1个和第7个X型网格需要进行进一步网格分割，其分割次轴均为Y轴，如图4所示为对所有X型网格按照次轴进行网格分割后的模型，其中5为X型网格，6为Y型网格。

第三步，针对第一步和第二步按照主轴和次轴已生成的X型、Y型网格，以最后剩下的一个坐标轴为侧轴继续进行网格分割，若某一网格包含的几何体大于N_X,max、N_Y,max或N_Z,max，则继续按照该轴对应的分割步长进行网格分割，分割方法和原理同上。

4.整个网格分割的基本流程如图5所示，最终输出结果包括以下内容：

整个几何模型的网格分割数据文件。

5.基于网格分割模型，进行快速蒙特卡罗粒子输运几何处理：

在进行蒙特卡罗粒子输运计算时，几何处理主要包括粒子空间位置定位和几何步长计算。基于网格分割模型进行快速蒙特卡罗粒子输运几何处理如图6的主要过程如下：

第一步，粒子空间位置定位。根据粒子坐标，确定包含该点的网格编号；在进行蒙特卡罗粒子输运模拟时，在给定粒子三维空间坐标数据P(X_p,Y_p,Z_p)时，首先根据以下公式，快速进行网格定位，查找该点所在的网格编号：

其中C_i为点P(X_p,Y_p,Z_p)在坐标轴i方向上的分量，|C_i-L_i|为给定点沿坐标轴i方向到模型边界的距离，表示取大于的最小正整数，则[ID(X),ID(Y),ID(Z)]即为点P(X_p,Y_p,Z_p)所在的网格编号；

假设网格[ID(X),ID(Y),ID(Z)]包含的几何体为G₁,G₂,...,G_M，然后依次对这M个几何体进行点包容性计算，从而便可快速查找到粒子所在几何体。

第二步，粒子几何步长计算。在确定粒子空间位置后，便可根据粒子与几何模型中几何体的位置关系，计算粒子的几何步长。由于给模型中所有几何体引入了包围盒，在进行几何步长计算时，若粒子径迹不与某一几何体的包围盒相交，则不必进行与该几何体的射线碰撞检测，从而实现加速粒子几何步长计算的功能。

6.计算实例，包含以下内容：

计算实例如图7所示是一个沸水反应堆BWR8×8单燃料组件问题，其中7为组件模型中的连接棒，8为组件模型中的水棒，9为组件模型中的燃料棒。整个燃料组件在X、Y、Z轴向的尺寸为L×W×H，其中L＝13.86cm，W＝13.86cm，H＝380cm。该组件总共由64根棒填充组成，包括2根水棒(黑色)，8根连接棒(阴影)，54根燃料棒(灰色)，每种棒的外径均为d＝1.23cm，每个燃料元件pin的边长为l＝1.7325cm。在进行蒙特卡罗粒子输运模拟计算时，假定给定一个空间点位置P(7.79625,7.79625,190)，要求基于该组件的空间网格模型实现快速几何处理功能。

如图7所示，以组件模型右下角顶点为笛卡尔XYZ直角坐标系原点O(0,0,0)，并给定以下网格分割参数：组件模型中填充几何体的总数目N＝64，一个网格可填充几何体的最大数目N_max＝20，一个X型网格可填充几何体的最大数目N_X,max＝5，一个Y型网格可填充几何体的最大数目N_Y,max＝2，一个Z型网格可填充几何体的最大数目N_Z,max＝1。

为了实现基于组件网格模型的快速蒙特卡罗粒子输运几何处理，主要包括以下步骤：

(1)计算组件模型中的所有几何体的轴向包围盒。根据轴向包围盒的定义，很容易知道每根燃料元件的轴向包围盒尺寸均为1.7325cm×1.7325cm×380cm，整个组件模型的轴向包围盒尺寸即为模型的几何尺寸13.86cm×13.86cm×380cm。

(2)确定网格分割的主轴、次轴和侧轴，并依次进行网格分割。

第一步，根据包围盒数据和网格分割参数确定主轴，并沿主轴对组件模型进行网格分割。当沿X轴进行网格分割时， $X_{\mathrm{step}}=\underset{J=1~N}{\mathrm{Min}}\left\{\right|X_{J,\max }-X_{J,\min }\left|\right\}=1.7325,$ 则代价函数：

$F\left(X\right)=\frac{\underset{J=1~N_X}{\mathrm{\&Sigma;}}{M}_{J,X}}{\underset{J=1~N_X}{\mathrm{\&Sigma;}}{N}_{J,X}}=\frac{8\&times;8}{8}=8(N_X=8)$

当沿Y轴进行网格分割时， $X_{\mathrm{step}}=\underset{J=1~N}{\mathrm{Min}}\left\{\right|X_{J,\max }-X_{J,\min }\left|\right\}=1.7325,$ 则代价函数：

$F\left(Y\right)=\frac{\underset{J=1~N_Y}{\mathrm{\&Sigma;}}{M}_{J,Y}}{\underset{J=1~N_Y}{\mathrm{\&Sigma;}}{N}_{J,Y}}=\frac{8\&times;8}{8}=8(N_Y=8)$

当沿Z轴进行网格分割时， $Z_{\mathrm{step}}=\underset{J=1~N}{\mathrm{Min}}\left\{\right|Z_{J,\max }-Z_{J,\min }\left|\right\}=380,$ 则代价函数：

$F\left(Z\right)=\frac{\underset{J=1~N_Z}{\mathrm{\&Sigma;}}{M}_{J,Z}}{\underset{J=1~N_Z}{\mathrm{\&Sigma;}}{N}_{J,Z}}=\frac{1\&times;64}{1}=64(N_Z=1)$

三个代价函数关系：F(X)＝F(Y)<F(Z)，根据代价函数最小原则，可以选X或Y轴为主轴进行网格分割，此处假定选择X轴为主轴进行网格分割。即按照最优步长X_step＝1.7325沿X轴进行分割后，得到8个X型网格，每个网格均含有8个几何体。

第二步，对包含几何体数目大于N_X,max＝5的每个X型网格，对其进一步确定次轴，并进行网格分割。由第一步可知，每个X型网格的包含的几何体数目(M_J,X＝8)均大于N_x,max，因此每个X型网格均需进一步网格分割。以第一个X型网格为例，此时有：

当沿Y轴进行网格分割时， $Y_{\mathrm{step}}=\underset{J=1~N}{\mathrm{Min}}\left\{\right|Y_{J,\max }-Y_{J,\min }\left|\right\}=1.7325,$ 则代价函数：

$F\left(Y\right)=\frac{\underset{J=1~N_Y}{\mathrm{\&Sigma;}}{M}_{J,Y}}{\underset{J=1~N_Y}{\mathrm{\&Sigma;}}{N}_{J,Y}}=\frac{8\&times;1}{8}=1(N_Y=8)$

当沿Z轴进行网格分割时， $Z_{\mathrm{step}}=\underset{J=1~N}{\mathrm{Min}}\left\{\right|Z_{J,\max }-Z_{J,\min }\left|\right\}=380,$ 则代价函数：

$F\left(Z\right)=\frac{\underset{J=1~N_Z}{\mathrm{\&Sigma;}}{M}_{J,Z}}{\underset{J=1~N_Z}{\mathrm{\&Sigma;}}{N}_{J,Z}}=\frac{1\&times;8}{1}=8(N_Z=1)$

两个代价函数关系：F(Y)<F(Z)，因此根据代价函数最小原则，可以选Y轴为次轴对该网格进行分割，即按照最优步长Y_step＝1.7325沿Y轴对第一个X型网格进行分割，分割后得到8个Y型网格，每个网格只包含1个几何体。同理，对剩下7个X型网格可进行同样分割。

第三步，对填充几何体数目大于N_Y,max＝2的每个Y型网格，对其进一步确定侧轴，并进行网格分割。在本计算实例中，由第一、二步按照主轴和次轴网格分割后，容易知道所有Y型网格只包含1个几何体，小于N_Y,max，故可以不再继续分割，即可以认为每个Y型网格以Z轴为侧轴只划分了1个Z型网格；此时得到的网格模型即为最优网格模型(如图7)。

(3)根据第(2)步得到的组件几何模型的最优网格模型，对于实例中给定的一空间点P(7.79625,7.79625,190)进行快速几何体定位和几何步长计算。

对于给定空间点P(7.79625,7.79625,190)，即X_p＝7.79625，Y_p＝7.79625，Z_p＝190。首先根据点P的X坐标确定主轴网格编号：

再根据点P的Y坐标，在第5个主轴网格中确定次轴网格编号：

对每个次轴网格而言，由于在Z轴上无进一步网格分割，只有1个网格，故侧轴网格编号为：

ID(Z)＝1

故点P在网格模型中的所在网格的编号为(5,5,1)，从图7中可以看到，该网格的包含的几何体为水棒，从而根据点P坐标与水棒几何面方程之间的关系，只需1次比较过程，便可确定P就在(5,5,1)网格所在的水棒几何中。而传统的蒙特卡罗粒子输运几何处理方法，需要对整个模型中的所有几何体进行比较，而基于网格模型只需定位网格后进行数次比较便可快速确定点P的所在的几何体。在得到点P所在几何体后，便可进行几何步长计算。

本发明未详细阐述的部分属于本领域公知技术。

尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

Claims (1)

1.一种基于空间网格分割的蒙特卡罗粒子输运快速几何处理方法，其特征在于，包括以下步骤：

步骤(1)计算几何模型中各个几何体的轴向包围盒：

分析待处理的几何模型，根据计算机图形学中的轴向包围盒AABB技术，为几何模型中的各个几何体构建轴向包围盒，轴向包围盒为各边平行于笛卡尔直角坐标轴的平行六面体，每个包围盒有且完全包含一个几何体，因此描述一个几何体的轴向包围盒只需六个标量即可，一个几何体的轴向包围盒可表示为：

AABB_J＝[(X_J,min,X_J,max)；(Y_J,min,Y_J,max)；(Z_J,min,Z_J,max)](J＝0～N)

其中J表示几何模型中几何体的编号，N表示整个模型中几何体的总数，AABB₀表示整个模型的轴向包围盒，AABB_J表示模型中第J个几何体的轴向包围盒，(X_J,min,X_J,max)表示几何体J沿X轴的尺寸范围从坐标值X_J,min到坐标值X_J,max，(Y_J,min,Y_J,max)表示几何体J沿Y轴的尺寸范围从坐标值Y_J,min到坐标值Y_J,max，(Z_J,min,Z_J,max)表示几何体J沿Z轴的尺寸范围从坐标值Z_J,min到坐标值Z_J,max；

步骤(2)基于代价函数最小原则，对给定几何模型进行空间网格分割，其步骤如下：

(21)根据步骤(1)得到的整个几何模型的轴向包围盒AABB₀，先假设对其沿X、Y、Z轴按照一定步长依次进行空间网格分割，其分割步长一般选取沿各轴方向的某一最小几何体的包围盒尺寸作为沿对应轴分割的步长，由如下公式获得：

$X_{step}=\underset{J=1~N}{Min}\left\{{\mathrm{AABB}}_{J,X}\right\}=\underset{J=1~N}{Min}\left\{\right|X_{J,max}-X_{J,min}\left|\right\},$

$Y_{step}=\underset{J=1~N}{Min}\left\{{\mathrm{AABB}}_{J,Y}\right\}=\underset{J=1~N}{Min}\left\{\right|Y_{J,max}-Y_{J,min}\left|\right\},$

$Z_{step}=\underset{J=1~N}{Min}\left\{{\mathrm{AABB}}_{J,Z}\right\}=\underset{J=1~N}{Min}\left\{\right|Z_{J,max}-Z_{J,min}\left|\right\},$

其中X_step、Y_step、Z_step为将整个几何模型沿X、Y、Z轴进行空间网格分割的步长；为从编号为1到N的几何体中选择沿X轴的包围盒最小尺寸，为同理，为从编号为1到N的几何体中选择沿Y轴的包围盒最小尺寸，为 为从编号为1到N的几何体中选择沿Z轴的包围盒最小尺寸，为

(22)假设以X_step、Y_step、Z_step为步长依次对整个几何模型沿X、Y、Z轴进行空间网格分割，得到三种空间网格分割方案，这三种空间网格分割方案得到的网格数目N_X、N_Y、N_Z分别为：

$N_X=\frac{|X_{0,max}-X_{0,min}|}{X_{step}}$

$N_Y=\frac{|Y_{0,max}-Y_{0,min}|}{Y_{step}}$

$N_Z=\frac{|Z_{0,max}-Z_{0,min}|}{Z_{step}}$

其中，X_0,min和X_0,max为整个模型沿X轴的最小尺寸和最大尺寸坐标值，Y_0,min和Y_0,max为整个模型沿Y轴的最小尺寸和最大尺寸坐标值，Z_0,min和Z_0,max为整个模型沿Z轴的最小尺寸和最大尺寸坐标值；三种网格分割方案的代价函数F(X)、F(Y)、F(Z)分别为：

$F\left(X\right)=\frac{\underset{J=1~N_X}{\&Sigma;}{M}_{J,X}}{\underset{J=1~N_X}{\&Sigma;}{N}_{J,X}}$

$F\left(Y\right)=\frac{\underset{J=1~N_Y}{\&Sigma;}{M}_{J,Y}}{\underset{J=1~N_Y}{\&Sigma;}{N}_{J,Y}}$

$F\left(Z\right)=\frac{\underset{J=1~N_Z}{\&Sigma;}{M}_{J,Z}}{\underset{J=1~N_Z}{\&Sigma;}{N}_{J,Z}}$

其中，F(X)为沿X轴分割得到的网格模型的代价函数，F(Y)为沿Y轴分割得到的网格模型的代价函数，F(Z)为沿Z轴分割得到的网格模型的代价函数；M_J,X、M_J,Y、M_J,Z表示分别沿X、Y、Z轴进行空间网格分割得到的网格模型中第J个网格中包含的几何体数目，N_J,X、N_J,Y、N_J,Z表示分别沿X、Y、Z轴进行空间网格分割得到的网格模型中第J个网格是否为空网格，即该网格是否含有几何体，它们的关系如下：

$N_{J,X}=\left\{\begin{array}{cc}1& if(M_{J,X}>0)\\ 0& if(M_{J,X}=0)\end{array}\right.$

$N_{J,Y}=\left\{\begin{array}{cc}1& if(M_{J,Y}>0)\\ 0& if(M_{J,Y}=0)\end{array}\right.$

$N_{J,Z}=\left\{\begin{array}{cc}1& if(M_{J,Z}>0)\\ 0& if(M_{J,Z}=0)\end{array}\right.$

由此可知，代价函数的意义在于衡量一个给定的空间网格分割模型，平均每个网格包含的几何体数目；

(23)选择代价函数最小的网格分割为最优的空间网格分割方案：

根据步骤(22)得到的三种空间网格分割方案，选取代价函数最小的Min{F(X),F(Y),F(Z)}方案作为步骤(2)最终的空间网格分割方案；

步骤(3)针对步骤(2)得到的空间网格分割得到的每个网格，再重复步骤(2)，并依次递归下去，直到待分割的网格中包含的几何体数目少于一定数目时，停止分割；步骤(3)和步骤(2)最大的区别在于：步骤(2)以整个几何模型为分割对象，而步骤(3)则针对某一个特定网格继续进行空间网格分割，但其本质是一样的；

步骤(4)基于空间网格分割得到的网格模型进行蒙特卡罗粒子输运快速几何处理，查找粒子所在的几何体：

(41)根据粒子坐标，确定包含该粒子坐标的点的网格编号；在进行蒙特卡罗粒子输运模拟时，在给定粒子三维空间坐标数据P(X_p,Y_p,Z_p)时，首先根据以下公式，快速进行网格定位，查找该粒子坐标的点所在的网格编号：

其中C_i为点P(X_p,Y_p,Z_p)在坐标轴i方向上的分量，L_i为整个模型或某一网格沿坐标轴i的最小边界坐标值，|C_i-L_i|为给定点沿坐标轴i方向到模型边界的距离，S_i为整个模型或某一网格沿坐标轴i进行网格分割的步长，表示取大于的最小正整数，则[ID(X),ID(Y),ID(Z)]即为点P(X_p,Y_p,Z_p)所在的网格编号；

(42)假设网格[ID(X),ID(Y),ID(Z)]包含的几何体的编号依次为G₁,G₂,...,G_M，然后依次对这M个几何体进行点包容性计算，从而便可快速查找到粒子所在的几何体。