CN107729667A - 基于四维空间理论的三维空间球体填充算法 - Google Patents

Abstract

本发明公开了一种基于四维空间理论的三维空间球体填充算法，该方法依次通过随机生成松散的四维超球体颗粒、超球体颗粒在重力作用下坠落并被三维容器捕捉、可选择地对超球体颗粒在四维空间中做振动处理和通过三维容器与超球体颗粒切割最终生成三维颗粒模型实现；该方法基于物理方法的基础上通过引入多维空间来克服传统算法的一些缺点，使堆积过程不受三维几何形态的约束，并且还适用于对已有的颗粒模型进行二次填充和动态调整，实现对任意形态容器的填充。

Description

基于四维空间理论的三维空间球体填充算法

技术领域

本发明涉及球体填充模拟算法技术领域，特别涉及一种基于四维空间理论的三维空间球体填充算法。

背景技术

如何使用大量的球体对固定空间进行填充(如将乒乓球装入木箱)以求得最大的填充密实度称为球体填充问题。球体填充问题看上去虽然简单，却是一个典型的跨学科问题，在数学、物理化学、材料科学、生物医学、冶金工程、机械制造及岩土工程等多个领域都有着广泛的应用背景，因此一直是相关领域研究的热点、重点问题之一。对球体填充问题的研究最早可追溯到1611年开普勒提出的关于等径球体最密填充的猜想，但目前为止，关于球体填充的研究依然还不完善，例如关于非等径球体混合填充的研究尚不全面，各类球体填充算法有待完善等。

上世纪六十年代以来，研究人员开发了各种球体随机填充算法，各有优缺点。目前基于球体随机填充方法主要可分为两大类，一类是基于几何的方法，这类算法主要侧重于球体堆积的结果；另一类是基于物理的方法，此类算法在计算机模拟中主要侧重于球体堆积的过程，也称之为动态模拟法。基于几何的颗粒排布方法有很多，大体分为序列添加和集合重排两类。集合重排有三种经典的方法，Jodrey和Tory的算法(JT算法)、Clarke和Wiley的算法(CW算法)、He方法。基于物理的方法大都基于分子动力学和离散单元法。1990年Lubachevsky和Stillinger提出了基于分子动力学的球体堆积算法用来模拟等径球体的堆积过程，后来Kansal等学者在该算法的基础上进行修改使该算法可以模拟非等径球体的堆积过程。近年来研究人员将离散单元法用于球体随机堆积问题的研究，离散单元法可考虑颗粒的平运和旋动，在颗粒的运动过程中则可考虑接触力、范德华力、重力、摩擦力等各类作用力的影响。在前人研究成果基础上，本申请提出一种基于空间四维理论对三维空间进行球体填充的算法。

发明内容

本发明的目的是提供一种在物理方法的基础上通过引入多维空间来克服传统算法，基于四维空间理论的三维空间球体填充算法。

传统基于重力堆积法对容器口部分颗粒的排布不平整及对含有复杂孔洞形态容器的颗粒排布困难等方面。为了解决这些问题，本申请提出的方法通过在第四个维度进行重力堆积，实现对三维空间的球体堆积。如图1所示，当将四个维度分别用x轴、y轴、z轴、zu轴表示。如果要用二维图完全表现四维空间，则需幅二维图；如果要用三维图表现四维空间，则需幅三维图，即4幅三维图即可完全表现出四维空间超球体颗粒的形态。其中，三维视图中超球体颗粒的灰度与其在三视图垂直的另一维度的位置相关联，以第一视图(Z-Y-ZU视图)为例进行说明：该视图表现出了超球体颗粒在Z-Y-ZU三维视角下的形态，超球体颗粒在垂直于Z-Y-ZU三维视图的另一维度──x轴的位置则由其在Z-Y-ZU三维视图中的色彩来表示。将超球体颗粒的度量化为色彩空间，系统指定x轴的范围为x₁到x₂，则超球体颗粒灰度值与其所在x维度的坐标位置相对应且呈现出线性变化的关系。由于对于三维空间形体在第四维视角看，其内部结构都是暴露在外的，因此堆积过程可以不受三维几何形态的约束。

为此，本发明技术方案如下：

一种基于四维空间理论的三维空间球体填充算法，其中：

情况一：当三维球体颗粒在三维容器中的填充密度＞0.5时，三维球体颗粒空间排布算法包括如下步骤：

S1、在计算机中模拟给定三维容器的三维结构图形，并在其上方构造四维超球体的生成区域；设定四维超球体的球径和数量，在所述生成区域随机生成松散排布的四维超球体颗粒，并获取每个四维超球体颗粒的球心坐标和球体半径；

