[go: up one dir, main page]

CN112084647B - Large-scale granular material internal stress and crushing simulation analysis method and device - Google Patents

Large-scale granular material internal stress and crushing simulation analysis method and device Download PDF

Info

Publication number
CN112084647B
CN112084647B CN202010917113.3A CN202010917113A CN112084647B CN 112084647 B CN112084647 B CN 112084647B CN 202010917113 A CN202010917113 A CN 202010917113A CN 112084647 B CN112084647 B CN 112084647B
Authority
CN
China
Prior art keywords
particle
boundary
force
coordinate system
integral
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202010917113.3A
Other languages
Chinese (zh)
Other versions
CN112084647A (en
Inventor
王桥
刘彪
周伟
岳强
黄诚斌
马刚
姬翔
田文祥
黄康桥
刘摇
张思凡
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Wuhan University WHU
Original Assignee
Wuhan University WHU
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Wuhan University WHU filed Critical Wuhan University WHU
Priority to CN202010917113.3A priority Critical patent/CN112084647B/en
Publication of CN112084647A publication Critical patent/CN112084647A/en
Application granted granted Critical
Publication of CN112084647B publication Critical patent/CN112084647B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/25—Design optimisation, verification or simulation using particle-based methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Data Mining & Analysis (AREA)
  • Computational Mathematics (AREA)
  • Pure & Applied Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Algebra (AREA)
  • Geometry (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Computing Systems (AREA)
  • Operations Research (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种大规模颗粒材料内部应力及破碎模拟分析方法和装置，属于颗粒破碎技术领域，本发明基于连续离散耦合基本思想，将边界元法和离散元法相结合进行大规模颗粒破碎研究，能够利用边界元法进行颗粒内部应力计算分析，并结合连续介质力学断裂理论，霍克布朗准则判断颗粒是否发生破碎。此外，在利用边界元进行内部应力模拟时，对于形状相似的颗粒集合体，例如：圆形颗粒，本发明只需计算一个颗粒的系数矩阵，其他相似颗粒通过坐标转换和系数缩放获得系数矩阵，大大提高了计算效率。本发明还进行了圆形堆石料的内部应力模拟，并证明其有效性。

Figure 202010917113

The invention discloses a large-scale particle material internal stress and crushing simulation analysis method and device, belonging to the technical field of particle crushing. Based on the basic idea of continuous discrete coupling, the invention combines boundary element method and discrete element method to conduct large-scale particle crushing research , can use the boundary element method to calculate and analyze the internal stress of the particles, and combine with the fracture theory of continuum mechanics and the Hawke-Brown criterion to determine whether the particles are broken. In addition, when using boundary elements to simulate internal stress, for particle aggregates with similar shapes, such as circular particles, the present invention only needs to calculate the coefficient matrix of one particle, and other similar particles obtain the coefficient matrix through coordinate transformation and coefficient scaling, The computational efficiency is greatly improved. The present invention also conducts the internal stress simulation of circular rockfill and proves its effectiveness.

Figure 202010917113

Description

大规模颗粒材料内部应力及破碎模拟分析方法和装置Method and device for simulation analysis of internal stress and crushing of large-scale granular materials

技术领域technical field

本发明属于颗粒破碎技术领域，更具体地，涉及一种大规模颗粒材料内部应力及破碎模拟分析方法和装置。The invention belongs to the technical field of particle crushing, and more particularly, relates to a method and device for simulating and analyzing the internal stress and crushing of large-scale granular materials.

背景技术Background technique

颗粒材料作为一种重要的建筑材料，在工程中得以广泛的应用，例如堆石坝的主要填筑材料为堆石体，以及工民建工程中常见的砂砾石材料。由于外力作用，颗粒材料内部经常处于高应力状态，导致颗粒发生破碎，进而使得颗粒集合体的级配发生改变，影响颗粒集合体的宏细观力学性质。因此研究颗粒破碎的力学机制是十分必要的。然而受到颗粒材料尺寸和位置的影响，物理实验很难还原颗粒材料的原位力学机制分析，因而颗粒材料的数值模拟方法广受欢迎。其中，连续离散耦合方法逐渐成为一种科学的颗粒破碎数值分析手段，连续离散耦合方法主要基于有限元法，计算成本非常高。As an important building material, granular materials are widely used in engineering. For example, the main filling material of rockfill dams is rockfill body, and the common sand and gravel materials in industrial and civil construction projects. Due to the action of external force, the interior of granular materials is often in a state of high stress, which leads to the fragmentation of the particles, which in turn changes the gradation of the particle aggregates and affects the macro-mechanical properties of the particle aggregates. Therefore, it is necessary to study the mechanical mechanism of particle breakage. However, due to the influence of the size and position of the granular material, it is difficult to reduce the in-situ mechanical mechanism analysis of the granular material by physical experiments, so the numerical simulation method of the granular material is very popular. Among them, the continuous discrete coupling method has gradually become a scientific numerical analysis method for particle crushing. The continuous discrete coupling method is mainly based on the finite element method, and the computational cost is very high.

目前所存在的离散元法，在模拟颗粒集合体相互作用时，一般将颗粒视为刚性体，无法计算颗粒内部应力大小，故而采用最大接触力原则或者设置应力阈值进行颗粒破碎判断准则，均由经验所得，理论依据薄弱。而对于现有的连续离散耦合方法，如有限元-离散元耦合法，比例边界有限元-离散元耦合法，能够利用连续介质力学中较为成熟的断裂理论进行颗粒破碎的判断，但计算效率一般较低。The existing discrete element method, when simulating the interaction of particle aggregates, generally regards particles as rigid bodies and cannot calculate the internal stress of particles. Therefore, the principle of maximum contact force or the setting of stress thresholds are used to judge particle breakage. Based on experience, the theoretical basis is weak. For the existing continuous discrete coupling methods, such as the finite element-discrete element coupling method and the proportional boundary finite element-discrete element coupling method, the relatively mature fracture theory in continuum mechanics can be used to judge the particle breakage, but the calculation efficiency is average. lower.

对于重力等体积力引起的域积分，需要进一步处理，否则将使得边界元丧失降维的优势，目前采用的域积分处理方法主要有：双互易法(Dual Reciprocity Method，DRM)，在DRM中，非齐次项可以用径向基函数(Radial basis function，RBF)等一系列函数来逼近，并应用第二个互易性将域积分转化为边界积分。只有域或边界上的点需要提供由非齐次项表述的信息。然而，DRM的精度在很大程度上取决于域点的分布和位置，以及用于近似非齐次项的函数类型。此外，复杂域中的点的排列可能不容易，特别是对于三维问题。类似于DRM的方法是多重互易法(Multiple Reciprocity Method，MRM)，其中互易定理通过一系列高阶基本解反复应用，将域积分转化为边界。另一种方法是特殊解方法(PropensityScore Matching，PSM)，其中构造了一个近似的特殊解，而不是在DRM中执行第二个互易性将域积分转化为边界积分。For the domain integral caused by the gravitational iso-body force, further processing is required, otherwise the boundary element will lose the advantage of dimensionality reduction. The currently used domain integral processing methods mainly include: Dual Reciprocity Method (DRM). , the inhomogeneous term can be approximated by a series of functions such as radial basis function (RBF), and the second reciprocity is applied to convert the domain integral to the boundary integral. Only points on domains or boundaries need to provide information represented by inhomogeneous terms. However, the accuracy of DRM is highly dependent on the distribution and location of domain points, as well as the type of function used to approximate the inhomogeneous terms. Furthermore, the arrangement of points in complex domains may not be easy, especially for three-dimensional problems. A method similar to DRM is the Multiple Reciprocity Method (MRM), in which the reciprocity theorem is applied iteratively through a series of higher-order fundamental solutions to convert domain integrals into boundaries. Another approach is the special solution method (PropensityScore Matching, PSM), in which an approximate special solution is constructed instead of performing a second reciprocity in the DRM to convert domain integrals to boundary integrals.

发明内容SUMMARY OF THE INVENTION

针对现有技术的以上缺陷或改进需求，本发明提出了一种大规模颗粒材料内部应力及破碎模拟分析方法和装置，具有极高的准确性和有效性。In view of the above defects or improvement needs of the prior art, the present invention proposes a large-scale granular material internal stress and crushing simulation analysis method and device, which has extremely high accuracy and effectiveness.

为实现上述目的，按照本发明的一个方面，提供了一种大规模颗粒材料内部应力及破碎模拟分析方法，包括：In order to achieve the above object, according to one aspect of the present invention, a method for simulating and analyzing the internal stress and crushing of large-scale granular materials is provided, including:

S1：利用颗粒相互作用求解器进行大规模颗粒集合体相互作用的模拟计算；S1: Use the particle interaction solver to simulate the interaction of large-scale particle aggregates;

S2：将颗粒相互作用求解器中某一时间步的计算结果导出，在颗粒内部应力求解器中，利用建模软件建立待求解颗粒的几何模型，输入从颗粒相互作用求解器中导出的待求解颗粒模型的材料参数、网格化分数、网格类型以及边界条件，然后利用边界元法输出待求解颗粒几何模型的内部应力；S2: Export the calculation result of a certain time step in the particle interaction solver. In the particle internal stress solver, use the modeling software to establish the geometric model of the particle to be solved, and input the to-be-solved particle derived from the particle interaction solver. Material parameters, mesh fraction, mesh type and boundary conditions of the particle model, and then use the boundary element method to output the internal stress of the particle geometry model to be solved;

S3：将每一个颗粒边界上受到的集中力等效为边界上的面力；S3: The concentrated force on the boundary of each particle is equivalent to the surface force on the boundary;

S4：建立控制方程，并由控制方程得到位移边界积分方程；S4: establish the control equation, and obtain the displacement boundary integral equation from the control equation;

S5：将位移边界积分方程中体积力导致的域积分转化为边界积分；S5: Convert the domain integral caused by the body force in the displacement boundary integral equation into a boundary integral;

S6：分别对每一个颗粒建立局部坐标系，并对相似颗粒计算转换矩阵；S6: respectively establish a local coordinate system for each particle, and calculate a transformation matrix for similar particles;

S7：对单颗粒边界积分方程进行离散，划分边界单元，并在域内布置节点；S7: Discretize the single particle boundary integral equation, divide the boundary elements, and arrange nodes in the domain;

S8：将每一个颗粒边界上的面力等效为边界单元上的面力；S8: Equivalent the surface force on each particle boundary to the surface force on the boundary element;

S9：计算单颗粒的系数矩阵，并对系数矩阵求逆；S9: Calculate the coefficient matrix of a single particle, and invert the coefficient matrix;

S10：求解边界节点的位移；S10: Solve the displacement of the boundary nodes;

S11：求解域内的节点的应力和应变等值；S11: Solve the equivalent stress and strain of nodes in the domain;

S12：根据内部应力值判断颗粒是否开裂，并利用颗粒替代法进行破碎颗粒的替代，并将代替颗粒的模型信息返回到颗粒相互作用求解器中。S12: Determine whether the particles are cracked according to the internal stress value, use the particle replacement method to replace the broken particles, and return the model information of the replaced particles to the particle interaction solver.

在一些可选的实施方案中，步骤S3包括：In some optional embodiments, step S3 includes:

对于作用在中心为O(x0,y0)颗粒上一点P(xp,yp)的集中力F(Fx,Fy)，可等效为均布力p，其中，

Figure BDA0002665400760000031
角度范围[θ1,θ2]为面力等效范围，且
Figure BDA0002665400760000032
(x0,y0)表示颗粒中心点的横纵坐标，(xp,yp)表示点P的横纵坐标，(Fx,Fy)表示集中力F的横纵轴的分解力，R表示颗粒半径。For a concentrated force F(F x , F y ) acting on a point P(x p , y p ) on a particle whose center is O(x 0 , y 0 ), it can be equivalent to a uniform force p, where,
Figure BDA0002665400760000031
The angular range [θ 1 , θ 2 ] is the equivalent range of the surface force, and
Figure BDA0002665400760000032
(x 0 , y 0 ) represents the horizontal and vertical coordinates of the particle center point, (x p , y p ) represents the horizontal and vertical coordinates of point P, (F x , F y ) represents the decomposition force of the horizontal and vertical axes of the concentration force F, R represents the particle radius.

在一些可选的实施方案中，由

Figure BDA0002665400760000033
建立控制方程，其中，x为研究域内的点，u为应力张量；a为加速度；ρ为颗粒材料密度；μ为剪切模量；λ＝2vμ/(1-2v)为拉美常数；v为泊松比，
Figure BDA0002665400760000034
和
Figure BDA0002665400760000035
为散度运算符；
Figure BDA0002665400760000036
为合力矢量，b(x)为等效体力，S表示颗粒面积，FI表示颗粒的集中力项，n表示集中力个数。In some optional embodiments, by
Figure BDA0002665400760000033
Establish the governing equation, where x is the point in the research domain, u is the stress tensor; a is the acceleration; ρ is the density of the granular material; μ is the shear modulus; λ=2vμ/(1-2v) is the Latin American constant; v is Poisson's ratio,
Figure BDA0002665400760000034
and
Figure BDA0002665400760000035
is the divergence operator;
Figure BDA0002665400760000036
is the resultant force vector, b(x) is the equivalent body force, S represents the particle area, F I represents the particle's concentrated force term, and n represents the number of concentrated forces.

在一些可选的实施方案中，步骤S5包括：In some optional embodiments, step S5 includes:

采用直线积分法进行域积分计算，其中，域积分表示为：D1＝∫ΩU(x,y)bdΩ(y)，基于直线积分法，由

Figure BDA0002665400760000041
将域积分转化为等效边界积分，其中，M表示积分线个数，Wi表示权重因子，Li表示积分线，U(x,y)表示格林函数基本解，b表示等效体力，dy1表示微分。The domain integral calculation is carried out by the straight-line integral method, where the domain integral is expressed as: D 1 =∫ Ω U(x,y)bdΩ(y), based on the straight-line integral method, by
Figure BDA0002665400760000041
Convert the domain integral to the equivalent boundary integral, where M is the number of integral lines, Wi is the weight factor, Li is the integral line, U(x, y) is the basic solution of Green's function, b is the equivalent physical strength, dy 1 means differential.

在一些可选的实施方案中，步骤S6包括：In some optional embodiments, step S6 includes:

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

Figure BDA0002665400760000042
Figure BDA0002665400760000043
为Jacobi矩阵的常系数，
Figure BDA0002665400760000044
|J|＝c2,
Figure BDA0002665400760000045
c为比例因子，逆矩阵为：
Figure BDA0002665400760000046
For a point x(x 1 ,x 2 ) in global coordinates and a point s i (s 1 ,s 2 ) in local coordinates, the transformation can be done by the following expressions:
Figure BDA0002665400760000042
Figure BDA0002665400760000043
are the constant coefficients of the Jacobi matrix,
Figure BDA0002665400760000044
|J|=c 2 ,
Figure BDA0002665400760000045
c is the scale factor, and the inverse matrix is:
Figure BDA0002665400760000046

在一些可选的实施方案中，对于形状相似的颗粒，只需要在局部坐标系内建立边界积分方程，其中，边界积分方程在局部坐标中可表达为：

Figure BDA0002665400760000047
s和
Figure BDA0002665400760000048
为局部坐标系中的点，
Figure BDA0002665400760000049
和
Figure BDA00026654007600000410
为局部坐标系中的位移和面力，Γ为问题研究域Ω的边界，
Figure BDA00026654007600000411
表示和点s位置相关的系数，
Figure BDA00026654007600000412
表示局部坐标系下的基本解函数，
Figure BDA00026654007600000413
表示局部坐标系下的边界，
Figure BDA00026654007600000414
表示局部坐标系下的基本解函数，
Figure BDA00026654007600000415
表示局部坐标系下点
Figure BDA00026654007600000416
的位移，
Figure BDA00026654007600000417
表示局部坐标系下的体力矢量，
Figure BDA00026654007600000418
表示局部坐标系下的求解域，
Figure BDA00026654007600000419
表示局部坐标系下点
Figure BDA00026654007600000420
和s的距离，δki表示克罗内克符号，
Figure BDA00026654007600000421
表示局部坐标系下r的偏导数，
Figure BDA00026654007600000422
表示局部坐标系下r的偏导数。In some alternative embodiments, for particles with similar shapes, only the boundary integral equation needs to be established in the local coordinate system, where the boundary integral equation can be expressed in the local coordinate as:
Figure BDA0002665400760000047
s and
Figure BDA0002665400760000048
is a point in the local coordinate system,
Figure BDA0002665400760000049
and
Figure BDA00026654007600000410
is the displacement and surface force in the local coordinate system, Γ is the boundary of the problem research domain Ω,
Figure BDA00026654007600000411
represents the coefficient related to the position of point s,
Figure BDA00026654007600000412
represents the basic solution function in the local coordinate system,
Figure BDA00026654007600000413
represents the boundary in the local coordinate system,
Figure BDA00026654007600000414
represents the basic solution function in the local coordinate system,
Figure BDA00026654007600000415
Represents a point in the local coordinate system
Figure BDA00026654007600000416
displacement,
Figure BDA00026654007600000417
represents the body force vector in the local coordinate system,
Figure BDA00026654007600000418
represents the solution domain in the local coordinate system,
Figure BDA00026654007600000419
Represents a point in the local coordinate system
Figure BDA00026654007600000420
the distance from s, δ ki denotes Kronecker notation,
Figure BDA00026654007600000421
represents the partial derivative of r in the local coordinate system,
Figure BDA00026654007600000422
represents the partial derivative of r in the local coordinate system.

在一些可选的实施方案中，步骤S8包括：In some optional embodiments, step S8 includes:

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

Figure BDA0002665400760000051
其中，
Figure BDA0002665400760000052
为集中力等效的范围，
Figure BDA0002665400760000053
和
Figure BDA0002665400760000054
为离散边界上的两个端点，单元所受等效面力为：
Figure BDA0002665400760000055
p0为初始等效面力，
Figure BDA0002665400760000056
β1和β2分别为单元集的起点和终点角度，p1和p2为在第一个和最后一个单元上残余力引起的压力，可由以下两组公式求得：
Figure BDA0002665400760000057
及
Figure BDA0002665400760000058
β3和β4分别为第一个单元的终点和最后一个单元的起点与集中力的夹角。For the unit Z j [s j ,s j+1 ], s j and s j+1 are the two endpoints of the unit, and j represents the unit number, if it satisfies:
Figure BDA0002665400760000051
in,
Figure BDA0002665400760000052
is the equivalent range of concentrated force,
Figure BDA0002665400760000053
and
Figure BDA0002665400760000054
are the two endpoints on the discrete boundary, and the equivalent surface force on the element is:
Figure BDA0002665400760000055
p 0 is the initial equivalent surface force,
Figure BDA0002665400760000056
β 1 and β 2 are the starting and ending angles of the element set, respectively, and p 1 and p 2 are the pressures due to residual forces on the first and last elements, which can be obtained from the following two sets of formulas:
Figure BDA0002665400760000057
and
Figure BDA0002665400760000058
β 3 and β 4 are the angles between the end point of the first element and the start point of the last element and the concentrated force, respectively.

在一些可选的实施方案中，由

Figure BDA0002665400760000059
确定每一个颗粒的边界节点，其中，
Figure BDA00026654007600000510
表示边界上点的位移，H-1表示H的逆矩阵，H表示边界积分里形成的系数矩阵，G表示边界积分里形成的系数矩阵，
Figure BDA00026654007600000511
表示局部坐标下节点的面力，
Figure BDA00026654007600000512
表示局部坐标下节点的体力系数矩阵，
Figure BDA00026654007600000513
表示局部坐标下节点的体力。In some optional embodiments, by
Figure BDA0002665400760000059
Determine the boundary nodes of each particle, where,
Figure BDA00026654007600000510
represents the displacement of the point on the boundary, H -1 represents the inverse matrix of H, H represents the coefficient matrix formed in the boundary integral, G represents the coefficient matrix formed in the boundary integral,
Figure BDA00026654007600000511
represents the surface force of the node in local coordinates,
Figure BDA00026654007600000512
represents the physical strength coefficient matrix of nodes in local coordinates,
Figure BDA00026654007600000513
Represents the physical force of the node in local coordinates.

在一些可选的实施方案中，由

Figure BDA00026654007600000514
确定颗粒域内节点的应变，由
Figure BDA00026654007600000515
确定应力值，其中，G′表示应变离散矩阵方程中对应位移的系数矩阵，H′表示应变离散矩阵方程中对应面力的系数矩阵，
Figure BDA00026654007600000516
表示应变离散矩阵方程中对应体力的系数矩阵，
Figure BDA00026654007600000517
m表示需要求解的内部点的数量，
Figure BDA00026654007600000518
和
Figure BDA00026654007600000519
为每个点在x，y方向上的应变值，
Figure BDA00026654007600000520
表示局部坐标系下的应力，
Figure BDA00026654007600000521
表示局部坐标系下的刚度矩阵，
Figure BDA00026654007600000522
表示局部坐标系下的应变，
Figure BDA00026654007600000523
δij、δkl、δik、δjl、δil及δjk表示克罗内克符号，v为泊松比。In some optional embodiments, by
Figure BDA00026654007600000514
Determine the strain at the nodes within the particle domain, given by
Figure BDA00026654007600000515
Determine the stress value, where G′ represents the coefficient matrix of the corresponding displacement in the strain discrete matrix equation, H′ represents the coefficient matrix of the corresponding surface force in the strain discrete matrix equation,
Figure BDA00026654007600000516
represents the coefficient matrix of the corresponding body force in the discrete matrix equation of strain,
Figure BDA00026654007600000517
m represents the number of interior points to be solved,
Figure BDA00026654007600000518
and
Figure BDA00026654007600000519
is the strain value of each point in the x, y direction,
Figure BDA00026654007600000520
represents the stress in the local coordinate system,
Figure BDA00026654007600000521
represents the stiffness matrix in the local coordinate system,
Figure BDA00026654007600000522
represents the strain in the local coordinate system,
Figure BDA00026654007600000523
δ ij , δ kl , δ ik , δ jl , δ il and δ jk represent Kronecker symbols, and v is Poisson's ratio.

按照本发明的另一方面，提供了一种大规模颗粒材料内部应力及破碎模拟分析装置，包括：粒相互作用求解器和颗粒内部应力求解器；According to another aspect of the present invention, a large-scale granular material internal stress and crushing simulation and analysis device is provided, including: a particle interaction solver and a particle internal stress solver;

所述颗粒相互作用求解器，用于进行大规模颗粒集合体相互作用的模拟计算；The particle interaction solver is used for simulation calculation of large-scale particle aggregate interaction;

所述颗粒内部应力求解器包括：输入模块，其用于接收操作者输入的待求解颗粒的材料参数、网格化分数、网格类型以及边界条件信息；求解模块，其用于采用上述任意一项大规模颗粒材料内部应力及破碎模拟分析方法得到待求解颗粒体内部的应力。The particle internal stress solver includes: an input module, which is used to receive the material parameters, meshing fraction, mesh type and boundary condition information of the particle to be solved input by the operator; a solution module, which is used to adopt any one of the above. A large-scale granular material internal stress and crushing simulation analysis method is used to obtain the internal stress of the particle to be solved.

总体而言，通过本发明所构思的以上技术方案与现有技术相比，能够取得下列有益效果：In general, compared with the prior art, the above technical solutions conceived by the present invention can achieve the following beneficial effects:

根据本发明所提供的一种大规模颗粒材料内部应力及破碎模拟分析方法和装置，能够容易地求解出砂砾石等颗粒材料的内部应力并进行颗粒破碎判断，进而把握颗粒材料的级配变化和宏细观力学性能上的改变，以确保工程结构的安全性和可靠性。According to the method and device for simulating and analyzing the internal stress and crushing of large-scale granular materials provided by the present invention, the internal stress of granular materials such as sand and gravel can be easily solved and the particle crushing judgment can be carried out, and then the gradation changes and macroscopic changes of granular materials can be grasped. Changes in meso-mechanical properties to ensure the safety and reliability of engineering structures.

附图说明Description of drawings

图1是本发明实施例提供的一种颗粒集合体内部应力及破碎模拟方法的流程示意图；1 is a schematic flowchart of a method for simulating the internal stress and crushing of a particle aggregate provided by an embodiment of the present invention;

图2是本发明实施例提供的一种模型及其边界条件示意图；2 is a schematic diagram of a model and its boundary conditions provided by an embodiment of the present invention;

图3是本发明实施例提供的一种模拟装置计算所得的内部应力云图。FIG. 3 is an internal stress nephogram calculated by a simulation device provided in an embodiment of the present invention.

具体实施方式Detailed ways

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。In order to make the objectives, technical solutions and advantages of the present invention clearer, the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention, but not to limit the present invention. In addition, the technical features involved in the various embodiments of the present invention described below can be combined with each other as long as they do not conflict with each other.

本发明将边界元法和离散元法相结合进行大规模颗粒破碎研究，能够利用边界元法进行颗粒内部应力计算分析，并结合连续介质力学断裂理论，霍克布朗准则判断颗粒是否发生破碎。此外，在利用边界元进行内部应力模拟时，对于形状相似的颗粒集合体，例如：圆形颗粒，本发明只需计算一个颗粒的系数矩阵，其他相似颗粒通过坐标转换和系数缩放获得系数矩阵，大大提高了计算效率。本发明还进行了圆形堆石料的内部应力模拟，并证明其有效性。本发明能够实现相似颗粒系数矩阵的相互转化，无需对每个相似的颗粒进行计算，故而进一步提高了计算效率。对于集中力，采用边界元直接处理存在一定困难，需要将其等效到边界单元上，本发明提出了一种“两步转换法”：第一步将集中力转换为边界上均匀分布的面力；第二步将边界上的面力转化为边界单元上的面力。采用该等效方法后，一方面消除了集中力造成的奇异性，另一方面对形状相同的单元只需要计算一次系数矩阵，不同边界条件无需重复进行边界积分来处理边界条件，可大大节省计算量。本发明采用直线积分法(Line Lim Method，LIM)处理边界元域积分法，在直线积分法中，只使用边界单元来离散边界，不需要体单元。因此，该方法可以看作是一种边界离散化方法。在直线积分中，域积分可以用直线上一维积分的和来计算，积分线可以很容易地由边界元和OCT树结构中的边界元来创建。The invention combines the boundary element method and the discrete element method for large-scale particle crushing research, and can use the boundary element method to calculate and analyze the internal stress of the particles, and combine the continuum mechanics fracture theory and the Hawke-Brown criterion to determine whether the particles are broken. In addition, when using boundary elements to simulate internal stress, for particle aggregates with similar shapes, such as circular particles, the present invention only needs to calculate the coefficient matrix of one particle, and other similar particles obtain the coefficient matrix through coordinate transformation and coefficient scaling, The computational efficiency is greatly improved. The present invention also conducts the internal stress simulation of circular rockfill and proves its effectiveness. The invention can realize the mutual conversion of the coefficient matrix of similar particles, and does not need to calculate each similar particle, so the calculation efficiency is further improved. For the concentrated force, it is difficult to use the boundary element to directly deal with it, and it needs to be equivalent to the boundary element. force; the second step converts the surface force on the boundary to the surface force on the boundary element. After adopting this equivalent method, on the one hand, the singularity caused by the concentrated force is eliminated, on the other hand, the coefficient matrix only needs to be calculated once for the elements with the same shape, and there is no need to repeat the boundary integration to deal with the boundary conditions for different boundary conditions, which can greatly save the calculation. quantity. The present invention adopts the line integral method (Line Lim Method, LIM) to process the boundary element domain integral method. In the linear integral method, only the boundary element is used to discretize the boundary, and the volume element is not needed. Therefore, this method can be regarded as a boundary discretization method. In line integrals, domain integrals can be computed with the sum of dimensional integrals on a straight line, and integral lines can be easily created from boundary elements and boundary elements in the OCT tree structure.

如图1所示，本发明实施例所提供的一种大规模颗粒材料内部应力及破碎模拟分析方法，包括以下步骤：As shown in FIG. 1 , a method for simulating and analyzing the internal stress and crushing of a large-scale granular material provided by an embodiment of the present invention includes the following steps:

步骤S1：利用颗粒相互作用求解器进行大规模颗粒集合体相互作用的模拟计算；Step S1: use a particle interaction solver to simulate the interaction of large-scale particle aggregates;

其中，颗粒相互作用求解器表示基于离散元法求解颗粒之间相互作用的模块。Among them, the particle interaction solver represents a module that solves the interaction between particles based on the discrete element method.

在本发明实施例中，如图2所示，为验证准确性及有效性，选取一个含有292颗圆形堆石料受到自上而下的压力的计算实例，设其弹性模量为1000Mpa，泊松比为0.25。In the embodiment of the present invention, as shown in Fig. 2, in order to verify the accuracy and validity, a calculation example containing 292 circular rockfills subjected to pressure from top to bottom is selected, and the elastic modulus is set to be 1000Mpa, and the poise The loose ratio is 0.25.

本发明实施例中的大规模优指在100个颗粒以上的规模。The large scale in the embodiment of the present invention preferably refers to the scale of more than 100 particles.

步骤S2：将颗粒相互作用求解器中某一时间步的计算结果导出，计算结果主要包括颗粒的几何信息，材料属性以及边界条件等；在颗粒内部应力求解器中，利用建模软件建立待求解颗粒的几何模型，输入从颗粒相互作用求解器中导出的待求解颗粒模型的材料参数、网格化分数、网格类型以及边界条件等信息，然后输出待求解颗粒几何模型的内部应力；Step S2: Export the calculation result of a certain time step in the particle interaction solver, and the calculation result mainly includes the geometric information of the particle, material properties and boundary conditions, etc.; in the particle internal stress solver, the modeling software is used to establish the solution to be solved. For the particle geometry model, input the material parameters, mesh fraction, mesh type, and boundary conditions of the particle model to be solved derived from the particle interaction solver, and then output the internal stress of the particle geometry model to be solved;

其中，颗粒内部应力求解器表示基于边界元法求解颗粒内部应力的模块。Among them, the particle internal stress solver represents a module that solves the internal stress of particles based on the boundary element method.

在本发明实施例中，利用边界元法来计算颗粒材料的应力应变分布，边界元法的基本方程为：In the embodiment of the present invention, the boundary element method is used to calculate the stress-strain distribution of the granular material, and the basic equation of the boundary element method is:

Figure BDA0002665400760000081
Figure BDA0002665400760000081

其中，Γ为问题研究域Ω的边界；x和y代表问题域Ω内的源点和场点；t为边界上的面力；U(x,y)和T(x,y)为基本解，可写为：Among them, Γ is the boundary of the problem domain Ω; x and y represent the source and field points in the problem domain Ω; t is the surface force on the boundary; U(x,y) and T(x,y) are the basic solutions , which can be written as:

Figure BDA0002665400760000082
Figure BDA0002665400760000082

Figure BDA0002665400760000083
Figure BDA0002665400760000083

其中，r代表源点和场点之间的距离；n为边界Γ的单位外法线向量。where r represents the distance between the source point and the field point; n is the unit outer normal vector of the boundary Γ.

步骤S3：将每一个颗粒边界上受到的集中力等效为边界上的面力；Step S3: Equating the concentrated force on the boundary of each particle as the surface force on the boundary;

在本发明实施例中，采用等效的方法将单颗粒边界上受到的集中力等效为边界上的面力，等效方法为：In the embodiment of the present invention, an equivalent method is used to equate the concentrated force on the boundary of a single particle as the surface force on the boundary, and the equivalent method is:

对于作用在中心为O(x0,y0)颗粒上一点P(xp,yp)的集中力F(Fx,Fy)，可等效为均布力p，其计算表达式是为：For the concentrated force F(F x ,F y ) acting on a point P(x p ,y p ) on a particle whose center is O(x 0 ,y 0 ), it can be equivalent to a uniform force p, and its calculation expression is for:

Figure BDA0002665400760000091
Figure BDA0002665400760000091

其中的角度范围[θ1,θ2]为面力等效范围，具体可利用下式进行计算：The angle range [θ 1 , θ 2 ] is the equivalent range of the surface force, which can be calculated by the following formula:

Figure BDA0002665400760000092
Figure BDA0002665400760000092

其中，(x0,y0)表示颗粒中心点的横纵坐标，(xp,yp)表示点P的横纵坐标，(Fx,Fy)表示集中力F的横纵轴的分解力，R表示颗粒半径。Among them, (x 0 , y 0 ) represents the horizontal and vertical coordinates of the particle center point, (x p , y p ) represents the horizontal and vertical coordinates of point P, and (F x , F y ) represents the decomposition of the horizontal and vertical axes of the concentration force F force, and R is the particle radius.

步骤S4：建立控制方程：Step S4: Establish the control equation:

Figure BDA0002665400760000093
Figure BDA0002665400760000093

其中，x为研究域内的点，u为应力张量；a为加速度；ρ为颗粒材料密度；μ为剪切模量；λ＝2vμ/(1-2v)为拉美常数；v为泊松比，

Figure BDA0002665400760000094
和
Figure BDA0002665400760000095
为散度运算符；
Figure BDA0002665400760000096
为合力矢量，b(x)为等效体力。S表示颗粒面积，FI表示颗粒的集中力项，n表示集中力个数。where x is the point in the research domain, u is the stress tensor; a is the acceleration; ρ is the density of the granular material; μ is the shear modulus; λ=2vμ/(1-2v) is the Latin American constant; v is the Poisson’s ratio ,
Figure BDA0002665400760000094
and
Figure BDA0002665400760000095
is the divergence operator;
Figure BDA0002665400760000096
is the resultant force vector, and b(x) is the equivalent physical force. S is the particle area, F I is the concentration term of the particle, and n is the number of the concentration.

步骤S5：由控制方程建立位移边界积分方程；Step S5: establish the displacement boundary integral equation from the control equation;

其中，位移边界积分方程用来求解颗粒的位移和颗粒边界上的应力。Among them, the displacement boundary integral equation is used to solve the displacement of the particle and the stress on the particle boundary.

步骤S6：将位移边界积分方程中体积力导致的域积分转化为边界积分；Step S6: Convert the domain integral caused by the body force in the displacement boundary integral equation into a boundary integral;

在本发明实施例中，对于由于重力等原因产生的域内积分，采用直线积分法进行计算：In the embodiment of the present invention, for the integral in the domain due to gravity and other reasons, the linear integral method is used to calculate:

对于域积分：For domain integrals:

D1＝∫ΩU(x,y)bdΩ(y) (7)D 1 = ∫Ω U(x,y)bdΩ(y) (7)

基于直线积分法，利用如下公式转化为等效边界积分：Based on the straight-line integral method, the following formula is used to convert to the equivalent boundary integral:

Figure BDA0002665400760000097
Figure BDA0002665400760000097

其中，M表示积分线个数，Wi表示权重因子，Li表示积分线，U(x,y)表示格林函数基本解，b表示等效体力，dy1表示微分。Among them, M represents the number of integral lines, Wi represents the weight factor, Li represents the integral line, U(x, y) represents the basic solution of Green's function, b represents the equivalent physical strength, and dy 1 represents the differential.

步骤S7：分别对每一个颗粒建立局部坐标系，并对相似颗粒计算转换矩阵；Step S7: establishing a local coordinate system for each particle respectively, and calculating a transformation matrix for similar particles;

本发明实施例中的相似颗粒指形状大小无限接近的颗粒。Similar particles in the embodiments of the present invention refer to particles with infinitely close shapes and sizes.

在本发明实施例中，对每一个颗粒建立局部坐标系，对相似的颗粒建立转换矩阵，其主要方法为：In the embodiment of the present invention, a local coordinate system is established for each particle, and a transformation matrix is established for similar particles, and the main methods are:

对于全局坐标中的点x(x1,x2)以及局部坐标中的点si(s1,s2)，可通过以下表达式进行转化：For a point x(x 1 ,x 2 ) in global coordinates and a point s i (s 1 ,s 2 ) in local coordinates, the transformation can be done by the following expressions:

Figure BDA0002665400760000101
Figure BDA0002665400760000101

其中，

Figure BDA0002665400760000102
为Jacobi矩阵的常系数，以及in,
Figure BDA0002665400760000102
are the constant coefficients of the Jacobi matrix, and

Figure BDA0002665400760000103
Figure BDA0002665400760000103

其中，

Figure BDA0002665400760000104
|J|＝c2,
Figure BDA0002665400760000105
c为比例因子。in,
Figure BDA0002665400760000104
|J|=c 2 ,
Figure BDA0002665400760000105
c is the scale factor.

逆矩阵为：The inverse matrix is:

Figure BDA0002665400760000106
Figure BDA0002665400760000106

步骤S8：对单颗粒边界积分方程进行离散，划分边界单元，并在域内布置节点；Step S8: discretize the single particle boundary integral equation, divide the boundary elements, and arrange nodes in the domain;

在本发明实施例中，对于形状相似的颗粒，只需要在局部坐标系内建立边界积分方程：In the embodiment of the present invention, for particles with similar shapes, only the boundary integral equation needs to be established in the local coordinate system:

边界积分方程在局部坐标中可表达为：The boundary integral equation can be expressed in local coordinates as:

Figure BDA0002665400760000111
Figure BDA0002665400760000111

s和

Figure BDA0002665400760000112
为局部坐标系中的点，
Figure BDA0002665400760000113
和
Figure BDA0002665400760000114
为局部坐标系中的位移和面力，
Figure BDA0002665400760000115
同公式(3)所示，以及s and
Figure BDA0002665400760000112
is a point in the local coordinate system,
Figure BDA0002665400760000113
and
Figure BDA0002665400760000114
are the displacement and surface force in the local coordinate system,
Figure BDA0002665400760000115
as shown in formula (3), and

Figure BDA0002665400760000116
Figure BDA0002665400760000116

其中，Γ为问题研究域Ω的边界，

Figure BDA0002665400760000117
表示和点s位置相关的系数，
Figure BDA0002665400760000118
表示局部坐标系下的基本解函数，
Figure BDA0002665400760000119
表示局部坐标系下的边界，
Figure BDA00026654007600001110
表示局部坐标系下的基本解函数，
Figure BDA00026654007600001111
表示局部坐标系下点
Figure BDA00026654007600001112
的位移，
Figure BDA00026654007600001113
表示局部坐标系下的体力矢量，
Figure BDA00026654007600001114
表示局部坐标系下的求解域，
Figure BDA00026654007600001115
表示局部坐标系下点
Figure BDA00026654007600001116
和s的距离，δki表示克罗内克符号，
Figure BDA00026654007600001117
表示局部坐标系下r的偏导数，
Figure BDA00026654007600001118
表示局部坐标系下r的偏导数。Among them, Γ is the boundary of the problem research domain Ω,
Figure BDA0002665400760000117
represents the coefficient related to the position of point s,
Figure BDA0002665400760000118
represents the basic solution function in the local coordinate system,
Figure BDA0002665400760000119
represents the boundary in the local coordinate system,
Figure BDA00026654007600001110
represents the basic solution function in the local coordinate system,
Figure BDA00026654007600001111
Represents a point in the local coordinate system
Figure BDA00026654007600001112
displacement,
Figure BDA00026654007600001113
represents the body force vector in the local coordinate system,
Figure BDA00026654007600001114
represents the solution domain in the local coordinate system,
Figure BDA00026654007600001115
Represents a point in the local coordinate system
Figure BDA00026654007600001116
the distance from s, δ ki denotes Kronecker notation,
Figure BDA00026654007600001117
represents the partial derivative of r in the local coordinate system,
Figure BDA00026654007600001118
represents the partial derivative of r in the local coordinate system.

其中，对于形状相似的颗粒，只需要求解一次边界积分方程的系数矩阵，并采用奇异值分解法求系数矩阵的逆。Among them, for particles with similar shapes, only the coefficient matrix of the first-order boundary integral equation needs to be solved, and the singular value decomposition method is used to obtain the inverse of the coefficient matrix.

步骤S9：将每一个颗粒边界上的面力等效为边界单元上的面力；Step S9: Equivalent the surface force on each particle boundary to the surface force on the boundary element;

在本发明实施例中，采用等效的方法将单颗粒边界上的面力等效为边界单元上的面力，等效方法为：In the embodiment of the present invention, the equivalent method is used to equalize the surface force on the single particle boundary to the surface force on the boundary element, and the equivalent method is:

对于单元Zj[sj,sj+1]，sj和sj+1为单元的两个端点，j表示单元序号，如果满足：For the unit Z j [s j ,s j+1 ], s j and s j+1 are the two endpoints of the unit, and j represents the unit number, if it satisfies:

Figure BDA00026654007600001119
Figure BDA00026654007600001119

其中，

Figure BDA00026654007600001120
为集中力等效的范围，
Figure BDA00026654007600001121
和
Figure BDA00026654007600001122
为离散边界上的两个端点，单元所受等效面力为：in,
Figure BDA00026654007600001120
is the equivalent range of concentrated force,
Figure BDA00026654007600001121
and
Figure BDA00026654007600001122
are the two endpoints on the discrete boundary, and the equivalent surface force on the element is:

Figure BDA0002665400760000121
Figure BDA0002665400760000121

其中，p0为初始等效面力，其计算表达式为：Among them, p 0 is the initial equivalent surface force, and its calculation expression is:

Figure BDA0002665400760000122
Figure BDA0002665400760000122

其中，β1和β2分别为单元集的起点和终点角度，p1和p2为在第一个和最后一个单元上残余力引起的压力，可由以下两组公式求得：where β 1 and β 2 are the starting and ending angles of the element set, respectively, and p 1 and p 2 are the pressures caused by residual forces on the first and last elements, which can be obtained from the following two sets of formulas:

Figure BDA0002665400760000123
Figure BDA0002665400760000123

Figure BDA0002665400760000124
Figure BDA0002665400760000124

其中，β3和β4分别为第一个单元的终点和最后一个单元的起点与集中力的夹角。Among them, β 3 and β 4 are the angles between the end point of the first element and the start point of the last element and the concentrated force, respectively.

步骤S10：计算单颗粒的系数矩阵，并对系数矩阵求逆；Step S10: calculating the coefficient matrix of the single particle, and inverting the coefficient matrix;

步骤S11：求解边界节点的位移；Step S11: solve the displacement of the boundary nodes;

步骤S12：求解域内的节点的应力和应变等值；Step S12: Solve the stress and strain equivalents of nodes in the domain;

在本发明实施例中，对于所有形状相似的颗粒，在只计算了局部坐标系下一个基础颗的系数矩阵及其逆矩阵后，所有相似颗粒边界节点和域内节点的未知量均可以通过矩阵相乘来计算，可大量节省计算时间：In the embodiment of the present invention, for all particles with similar shapes, after only calculating the coefficient matrix of the next basic particle in the local coordinate system and its inverse matrix, the unknowns of all similar particle boundary nodes and nodes in the domain can be calculated through the matrix phase. Multiply to calculate, which can save a lot of calculation time:

对于每一个颗粒的边界节点可通过以下公式计算：The boundary node for each particle can be calculated by the following formula:

Figure BDA0002665400760000125
Figure BDA0002665400760000125

其中，

Figure BDA0002665400760000126
表示边界上点的位移，H-1表示H的逆矩阵，H表示边界积分里形成的系数矩阵，G表示边界积分里形成的系数矩阵，
Figure BDA0002665400760000127
表示局部坐标下节点的面力，
Figure BDA0002665400760000128
表示局部坐标下节点的体力系数矩阵，
Figure BDA0002665400760000129
表示局部坐标下节点的体力。in,
Figure BDA0002665400760000126
represents the displacement of the point on the boundary, H -1 represents the inverse matrix of H, H represents the coefficient matrix formed in the boundary integral, G represents the coefficient matrix formed in the boundary integral,
Figure BDA0002665400760000127
represents the surface force of the node in local coordinates,
Figure BDA0002665400760000128
represents the physical strength coefficient matrix of nodes in local coordinates,
Figure BDA0002665400760000129
Represents the physical force of the node in local coordinates.

由于矩阵H-1，G和

Figure BDA00026654007600001210
只需要计算一次，公式(19)的效率可以非常高而且适用于大规模颗粒计算。由于H是奇异矩阵，因此其逆矩阵可以通过奇异值分解(Singular ValueDecomposition，SVD)进行计算，如下所示：Since the matrix H -1 , G and
Figure BDA00026654007600001210
With only one calculation required, equation (19) can be very efficient and suitable for large-scale particle calculations. Since H is a singular matrix, its inverse matrix can be calculated by Singular Value Decomposition (SVD) as follows:

H-1＝V2∑+V1 T (20)H -1 = V 2 ∑ + V 1 T (20)

其中，V1和V2左奇异向量和右奇异向量，并且∑+为∑的广义逆，∑可写为：where V 1 and V 2 left singular vectors and right singular vectors, and ∑ + is the generalized inverse of ∑, ∑ can be written as:

Figure BDA0002665400760000131
Figure BDA0002665400760000131

其中，∑0为前N-3个奇异值(从大至小排列)，最小的三个设为零。因此可以防止颗粒的刚体位移并且获得更高的精度，N表示矩阵行数。Among them, Σ 0 is the first N-3 singular values (arranged from largest to smallest), and the smallest three are set to zero. Therefore, rigid body displacement of the particles can be prevented and higher accuracy can be obtained, N represents the number of matrix rows.

此外，对于颗粒域内节点的应变可以采用如下方式计算：In addition, the strain of nodes in the particle domain can be calculated as follows:

Figure BDA0002665400760000132
Figure BDA0002665400760000132

其中，G′表示应变离散矩阵方程中对应位移的系数矩阵，H′表示应变离散矩阵方程中对应面力的系数矩阵，

Figure BDA0002665400760000133
表示应变离散矩阵方程中对应体力的系数矩阵。Among them, G′ represents the coefficient matrix of the corresponding displacement in the strain discrete matrix equation, H′ represents the coefficient matrix of the corresponding surface force in the strain discrete matrix equation,
Figure BDA0002665400760000133
Represents the matrix of coefficients for the corresponding body forces in the discrete matrix equation of strain.

其中：in:

Figure BDA0002665400760000134
Figure BDA0002665400760000134

其中，m表示需要求解的内部点的数量，

Figure BDA0002665400760000135
和
Figure BDA0002665400760000136
为每个点在x，y方向上的应变值。where m represents the number of interior points to be solved,
Figure BDA0002665400760000135
and
Figure BDA0002665400760000136
is the strain value of each point in the x, y direction.

在计算了应变值之后，应力值可以通过如下方式计算：After calculating the strain value, the stress value can be calculated as follows:

Figure BDA0002665400760000137
Figure BDA0002665400760000137

其中，

Figure BDA0002665400760000138
表示局部坐标系下的应力，
Figure BDA0002665400760000139
表示局部坐标系下的刚度矩阵，
Figure BDA00026654007600001310
表示局部坐标系下的应变。in,
Figure BDA0002665400760000138
represents the stress in the local coordinate system,
Figure BDA0002665400760000139
represents the stiffness matrix in the local coordinate system,
Figure BDA00026654007600001310
represents the strain in the local coordinate system.

以及as well as

Figure BDA0002665400760000141
Figure BDA0002665400760000141

其中，δij、δkl、δik、δjl、δil及δjk表示克罗内克符号，v为泊松比。Among them, δ ij , δ kl , δ ik , δ jl , δ il and δ jk represent Kronecker symbols, and v is Poisson's ratio.

如图3所示，图3为颗粒集合体的内部应力云图，其中，图3中(a)表示应力sxx的应力云图，图3中(b)表示应力syy的应力云图。As shown in FIG. 3 , FIG. 3 is the internal stress nephogram of the particle aggregate, wherein (a) in FIG. 3 represents the stress nephogram of stress sxx, and (b) in FIG. 3 represents the stress nephogram of stress syy.

步骤S13：根据应力值判断颗粒是否开裂，并利用颗粒替代法进行破碎颗粒的替代，并将代替颗粒的模型信息返回到颗粒相互作用求解器中。Step S13: Determine whether the particles are cracked according to the stress value, use the particle replacement method to replace the broken particles, and return the model information of the replaced particles to the particle interaction solver.

本申请还提供一种大规模颗粒材料内部应力及破碎模拟分析装置，包括：颗粒相互作用求解器和颗粒内部应力求解器；The application also provides a large-scale particle material internal stress and crushing simulation analysis device, including: a particle interaction solver and a particle internal stress solver;

所述颗粒相互作用求解器，用于进行大规模颗粒集合体相互作用的模拟计算；The particle interaction solver is used for simulation calculation of large-scale particle aggregate interaction;

所述颗粒内部应力求解器包括：输入模块，其用于接收操作者输入的待求解颗粒的材料参数、网格化分数、网格类型以及边界条件信息；求解模块，其用于采用上述大规模颗粒材料内部应力及破碎模拟分析方法得到待求解颗粒体内部的应力。The particle internal stress solver includes: an input module, which is used to receive the material parameters, meshing fraction, mesh type and boundary condition information of the particle to be solved input by the operator; The internal stress and crushing simulation analysis method of granular material is used to obtain the internal stress of the particle to be solved.

其中，具体实施方式可以参考上述方法实施例的描述，本发明实施例将不再复述。For the specific implementation manner, reference may be made to the descriptions of the foregoing method embodiments, which will not be repeated in the embodiments of the present invention.

需要指出，根据实施的需要，可将本申请中描述的各个步骤/部件拆分为更多步骤/部件，也可将两个或多个步骤/部件或者步骤/部件的部分操作组合成新的步骤/部件，以实现本发明的目的。It should be pointed out that, according to the needs of implementation, the various steps/components described in this application may be split into more steps/components, or two or more steps/components or partial operations of steps/components may be combined into new steps/components to achieve the purpose of the present invention.

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。Those skilled in the art can easily understand that the above are only preferred embodiments of the present invention, and are not intended to limit the present invention. Any modifications, equivalent replacements and improvements made within the spirit and principles of the present invention, etc., All should be included within the protection scope of the present invention.

Claims (4)

1.一种大规模颗粒材料内部应力及破碎模拟分析方法，其特征在于，包括：1. a large-scale granular material internal stress and crushing simulation analysis method, is characterized in that, comprises: S1：利用颗粒相互作用求解器进行大规模颗粒集合体相互作用的模拟计算；S1: Use the particle interaction solver to simulate the interaction of large-scale particle aggregates; S2：将颗粒相互作用求解器中某一时间步的计算结果导出，在颗粒内部应力求解器中，利用建模软件建立待求解颗粒的几何模型，输入从颗粒相互作用求解器中导出的待求解颗粒模型的材料参数、网格化分数、网格类型以及边界条件，然后利用边界元法输出待求解颗粒几何模型的内部应力；S2: Export the calculation result of a certain time step in the particle interaction solver. In the particle internal stress solver, use the modeling software to establish the geometric model of the particle to be solved, and input the to-be-solved particle derived from the particle interaction solver. Material parameters, mesh fraction, mesh type and boundary conditions of the particle model, and then use the boundary element method to output the internal stress of the particle geometry model to be solved; S3：将每一个颗粒边界上受到的集中力等效为边界上的面力；S3: The concentrated force on the boundary of each particle is equivalent to the surface force on the boundary; S4：建立控制方程，并由控制方程得到位移边界积分方程；S4: establish the control equation, and obtain the displacement boundary integral equation from the control equation; S5：将位移边界积分方程中体积力导致的域积分转化为边界积分；S5: Convert the domain integral caused by the body force in the displacement boundary integral equation into a boundary integral; S6：分别对每一个颗粒建立局部坐标系，并对相似颗粒计算转换矩阵；S6: respectively establish a local coordinate system for each particle, and calculate a transformation matrix for similar particles; S7：对单颗粒边界积分方程进行离散，划分边界单元，并在域内布置节点；S7: Discretize the single particle boundary integral equation, divide the boundary elements, and arrange nodes in the domain; S8：将每一个颗粒边界上的面力等效为边界单元上的面力；S8: Equivalent the surface force on each particle boundary to the surface force on the boundary element; S9：计算单颗粒的系数矩阵，并对系数矩阵求逆；S9: Calculate the coefficient matrix of a single particle, and invert the coefficient matrix; S10：求解边界节点的位移；S10: Solve the displacement of the boundary node; S11：求解域内的节点的应力和应变等值；S11: Solve the equivalent stress and strain of nodes in the domain; S12：根据内部应力值判断颗粒是否开裂，并利用颗粒替代法进行破碎颗粒的替代，并将代替颗粒的模型信息返回到颗粒相互作用求解器中；S12: Determine whether the particles are cracked according to the internal stress value, use the particle replacement method to replace the broken particles, and return the model information of the replaced particles to the particle interaction solver; 步骤S3包括：Step S3 includes: 对于作用在中心为O(x0,y0)颗粒上一点P(xp,yp)的集中力F(Fx,Fy)，可等效为均布力p，其中，
Figure FDA0003500888510000021
角度范围[θ1,θ2]为面力等效范围，且
Figure FDA0003500888510000022
(x0,y0)表示颗粒中心点的横纵坐标，(xp,yp)表示点P的横纵坐标，(Fx,Fy)表示集中力F的横纵轴的分解力，R表示颗粒半径；For a concentrated force F(F x , F y ) acting on a point P(x p , y p ) on a particle whose center is O(x 0 , y 0 ), it can be equivalent to a uniform force p, where,
Figure FDA0003500888510000021
The angular range [θ 1 , θ 2 ] is the equivalent range of the surface force, and
Figure FDA0003500888510000022
(x 0 , y 0 ) represents the horizontal and vertical coordinates of the particle center point, (x p , y p ) represents the horizontal and vertical coordinates of point P, (F x , F y ) represents the decomposition force of the horizontal and vertical axes of the concentration force F, R represents the particle radius; 由
Figure FDA0003500888510000023
建立控制方程，其中，x为研究域内的点，u为应力张量；a为加速度；ρ为颗粒材料密度；μ为剪切模量；λ＝2vμ/(1-2v)为拉美常数；v为泊松比，▽x()和▽x·()为散度运算符；
Figure FDA0003500888510000024
为合力矢量，b(x)为等效体力，S表示颗粒面积，FI表示颗粒的集中力项，n表示集中力个数；Depend on
Figure FDA0003500888510000023
Establish the governing equation, where x is the point in the research domain, u is the stress tensor; a is the acceleration; ρ is the density of the granular material; μ is the shear modulus; λ=2vμ/(1-2v) is the Latin American constant; v is Poisson’s ratio, ▽ x () and ▽ x ·() are divergence operators;
Figure FDA0003500888510000024
is the resultant force vector, b(x) is the equivalent physical force, S represents the particle area, F I represents the particle's concentrated force term, and n represents the number of concentrated forces; 对于形状相似的颗粒，只需要在局部坐标系内建立边界积分方程，其中，边界积分方程在局部坐标中可表达为：
Figure FDA0003500888510000025
s和
Figure FDA0003500888510000026
为局部坐标系中的点，
Figure FDA0003500888510000027
和
Figure FDA0003500888510000028
为局部坐标系中的位移和面力，Γ为问题研究域Ω的边界，
Figure FDA0003500888510000029
表示和点s位置相关的系数，
Figure FDA00035008885100000210
表示局部坐标系下的基本解函数，
Figure FDA00035008885100000211
表示局部坐标系下的边界，
Figure FDA00035008885100000212
表示局部坐标系下的基本解函数，
Figure FDA00035008885100000213
表示局部坐标系下点
Figure FDA00035008885100000214
的位移，
Figure FDA00035008885100000215
表示局部坐标系下的体力矢量，
Figure FDA00035008885100000216
表示局部坐标系下的求解域，
Figure FDA00035008885100000217
表示局部坐标系下点
Figure FDA00035008885100000218
和s的距离，δki表示克罗内克符号，
Figure FDA00035008885100000219
表示局部坐标系下r的偏导数，
Figure FDA00035008885100000220
表示局部坐标系下r的偏导数；For particles with similar shapes, it is only necessary to establish a boundary integral equation in the local coordinate system, where the boundary integral equation can be expressed in local coordinates as:
Figure FDA0003500888510000025
s and
Figure FDA0003500888510000026
is a point in the local coordinate system,
Figure FDA0003500888510000027
and
Figure FDA0003500888510000028
is the displacement and surface force in the local coordinate system, Γ is the boundary of the problem research domain Ω,
Figure FDA0003500888510000029
represents the coefficient related to the position of point s,
Figure FDA00035008885100000210
represents the basic solution function in the local coordinate system,
Figure FDA00035008885100000211
represents the boundary in the local coordinate system,
Figure FDA00035008885100000212
represents the basic solution function in the local coordinate system,
Figure FDA00035008885100000213
Represents a point in the local coordinate system
Figure FDA00035008885100000214
displacement,
Figure FDA00035008885100000215
represents the body force vector in the local coordinate system,
Figure FDA00035008885100000216
represents the solution domain in the local coordinate system,
Figure FDA00035008885100000217
Represents a point in the local coordinate system
Figure FDA00035008885100000218
the distance from s, δ ki denotes Kronecker notation,
Figure FDA00035008885100000219
represents the partial derivative of r in the local coordinate system,
Figure FDA00035008885100000220
represents the partial derivative of r in the local coordinate system; 步骤S8包括：Step S8 includes: 对于单元Zj[sj,sj+1]，sj和sj+1为单元的两个端点，j表示单元序号，如果满足：
Figure FDA0003500888510000031
其中，
Figure FDA0003500888510000032
为集中力等效的范围，
Figure FDA0003500888510000033
和
Figure FDA0003500888510000034
为离散边界上的两个端点，单元所受等效面力为：
Figure FDA0003500888510000035
p0为初始等效面力，
Figure FDA0003500888510000036
β1和β2分别为单元集的起点和终点角度，p1和p2为在第一个和最后一个单元上残余力引起的压力，可由以下两组公式求得：
Figure FDA0003500888510000037
及
Figure FDA0003500888510000038
β3和β4分别为第一个单元的终点和最后一个单元的起点与集中力的夹角；For unit Z j [s j ,s j+1 ], s j and s j+1 are the two endpoints of the unit, and j represents the unit number, if:
Figure FDA0003500888510000031
in,
Figure FDA0003500888510000032
is the equivalent range of concentrated force,
Figure FDA0003500888510000033
and
Figure FDA0003500888510000034
are the two endpoints on the discrete boundary, and the equivalent surface force on the element is:
Figure FDA0003500888510000035
p 0 is the initial equivalent surface force,
Figure FDA0003500888510000036
β 1 and β 2 are the starting and ending angles of the element set, respectively, and p 1 and p 2 are the pressures due to residual forces on the first and last elements, which can be obtained from the following two sets of formulas:
Figure FDA0003500888510000037
and
Figure FDA0003500888510000038
β 3 and β 4 are the angle between the end point of the first element and the start point of the last element and the concentrated force, respectively; 由
Figure FDA0003500888510000039
确定每一个颗粒的边界节点，其中，
Figure FDA00035008885100000310
表示边界上点的位移，H-1表示H的逆矩阵，H表示边界积分里形成的系数矩阵，G表示边界积分里形成的系数矩阵，
Figure FDA00035008885100000311
表示局部坐标下节点的面力，
Figure FDA00035008885100000312
表示局部坐标下节点的体力系数矩阵，
Figure FDA00035008885100000313
表示局部坐标下节点的体力；Depend on
Figure FDA0003500888510000039
Determine the boundary nodes of each particle, where,
Figure FDA00035008885100000310
represents the displacement of the point on the boundary, H -1 represents the inverse matrix of H, H represents the coefficient matrix formed in the boundary integral, G represents the coefficient matrix formed in the boundary integral,
Figure FDA00035008885100000311
represents the surface force of the node in local coordinates,
Figure FDA00035008885100000312
represents the physical strength coefficient matrix of nodes in local coordinates,
Figure FDA00035008885100000313
Represents the physical force of the node in local coordinates; 由
Figure FDA00035008885100000314
确定颗粒域内节点的应变，由
Figure FDA00035008885100000315
确定应力值，其中，G′表示应变离散矩阵方程中对应位移的系数矩阵，H′表示应变离散矩阵方程中对应面力的系数矩阵，
Figure FDA00035008885100000316
表示应变离散矩阵方程中对应体力的系数矩阵，
Figure FDA00035008885100000317
m表示需要求解的内部点的数量，
Figure FDA00035008885100000318
和
Figure FDA00035008885100000319
为每个点在x，y方向上的应变值，
Figure FDA00035008885100000320
表示局部坐标系下的应力，
Figure FDA00035008885100000321
表示局部坐标系下的刚度矩阵，
Figure FDA00035008885100000322
表示局部坐标系下的应变，
Figure FDA00035008885100000323
δij、δkl、δik、δjl、δil及δjk表示克罗内克符号，v为泊松比。Depend on
Figure FDA00035008885100000314
Determine the strain at the nodes within the particle domain, given by
Figure FDA00035008885100000315
Determine the stress value, where G′ represents the coefficient matrix of the corresponding displacement in the strain discrete matrix equation, H′ represents the coefficient matrix of the corresponding surface force in the strain discrete matrix equation,
Figure FDA00035008885100000316
represents the coefficient matrix of the corresponding body force in the discrete matrix equation of strain,
Figure FDA00035008885100000317
m represents the number of interior points to be solved,
Figure FDA00035008885100000318
and
Figure FDA00035008885100000319
is the strain value of each point in the x, y direction,
Figure FDA00035008885100000320
represents the stress in the local coordinate system,
Figure FDA00035008885100000321
represents the stiffness matrix in the local coordinate system,
Figure FDA00035008885100000322
represents the strain in the local coordinate system,
Figure FDA00035008885100000323
δ ij , δ kl , δ ik , δ jl , δ il and δ jk represent Kronecker symbols, and v is Poisson's ratio. 2.根据权利要求1所述的方法，其特征在于，步骤S5包括：2. The method according to claim 1, wherein step S5 comprises: 采用直线积分法进行域积分计算，其中，域积分表示为：D1＝∫ΩU(x,y)bdΩ(y)，基于直线积分法，由
Figure FDA0003500888510000041
将域积分转化为等效边界积分，其中，M表示积分线个数，Wi表示权重因子，Li表示积分线，U(x,y)表示格林函数基本解，b表示等效体力，dy1表示微分。The domain integral calculation is carried out by the straight-line integral method, where the domain integral is expressed as: D 1 =∫ Ω U(x,y)bdΩ(y), based on the straight-line integral method, by
Figure FDA0003500888510000041
Convert the domain integral to the equivalent boundary integral, where M is the number of integral lines, Wi is the weight factor, Li is the integral line, U(x, y) is the basic solution of Green's function, b is the equivalent physical strength, dy 1 means differential. 3.根据权利要求1所述的方法，其特征在于，步骤S6包括：3. The method according to claim 1, wherein step S6 comprises: 对于全局坐标中的点x(x1,x2)以及局部坐标中的点si(s1,s2)，可通过以下表达式进行转化：
Figure FDA0003500888510000042
Figure FDA0003500888510000043
为Jacobi矩阵的常系数，
Figure FDA0003500888510000044
|J|＝c2,
Figure FDA0003500888510000045
c为比例因子，逆矩阵为：
Figure FDA0003500888510000046
For a point x(x 1 ,x 2 ) in global coordinates and a point s i (s 1 ,s 2 ) in local coordinates, the transformation can be done by the following expressions:
Figure FDA0003500888510000042
Figure FDA0003500888510000043
are the constant coefficients of the Jacobi matrix,
Figure FDA0003500888510000044
|J|=c 2 ,
Figure FDA0003500888510000045
c is the scale factor, and the inverse matrix is:
Figure FDA0003500888510000046
 4.一种大规模颗粒材料内部应力及破碎模拟分析装置，其特征在于，包括：粒相互作用求解器和颗粒内部应力求解器；4. A large-scale particle material internal stress and crushing simulation analysis device, characterized in that, comprising: a particle interaction solver and a particle internal stress solver; 所述颗粒相互作用求解器，用于进行大规模颗粒集合体相互作用的模拟计算；The particle interaction solver is used for the simulation calculation of large-scale particle aggregate interaction; 所述颗粒内部应力求解器包括：输入模块，其用于接收操作者输入的待求解颗粒的材料参数、网格化分数、网格类型以及边界条件信息；求解模块，其用于采用权利要求1至3任意一项大规模颗粒材料内部应力及破碎模拟分析方法得到待求解颗粒体内部的应力。The particle internal stress solver includes: an input module, which is used for receiving the material parameters, meshing fraction, mesh type and boundary condition information of the particle to be solved input by an operator; a solving module, which is used for adopting claim 1 Any one of the internal stress and crushing simulation analysis methods of large-scale granular materials to 3 can obtain the internal stress of the particle body to be solved.
CN202010917113.3A 2020-09-03 2020-09-03 Large-scale granular material internal stress and crushing simulation analysis method and device Active CN112084647B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010917113.3A CN112084647B (en) 2020-09-03 2020-09-03 Large-scale granular material internal stress and crushing simulation analysis method and device

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010917113.3A CN112084647B (en) 2020-09-03 2020-09-03 Large-scale granular material internal stress and crushing simulation analysis method and device

Publications (2)

Publication Number Publication Date
CN112084647A CN112084647A (en) 2020-12-15
CN112084647B true CN112084647B (en) 2022-04-01

Family

ID=73731467

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010917113.3A Active CN112084647B (en) 2020-09-03 2020-09-03 Large-scale granular material internal stress and crushing simulation analysis method and device

Country Status (1)

Country Link
CN (1) CN112084647B (en)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117540600B (en) * 2023-11-17 2025-04-15 水电水利规划设计总院 A macroscopic step-by-step coupled multi-scale stress-deformation calculation method for rockfill dams
CN117540601A (en) * 2023-11-17 2024-02-09 水电水利规划设计总院 A method for solving the crushing path of rockfill material

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7265752B2 (en) * 2004-01-09 2007-09-04 Microsoft Corporation Multi-chart geometry images
CN102628861B (en) * 2012-04-11 2014-10-22 武汉大学 Method for simulating temperature cracking value of mass concrete
CN106980735B (en) * 2017-04-06 2020-07-28 西南科技大学 Numerical Simulation Method for Thermal Fracture of Brittle Materials
CN109284523A (en) * 2018-07-19 2019-01-29 同济大学 A simulation method for progressive failure and solid-liquid phase transition behavior of geotechnical media
CN110598293B (en) * 2019-09-03 2020-05-05 上海交通大学 Method for predicting fracture damage behavior of micro-nano fiber reinforced composite material

Also Published As

Publication number Publication date
CN112084647A (en) 2020-12-15

Similar Documents

Publication Publication Date Title
Schneider et al. A large-scale comparison of tetrahedral and hexahedral elements for solving elliptic PDEs with the finite element method
CN110555263B (en) Level set topology optimization method for curved shell structure optimization design
Vincent et al. Facilitating the Adoption of Unstructured High-Order MethodsAmongst a Wider Community of Fluid Dynamicists
CN112084647B (en) Large-scale granular material internal stress and crushing simulation analysis method and device
CN103838913B (en) The Finite Element of the curved bridge of curved box girder
CN106354918B (en) The construction method of fluid structurecoupling problem numerical simulation in a kind of hydraulic fracturing
Gharti et al. Simulation of multistage excavation based on a 3D spectral-element method
CN102298660A (en) Universal method of boundary modeling based on distinct element method
Feng et al. An edge/face-based smoothed radial point interpolation method for static analysis of structures
CN108170898A (en) A kind of jointed rock slope reliability analysis Lower Bound Limit
CN110532727A (en) It can be used for the method for numerical simulation of common non-newtonian fluid
CN104281730A (en) Great-rotating-deformation plate shell structure dynamic response finite element analysis method
Tong et al. Three-dimensional numerical manifold method for heat conduction problems with a simplex integral on the boundary
Liang et al. Extended material point method for the three‐dimensional crack problems
Radu et al. A doubly nonlocal Laplace operator and its connection to the classical Laplacian
Tang et al. An efficient adaptive analysis procedure using the edge-based smoothed point interpolation method (ES-PIM) for 2D and 3D problems
Tate et al. Fixed-topology Lorentzian triangulations: quantum Regge calculus in the Lorentzian domain
Shang et al. An efficient 4‐node facet shell element for the modified couple stress elasticity
Li et al. A staggered asynchronous step integration algorithm for hybrid finite-element and discrete-element modeling
CN117540600B (en) A macroscopic step-by-step coupled multi-scale stress-deformation calculation method for rockfill dams
He et al. Simulation of solids with multiple rectangular inhomogeneities using non-uniform eigenstrain formulation of BIEs
CN114818523A (en) A FDM-DDA-Based Calculation Method for Unsaturated Transient Fluid-Structure Interaction
CN107220212A (en) A kind of boundary Element method of two-dimentional non-compact border sound scattering
Lv et al. Extended multiscale finite element method based on polyhedral coarse grid elements for heterogeneous materials and structures
Schneidera et al. A Large Scale Comparison of Tetrahedral and Hexahedral Elements for Finite Element Analysis

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant