Rock mass single-bond two-parameter dynamic damage constitutive model construction and implementation method thereof

Abstract

The invention provides a construction and implementation method of a rock mass single-bond two-parameter dynamic damage constitutive model, relates to the technical field of tunnel and underground engineering rock mass mechanical calculation, and solves the attenuation effect of a non-local PD force by introducing a hexagonal polynomial attenuation function to reflect the influence of an internal length on a non-local long-range force; a new constitutive force function and a new failure criterion are constructed by introducing a damage function of a key; the influence of the key deformation rate on the critical elongation is considered to reflect the rate effect of the rock material; the defect of a fixed Poisson's ratio is overcome by introducing tangential rigidity into the bond, and meanwhile, the heterogeneity of the rock is described by utilizing normal rigidity and tangential rigidity obeying Weibull distribution, so that the defect that the heterogeneity of the rock material cannot be reflected by traditional bond-based near-field dynamics is overcome.

Description

Technical Field

[0001] The invention relates to the field of geotechnical engineering, and in particular to a construction method of a rock single-bond dual-parameter dynamic damage constitutive model and an implementation method thereof. Background Art

[0002] In underground engineering construction, rock mass is often subjected to dynamic loads such as blasting, drilling, and mining earthquakes. These dynamic loads may accelerate rock mass damage, causing cracks to expand and penetrate faster, thereby inducing rock mass instability and damage, affecting the stability of the rock mass, and threatening the lives and property safety of construction workers. Therefore, studying the mechanical properties and fracture mechanism of rock under dynamic loads can provide the most realistic and effective theoretical basis for the stability research of underground engineering.

[0003] As an emerging non-local theory, peridynamics uses integral form to solve the motion equations of solid mechanics, solving the problem that traditional methods have no definition in solving discontinuities. In peridynamics, there is no need to introduce additional fracture criteria, and cracks can be spontaneously simulated to start, expand, and bifurcate in any direction, which is widely used in the study of material fracture behavior.

[0004] The PD model originally proposed is called the Bond-Based Peridynamics (BBPD) model, which only considers the point force between two material points, so the Poisson's ratio is a constant. The bond-based peridynamics model describes brittle fracture and is not suitable for simulating the nonlinear mechanical behavior and progressive failure of rock materials. In addition, the traditional bond-based peridynamics model ignores the influence of internal length on nonlocal long-range forces, which affects the calculation accuracy. The mechanical properties and deformation characteristics of rock materials under impact loads show obvious rate effects. The traditional peridynamic model ignores the influence of strain rate on bond strength and cannot accurately reflect the deformation characteristics of rock materials under dynamic loads. As a natural medium, heterogeneity is its typical characteristic, while the traditional bond-based peridynamics model ignores the heterogeneity of rock. At present, some scholars have improved the bond-based peridynamics model and broken through the limitation of fixed Poisson's ratio, but have not considered the influence of strain rate on the dynamic mechanical properties of rock materials, nor the heterogeneity of rock. Summary of the invention

[0005] In view of the deficiencies in the prior art, the purpose of the present invention is to provide a single-bond dual-parameter rock dynamic damage constitutive model construction method and implementation method, which can effectively simulate the entire process of crack initiation, expansion and penetration of rock materials under impact loads.

[0006] In order to achieve the above object, the present invention adopts the following technical solution. The present invention provides a method for constructing a dynamic damage constitutive model of a peridynamic rock material, comprising the following steps:

[0007] S10 introduces an influence function that can reflect the change of the interaction intensity of material points with the distance between material points into the single-bond dual-parameter bond-based peridynamic constitutive force function, thus improving the constitutive force function.

[0008] S20 derives the strain energy density of the improved constitutive force function;

[0009] S30 The strain energy density of the model in the single-bond two-parameter bond-based peridynamics is equal to the strain energy density of the continuum mechanics, and the corrected normal stiffness and tangential stiffness are obtained;

[0010] S40 uses Weibull distribution to assign values ​​to the tangential stiffness and normal stiffness described in step S13.

[0011] S50 introduces the damage function of the bond and constructs a new fracture criterion;

[0012] S60 considers the influence of strain rate on the dynamic strength of the bond and constructs a dynamic failure criterion.

[0013] S70 substitutes the tangential stiffness and normal stiffness obtained in step S14 and the damage function of the key in step S15 into the improved constitutive force function in step S11 to complete the construction of the dynamic damage constitutive model of rock materials.

[0014] The expression form of the influence function described in step S10 is:

[0015]

[0016] Among them, δ is the field radius and ξ is the initial distance between material points.

[0017] The step S20 comprises the following steps:

[0018] S21 determines the expression form of the constitutive force function of the improved single-bond two-parameter bond-based peridynamic model.

[0019] For two-dimensional problems, the constitutive force function is expressed as:

[0020]

[0021] Among them, δ is the field radius, ξ is the initial distance between material points, c n (0,δ) is the axial stiffness, c t (0,δ) is the tangential stiffness, η n is the displacement vector η along The weight, η t is the displacement vector η along The weight, represents the unit vector along the initial direction of the bond, represents the unit vector perpendicular to the initial direction of the bond.

[0022] For three-dimensional problems, the constitutive force function is expressed as:

[0023]

[0024] Among them, δ is the field radius, ξ is the initial distance between material points, and c n (0,δ) is the axial stiffness, c t (0,δ) is the tangential stiffness, η n is the displacement vector η along The weight, η t1 is the displacement vector η along The weight, η t2 is the displacement vector η along The weight, represents the unit vector along the initial direction of the bond, and represents the unit vector perpendicular to the initial direction of the bond, and and Perpendicular to each other.

