CN113221200B - 一种适用于堆芯颗粒分布不确定性分析的三维高效随机排布方法 - Google Patents

Abstract

本发明涉及核反应堆建模与仿真领域，涉及适用于堆芯颗粒分布不确定性分析的三维高效随机排布方法，步骤如下：一：获取计算参数信息；二：在几何区域内随机生成出所有颗粒位置；三：给几何区域划分网格，把所有颗粒归并到各个网格当中；四：判断颗粒之间是否存在重叠；五：消除颗粒间重叠；六：判断颗粒与边界之间是否存在重叠；七：消除颗粒与边界的重叠；八：重复步骤四到步骤七，直到所有颗粒都消除了重叠情况；步骤九：判断颗粒和几何尺寸，若没达到规定尺寸，则增长颗粒半径，减小几何尺寸，并更新网格内颗粒信息，回到步骤三继续进行。本发明只要获取了堆芯的几何尺寸和燃料颗粒数，即可得出最终的全堆芯模型。

Description

一种适用于堆芯颗粒分布不确定性分析的三维高效随机排布 方法

技术领域

本发明涉及核反应堆建模与仿真领域，具体涉及适用于堆芯颗粒分布不确定性分析的三维高效随机排布方法。

背景技术

三维高效随机排布方法通过计算机模拟出颗粒在几何容器内部的真实填充情况，可以快速方便的来研究各种几何形状容器内部的颗粒的随机填充位置、填充率与空隙率、颗粒堆积结构等问题。此外，还可为后续的模拟计算提供相应参数，如颗粒的球流运动，颗粒排布的不确定性分析等。这尤其在核反应堆中燃料颗粒的随机排布不确定性分析中具体重要技术意义。

目前核反应堆建模与计算软件中MCNP(Benoit Forestiera,Joachim Missa,Franck Bernard.Criticality calculations on pebble-bed HTR-PROTEUSconfiguration as a validation for the pseudo-scattering tracking methodimplemented in the MORET5Monte Carlo code.International Conference on thePhysics of Reactors,Interlaken,Switzerland,September 14-19,2008.),Serpent(Ville Rintala,Heikki Suikkanen,Jaakko

Modeling of realistic pebblebed reactor geometries using the Serpent Monte Carlo code,Annals of NuclearEnergy,v77,p.223-230,2015.),VSOP，Dragon都有产生燃料颗粒的随机坐标的方法，而颗粒随机排布方法主要有蒙特卡罗移除方法(Monte Carlo Rejection Method)，系统扩展法(Expending System Method)，离散单元法(Discrete Elements Method)等。

计算机模拟仿真随机球排布方法一般分为随机松散排布(Random LoosePacking,RLP)和随机紧密排布(Random Close Packing,RCP)(W.S.Jodrey,E.M.Tory,Computer simulation of close random packing of equal spheres,PHYSICAL REVIEW,v32，p.4，1985.)；随机松散排布一般用于没有重力作用的弥散模型，这些模型一般颗粒填充率较低，易于实现；而随机紧密排布模型一般是填充率大于60％的模型。对于随机紧密排布，目前分为基于物理的方法与基于几何的模拟方法。其中，基于物理的模拟方法目前主要有离散单元法(Discrete Element Method，DEM)和分子动力学法(Molecular Dynamics,MD)。离散单元法主要是划分单元，单元内计算颗粒与颗粒之间的各种相互作用力和重力作用，按照牛顿运动定理计算颗粒的运动方向与速度大小，是一种严谨的方法，目前广泛应用于土木和材料工程邻域。基于几何的方法主要有序列添加法和集合重排法两类，序列添加法首先放置一定量颗粒在容器中，考虑重力所以颗粒一般都放置在容器底部，然后每次迭代生成一个颗粒按照一个给定的方向运动(一般为重力方向)，直到颗粒碰到边界或其他颗粒，达到稳定位置后停止，如果生成的颗粒与边界或者与其他颗粒重叠，则重新继续生成，比如蒙特卡洛移除法(Monte Carol Rejection Method,MC Rejection)(G.J.Auwerda,J.L.Kloosterman,Comparison of Experiments and Calculations of Void FractionDistributions in Randomly Stacked Pebble Beds,In:PHYSOR-2010,Pittsburgh,Pennsylvania,USA,May 9-14,2010.)；集合重排法则是在几何容器内一次性投入所有颗粒的坐标点，不考虑颗粒实际的重叠情况，再根据几何判断条件挪动颗粒，直到所有颗粒消除重叠位置，主要有代表性方法的就是扩展系统法(Expending System Method,ESM)，该方法在集合重排的基础上，增加了增长球半径的功能，使得排列可以更为紧密，最高可达到64％填充率(E.I.Zubko1 and Yu.E.Zubko,MAXIMUM POSSIBLE DENSITIES OF RANDOM SPHEREPACKING WITHIN THE COMPOSITE,Powder Metallurgy and Metal Ceramics,v58,p.133,2019)。

现有的产生随机颗粒的方法，通常有以下几点缺陷：颗粒的随机性不强，或者不能进行随机排布只能进行规则排布；随机紧密排布的颗粒填充率不够高；随机排布结果运行计算时间长；随机颗粒排布不能反映真实的模型情况(有颗粒重叠或者颗粒悬空)；不能在几何尺寸，随机颗粒数量确定的情况下按照一定填充率随机填充；不能输出颗粒随机位置来进行不确定性分析。

对于当前核反应堆建模计算软件当中，DRAGON程序只能进行颗粒的规则排布，对于随机排布则存在颗粒重叠问题，不能反应真实情况。对于Serpent和MCNP软件，都采用的蒙特卡罗方法产生的随机颗粒坐标位置，根据蒙卡的产生原理，对于随机紧密排布的堆积情况，如HTR-10，PBR-250球床中燃料球的堆积，往往达不到相应的填充率要求，而且因为蒙卡方法随机位置的产生，程序的计算时间长。基于几何的随机排布方法比如ESM方法可以得到较高的填充率但是没有考虑重力因素的影响可能会在容器顶部存在颗粒悬空的非真实情况，而基于物理的随机排布方法如DEM方法可以真实反应颗粒的填充，但是因为需要输入颗粒质量，摩擦系数等一系列物理参数，在每一步中都需要计算容器内每个颗粒的受力情况而使得计算会占用大量的计算机内存，计算时间也大大增加。

对于基于物理的随机排布方法，比如DEM方法，QDM方法(Yanheng Li,WeiJi.Stability and convergence analysis of a dynamics-based collective methodfor random sphere packing.Journal of Computational Physics.v250,p.373-387，2013.)，因为考虑到了重力、颗粒间的相互作用力，颗粒间的摩擦力、颗粒和容器壁的摩擦力等等，所以主要的问题在于计算机的运行时间过长，而且颗粒数越多，填充率越高的模型，计算时间往往成指数上涨。因此为了解决时间上的问题，主要会从程序结构和计算机性能上着手；程序的优化和计算机的性能则需要花费较多的资金。

发明内容

针对已有技术的缺陷和不足，本发明人探索采用了基于几何的随机排布方法，不考虑受力的影响，只有几何上的判断而使得程序变得十分简便，从而使得计算效率大大提高，但是也因为没有重力作用，在填充率较低情况下会出现颗粒悬空情况，为了能够模拟出真实排布情况往往会在颗粒的几何排布方法上进行研究，比如ESM方法的每一迭代步增大一定颗粒尺寸来减少颗粒间的空隙，从而达到更高的填充率。但是ESM方法往往填充率和模型真实性不能完全保证，为了保证真实性通常最终模型只会选择几何中间一段而去掉模型边界两端的一些颗粒，这样最终模型的颗粒数和填充率会和预期值产生一定偏差从而达不到精确的设计要求。

进而，本发明目的就是为了能够计算真实堆芯中燃料颗粒或燃料球的随机排布造成的不确定性影响，需要保证在规定几何尺寸下的颗粒数量和填充率要求，而且在保证真实度的情况下还有不错的计算效率，在保证真实排布的情况下还采用了一种较快的搜索颗粒方法来加速程序的运行速度。

本发明是一种改进的ESM方法，属于基于几何的随机排布方法。一方面，对于几何随机排布方法的缺陷，比如蒙卡方法的低填充率问题，该方法通过不断增加颗粒的直径，再根据颗粒间的位置判断不断挪动颗粒位置使之达到更为紧密的填充。另一方面，对于不能反应真实排布问题，有颗粒悬空的问题，该方法通过压缩手段实现，即先在比目标容器更大的几何容器中产生出所有颗粒，通常在高度尺寸上放大，然后在每一次迭代之后，增大一定球半径的同时减小一定的容积几何尺寸，最终高度和球半径都变为目标模型的尺寸后结束运行。这样相当于在边界处给予了一个“挤压”使得填充更加紧密而且不会出现颗粒悬空的情况。再一方面，对于计算时间的优化，本发明使用了一种网格近邻表方法来加快颗粒之间的比较，从而大大加快了颗粒数多填充率高情况下的计算效率。

本发明提出了一种基于网格近邻表的ESM方法，目的是为了生成真实的反应堆内燃料球和燃料颗粒的随机排布，从而研究燃料球、燃料颗粒随机排布对堆芯计算所造成的不确定性。本发明方法的实现步骤如下：

步骤一：获取计算参数信息，包括颗粒总数(即几何区域内颗粒的数量，如燃料球数、燃料颗粒数)，最终边界几何尺寸(对于圆柱、立方体几何只需要压缩几何高度即可，对于球型几何则需压缩球型几何半径)H_end，最终颗粒半径R_peb和迭代步数N。

步骤二：在几何区域内随机生成出所有颗粒位置，并给颗粒编号，不考虑重叠情况；

步骤三：给几何区域划分网格，把所有颗粒归并到各个网格当中，网格中记录颗粒编号；

步骤四：判断颗粒之间是否存在重叠(颗粒间距离小于颗粒直径)；

步骤五：消除颗粒间重叠；

步骤六：判断颗粒与边界之间是否存在重叠(颗粒与边界距离小于颗粒半径)；

步骤七：消除颗粒与边界的重叠；

步骤八：重复步骤四到步骤七，直到所有颗粒都消除了重叠情况；

步骤九：判断颗粒和几何尺寸，若没达到规定尺寸，则增长颗粒半径，减小几何尺寸，并更新网格内颗粒信息，然后回到步骤三继续进行；若已达到规定的颗粒和几何尺寸，步骤结束，得出最终的随机填充颗粒坐标位置，颗粒随机排布完成。

通过本发明，对于核反应堆内的燃料颗粒随机排布，只要获取了堆芯的几何尺寸、燃料颗粒尺寸和燃料颗粒数，通过上述步骤即可得出最终的全堆芯模型。本发明是在传统的ESM方法上，加入了网格近邻表法与几何压缩法。对于网格近邻表方法的整数特性可以快速搜索附近颗粒，大大加快了颗粒间的判断比较时间；对于几何压缩法，本发明使用了同半径增长一致的对数插值公式来压缩几何尺寸，从而在不用减少颗粒的情况下达到真实的填充模型，使得确定颗粒数目和几何尺寸的真实堆芯模型得以建立。

附图说明

图1本发明的ESM方法流程图。

图2网格近邻表方法示意图。

图3填充率随半径变化示意图。

图4HTR-10球床颗粒排布模型。

图5燃料球内TRISO颗粒随机排布模型。

图6该方法在RMC内产生的HTR-10初装堆模型。

具体实施方式

下面通过具体实施例对本发明作进一步的说明，但并不构成对本发明的限制。

实施例一

本发明方法产生真实堆芯颗粒随机排布流程如图1所示。

具体步骤如下：

(1)获取堆芯内随机排布的颗粒总数N_peb，最终颗粒半径R_peb，边界几何尺寸和总迭代步数N这些必要的参数信息，在该方法中颗粒半径和几何尺寸会随着迭代步数而发生变化，第i步的颗粒半径表示为R_i，几何高度表示为H_i。通常最初颗粒半径R_ini定义为最终颗粒半径R_peb的2/3，最初几何高度H_ini为最终几何高度H_end的3/2。

(2)通过基于线性同余法的rand随机数生成函数在几何区域内生成N_peb个点作为各随机颗粒的球心，不考虑颗粒之间的重叠情况；

(3)判断球与球之间的位置关系，本发明所采用的基于网格近邻表的ESM方法中，首先在颗粒填充的几何区域内划分网格，网格宽度一般为一倍球直径或两倍球直径，使用者可以自行设定，但必须大于等于一倍球直径。采用笛卡尔坐标系，在X，Y，Z方向分别划分网格，使得几何空间被分割为一个个立方体晶格。定义如果球心坐标在晶格之内，该球就属于该晶格，因为晶格宽度大于等于一倍球直径，所以一个球只可能在一个晶格之中，因此可以用晶格坐标存储球的信息。具体存储方式如下：

首先给每一个球编号(1,2,3……N_peb)，每个晶格也可以由一组(m,n,l)整数编号表示，其中m，n，l为对应方向上网格宽度的倍数，例如一个编号为k(1≤k≤N_peb)的球的球心坐标为(7.625，25.395，11.431)，假设球当前直径D＝2R_i＝6，X，Y，Z方向上网格划分宽度都为一倍球直径，则该球所属的晶格为(2，5，2)，以此类推，所有球都可以归到各个由三个整数坐标表示的晶格中，随后通过比较该晶格内的和其周围26个晶格内的球位置进行判断即可，如图2所示。

颗粒间的具体判断流程为，比较该球和其周围26个晶格中的球的球心之间距离，定义球心之间距离为l，球当前直径为D＝2R_i，采用两点间距离公式计算出两球球心之间的距离l，如果l<D，则说明球心之间的距离小于球直径，那么两球处于重叠情况。

对于传统ESM的方法，由于没有加入几何区域划分的网格法，会在所有球中进行两两之间的距离判断，虽然距离判断方法和本发明方法一致。但是在在球数较多的情况下,传统的ESM方法在进行颗粒之间的距离比较时往往会耗费大量的时间，而本发明方法在传统的ESM上使用网格近邻表方法加速了这一过程，使得计算效率大大提高。

(4)消除颗粒间的重叠，在上一步中计算了颗粒间球心的距离l，当两球球心之间的距离l小于当前球直径D＝2R_i，则判定两球处于重叠情况，而在真实的颗粒排布模型中，特别是在反应堆中，颗粒都为刚性球体，不可能出现这样的颗粒重叠现象，所以需要对重叠的颗粒进行消除重叠操作，具体方法为当判断l<D时，使两球在其球心连线方向上分别相背移动(D–l)/2的距离使其变为相切状态。由于每次只能进行两个重叠颗粒之间的移动，当模型内颗粒数量很多，紧密堆积(填充率>60％)的情况时，一个颗粒可能和多个颗粒发生重叠现象，通过一次消除重叠的移动操作不一定能消除该颗粒的重叠现象，因此每次消除重叠移动步骤后，都需要返回上一步的颗粒间的位置重新再进行比较。

(5)在上一步颗粒进行了消除重叠移动后，如果再次判断颗粒间的距离已经完全消除了重叠后，下一步会进行颗粒和边界间的位置判断。通过点到直线距离公式计算出颗粒球心和边界之间的距离L(如果边界为圆边，则L为球心到圆边切线的距离)。当颗粒与边界之间的距离小于当前颗粒半径，即L<R_i时，判断为颗粒和边界发生了重叠现象，此时该方法通过使颗粒沿着其边线法线方向向边线内移动R_i-L的距离，移动后颗粒与边界成为相切状态。同样，在进行颗粒的消除重叠的移动操作后，颗粒位置发生了变化，所以需要返回颗粒间的位置比较步骤重新比较位置。

在传统的ESM方法中，因为程序每次只能进行两个球之间的顺序位置判断，因此可能会出现球与球，球与边界多次重叠和消除重叠移动的情况。例如球A与球B判断完移动消除重叠后，球A因为发生了移动可能和之前没有重叠的球C产生重叠，球A与球C判断移动消除重叠后，球A可能又会与之前消除了重叠的球B发生重叠。当两个球球心连线与坐标轴接近平行且边界垂直于坐标轴时，可能会出现球与球、球与边界之间来回不断的判断与移动操作，从而最终可能导致死循环的发生。而本发明采用的是一种网格近邻表方法，通过给每个球和网格进行整数值定位，通过每次颗粒判断移动操作后重新更新网格内的球编号信息来避免死循环的发生，而且在来回不断的颗粒间判断操作中提高了计算的效率。

(6)当所有颗粒都进行了比较并且没有了重叠情况后，判断颗粒当前半径R_i是否达到输入的颗粒最终半径R_peb，几何高度H_i是否达到最终高度H_end；如果没有，则通过对数插值公式(1)增长颗粒的半径，通过公式(2)减小几何高度，随后重新进入步骤(3)进行下一步迭代计算。

