大规模颗粒材料内部应力及破碎模拟分析方法和装置与流程

本发明属于颗粒破碎技术领域，更具体地，涉及一种大规模颗粒材料内部应力及破碎模拟分析方法和装置。

背景技术：

颗粒材料作为一种重要的建筑材料，在工程中得以广泛的应用，例如堆石坝的主要填筑材料为堆石体，以及工民建工程中常见的砂砾石材料。由于外力作用，颗粒材料内部经常处于高应力状态，导致颗粒发生破碎，进而使得颗粒集合体的级配发生改变，影响颗粒集合体的宏细观力学性质。因此研究颗粒破碎的力学机制是十分必要的。然而受到颗粒材料尺寸和位置的影响，物理实验很难还原颗粒材料的原位力学机制分析，因而颗粒材料的数值模拟方法广受欢迎。其中，连续离散耦合方法逐渐成为一种科学的颗粒破碎数值分析手段，连续离散耦合方法主要基于有限元法，计算成本非常高。

目前所存在的离散元法，在模拟颗粒集合体相互作用时，一般将颗粒视为刚性体，无法计算颗粒内部应力大小，故而采用最大接触力原则或者设置应力阈值进行颗粒破碎判断准则，均由经验所得，理论依据薄弱。而对于现有的连续离散耦合方法，如有限元-离散元耦合法，比例边界有限元-离散元耦合法，能够利用连续介质力学中较为成熟的断裂理论进行颗粒破碎的判断，但计算效率一般较低。

对于重力等体积力引起的域积分，需要进一步处理，否则将使得边界元丧失降维的优势，目前采用的域积分处理方法主要有：双互易法(dualreciprocitymethod，drm)，在drm中，非齐次项可以用径向基函数(radialbasisfunction，rbf)等一系列函数来逼近，并应用第二个互易性将域积分转化为边界积分。只有域或边界上的点需要提供由非齐次项表述的信息。然而，drm的精度在很大程度上取决于域点的分布和位置，以及用于近似非齐次项的函数类型。此外，复杂域中的点的排列可能不容易，特别是对于三维问题。类似于drm的方法是多重互易法(multiplereciprocitymethod，mrm)，其中互易定理通过一系列高阶基本解反复应用，将域积分转化为边界。另一种方法是特殊解方法(propensityscorematching，psm)，其中构造了一个近似的特殊解，而不是在drm中执行第二个互易性将域积分转化为边界积分。

技术实现要素：

针对现有技术的以上缺陷或改进需求，本发明提出了一种大规模颗粒材料内部应力及破碎模拟分析方法和装置，具有极高的准确性和有效性。

为实现上述目的，按照本发明的一个方面，提供了一种大规模颗粒材料内部应力及破碎模拟分析方法，包括：

s1：利用颗粒相互作用求解器进行大规模颗粒集合体相互作用的模拟计算；

s2：将颗粒相互作用求解器中某一时间步的计算结果导出，在颗粒内部应力求解器中，利用建模软件建立待求解颗粒的几何模型，输入从颗粒相互作用求解器中导出的待求解颗粒模型的材料参数、网格化分数、网格类型以及边界条件，然后利用边界元法输出待求解颗粒几何模型的内部应力；

s3：将每一个颗粒边界上受到的集中力等效为边界上的面力；

s4：建立控制方程，并由控制方程得到位移边界积分方程；

s5：将位移边界积分方程中体积力导致的域积分转化为边界积分；

s6：分别对每一个颗粒建立局部坐标系，并对相似颗粒计算转换矩阵；

s7：对单颗粒边界积分方程进行离散，划分边界单元，并在域内布置节点；

s8：将每一个颗粒边界上的面力等效为边界单元上的面力；

s9：计算单颗粒的系数矩阵，并对系数矩阵求逆；

s10：求解边界节点的位移；

s11：求解域内的节点的应力和应变等值；

s12：根据内部应力值判断颗粒是否开裂，并利用颗粒替代法进行破碎颗粒的替代，并将代替颗粒的模型信息返回到颗粒相互作用求解器中。

在一些可选的实施方案中，步骤s3包括：

对于作用在中心为o(x0,y0)颗粒上一点p(xp,yp)的集中力f(fx,fy)，可等效为均布力p，其中，角度范围[θ1,θ2]为面力等效范围，且(x0,y0)表示颗粒中心点的横纵坐标，(xp,yp)表示点p的横纵坐标，(fx,fy)表示集中力f的横纵轴的分解力，r表示颗粒半径。

在一些可选的实施方案中，由建立控制方程，其中，x为研究域内的点，u为应力张量；a为加速度；ρ为颗粒材料密度；μ为剪切模量；λ＝2vμ/(1-2v)为拉美常数；v为泊松比，和为散度运算符；为合力矢量，b(x)为等效体力，s表示颗粒面积，fⁱ表示颗粒的集中力项，n表示集中力个数。