S2、使在所述生成区域随机生成松散排布的四维超球体颗粒在第四维重力作用下发生坠落并被下方三维容器捕捉，当三维容器的填充密度达到0.5时四维超球体颗粒停止下落，计算获取此时被三维容器捕捉的各个四维超球体颗粒的最终球心坐标和球体半径；

S3、对被三维容器捕捉的全部四维超球体颗粒施加一个随机向上的力，使全部四维超球体颗粒在四维空间中发生振动而重新排布，同时原四维超球体的生成区域内的四维超球体颗粒继续下落，直至被三维容器捕捉的全部四维超球体颗粒的填充密度达到预设值，计算获取被三维容器捕捉的各个四维超球体颗粒的填充密度达到预设值时的最终球心坐标和球体半径；

S4、根据被三维容器捕捉的全部四维超球体颗粒的最终位置坐标得出三维容器与超球体颗粒最终生成三维球体颗粒模型，并获取该三维球体颗粒模型中各个三维球体颗粒的球心坐标和球体半径；

情况二：当三维球体颗粒在三维容器中的填充密度≤0.5时，三维球体颗粒空间排布算法包括如下步骤：

S1、在计算机中模拟给定三维容器的三维结构图形，并在其上方构造四维超球体的生成区域；设定四维超球体的球径和数量，在所述生成区域随机生成松散排布的四维超球体颗粒，并获取每个四维超球体颗粒的球心坐标和球体半径；

S2、使在所述生成区域随机生成松散排布的四维超球体颗粒全部在第四维重力作用下发生坠落并被下方三维容器捕捉并计算获取当三维容器的填充密度达到预设填充密度时被三维容器捕捉的各个四维超球体颗粒的最终球心坐标和球体半径；

S4、根据被三维容器捕捉的全部四维超球体颗粒的最终位置坐标得出三维容器与超球体颗粒最终生成三维颗粒模型，并获取该三维颗粒模型中各个三维球体颗粒的球心坐标和球体半径。

具体而言，该算法的基本原理如图2所示，在四维空间中的三维容器在第四个维度看，可表示为四维超平面上的平面区域。本申请提出的颗粒排布基本思想是让在四维空间中随机分布的超球体颗粒在第四维方向重力的作用下向三维容器所在的超平面坠落，容器对超球体颗粒的质心具有粘结作用，这样容器就可以捕捉落入其内的超球体颗粒，随着不断的球体颗粒坠落，在重力势能的作用下三维容器内的颗粒将会越来越密实。最终，通过对捕捉三维容器中的超球体颗粒进行超平面切割操作来生成对应的三维颗粒模型。

进一步地，步骤S1的具体实施方法如下：

S101、在计算机中绘制出给定的三维容器的三维结构图形，并在第四维度以容器为第四维的原点，在其上方构造一个四维包盒：[x₁，x₂，y₁，y₂，z₁，z₂，zu₁ ^4D，zu₂ ^4D]，其中，zu₁ ^4D为包盒第四维度底部坐标,zu₂ ^4D为包盒第四维度顶部坐标；

S102、对用于填充三维容器的四维超球体颗粒的参数进行设定，参数由两个数列[R_i]和[N_i]来控制，其中，R_i代表的是第i组球体颗粒的半径，N_i代表的是对应第i组球体颗粒的数目；

S103、以步骤S102设定的最小超球体颗粒的直径将步骤S101构造的四维包盒划分成一个个四维单元体，然后计算机不断随机选取任意一个四维单元体的中心点为球心生成四维超球体颗粒，直至完成将步骤S102设定的全部四维超球体颗粒在四维包盒内的稀疏生成；同时获取在四维包盒内生成的每个四维超球体颗粒的球心坐标和球体半径。

进一步地，步骤S2的具体实施方法如下：

S201、采用有限元网格划分的方法将步骤S1绘制的三维容器的三维结构图形的外表面划分为N个三角形以表征三维容器的表面轮廓，然后通过获取每个三角形的顶点坐标，并将在四维包盒内生成的每个四维超球体颗粒的球心坐标与三角形顶点坐标进行比较，判断四维超球体颗粒在第一、第二、第三维度的坐标(x，y，z)是否在N个三角形顶点坐标构成的三维容器外表面轮廓内，筛选出在第四维度重力作用下自四维包盒内下落的全部四维超球体颗粒中经过下方三维容器并被捕捉到的填充颗粒；

其中，引入有限元网格划分的方法对三维容器进行表征实际上是首先三维结构图形的外表面划分为N个三角形，然后获取所有三角形的所有顶点坐标，用三角形的顶点坐标表征划分的N个三角形的信息，进一步通过N个三角形的信息表征了三维容器的外表面轮廓；这样，通过量化就可以实现判断四维超球体颗粒的在第四维度的下落轨迹是否经过三维容器，即判断出自四维包盒内下落的全部四维超球体颗粒中哪些是进入到三维容器内的填充颗粒。

S202、对经步骤S201筛选出的填充颗粒按照设定的间隔时间计算出经过每段间隔时间Δt后位于三维容器内的各个四维超球体颗粒的实时球心坐标和速度；其中，具体计算公式如下：

某四维超球体颗粒i受到的合力F：F：F＝m_ig+F_AB+∑F^c+∑F^d，其中，m_i为i球的质量，g为重力加速度，F_AB为i球受到的粘结力，∑F^c为i球受到的其它超球体颗粒的总接触力，∑F^d为i球受到的其它超球体颗粒的总阻尼力；

i球受到的粘结力F_AB：F_AB＝∣min(Zu·k_h,0)∣·n_AB，其中，Zu表示在下落过程中四维超球体颗粒的球心在第四维度的实时坐标，n_AB为从i球质心指向超平面的单位向量，k_h为设定的粘结系数；

任意其它超球体颗粒w对i的接触力F_iw ^c：其中，Δδ为接触变形量，k_δ为刚度系数；n_wi为从超球体w的中心点指向超球体i中心点的单位向量；

任意其它超球体颗粒w对i的阻尼力F_iw ^d计算为：其中，ξ为阻尼系数，v_wi为超球体颗粒i对w的相对速度；

然后，根据设定的间隔时间Δt通过迭代的方法进行计算出每颗四维超球体颗粒的实时位置和速度，计算公式如下：

v^(t+Δt/2)＝v^(t-Δt/2)+a·Δt，其中，v^(t+Δt/2)为超球体颗粒在t+Δt/2时刻的速度，v^(t-Δt/2)为超球体颗粒在t-Δt/2时刻的速度；

s^(t+Δt)＝s^(t)+v^(t+Δt/2)·Δt，其中，s^(t+Δt)为超球体颗粒在t+Δt时刻的位移，s^(t)为超球体颗粒在t时刻的位移；

在计算出各个四维超球体颗粒的实时位置和速度的同时，计算出当前状态下三维容器的填充密度；当预设的填充密度≤0.5时，上述计算过程在三维容器的填充密度达到预设填充密度时结束；当预设的填充密度＞0.5时，上述计算过程在三维容器的填充密度达到0.5时结束。

进一步地，步骤S3的具体实施方法如下：

S301、对位于三维容器中的全部超球体颗粒在第四维度(zu轴表示)施加力一个竖直向上的振动力F'_zu：F'_zu(t)＝max(0,ζA_zuSin(ω_zu(t-t₀)))，其中，A_zu和ω_zu分别是振动幅度和振动频率，t₀表示在zu方向上的振动开始时间，t表示在zu方向上的振动结束时间，ζ是一个0到1的随机数；

S302、对受到振动力F'_zu的超球体颗粒和自四维包盒内新下落的超球体颗粒按照步骤S202设定的间隔时间Δt计算出经过每段间隔时间Δt后，位于三维容器内的各个四维超球体颗粒的实时球心坐标和速度。

进一步地，步骤S4的具体实施方法如下：设任意被三维容器捕捉的超球体颗粒质点坐标为(x'，y'，z'，zu')，通过公式：将三维容器与四维超球体切割并生成三维颗粒模型，其中，R^3D为三维球体颗粒半径，R^4D为四维超球体颗粒半径，zu'为超球体颗粒第四维度坐标；最终得到该任意三维球体颗粒模型的坐标为(x'，y'，z')，半径为R^3D。

与传统的计算机模拟球体堆积相比较，该基于四维空间理论的三维空间球体填充算法中球体颗粒在堆积过程中不同尺寸三维球体颗粒的填充比例基本不会发生变化，在保证与设定值基本一致的基础上实现了较高的堆积密度；此外，本专利所提出的堆积方法便于对不规则或具有孔洞状态的三维填充区域进行填充，对传统的基于三维重力球体排布方法是一个有力的补充；此外，传统基于物理模拟球体堆积法受限于三维空间，得到的堆积结果往往是局部最优解而不是全局最优解，基于四维空间理论的三维球体颗粒排布方法能够突破三维空间的几何障碍限制，更加利于找到全局最优解。

附图说明

图1为多颗四维超球体颗粒堆积状态的四视图；

图2为本发明的基于四维空间的三维颗粒排布算法原理示意图；

图3为本发明的基于四维空间理论的三维空间球体填充算法的流程图；

