CN111611695A - Automatic calibration method for discrete element linear stiffness parameter in simulation of rock and soil material - Google Patents

Automatic calibration method for discrete element linear stiffness parameter in simulation of rock and soil material Download PDF

Info

Publication number
CN111611695A
CN111611695A CN202010393302.5A CN202010393302A CN111611695A CN 111611695 A CN111611695 A CN 111611695A CN 202010393302 A CN202010393302 A CN 202010393302A CN 111611695 A CN111611695 A CN 111611695A
Authority
CN
China
Prior art keywords
sample
parameter
parameters
simulation
particle
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.)
Pending
Application number
CN202010393302.5A
Other languages
Chinese (zh)
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.)
Taiyuan University of Technology
Original Assignee
Taiyuan University of Technology
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 Taiyuan University of Technology filed Critical Taiyuan University of Technology
Priority to CN202010393302.5A priority Critical patent/CN111611695A/en
Publication of CN111611695A publication Critical patent/CN111611695A/en
Pending legal-status Critical Current

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
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • 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)
  • Theoretical Computer Science (AREA)
  • Geometry (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Graphics (AREA)
  • Software Systems (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Investigating Strength Of Materials By Application Of Mechanical Stress (AREA)

Abstract

The invention discloses an automatic calibration method for discrete element linear stiffness parameters in simulation of rock and soil materials. The method utilizes an analytic formula as an initial estimation value of linear elastic contact parameters, and realizes automatic training of the contact parameters in the discrete element model by taking the particle stiffness parameters as independent variables and numerical simulation actual results as dependent variables. By the method, the target parameters can be obtained through automatic training, and the time required by the traditional manual checking method is greatly reduced; the macroscopic deformation analytical solution based on the isotropic discrete sample is used as the initial value of parameter training, so that the problems of local optimal trap, overlong training time and the like in the parameter training process are effectively avoided; the elastic parameters are calibrated through the response of the sample under the condition of small strain, and the parameter calibration efficiency is high. The method is simple in implementation process, and is suitable for rapidly calibrating mesoscopic deformation parameters of various discrete materials such as sand, rockfill and the like and bonded discrete materials such as deposited rocks and ceramics during simulation research by adopting discrete elements.

Description

Automatic calibration method for discrete element linear stiffness parameter in simulation of rock and soil material
Technical Field
The invention belongs to the field of geotechnical engineering, and particularly relates to an automatic calibration method for discrete element linear stiffness parameters in simulation of geotechnical materials.
Background
As the earth materials such as sandy soil, rockfill and the like have natural discrete characteristics, the discrete element method taking particles as basic units is widely applied to the analysis and research of the problems of the earth materials and the geotechnical engineering. However, since the mesoscale parameters adopted by the discrete element algorithm are not easily measured directly by physical tests, most of the current discrete element simulations for geotechnical engineering firstly simulate some conventional geotechnical tests (such as triaxial tests) and continuously debug the discrete element mesoscale parameters manually, until the simulation object is basically consistent with the physical tests on the target macroscopic indexes, the used simulation parameters are not taken as a set of reliable input parameters. The traditional manual calibration method is time-consuming and labor-consuming, and a higher precision is difficult to achieve.
The linear model and the hertzian model are two basic contact models that are widely used in discrete element simulations. The input parameters used by each model are different. The linear contact model adopts normal contact stiffness and tangential contact stiffness as parameters for controlling deformation; the hertzian model uses the particle shear modulus and poisson's ratio as parameters to control deformation. Due to the difference of microscopic deformation mechanisms, the two methods have difference in calibration strategies. The method mainly solves the problem of parameter calibration of microscopic deformation parameters (namely, normal contact stiffness and tangential contact stiffness) in a linear model.
In recent years, the problem of calibrating discrete element parameters has attracted much attention.
The Chinese invention patent (application number: 201810632432.2, patent name: a method for rapidly determining discrete element model parameters of rock brittle materials) provides a method for determining discrete element model parameters of rock brittle materials. The method needs to simulate enough number of samples to perform parameter sensitivity analysis to estimate the value range of the parameters and provide a basis for optimizing the parameter value by adopting a statistical method in the next step. The method provides a set of available reference processes for calibrating the discrete element parameters, and is still a statistical fitting method based on sufficient sample empirical data. The method has the following disadvantages: (1) a sufficient number of numerical samples need to be simulated to provide basis for sensitivity analysis and statistical fitting, and the calculation amount is large; (2) after numerical value output, main microscopic influence parameters behind a single macroscopic index are determined from sensitivity analysis, and then parameter values are optimized through Latin hypercube sampling and ant colony algorithm, so that physical relation of discrete materials behind the macroscopic and microscopic parameters is ignored by depending on a statistical method too, and the implementation process is relatively complex.
The Chinese invention patent (application number: 201811477615.8, patent name: microscopic parameter calibration method in discrete element simulation) provides a microscopic parameter calibration method in discrete element simulation. The method is also a traditional manual trial and error checking method in the aspect of calibrating the particle rigidity parameters, and cannot improve the calibration efficiency substantially.
The Chinese invention patent (application number: 201910080559.2, patent name: an automatic discrete meta-material training method) provides an automatic discrete meta-material parameter training method. The method takes the analytic relation of macro-micro parameters of granules based on regular arrangement as a material training initial value, obtains the next updated micro parameters by multiplying the current micro parameters by the ratio of the actual macro parameters to the target macro parameters, and continuously trains until the target value is reached. The method realizes automatic training and makes obvious progress compared with the traditional method. However, the following disadvantages also exist: (1) all the microscopic parameters are trained simultaneously, and all the parameters interfere with each other in the process, so that the calibration efficiency is influenced; (2) the updating strategy is derived from an empirical assumption, lacks a strict mathematical basis and cannot ensure sufficient convergence efficiency and convergence accuracy; (3) the input mesoscopic initial value is derived from the deformation relation of the particles with single particle size and regularly arranged, and the difference with the actual situation that the texture and the particle size of the particles in the actual geotechnical material are randomly distributed is larger. The initial estimation of the macroscopic response of the actual dispersion is too deviated, so that the training time is too long, and the risk of the material facing a local optimal trap in the training process is obviously increased; (4) the parameter update process ignores the macro-meso parameter relationship of the dispersoid, for example: the normal stiffness and the tangential stiffness simultaneously influence the elastic modulus and the Poisson ratio in macroscopic parameters, and according to a parameter updating strategy used by the method, the normal stiffness or the tangential stiffness is updated only by the elastic modulus or the Poisson ratio, which is actually not good; (5) the differences between the two-dimensional model and the three-dimensional model in the parameter training strategy and process are not distinguished. For example: the same initial value estimation is used for the two-dimensional model and the three-dimensional model, and the model cannot achieve the optimal convergence efficiency.
Disclosure of Invention
The invention aims to provide a set of automatic linear contact rigidity calibration method which is suitable for two-dimensional and three-dimensional simulation and has high efficiency and stable performance when discrete elements simulate sand, rockfill and other discrete bodies.
The object of the invention is achieved by the following steps: the invention provides an automatic calibration method of discrete element linear stiffness parameters in simulation of rock and soil materials, which comprises the following steps:
determining the Young modulus and Poisson's ratio of a rock-soil material to be simulated in an indoor true triaxial or biaxial test mode;
secondly, calculating initial estimated values of the normal rigidity and the tangential rigidity of the particle contact according to a uniform dispersion macro-microscopic deformation parameter analytical formula deduced under a small strain condition;
thirdly, establishing a discrete element biaxial or triaxial numerical test by using the calculated stiffness parameter initial estimation value;
setting the friction coefficient among the sample particles to be a value of 1.0 or more before loading, implementing numerical simulation of a small-strain triaxial test, and stopping loading until the axial strain reaches a preset threshold value to obtain the elastic modulus and the Poisson ratio of the current simulation sample;
fifthly, constructing an error function L;
sixthly, judging whether the size of the error function L meets an error allowable value or not; if the rigidity parameter meets the requirement, the rigidity parameter used by the current model is a calibration parameter, and the calibration is finished; if not, continuing to execute calibration;
updating the rigidity parameters based on an improved gradient descent method;
recalculating the numerical model by using the updated stiffness parameters, updating the size of the error function and evaluating whether the size of the error function meets an error allowable value or not; if not, the process of the steps seven and eight is continuously executed until the error function is lower than the error allowable value.
Wherein the initial estimate of particle contact normal
Figure BDA0002486424920000021
And an initial estimate of tangential stiffness
Figure BDA0002486424920000022
The calculation formula of (2) is as follows:
two-dimensional discrete element model:
Figure BDA0002486424920000023
three-dimensional discrete element model:
Figure BDA0002486424920000024
wherein E is the Young modulus with small strain obtained by the test, V is the Poisson' S ratio obtained by the test, S is the area of the two-dimensional sample, V is the volume of the three-dimensional sample, and N iscIs the number of contacts in the sample; when the sample is in single particle size distribution, r is the particle radius; when the sample has a multi-particle size distribution, r is the median particle size of the sample.
The method for calculating the ratio of the area/volume to the contact number of the discrete rock-soil material sample to be simulated comprises the following steps:
Figure BDA0002486424920000031
wherein 2D and 3D represent 2-dimensional and 3-dimensional models, respectively; phi is porosity;
Figure BDA0002486424920000032
is the average coordination number of the particles.
Wherein the initial friction coefficient f at the time of particle generation is usedcThe porosity and average coordination number empirical values of the final sample are estimated as follows:
Figure BDA0002486424920000033
。
wherein the median particle size used to calculate the initial stiffness of the model for the multi-particle dispersion samples.
After the true triaxial or biaxial simulation is finished, the stress-strain curve of the sample in the loading process should be checked to ensure that the sample only generates recoverable elastic deformation during the loading process.
Wherein the error function L is:
Figure BDA0002486424920000034
wherein, β1And β2Is two value ranges of [0,1 ]]And satisfy β1+β2＝1，E*And ν is the calculated elastic modulus and poisson's ratio, respectively, of the numerical samples.
Wherein, if the particle contact stiffness parameter is updated for the first time, updating according to the following mode:
Figure BDA0002486424920000035
wherein, the symbol is assigned αnAnd αsIs coefficient, and the value is-0.5-0.5.
If the current particle contact stiffness parameter is updated for the second time or more, calculating an error function L to k according to a difference methodnAnd ksPartial derivatives of
Figure BDA0002486424920000036
And
Figure BDA0002486424920000037
if order is given by knAnd ksThe error function for the argument is f, let:
f(kn,ks)＝L(E*,ν*) (7)
the two partial derivatives are estimated differentially as follows:
Figure BDA0002486424920000038
。
if the current particle contact stiffness parameter is updated for the second time or more, updating according to an improved gradient descent method, wherein the updating formula is as follows:
Figure BDA0002486424920000041
。
compared with the prior art, the invention has the technical effects that: (1) the microscopic linear stiffness parameters in the particle units in the discrete element simulation process can be automatically and quickly calibrated, and the analysis and research of discrete elements on geotechnical engineering problems are greatly promoted; (2) the calibration process is iterated based on the improved gradient descent method, the calibration convergence is good, and the efficiency is high; (3) the deformation parameters are calibrated based on a small-strain triaxial or biaxial loading test, the simulation time is short, and the calibration result is reliable.
In conclusion, the method is efficient, stable, clear in implementation process and convenient to program, can provide technical support for calibrating the linear stiffness parameters of the discrete elements, and promotes the analysis and research of the discrete element technology on geotechnical engineering problems; the method is not only suitable for single-particle-size particle dispersions and multi-particle-size particle dispersions, but also has better calibration capability for two-dimensional and three-dimensional samples. Within several iteration steps, the prediction precision of the parameters can reach within 1 percent.
Drawings
FIG. 1 is a schematic diagram of an automatic calibration technique for discrete element linear contact stiffness parameters;
FIG. 2 is a schematic diagram of a model of the two-dimensional single-particle-size random particle sample;
FIG. 3 is a schematic diagram of a model of the two-dimensional multi-particle-size random particle sample;
FIG. 4 is a schematic diagram of a three-dimensional single-particle random particle sample.
Detailed Description
For a more clear understanding of the technical features, objects and effects of the present invention, embodiments of the present invention will now be described in detail with reference to the accompanying drawings.
As shown in FIG. 1, the invention provides an automatic calibration method of discrete element linear stiffness parameters in simulation of geotechnical materials, which comprises the following steps:
step (1): and determining the macroscopic parameters of the rock and soil sample to be simulated by means of indoor true triaxial or biaxial tests and the like. The macroscopic parameters used were young's modulus E at small strain and poisson ratio υ.
Step (2): inputting the target macroscopic parameters obtained by the test into a uniform dispersion macroscopic deformation parameter analytical formula (the formula is shown below, and the derivation process is shown later) derived under the condition of small strain, and respectively determining initial estimated values of the normal stiffness and the tangential stiffness of the particles in contact under the conditions of two-dimensional simulation and three-dimensional simulation
Figure BDA0002486424920000044
And
Figure BDA0002486424920000045
two-dimensional discrete element model:
Figure BDA0002486424920000042
three-dimensional discrete element model:
Figure BDA0002486424920000043
wherein S is the area of the two-dimensional sample, V is the volume of the three-dimensional sample, and NcThe total number of contacts participating in load transfer in the sample; when the single particle size distribution of the sample is present,r is the particle radius; when the sample has a multi-particle size distribution, r is the median particle size of the sample. The ratio of area/volume to number of contacts for the sample was:
Figure BDA0002486424920000051
wherein 2D and 3D represent 2-dimensional and 3-dimensional models, respectively; phi is porosity;
Figure BDA0002486424920000052
is the average coordination number of the particles. For discrete dispersoid samples, the initial coefficient of friction at which the particles are produced tends to have a large effect on the final degree of compaction of the particle body. When the initial coefficient of friction is fcThe values of porosity and average coordination number of the final sample can then be estimated according to the following empirical formula:
Figure BDA0002486424920000053
and (3): and (3) establishing a discrete element biaxial or true triaxial numerical test, and inputting the particle contact stiffness parameter calculated in the step (2) as an initial parameter into a model for numerical calculation, wherein the numerical test adopts the confining pressure same as that of an indoor test.
And (4): before loading, the friction coefficient among sample particles is set to be 1.0, and the numerical simulation of the small-strain triaxial test is implemented until the axial strain reaches a small value, and the loading is stopped. (e.g., axial strain 0.05%). After the simulation is complete, the stress-strain curve of the specimen during loading should be examined to ensure that only recoverable elastic deformation of the specimen occurs during this loading. The inter-particle friction coefficient was set to 1.0 in order that no sliding deformation (plastic deformation) occurred between particles when the sample was deformed in a small axial direction as a whole. Calculating the modulus of elasticity E obtained by numerical sample simulation*And poisson's ratio ν. The calculation method comprises the following steps:
Figure BDA0002486424920000054
Figure BDA0002486424920000055
and (5): constructing an error function L:
Figure BDA0002486424920000056
wherein, β1And β2Is two value ranges of [0,1 ]]And satisfy β1+β2When modulus of elasticity and poisson's ratio need to be calibrated simultaneously, β1＝β20.5, when calibrated to macroscopic parameters targeting only the modulus of elasticity, β1＝1，β1＝0.
And (6): and judging whether the size of the error function L meets an error allowable value or not. If the rigidity parameter meets the requirement, the rigidity parameter used by the current model is a calibration parameter, and the calibration is finished; if not, continuing to execute the calibration step. The specific process is as follows: and (4) substituting the elastic modulus E and the Poisson ratio v obtained by the simulation in the step (4) into the error function constructed in the step (5), and judging whether the error function is smaller than an error allowable value (for example: 0.01%). If so, the used microscopic parameters are reasonable estimation parameters for simulating the test; if not, updating the parameters according to the description of the step (7).
And (7): and updating the rigidity parameter. If the particle contact stiffness parameter is updated for the first time, updating according to the following mode:
Figure BDA0002486424920000061
wherein, the symbol is assigned αnAnd αsIs a coefficient, generally takes a value between-0.5 and 0.5.
If the current particle contact stiffness parameter is updated for the second time or more, updating according to an improved gradient descent method, wherein the updating formula is as follows:
Figure BDA0002486424920000062
where η is the learning rate,
Figure BDA0002486424920000063
and
Figure BDA0002486424920000064
are error functions L to k respectivelynAnd ksPartial derivatives of (a). If order is given by knAnd ksThe error function for the argument is f, let:
f(kn,ks)＝L(E*,ν*) (7)
the two partial derivatives can be estimated differentially:
Figure BDA0002486424920000065
and (8): and recalculating the numerical model by adopting the updated rigidity parameters, updating the size of the error function and evaluating whether the size of the error function meets an error allowable value. If not, continuing to execute the processes in the steps (7) and (8) until the error function is lower than the error allowable value. The specific process is as follows: and (4) carrying out discrete element biaxial or true triaxial compression test calculation on the stiffness parameters obtained by updating in the step (7), calculating the size of an error function according to the newly calculated elastic modulus E and Poisson ratio v, and judging whether the error function is smaller than an error allowable value or not. If the stiffness parameter is smaller than the calibration parameter, ending the calculation, wherein the stiffness parameter of the current model is the calibration parameter; and (5) if the error requirement is not met, updating the parameter values according to the step (7) and recalculating until the error meets the requirement.
Attached: uniform dispersion macro-micro deformation parameter analytic formula derivation process under small strain condition
(1) Basic assumptions and inferences
The basic assumption is that: the homogeneous bulk sample deforms uniformly (uniform strain) throughout the sample space, and the bulk sample deforms to a small deformation (recoverable elastic deformation).
Deducing the idea: using continuous medium mechanics theory, the dispersion samples were homogenized and equivalent to the corresponding continuum. The particle system is consistent in strain energy with the corresponding equivalent continuum.
Note: the following derivation process is done based on tensor operations. In tensor calculation, the x axis, the y axis and the z axis in a cartesian coordinate system are generally represented by subscripts 1,2 and 3; the component notation of tensor x is xi(first order) or xij(second order), where i and j may be any of 1,2, 3.
(2) And (3) strain energy calculation:
the strain energy U stored in the particle system is:
Figure BDA0002486424920000066
wherein: n is a radical ofcThe total number of contacts participating in load transfer in the overall particulate system;
Figure BDA0002486424920000067
and
Figure BDA0002486424920000068
respectively, the normal and tangential relative displacements of the kth contact in the particle system, k being 1, … Nc；
Figure BDA0002486424920000069
And
Figure BDA00024864249200000610
normal and tangential contact forces, respectively, imparted by contact k;
Figure BDA0002486424920000071
and
Figure BDA0002486424920000072
refers to the slight deformation of the contact k in the normal and tangential directions, respectively.
For the linear contact model, before the contact does not slip, the relationship between the contact force and the relative displacement between the particles is:
Figure BDA0002486424920000073
wherein: k is a radical ofnAnd ksIs the normal and tangential stiffness of the particle. Strain energy U stored in kth linear contactkComprises the following steps:
Figure BDA0002486424920000074
(3) particle displacement versus strain
Assuming that the particles A and B are any two particles in contact with each other in the sample, and the displacement of the particles is made equal to the displacement of the equivalent continuum at the corresponding point, the relative displacement between the particles
Figure BDA0002486424920000075
Corresponding strain on equivalent continuum
Figure BDA0002486424920000076
The relationship of (1) is:
Figure BDA0002486424920000077
wherein:
Figure BDA0002486424920000078
and
Figure BDA0002486424920000079
particles A and B, respectively, along xjCoordinate of direction, LkIs the distance between the particles a and B,
Figure BDA00024864249200000710
is the equivalent strain between particles a and B.
Figure BDA00024864249200000711
The jth component of the kth contact.
Normal relative displacement between particles
Figure BDA00024864249200000712
Can be expressed as:
Figure BDA00024864249200000713
tangential relative displacement between particles
Figure BDA00024864249200000714
Can be expressed as:
Figure BDA00024864249200000715
by substituting equations (14), (16) and (17) into equation (12), the strain energy of the particle system can be expressed as:
Figure BDA00024864249200000716
(4) equivalent elastic stiffness matrix for particle systems
The strain energy density of the entire particle system can be obtained using the strain energy divided by the area (two-dimensional) or volume (three-dimensional) of the particle system:
Figure BDA00024864249200000717
wherein S is the area of the two-dimensional sample and V is the volume of the three-dimensional sample. To simplify the derivation process, the following derivation will be done based on a three-dimensional particle system, and the results of the two-dimensional model can be obtained by replacing the volume V with the area S and ignoring the 3 rd index of each physical tensor.
Based on the assumption that the deformation is uniform within the dispersion (strain is equal everywhere), the local strain at any location in the particle system is equal to the overall strain of the system:
Figure BDA00024864249200000718
according to the theory of elastic mechanics, the stress tensor of a continuum can be obtained by partial differentiation of the strain tensor with respect to strain energy density:
Figure BDA0002486424920000081
elastic stiffness matrix C of particle systemijmnCan be further improved by using the stress tensor σijTo strain tensormnTaking partial differential to obtain:
Figure BDA0002486424920000082
wherein,inis the kronecker function.
(5) Analytic relationship of macro-micro deformation parameters in uniform particle system
In uniform particle systems with uniform particle size and isotropic internal texture, the elastic stiffness matrix C of the particle systemijmnCan be further simplified into:
Figure BDA0002486424920000083
in elastomechanics, the generally isotropic elastic continuum stiffness matrix CijmnCommonly expressed as:
Figure BDA0002486424920000084
wherein E and v are the elastic modulus and Poisson's ratio of a generally isotropic elastic continuum.
By comparing the elastic stiffness matrix of the particle system (equation (23)) with the stiffness matrix of a generally isotropic elastic continuum (equation (24)), the relationship between macro and micro deformation parameters of the particle system at small strains can be determined.
Figure BDA0002486424920000085
Figure BDA0002486424920000086
Equation (25) and equation (26) are combined, and the estimated equation for the particle stiffness parameter can be derived as follows:
Figure BDA0002486424920000087
Figure BDA0002486424920000088
this formula is derived based on isotropic homogeneous particle systems of the same radius, for particle systems with arbitrary particle sizes, the internal texture inevitably has a certain degree of anisotropic character. Therefore, the present analytical formula is used in the present invention only for the estimation of the true particle stiffness parameter. The accurate result will be obtained by the iterative solution auto-training proposed by this patent.
Detailed description of the invention
The method is applied to the calibration of a two-dimensional single-particle-diameter random particle sample (as shown in figure 2), and the test parameters are as follows:
Figure BDA0002486424920000091
TABLE 1 two-dimensional single-particle-size randomly arranged particle model parameters
The simulation steps are as follows:
step (1): and determining the macroscopic parameters of the geotechnical material to be simulated in an indoor test mode. The macroscopic parameters of the geotechnical material in the embodiment are as follows: young modulus E was 10GPa and Poisson ratio upsilon was 0.2.
Step (2): inputting the target macroscopic parameters obtained by the test into a formula, and respectively determining initial estimated values of the normal stiffness and the tangential stiffness of the particle contact in the two-dimensional model
Figure BDA0002486424920000092
And
Figure BDA0002486424920000093
Figure BDA0002486424920000094
wherein S is the area of the two-dimensional sample, NcThe number of contacts in the sample, r the median diameter of the granules, and phi the porosity;
Figure BDA0002486424920000096
is the average coordination number of the particles.
And (3): and (3) establishing a discrete element biaxial numerical test, putting the particle contact rigidity parameter calculated in the step (2) into a model as an initial parameter for numerical calculation, wherein the same confining pressure as that of an indoor test is adopted in the numerical test and is 1 MPa.
And (4): before loading, the friction coefficient among sample particles is set to be 1.0, and the numerical simulation of the small-strain triaxial test is implemented until the axial strain reaches a smaller value of 0.05 percent, and the loading is stopped. The stress-strain curve during the loading of the specimen was examined to confirm that only recoverable elastic deformation of the specimen occurred during this loading. Calculating the modulus of elasticity E obtained by numerical sample simulation*And poisson's ratio ν.
And (5): constructing an error function L:
Figure BDA0002486424920000095
wherein, β1＝β2＝0.5。
And (6): and judging whether the size of the error function L meets an error allowable value or not. The first set of error functions is found to be 0.0671 in size and does not meet the tolerance requirements.
And (7): and updating the rigidity parameter. If the particle contact stiffness parameter is updated for the first time, updating according to the following mode:
Figure BDA0002486424920000101
wherein, the symbol is assigned αnAnd αsThe coefficients are all 0.22.
If the current particle contact stiffness parameter is updated for the second time or more, updating according to an improved gradient descent method, wherein the updating formula is as follows:
Figure BDA0002486424920000102
where η is the learning rate,
Figure BDA0002486424920000103
and
Figure BDA0002486424920000104
are error functions L to k respectivelynAnd ksPartial derivatives of (a). If order is given by knAnd ksThe error function for the argument is f, let:
f(kn,ks)＝L(E*,ν*) (7)
the two partial derivatives can be estimated differentially:
Figure BDA0002486424920000105
and (8): and recalculating the numerical model by adopting the updated rigidity parameters, updating the size of the error function and evaluating whether the size of the error function meets an error allowable value. If not, continuing to execute the processes in the steps (7) and (8) until the error function is lower than the error allowable value. The specific process is as follows: and (4) carrying out discrete element biaxial or true triaxial compression test calculation on the stiffness parameters obtained by updating in the step (7), calculating the size of an error function according to the newly calculated elastic modulus E and Poisson ratio v, and judging whether the error function is smaller than an error allowable value or not. If the stiffness parameter is smaller than the calibration parameter, ending the calculation, wherein the stiffness parameter of the current model is the calibration parameter; and (5) if the error requirement is not met, updating the parameter values according to the step (7) and recalculating. The model iterates to the 9 th timeAnd when the error function is 2.80E-05, the error meets the requirement, the calibration process is finished, and the obtained optimal stiffness parameter estimation value is as follows: k is a radical ofn＝9605966N/m,ks3687070N/m. The parameter update procedure in the record calibration is shown in table 2:
Figure BDA0002486424920000106
TABLE 2 two-dimensional single-grain-diameter random-grain sample calibration process
Detailed description of the invention
The method of the invention is applied to the calibration of a two-dimensional multi-particle-size random particle sample (as shown in FIG. 3), and the test parameters are shown in Table 3:
Figure BDA0002486424920000111
TABLE 3 two-dimensional multiple-particle-size randomly arranged particle model parameters
The simulation steps are as follows:
step (1): and determining the macroscopic parameters of the geotechnical material to be simulated in an indoor test mode. The macroscopic parameters of the geotechnical material in the embodiment are as follows: young modulus E was 10GPa and Poisson ratio upsilon was 0.2.
Step (2): inputting the target macroscopic parameters obtained by the test into a formula, and respectively determining initial estimated values of the normal stiffness and the tangential stiffness of the particle contact in the two-dimensional model
Figure BDA0002486424920000112
And
Figure BDA0002486424920000113
Figure BDA0002486424920000114
wherein S is the area of the two-dimensional sample, NcThe number of contacts in the sample, r the median diameter of the granules, and phi the porosity;
Figure BDA0002486424920000115
is the average coordination number of the particles.
And (3): and (3) establishing a discrete element biaxial numerical test (as shown in figure 3), putting the particle contact rigidity parameter calculated in the step (2) into a model as an initial parameter for numerical calculation, and adopting the same confining pressure as that of an indoor test in the numerical test, wherein the confining pressure is 1 MPa.
And (4): before loading, the friction coefficient among sample particles is set to be 1.0, and the numerical simulation of the small-strain triaxial test is implemented until the axial strain reaches a smaller value of 0.05 percent, and the loading is stopped. The stress-strain curve during the loading of the specimen was examined to confirm that only recoverable elastic deformation of the specimen occurred during this loading. Calculating the modulus of elasticity E obtained by numerical sample simulation*And poisson's ratio ν.
And (5): constructing an error function L:
Figure BDA0002486424920000116
wherein, β1＝β2＝0.5。
And (6): and judging whether the size of the error function L meets an error allowable value or not. The first set of error functions is obtained with a magnitude of 0.0253, which does not meet the tolerance requirements.
And (7): and updating the rigidity parameter. If the particle contact stiffness parameter is updated for the first time, updating according to the following mode:
Figure BDA0002486424920000117
wherein, the symbol is assigned αnAnd αsThe coefficients are all 0.2.
If the current particle contact stiffness parameter is updated for the second time or more, updating according to an improved gradient descent method, wherein the updating formula is as follows:
Figure BDA0002486424920000121
where η is the learning rate,
Figure BDA0002486424920000122
and
Figure BDA0002486424920000123
are error functions L to k respectivelynAnd ksPartial derivatives of (a). If order is given by knAnd ksThe error function for the argument is f, let:
f(kn,ks)＝L(E*,ν*) (7)
the two partial derivatives can be estimated differentially:
Figure BDA0002486424920000124
and (8): and recalculating the numerical model by adopting the updated rigidity parameters, updating the size of the error function and evaluating whether the size of the error function meets an error allowable value. If not, the process of the steps seven and eight is continuously executed until the error function is lower than the error allowable value. The specific process is as follows: and (4) carrying out discrete element biaxial or triaxial compression test calculation on the stiffness parameters obtained by updating in the step (7), calculating the size of an error function according to the newly calculated elastic modulus E and Poisson ratio v, and judging whether the error function is smaller than an error allowable value. If the stiffness parameter is smaller than the calibration parameter, ending the calculation, wherein the stiffness parameter of the current model is the calibration parameter; and (5) if the error requirement is not met, updating the parameter values according to the step (7) and recalculating. When the model is iterated to the 7 th time, the size of an error function is 6.65E-05, the error meets the requirement, the calibration process is finished, and the obtained optimal stiffness parameter estimation value is as follows: kn is 10533175N/m, ks is 3546570N/m. The procedure for recording parameter updates in the calibration is shown in table 4:
Figure BDA0002486424920000125
TABLE 4 two-dimensional multiple-grain-size random grain sample calibration procedure
Detailed description of the invention
The method is applied to the calibration of a three-dimensional single-particle-diameter random particle sample (as shown in FIG. 4), and the test parameters are as shown in Table 5:
Figure BDA0002486424920000131
TABLE 5 three-dimensional multiple-particle-size randomly arranged particle model parameters
The simulation steps are as follows:
step (1): and determining the macroscopic parameters of the geotechnical material to be simulated in an indoor test mode. The macroscopic parameters of the geotechnical material in the embodiment are as follows: young modulus E was 10GPa and Poisson ratio upsilon was 0.2.
Step (2): inputting the target macroscopic parameters obtained by the test into a formula, and respectively determining initial estimated values of the normal stiffness and the tangential stiffness of the particle contact in the three-dimensional model
Figure BDA0002486424920000132
And
Figure BDA0002486424920000133
Figure BDA0002486424920000134
wherein V is the volume of the three-dimensional sample, NcThe number of contacts in the sample, r the median diameter of the granules, and phi the porosity;
Figure BDA0002486424920000135
is the average coordination number of the particles.
And (3): and (3) establishing a discrete element biaxial numerical test (as shown in figure 3), putting the particle contact rigidity parameter calculated in the step (2) into a model as an initial parameter for numerical calculation, and adopting the same confining pressure as that of an indoor test in the numerical test, wherein the confining pressure is 1 MPa.
And (4): before loading, the coefficient of friction between the sample particles was set to 1.0, and small strain triaxial was performedThe experiment numerical simulation is carried out, and the loading is stopped until the axial strain reaches a smaller value of 0.05 percent. The stress-strain curve during the loading of the specimen was examined to confirm that only recoverable elastic deformation of the specimen occurred during this loading. Calculating the modulus of elasticity E obtained by numerical sample simulation*And poisson's ratio ν.
And (5): constructing an error function L:
Figure BDA0002486424920000136
if the elastic modulus is only used as the check target, β1＝1，β2＝0。
And (6): and judging whether the size of the error function L meets an error allowable value or not. The first set of error functions is obtained with a magnitude of 0.0253, which does not meet the tolerance requirements.
And (7): and updating the rigidity parameter. If the particle contact stiffness parameter is updated for the first time, updating according to the following mode:
Figure BDA0002486424920000137
wherein, the symbol is assigned αnAnd αsThe values are respectively 0.1 and-0.5.
If the current particle contact stiffness parameter is updated for the second time or more, updating according to an improved gradient descent method, wherein the updating formula is as follows:
Figure BDA0002486424920000138
where η is the learning rate,
Figure BDA0002486424920000141
and
Figure BDA0002486424920000142
are error functions L to k respectivelynAnd ksPartial derivatives of (a). If order is given by knAnd ksThe error function for the argument is f, let:
f(kn,ks)＝L(E*,ν*) (7)
the two partial derivatives can be estimated differentially:
Figure BDA0002486424920000143
and (8): and recalculating the numerical model by adopting the updated rigidity parameters, updating the size of the error function and evaluating whether the size of the error function meets an error allowable value. If not, the process of the steps seven and eight is continuously executed until the error function is lower than the error allowable value. The specific process is as follows: and (4) carrying out discrete element biaxial or triaxial compression test calculation on the stiffness parameters obtained by updating in the step (7), calculating the size of an error function according to the newly calculated elastic modulus E and Poisson ratio v, and judging whether the error function is smaller than an error allowable value. If the stiffness parameter is smaller than the calibration parameter, ending the calculation, wherein the stiffness parameter of the current model is the calibration parameter; and (5) if the error requirement is not met, updating the parameter values according to the step (7) and recalculating. When the model is iterated to the 7 th time, the size of an error function is 6.65E-05, the error meets the requirement, the calibration process is finished, and the obtained optimal stiffness parameter estimation value is as follows: kn is 6087605N/m, ks is 2056510N/m. The procedure for recording parameter updates in the calibration is shown in table 6:
Figure BDA0002486424920000144
TABLE 6 calibration process of three-dimensional single-particle-diameter random particle sample
In conclusion, the method is efficient, stable, clear in implementation process and convenient to program, can provide technical support for calibrating the linear stiffness parameters of the discrete elements, and promotes the analysis and research of the discrete element technology on geotechnical engineering problems; the method is not only suitable for single-particle-size particle dispersions and multi-particle-size particle dispersions, but also has better calibration capability for two-dimensional and three-dimensional samples. Within several iteration steps, the prediction precision of the parameters can reach within 1 percent.
Compared with the prior art, the invention has the technical effects that: (1) the microscopic linear stiffness parameters in the particle units in the discrete element simulation process can be automatically and quickly calibrated, and the analysis and research of discrete elements on geotechnical engineering problems are greatly promoted; (2) the calibration process is iterated based on the improved gradient descent method, the calibration convergence is good, and the efficiency is high; (3) the deformation parameters are calibrated based on a small-strain triaxial or biaxial loading test, the simulation time is short, and the calibration result is reliable.
While the present invention has been described with reference to the embodiments shown in the drawings, the present invention is not limited to the embodiments, which are illustrative and not restrictive, and it will be apparent to those skilled in the art that various changes and modifications can be made therein without departing from the spirit and scope of the invention as defined in the appended claims.

