CN112084647A - 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 PDFInfo
- Publication number
- CN112084647A CN112084647A CN202010917113.3A CN202010917113A CN112084647A CN 112084647 A CN112084647 A CN 112084647A CN 202010917113 A CN202010917113 A CN 202010917113A CN 112084647 A CN112084647 A CN 112084647A
- Authority
- CN
- China
- Prior art keywords
- particle
- boundary
- coordinate system
- representing
- force
- 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.)
- Granted
Links
- 238000004088 simulation Methods 0.000 title claims abstract description 30
- 238000004458 analytical method Methods 0.000 title claims abstract description 17
- 239000008187 granular material Substances 0.000 title claims description 14
- 239000002245 particle Substances 0.000 claims abstract description 141
- 239000011159 matrix material Substances 0.000 claims abstract description 62
- 238000000034 method Methods 0.000 claims abstract description 51
- 238000004364 calculation method Methods 0.000 claims abstract description 25
- 239000000463 material Substances 0.000 claims abstract description 20
- 238000006243 chemical reaction Methods 0.000 claims abstract description 10
- 238000011160 research Methods 0.000 claims abstract description 4
- 230000003993 interaction Effects 0.000 claims description 27
- 238000006073 displacement reaction Methods 0.000 claims description 22
- 230000010354 integration Effects 0.000 claims description 18
- 239000000243 solution Substances 0.000 claims description 11
- 239000013598 vector Substances 0.000 claims description 8
- 239000003637 basic solution Substances 0.000 claims description 7
- 239000011236 particulate material Substances 0.000 claims description 5
- 230000009471 action Effects 0.000 claims description 4
- 230000001133 acceleration Effects 0.000 claims description 3
- 238000004613 tight binding model Methods 0.000 claims description 3
- 230000004304 visual acuity Effects 0.000 claims description 3
- 238000005229 chemical vapour deposition Methods 0.000 claims 1
- 229910044991 metal oxide Inorganic materials 0.000 claims 1
- 150000004706 metal oxides Chemical class 0.000 claims 1
- 238000010168 coupling process Methods 0.000 abstract description 5
- 238000011439 discrete element method Methods 0.000 abstract description 4
- 230000008878 coupling Effects 0.000 abstract 1
- 238000005859 coupling reaction Methods 0.000 abstract 1
- 230000006870 function Effects 0.000 description 7
- 238000013459 approach Methods 0.000 description 3
- 239000000203 mixture Substances 0.000 description 3
- 230000008901 benefit Effects 0.000 description 2
- 239000004566 building material Substances 0.000 description 2
- 230000008859 change Effects 0.000 description 2
- 238000000354 decomposition reaction Methods 0.000 description 2
- 230000005484 gravity Effects 0.000 description 2
- 230000006872 improvement Effects 0.000 description 2
- 230000007246 mechanism Effects 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 239000004576 sand Substances 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 230000009977 dual effect Effects 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 238000011065 in-situ storage Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000008569 process Effects 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
- 230000001131 transforming effect Effects 0.000 description 1
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
The invention discloses a large-scale particle material internal stress and crushing simulation analysis method and device, belongs to the technical field of particle crushing, and is based on the basic idea of continuous discrete coupling, a boundary element method and a discrete element method are combined to carry out large-scale particle crushing research, the internal stress of particles can be calculated and analyzed by using the boundary element method, and whether the particles are crushed or not is judged by combining the mechanical fracture theory of continuous media and the HokkBrown criterion. Further, in the case of internal stress simulation using boundary elements, for a collection of particles having similar shapes, for example: the invention only needs to calculate the coefficient matrix of one particle, and other similar particles obtain the coefficient matrix through coordinate conversion and coefficient scaling, thereby greatly improving the calculation efficiency. The invention also performs internal stress simulation of the circular rockfill material and proves the effectiveness of the internal stress simulation.
Description
Technical Field
The invention belongs to the technical field of particle crushing, and particularly relates to a large-scale particle material internal stress and crushing simulation analysis method and device.
Background
The granular material is used as an important building material and widely applied to engineering, for example, the main filling material of a rock-fill dam is a rock-fill body, and the material of sand and gravel commonly seen in civil engineering is used as a building material. Due to the action of external force, the interior of the particle material is often in a high stress state, so that the particles are broken, the gradation of the particle aggregate is changed, and the macro-micro mechanical property of the particle aggregate is influenced. It is therefore essential to study the mechanical mechanism of particle breakage. However, due to the influence of the size and the position of the granular material, in-situ mechanical mechanism analysis of the granular material is difficult to reduce through physical experiments, so that the numerical simulation method of the granular material is popular. The continuous discrete coupling method is a scientific particle breakage numerical analysis means, is mainly based on a finite element method, and has very high calculation cost.
The existing discrete element method generally treats particles as rigid bodies when simulating the interaction of particle aggregates, and cannot calculate the internal stress of the particles, so that the judgment criterion of particle breakage by adopting a maximum contact force principle or setting a stress threshold is obtained by experience, and the theoretical basis is weak. The existing continuous discrete coupling method, such as a finite element-discrete element coupling method and a proportional boundary finite element-discrete element coupling method, can utilize a mature fracture theory in continuous medium mechanics to judge the particle breakage, but the calculation efficiency is generally low.
For the domain integration caused by gravity and the like, further processing is needed, otherwise, the boundary element loses the advantage of dimensionality reduction, and the currently adopted domain integration processing method mainly comprises the following steps: dual Reciprocity Method (DRM), in which a heterogeneous term can be approximated by a series of functions such as Radial Basis Function (RBF), and a second Reciprocity is applied to convert a domain integral to a boundary integral. Only points on the domain or boundary need to provide information represented by the non-homogeneous term. However, the accuracy of DRM depends to a large extent on the distribution and location of the domain points, as well as the type of function used to approximate the non-homogeneous term. Furthermore, the arrangement of points in a complex domain may not be easy, especially for three-dimensional problems. A DRM-like approach is the Multiple Reciprocity Method (MRM), in which the Reciprocity theorem is applied iteratively through a series of high-order fundamental solutions, transforming the domain integrals into boundaries. Another approach is the special solution approach (PSM), in which an approximate special solution is constructed instead of performing a second reciprocity in DRM to convert the domain integrals into boundary integrals.
Disclosure of Invention
Aiming at the defects or the improvement requirements of the prior art, the invention provides a large-scale simulation analysis method and device for the internal stress and crushing of the granular material, and the method and device have extremely high accuracy and effectiveness.
To achieve the above object, according to one aspect of the present invention, there is provided a large-scale simulation analysis method for internal stress and fracture of a particulate material, comprising:
s1: carrying out simulation calculation of large-scale particle aggregate interaction by using a particle interaction solver;
s2: deriving a calculation result of a certain time step in a particle interaction solver, establishing a geometric model of particles to be solved by utilizing modeling software in a particle internal stress solver, inputting material parameters, gridding fractions, grid types and boundary conditions of the particle model to be solved derived from the particle interaction solver, and outputting the internal stress of the geometric model of the particles to be solved by utilizing a boundary element method;
s3: the concentration force on each particle boundary is equivalent to the surface force on the boundary;
s4: establishing a control equation, and obtaining a displacement boundary integral equation by the control equation;
s5: converting domain integration caused by volume force in a displacement boundary integral equation into boundary integration;
s6: respectively establishing a local coordinate system for each particle, and calculating a conversion matrix for similar particles;
s7: dispersing a single-grain boundary integral equation, dividing boundary units, and arranging nodes in a domain;
s8: equating the surface force on each grain boundary to the surface force on the boundary unit;
s9: calculating a coefficient matrix of the single particle, and inverting the coefficient matrix;
s10: solving the displacement of the boundary node;
s11: solving the stress and strain equivalence of nodes in the domain;
s12: and judging whether the particles crack or not according to the internal stress value, replacing the broken particles by using a particle replacement method, and returning the model information of the replaced particles to the particle interaction solver.
In some alternative embodiments, step S3 includes:
for action centered on O (x)0,y0) A point P (x) on the particlep,yp) Concentration force F (F)x,Fy) It can be equivalent to a uniform force p, wherein,
angular range [ theta ]1,θ2]Is the equivalent range of surface force, and
(x0,y0) The abscissa and ordinate of the center point of the particle, (x)p,yp) Represents the abscissa and ordinate of the point P, (F)x,Fy) The resolving power is shown transverse to the longitudinal axis of the concentration force F, and R is the particle radius.
In some alternative embodiments, the composition is prepared by
Establishing a control equation, wherein x is a point in a research domain, and u is a stress tensor; a is the acceleration; ρ is the density of the particulate material; μ is the shear modulus; λ ═ 2 ν μ/(1-2v) for the lame constant; v is the poisson's ratio and,
and
is the divergence operator;
is the resultant force vector, b (x) is the equivalent physical force, S represents the particle area, FIThe term of the concentration force of the particles is shown, and n represents the number of the concentration forces.
In some alternative embodiments, step S5 includes:
and performing domain integral calculation by adopting a straight line integral method, wherein the domain integral is expressed as: d1＝∫ΩU (x, y) bd omega (y) based on a linear integration method
Converting the domain integral into equivalent boundary integral, wherein M represents the number of integral lines, WiRepresents a weight factor, LiRepresents the integral line, U (x, y) represents the fundamental solution of the Green function, b represents the equivalent physical strength，dy1The differential is indicated.
In some alternative embodiments, step S6 includes:
for point x (x) in global coordinates1,x2) And a point s in local coordinatesi(s1,s2) The conversion can be performed by the following expression:
being the constant coefficients of the Jacobi matrix,
|J|＝c2,
c is a scale factor and the inverse matrix is:
in some alternative embodiments, for similarly shaped particles, only a boundary integral equation needs to be established in the local coordinate system, wherein the boundary integral equation can be expressed in the local coordinates as:
s and
are points in the local coordinate system and are,
and
the boundaries of the domain omega are studied for the problem, for displacements and surface forces in the local coordinate system,
the coefficient relating to the position of the point s is represented,
representing the basic solution function in a local coordinate system,
which represents the boundary under the local coordinate system,
representing the basic solution function in a local coordinate system,
representing points in a local coordinate system
Is detected by the displacement of (a) a,
representing the physical force vector in a local coordinate system,
representing the solution domain in the local coordinate system,
representing points in a local coordinate system
And the distance of s is such that,kirepresents a symbol of the form of a kronecker,
representing the partial derivative of r in the local coordinate system,
representing the partial derivative of r in the local coordinate system.
In some alternative embodiments, step S8 includes:
for unit Zj[sj,sj+1]，sjAnd sj+1J represents the unit number for both ends of the unit, if:
wherein,
in order to be the range of equivalence of the concentration force,
and
for two endpoints on a discrete boundary, the equivalent surface force experienced by the cell is:
p0in order to be the initial equivalent surface force,
β1and beta2Respectively start and end angles, p, of the cell set1And p2The pressure caused by the residual force on the first and last cells can be determined by the following two equations:
and
β3and beta4The included angles between the end point of the first unit and the start point of the last unit and the concentrated force are respectively.
In some alternative embodiments, the composition is prepared by
The boundary nodes for each grain are determined, wherein,
representing boundariesDisplacement of upper point, H-1An inverse matrix representing H, a coefficient matrix formed in the boundary integral, G a coefficient matrix formed in the boundary integral,
representing the surface force of the node at local coordinates,
a physical coefficient matrix representing the nodes at local coordinates,
representing the physical strength of the node in local coordinates.
In some alternative embodiments, the composition is prepared by
Determining the strain of a node in the particle domain from
Determining stress values, wherein G 'represents a coefficient matrix corresponding to displacement in the strain discrete matrix equation, H' represents a coefficient matrix corresponding to surface force in the strain discrete matrix equation,
a matrix of coefficients representing the corresponding physical forces in the strain discrete matrix equation,
m represents the number of interior points that need to be solved,
and
for each point strain value in the x, y direction,
representing stress in a local coordinate system，
Representing a stiffness matrix in a local coordinate system,
which represents the strain in the local coordinate system,
ij、kl、ik、jl、ilandjkrepresenting the kronecker notation, v is the poisson ratio.
According to another aspect of the present invention, there is provided a large scale grain material internal stress and fracture simulation analysis apparatus, comprising: a particle interaction solver and a particle internal stress solver;
the particle interaction solver is used for carrying out simulation calculation on the interaction of the large-scale particle aggregate;
the intra-particle stress solver comprises: the input module is used for receiving the material parameters, the gridding fraction, the grid type and the boundary condition information of the particles to be solved, which are input by an operator; and the solving module is used for obtaining the stress inside the particle body to be solved by adopting any one of the large-scale particle material internal stress and the crushing simulation analysis method.
In general, compared with the prior art, the above technical solution contemplated by the present invention can achieve the following beneficial effects:
according to the large-scale granular material internal stress and crushing simulation analysis method and device provided by the invention, the internal stress of granular materials such as sand gravel and the like can be easily solved, the judgment of the crushing of the granules can be carried out, and the gradation change and the change of macro-micro mechanical properties of the granular materials can be further grasped, so that the safety and the reliability of an engineering structure can be ensured.
Drawings
FIG. 1 is a schematic flow chart of a simulation method for internal stress and fracture of a particle assembly according to an embodiment of the present invention;
FIG. 2 is a schematic diagram of a model and its boundary conditions according to an embodiment of the present invention;
fig. 3 is a cloud image of internal stresses calculated by a simulation apparatus according to an embodiment of the present invention.
Detailed Description
In order to make the objects, technical solutions and advantages of the present invention more apparent, the present invention is described in further detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and are not intended to limit the invention. In addition, the technical features involved in the embodiments of the present invention described below may be combined with each other as long as they do not conflict with each other.
According to the invention, a boundary element method and a discrete element method are combined to carry out large-scale particle crushing research, the internal stress of the particles can be calculated and analyzed by using the boundary element method, and whether the particles are crushed or not is judged by combining a continuous medium mechanical fracture theory and a HockBrown criterion. Further, in the case of internal stress simulation using boundary elements, for a collection of particles having similar shapes, for example: the invention only needs to calculate the coefficient matrix of one particle, and other similar particles obtain the coefficient matrix through coordinate conversion and coefficient scaling, thereby greatly improving the calculation efficiency. The invention also performs internal stress simulation of the circular rockfill material and proves the effectiveness of the internal stress simulation. The method can realize the mutual conversion of the similar particle coefficient matrixes, and each similar particle does not need to be calculated, so the calculation efficiency is further improved. For the concentration force, the direct processing by adopting the boundary element has certain difficulty, and the boundary element needs to be equivalent to the boundary unit, and the invention provides a two-step conversion method: the first step is to convert the concentration force into surface force evenly distributed on the boundary; the second step translates the surface force on the boundary into a surface force on the boundary element. After the equivalent method is adopted, on one hand, the singularity caused by the concentration force is eliminated, on the other hand, the coefficient matrix is only needed to be calculated once for the units with the same shape, and the boundary conditions are not needed to be processed by repeatedly carrying out boundary integration under different boundary conditions, so that the calculation amount can be greatly saved. The invention adopts a Line Lim Method (LIM) to process a boundary element domain integration Method, in the linear integration Method, only a boundary unit is used for dispersing the boundary, and a volume unit is not needed. Therefore, the method can be regarded as a boundary discretization method. In line integration, the domain integral can be computed as the sum of the one-dimensional integrals on the line, and the integral line can be easily created from the boundary elements and the boundary elements in the OCT tree structure.
As shown in fig. 1, a large-scale simulation analysis method for internal stress and fracture of a granular material provided by an embodiment of the present invention includes the following steps:
step S1: carrying out simulation calculation of large-scale particle aggregate interaction by using a particle interaction solver;
the particle interaction solver represents a module for solving the interaction between particles based on a discrete element method.
In the embodiment of the present invention, as shown in fig. 2, to verify the accuracy and effectiveness, a calculation example of 292 circular rockfill materials subjected to a top-down pressure is selected, and the elastic modulus is set to be 1000Mpa, and the poisson ratio is set to be 0.25.
The large scale in the examples of the present invention preferably refers to a scale of more than 100 particles.
Step S2: exporting a calculation result of a certain time step in the particle interaction solver, wherein the calculation result mainly comprises geometric information, material properties, boundary conditions and the like of particles; in a particle internal stress solver, establishing a geometric model of particles to be solved by utilizing modeling software, inputting information such as material parameters, gridding fractions, grid types and boundary conditions of the particle model to be solved, which is derived from the particle interaction solver, and then outputting internal stress of the particle geometric model to be solved;
the particle internal stress solver represents a module for solving the particle internal stress based on a boundary element method.
In the embodiment of the present invention, the stress-strain distribution of the particulate material is calculated by using a boundary element method, and the basic equation of the boundary element method is as follows:
wherein the boundary of the domain Ω is investigated for the problem; x and y represent source and field points within the problem domain Ω; t is the surface force on the boundary; u (x, y) and T (x, y) are basic solutions, which can be written as:
wherein r represents the distance between the source point and the field point; n is the unit external normal vector of the boundary.
Step S3: the concentration force on each particle boundary is equivalent to the surface force on the boundary;
in the embodiment of the invention, the concentration force on the single particle boundary is equivalent to the surface force on the boundary by adopting an equivalent method, and the equivalent method comprises the following steps:
for action centered on O (x)0,y0) A point P (x) on the particlep,yp) Concentration force F (F)x,Fy) Equivalent to a uniform force p, the calculation expression is:
angle range [ theta ] therein1,θ2]For the area force equivalent range, the following formula can be used for calculation:
wherein (x)0,y0) The abscissa and ordinate of the center point of the particle, (x)p,yp) Represents the abscissa and ordinate of the point P, (F)x,Fy) Indicating a concentration forceThe resolving power of F transverse to the longitudinal axis, R represents the particle radius.
Step S4: establishing a control equation:
wherein x is a point within the study domain and u is the stress tensor; a is the acceleration; ρ is the density of the particulate material; μ is the shear modulus; λ ═ 2 ν μ/(1-2v) for the lame constant; v is the poisson's ratio and,
and
is the divergence operator;
is the resultant force vector, and b (x) is the equivalent physical force. S represents the particle area, FIThe term of the concentration force of the particles is shown, and n represents the number of the concentration forces.
Step S5: establishing a displacement boundary integral equation by a control equation;
where the displacement boundary integral equation is used to solve for the displacement of the grain and the stress on the grain boundary.
Step S6: converting domain integration caused by volume force in a displacement boundary integral equation into boundary integration;
in the embodiment of the invention, the intra-domain integral generated due to gravity and the like is calculated by adopting a linear integration method:
for domain integration:
D1＝∫ΩU(x,y)bdΩ(y) (7)
based on a linear integration method, the method is converted into equivalent boundary integration by using the following formula:
wherein M represents the number of integral lines, WiRepresents a weight factor, LiRepresents the integral line, U (x, y) represents the fundamental solution of the Green function, b represents the equivalent physical power, dy1The differential is indicated.
Step S7: respectively establishing a local coordinate system for each particle, and calculating a conversion matrix for similar particles;
similar particles in embodiments of the present invention refer to particles of infinitely close size.
In the embodiment of the invention, a local coordinate system is established for each particle, and a transformation matrix is established for similar particles, and the method mainly comprises the following steps:
for point x (x) in global coordinates1,x2) And a point s in local coordinatesi(s1,s2) The conversion can be performed by the following expression:
wherein,
is the constant coefficient of Jacobi matrix, and
wherein,
|J|＝c2,
c is a scale factor.
The inverse matrix is:
step S8: dispersing a single-grain boundary integral equation, dividing boundary units, and arranging nodes in a domain;
in the embodiment of the present invention, for particles with similar shapes, only a boundary integral equation needs to be established in the local coordinate system:
the boundary integral equation can be expressed in local coordinates as:
s and
are points in the local coordinate system and are,
and
for displacements and surface forces in the local coordinate system,
as shown in formula (3), and
in which the boundaries of the domain omega are investigated for the problem,
the coefficient relating to the position of the point s is represented,
representing the basic solution function in a local coordinate system,
which represents the boundary under the local coordinate system,
representing the basic solution function in a local coordinate system,
representing points in a local coordinate system
Is detected by the displacement of (a) a,
representing the physical force vector in a local coordinate system,
representing the solution domain in the local coordinate system,
representing points in a local coordinate system
And the distance of s is such that,kirepresents a symbol of the form of a kronecker,
representing the partial derivative of r in the local coordinate system,
representing the partial derivative of r in the local coordinate system.
For particles with similar shapes, only a coefficient matrix of a first-order boundary integral equation needs to be solved, and the inverse of the coefficient matrix is solved by adopting a singular value decomposition method.
Step S9: equating the surface force on each grain boundary to the surface force on the boundary unit;
in the embodiment of the invention, an equivalent method is adopted to equate the surface force on the single grain boundary to the surface force on the boundary unit, and the equivalent method is as follows:
for unit Zj[sj,sj+1]，sjAnd sj+1Is the two ends of a cell, j represents the cell number, if fullFoot:
wherein,
in order to be the range of equivalence of the concentration force,
and
for two endpoints on a discrete boundary, the equivalent surface force experienced by the cell is:
wherein p is0The calculation expression for the initial equivalent surface force is as follows:
wherein, beta1And beta2Respectively start and end angles, p, of the cell set1And p2The pressure caused by the residual force on the first and last cells can be determined by the following two equations:
wherein, beta3And beta4The included angles between the end point of the first unit and the start point of the last unit and the concentrated force are respectively.
Step S10: calculating a coefficient matrix of the single particle, and inverting the coefficient matrix;
step S11: solving the displacement of the boundary node;
step S12: solving the stress and strain equivalence of nodes in the domain;
in the embodiment of the invention, for all the grains with similar shapes, after the coefficient matrix and the inverse matrix of the next basic grain in the local coordinate system are only calculated, the unknown quantities of all the boundary nodes of the similar grains and the nodes in the domain can be calculated by matrix multiplication, so that the calculation time can be greatly saved:
the boundary node for each grain can be calculated by the following formula:
wherein,
representing the displacement of points on the boundary, H-1An inverse matrix representing H, a coefficient matrix formed in the boundary integral, G a coefficient matrix formed in the boundary integral,
representing the surface force of the node at local coordinates,
a physical coefficient matrix representing the nodes at local coordinates,
representing the physical strength of the node in local coordinates.
Due to the matrix H-1G and
the efficiency of equation (19) can be very high and is suitable for large scale particle calculation, requiring only one calculation. Since H is a singular matrix, its inverse can passSingular Value Decomposition (SVD) was calculated as follows:
H-1＝V2∑+V1 T (20)
wherein, V1And V2Left and right singular vectors, and ∑+For the generalized inverse of Σ, Σ can be written as:
therein, sigma0For the first N-3 singular values (arranged from large to small), the smallest three are set to zero. Rigid body displacement of the particles can be prevented and higher accuracy is obtained, N representing the number of matrix rows.
Furthermore, the strain for a node within the granular domain may be calculated as follows:
wherein G 'represents a coefficient matrix corresponding to displacement in the strain discrete matrix equation, H' represents a coefficient matrix corresponding to surface force in the strain discrete matrix equation,
and the coefficient matrix represents the corresponding physical force in the strain discrete matrix equation.
Wherein:
where m represents the number of interior points that need to be solved,
and
strain values in the x, y direction for each point.
After the strain value is calculated, the stress value may be calculated as follows:
wherein,
which represents the stress in a local coordinate system,
representing a stiffness matrix in a local coordinate system,
representing the strain in the local coordinate system.
And
wherein,ij、kl、ik、jl、ilandjkrepresenting the kronecker notation, v is the poisson ratio.
As shown in fig. 3, fig. 3 is a stress cloud of the inside of the particle assembly, in which fig. 3 (a) shows a stress cloud of stress sxx, and fig. 3 (b) shows a stress cloud of stress syy.
Step S13: and judging whether the particles crack or not according to the stress value, replacing the broken particles by using a particle replacement method, and returning the model information of the replaced particles to the particle interaction solver.
The application also provides a large-scale granular material internal stress and broken simulation analysis device, includes: a particle interaction solver and a particle internal stress solver;
the particle interaction solver is used for carrying out simulation calculation on the interaction of the large-scale particle aggregate;
the intra-particle stress solver comprises: the input module is used for receiving the material parameters, the gridding fraction, the grid type and the boundary condition information of the particles to be solved, which are input by an operator; and the solving module is used for obtaining the stress inside the particle body to be solved by adopting the large-scale particle material internal stress and crushing simulation analysis method.
The detailed description may refer to the description of the method embodiments, and the embodiments of the present invention will not be repeated.
It should be noted that, according to the implementation requirement, each step/component described in the present application can be divided into more steps/components, and two or more steps/components or partial operations of the steps/components can be combined into new steps/components to achieve the purpose of the present invention.
It will be understood by those skilled in the art that the foregoing is only a preferred embodiment of the present invention, and is not intended to limit the invention, and that any modification, equivalent replacement, or improvement made within the spirit and principle of the present invention should be included in the scope of the present invention.
Claims (10)
1. A large-scale particle material internal stress and crushing simulation analysis method is characterized by comprising the following steps:
s1: carrying out simulation calculation of large-scale particle aggregate interaction by using a particle interaction solver;
s2: deriving a calculation result of a certain time step in a particle interaction solver, establishing a geometric model of particles to be solved by utilizing modeling software in a particle internal stress solver, inputting material parameters, gridding fractions, grid types and boundary conditions of the particle model to be solved derived from the particle interaction solver, and outputting the internal stress of the geometric model of the particles to be solved by utilizing a boundary element method;
s3: the concentration force on each particle boundary is equivalent to the surface force on the boundary;
s4: establishing a control equation, and obtaining a displacement boundary integral equation by the control equation;
s5: converting domain integration caused by volume force in a displacement boundary integral equation into boundary integration;
s6: respectively establishing a local coordinate system for each particle, and calculating a conversion matrix for similar particles;
s7: dispersing a single-grain boundary integral equation, dividing boundary units, and arranging nodes in a domain;
s8: equating the surface force on each grain boundary to the surface force on the boundary unit;
s9: calculating a coefficient matrix of the single particle, and inverting the coefficient matrix;
s10: solving the displacement of the boundary node;
s11: solving the stress and strain equivalence of nodes in the domain;
s12: and judging whether the particles crack or not according to the internal stress value, replacing the broken particles by using a particle replacement method, and returning the model information of the replaced particles to the particle interaction solver.
2. The method according to claim 1, wherein step S3 includes:
for action centered on O (x)0,y0) A point P (x) on the particlep,yp) Concentration force F (F)x,Fy) It can be equivalent to a uniform force p, wherein,
angular range [ theta ]1,θ2]Is the equivalent range of surface force, and
(x0,y0) The abscissa and ordinate of the center point of the particle, (x)p,yp) Represents the abscissa and ordinate of the point P, (F)x,Fy) The resolving power is shown transverse to the longitudinal axis of the concentration force F, and R is the particle radius.
3. The method of claim 1, wherein the step of removing the metal oxide is performed by a chemical vapor deposition processIn that, by
Establishing a control equation, wherein x is a point in a research domain, and u is a stress tensor; a is the acceleration; ρ is the density of the particulate material; μ is the shear modulus; λ ═ 2 ν μ/(1-2v) for the lame constant; v is Poisson ratio +x() And +x() is the divergence operator;
is the resultant force vector, b (x) is the equivalent physical force, S represents the particle area, FIThe term of the concentration force of the particles is shown, and n represents the number of the concentration forces.
4. The method according to claim 1, wherein step S5 includes:
and performing domain integral calculation by adopting a straight line integral method, wherein the domain integral is expressed as: d1＝∫ΩU (x, y) bd omega (y) based on a linear integration method
Converting the domain integral into equivalent boundary integral, wherein M represents the number of integral lines, WiRepresents a weight factor, LiRepresents the integral line, U (x, y) represents the fundamental solution of the Green function, b represents the equivalent physical power, dy1The differential is indicated.
5. The method according to claim 1, wherein step S6 includes:
for point x (x) in global coordinates1,x2) And a point s in local coordinatesi(s1,s2) The conversion can be performed by the following expression:
is JacoThe constant coefficients of the bi-matrix are,
|J|＝c2,
c is a scale factor and the inverse matrix is:
6. the method of claim 1, wherein for similarly shaped particles, only a boundary integral equation needs to be established in the local coordinate system, wherein the boundary integral equation can be expressed in the local coordinates as:
s and
are points in the local coordinate system and are,
and
the boundaries of the domain omega are studied for the problem, for displacements and surface forces in the local coordinate system,
the coefficient relating to the position of the point s is represented,
representing the basic solution function in a local coordinate system,
which represents the boundary under the local coordinate system,
representing the basic solution function in a local coordinate system,
representing points in a local coordinate system
Is detected by the displacement of (a) a,
representing the physical force vector in a local coordinate system,
representing the solution domain in the local coordinate system,
representing points in a local coordinate system
And the distance of s is such that,kirepresents a symbol of the form of a kronecker,
representing the partial derivative of r in the local coordinate system,
representing the partial derivative of r in the local coordinate system.
7. The method according to claim 1, wherein step S8 includes:
for unit Zj[sj,sj+1]，sjAnd sj+1J represents the unit number for both ends of the unit, if:
wherein,
in order to be the range of equivalence of the concentration force,
and
for two endpoints on a discrete boundary, the equivalent surface force experienced by the cell is:
p0in order to be the initial equivalent surface force,
β1and beta2Respectively start and end angles, p, of the cell set1And p2The pressure caused by the residual force on the first and last cells can be determined by the following two equations:
and
β3and beta4The included angles between the end point of the first unit and the start point of the last unit and the concentrated force are respectively.
8. The method of claim 1, wherein the method is performed by
The boundary nodes for each grain are determined, wherein,
representing points on the boundaryDisplacement, H-1An inverse matrix representing H, a coefficient matrix formed in the boundary integral, G a coefficient matrix formed in the boundary integral,
representing the surface force of the node at local coordinates,
a physical coefficient matrix representing the nodes at local coordinates,
representing the physical strength of the node in local coordinates.
9. The method of claim 1, wherein the method is performed by
Determining the strain of a node in the particle domain from
Determining stress values, wherein G 'represents a coefficient matrix corresponding to displacement in the strain discrete matrix equation, H' represents a coefficient matrix corresponding to surface force in the strain discrete matrix equation,
a matrix of coefficients representing the corresponding physical forces in the strain discrete matrix equation,
m represents the number of interior points that need to be solved,
and
for each point in the x, y directionThe value of the strain is determined,
which represents the stress in a local coordinate system,
representing a stiffness matrix in a local coordinate system,
which represents the strain in the local coordinate system,
ij、kl、ik、jl、ilandjkrepresenting the kronecker notation, v is the poisson ratio.
10. A large-scale granular material internal stress and crushing simulation analysis device is characterized by comprising: a particle interaction solver and a particle internal stress solver;
the particle interaction solver is used for carrying out simulation calculation on the interaction of the large-scale particle aggregate;
the intra-particle stress solver comprises: the input module is used for receiving the material parameters, the gridding fraction, the grid type and the boundary condition information of the particles to be solved, which are input by an operator; a solving module for obtaining the internal stress of the particle body to be solved by adopting the large-scale particle material internal stress and crushing simulation analysis method of any one of claims 1 to 9.
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 true CN112084647A (en) | 2020-12-15 |
| CN112084647B 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) |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20050151733A1 (en) * | 2004-01-09 | 2005-07-14 | Microsoft Corporation | Multi-chart geometry images |
| CN102628861A (en) * | 2012-04-11 | 2012-08-08 | 武汉大学 | Method for simulating temperature cracking value of mass concrete |
| CN106980735A (en) * | 2017-04-06 | 2017-07-25 | 西南科技大学 | The method for numerical simulation of fragile material thermal fracture |
| CN109284523A (en) * | 2018-07-19 | 2019-01-29 | 同济大学 | A kind of rock soil medium Progressive failure, class solid-liquid phase change behavior analogy method |
| CN110598293A (en) * | 2019-09-03 | 2019-12-20 | 上海交通大学 | Method for predicting fracture damage behavior of micro-nano fiber reinforced composite material |
-
2020
- 2020-09-03 CN CN202010917113.3A patent/CN112084647B/en active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20050151733A1 (en) * | 2004-01-09 | 2005-07-14 | Microsoft Corporation | Multi-chart geometry images |
| CN102628861A (en) * | 2012-04-11 | 2012-08-08 | 武汉大学 | Method for simulating temperature cracking value of mass concrete |
| CN106980735A (en) * | 2017-04-06 | 2017-07-25 | 西南科技大学 | The method for numerical simulation of fragile material thermal fracture |
| CN109284523A (en) * | 2018-07-19 | 2019-01-29 | 同济大学 | A kind of rock soil medium Progressive failure, class solid-liquid phase change behavior analogy method |
| CN110598293A (en) * | 2019-09-03 | 2019-12-20 | 上海交通大学 | Method for predicting fracture damage behavior of micro-nano fiber reinforced composite material |
Non-Patent Citations (3)
| Title |
|---|
| WEIZHOU等: "Formulations of displacement discontinuity method for crack problems based on boundary element method", 《ENGINEERING ANALYSIS WITH BOUNDARY ELEMENTS》 * |
| 李洪超: "岩石RHT模型理论及主要参数确定方法研究", 《中国博士学位论文全文数据库 工程科技Ⅱ辑》 * |
| 邹宇雄等: "抗转动对颗粒材料组构特性的影响研究", 《岩土力学》 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN112084647B (en) | 2022-04-01 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| González et al. | Recent advances on the use of separated representations | |
| Peirce et al. | A spectral multipole method for efficient solution of large‐scale boundary element models in elastostatics | |
| Jin et al. | An edge‐based strain smoothing particle finite element method for large deformation problems in geotechnical engineering | |
| Benedetti et al. | A fast dual boundary element method for 3D anisotropic crack problems | |
| Yang et al. | A non-matching finite element-scaled boundary finite element coupled method for linear elastic crack propagation modelling | |
| Schneider et al. | A Large-Scale Comparison of Tetrahedral and Hexahedral Elements for Solving Elliptic PDEs with the Finite Element Method | |
| Huang et al. | A phase-field cohesive zone model integrated with cell-based smoothed finite element method for quasi-brittle fracture simulations of concrete at mesoscale | |
| CN104239277A (en) | Rapid discrete-element numerical-value calculating method applicable to GPU (Graphics Processing Unit) scalar-matrix operation | |
| Huang et al. | A CT image-driven computational framework for investigating complex 3D fracture in mesoscale concrete | |
| Yang et al. | Development of ABAQUS UEL/VUEL subroutines for scaled boundary finite element method for general static and dynamic stress analyses | |
| Zhang et al. | Numerical manifold method based on isogeometric analysis | |
| Goodarzi et al. | Modelling slope failure using a quasi-static MPM with a non-local strain softening approach | |
| Wu et al. | Arbitrary polygon mesh for elastic and elastoplastic analysis of solids using smoothed finite element method | |
| Bathe | The finite element method with “overlapping finite elements” | |
| CN112084647B (en) | Large-scale granular material internal stress and crushing simulation analysis method and device | |
| CN106407620A (en) | ABAQUS-based engineering structure response surface stochastic finite element analysis processing method | |
| CN104992046A (en) | Computing system and method of fluid mechanics | |
| Ciancio et al. | A method for the calculation of inter-element stresses in 3D | |
| Daneshmand et al. | Static and dynamic analysis of 2D and 3D elastic solids using the modified FGFEM | |
| Xia et al. | A parallel, implicit reconstructed discontinuous Galerkin method for the compressible flows on 3D arbitrary grids | |
| Liu et al. | Circumventing volumetric locking in stabilized smoothed particle finite element method and its application to dynamic large deformation problems | |
| Liu et al. | Smoothed finite elements large deformation analysis | |
| Rossikhin et al. | Vibrations of Suspension Bridges: Fractional Derivative Model | |
| Okada et al. | Application of s-version finite element method to two-dimensional fracture mechanics problems | |
| Ptaszny et al. | Stress analysis of linear elastic structures by the fast multipole boundary element method |
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 |