图4为步骤S1的初始稀疏超球体颗粒生成过程的流程图；

图5为步骤S1的初始稀疏四维超球体生成过程的二维原理示意图；

图6为步骤S2的多颗四维超球体颗粒受重力下降过程的四视图；

图7为步骤S2的多颗四维超球体颗粒填充的过程的示意图；

图8为步骤S2的多颗四维超球体颗粒的填充计算过程的流程图。

具体实施方式

下面结合附图及具体实施例对本发明做进一步的说明，但下述实施例绝非对本发明有任何限制。

采用该基于四维空间理论的三维空间球体填充算法对多组不同粒径的三维球体颗粒在一个给定容积的三维被填充容器内按照预设的填充密度进行排布的排布方式进行模拟并获得合理的排布方式。

如图3所示，以预设的填充密度＞0.5的情况为例进行说明。其中，填充密度的定义为三维容器内全部超球体颗粒的总体积与三维容器容积的比值。

该填充算法具体包括如下步骤：

S1、在计算机中模拟给定三维容器的三维结构图形，并在其上方构造四维超球体的生成区域；设定四维超球体的球径和数量，在所述生成区域随机生成松散排布的四维超球体颗粒，并获取每个四维超球体颗粒的球心坐标和球体半径；

S2、使在所述生成区域随机生成松散排布的四维超球体颗粒在第四维重力作用下发生坠落并被下方三维容器捕捉，并计算获取当三维容器的填充密度达到0.5时被三维容器捕捉的各个四维超球体颗粒的最终球心坐标和球体半径；

S3、对被三维容器捕捉的全部四维超球体颗粒施加一个随机向上的力，使全部四维超球体颗粒在四维空间中发生振动而重新排布，同时原四维超球体的生成区域内的四维超球体颗粒继续下落，直至被三维容器捕捉的全部四维超球体颗粒的填充密度达到预设值，计算出被三维容器捕捉的各个四维超球体颗粒的填充密度达到预设值时的最终球心坐标和球体半径；

S4、根据被三维容器捕捉的全部四维超球体颗粒的最终位置坐标得出三维容器与超球体颗粒最终生成三维球体颗粒模型，并获取该三维球体颗粒模型中各个三维球体颗粒的球心坐标和球体半径。

具体地，步骤S1的具体实现过程如下：

S101、在计算机中绘制出给定的三维容器的三维结构图形，并在第四维度以容器为第四维的原点，在其上方构造一个四维包盒：[x₁，x₂，y₁，y₂，z₁，z₂，zu₁ ^4D，zu₂ ^4D]，其中，zu₁ ^4D为包盒第四维度底部坐标,zu₂ ^4D为包盒第四维度顶部坐标，那么，包盒沿第四维度的厚度即为zu₂ ^4D-zu₁ ^4D；

S102、根据对三维容器的填充要求，对用于填充三维容器的四维超球体颗粒的参数进行设定，参数由两个数列[R_i]和[N_i]来控制，其中，R_i代表的是第i组球体颗粒的半径，N_i代表的是对应第i组球体颗粒的数目；通过这两组参数的设定实现步骤S2中填充在给定三维容器中的四维超球体颗粒按照预设颗粒级配(即颗粒粒径的大小及其在总体中所占的比例)模型的生成；

S103、以步骤S102设定的最小超球体颗粒的直径(即min(2R_i))将步骤S101构造的四维包盒划分成一个个四维单元体，然后计算机不断随机选取任意一个四维单元体的中心点为球心生成四维超球体颗粒，直至完成将步骤S102设定的全部四维超球体颗粒在四维包盒内的稀疏生成；同时获取在四维包盒内生成的每个四维超球体颗粒的球心坐标。

其中，在步骤S103中，计算机不断随机选取四维单元体的中心点生成超球体，当某个四维单元体的中心点被四维超球体颗粒占据时，该中心点称为“占据点”，将被系统在后续的随机选点中自动排除，不再被计算机选取生成超球体颗粒；同时，系统不允许生成的超球体之间存在重合，因此当四维单元体中心点处生成的超球体颗粒因周围超球体颗粒的存在被限制，即最大尺寸(该中心点到最近一颗超球体颗粒的距离)小于min(R_i)时，该中心点称为“抑制点”，也会被系统舍弃。图5所示为超球体颗粒在四维包盒内生成的二维原理示意图，图中的“●”代表未生成超球体颗粒的四维单元体的中心，代表“占据点”，代表“抑制点”。通过此过程最终完成在四维超球体颗粒在定义范围内的稀疏生成，其具体流程如图5所示。四维超球体颗粒最终生成的半径和数量符合用户定义参数[R_i]和[N_i]。需要说明的是，为了保证能够生成足够数量的四维超球体，因此四维包盒坐标范围：[x₁，x₂，y₁，y₂，z₁，z₂，zu₁ ^4D，zu₂ ^4D]中，第四维度的范围zu₁ ^4D～zu₂ ^4D之间的距离足够大，能够保证所有四维超球体颗粒的生成。