[0025] S22 Determination of the micropotential energy function of the improved single-bond two-parameter bond-based peridynamics.

[0026] For two-dimensional problems

[0027]

[0028] Among them, δ is the field radius, ξ is the initial distance between material points, and c n (0,δ) is the axial stiffness, c t (0,δ) is the tangential stiffness, η n / ξ and η t / ξ respectively represent the bond and Elongation in direction.

[0029] For three-dimensional problems

[0030]

[0031] Among them, δ is the field radius, ξ is the initial distance between material points, and c n (0,δ) is the axial stiffness, c t (0,δ) is the tangential stiffness, η n / ξ,η t1 / ξ and η t2 / ξ respectively represent the bond and Elongation in direction.

[0032] S23 integrates the micro potential energy function in step S12.2 to obtain the strain energy density.

[0033] For three-dimensional problems, the strain energy density is expressed as

[0034]

[0035] Among them, δ is the field radius, c n (0,δ) is the axial stiffness, c t (0,δ) is the tangential stiffness. ε ij (i=1,2,3; j=1,2,3) represents strain.

[0036] For a two-dimensional problem, the strain energy density is expressed as

[0037]

[0038] Among them, h is the thickness, δ is the field radius, c n (0,δ) is the axial stiffness, c t (0,δ) is the tangential stiffness, ε ij (i=1,2; j=1,2) represents strain.

[0039] The normal stiffness and tangential stiffness in step S30 can be obtained by derivation:

[0040] Under plane stress

[0041]

[0042] Among them, h is the thickness, δ is the field radius, E is the elastic modulus, and ν is the Poisson's ratio.

[0043] In plane strain state

[0044]

[0045] Among them, h is the thickness, δ is the field radius, E is the elastic modulus, and ν is the Poisson's ratio.

[0046] For three-dimensional problems

[0047]

[0048] Among them, δ is the domain radius, E is the elastic modulus, and ν is the Poisson's ratio.

[0049] The normal stiffness and tangential stiffness considering the Weibull distribution described in step S40 can be expressed as:

[0050]

[0051] Among them, c n(0,δ) is the axial stiffness, c t (0,δ) is the tangential stiffness, r(x) represents a random parameter uniformly distributed between (0,1], β1 and β2 represent c nw and c tw Uniformity level index.

[0052] The damage function of the key described in step S50 can be expressed as:

[0053]

[0054] in, represents the critical elongation of bond elasticity, represents the critical elongation of key elastic compression, represents the critical elongation of bond tension, represents the critical elongation of key compression, and m represents the softening coefficient, which is used to control the attenuation of key force.

[0055] The failure criterion after introducing the damage function can be expressed as:

[0056] Critical elongation of bond tension

[0057]

[0058] Key compression critical elongation

[0059]

[0060] Among them, α and β represent the microscopic parameters in the dynamic damage constitutive model of rock materials, α represents the microscopic parameter matching the fracture energy, and β represents the microscopic parameter matching the uniaxial compression absorption energy of rock materials.

[0061] When considering the effect of strain rate on the dynamic strength of the bond in step S60, the critical elongation of bond elastic tension is and elastic compression critical elongation It can be expressed as:

[0062]

[0063]

[0064] in, represents the tensile deformation rate, Indicates the compression deformation rate. represents the quasi-static strain rate under tension, represents the quasi-static strain rate under compression. Take 3×10 -6 s -1 ,30×10-6 s -1 .

[0065] For the plane stress problem, the dynamic damage constitutive model of rock material described in step S70 can be expressed as:

[0066]

[0067] Among them, E represents the elastic modulus; δ represents the field radius; ν is the Poisson's ratio; ξ represents the initial distance between material points; represents the critical elongation of bond elastic tension, and h represents the thickness; represents the critical elongation of key elastic compression, represents the critical elongation of bond tension, represents the critical elongation of the key compression, m represents the softening coefficient, which is used to control the attenuation of the key force, and c n (0,δ) is the axial stiffness, c t (0,δ) is the tangential stiffness, r(x) represents a random parameter uniformly distributed between (0,1], β1 and β2 represent c nw and c tw Uniformity level index; η n is the displacement vector η along The weight, η t is the displacement vector η along The weight, represents the unit vector along the initial direction of the bond, represents the unit vector perpendicular to the initial direction of the bond.

[0068] For the plane strain problem, the dynamic damage constitutive model of rock material described in step S17 can be expressed as

[0069]

[0070] Among them, E represents the elastic modulus; δ represents the field radius; ν is the Poisson's ratio; ξ represents the initial distance between material points; represents the critical elongation of bond elastic tension, and h represents the thickness; represents the critical elongation of key elastic compression, represents the critical elongation of bond tension, represents the critical elongation of the key compression, m represents the softening coefficient, which is used to control the attenuation of the key force, and c n (0,δ) is the axial stiffness, c t (0,δ) is the tangential stiffness, r(x) represents a random parameter uniformly distributed between (0,1], β1 and β2 represent c nw and c tw Uniformity level index; η n is the displacement vector η along The weight, η t is the displacement vector η along The weight, represents the unit vector along the initial direction of the bond, represents the unit vector perpendicular to the initial direction of the bond.

[0071] For three-dimensional problems, the dynamic damage constitutive model of rock materials described in step S17 can be expressed as

[0072]

[0073] Among them, E represents the elastic modulus; δ represents the field radius; ν is the Poisson's ratio; ξ represents the initial distance between material points; It represents the critical elongation of bond elastic tension; represents the critical elongation of key elastic compression, represents the critical elongation of bond tension, represents the critical elongation of the key compression, m represents the softening coefficient, which is used to control the attenuation of the key force, and c n (0,δ) is the axial stiffness, c t (0,δ) is the tangential stiffness, r(x) represents a random parameter uniformly distributed between (0,1], β1 and β2 represent c nw and c tw Uniformity level index; η n is the displacement vector η along The weight, η t1 is the displacement vector η along The weight, η t2 is the displacement vector η along The weight, represents the unit vector along the initial direction of the bond, and represents the unit vector perpendicular to the initial direction of the bond, and and Perpendicular to each other.

[0074] The present invention also provides a method for implementing a single-bond dual-parameter rock dynamic damage constitutive model, comprising the following steps:

[0075] S101: Parameter initialization, determining model geometry parameters, and defining material parameters;

[0076] S102: Discretize the solution domain and generate material point coordinates;

[0077] S103: Determine the points within the material point domain, and initialize the keys within the material point domain;

[0078] S104: Assign values ​​of the axial stiffness and tangential stiffness of each material point according to Weibull distribution, and take the average value to calculate the axial stiffness and normal stiffness of each key;

[0079] S105: Determine the time step;

[0080] S106: Apply velocity and displacement boundary conditions;

[0081] S107: Calculate the initial orientation angle, axial deformation and normal deformation of the key.

[0082] S108: Calculate volume correction factor;

[0083] S109: Calculate the tensile and compressive deformation rates of the key

[0084] S110: Determine the critical elongation of the key;

[0085] S111: Calculate the constitutive forces of material points;

[0086] S112: Calculate the elongation of the bond and determine whether the elongation is greater than the critical elongation or less than the critical compression rate. If "yes", the bond between the material points is broken and the constitutive force between the material points is 0. If "no", calculate the constitutive force according to step 11.

[0087] S113: Use forward difference technique to solve displacement and velocity.

[0088] S114: Output the calculation results and draw a damage cloud map.

[0089] Compared with the prior art, the present invention has the following beneficial effects:

[0090] (1) A sixth-order polynomial attenuation function is introduced to reflect the influence of internal length on nonlocal long-range force, solve the attenuation effect of nonlocal PD force, and improve the calculation accuracy.

[0091] (2) The traditional bond-based peridynamic model describes brittle fracture rather than quasi-brittle fracture. When the bond elongation exceeds the critical elongation, the force between the material points disappears and the bond force suddenly drops to 0. Obviously, this failure criterion cannot describe the nonlinear mechanical behavior and progressive failure process of rock materials. To this end, the damage function of the bond is introduced to construct a new constitutive force function and failure criterion.

[0092] (3) The strength of rock materials has obvious strain rate effect. The traditional bond-based peridynamic model is not suitable for simulating the dynamic fracture process of rock under impact load. The rate effect of rock materials can be reflected by considering the influence of bond deformation rate on critical elongation.

[0093] (4) The defect of fixed Poisson's ratio is overcome by introducing tangential stiffness in the bond. At the same time, the normal stiffness and tangential stiffness obeying Weibull distribution are used to describe the heterogeneity of rock, which makes up for the deficiency of traditional bond-based peridynamics that cannot reflect the heterogeneity of rock materials. BRIEF DESCRIPTION OF THE DRAWINGS

[0094] Figure 1 Schematic diagram of the motion principle of material points in the near-field dynamics in the construction method of the single-bond dual-parameter dynamic damage constitutive model of rock mass provided in the embodiment

[0095] Figure 2 The relationship between the bond force and elongation of the PMB material in the construction method of the single-bond dual-parameter dynamic damage constitutive model of rock mass provided in the embodiment;

[0096] Figure 3 The relationship between the bond force and elongation of the rock material in the construction method of the single-bond dual-parameter dynamic damage constitutive model of the rock mass provided in the embodiment;

[0097] Figure 4 A calculation flow chart of a single-bond dual-parameter rock material dynamic damage constitutive model in a rock mass detection method provided in an embodiment;

[0098] Figure 5 A rock crack propagation path diagram under impact load in a rock mass detection method provided in an embodiment;

[0099] Figure 6 A numerical model diagram of rock dynamic impact in the rock mass detection method provided in the embodiment;

[0100] Figure 7 A random distribution diagram of bond stiffness in a rock mass detection method provided in an embodiment. DETAILED DESCRIPTION

[0101] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0102] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the accompanying drawings in the embodiments of the present invention.

[0103] Example 1

[0104] like Figure 1As shown in the bond-based peridynamics model, there is an interaction between the material point x and the material point x′ in its field. According to Newton's second law, the motion equation of the material point x at any time can be expressed as:

[0105]

[0106] ρ(x) represents the mass density of material point x, u(x,t) and u(x′,t) represent the displacements of material points x and x′ at time t, respectively. represents the acceleration of the material point x, H x ={x′∈R,||x′-x||<δ} represents the set of other material points within the domain of material point x, f represents the point force function (also known as the constitutive force function), b(x, t) represents the body force density, V x ' is the volume element of the material point x';

[0107] The traditional bond-based peridynamic model has defects in simulating the fracture of rock materials under impact loads. As a quasi-brittle material, rock exhibits linear and nonlinear mechanical behaviors under external loads. Figure 2 , 3 As shown, in the traditional bond-based peridynamic model, when the elongation of the bond exceeds the critical elongation, the force between the material points disappears and the bond force suddenly drops to 0. Obviously, this failure criterion cannot describe the nonlinear mechanical behavior and progressive failure process of rock materials. The traditional bond-based peridynamic model only considers the point force between two material points, so the Poisson's ratio is a constant. Moreover, the traditional bond-based peridynamic model does not consider the influence of strain rate on bond strength. Accordingly, the present invention proposes a method for constructing a single-bond dual-parameter rock dynamic damage constitutive model, comprising the following steps:

