CN110414116A - 一种颗粒材料的颗粒状态分析方法、装置及设备 - Google Patents

Abstract

一种颗粒材料的颗粒状态分析方法包括：获取待分析的颗粒材料中的颗粒的几何形态，通过傅里叶级数表征所述颗粒材料的数学描述，所述数学描述包括颗粒形状描述符和颗粒位置描述符；根据所述颗粒材料的数学描述，判断颗粒之间是否接触；如果颗粒之间有接触，则根据所述颗粒材料的数学描述，计算颗粒之间的接触特征；根据所述颗粒之间的接触特征，计算颗粒之间的接触力，根据所计算的接触力更新颗粒的运动状态，得到颗粒材料的宏观力学响应。通过该数学描述方式，能够适用于任意形状颗粒，而且颗粒形状描述符保持常量，独立于颗粒的位置描述符，在颗粒发生移动或旋转时，不需要更新颗粒的形状描述符，运算量较低，计算效率较高。

Description

一种颗粒材料的颗粒状态分析方法、装置及设备

技术领域

本申请属于材料分析领域，尤其涉及一种颗粒材料的颗粒状态分析方法、装置及设备。

背景技术

离散元是一种基于颗粒的数值模型，特别适用于描述块状颗粒材料的力学行为。它最初是由Cundall(康德尔)和Strack(斯特拉克)提出的，用于岩土材料的分析。此后，离散元被广泛应用于各种颗粒材料的建模，并模拟各种不同工程分支中的固体处理、粉末流动等各种问题。离散元对颗粒材料中的所有颗粒进行显式建模，直接捕获颗粒之间的相互作用力，并跟踪每个颗粒的运动。颗粒材料的宏观力学行为即表现为所有组成颗粒的相互作用和运动的集合。

通常，离散元模型中有两种类型的基本元素：颗粒和边界。基本元素通常被假定为是刚性的，但可以相互重叠。颗粒是具有闭合表面的物体，它可以用一个简单的几何形状(例如，球体或椭球)来表示，也可以用多个简单几何形状的组合来表示。在现有的离散元方法中，颗粒只能通过简单的几何形状来描述，对任意不规则颗粒不具有普适性。虽然采用简单几何形体堆叠的方式可一定程度描述颗粒的不规则形状，但精度不高，而且导致运算量较大。

发明内容

有鉴于此，本申请实施例提供了一种颗粒材料的颗粒状态分析方法、装置及设备，以解决现有技术中采用简单几何形体堆叠的方式的精度不高，而且运算量较大的问题。

本申请实施例的第一方面提供了一种颗粒材料的颗粒状态分析方法，所述颗粒材料的颗粒状态分析方法包括：

获取待分析的颗粒材料中的颗粒的几何形态，通过傅里叶级数表征所述颗粒材料的数学描述，所述数学描述包括颗粒形状描述符和颗粒位置描述符；

根据所述颗粒材料的数学描述，判断颗粒之间是否接触；

如果颗粒之间有接触，则根据所述颗粒材料的数学描述，计算颗粒之间的接触特征；

根据所述颗粒之间的接触特征，计算颗粒之间的接触力，根据所计算的接触力更新颗粒的运动状态；

根据所获取的颗粒当前的状态以及颗粒在下一时刻的状态，通过时间步的迭代，获取颗粒在一时间段的力学响应与状态演化。

结合第一方面，在第一方面的第一种可能实现方式中，所述获取待分析的颗粒材料中的颗粒的几何形态，通过傅里叶级数表征所述颗粒材料的数学描述，所述数学描述包括颗粒形状描述符和颗粒位置描述符的步骤包括：

在局部坐标系中，颗粒表面上的任何点通过其相对于原点的极半径r和相对于x轴的极角θ来确定，将颗粒表面上的任何点A的局部坐标(x_A，y_A)可表示为x_a＝rcosθ,y_a＝rsinθ；

将极半径r通过傅里叶展开的方式表示为： 其中，a₀，a_n，b_n为颗粒形态描述符，n为傅里叶级数；

在全局坐标系中，颗粒表面上的任何点位置表示为X_A＝X_C+r(θ)cosθ，Y_A＝Y_C+r(θ)sinθ，其中，X_C、Y_C为颗粒形心位置的颗粒位置描述符。