步骤S2的具体实现过程如下：

S201、采用有限元网格划分的方法将步骤S1绘制的三维容器的三维结构图形的外表面划分为N个三角形以表征三维容器的表面轮廓，然后通过获取每个三角形的顶点坐标，并将在四维包盒内生成的每个四维超球体颗粒的球心坐标与三角形顶点坐标进行比较，判断四维超球体颗粒在第一、第二、第三维度的坐标(x，y，z)是否在N个三角形顶点坐标构成的三维容器外表面轮廓内，筛选出在第四维度重力作用下自四维包盒内下落的全部四维超球体颗粒中经过下方三维容器并被捕捉到的填充颗粒；

具体来说，有限元网格划分的方法即通过将三维被填充容器的外表面划分成N个三角形；当划分的三角形的数量足够多时，即可实现用N个三角形表征三维被填充区域外表面，而每个三角形均由三个顶点连线形成，因此通过获得所有三角形的顶点坐标即可实现对容器的表征，同时也实现了将待填充的三维容器坐标化。进一步即可通过判断四维超球体颗粒在第一、第二、第三维度的坐标(x，y，z)是否在N个三角形顶点坐标构成的三维容器外表面轮廓内，筛选出能够被下方三维容器捕捉到的填充颗粒。

如图6所示为生成区域内四维超球体颗粒下降过程中的四视图。超球体颗粒在第四维度重力的作用下垂直三维容器所在的超平面纷纷下降。由于在第四维度，待填充的三维容器视为一个平面，即上文提到的超平面，为保证所有四维填充颗粒能够粘结在超平面上，实现对三维容器的填充，同时便于对四维超球体的后续受力进行分析，因此设定三维容器所在区域对填充颗粒的质心均具有粘结作用，即任意超球体颗粒在第四维重力作用下降落，如果落点在三维容器所在的范围内，三维容器所在的超平面对四维超球体具有粘结力F_AB，其对应公式为：F_AB＝∣min(Zu·k_h,0)∣·n_AB，其中，Zu表示在下落过程中四维超球体颗粒的球心在第四维度的实时坐标，n_AB为从i球质心指向超平面的单位向量，k_h为设定的粘结系数，其具体值一般设定在10⁷～10⁸N/m的范围内；如图7所示，在粘结力F_AB和第四维重力的平衡作用下超球体颗粒最终会在三维空间视图下将处于“悬浮”状态，这样三维容器的填充区域即通过粘结力完成对四维超球体颗粒的“捕捉”，逐渐实现对指定三维空间的填充。

具体填充过程流程图如图8所示。

S202、对经步骤S201筛选出的填充颗粒按照设定的间隔时间计算出经过每段间隔时间后位于三维容器内的各个四维超球体颗粒的实时位置和速度，并计算出当前状态下三维容器的填充密度。

由于在球体被“捕捉”进入填充空间的过程中，填充空间内的超球体不断发生接触、碰撞、移动，超球体的位置不断更新变化最终达到静态的力学平衡。假定所有超球体颗粒为刚体，那么对于任意超球体颗粒i来说，其与其它超球体颗粒之间的作用力包括接触力F_c和阻尼力F_d；因此，假设存在n个超球体颗粒对超球体颗粒i有力学作用，对接触力F_c、阻尼力F_d进行n个求和就是n个超球体颗粒对超球体颗粒i的总力；具体地，

i球的运动规律满足牛顿第二定律，即：F＝m_ia，其中，F为i球受到的总力，m_i为i球的质量，a为i球的加速度；

进一步对i球受力F进行展开，F＝m_ig+F_AB+∑F^c+∑F^d＝F_x+F_y+F_z+F_zu，其中，g为重力加速度，F_AB为i球受到的粘结力，∑F^c为i球受到的其它超球体颗粒的总接触力，∑F^d为i球受到的其它超球体颗粒的总阻尼力，F_x、F_y、F_z和F_zu分别为i球受到的总力分别在x轴、y轴、z轴和zu轴上的分力；

对应地，加速度a为：a＝F/m_i＝a_x+a_y+a_z+a_zu，在四维空间中可以分解成以下四个方程：

在与超球体颗粒接触、有力学作用的n个超球体中，我们任意选取一个其它的超球体颗粒w对i球的力学作用进行分析。