[0108] S10 introduces an influence function that can reflect the change of the interaction intensity of material points with the distance between material points into the single-bond dual-parameter bond-based peridynamic constitutive force function, thus improving the constitutive force function.

[0109] S20 derives the strain energy density of the improved constitutive force function;

[0110] S30 The strain energy density of the model in the single-bond two-parameter bond-based peridynamics is equal to the strain energy density of the continuum mechanics, and the corrected normal stiffness and tangential stiffness are obtained;

[0111] S40 uses Weibull distribution to assign values ​​to the tangential stiffness and normal stiffness described in step S30.

[0112] S50 introduces the damage function of the bond and constructs a new fracture criterion;

[0113] S60 considers the influence of strain rate on the dynamic strength of the bond and constructs a dynamic failure criterion.

[0114] S70 substitutes the tangential stiffness and normal stiffness obtained in step S14 and the damage function of the key in step S50 into the improved constitutive force function in step S10 to complete the construction of the dynamic damage constitutive model of rock materials.

[0115] The expression form of the influence function described in step S10 is:

[0116]

[0117] Among them, δ is the field radius and ξ is the initial distance between material points.

[0118] The step S20 comprises the following steps:

[0119] S21 determines the expression form of the constitutive force function of the improved single-bond two-parameter bond-based peridynamic model.

[0120] For two-dimensional problems, the constitutive force function is expressed as:

[0121]

[0122] Among them, δ is the field radius, ξ is the initial distance between material points, c n (0,δ) is the axial stiffness, c t (0,δ) is the tangential stiffness, η n is the displacement vector η along The weight, η t is the displacement vector η along The weight, represents the unit vector along the initial direction of the bond, represents the unit vector perpendicular to the initial direction of the bond.

[0123] For three-dimensional problems, the constitutive force function is expressed as:

[0124]

[0125] Among them, δ is the field radius, ξ is the initial distance between material points, c n (0,δ) is the axial stiffness, c t (0,δ) is the tangential stiffness, η n is the displacement vector η along The weight, η t1 is the displacement vector η along The weight, η t2 is the displacement vector η along The weight, represents the unit vector along the initial direction of the bond, and represents the unit vector perpendicular to the initial direction of the bond, and and Perpendicular to each other.

[0126] S22 Determination of the micropotential energy function of the improved single-bond two-parameter bond-based peridynamics.

[0127] For two-dimensional problems

[0128]

[0129] Among them, δ is the field radius, ξ is the initial distance between material points, and c n (0,δ) is the axial stiffness, c t (0,δ) is the tangential stiffness, η n / ξ and η t / ξ respectively represent the bond and Elongation in direction.

[0130] For three-dimensional problems

[0131]

[0132] Among them, δ is the field radius, ξ is the initial distance between material points, and c n (0,δ) is the axial stiffness, c t (0,δ) is the tangential stiffness, η n / ξ,η t1 / ξ and η t2 / ξ respectively represent the bond and Elongation in direction.

[0133] S23 integrates the micro potential energy function in step S12.2 to obtain the strain energy density.

[0134] For three-dimensional problems, the strain energy density is expressed as

[0135]

[0136] Among them, δ is the field radius, c n (0,δ) is the axial stiffness, c t (0,δ) is the tangential stiffness. ε ij (i=1,2,3; j=1,2,3) represents strain.

[0137] For a two-dimensional problem, the strain energy density is expressed as

[0138]

[0139] Among them, h is the thickness, δ is the field radius, c n (0,δ) is the axial stiffness, c t (0,δ) is the tangential stiffness, ε ij (i=1,2; j=1,2) represents strain.

[0140] The normal stiffness and tangential stiffness in step S30 can be obtained by derivation:

[0141] Under plane stress

[0142]

[0143] Among them, h is the thickness, δ is the field radius, E is the elastic modulus, and ν is the Poisson's ratio.

[0144] In plane strain state

[0145]

[0146] Among them, h is the thickness, δ is the field radius, E is the elastic modulus, and ν is the Poisson's ratio.

[0147] For three-dimensional problems

[0148]

[0149] Among them, δ is the domain radius, E is the elastic modulus, and ν is the Poisson's ratio.

[0150] The normal stiffness and tangential stiffness considering the Weibull distribution described in step S40 can be expressed as:

[0151]

[0152] Among them, c n (0,δ) is the axial stiffness, c t (0,δ) is the tangential stiffness, r(x) represents a random parameter uniformly distributed between (0,1], β1 and β2 represent c nw and c tw Uniformity level index.

[0153] The damage function of the key described in step S50 can be expressed as:

[0154]

[0155] in, represents the critical elongation of bond elasticity, represents the critical elongation of key elastic compression, represents the critical elongation of bond tension, represents the critical elongation of key compression, and m represents the softening coefficient, which is used to control the attenuation of key force.

[0156] The failure criterion after introducing the damage function can be expressed as:

[0157] Critical elongation of bond tension

[0158]

[0159] Key compression critical elongation

[0160]

[0161] Among them, α and β represent the microscopic parameters in the dynamic damage constitutive model of rock materials, α represents the microscopic parameter matching the fracture energy, and β represents the microscopic parameter matching the uniaxial compression absorption energy of rock materials.

[0162] When considering the influence of the bond deformation rate on the bond dynamic strength in step S60, the bond elastic tensile critical elongation and elastic compression critical elongation It can be expressed as:

[0163]

[0164] in, represents the tensile deformation rate, Indicates the compression deformation rate. represents the quasi-static strain rate under tension, represents the quasi-static strain rate under compression. Take 3×10 -6 s -1 ,30×10 -6 s -1 .