在一些可选的实施方案中，步骤s5包括：

采用直线积分法进行域积分计算，其中，域积分表示为：d1＝∫ωu(x,y)bdω(y)，基于直线积分法，由将域积分转化为等效边界积分，其中，m表示积分线个数，wi表示权重因子，li表示积分线，u(x,y)表示格林函数基本解，b表示等效体力，dy1表示微分。

在一些可选的实施方案中，步骤s6包括：

对于全局坐标中的点x(x1,x2)以及局部坐标中的点si(s1,s2)，可通过以下表达式进行转化：为jacobi矩阵的常系数，|j|＝c²,c为比例因子，逆矩阵为：

在一些可选的实施方案中，对于形状相似的颗粒，只需要在局部坐标系内建立边界积分方程，其中，边界积分方程在局部坐标中可表达为：s和为局部坐标系中的点，和为局部坐标系中的位移和面力，γ为问题研究域ω的边界，表示和点s位置相关的系数，表示局部坐标系下的基本解函数，表示局部坐标系下的边界，表示局部坐标系下的基本解函数，表示局部坐标系下点的位移，表示局部坐标系下的体力矢量，表示局部坐标系下的求解域，表示局部坐标系下点和s的距离，δki表示克罗内克符号，表示局部坐标系下r的偏导数，表示局部坐标系下r的偏导数。

在一些可选的实施方案中，步骤s8包括：

对于单元zj[s^j,s^j+1]，s^j和s^j+1为单元的两个端点，j表示单元序号，如果满足：其中，为集中力等效的范围，和为离散边界上的两个端点，单元所受等效面力为：p0为初始等效面力，β1和β2分别为单元集的起点和终点角度，p1和p2为在第一个和最后一个单元上残余力引起的压力，可由以下两组公式求得：及β3和β4分别为第一个单元的终点和最后一个单元的起点与集中力的夹角。

在一些可选的实施方案中，由确定每一个颗粒的边界节点，其中，表示边界上点的位移，h^-1表示h的逆矩阵，h表示边界积分里形成的系数矩阵，g表示边界积分里形成的系数矩阵，表示局部坐标下节点的面力，表示局部坐标下节点的体力系数矩阵，表示局部坐标下节点的体力。

在一些可选的实施方案中，由确定颗粒域内节点的应变，由确定应力值，其中，g′表示应变离散矩阵方程中对应位移的系数矩阵，h′表示应变离散矩阵方程中对应面力的系数矩阵，表示应变离散矩阵方程中对应体力的系数矩阵，m表示需要求解的内部点的数量，和为每个点在x，y方向上的应变值，表示局部坐标系下的应力，表示局部坐标系下的刚度矩阵，表示局部坐标系下的应变，δij、δkl、δik、δjl、δil及δjk表示克罗内克符号，v为泊松比。

按照本发明的另一方面，提供了一种大规模颗粒材料内部应力及破碎模拟分析装置，包括：粒相互作用求解器和颗粒内部应力求解器；

所述颗粒相互作用求解器，用于进行大规模颗粒集合体相互作用的模拟计算；

所述颗粒内部应力求解器包括：输入模块，其用于接收操作者输入的待求解颗粒的材料参数、网格化分数、网格类型以及边界条件信息；求解模块，其用于采用上述任意一项大规模颗粒材料内部应力及破碎模拟分析方法得到待求解颗粒体内部的应力。

总体而言，通过本发明所构思的以上技术方案与现有技术相比，能够取得下列有益效果：

根据本发明所提供的一种大规模颗粒材料内部应力及破碎模拟分析方法和装置，能够容易地求解出砂砾石等颗粒材料的内部应力并进行颗粒破碎判断，进而把握颗粒材料的级配变化和宏细观力学性能上的改变，以确保工程结构的安全性和可靠性。

附图说明

图1是本发明实施例提供的一种颗粒集合体内部应力及破碎模拟方法的流程示意图；

图2是本发明实施例提供的一种模型及其边界条件示意图；

图3是本发明实施例提供的一种模拟装置计算所得的内部应力云图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