Claims (10)

1. An automatic calibration method for discrete element linear stiffness parameters in simulation of geotechnical materials is characterized by comprising the following steps:
determining the Young modulus and Poisson's ratio of a rock-soil material sample to be simulated in an indoor true triaxial or biaxial test mode;
secondly, calculating initial estimated values of the normal rigidity and the tangential rigidity of the particle contact according to a uniform dispersion macro-microscopic deformation parameter analytical formula deduced under a small strain condition;
thirdly, establishing a discrete element biaxial or triaxial numerical test of the rock-soil material sample to be simulated by using the calculated initial estimated value of the stiffness parameter;
setting the friction coefficient among the test sample particles to be a value of 1.0 or more before loading, implementing numerical simulation of a small-strain triaxial test, and stopping loading until the axial strain reaches a preset threshold value to obtain the elastic modulus and the Poisson ratio of the current simulation sample;
fifthly, constructing an error function L;
sixthly, judging whether the size of the error function L meets an error allowable value or not; if the rigidity parameter meets the requirement, the rigidity parameter used by the current model is a calibration parameter, and the calibration is finished; if not, continuing to execute calibration;
updating the rigidity parameters based on an improved gradient descent method;
recalculating the numerical model by using the updated stiffness parameters, updating the size of the error function and evaluating whether the size of the error function meets an error allowable value or not; if not, the process of the steps seven and eight is continuously executed until the error function is lower than the error allowable value.
2. The automatic calibration method for the linear stiffness parameters of the discrete elements during simulation of the geotechnical materials according to claim 1, wherein the initial particle contact normal estimation value
Figure FDA0002486424910000014
And an initial estimate of tangential stiffness
Figure FDA0002486424910000015
The calculation formula of (2) is as follows:
two-dimensional discrete element model:
Figure FDA0002486424910000011
three-dimensional discrete element model:
Figure FDA0002486424910000012
wherein E is the Young modulus with small strain obtained by the test, V is the Poisson' S ratio obtained by the test, S is the area of the two-dimensional sample, V is the volume of the three-dimensional sample, and N iscIs the number of contacts in the sample; when the sample is in single particle size distribution, r is the particle radius; when the sample has a multi-particle size distribution, r is the median particle size of the sample.
3. The automatic calibration method for the linear stiffness parameters of the discrete elements during the simulation of the geotechnical materials according to claim 1, wherein for discrete samples of the geotechnical materials to be simulated, the calculation method for the ratio of the area/volume to the contact number is as follows:
Figure FDA0002486424910000013
wherein 2D and 3D represent 2-dimensional and 3-dimensional models, respectively; phi is porosity;
Figure FDA0002486424910000027
is the average coordination number of the particles.
4. The automatic calibration method for the linear stiffness parameter of the discrete element during simulation of the geotechnical material according to claim 1, wherein the initial friction coefficient f during particle generation is adoptedcThe porosity and average coordination number empirical values of the final sample are estimated as follows:
Figure FDA0002486424910000021
。
5. the automatic calibration method for the linear stiffness parameters of the discrete elements during the simulation of the geotechnical materials according to claim 1, wherein the initial stiffness of the model is calculated by using the median particle size for the multi-particle size dispersion sample.
6. The automatic calibration method for the linear stiffness parameters of the discrete elements during simulation of the geotechnical materials according to claim 1, wherein after the true triaxial or biaxial simulation is finished, a stress-strain curve of the sample during loading is checked to ensure that the sample is only elastically deformed in a recoverable way during the loading period.
7. The automatic calibration method for the linear stiffness parameter of the discrete element during the simulation of the geotechnical material according to claim 1, wherein an error function L is as follows:
Figure FDA0002486424910000022
wherein, β1And β2Is two value ranges of [0,1 ]]And satisfy β1+β2＝1，E*And ν is the calculated elastic modulus and poisson's ratio, respectively, of the numerical samples.
8. The automatic calibration method for the linear stiffness parameter of the discrete element during simulation of the geotechnical material according to claim 1, wherein: if the particle contact stiffness parameter is updated for the first time, updating according to the following mode:
Figure FDA0002486424910000023
wherein, the symbol is assigned αnAnd αsIs coefficient, and the value is-0.5-0.5.
9. The automatic calibration method for the linear stiffness parameter of the discrete element during simulation of the geotechnical material according to claim 1, wherein if the current particle contact stiffness parameter is updated for the second time or more, the error function L to k is calculated according to a difference methodnAnd ksPartial derivatives of
Figure FDA0002486424910000024
And
Figure FDA0002486424910000025
if order is given by knAnd ksThe error function for the argument is f, let:
f(kn,ks)＝L(E*,ν*) (7)
the two partial derivatives are estimated differentially as follows:
Figure FDA0002486424910000026
。
10. the automatic calibration method for the linear stiffness parameter of the discrete element during the simulation of the geotechnical material according to claim 1, wherein if the current particle contact stiffness parameter is updated for the second time or more, the current particle contact stiffness parameter is updated according to an improved gradient descent method, and the update formula is as follows:
Figure FDA0002486424910000031
。
CN202010393302.5A 2020-05-11 2020-05-11 Automatic calibration method for discrete element linear stiffness parameter in simulation of rock and soil material Pending CN111611695A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010393302.5A CN111611695A (en) 2020-05-11 2020-05-11 Automatic calibration method for discrete element linear stiffness parameter in simulation of rock and soil material

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010393302.5A CN111611695A (en) 2020-05-11 2020-05-11 Automatic calibration method for discrete element linear stiffness parameter in simulation of rock and soil material