式中R_i，H_i分别为第i次迭下的颗粒半径和几何高度；R_ini，H_ini分别为最初颗粒半径和最初几何高度；R_peb，H_end分别为最终球半径和最终几何高度，i为当前迭代次数，N为总迭代次数，通常

不断进行迭代计算，直到判断出R_i＝R_peb；H_i＝H_end时，结束计算，颗粒随机排布完成。

本发明因为其实用性和高效性，可以产生基于RMC程序的颗粒坐标来进行燃料颗粒随机位置的不确定性分析。

验证实施例

因为本发明方法没有加入重力作用，方法验证是为了验证最终生成的颗粒模型是否符合真实的颗粒堆积情况，从而才能说明应用该方法生成的HTR-10等堆芯模型是可以模拟真实情况的。而对比方法是在验证了方法的正确性之后，通过和其他方法的对比体现出该方法的优势。以下以G.J.Auwerda,J.L.Kloosterman,D.Lathouwers,et al.,Packingmicrostructure and local density variations of experimental and computationalpebble beds.In:PHYSOR-2012.Knoxville,TN.2012.这篇论文中的数据，通过不同方法模拟出真实搭建的模型，本发明方法也按文中提供的数据构造出同样的几何条件填充模型，从耗时上说明了该方法的高效性，从填充率上说明了该方法的可操作性和准确性。

将本发明方法与真实实验测量的结果进行了对比，如图3所示。图3为颗粒填充率随几何半径从边缘到中心的变化曲线，数据点为Gert等人在代尔夫特工业大学搭建的一个真实颗粒实验模型装置(PebBeX facility)测量得到的(如上述论文)。真实模型采用的是亚克力作为材料制作的一个圆柱型容器，不断往容器内投放同样亚克力材料的小球，直到小球填满整个容器为止。圆柱容器的几何直径和填充小球的几何直径比D/d＝18,Gert采用X光线扫描的方式探测模型内部的颗粒填充情况去获得数据，从而得到了图3中的数据点。而黑色的数据曲线为本发明方法产生和实验模型一样尺寸和球数的随机颗粒填充模型的结果，可以看出填充率随几何半径的变化结果与Gert在真实颗粒排布中测量的数据结果吻合，可以说明本发明方法产生模型的真实性。

在表1中，进行了在球数和几何条件都一致的情况下，各种方法计算时间的对比,其中圆柱容器几何的内直径为229mm，高度为235mm，填充的小球直径为12.7mm,填充球数为5457±10个。

其中，实验结果是Gert搭建的真实模型的结果，DEM为严谨的物理计算方法(Yanheng Li，Wei Ji.A collective dynamics-based method for initial pebblepacking in pebble flow simulations.Nuclear Engineering and Design.v250,p:229-236,2012.)，在DEM方法中，需要获取实验模型中亚克力材料的几何容器和填充的球的物理参数，如重力加速度，球的质量、密度，摩擦系数，杨氏模量等等物理信息。随后根据牛顿第二定律计算每个时间点下各个小球的受力情况，还原球真实的填充过程和运动情况，最终得出填充的模型。传统ESM为纯几何的方法(P.Mrafko,Homogeneous and Isotropic HardSphere Model of Amorphous Metals,Journal de Physique Colloques,v44,pp.C8-322-C8-325,1980)，相对于本发明方法没有加入网格近邻表和几何压缩法。在相同模型条件的计算结果中可以看出，本发明方法在球数和填充率严格控制的情况下保证了模型的准确性而且计算时间最少。在利用抽样统计方法进行不确定性分析中，大量样本的计算需要花费较多时间，进而通过本发明方法可以更加有效率的产生多组随机颗粒坐标样本来进行反应堆内燃料球随机排布和随机填充率的不确定性研究。

表1各方法对比数据表

本发明人还研究发现，如果没有采用网格近邻表方法的ESM方法，在2D低填充率情况下，还可以快速获得准确的结果，但是在3D的高填充率、高D/d比情况下，计算时间显著增加，而且可能会出现死循环的情况。因此，本发明加入了网格近邻表方法的ESM方法，直接效果就是加快了程序颗粒之间对比的速率，而更深层次的效果在于因为网格的整数定义特性，使得程序在计算效率提高的同时变得更加的可观与简化，通过不同情况下更改网格的宽度可以达到一个最优的效率，而且解决了死循环的问题从而提升了程序的稳定性。更因为计算效率的显著提高，用抽样方法进行反应堆内随机颗粒排布的不确定性分析可以变的更加的准确高效。

应用实施例一：

采用本发明方法生成了真实HTR-10初装堆堆芯的高保真燃料填充模型，用来解决燃料球在堆芯内的随机堆积。

第一步，获取反应堆的堆芯尺寸和燃料球、石墨球尺寸与数目等参数，其中具体有堆芯活性区圆柱直径180cm，高度设置为123cm，下圆台高度36.946cm，下方卸料管直径25cm，长度33cm，燃料球与石墨球直径均为6cm，总球数N_peb总共有16890个，迭代步数N设置为3步。

第二步，按照相应的几何参数在设置的几何区域内生成出16890个颗粒的随机坐标位置，网格宽度设置为当前一倍球直径D_i，通过把所有颗粒信息统计到区域的各个网格当中。

第三步，通过网格近邻表方法进行颗粒间的位置判断，即使用两点之间距离公式计算两颗粒球心之间距离l是否小于当前球直径D_i。

第四步，如果判断两球之间距离l<D_i，则两球沿着球心沿线方向相背移动(D-l)/2的距离成为相切状态，通过该方法消除所有颗粒间的重叠情况。

第五步，判断颗粒与边界的重叠情况，判断颗粒球心与边界的距离L是否小于R_i，如果L<R_i，则颗粒往边界内部移动R_i-L的距离使得颗粒与边界相切。

第六步，所有颗粒消除和边界的重叠后回到第四步重新进行颗粒间的判断，直到所有颗粒都消除了重叠情况，通过对数插值公式(1)与公式(2)增加颗粒半径和减小活性区大圆柱的高度。

第七步，判断迭代步数是否达到最终迭代步数3，颗粒半径R_i是否达到最终半径3cm，大圆柱几何高度是否达到最终高度123cm，如果没有，继续从第三步开始迭代计算，并且同时更新网格内颗粒的信息。如果达到最终几何半径和高度，得出最终的随机填充颗粒坐标位置，通过每个球的位置坐标和半径，用绘图软件画出图像如图4所示。

有了颗粒的坐标位置，即可以把颗粒的坐标位置写成RMC建模软件的输入卡，输入卡中有颗粒的编号(从1到16890)，颗粒的半径3cm，颗粒的三维坐标位置(x,y,z)。RMC通过读取颗粒坐标位置输入卡和其相应的几何与材料信息即可生成出颗粒随机排布的真实HTR-10堆芯模型，如图6所示。因为本发明方法采用了线性同余法随机函数产生颗粒随机坐标位置，因此每次生成的颗粒随机填充的位置都不一样，通过不断生成颗粒坐标位置得到多组燃料颗粒随机坐标样本，再通过RMC的中子输运计算算出每个样本模型的k_eff值，从而可以分析出燃料颗粒的随机排布对k_eff值的不确定性。此外，在燃料颗粒数不变的情况下采用更改填充几何的最终高度来更改燃料的填充率，同样多次产生多组不同填充率下的随机坐标样本，通过抽样统计的方法研究填充率的变化对k_eff值的不确定性。

应用实施例二：

应用本发明方法产生了球型边界颗粒随机排布，用来模拟反应堆燃料球内TRISO燃料颗粒的随机弥散，如图5所示。模型边界为球型边界情况，最终燃料球半径为3cm，内部TRISO燃料颗粒最终半径为0.91mm，在HTR-10反应堆中，燃料球内TRISO颗粒个数为8335。获取到这些几何信息后，设置迭代步数N为3既可以通过以上步骤产生出燃料球内燃料颗粒的随机弥散排布模型，需要注意的是填充的几何边界变为了球型边界，所以最终的几何尺寸为球型边界的最终半径。通过多次使用产生即可得到多个不同排布的颗粒弥散模型。该模型可以真实的模拟反应堆中燃料球内燃料颗粒的随机排布，从而研究每个燃料球内燃料颗粒的随机弥散对反应堆物理计算造成的不确定性。

Claims (8)

1.一种适用于堆芯颗粒分布不确定性分析的三维高效随机排布方法，其特征在于，包括步骤如下：

步骤一：获取计算参数信息；