结合第一方面，在第一方面的第二种可能实现方式中，所述根据所述颗粒材料的数学描述，判断颗粒之间是否接触的步骤包括：

对于任意两个颗粒p，q，其颗粒形心位置分别为C_p、C_q，点A_p是颗粒p上的任意一点，向量C_pA_p的极角为v，直线C_qA_p与颗粒q的表面相交于A_q，点Ap与点Aq的距离表示为其中L表示其下标两点之间的距离；

计算所述距离D的最小值，如果所述最小值小于或等于0，则颗粒p、q相接触。

结合第一方面的第二种可能实现方式，在第一方面的第三种可能实现方式中，所述计算所述距离D的最小值的步骤包括：

计算线段C_qA_p的长度：其中，点A_p的坐标表示为： 角标p、q表示颗粒p和q，为傅里叶基函数；

计算线段C_qA_q的长度：其中，点Ap的坐标表示为：

根据线段CqAp和线段CqAq的长度，计算距离D为：

通过牛顿迭代法计算所述距离的最小值。

结合第一方面，在第一方面的第四种可能实现方式中，所述根据所述颗粒材料的数学描述，计算颗粒之间的接触特征的步骤包括：

根据计算的距离确定两个颗粒的表面的接触点；

根据所述接触点确定所述接触点对应的接触中心；

根据所述接触中心确定接触分支向量，所述接触分支向量即为接触中心指向颗粒形心的向量；

获取所述接触点的法向和切向，根据通过所述接触中心、方向为接触法向的直线与两个颗粒表面轮廓的交点，确定接触重叠点；

根据所述接触重叠点之间的距离计算接触重叠量；

计算接触相对切向位移变化量Δδ_t＝v_tΔt，其中， 为相对切向速度，v为颗粒的速度，w为颗粒的角速度，ln为接触分支向量长度，上标p、q指示颗粒p和q，颗粒的状态包括位置、倾角、速度和角速度，颗粒的速度、角速度由颗粒当前的状态决定，。

结合第一方面，在第一方面的第五种可能实现方式中，所述根据所述颗粒之间的接触特征，计算颗粒之间的接触力，根据所计算的接触力更新颗粒的运动状态的步骤包括：

根据公式计算接触力和接触弯矩：F_n＝k_nδ_n，其中F_n为法向接触力，F_t为切向接触力；M_c＝0接触弯矩，δ_n为法向位移量，k_n、k_t为接触刚度，μ_c为接触摩擦系数；

根据所述接触力，结合颗粒重力和外部荷载，计算颗粒所受全力，根据所受全力计算颗粒的加速度、速度、移动位置。

本申请实施例的第二方面提供了一种颗粒材料的颗粒状态分析装置，所述颗粒材料的颗粒状态分析装置包括：

数学描述单元，用于获取待分析的颗粒材料中的颗粒的几何形态，通过傅里叶级数表征所述颗粒材料的数学描述，所述数学描述包括颗粒形状描述符和颗粒位置描述符；

接触判断单元，用于根据所述颗粒材料的数学描述，判断颗粒之间是否接触；

接触特征计算单元，用于如果颗粒之间有接触，则根据所述颗粒材料的数学描述，计算颗粒之间的接触特征；

运动状态计算单元，用于根据所述颗粒之间的接触特征，计算颗粒之间的接触力，根据所计算的接触力更新颗粒的运动状态，得到颗粒材料的宏观力学响应。

结合第二方面，在第二方面的第一种可能实现方式中，所述数学描述单元包括：

表面点表示子单元，用于在局部坐标系中，颗粒表面上的任何点通过其相对于原点的极半径r和相对于x轴的极角θ来确定，将颗粒表面上的任何点A的局部坐标(x_A，y_A)可表示为x_a＝rcosθ,y_a＝rsinθ；

颗粒形态描述符获取子单元，用于将极半径r通过傅里叶展开的方式表示为：其中，a₀，a_n，b_n为颗粒形态描述符，n为傅里叶级数；