设超球体颗粒i中心点坐标为(x_i，y_i，z_i，zu_i)，超球体颗粒w中心点坐标为(x_w，y_w，z_w，zu_w)，因此，

(1)超球体颗粒w对i的接触力F_iw ^c计算为：其中，Δδ为接触变形量；k_δ为刚度系数；n_wi为从超球体w的中心点指向超球体i中心点的单位向量；

其中，接触变形量Δδ的计算公式为：

其中，R_i为超球体颗粒i的半径，R_w为超球体颗粒w的半径；

(2)超球体颗粒w对i的阻尼力F_iw ^d计算为：其中，ξ为阻尼系数，v_wi为超球体颗粒i对w的相对速度；

在上述对接触力F_iw ^c和F_iw ^d的计算中，Δδ、k_δ和ξ的实际取值均与实际填充在三维容器内的颗粒的实际性质相关。

由于四维超球体颗粒为独立的单元体，每个球体的受力情况和运动速度、位移等信息可以使用迭代的方法进行计算，具体计算过程：

首先，设定更新三维容器内的超球体颗粒实时位置和速度的间隔时间，即时间步长Δt；

时间步长Δt基于以下两个假设进行设定：假设1：Δt足够小，因此在一个Δt范围内，颗粒只和与其接触的颗粒存在力的关系，其他颗粒对该颗粒没有影响；假设2：在所选取的Δt范围内，颗粒的速度和加速度保持不变；

因此，在每个时间步长内根据公式：v^(t+Δt/2)＝v^(t-Δt/2)+a·Δt和公式：s^(t+Δt)＝s^(t)+v^(t+Δt/2)·Δt完成对计算单元体的速度和位移的更新，其中，v^(t+Δt/2)为超球体颗粒在t+Δt/2时刻的速度，v^(t-Δt/2)为超球体颗粒在t-Δt/2时刻的速度，s^(t+Δt)为超球体颗粒在t+Δt时刻的位移，s^(t)为超球体颗粒在t时刻的位移。

当预设的填充密度≤0.5时，这一计算过程在三维容器的填充密度达到预设填充密度时结束；当预设的填充密度＞0.5时，这一计算过程在三维容器的填充密度达到0.5时结束。

步骤S3的具体实现过程如下：

四维超球体颗粒在重力的作用下下降，对系统定义的三维容器进行填充。当用户定义三维容器的被填充密度≤0.5时，该填充过程可以通过三维容器对超球体颗粒的“捕捉”直接实现；而当用户定义三维容器的被填充密度＞0.5时，则需要经过步骤S3的振动处理，使三维容器内球体分布更加均匀、球体堆积状态更加稳定，以此提高填充密度。

具体地，对位于三维容器中的全部超球体颗粒在第四维度(zu轴表示)施加力一个竖直向上的振动力F'_zu，作用时间为t-t₀，对颗粒振动进行模拟。在振动模拟过程中，第四维度的粘结力、重力依然存在，因此超球体颗粒在第四维重力作用下正常坠落。对颗粒体进行模拟振动，通过超球体颗粒之间相互碰撞产生力的作用，超球体颗粒在力的作用下产生新的速度和位移，位置不断进行更新直至振动结束，在粘结力F_AB与重力的共同作用下三维容器内的超球体颗粒最终回到静态平衡状态。

这一过程的具体计算方法如下：

S301、对位于三维容器中的全部超球体颗粒在第四维度(zu轴表示)施加力一个竖直向上的振动力F'_zu：F'_zu(t)＝max(0,ζA_zuSin(ω_zu(t-t₀)))，其中，A_zu和ω_zu分别是振动幅度和振动频率，t₀表示在zu方向上的振动开始时间，t表示在zu方向上的振动结束时间，ζ是一个0到1的随机数；

S302、对受到振动力F'_zu的超球体颗粒和自四维包盒内新下落的超球体颗粒按照步骤S202设定的间隔时间Δt计算出经过每段间隔时间Δt后，位于三维容器内的各个四维超球体颗粒的实时球心坐标和速度。

具体计算方法如下：对于原三维容器内的全部四维超球体颗粒来说，其在受振动力作用时间段t₀～t内受到的合力F₁:F₁＝m_ig+F_AB+F′zu∑F^c+∑F^d；由于受到振动力，原三维容器内的四维超球体颗粒向上运动，而当振动力消失时，这些四维超球体颗粒的受力发生改变，所受到的合力变为F₂：F₂＝m_ig+F_AB+∑F^c+∑F^d；同时，原四维超球体的生成区域内的四维超球体颗粒继续下落，这些四维超球体颗粒所受到的合力同样可以用F₂的计算公式：F₂＝m_ig+F_AB+∑F^c+∑F^d进行计算；按照步骤S202设定的间隔时间Δt计算出经过每段间隔时间Δt后，位于三维容器内的各个四维超球体颗粒的实时球心坐标和速度，同时计算出每段间隔时间Δt后被三维容器捕捉的各个四维超球体颗粒的填充密度。

如果一次振动依然不能达到预设的的填充密度时，则重新对三维容器内的所有四维超球体颗粒一个新的振动力F'_zu，即通过重复步骤S302以对四维超球体颗粒进行重新排布，直到达到预设的填充密度。由于这一过程实现了超球体颗粒在振动过程中以第四维度为路径重新在三维容器中确定自身的合理位置，突破了三维振动过程中的几何限制，为得到最优堆积模型提供了可能。

步骤S4的具体实现过程如下：

三维容器与四维超球体颗粒进行切割生成三维球体颗粒模型即为最终的三维填充球体颗粒模型。

具体地，设任意被三维容器的粘结力F_AB捕捉的超球体颗粒质点坐标为(x'，y'，z'，zu')，通过公式：

将三维容器与四维超球体切割并生成三维颗粒模型，其中，R^3D为三维球体颗粒半径，R^4D为四维超球体颗粒半径，zu'为超球体颗粒第四维度坐标。最终，切割生成的该任意三维球体颗粒模型的坐标为(x'，y'，z')，半径为R^3D。

需要说明的是，当得到的四维球体颗粒半径小于给定颗粒最小半径(即min(R_i))时，则将该四维球体颗粒排除在最终模型中。

由于该四维填充算法从更高一个维度进行重力填充的优点在于不受三维几何的约束，因此除实现上述填充过程以外，还可以对已有的颗粒模型进行二次填充，加入更细小的颗粒，或者对模型进行动态调整，实现对任意形态的容器填充。

Claims (5)

1.一种基于四维空间理论的三维空间球体填充算法，其特征在于，

i)当三维球体颗粒在三维容器中的填充密度＞0.5时，三维球体颗粒空间排布算法包括如下步骤：

S1、在计算机中模拟给定三维容器的三维结构图形，并在其上方构造四维超球体的生成区域；设定四维超球体的球径和数量，在所述生成区域随机生成松散排布的四维超球体颗粒，并获取每个四维超球体颗粒的球心坐标和球体半径；

S2、使在所述生成区域随机生成松散排布的四维超球体颗粒全部在第四维重力作用下发生坠落并被下方三维容器捕捉，并计算获取当三维容器的填充密度达到0.5时被三维容器捕捉的各个四维超球体颗粒的最终球心坐标和球体半径；

S3、对被三维容器捕捉的全部四维超球体颗粒施加一个随机向上的力，使四维超球体颗粒在四维空间中发生振动而重新排布，同时原四维超球体的生成区域内的四维超球体颗粒继续下落，直至被三维容器捕捉的全部四维超球体颗粒的填充密度达到预设值，计算获取被三维容器捕捉的全部四维超球体颗粒的填充密度达到预设值时的球心坐标和球体半径；

S4、根据被三维容器捕捉的全部四维超球体颗粒的最终位置坐标得出三维容器与超球体颗粒最终生成三维球体颗粒模型，并获取该三维球体颗粒模型中各个三维球体颗粒的球心坐标和球体半径；

ii)当三维球体颗粒在三维容器中的填充密度≤0.5时，三维球体颗粒空间排布算法仅需进行上述步骤的S1、S2和S4即可，其中在步骤S2中，计算获取当三维容器的填充密度达到预设填充密度时被三维容器捕捉的各个四维超球体颗粒的最终球心坐标和球体半径。

2.根据权利要求1所述的基于四维空间理论的三维空间球体填充算法，其特征在于，所述步骤S1的具体步骤为：

S101、在计算机中绘制出给定的三维容器的三维结构图形，并在第四维度以容器为第四维的原点，在其上方构造一个四维包盒：[x₁，x₂，y₁，y₂，z₁，z₂，zu₁ ^4D，zu₂ ^4D]，其中，zu₁ ^4D为包盒第四维度底部坐标,zu₂ ^4D为包盒第四维度顶部坐标；

S102、对用于填充三维容器的四维超球体颗粒的参数进行设定，参数由两个数列[R_i]和[N_i]来控制，其中，R_i代表的是第i组球体颗粒的半径，N_i代表的是对应第i组球体颗粒的数目；