[0165] The deformation rate of the bond is similar to the concept of strain rate in continuum mechanics and can be expressed as:

[0166]

[0167] Where ξ is the initial distance between material points, and represent the velocities of material point j and material point i respectively.

[0168] For the plane stress problem, the dynamic damage constitutive model of rock material described in step S70 can be expressed as:

[0169]

[0170] Among them, E represents the elastic modulus; δ represents the field radius; ν is the Poisson's ratio; ξ represents the initial distance between material points; represents the critical elongation of bond elastic tension, and h represents the thickness; represents the critical elongation of key elastic compression, represents the critical elongation of bond tension, represents the critical elongation of the key compression, m represents the softening coefficient, which is used to control the attenuation of the key force, and c n (0,δ) is the axial stiffness, c t (0,δ) is the tangential stiffness, r(x) represents a random parameter uniformly distributed between (0,1], β1 and β2 represent c nw and c tw Uniformity level index; η n is the displacement vector η along The weight, η t is the displacement vector η along The weight, represents the unit vector along the initial direction of the bond, represents the unit vector perpendicular to the initial direction of the bond.

[0171] For the plane strain problem, the dynamic damage constitutive model of rock material described in step S17 can be expressed as

[0172]

[0173] Among them, E represents the elastic modulus; δ represents the field radius; ν is the Poisson's ratio; ξ represents the initial distance between material points; represents the critical elongation of bond elastic tension, and h represents the thickness; represents the critical elongation of key elastic compression, represents the critical elongation of bond tension, represents the critical elongation of the key compression, m represents the softening coefficient, which is used to control the attenuation of the key force, and c n (0,δ) is the axial stiffness, c t (0,δ) is the tangential stiffness, r(x) represents a random parameter uniformly distributed between (0,1], β1 and β2 represent c nw and c tw Uniformity level index; η n is the displacement vector η along The weight, η t is the displacement vector η along The weight, represents the unit vector along the initial direction of the bond, represents the unit vector perpendicular to the initial direction of the bond.

[0174] For a three-dimensional problem, the dynamic damage constitutive model of rock materials described in step S70 can be expressed as:

[0175]

[0176] Among them, E represents the elastic modulus; δ represents the field radius; ν is the Poisson's ratio; ξ represents the initial distance between material points; It represents the critical elongation of bond elastic tension; represents the critical elongation of key elastic compression, represents the critical elongation of bond tension, represents the critical elongation of the key compression, m represents the softening coefficient, which is used to control the attenuation of the key force, and c n (0,δ) is the axial stiffness, c t (0,δ) is the tangential stiffness, r(x) represents a random parameter uniformly distributed between (0,1], β1 and β2 represent c nw and c tw Uniformity level index; η n is the displacement vector η along The weight, η t1 is the displacement vector η along The weight, η t2 is the displacement vector η along The weight, represents the unit vector along the initial direction of the bond, and represents the unit vector perpendicular to the initial direction of the bond, and and Perpendicular to each other.

[0177] Example 2

[0178] like Figure 4 As shown, the present invention also provides an implementation method of a single-bond dual-parameter rock dynamic damage constitutive model, comprising the following steps:

[0179] S101: Parameter initialization, determining model geometry parameters, and defining material parameters;

[0180] Taking the three-dimensional problem as an example, the model length is defined as L, the width is W, the thickness is H, the elastic modulus of the material is E, the density is ρ, the Poisson's ratio is υ, and the tensile strength is f t , tensile strength f c .

[0181] S102: Discretize the solution domain, determine the material point spacing, and generate material point coordinates;

[0182] The model is L long, W wide, and H thick. The number of material points distributed along the length direction is m, the number distributed along the width direction is n, the number distributed along the thickness direction is p, and the spacing between material points is

[0183] The length direction is parallel to the x-axis, the width direction is parallel to the y-axis, and the thickness direction is parallel to the z-axis. The coordinates of the material point can be expressed as:

[0184]

[0185] Wherein, i = 1, 2, ..., m, represents the number of material points arranged along the length direction. j = 1, 2, ..., n, represents the number of material points arranged along the width direction. k = 1, 2, ..., p, represents the number of material points arranged along the thickness direction.

[0186] S103: Determine the points within the material point domain, and initialize the keys within the material point domain;

[0187] Traverse all material points, calculate the distance between material points x and x′, and determine whether the distance between material points satisfies H={x′∈R,||x′-x||<δ}. If so, number the material points. Construct the key in the material point domain and initialize it, setting μ(x,t,ξ)=1.

[0188] Among them, μ(x, t, ξ) is a scalar function used to judge the integrity of the bond. If u(x, t, ξ) = 1, the bond is broken, and if u(x, t, ξ) = 0, the bond is broken.

[0189] Damage at one point It can be calculated by the ratio of the number of broken bonds to the total number of bonds within the material point area:

[0190]

[0191] S104: Assign values ​​of the axial stiffness and tangential stiffness of each material point according to Weibull distribution, and take the average value to calculate the axial stiffness and normal stiffness of each key;

[0192] Assume that the axial stiffness follows a Weibull distribution:

[0193]

[0194] Among them, c n (0,δ) represents the axial stiffness, represents the average value of the axial stiffness, and β1 is the uniformity level index of the axial stiffness obeying the Weibull distribution.

[0195] The axial stiffness that obeys the Weibull distribution can be obtained by inverse transformation:

[0196]

[0197] Among them, c n (0,δ) represents the axial stiffness, represents the average value of the axial stiffness, and β1 is the uniformity level index of the axial stiffness obeying the Weibull distribution.