颗粒位置描述符获取子单元，用于在全局坐标系中，颗粒表面上的任何点位置表示为X_A＝X_C+r(θ)cosθ，Y_A＝Y_C+r(θ)sinθ，其中，X_C、Y_C为颗粒形心位置的颗粒位置描述符。

本申请实施例的第三方面提供了一种颗粒材料的颗粒状态分析设备，包括存储器、处理器以及存储在所述存储器中并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现如第一方面任一项所述颗粒材料的颗粒状态分析方法的步骤。

本申请实施例的第四方面提供了一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序被处理器执行时实现如第一方面任一项所述颗粒材料的颗粒状态分析方法的步骤。

本申请实施例与现有技术相比存在的有益效果是：获取颗粒材料中的颗粒的几何形态后，通过傅里叶级数表征颗粒材料的数学描述，根据所述数学描述判断颗粒是否接触，如果接触后则根据数学描述确定颗粒之间的接触特征，并根据所述接触特征计算颗粒之间的接触力，根据计算的接触力更新颗粒的运动状态，从而得到颗粒材料的宏观力学响应。由于本申请能够直接对不规则粒子进行建模，能够适用于任意形状颗粒，而且颗粒形状描述符保持常量，独立于颗粒的位置描述符，在颗粒发生移动或旋转时，不需要更新颗粒的形状描述符，相对于简单几何锥体堆叠的方式，特别是在高精度颗粒形态的条件下，运算量较低，计算效率较高。

附图说明

为了更清楚地说明本申请实施例中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本申请的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1是本申请实施例提供的一种颗粒材料的颗粒状态分析方法的实现流程示意图；

图2是本申请实施例提供的一种确定颗粒的数学描述方法的实现流程示意图；

图3是本申请实施例提供的一种判断颗粒是否接触方法的实现流程示意图；

图4是本申请实施例提供的一种判断颗粒是否接触的示意图；

图5是本申请实施例提供的一种计算接触特征方法的实现流程示意图；

图6是本申请实施例提供的接触特征的示意图；

图7是本申请实施例提供的一种确定颗粒材料的宏观力学响应的实现流程示意图；

图8是本申请实施例提供的一种颗粒材料的颗粒状态分析装置的示意图；

图9是本申请实施例提供的一种颗粒材料的颗粒状态分析设备的示意图。

具体实施方式

以下描述中，为了说明而不是为了限定，提出了诸如特定系统结构、技术之类的具体细节，以便透彻理解本申请实施例。然而，本领域的技术人员应当清楚，在没有这些具体细节的其它实施例中也可以实现本申请。在其它情况中，省略对众所周知的系统、装置、电路以及方法的详细说明，以免不必要的细节妨碍本申请的描述。

为了说明本申请所述的技术方案，下面通过具体实施例来进行说明。

图1为本申请实施例提供的一种颗粒材料的颗粒状态分析方法的实现流程示意图，详述如下：

在步骤S101中，获取待分析的颗粒材料中的颗粒的几何形态，通过傅里叶级数表征所述颗粒材料的数学描述，所述数学描述包括颗粒形状描述符和颗粒位置描述符；

根据待分析的颗粒材料，可以通过测量的方式，或者通过机器识别的方式，获取所述颗粒材料中的颗粒的几何形态，所述几何形态可以包括颗粒表面上的点的极半径、以及颗粒表面上的点相对于形心的位置等。获取所述颗粒材料的颗粒的数学描述的步骤可以如图2所示，包括：

在步骤S201中，在局部坐标系中，颗粒表面上的任何点通过其相对于原点的极半径r和相对于x轴的极角θ来确定，将颗粒表面上的任何点A的局部坐标(x_A，y_A)可表示为x_a＝rcosθ,y_a＝rsinθ；

所述局部坐标系，以颗粒的中心为坐标原点，颗粒的旋转、平移等操作都是围绕局部坐标系进行的，这时，当颗粒模型进行旋转或平移等操作时，局部坐标系也执行相应的旋转或平移操作。

在步骤S202中，将极半径r通过傅里叶展开的方式表示为： 其中，a₀，a_n，b_n为颗粒形态描述符，n为傅里叶级数；

其中，r(θ)也可表示为：r(θ)＝T_FS(θ)*C_FS，其中， 为傅里叶基函数，C_FS＝[a₀，a₁，b₁，…，b_N，b_N]为傅里叶系数。

在步骤S203中，在全局坐标系中，颗粒表面上的任何点位置表示为X_A＝X_C+r(θ)cosθ，Y_A＝Y_C+r(θ)sinθ，其中，X_C、Y_C为颗粒形心位置的颗粒位置描述符。

在步骤S102中，根据所述颗粒材料的数学描述，判断颗粒之间是否接触；

可以在两个颗粒的表面任取一个点，判断两个点之间的距离是否小于或等于零来判断两个颗粒是否接触，具体可以如图3所示，包括：

在步骤S301中，对于任意两个颗粒p，q，其颗粒形心位置分别为C_p、C_q，点A_p是颗粒p上的任意一点，向量C_pA_p的极角为v，直线C_qA_p与颗粒q的表面相交于A_q，点A_p与点A_q的距离表示为其中L表示其下标两点之间的距离；

在步骤S302中，计算所述距离D的最小值，如果所述最小值小于或等于0，则颗粒p、q相接触。

其中，计算所述距离D可以分别两个颗粒表面上的点与形心位置的距离，并通过求导的方式确定最小距离，具体可以包括：

计算线段C_qA_p的长度：其中，点Ap的坐标表示为： 角标p、q表示颗粒p和q，为傅里叶基函数，θ₀为颗粒初始倾角。

计算线段C_qA_q的长度：其中，点A_p的坐标表示为：

根据线段C_qA_p和线段C_qA_q的长度，计算距离D为：

通过牛顿迭代法计算所述距离的最小值。

其中，所述距离D的最小值可以通过获得，其可以通过牛顿迭代方法求解。牛顿迭代计算时，关键的导数计算式为：

如图4所示，当线段C_qA_p和线段C_qA_q的差值大于零时，颗粒没有发生接触，当线段C_qA_p和线段C_qA_q的差值等于零或小于零时，则表示颗粒之间有发生接触。

在步骤S103中，如果颗粒之间有接触，则根据所述颗粒材料的数学描述，计算颗粒之间的接触特征；

在判断颗粒之间有接触后，可以根据所述颗粒材料的数学描述，计算所述颗粒之间的接触特征，具体可以如图5所示，包括：

在步骤S501中，根据计算的距离确定两个颗粒的表面的接触点；

可以在计算距离D＝0时得到两个颗粒p、q表面的接触点I₁和I₂。

在步骤S502中，根据所述接触点确定所述接触点对应的接触中心；

接触点I₁和I₂的中点，即可确定接触中心M。

在步骤S503中，根据所述接触中心确定接触分支向量，所述接触分支向量即为接触中心指向颗粒形心的向量；

所述接触分支向量为接触中心M指向颗粒形心C_p和C_q的向量，即MC_p、MC_q。

在步骤S504中，获取所述接触点的法向和切向，根据通过所述接触中心、方向为接触法向的直线与两个颗粒表面轮廓的交点，确定接触重叠点；

如图6所示，计算接触法向n，即为接触点I₁和I₂所表示的向量I₁I₂的垂直方向，计算接触切向t，即为接触点I₁和I₂所表示的向量I₁I₂的方向，计算接触重叠点N₁和N₂，即为通过接触中心M、方向为接触法向n的直线与颗粒p、q的表面轮廓的交点。

在步骤S505中，根据所述接触重叠点之间的距离计算接触重叠量；

计算接触重叠量δn，即为重叠点N₁和N₂之间的距离。

在步骤S506中，计算接触相对切向位移变化量Δδ_t＝v_tΔt，其中， 为相对切向速度，v颗粒的速度，w为颗粒的角速度，l_n为接触分支向量长度，上标p、q指示颗粒p和q，且颗粒的状态包括位置、倾角、速度和角速度，颗粒的速度、角速度由颗粒当前的状态决定。

在步骤S104中，根据所述颗粒之间的接触特征，计算颗粒之间的接触力，根据所计算的接触力更新颗粒的运动状态。

得到颗粒材料的宏观力学响应的具体步骤可以如图7所示，包括：

在步骤S701中，根据公式计算接触力和接触弯矩：F_n＝k_nδ_n， 其中F_n为法向接触力，F_t为切向接触力；M_c＝0接触弯矩，δ_n为法向位移量，k_n、k_t为接触刚度，μ_c为接触摩擦系数；

在步骤S702中，根据所述接触力，结合颗粒重力和外部荷载，计算颗粒所受全力，根据所受全力计算颗粒的加速度、速度、移动位置；

在步骤S105中，根据所获取的颗粒当前的状态以及颗粒在下一时刻的状态，通过时间步的迭代，获取颗粒在一时间段的力学响应与状态演化。

其中，在每个颗粒的加速度、速度、移动位置确定后，可以根据所有颗粒状态，通过时间步的迭代计算，得到材料的宏观力学响应，从而能够有效的反应颗粒材料的真实反应，便于工作人员在施工过程中对固体处理、粉末流动等问题的处理。

应理解，上述实施例中各步骤的序号的大小并不意味着执行顺序的先后，各过程的执行顺序应以其功能和内在逻辑确定，而不应对本申请实施例的实施过程构成任何限定。

图8为本申请实施例提供的一种颗粒材料的颗粒状态分析装置的结构示意图，详述如下：

所述颗粒材料的颗粒状态分析装置包括：

数学描述单元801，用于获取待分析的颗粒材料中的颗粒的几何形态，通过傅里叶级数表征所述颗粒材料的数学描述，所述数学描述包括颗粒形状描述符和颗粒位置描述符；

接触判断单元802，用于根据所述颗粒材料的数学描述，判断颗粒之间是否接触；

接触特征计算单元803，用于如果颗粒之间有接触，则根据所述颗粒材料的数学描述，计算颗粒之间的接触特征；

运动状态计算单元804，用于根据所述颗粒之间的接触特征，计算颗粒之间的接触力，根据所计算的接触力更新颗粒的运动状态；

响应获取单元805，用于根据所获取的颗粒当前的状态以及颗粒在下一时刻的状态，通过时间步的迭代，获取颗粒在一时间段的力学响应与状态演化。

优选的，所述数学描述单元包括：

表面点表示子单元，用于在局部坐标系中，颗粒表面上的任何点通过其相对于原点的极半径r和相对于x轴的极角θ来确定，将颗粒表面上的任何点A的局部坐标(x_A，y_A)可表示为x_a＝rcosθ,y_a＝rsinθ；

颗粒形态描述符获取子单元，用于将极半径r通过傅里叶展开的方式表示为：其中，a₀，a_n，b_n为颗粒形态描述符，n为傅里叶级数；

颗粒位置描述符获取子单元，用于在全局坐标系中，颗粒表面上的任何点位置表示为X_A＝X_C+r(θ)cosθ，Y_A＝Y_C+r(θ)sinθ，其中，X_C、Y_C为颗粒形心位置的颗粒位置描述符。

图8所述颗粒材料的颗粒状态分析装置，与图1所述的颗粒材料的颗粒状态分析方法对应。

图9是本申请一实施例提供的颗粒材料的颗粒状态分析设备的示意图。如图9所示，该实施例的颗粒材料的颗粒状态分析设备9包括：处理器90、存储器91以及存储在所述存储器91中并可在所述处理器90上运行的计算机程序92，例如颗粒材料的颗粒状态分析程序。所述处理器90执行所述计算机程序92时实现上述各个颗粒材料的颗粒状态分析方法实施例中的步骤，例如图1所示的步骤101至103。或者，所述处理器90执行所述计算机程序92时实现上述各装置实施例中各模块/单元的功能，例如图5所示模块501至503的功能。

示例性的，所述计算机程序92可以被分割成一个或多个模块/单元，所述一个或者多个模块/单元被存储在所述存储器91中，并由所述处理器90执行，以完成本申请。所述一个或多个模块/单元可以是能够完成特定功能的一系列计算机程序指令段，该指令段用于描述所述计算机程序92在所述颗粒材料的颗粒状态分析设备9中的执行过程。例如，所述计算机程序92可以被分割成，各单元具体功能如下：

数学描述单元，用于获取待分析的颗粒材料中的颗粒的几何形态，通过傅里叶级数表征所述颗粒材料的数学描述，所述数学描述包括颗粒形状描述符和颗粒位置描述符；

接触判断单元，用于根据所述颗粒材料的数学描述，判断颗粒之间是否接触；

接触特征计算单元，用于如果颗粒之间有接触，则根据所述颗粒材料的数学描述，计算颗粒之间的接触特征；

运动状态计算单元，用于根据所述颗粒之间的接触特征，计算颗粒之间的接触力，根据所计算的接触力更新颗粒的运动状态，得到颗粒材料的宏观力学响应。

所述颗粒材料的颗粒状态分析设备9可以是桌上型计算机、笔记本、掌上电脑及云端服务器等计算设备。所述颗粒材料的颗粒状态分析设备可包括，但不仅限于，处理器90、存储器91。本领域技术人员可以理解，图9仅仅是颗粒材料的颗粒状态分析设备9的示例，并不构成对颗粒材料的颗粒状态分析设备9的限定，可以包括比图示更多或更少的部件，或者组合某些部件，或者不同的部件，例如所述颗粒材料的颗粒状态分析设备还可以包括输入输出设备、网络接入设备、总线等。

所称处理器90可以是中央处理单元(Central Processing Unit，CPU)，还可以是其他通用处理器、数字信号处理器(Digital Signal Processor，DSP)、专用集成电路(Application Specific Integrated Circuit，ASIC)、现成可编程门阵列(Field-Programmable Gate Array，FPGA)或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件等。通用处理器可以是微处理器或者该处理器也可以是任何常规的处理器等。

所述存储器91可以是所述颗粒材料的颗粒状态分析设备9的内部存储单元，例如颗粒材料的颗粒状态分析设备9的硬盘或内存。所述存储器91也可以是所述颗粒材料的颗粒状态分析设备9的外部存储设备，例如所述颗粒材料的颗粒状态分析设备9上配备的插接式硬盘，智能存储卡(Smart Media Card,SMC)，安全数字(Secure Digital,SD)卡，闪存卡(Flash Card)等。进一步地，所述存储器91还可以既包括所述颗粒材料的颗粒状态分析设备9的内部存储单元也包括外部存储设备。所述存储器91用于存储所述计算机程序以及所述颗粒材料的颗粒状态分析设备所需的其他程序和数据。所述存储器91还可以用于暂时地存储已经输出或者将要输出的数据。

所属领域的技术人员可以清楚地了解到，为了描述的方便和简洁，仅以上述各功能单元、模块的划分进行举例说明，实际应用中，可以根据需要而将上述功能分配由不同的功能单元、模块完成，即将所述装置的内部结构划分成不同的功能单元或模块，以完成以上描述的全部或者部分功能。实施例中的各功能单元、模块可以集成在一个处理单元中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个单元中，上述集成的单元既可以采用硬件的形式实现，也可以采用软件功能单元的形式实现。另外，各功能单元、模块的具体名称也只是为了便于相互区分，并不用于限制本申请的保护范围。上述系统中单元、模块的具体工作过程，可以参考前述方法实施例中的对应过程，在此不再赘述。

在上述实施例中，对各个实施例的描述都各有侧重，某个实施例中没有详述或记载的部分，可以参见其它实施例的相关描述。

本领域普通技术人员可以意识到，结合本文中所公开的实施例描述的各示例的单元及算法步骤，能够以电子硬件、或者计算机软件和电子硬件的结合来实现。这些功能究竟以硬件还是软件方式来执行，取决于技术方案的特定应用和设计约束条件。专业技术人员可以对每个特定的应用来使用不同方法来实现所描述的功能，但是这种实现不应认为超出本申请的范围。

在本申请所提供的实施例中，应该理解到，所揭露的装置/终端设备和方法，可以通过其它的方式实现。例如，以上所描述的装置/终端设备实施例仅仅是示意性的，例如，所述模块或单元的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方式，例如多个单元或组件可以结合或者可以集成到另一个系统，或一些特征可以忽略，或不执行。另一点，所显示或讨论的相互之间的耦合或直接耦合或通讯连接可以是通过一些接口，装置或单元的间接耦合或通讯连接，可以是电性，机械或其它的形式。