S103、以步骤S102设定的最小超球体颗粒的直径将步骤S101构造的四维包盒划分成一个个四维单元体，然后计算机不断随机选取任意一个四维单元体的中心点为球心生成四维超球体颗粒，直至完成将步骤S102设定的全部四维超球体颗粒在四维包盒内的稀疏生成；获取在四维包盒内生成的每个四维超球体颗粒的球心坐标和球体半径。

3.根据权利要求1所述的基于四维空间理论的三维空间球体填充算法，其特征在于，所述步骤S2的具体步骤为：

S201、采用有限元网格划分的方法将步骤S1绘制的三维容器的三维结构图形的外表面划分为N个三角形以表征三维容器的表面轮廓，然后将获取每个三角形的顶点坐标，并将在四维包盒内生成的每个四维超球体颗粒的球心坐标与三角形顶点坐标进行比较，判断四维超球体颗粒在第一、第二、第三维度的坐标(x，y，z)是否在N个三角形顶点坐标构成的三维容器外表面轮廓内，筛选出在第四维度重力作用下自四维包盒内下落的全部四维超球体颗粒中经过下方三维容器并被捕捉到的填充颗粒；

S202、对经步骤S201筛选出的填充颗粒按照设定的间隔时间Δt计算出经过每段间隔时间Δt后位于三维容器内的各个四维超球体颗粒的实时球心坐标和速度，计算公式如下：

某四维超球体颗粒i受到的合力F：F＝m_ig+F_AB+∑F^c+∑F^d，其中，m_i为i球的质量，g为重力加速度，F_AB为i球受到的粘结力，∑F^c为i球受到的其它超球体颗粒的总接触力，∑F^d为i球受到的其它超球体颗粒的总阻尼力；

i球受到的粘结力F_AB：F_AB＝∣min(Zu·k_h,0)∣·n_AB，其中，Zu表示在下落过程中四维超球体颗粒的球心在第四维度的实时坐标，n_AB为从i球质心指向超平面的单位向量，k_h为设定的粘结系数；

任意其它超球体颗粒w对i的接触力F_iw ^c：其中，Δδ为接触变形量，k_δ为刚度系数；n_wi为从超球体w的中心点指向超球体i中心点的单位向量；

任意其它超球体颗粒w对i的阻尼力F_iw ^d计算为：其中，ξ为阻尼系数，v_wi为超球体颗粒i对w的相对速度；

然后，根据设定的间隔时间Δt通过迭代的方法进行计算出每颗四维超球体颗粒的实时位置和速度，计算公式如下：

v^(t+Δt/2)＝v^(t-Δt/2)+a·Δt，其中，v^(t+Δt/2)为超球体颗粒在t+Δt/2时刻的速度，v^(t-Δt/2)超球体颗粒在t-Δt/2时刻的速度；

s^(t+Δt)＝s^(t)+v^(t+Δt/2)·Δt，其中，s^(t+Δt)为超球体颗粒在t+Δt时刻的位移，s^(t)为超球体颗粒在t时刻的位移；

在计算出各个四维超球体颗粒的实时位置和速度的同时，计算出当前状态下三维容器的填充密度；当预设的填充密度≤0.5时，上述计算过程在三维容器的填充密度达到预设填充密度时结束；当预设的填充密度＞0.5时，上述计算过程在三维容器的填充密度达到0.5时结束。

4.根据权利要求3所述的基于四维空间理论的三维空间球体填充算法，其特征在于，所述步骤S3的具体步骤为：

S301、对位于三维容器中的全部超球体颗粒在第四维度(zu轴表示)施加一个竖直向上的振动力F'_zu：F'_zu(t)＝max(0,ζA_zuSin(ω_zu(t-t₀)))，其中，A_zu和ω_zu分别是振动幅度和振动频率，t₀表示在zu方向上的振动开始时间，t表示在zu方向上的振动结束时间，ζ是一个0到1的随机数；

S302、对受到振动力F'_zu的超球体颗粒和自四维包盒内新下落的超球体颗粒按照步骤S202设定的间隔时间Δt计算出经过每段间隔时间Δt后，位于三维容器内的各个四维超球体颗粒的实时球心坐标和速度。

5.根据权利要求1所述的基于四维空间理论的三维空间球体填充算法，其特征在于，所述步骤S4的具体步骤为：设任意被三维容器捕捉的超球体颗粒质点坐标为(x'，y'，z'，zu')，通过公式：将三维容器与四维超球体切割并生成三维颗粒模型，其中，R^3D为三维球体颗粒半径，R^4D为四维超球体颗粒半径，zu'为超球体颗粒第四维度坐标；最终得到该任意三维球体颗粒模型的坐标为(x'，y'，z')，半径为R^3D。