Publications (1)

Publication Number Publication Date
CN111611695A true CN111611695A (en) 2020-09-01

Family

ID=72200225

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010393302.5A Pending CN111611695A (en) 2020-05-11 2020-05-11 Automatic calibration method for discrete element linear stiffness parameter in simulation of rock and soil material

Country Status (1)

Country Link
CN (1) CN111611695A (en)

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112098209A (en) * 2020-09-15 2020-12-18 河北工业大学 Rock-soil particle damage localization identification method
CN112163328A (en) * 2020-09-18 2021-01-01 武汉大学 Geotechnical particle material constitutive modeling method based on deep learning and data driving
CN112765895A (en) * 2021-01-28 2021-05-07 南京大学 Machine learning-based automatic modeling method for discrete elements of rock and soil materials
CN113191004A (en) * 2021-05-09 2021-07-30 甘肃省地震局 Loess soil body internal rigidity calculation method
CN115828747A (en) * 2022-12-01 2023-03-21 山东大学 Fractured rock mass parameter intelligent calibration method and system considering particle interlocking effect

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20130289953A1 (en) * 2012-01-24 2013-10-31 The University Of Akron Self-optimizing, inverse analysis method for parameter identification of nonlinear material constitutive models
CN109030202A (en) * 2018-06-19 2018-12-18 湘潭大学 A kind of method of quick determining rock fragile materials discrete element analysis parameter
CN109543337A (en) * 2018-12-05 2019-03-29 陕西理工大学 A kind of rill evolution scaling method in discrete element simulation
CN109612885A (en) * 2019-01-08 2019-04-12 东北大学 A kind of mineral grain model parameter scaling method based on distinct element method
CN109815599A (en) * 2019-01-28 2019-05-28 南京大学 A kind of automatic training method of discrete element material

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20130289953A1 (en) * 2012-01-24 2013-10-31 The University Of Akron Self-optimizing, inverse analysis method for parameter identification of nonlinear material constitutive models
CN109030202A (en) * 2018-06-19 2018-12-18 湘潭大学 A kind of method of quick determining rock fragile materials discrete element analysis parameter
CN109543337A (en) * 2018-12-05 2019-03-29 陕西理工大学 A kind of rill evolution scaling method in discrete element simulation
CN109612885A (en) * 2019-01-08 2019-04-12 东北大学 A kind of mineral grain model parameter scaling method based on distinct element method
CN109815599A (en) * 2019-01-28 2019-05-28 南京大学 A kind of automatic training method of discrete element material

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
TONGMING QU等: "Calibration of linear contact stiffnesses in discrete element models using a hybrid analytical-computational framework", 《POWDER TECHNOLOGY》 *

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112098209A (en) * 2020-09-15 2020-12-18 河北工业大学 Rock-soil particle damage localization identification method
CN112163328A (en) * 2020-09-18 2021-01-01 武汉大学 Geotechnical particle material constitutive modeling method based on deep learning and data driving
CN112765895A (en) * 2021-01-28 2021-05-07 南京大学 Machine learning-based automatic modeling method for discrete elements of rock and soil materials
CN112765895B (en) * 2021-01-28 2023-10-17 南京大学 Automatic modeling method for discrete elements of rock and soil materials based on machine learning
CN113191004A (en) * 2021-05-09 2021-07-30 甘肃省地震局 Loess soil body internal rigidity calculation method
CN113191004B (en) * 2021-05-09 2022-08-30 甘肃省地震局 Method for calculating internal rigidity of loess soil body
CN115828747A (en) * 2022-12-01 2023-03-21 山东大学 Fractured rock mass parameter intelligent calibration method and system considering particle interlocking effect
CN115828747B (en) * 2022-12-01 2023-09-26 山东大学 Intelligent calibration method and system for fracture rock parameters by considering particle interlocking effect