所述作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部单元来实现本实施例方案的目的。

另外，在本申请各个实施例中的各功能单元可以集成在一个处理单元中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个单元中。上述集成的单元既可以采用硬件的形式实现，也可以采用软件功能单元的形式实现。

所述集成的模块/单元如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本申请实现上述实施例方法中的全部或部分流程，也可以通过计算机程序来指令相关的硬件来完成，所述的计算机程序可存储于一计算机可读存储介质中，该计算机程序在被处理器执行时，可实现上述各个方法实施例的步骤。。其中，所述计算机程序包括计算机程序代码，所述计算机程序代码可以为源代码形式、对象代码形式、可执行文件或某些中间形式等。所述计算机可读介质可以包括：能够携带所述计算机程序代码的任何实体或装置、记录介质、U盘、移动硬盘、磁碟、光盘、计算机存储器、只读存储器(ROM，Read-Only Memory)、随机存取存储器(RAM，Random Access Memory)、电载波信号、电信信号以及软件分发介质等。需要说明的是，所述计算机可读介质包含的内容可以根据司法管辖区内立法和专利实践的要求进行适当的增减，例如在某些司法管辖区，根据立法和专利实践，计算机可读介质不包括是电载波信号和电信信号。

以上所述实施例仅用以说明本申请的技术方案，而非对其限制；尽管参照前述实施例对本申请进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本申请各实施例技术方案的精神和范围，均应包含在本申请的保护范围之内。

Claims (10)

1.一种颗粒材料的颗粒状态分析方法，其特征在于，所述颗粒材料的颗粒状态分析方法包括：

获取待分析的颗粒材料中的颗粒的几何形态，通过傅里叶级数表征所述颗粒材料的数学描述，所述数学描述包括颗粒形状描述符和颗粒位置描述符；

根据所述颗粒材料的数学描述，判断颗粒之间是否接触；

如果颗粒之间有接触，则根据所述颗粒材料的数学描述，计算颗粒之间的接触特征；

根据所述颗粒之间的接触特征，计算颗粒之间的接触力，根据所计算的接触力更新颗粒的运动状态；

根据所获取的颗粒当前的状态以及颗粒在下一时刻的状态，通过时间步的迭代，获取颗粒在一时间段的力学响应与状态演化。

2.根据权利要求1所述的颗粒材料的颗粒状态分析方法，其特征在于，所述获取待分析的颗粒材料中的颗粒的几何形态，通过傅里叶级数表征所述颗粒材料的数学描述，所述数学描述包括颗粒形状描述符和颗粒位置描述符的步骤包括：

在局部坐标系中，颗粒表面上的任何点通过其相对于原点的极半径r和相对于x轴的极角θ来确定，将颗粒表面上的任何点A的局部坐标(x_A，y_A)可表示为x_a＝rcosθ，y_a＝rsinθ；

将极半径r通过傅里叶展开的方式表示为： 其中，a₀，a_n，b_n为颗粒形态描述符，n为傅里叶级数；

在全局坐标系中，颗粒表面上的任何点位置表示为X_A＝X_C+r(θ)cosθ，Y_A＝Y_C+r(θ)sinθ，其中，X_C、Y_C为颗粒形心位置的颗粒位置描述符。

3.根据权利要求1所述的颗粒材料的颗粒状态分析方法，其特征在于，所述根据所述颗粒材料的数学描述，判断颗粒之间是否接触的步骤包括：

对于任意两个颗粒p，q，其颗粒形心位置分别为C_p、C_q，点A_p是颗粒p上的任意一点，向量C_pA_p的极角为v，直线C_qA_p与颗粒q的表面相交于A_q，点A_p与点A_q的距离表示为其中L表示其下标两点之间的距离；

计算所述距离D的最小值，如果所述最小值小于或等于0，则颗粒p、q相接触。

4.根据权利要求3所述的颗粒材料的颗粒状态分析方法，其特征在于，所述计算所述距离D的最小值的步骤包括：