本发明将边界元法和离散元法相结合进行大规模颗粒破碎研究，能够利用边界元法进行颗粒内部应力计算分析，并结合连续介质力学断裂理论，霍克布朗准则判断颗粒是否发生破碎。此外，在利用边界元进行内部应力模拟时，对于形状相似的颗粒集合体，例如：圆形颗粒，本发明只需计算一个颗粒的系数矩阵，其他相似颗粒通过坐标转换和系数缩放获得系数矩阵，大大提高了计算效率。本发明还进行了圆形堆石料的内部应力模拟，并证明其有效性。本发明能够实现相似颗粒系数矩阵的相互转化，无需对每个相似的颗粒进行计算，故而进一步提高了计算效率。对于集中力，采用边界元直接处理存在一定困难，需要将其等效到边界单元上，本发明提出了一种“两步转换法”：第一步将集中力转换为边界上均匀分布的面力；第二步将边界上的面力转化为边界单元上的面力。采用该等效方法后，一方面消除了集中力造成的奇异性，另一方面对形状相同的单元只需要计算一次系数矩阵，不同边界条件无需重复进行边界积分来处理边界条件，可大大节省计算量。本发明采用直线积分法(linelimmethod，lim)处理边界元域积分法，在直线积分法中，只使用边界单元来离散边界，不需要体单元。因此，该方法可以看作是一种边界离散化方法。在直线积分中，域积分可以用直线上一维积分的和来计算，积分线可以很容易地由边界元和oct树结构中的边界元来创建。

如图1所示，本发明实施例所提供的一种大规模颗粒材料内部应力及破碎模拟分析方法，包括以下步骤：

步骤s1：利用颗粒相互作用求解器进行大规模颗粒集合体相互作用的模拟计算；

其中，颗粒相互作用求解器表示基于离散元法求解颗粒之间相互作用的模块。

在本发明实施例中，如图2所示，为验证准确性及有效性，选取一个含有292颗圆形堆石料受到自上而下的压力的计算实例，设其弹性模量为1000mpa，泊松比为0.25。

本发明实施例中的大规模优指在100个颗粒以上的规模。

步骤s2：将颗粒相互作用求解器中某一时间步的计算结果导出，计算结果主要包括颗粒的几何信息，材料属性以及边界条件等；在颗粒内部应力求解器中，利用建模软件建立待求解颗粒的几何模型，输入从颗粒相互作用求解器中导出的待求解颗粒模型的材料参数、网格化分数、网格类型以及边界条件等信息，然后输出待求解颗粒几何模型的内部应力；

其中，颗粒内部应力求解器表示基于边界元法求解颗粒内部应力的模块。

在本发明实施例中，利用边界元法来计算颗粒材料的应力应变分布，边界元法的基本方程为：

其中，γ为问题研究域ω的边界；x和y代表问题域ω内的源点和场点；t为边界上的面力；u(x,y)和t(x,y)为基本解，可写为：

其中，r代表源点和场点之间的距离；n为边界γ的单位外法线向量。

步骤s3：将每一个颗粒边界上受到的集中力等效为边界上的面力；

在本发明实施例中，采用等效的方法将单颗粒边界上受到的集中力等效为边界上的面力，等效方法为：

对于作用在中心为o(x0,y0)颗粒上一点p(xp,yp)的集中力f(fx,fy)，可等效为均布力p，其计算表达式是为：

其中的角度范围[θ1,θ2]为面力等效范围，具体可利用下式进行计算：

其中，(x0,y0)表示颗粒中心点的横纵坐标，(xp,yp)表示点p的横纵坐标，(fx,fy)表示集中力f的横纵轴的分解力，r表示颗粒半径。

步骤s4：建立控制方程：

其中，x为研究域内的点，u为应力张量；a为加速度；ρ为颗粒材料密度；μ为剪切模量；λ＝2vμ/(1-2v)为拉美常数；v为泊松比，和为散度运算符；为合力矢量，b(x)为等效体力。s表示颗粒面积，fⁱ表示颗粒的集中力项，n表示集中力个数。

步骤s5：由控制方程建立位移边界积分方程；

其中，位移边界积分方程用来求解颗粒的位移和颗粒边界上的应力。

步骤s6：将位移边界积分方程中体积力导致的域积分转化为边界积分；

在本发明实施例中，对于由于重力等原因产生的域内积分，采用直线积分法进行计算：

对于域积分：

d1＝∫ωu(x,y)bdω(y)(7)

基于直线积分法，利用如下公式转化为等效边界积分：

其中，m表示积分线个数，wi表示权重因子，li表示积分线，u(x,y)表示格林函数基本解，b表示等效体力，dy1表示微分。

步骤s7：分别对每一个颗粒建立局部坐标系，并对相似颗粒计算转换矩阵；

本发明实施例中的相似颗粒指形状大小无限接近的颗粒。

在本发明实施例中，对每一个颗粒建立局部坐标系，对相似的颗粒建立转换矩阵，其主要方法为：

对于全局坐标中的点x(x1,x2)以及局部坐标中的点si(s1,s2)，可通过以下表达式进行转化：

其中，为jacobi矩阵的常系数，以及

其中，|j|＝c²,c为比例因子。

逆矩阵为：

步骤s8：对单颗粒边界积分方程进行离散，划分边界单元，并在域内布置节点；

在本发明实施例中，对于形状相似的颗粒，只需要在局部坐标系内建立边界积分方程：