Similar Documents

Publication Publication Date Title
CN111611695A (en) Automatic calibration method for discrete element linear stiffness parameter in simulation of rock and soil material
Shi et al. Calibration of micro-scaled mechanical parameters of granite based on a bonded-particle model with 2D particle flow code
Yu et al. An intelligent displacement back-analysis method for earth-rockfill dams
Armaghani et al. Indirect measure of shale shear strength parameters by means of rock index tests through an optimized artificial neural network
Kahraman et al. The usability of Cerchar abrasivity index for the prediction of UCS and E of Misis Fault Breccia: regression and artificial neural networks analysis
Kim et al. Prediction of subgrade resilient modulus using artificial neural network
CN111610091A (en) Automatic calibration method for discrete element Hertz contact parameter during simulation of geotechnical material
CN112163328B (en) Geotechnical particle material constitutive modeling method based on deep learning and data driving
Najjar et al. Simulating the stress–strain behavior of Georgia kaolin via recurrent neuronet approach
CN105259331A (en) Uniaxial strength forecasting method for jointed rock mass
JP6789274B2 (en) How to build a meso force chain particle model based on BPM theory
Kayadelen Estimation of effective stress parameter of unsaturated soils by using artificial neural networks
CN115508206B (en) Joint rock slope rock mass intensity parameter probability inversion method
Ma et al. Crack evolution and acoustic emission characteristics of rock specimens containing random joints under uniaxial compression
CN110929412B (en) Dynamic joint friction coefficient attenuation calculation method based on DDA theory
CN110457746B (en) BP neural network-based structural plane peak shear strength prediction model construction method
CN115641702B (en) Single landslide early warning and forecasting method based on space-time combination
Garg et al. A new simulation approach of genetic programming in modelling of soil water retention property of unsaturated soil
CN112884739B (en) Deep learning network-based method for rapidly detecting filling compactness of rock-fill body
CN110750901B (en) Discrete element model-based soil disturbance range judgment method
Li et al. Parameter Estimation of Particle Flow Model for Soils Using Neural Networks.
Kim et al. Preliminary study on PFC3D microparameter calibration using optimization of an artificial neural network
Moon et al. Non-equibiaxial residual stress evaluation methodology using simulated indentation behavior and machine learning
CN109960856A (en) A kind of false triaxial shear test discrete element numerical simulation method based on cylindrical body servo condition
Duffey et al. Non-normal impact of earth penetrators

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
WD01 Invention patent application deemed withdrawn after publication
WD01 Invention patent application deemed withdrawn after publication

Application publication date: 20200901