计算线段CqAp的长度：其中，点A_p的坐标表示为： 角标p、q表示颗粒p和q，为傅里叶基函数，θ₀为颗粒的初始倾角。

计算线段C_qA_q的长度：其中，点Ap的坐标表示为：

根据线段C_qA_p和线段C_qA_q的长度，计算距离D为： 通过牛顿迭代法计算所述距离的最小值。

5.根据权利要求1所述的颗粒材料的颗粒状态分析方法，其特征在于，所述根据所述颗粒材料的数学描述，计算颗粒之间的接触特征的步骤包括：

根据计算的距离确定两个颗粒的表面的接触点；

根据所述接触点确定所述接触点对应的接触中心；

根据所述接触中心确定接触分支向量，所述接触分支向量即为接触中心指向颗粒形心的向量；

获取所述接触点的法向和切向，根据通过所述接触中心、方向为接触法向的直线与两个颗粒表面轮廓的交点，确定接触重叠点；

根据所述接触重叠点之间的距离计算接触重叠量；

计算接触相对切向位移变化量Δδ_t＝v_tΔt，其中， 为相对切向速度，v为颗粒的速度，w为颗粒的角速度，l_n为接触分支向量长度，上标p、q指示颗粒p和q，且颗粒的状态包括位置、倾角、速度和角速度，颗粒的速度、角速度由颗粒当前的状态决定。

6.根据权利要求1所述的颗粒材料的颗粒状态分析方法，其特征在于，所述根据所述颗粒之间的接触特征，计算颗粒之间的接触力，根据所计算的接触力更新颗粒的运动状态的步骤包括：

根据公式计算接触力和接触弯矩：F_n＝k_nδ_n，其中F_n为法向接触力，F_t为切向接触力；M_c＝0接触弯矩，δ_n为法向位移量，k_n、k_t为接触刚度，μ_c为接触摩擦系数；

根据所述接触力，结合颗粒重力和外部荷载，计算颗粒所受全力，根据所受全力计算颗粒的加速度、速度、移动位置。

7.一种颗粒材料的颗粒状态分析装置，其特征在于，所述颗粒材料的颗粒状态分析装置包括：

数学描述单元，用于获取待分析的颗粒材料中的颗粒的几何形态，通过傅里叶级数表征所述颗粒材料的数学描述，所述数学描述包括颗粒形状描述符和颗粒位置描述符；

接触判断单元，用于根据所述颗粒材料的数学描述，判断颗粒之间是否接触；

接触特征计算单元，用于如果颗粒之间有接触，则根据所述颗粒材料的数学描述，计算颗粒之间的接触特征；

运动状态计算单元，用于根据所述颗粒之间的接触特征，计算颗粒之间的接触力，根据所计算的接触力更新颗粒的运动状态；

响应获取单元，用于根据所获取的颗粒当前的状态以及颗粒在下一时刻的状态，通过时间步的迭代，获取颗粒在一时间段的力学响应与状态演化。

8.根据权利要求7所述的颗粒材料的颗粒状态分析装置，其特征在于，所述数学描述单元包括：

表面点表示子单元，用于在局部坐标系中，颗粒表面上的任何点通过其相对于原点的极半径r和相对于x轴的极角0来确定，将颗粒表面上的任何点A的局部坐标(x_A，y_A)可表示为x_a＝rcosθ，y_a＝rsinθ；

颗粒形态描述符获取子单元，用于将极半径r通过傅里叶展开的方式表示为：其中，a₀，a_n，b_n为颗粒形态描述符，n为傅里叶级数；

颗粒位置描述符获取子单元，用于在全局坐标系中，颗粒表面上的任何点位置表示为X_A＝X_C+r(θ)cosθ，Y_A＝Y_C+r(θ)sinθ，其中，X_C、Y_C为颗粒形心位置的颗粒位置描述符。

9.一种颗粒材料的颗粒状态分析设备，包括存储器、处理器以及存储在所述存储器中并可在所述处理器上运行的计算机程序，其特征在于，所述处理器执行所述计算机程序时实现如权利要求1至6任一项所述颗粒材料的颗粒状态分析方法的步骤。

10.一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现如权利要求1至6任一项所述颗粒材料的颗粒状态分析方法的步骤。