边界积分方程在局部坐标中可表达为：

s和为局部坐标系中的点，和为局部坐标系中的位移和面力，同公式(3)所示，以及

其中，γ为问题研究域ω的边界，表示和点s位置相关的系数，表示局部坐标系下的基本解函数，表示局部坐标系下的边界，表示局部坐标系下的基本解函数，表示局部坐标系下点的位移，表示局部坐标系下的体力矢量，表示局部坐标系下的求解域，表示局部坐标系下点和s的距离，δki表示克罗内克符号，表示局部坐标系下r的偏导数，表示局部坐标系下r的偏导数。

其中，对于形状相似的颗粒，只需要求解一次边界积分方程的系数矩阵，并采用奇异值分解法求系数矩阵的逆。

步骤s9：将每一个颗粒边界上的面力等效为边界单元上的面力；

在本发明实施例中，采用等效的方法将单颗粒边界上的面力等效为边界单元上的面力，等效方法为：

对于单元zj[s^j,s^j+1]，s^j和s^j+1为单元的两个端点，j表示单元序号，如果满足：

其中，为集中力等效的范围，和为离散边界上的两个端点，单元所受等效面力为：

其中，p0为初始等效面力，其计算表达式为：

其中，β1和β2分别为单元集的起点和终点角度，p1和p2为在第一个和最后一个单元上残余力引起的压力，可由以下两组公式求得：

其中，β3和β4分别为第一个单元的终点和最后一个单元的起点与集中力的夹角。

步骤s10：计算单颗粒的系数矩阵，并对系数矩阵求逆；

步骤s11：求解边界节点的位移；

步骤s12：求解域内的节点的应力和应变等值；

在本发明实施例中，对于所有形状相似的颗粒，在只计算了局部坐标系下一个基础颗的系数矩阵及其逆矩阵后，所有相似颗粒边界节点和域内节点的未知量均可以通过矩阵相乘来计算，可大量节省计算时间：

对于每一个颗粒的边界节点可通过以下公式计算：

其中，表示边界上点的位移，h^-1表示h的逆矩阵，h表示边界积分里形成的系数矩阵，g表示边界积分里形成的系数矩阵，表示局部坐标下节点的面力，表示局部坐标下节点的体力系数矩阵，表示局部坐标下节点的体力。

由于矩阵h^-1，g和只需要计算一次，公式(19)的效率可以非常高而且适用于大规模颗粒计算。由于h是奇异矩阵，因此其逆矩阵可以通过奇异值分解(singularvaluedecomposition，svd)进行计算，如下所示：

h^-1＝v2∑⁺v1^t(20)

其中，v1和v2左奇异向量和右奇异向量，并且∑⁺为∑的广义逆，∑可写为：

其中，∑0为前n-3个奇异值(从大至小排列)，最小的三个设为零。因此可以防止颗粒的刚体位移并且获得更高的精度，n表示矩阵行数。

此外，对于颗粒域内节点的应变可以采用如下方式计算：

其中，g′表示应变离散矩阵方程中对应位移的系数矩阵，h′表示应变离散矩阵方程中对应面力的系数矩阵，表示应变离散矩阵方程中对应体力的系数矩阵。

其中：

其中，m表示需要求解的内部点的数量，和为每个点在x，y方向上的应变值。

在计算了应变值之后，应力值可以通过如下方式计算：

其中，表示局部坐标系下的应力，表示局部坐标系下的刚度矩阵，表示局部坐标系下的应变。

以及

其中，δij、δkl、δik、δjl、δil及δjk表示克罗内克符号，v为泊松比。

如图3所示，图3为颗粒集合体的内部应力云图，其中，图3中(a)表示应力sxx的应力云图，图3中(b)表示应力syy的应力云图。

步骤s13：根据应力值判断颗粒是否开裂，并利用颗粒替代法进行破碎颗粒的替代，并将代替颗粒的模型信息返回到颗粒相互作用求解器中。

本申请还提供一种大规模颗粒材料内部应力及破碎模拟分析装置，包括：颗粒相互作用求解器和颗粒内部应力求解器；

所述颗粒相互作用求解器，用于进行大规模颗粒集合体相互作用的模拟计算；

所述颗粒内部应力求解器包括：输入模块，其用于接收操作者输入的待求解颗粒的材料参数、网格化分数、网格类型以及边界条件信息；求解模块，其用于采用上述大规模颗粒材料内部应力及破碎模拟分析方法得到待求解颗粒体内部的应力。

其中，具体实施方式可以参考上述方法实施例的描述，本发明实施例将不再复述。

需要指出，根据实施的需要，可将本申请中描述的各个步骤/部件拆分为更多步骤/部件，也可将两个或多个步骤/部件或者步骤/部件的部分操作组合成新的步骤/部件，以实现本发明的目的。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。