步骤二：在几何区域内随机生成出所有颗粒位置，并给颗粒编号，不考虑重叠情况；

步骤三：给几何区域划分网格，把所有颗粒归并到各个网格当中，网格中记录颗粒编号；

步骤四：判断颗粒之间是否存在重叠；

步骤五：消除颗粒间重叠；

步骤六：判断颗粒与边界之间是否存在重叠；

步骤七：消除颗粒与边界的重叠；

步骤八：重复步骤四到步骤七，直到所有颗粒都消除了重叠情况；

步骤九：判断颗粒和几何尺寸，若没达到规定尺寸，则增长颗粒半径，减小几何尺寸，并更新网格内颗粒信息，然后回到步骤三继续进行；若已达到规定的颗粒和几何尺寸，步骤结束；

其中，步骤四中判断颗粒之间是否存在重叠的方法如下：首先在颗粒填充的几何区域内划分网格，网格宽度为一倍球直径至两倍球直径；采用笛卡尔坐标系，在X，Y，Z方向分别划分网格，使得几何空间被分割为一个个立方体晶格；其中定义如果球心坐标在晶格之内，该球就属于该晶格，并且通过基于线性同余法的rand随机数生成函数在几何区域内生成Npeb个点作为各随机颗粒的球心，不考虑颗粒之间的重叠情况；比较该球和其周围26个晶格中的球的球心之间距离，定义球心之间距离为l，球当前直径为D＝2Ri，采用两点间距离公式计算出两球球心之间的距离l，如果l<D，则说明球心之间的距离小于球直径，那么两球处于重叠情况。

2.如权利要求1所述的三维高效随机排布方法，其特征在于，步骤一中参数信息包括颗粒总数，最终边界几何尺寸H_end，最终颗粒半径R_peb和迭代步数N。

3.如权利要求2所述的三维高效随机排布方法，其特征在于，步骤一中最初颗粒半径R_ini定义为最终颗粒半径R_peb的2/3，最初边界几何高度H_ini为最终几何高度H_end的3/2。

4.如权利要求1所述的三维高效随机排布方法，其特征在于，步骤四中，当判断得出颗粒间距小于颗粒直径，则认为颗粒间存在重叠。

5.如权利要求1所述的三维高效随机排布方法，其特征在于，步骤五中，当步骤四中计算确定的颗粒间球心的距离l小于当前球直径D(D＝2Ri)，则判定两球处于重叠情况，所述消除颗粒间重叠操作的具体方法为当判断l<D时，使两球在其球心连线方向上分别相背移动(D–l)/2的距离使其变为相切状态；每次消除重叠移动步骤后，都需要返回步骤四进行颗粒间是否重叠的再判断。

6.如权利要求1所述的三维高效随机排布方法，其特征在于，如果判断颗粒间的距离已经完全消除了重叠后，下一步会进行步骤六进行颗粒和边界间的位置判断，当判断得出颗粒与边界距离小于颗粒半径，则认为颗粒与边界存在重叠，具体方法是，通过点到直线距离公式计算出颗粒球心和边界之间的距离L；当颗粒与边界之间的距离小于当前颗粒半径，即L<Ri时，判断为颗粒和边界发生了重叠现象。

7.如权利要求1所述的三维高效随机排布方法，其特征在于，如果步骤六判断得出颗粒和边界发生了重叠，则通过使颗粒沿着其边线法线方向向边线内移动Ri-L的距离，移动后颗粒与边界成为相切状态；返回颗粒步骤三。

8.如权利要求1所述的三维高效随机排布方法，其特征在于，当完成步骤八，即所有颗粒都进行了比较并且没有重叠情况后，判断颗粒当前半径R_i是否达到输入的颗粒最终半径R_peb，几何高度H_i是否达到最终高度H_end；如果没有，则通过对数插值公式(1)增长颗粒的半径，通过公式(2)减小几何高度，随后重新进入步骤三进行下一步迭代计算，不断进行迭代计算，直到判断出R_i＝R_peb；H_i＝H_end时，结束计算，得出最终的随机填充颗粒坐标位置，颗粒随机排布完成：

式中R_i，H_i分别为第i次迭下的颗粒半径和几何高度；R_ini，H_ini分别为最初颗粒半径和最初几何高度；R_peb，H_end分别为最终球半径和最终几何高度，i为当前迭代次数，N为总迭代次数，

Images