三维随机球包装代码：基于Matlab的软件X 18（2022）101051

软件X 18（2022）101051原始软件出版物3DRSP：基于Matlab的三维随机球包装代码特拉维斯·J阿列克谢·布莱克切维亚科夫萨斯喀彻温大学数学与统计系，Saskatoon，SK，Canada S7N 5E6ar t i cl e i nf o文章历史记录：2021年10月29日收到2022年3月15日收到修订版，2022年保留字：随机球填充三维域概率分布Matlaba b st ra ct提出了一种基于svertMatlabsvert的计算程序，用半径为随机挑选从任何指明概率分布支持由svertMatlabsvert .用球面填充单位立方体、平行六面体和两个半球表面之间的凹域的通用程序序列和实例，从威布尔分布和伽玛分布，提出。考虑了在粉末床3D打印过程中由激光束加热的颗粒之间的键形成的数值模拟的示例应用©2022作者（S）。由爱思唯尔公司出版这是CC BY-NC-ND下的开放获取文章许可证（http://creativecommons.org/licenses/by-nc-nd/4.0/）。代码元数据当前代码版本1.0用于此代码版本的代码/存储库的永久链接https://github.com/ElsevierSoftwareX/SOFTX-D-21-00213可复制胶囊的永久链接法律代码许可证公共领域使用的代码版本控制系统无使用的软件代码语言、工具和服务Matlab 2021a编译要求、操作环境依赖性如果可用，链接到开发人员文档/手册问题支持电子邮件shevyakov@math.usask.ca1. 介绍所提出的基于Matlab的代码填充一个给定的三维域与球的半径遵循任何规定的随机分布的球体在Matlab中可用。创建此代码的动机是构建粉末床3D打印过程的离散元素方法（DEM）模型。这种模型的初始步骤是确定粉末床域中单个颗粒堆积的初始位置，以及哪些颗粒接触。这就提出了一个没有简单解决办法的问题。事实上，找到不相等球体的最佳包装是一项具有挑战性的任务[1]。以前的DEM模型（例如，参考文献[2，3]）依赖于粒子动力学方法来确定重力和粒子间力的位置;相互作用和*通讯作者。电子邮件地址：tjb235@usask.ca（Travis J. Black），shevyakov@math.usask.ca（Alexei F.Cheviakov）。https://doi.org/10.1016/j.softx.2022.101051使用牛顿力学确定球体的最终位置。这需要对于大量的粒子，这样的计算可能是非常昂贵的计算。由于我们唯一关心的是颗粒上的最终位置，因此当前代码被设计为绕过粉末沉降的动力学，并找到一种基于几何考虑来确定颗粒的静止位置的方法。随机球填充具有广泛的应用，包括DEM建模、颗粒动力学、用于治疗脑肿瘤的放射外科[1]、最优填充问题等。其他球填充方法具有类似的目的[4，5]，但是它们没有满足我们的需求大量文献致力于相同尺寸的球体的随机堆积的相关但不同的问题（参见，例如，参考文献[6][7]及参考文献。下面提出的程序的目标不是找到最佳填充，而是不相等球体的填充，其紧密地模拟金属粉末颗粒的实际填充。2352-7110/©2022作者。由爱思唯尔公司出版。这是一篇开放获取的文章，使用CC BY-NC-ND许可证（http://creativecommons.org/licenses/by-nc-nd/4.0/）。可在ScienceDirect上获得目录列表SoftwareX期刊主页：www.elsevier.com/locate/softxTravis J. Black和Alexei F. 切维亚科夫软件X 18（2022）1010512VVVV∈=× ×粉末床印刷工艺。在最简单的设置中，要填充的总域可以表示为平行六面体形砖块的并集，最初，砖块填充有球体的分布，然后用作构建块以达到所需的体积。砖块的对称性被用来跟踪粒子的接触。开发了两种相关的方法，并在当前基于Matlab的软件包中提供给每一个都由一个负责填充单个砖块的部分和一个构建体积的部分组成。对于方法1，在第2节中讨论了填充单元块的方法、任意平行六面体形状的体积、主球体拟合函数、典型程序序列和三个运行示例。第一种方法随机填充单元块的边缘，然后是面，然后是体积，这导致单元块填充有非对称球形分布。总域可以与单个砖块重合，或者由x，y和z方向上的几个砖块组成;在后一种情况下，总域由单元砖块围绕其面的反射构成，提供边界球接触。对于第二种方法，第3节提供了详细信息和示例。在方法2中，与方法1不同，对于单元砖，在每个方向上仅填充一个边缘，并且每个边缘制作四个平行副本。同样，每个平面中只有三个面填充了随机球体，并复制到相对的面上。这产生了具有相同的相对面（但是非对称体积填充）的单元块，并且因此可以使用单元块的直接复制-粘贴来填充更大的总域。方法2与方法1的另一个区别是，在方法2中，每个面上的球体的中心位于砖面本身上，而不是完全位于单元砖的内部并接触砖面。示例3（第4节）说明了几何方法的应用，以创建形状更复杂的域的球形填充：减去在x方向上两个相对面的中间居中的两个半球的平行六面体.在第5节中考虑了包含考虑粉末床3D打印过程中离散激光诱导的基于热的结合的简单模型的物理示例，参考文献10。[8]的一项建议。5.1节描述了物理原理和常数，5.2节给出了一个小正方形的模拟打印结果。在当前工作中使用的示例中，采用了球面半径的Weibull 和Gamma概率分布[9，10]。选择威布尔和伽马分布参数以对应于涉及钢球的粉末床增材制造。该软件支持Matlab中提供的所有概率分布。最后，在第二节中进行了总结讨论第六条.2. 用球体填充域：方法1第一种方法，以及稍后描述的第二种方法，可以用来填充具有给定随机分布的球体的任何平行六面体形状的域。该域可以完全随机的方式填充，或者为了更快的计算，可以将其细分为更小的标准平行六面体（“单位砖”）。在后一种情况下，单个砖块将随机填充球体，并且砖块可以在x，y和z方向上复制和连接，如下所述，任何规定的次数，以创建.在前一种情况下，整个域被视为单个单元块。在球形填充期间，记录随机球的位置和半径，以及接触球之间的成对连接，以及与位于V的每个面上的球相对应的球索引集。2.1. 方法一：填充单元砖用给定的球体分布填充单元砖的第一种方法的初始步骤包括在砖的每个角中放置平均半径的球体，然后用从相同分布随机抽取的具有尺寸的接触球体填充每个相邻角之间的边缘，然后以相同的方式填充完成所有面后，将填充剩余体积当放置新球体时，将记录与新球体接触的相邻球体。这样，当体积被填充时，有一个包含所有接触的球体对的列表。当放置新球体时，与给定球体接触的其他球体的列表被保持，基于恒定的无量纲参数ε，其指定要被认为接触的两个颗粒的可接受的分离/重叠。被认为是“接近”的粒子 当需要放置一个新球体时，为了确定它的位置，程序会遍历可能的“父粒子”列表，如下所述，新球体被放置成与可能的父粒子接触（当我们在3D中工作时，这个想法在图2D中显示）。①的人。这是通过求解方程组来实现的=Rn+R1，=Rn+R2，=Rn+R3，其中xiR3，i1， 2， 3是三元组父球体，xn是新球体中心的未知位置，Ri表示相应的半径。一旦新球体被放置在其假定的位置，将运行检查以查看它是否与任何其他球体重叠。如果是这样，则丢弃当前推定位置，并且将该球体与下一组可能的父球体进行匹配。一旦放置了球体并且不与任何其他球体相交，就会存储其位置和接触，并从给定的分布中随机绘制下一个球体如果给定的新球体与任何父球体都不匹配，则将丢弃该球体，并随机绘制新球体半径。我们注意到，在方法1中，与下面的方法2不同，单元砖的每个边和面都覆盖有不同的随机球体组。2.2. 方法1：填充平行六面体体积当单元砖被填充时，它可以“原样”使用以在相对少量的计算时间内，通过利用单元砖的对称性。采用体积填充的方法1来生成尺寸为m的平行六面体np，其中m、n和p分别是x、y和z方向上的单元砖的数量。在方法1中，关于与单元砖的面接触的球体的信息用于通过相对于其面平面对称地反射单元砖来填充期望的体积，使得面上的球体将与单元砖的对称副本中的重复该过程直到满足所需长度，然后在垂直方向上重复该过程。因此，代替用球体填充空间的耗时过程，简单地通过原始单元砖中的球体坐标的对称变换来计算附加砖中的球体坐标，并且直接复制半径和接触信息。Travis J. Black和Alexei F. 切维亚科夫软件X 18（2022）1010513•¯（）下一页×V•×ׯׯV≥λKλFig. 1. 与父球体接触的新球体：2D动画。2.3. 方法一：主球面拟合函数Matlab函数Method1GenerateSpheres.m实现了如上所述的单元砖生成和体积填充的方法1。输入和输出参数按出现顺序给出。输入参数：ProbabilityDistr：球形粒子半径的Matlab概率分布。这个对象可以使用Matlab的makedist函数创建对于给定的分布，我们称r为平均球半径。FaceGoal：与单位砖面接触的球体所覆盖的单位砖面面积的分数。BodyGoal：球体填充的单位砖体积的分数。• SphereContactParameter：接触参数ε，在[0，1]范围内。如果两个球形粒子UnitBrickNSpheres：每个单元砖中的球体数量。• 职位：矩阵x1· · ·xn· · ·y1· ··yn· · ·z1· · ·zn···第一行存储V中所有球体的x坐标，第二行存储所有球体的y坐标，第三行存储所有球体的z坐标。 因此，位置矩阵的第n列给出了第n个球体的坐标。半径：a 1FinalNSpheres矩阵，存储整个域中所有球体的半径。第n个条目是第n个球体的半径。Contacts：跟踪哪些粒子处于接触状态。此矩阵由两列组成;同一行是两个球体的索引对，球体接触参数r被认为是接触。·彼此，他们是[客户端]contact.• ListXmin、ListYmin、ListZmin、ListXmax、ListYmax、ParentParameter： 0， 1。如果两个球形粒子在ParentParameterr内，它们被认为是进一步粒子的潜在BrickSideLength：一个由三个值组成的数组，对应于单元砖边沿x、y和z的绝对长度（物理单位）。BrickNumbers：三个整数值的数组，指定构建总体积V所需的单元砖在x、y和z方向上的副本数量。函数Method1GenerateSpheres.m使用预定的概率分布来填充在x、y和z方向上具有物理尺寸的球体给出ListZmax：单列矩阵，每个矩阵存储所有索引位置矩阵中的球形粒子的与总域的相应边界接触对应于最小x，最小y，最小z，最大x，最大y，和最大z2.4. 方法一：典型程序序列和运行实例作为第一种体积填充方法的运行示例，我们选择Weibull分布[9，10]作为球体半径，由概率密度函数（PDF）f（r;λ，k）=k（r）k−1e−（r/λ）k，（2.2）Li=BrickSideLnumbers（1）×BrickNumbers（i），i=1，2，3.（2.1）输出参数：• FinalNSpheres：总域V中的球体的最终总数。其中r0是描述球体半径的维度随机变量，λ >0是以与随机变量r相同的长度单位测量的尺度参数，k>0是无量纲形状参数。分布（2.2）具有平均值r<$=λΓ（1+1），（ 2.3）······Travis J. Black和Alexei F. 切维亚科夫软件X 18（2022）1010514× ×=VVVV× ×=V其中，r是gamma函数。对于当前示例，我们选择对应于316L不锈钢粉末的粉末床3D打印过程的随机球体参数[9]，表1对于实施例1A和1B，对于由不同边长的221单位立方体构建的 总 域 V ， 对于球体接触参数的两个不同值，单位立方体中球体N的样本桌面计算时间T。λ= 15。7µ m，k= 3。五十五岁，14岁。14µm。（2.4）立方体边长球体接触参数0.1 0.2我们注意到，在文献中，特别是致力于增材制造的文献中，经常使用球形颗粒的直径而例如，基于直径的值λ31岁参考文献[9]中使用4µm（参见参考文献[9]表1，钢样品S2）。matlab脚本下文附录A中列出的示例1A_Method1_Generate_and_Plot.m规定了概率分布（2.2），并定义了运行的主要参数，包括要用球体填充的域的特性，以及控制球体大小（2.4）和表面和体积填充比的变量。附录A中还列出了脚本中使用的一些命令和运行参数的说明。该脚本调用主体积填充函数Method1GenerateSpheres.m，保存数据，并绘制生成的图形。实施例1A. 在方法1的第一个实施例运行中，使用以下输入参数• ProbabilityDistr：Weibull（2.2），（2.4）.• FaceGoal：0.8.• BodyGoal：0.55。• SphereContactParameter：0.2.• ParentParameter：0.5.• BrickSideLn：[1; 1; 1]*std_length，其中std_length=15 r<$.• BrickNumbers：[2; 2; 1].我们注意到，面部目标参数的参考值可以被计算为球体投影面积与边长为2r的正方形的面积的比率。类似地，身体目标被估计为球体与外接砖体积比的比率，其产生FaceGoal<$（πr2）/（2 r）2= π/4 <$0。78，（2.5）BodyGoal<$（4 πr 3/3）/（2 r）3 <$0. 52.可以通过实验确定特定应用程序的FaceGoal和BodyGoal球体生成脚本Example1A_Method1_Generate_Plot.m还用于生成图并保存当前示例的图形文件。图图2（a，c，e）示出了单元砖、其内部结构以及与实施例1A的给定分布（2.2）的概率密度相比的实际颗粒尺寸的直方图实施例1B. 在这里，我们使用与示例1A相同的设置，具有更大的单位砖边长：std_length=30r<$，因此，每单位立方体有四倍多的球体。这个例子的单位立方体和球体大小直方图如图所示。2（b、d、f）。特别地，图1B中的实际球体尺寸分布直方图。 2（f）比实施例1A更接近给定的威布尔分布（图1A）。2（e））;这是由于增加了将随机大小的球体拟合到单元比例1A中的立方体大（相对于r′）。15r'（实施例1A）T=6 min，N=450T= 2 min，N=41230r'（实施例1B）T=187分钟， N=3075T=40 min， N= 2867示例1A和1B：计算时间。计算在Matlab 2021 a上进行，使用戴尔工作站，配备两个至强处理器、32个逻辑处理器和128 GB内存。 下表1中列出的计算时间取决于单元砖的尺寸大小和单元砖的数量沿轴的BrickNumbers以形成总域.计算时间对球接触参数ε有很强的依赖性。(All计算还在基于Intel i7的笔记本电脑上进行了测试，该笔记本电脑具有一个物理处理器和四个逻辑处理器以及16 GB内存，平均计算时间延长了1.5到3倍）。我们还注意到，绘图可能需要相对较长的时间，类似于球体填充的计算时间，这是由于大量的球形粒子，每个球形粒子在Matlab球体绘图例程中由具有多个面的图形对象表示。实施例1C：非立方体单元砖。 在当前示例中，我们使用类似于示例1的Example1A_SpherePackingMethod1_Generate.m的函数来调用具有不同参数的主球体拟合函数Method1GenerateSpheres.m：ProbabilityDistr：球体半径的威布尔分布（2.2）、（2.4）• FaceGoal：0.8，• 身体目标：0.55，• 球形接触参数：0.2，• ParentParameter：0.5，• BrickSideLn： [1.2; 1.7; 1]*std_length，其中std_length=15r<$，• BrickNumbers：[4; 4; 2].因此，总域大小为44232个单位的砖是用球体填充的。每个单元砖都是非立方体的，尺寸长度在BrickSideLounge中指定。图4显示了最终的单元砖和球体接触对应于最小和最大x值的边，以及域的总构建。在相同的硬件/软件配置上，使用上面指定的32个单元砖对总域进行方法1球体填充的总计算花费了大约10分钟。为了比较，如果相同的体积由四个较大的砖构建，即，• BrickSideLn： [1.2; 1.7; 1]*std_length，其中std_length=30r<$，• BrickNumbers： [2; 2; 1]，计算时间增加到大约160分钟。3. 用球体填充域：方法2图图3示出了总体积V的构造，2× 2× 1单位立方体和连接偶的连通图位于水平板z=（0. 6 ± 0。3）×std_length。类似于第一种方法，第二种方法用于用一个或多个单元填充平行六面体形状的域V·Travis J. Black和Alexei F. 切维亚科夫软件X 18（2022）1010515¯图二. （a）实施例1A：边长为15 r的单位立方体，使用方法1用球体填充，威布尔分布（2.2），（2.4）。(c)中间1 / 3的球体立方体，显示填充立方体的内部结构。(e)所得样品中的实际球体半径分布，与给定的威布尔分布（2.2）、（2.4）进行比较。（b，d，f）相同的图例1B：边长为30r'的单位立方体。Travis J. Black和Alexei F. 切维亚科夫软件X 18（2022）1010516×V图3.第三章。（a）实施例1B：两个单元立方体连接形成2 ×1行顶视图. (b)两行（四个单位立方体）连接形成最终的域V。(c)实施例1B的单位立方体的水平板子域的连通性图（球体按比例缩小到其实际尺寸的40%以可视化连接）。砖块里装满了随机的球体，球体的半径遵循给定的概率分布。第二种方法与第一种方法的主要区别之一是，在方法2中，每个面上的球体的中心位于面本身上，而不是完全位于单元砖的内部并接触砖面（方法1）。这会导致单位砖块的外观略有不同。另一个区别在于边和面的对称性当填充单元砖时，方法2首先将相等的球体放置到所有角上;然后，与第一种方法不同，然后在每个（x、y和z）方向上填充，并且制作每个边缘的四个平行副本，以填充单元砖的所有12个边缘。类似地，只有三个面（在xy、yz和xz平面中）填充有随机球体，并且这些相被复制到相对的面上。最后，将单位砖体积随机填充球体。作为该过程的结果，获得具有对称面但非对称体积填充的单元砖为了在由多个单元砖组成的总域中构建球体填充，单元砖在x，y和z方向上直接复制所需的次数;当单元砖连接在一起时，面对称性提供了无缝配合3.1. 方法2：主球面拟合函数与 方 法 1 并 行 ， 方 法 2 在 Matlab 中 实 现 ，在 名 为Method2GenerateSpheres.m的函数中，该函数具有与方法1函数相同的参数（请参见第2.3节）。该函数使用任何规定的概率分布，填充单元砖，然后是从单元砖构造的域，如指定的。方法2球体中的所有输入和输出参数Travis J. Black和Alexei F. 切维亚科夫软件X 18（2022）1010517=≥K，见图4。（a）实施例1C：边长为15 r<$× [1. 二，一。7， 1]使用方法1和威布尔分布（2.2），（2.4）填充球体。(b)与砖边相切的球体对应于最小和最大x。（c）总域V：4× 4× 2个单元砖。拟合函数也与方法1的拟合函数一致，使得两种方法完全可互换，但导致不同的填充。3.2. 方法2：典型程序序列和运行示例实施例2. 在方法2的当前示例中，我们为球体大小选择了不同的概率分布：由概率密度函数（PDF）对于当前示例，我们选择与上面示例1类似的随机球体参数，以便平均半径近似匹配（2.4）θ= 7微米，k = 2，r<$= 14微米。（3.3） 脚本Example2_Method2_Generate_and_Plot.m（附录B）用于调用球体填充函数Method2GenerateSpheres. m，输入参数如下所示。f（r;k，θ）1Γ（k）xk−1e−x/θθ（3.1）• ProbabilityDistr：Gamma（3.1），（3.3）.• FaceGoal：1.0.其中r0是描述球体半径的维度随机变量，θ >0是以与随机变量r相同的单位测量的尺度参数，k>0是无量纲形状参数。分布（3.1）具有平均值• BodyGoal：0.9。• SphereContactParameter：0.1.• ParentParameter：0.5.• BrickSideLn：[1; 1; 1]*std_length，其中r<$=kθ。（第3.2节）Travis J. Black和Alexei F. 切维亚科夫软件X 18（2022）1010518std_length=30r<$（见式（3.3））。• BrickNumbers：[2;2;1] .Travis J. Black和Alexei F. 切维亚科夫软件X 18（2022）1010519¯× ×V××× ×=-V图五. （ a）实施例2：边长为30 r的单元砖（立方体），使用具有伽马分布（3.1）、（3.3）的方法2填充有球体（单元立方体被剪切到其尺寸，在面上示出球体的中心）。(b)与（a）相同，其中完整的单位立方体显示在较大的域中。(c)理论上的概率函数伽马分布（3.1）、（3.3）与单位立方体中实际颗粒尺寸的直方图。(d)顶部和底部的总体积V由221单位砖包含单位立方体的水平面的八个相同副本在这里，整个域由2 2 1个单位砖（在这个例子中是立方体）构成。我们注意到FaceGoal和BodyGoal的值比方法1示例中使用的值更高。这是自然的，因为对于方法2，边界球的中心是表2示例2的单位立方体中的球体N的数量为T的示例桌面计算时间，具有Gamma分布（3.1），（3.3），对于由2 2构建的总域V1个不同边长的单位砖，用于两个不同的球体接触参数值。位于表面上，这导致更高的相对面积，单位砖中球体所占的体积面部和身体目标的最佳值可以由用户通过实验来选择立方体边长15岁球体接触参数0.1 0.2T=7 s，N=160T= 11 s， N=178这取决于特定的应用。图 5显示了单位立方体与球体中心的面（比较图。 2（a））和总体积的z-边界由xy0平面中的单位立方体面的相同重复副本组成。表2包含与方法1中执行的计算类似并在表1中列出的计算的运行时间。有趣的是，对于方法2，与方法1不同（cf. 表1）中，较小球体接触参数0.1的计算时间与接触参数值0.2（当30rT=14 min，N= 1757T=15 min，N= 1771在方法2中使用实施例1的威布尔分布代替伽马分布）。我们还使用与方法1相同的威布尔分布（2.2）、（2.4）进行了方法2的运行（未显示结果图）。在这种情况下，使用方法2，计算时间比使用方法1完成的时间更快，Travis J. Black和Alexei F. 切维亚科夫软件X 18（2022）10105110V•¯V==×表3示例2的单位立方体中的球体N的数量为T的示例桌面计算时间，具有威布尔分布（2.2）、（2.4），对于由以下各项构建的总域V2×2×1不同边长的单位立方体，对于两个不同的实施例3A. 在第一个样本运行中，我们使用Weibull分布（2.2），（2.4）填充基于边长为立方体形状的域球体接触参数（参见方法1的表1）。std_length=30r<$（4.1）立方体边长球体接触参数0.1 0.2（对于威布尔分布，r由公式2.3给出），15岁T=4 min，N= 591T=4 min，N= 594相等半径的两个半球的牵引力。5× std_length。30rT=32 min， N= 3332T=38 min，N= 3606计算还导致单位立方体中更多的球体（比较方法1的表1和方法2的表34. 示例3：形状更复杂的域下一个示例构建在与方法1和2的球体填充，延伸到一个非平行六面体形状的域。在这个例子中，域类似于方法1和方法2中的单个单元砖，减去两个半球，两个半球的中心位于对应于x0和最大x的两个相对面的中间。这些球体的半径必须小于或等于砖尺寸最小值的1/2。为此示例创建了一个特殊的球体生成函数Example3GenerateSpheres.m。此函数的输入和输出参数概述如下。Example3GenerateSpheres.m的参数集小于方法1和方法2中球生成函数的参数集，但含义保持不变（参见第2.3节）。输入参数：ProbabilityDistr：球形粒子半径的Matlab概率分布。平均球体半径由r表示。BrickSideLength：一个由三个值组成的数组，对应于单元砖边沿x、y和z的绝对长度（物理单位）。HemisphereRadii：一个数组，包含两个值，对应于两个半球的半径（物理单位），中心位于两个相对面的中间，对应于x0和最大x，定义了计算域。FaceGoal：域边界面面积的分数，由与之接触的球体的投影• BodyGoal：球体填充的域体积的分数。• sphereContactParameter：接触参数，在脚本Example3A_Cube_Hemisph_Generate_and_Plot.m（附录C）使用以下参数调用球体放置例程Example3GenerateSpheres.m。• ProbabilityDistr：Weibull（2.2），（2.4）.• BrickSideLtd.：[1; 1; 1]* 标准长度。• 半球半径：[0.5; 0.5]* 标准长度。• FaceGoal：0.4.• BodyGoal：0.4。• SphereContactParameter：0.1.本例在工作站配置中花费49分钟完成（在笔记本电脑配置中花费138分钟）。由此产生的球形排列如图所示。第6（a）段。实 施 例 3B. 在 第 二 个 样 本 运 行 中 ， 使 用 相 同 的 Weibull 分 布（2.2）、（2.4）和典型域大小（4.1）以填充基于非立方块的域 ， 其 中 减 去 两 个 半 径 不 相 等 的 半 球 。 基 于 示 例 3A 的Example3A_Brick_Hemisph_Generate_and_Plot.m的脚本调用具有不同域和计算参数的球体放置例程，如下所示。• ProbabilityDistr：Weibull（2.2），（2.4）.• BrickSideLtd.：[1.3; 1; 1]*std_length.• 半球半径：[0.2; 0.4]* 标准长度。• FaceGoal：0.4.• BodyGoal：0.4。• SphereContactParameter：0.2.由此产生的球形填充和相关图表如图所示。6（b、c、d）。在工作站配置中完成此计算需要41分钟（在笔记本电脑配置中为144分钟）。得到的球体尺寸分布直方图（图 6（d））显示出与给定的概率密度函数很好的一致性。5. 一个物理实例：粉末床3D打印过程我们现在使用球体填充方法2来模拟增材制造过程中的基于热的结合，该增材制造过程使用金属粉末床和引导的激光束来加热球形颗粒，从而在它们之间形成结合，如参考文献[8]中所述。316L不锈钢粉末的粒度分布近似于Weibull分布（2.2），（2.4）[9，10]。[0， 1]。 如果两个球面粒子在¯SphereContactParameterr，它们被认为是接触的。输出参数：• NSpheres：放置到域V中的随机球体的总数。5.1. 热源和粉床模型对于热源模型，物理假设如下。从激光到第i个球形粒子的热通量由下式给出：位置：随机球体的坐标矩阵（见2.3节）。qi激光器ri3=Qr3，• 半径：一个向量，存储了全域V中所有随机球体的半径（见2.3节）。接触：包含接触的球体对的矩阵（见2.3节）。其中Q是激光的总功率，ri是颗粒的半径，并且ri是激光束的半径• 对流产生的热通量为qiconv=kb（TR−Ti），·····ℓ·Travis J. Black和Alexei F. 切维亚科夫软件X 18（2022）10105111·=+∑不我我IJ不J我我我miCp（）下一页见图6。 （a）实施例3A：边长（4.1）减去两个相等半球的单位立方体，填充随机球体和威布尔分布（2.2），（2.4）。(b)例如图3B：一个非立方域减去两个不相等的半球的类似计算。(c)示例的计算域和第一个父球体父对象3B. （d）实施例3B中的实际球体半径分布，与给定的威布尔分布相比其中kb是传热系数，TR是周围空气的温度。• 两个粒子之间的热通量由下式给出：利用离散时间步长t1，t2，t1，. . . ，第i个粒子的温度更新表示为：QtQ =kT− T其中kt是导热系数，Ti、Tj表示粒子i和j的温度。进入第i个粒子的总热通量因此变为Nqi=qi激光器 +qiconv +qi j。j=1Tt+Tt =Tt +我不知道，其中Tt是前一个时间步的粒子qt是总的初始能流，mi是第i个物体的质量球形颗粒，Cp表示材料比热。如果颗粒i和j接触，并且都高于烧结温度Ts，则在它们之间形成结合·、Travis J. Black和Alexei F. 切维亚科夫软件X 18（2022）10105112我·==×=×·=∼ ××·=·见图7。一个正方形的模拟印刷品（第5节）。(a)打印域（大小以µm为单位）。(b)激光束路径。（c）激光的第二次通过（温度标度，单位为开尔文）。（d）激光的最后一次通过使用以下参数的模拟打印的以下样本值[9]的文件。• 激光器总功率：Q=100 W。• 球体接触参数：0.1，• ParentParameter：0.5，• BrickSideLn： [1; 1; 0.12]* std_length，其中• 空气的导热系数：τ= 0。262 W/（m· K）。std_length=30r<$，一钢的导热系数：τs15 W/（m K）。空气边界处的传热系数：k bπτaravg/2 2。9090 10−6 W/ K。钢的传热系数：kt πτs ravg/23。3309 10−4 W/ K。• 钢的比热容C p= 0。5 J/（g· K）。• 钢密度ρ=8× 10−12 g/（µm）3。• 球形粒子的质量m i=（4π/3）ρr 3。• 激光半径r=50µ m。• 环境空气温度：TR=300 K。• 烧结温度：Ts=1000 K。5.2. 一个正方形的模拟印刷品使用方法2（第3节）生成样品粉末床，参数• ProbabilityDistr：Weibull（2.2），（2.4），• FaceGoal：1.0，• BodyGoal：0.9，• BrickNumbers：[6; 6; 1].计算大约花了五分钟， 2545 254550块89µ m，平均颗粒半径14。14µ m（图7（a））。选择图7（b）中所示的激光路径来模拟在半导体层中间的小正方形的印刷。通过在粒子之间建立热诱导连接来控制域。图图7（c，d）示出了第二次和最后一次激光通过时的温度图。图8描绘了在模拟印刷过程的结果中结合的颗粒和结合图的俯视图。当前的软件对于采用离散方法的模型和应用程序是有用的，其中离散颗粒之间的相互作用（可以通过球体建模）变得重要。对于这样的应用，感兴趣的三维域需要填充有从规定的概率分布随机抽取的尺寸的大量接触的非重叠球体。目前的代码采用几何原理来放置初始球体，然后插入新的球体，以便它们与先前的球体以随机方式接触。用户指定的参数包括域维度，Travis J. Black和Alexei F. 切维亚科夫软件X 18（2022）10105113∼∼图8.第八条。 一个正方形的 模 拟 打 印 （第5节）。（a）在模拟打印过程中结合的颗粒（b）颗粒之间形成的结合所需的概率分布，以及所需的表面和球在域的表面和本体中的体积分数。输出包括球体中心位置、球体半径和由彼此接触的球体对定义的连通性图。使用标准的Matlab绘图功能可以很容易地可视化所实例包括在内。6. 讨论和影响目前的手稿描述了一个Matlab实现的几何为基础的算法填充一个三维域与球的半径随机分布根据一个给定的概率分布。两种方法，方法1和方法2（第2节）使用一个或多个标准“单元砖”的副本（由用户指定的尺寸形成），用球体这两种方法是实现为可自由互换的Matlab函数，具有相同的输入和输出参数。不同的方法在单位域的面上的球体的放置原则和对称性。包括使用各种球半径分布和不同域形状的例子。所提出的代码的输出包括球心位置、球半径和由球对定义的连通图相互接触。使用标准的Matlab绘图功能可以很容易地可视化所得到的球体填充。第4节中的一个例子说明了涉及半球边界的更复杂域形状的球形填充技术。所提出的软件将是有用的模型和applications，采用离散的方法，其中大量的随机尺寸的离散颗粒之间的相互作用，建模的球体，变得重要。应用包括增材制造（粉末床3D打印）中离散金属部件之间的热致连接模型，这是这项工作的主要动机（见第5节）。可直接受益于使用的其他应用领域该软件的各种模型的过滤，渗流和扩散的粒状介质中出现，例如，在地质科学、生命科学、过滤相关的工业应用和废物管理。在这项工作中提出的代码可以有效地适应生成不同尺寸和形状的域的球形填充物，使用Matlab中可用的任何常见粒度概率分布。所提出的算法在域方面的通用性形状和尺寸以及颗粒尺寸概率分布使得这些算法与生成球体（或其它对象）的随机或非随机紧密堆积的那些算法相当不同。具有相同的尺寸（参见，例如，[6、7]）。后一种设计的最优性的一个常用度量是重复物体所占的体积分数φ。众所周知，全同球体的无规密堆积可以达到φ 0。64，而面心立方晶格对应于φ 0。74（[6]和其中的参考文献）。 然而，在目前的工作中提出的代码，主要是针对球半径不相等的问题，下面是一个给定的非奇异（通常是连续的）概率分布模型的一个特定的物理情况。因此，所提供的代码可以与射击方法结合使用，以估计在任何给定的设置中由球形颗粒占据的接近最大体积分数，而不需要执行实际的实验。它也可以用来近似解决在直接测量不可行的情况下估计颗粒尺寸分布的逆问题。竞合利益作者声明，他们没有已知的竞争性财务利益或个人关系，可能会影响本文报告的工作致谢作者感谢匿名推荐人提供的宝贵建议和参考，并感谢加拿大NSERC通过发现基金RGPIN-2019-05570和USRA奖学金提供的支持。Travis J. Black和Alexei F. 切维亚科夫软件X 18（2022）10105114∗∗∗∗附录A.球体填充方法1的数据生成和绘图脚本，示例1A• 主文件：Example1A_Method1_Generate_and_Plot.m• 输出量：- 菲格 2（b、d、f）。– 图与图相似。图3（a，b）为实施例1B。– Matlab数据文件Example1A_Method1_Results.mat• 所需的其他文件（在同一个文件夹中）：– Method1GenerateSpheres.m：主要的方法1球体填充例程– Method1GenerateSpheres.m的辅助例程文件：M1位置2Xmax.m∗M1Position2Xmin.mM1位置2Ymax.m∗M1位置2Ymin.m∗M1Position2Zmax.m∗M1Position2Zmin.m∗M1位置3.m搜索M1Search2D.mM1Search3D.m脚本中使用的一些命令和运行参数的说明如下。• rng（0）：将Matlab伪随机种子设置为零（默认值）。• Weibull_scale，Weibull_shape：公式2.4中的分布参数λ，k。• ProbabilityDistribution： Matlab分布 可变 为 的 PDF（2.2），（2.4）.¯• cube_side_length：单位