[0198] Then the axial stiffness of the bond between material points x and x′ can be expressed as:

[0199]

[0200] Assume that the tangential stiffness follows a Weibull distribution:

[0201]

[0202] Among them, c t (0,δ) represents the tangential stiffness, represents the average value of the tangential stiffness, and β2 is the uniformity level index of the tangential stiffness following the Weibull distribution.

[0203] Then the tangential stiffness obeying the Weibull distribution can be expressed as:

[0204]

[0205] Among them, c t (0,δ) represents the tangential stiffness, represents the average value of the tangential stiffness, and β2 is the uniformity level index of the tangential stiffness following the Weibull distribution.

[0206] The tangential stiffness of the bond between material points x and x′ is obtained by averaging the axial stiffness at the two points:

[0207]

[0208] S105: Determine the time step required for the calculation model.

[0209] S106: Apply velocity and displacement boundary conditions;

[0210] S107: Calculate the initial orientation angle, axial deformation and normal deformation of the key.

[0211] S108: Calculate volume correction factor;

[0212] For material points at the boundary of the domain, their volumes need to be corrected. The volume correction factor can be calculated by the following formula:

[0213]

[0214] Wherein, r=Δx / 2, Δx represents the distance between material points.

[0215] S109: Calculate the strain rate of the key.

[0216] The strain rate of the bond can be expressed as:

[0217]

[0218] S110: Update rate related parameters. Calculate the critical elongation of the bond elasticity according to the current bond strain rate. Key elastic compression critical elongation Critical elongation of bond tension Key compression critical elongation

[0219] S111: Calculate the constitutive force of material points; taking a three-dimensional problem as an example, the constitutive force of the bond between material points at any time can be expressed as:

[0220]

[0221] The constitutive force of any material point can be expressed by the sum of the bond forces of all bonds within its domain:

[0222]

[0223] Wherein, N represents the number of material points within the domain of material point x.

[0224] S112: Calculate the elongation of the bond and determine whether the elongation is greater than the critical elongation or less than the critical compression rate. If "yes", the bond between the material points is broken and the constitutive force between the material points is 0. If "no", calculate the constitutive force according to step A11.

[0225] The elongation of the bond can be calculated by the following formula:

[0226]

[0227] S113: Use forward difference technique to solve displacement and velocity.

[0228] The velocity and displacement can be solved by the following equation:

[0229]

[0230] S114: Output the calculation results and draw a damage cloud map.

[0231] Figure 5The present invention is used to simulate the calculation effect diagram of the dynamic fracture of rock under impact load. The characteristic is that based on the simulation method of the dynamic fracture process of rock materials constructed by the present invention, a numerical model of rock under impact load is established. The implementation method of the dynamic fracture process of rock materials constructed by the present invention can simulate the dynamic fracture process and failure characteristics of rock under impact load. Due to the introduction of randomly distributed defects, the number of random microcracks is greater. The introduction of randomly distributed defects induces stress concentration, making the crack initiation position more random. As the load increases, microcracks begin to expand, the microcracks are interconnected, cracks cross each other, and the specimen is divided into multiple blocks, which are finally crushed.

[0232] Furthermore, the specific implementation steps of the calculation effect include:

[0233] B1 parameter initialization, determine the model geometry parameters, and define the material parameters.

[0234] In this embodiment, if Figure 6 As shown, the length of the rock specimen model is L = 50 mm, the diameter is D = 50 mm, the elastic modulus is E = 32.13 GPa, and the mass density is ρ = 2560 kg / m 3 , Poisson's ratio ν = 0.25, tensile strength f t =7.46MPa, compression strength f c =107MPa. The incident rod is 2000mm long, the transmission rod is 2000mm long, and the elastic modulus is E=210GPa.

[0235] In this embodiment, the total number of material points of the specimen is 401920, and the specimen is divided into 80 material points along the x-axis, 80 material points along the y-axis, and 80 material points along the z-axis. The material point spacing of the specimen is Δx = 0.625 mm, and the field radius is δ = 1.875 mm. The total number of material points of the incident rod is 2000, and the incident rod is simplified into a one-dimensional model, and 2000 material points are divided along the x-axis. The material point spacing of the incident rod is Δx = 1 mm, and the field radius is δ = 3 mm. The total number of material points of the transmission rod is 2000, and the transmission rod is simplified into a one-dimensional model, and 2000 material points are divided along the x-axis. The material point spacing of the transmission rod is Δx = 1 mm, and the field radius is δ = 3 mm.

[0236] B2 discretizes the solution domain and generates material point coordinates.

[0237] In this embodiment, the coordinates of each material point along the x direction can be expressed as follows, sorted from left to right:

[0238] x=-25+(i-1)+0.3125

[0239] Wherein, i=1, 2, ..., 80, indicating the number of lattices arranged along the x direction.

[0240] Arranged from bottom to top, the coordinate of each material point along the y direction can be expressed as;

[0241] y=-25+(j-1)+0.3125

[0242] Wherein, j=1, 2, ..., 80, indicating the number of lattices arranged along the y direction.

[0243] Sorting from front to back, the coordinate of each material point along the z direction can be expressed as

[0244] z=-25+(k-1)+0.3125

[0245] Wherein, k = 1, 2, ..., 80, indicating the number of lattices arranged along the z direction.

[0246] B3: Determine the points within the material point domain and initialize the keys within the material point domain;

[0247] Traverse all material points, calculate the distance between material points x and x′, and determine whether the distance between material points satisfies H={x′∈R,||x′-x||<δ}. If so, number the material points. Construct the key in the material point domain and initialize it, setting μ(x,t,ξ)=1.

[0248] B4: Assign the normal stiffness and tangential stiffness of each material point according to the Weibull distribution.

[0249] like Figure 7 As shown, in this embodiment, the average value of the normal stiffness is The shape parameter β1 of the normal stiffness obeying the Weibull distribution is 2.

[0250] B5: Apply boundary conditions.

[0251] Apply stress boundary conditions at the ends of the incident rod.

[0252] B6: Calculate the initial orientation angle, axial deformation and normal deformation of the key.

[0253] B7: Calculate the strain rate of the key and update rate related parameters.

[0254] The strain rate of the bond can be expressed as:

[0255]

[0256] B8: Calculate the constitutive force of the material point and solve the velocity and displacement of the material point.

[0257] The constitutive forces of all bonds in the near field of each material point are calculated and summed to obtain the constitutive force of a single material point. The velocity and displacement are solved by using the explicit forward and backward difference formulas in the time domain. The time step is 2.0×10 -7 s.

[0258] B9: Calculate the damage of material points.

[0259] Calculate the bond elongation in the field of each material point and determine whether the bond elongation is greater than the critical elongation. If so, the bonds between the material points are broken, that is, u(x, t, ξ) = 0. Record the number of broken bonds and obtain the damage of the material point by the ratio of the number of broken bonds to the total number of bonds in the near field.

[0260] The above are only specific embodiments of the present invention, but the protection scope of the present invention is not limited thereto. Any person skilled in the art who is familiar with the technical field can easily think of changes or substitutions within the technical scope disclosed by the present invention, which should be included in the protection scope of the present invention. Therefore, the protection scope of the present invention should be based on the protection scope of the claims.

Claims

1. A method for constructing a single-bond dual-parameter dynamic damage constitutive model for rock mass, characterized in that: The method comprises the following steps: S10 introduces an influence function that can reflect the change of the interaction intensity of material points with the distance between material points into the single-bond dual-parameter bond-based peridynamic constitutive force function, thus improving the constitutive force function. S20 derives the strain energy density of the improved constitutive force function; S30 The strain energy density of the model in the single-bond two-parameter bond-based peridynamics is equal to the strain energy density of the continuum mechanics, and the corrected normal stiffness and tangential stiffness are obtained; S40 uses Weibull distribution to assign values ​​to the tangential stiffness and normal stiffness described in step S30; S50 introduces the damage function of the bond and constructs a new fracture criterion; S60 considers the influence of strain rate on the dynamic strength of the bond and constructs the dynamic failure criterion; S70 substitutes the tangential stiffness and normal stiffness obtained in step S40 and the damage function of the key in step S50 into the improved constitutive force function in step S10 to complete the construction of the dynamic damage constitutive model of rock materials.

2. The method for constructing a single-bond dual-parameter dynamic damage constitutive model of rock mass according to claim 1, characterized in that: The expression of the influence function of step S10 is: Among them, δ is the field radius and ξ is the initial distance between material points.

3. The method for constructing a single-bond dual-parameter dynamic damage constitutive model of rock mass according to claim 1, characterized in that: Step S20 includes the following steps: S21. Determine the expression of the constitutive force function of the improved single-bond dual-parameter bond-based peridynamic model: For two-dimensional problems, the constitutive force function is expressed as: Among them, δ is the field radius, ξ is the initial distance between material points, c n (0,δ) is the axial stiffness, c t (0,δ) is the tangential stiffness, η n is the displacement vector η along The weight, η t is the displacement vector η along The weight, represents the unit vector along the initial direction of the bond, represents the unit vector perpendicular to the initial direction of the bond; For three-dimensional problems, the expression of the constitutive force function is: Among them, δ is the field radius, ξ is the initial distance between material points, c n (0,δ) is the axial stiffness, c t (0,δ) is the tangential stiffness, η n is the displacement vector η along The weight, η t1 is the displacement vector η along The weight, η t2 is the displacement vector η along The weight, represents the unit vector along the initial direction of the bond, and represents the unit vector perpendicular to the initial direction of the bond, and and perpendicular to each other; S22. Determine the micropotential energy function of the improved single-bond two-parameter bond-based peridynamics; For two-dimensional problems, the expression of the micro potential energy function is: Among them, δ is the field radius, ξ is the initial distance between material points, c n (0,δ) is the axial stiffness, c t (0,δ) is the tangential stiffness, η n / ξ and η t / ξ respectively represent the bond and Elongation in the direction; For three-dimensional problems, the expression of the micro potential energy function is: Among them, δ is the field radius, ξ is the initial distance between material points, c n (0,δ) is the axial stiffness, c t (0,δ) is the tangential stiffness, η n / ξ,η t1 / ξ and η t2 / ξ respectively represent the bond and Elongation in the direction; S23. Integrate the micro potential energy function in step S20 to obtain the strain energy density; For three-dimensional problems, the expression for strain energy density is: Among them, δ is the field radius, c n (0,δ) is the axial stiffness, c t (0,δ) is the tangential stiffness; ε ij (i=1,2,3; j=1,2,3) represents strain; For two-dimensional problems, the expression for strain energy density is: Among them, h is the thickness, δ is the field radius, c n (0,δ) is the axial stiffness, c t (0,δ) is the tangential stiffness, ε ij (i=1,2; j=1,2) represents strain.

4. The method for constructing a single-bond dual-parameter dynamic damage constitutive model of rock mass according to claim 1, characterized in that: The normal stiffness and tangential stiffness in step S30 are derived as follows: Under plane stress Among them, h is the thickness, δ is the field radius, E is the elastic modulus, and ν is the Poisson's ratio; In plane strain state Among them, h is the thickness, δ is the field radius, E is the elastic modulus, and ν is the Poisson's ratio; For three-dimensional problems Among them, δ is the domain radius, E is the elastic modulus, and ν is the Poisson's ratio.

5. The method for constructing a single-bond dual-parameter dynamic damage constitutive model of rock mass according to claim 1, characterized in that: The expressions of the normal stiffness and tangential stiffness of the Weibull distribution in step S40 are: Among them, c n (0,δ) is the axial stiffness, c t (0,δ) is the tangential stiffness, r(x) represents a random parameter uniformly distributed between (0,1], β1 and β2 represent c nw and c tw Uniformity level index.

6. The method for constructing a single-bond dual-parameter dynamic damage constitutive model of rock mass according to claim 1, characterized in that: The expression of the damage function in step S50 is: in, represents the critical elongation of bond elasticity, represents the critical elongation of key elastic compression, represents the critical elongation of bond tension, represents the critical elongation of key compression, and m represents the softening coefficient; The failure criterion after introducing the damage function can be expressed as: Critical elongation of bond tension Key compression critical elongation Among them, α and β represent the microscopic parameters in the dynamic damage constitutive model of rock materials, α represents the microscopic parameter matching the fracture energy, and β represents the microscopic parameter matching the uniaxial compression absorption energy of rock materials.

7. The method for constructing a single-bond dual-parameter dynamic damage constitutive model of rock mass according to claim 1, characterized in that: In step S60, when considering the influence of the bond deformation rate on the bond dynamic strength, the bond elastic tensile critical elongation and elastic compression critical elongation It can be expressed as: in, represents the tensile deformation rate, Indicates the compression deformation rate; represents the quasi-static strain rate under tension, represents the quasi-static strain rate under compression, Take 3×10 -6 s -1 ,30×10 -6 s -1 .

8. The method for constructing a single-bond dual-parameter dynamic damage constitutive model of rock mass according to claim 1, characterized in that: In step S70, the expression of the dynamic damage constitutive model of rock material is: Among them, E represents the elastic modulus; δ represents the field radius; ν is the Poisson's ratio; ξ represents the initial distance between material points; represents the critical elongation of bond elastic tension, and h represents the thickness; represents the critical elongation of key elastic compression, represents the critical elongation of bond tension, represents the critical elongation of the key compression, m represents the softening coefficient, which is used to control the attenuation of the key force, and c n (0,δ) is the axial stiffness, c t (0,δ) is the tangential stiffness, r(x) represents a random parameter uniformly distributed between (0,1], β1 and β2 represent c nw and c tw Uniformity level index; η n is the displacement vector η along The weight, η t is the displacement vector η along The weight, represents the unit vector along the initial direction of the bond, represents the unit vector perpendicular to the initial direction of the bond; For the plane strain problem, the dynamic damage constitutive model of rock material described in step S17 can be expressed as: Among them, E represents the elastic modulus; δ represents the field radius; ν is the Poisson's ratio; ξ represents the initial distance between material points; represents the critical elongation of bond elastic tension, and h represents the thickness; represents the critical elongation of key elastic compression, represents the critical elongation of bond tension, represents the critical elongation of the key compression, m represents the softening coefficient, which is used to control the attenuation of the key force, and c n (0,δ) is the axial stiffness, c t (0,δ) is the tangential stiffness, r(x) represents a random parameter uniformly distributed between (0,1], β1 and β2 represent c nw and c tw Uniformity level index; η n is the displacement vector η along The weight, η t is the displacement vector η along The weight, represents the unit vector along the initial direction of the bond, represents the unit vector perpendicular to the initial direction of the bond; For three-dimensional problems, the dynamic damage constitutive model of rock materials described in step S17 can be expressed as: Among them, E represents the elastic modulus; δ represents the field radius; ν is the Poisson's ratio; ξ represents the initial distance between material points; It represents the critical elongation of bond elastic tension; represents the critical elongation of key elastic compression, represents the critical elongation of bond tension, represents the critical elongation of the key compression, m represents the softening coefficient, which is used to control the attenuation of the key force, and c n (0,δ) is the axial stiffness, c t (0,δ) is the tangential stiffness, r(x) represents a random parameter uniformly distributed between (0,1], β1 and β2 represent c nw and c tw Uniformity level index; η n is the displacement vector η along The weight, η t1 is the displacement vector η along The weight, η t2 is the displacement vector η along The weight, represents the unit vector along the initial direction of the bond, and represents the unit vector perpendicular to the initial direction of the bond, and and Perpendicular to each other.

9. A constitutive model implementation method, used for the rock mass single-bond dual-parameter dynamic damage constitutive model according to any one of claims 1 to 8, characterized in that: S101: Parameter initialization, determining model geometry parameters, and defining material parameters; S102: Discretize the solution domain and generate material point coordinates; S103: Determine the points within the material point domain, and initialize the keys within the material point domain; S104: Assign values ​​of the axial stiffness and tangential stiffness of each material point according to Weibull distribution, and take the average value to calculate the axial stiffness and normal stiffness of each key; S105: Determine the time step; S106: Apply velocity and displacement boundary conditions; S107: Calculate the initial orientation angle, axial deformation and normal deformation of the key; S108: Calculate volume correction factor; S109: Calculate the tensile and compressive deformation rates of the key S110: Determine the critical elongation of the key; S111: Calculate the constitutive forces of material points; S112: Calculate the elongation of the bond and determine whether the elongation is greater than the critical elongation or less than the critical compression rate. If "yes", the bond between the material points is broken and the constitutive force between the material points is 0. If "no", calculate the constitutive force according to step 111; S113: Use forward difference technique to solve displacement and velocity; S114: Output the calculation results and draw a damage cloud map.

Images