Modeling and calculating method of phanerocrystalline rock with bedding structure based on neper-dem

Abstract

The application discloses a Neper-DEM-based modeling and calculation method of phanerocrystalline rocks with bedding structure, relates to the numerical calculation field of tunnel engineering technology, and solves the problem that the prior art lacks a simulation and calculation method of phanerocrystalline rocks such as gneiss and granulitic gneiss with transverse isotropy; the application comprises the following steps: constructing a polycrystal structure with crystal characteristics in Neper, and then generating the polycrystal structure into a geometric parameter text of a crystal grain; constructing a discrete element rock sample model in DEM; performing numerical simulation test on the discrete element rock sample model to obtain a stress-strain curve and a sample failure mode, wherein the numerical simulation is realized by adopting a Brazilian split test and a triaxial test; calibrating discrete element sample calculation parameters according to the numerical simulation result, and verifying the effectiveness of the modeling and calculation method of the discrete element rock sample model; and the application realizes the construction and assignment of the crystal structure and the bedding structure at the same time in the modeling method, so that the rock structure is closer to the real rock structure.

Description

Technical Field

[0001] The present invention relates to the field of numerical calculation of tunnel engineering technology, and in particular to a modeling and calculation method of phaneritic rock with bedding structure based on Neper-DEM. Background Art

[0002] As a natural material, rocks develop varying degrees of anisotropy and heterogeneity, such as microcracks and defects, under long-term geological conditions. Rock anisotropy refers to the directional variations in the rock's physical and mechanical properties, primarily influenced by its internal structure, such as bedding, foliation, fractures, and mineral arrangement. For example, sedimentary rocks (such as shale and slate) develop a layered structure through deposition at different times and under different environments. Metamorphic rocks (such as diabase and slate) also develop anisotropic characteristics during metamorphism.

[0003] Directed pressure during metamorphism causes the mineral grains or layers in the rock to align in parallel, forming structural planes with directional characteristics. Due to this layered or folded structure, these rocks are considered transversely isotropic rocks, that is, these rocks have the same physical properties in two perpendicular directions, but exhibit different properties in a third direction (usually the direction of the bedding plane or major fracture). In addition to the anisotropy caused by the bedding structure, heterogeneous mineral composition also affects the mechanical properties of the rock. As a typical mineral-structured rock, granite has a complex microstructure, such as mineral content, shape, size and spatial distribution, as well as defects in the form of voids, microcracks, weak grain boundaries and fractures, which control the generation, propagation and aggregation of cracks and fractures.

[0004] Many rocks in nature exhibit both stratification and a distinct polycrystalline structure. These rocks are typically metamorphic rocks, though some igneous rocks may also exhibit similar characteristics. Gneiss, a metamorphic rock, exhibits a distinct schistose structure, where light and dark mineral components are banded and oriented. Mineral crystals in gneiss are typically large, with minerals such as quartz and feldspar forming distinct crystals. Biotite and hornblende in the dark bands can also have large grains. Schist also exhibits distinct foliated stratification, where the oriented arrangement of minerals creates a parallel structure. Mineral crystals are typically large, particularly mica, chlorite, and garnet, which can form distinct lamellae or scale-like crystals. Granite, a typical igneous rock, can also exhibit transversely isotropy. For example, if fractures are prevalent in a particular direction, the tensile and compressive strengths of the granite may differ perpendicular to the fractures compared to those along the fractures. Metamorphism, sedimentation, or intrusive bedding can also produce oriented mineral crystals, particularly in granitic gneiss. These phenomena can cause granite to exhibit transverse isotropy in microscopic or macroscopic structures. The crystal structure and bedding structure of rocks have a significant impact on rock mass stability, stress distribution, deformation characteristics, crack propagation, and failure modes. These rocks are both anisotropic and heterogeneous, and may exhibit larger deformations perpendicular to bedding or cracks, while smaller deformations parallel to the bedding or cracks. Due to the presence of non-uniformly distributed minerals, large local deformations may also occur. In terms of failure modes, cracks are often more likely to propagate along bedding, mineral orientation, and weaker mineral structure interfaces, which may lead to a non-uniform distribution of cracks and faults, affecting the stability of the rock mass.

[0005] Discrete element software can more clearly and directly observe the deformation and failure characteristics of such rocks and study their crack propagation mechanisms. The discrete element method can simulate the generation, expansion, and merging of cracks at the particle scale, which is more advantageous than traditional continuous medium methods (such as the finite element method) and is particularly suitable for studying the failure mechanism of rock materials. The discrete element method can also carefully consider the interactions between particles, including adhesion and friction, making the simulation results closer to the actual mechanical behavior of the material, and can accurately simulate the local failure process of the material, analyze the initiation and propagation path of micro cracks in the material, and reveal the material's failure mechanism.

[0006] Existing specimen modeling and calculation technologies have the following deficiencies: (1) They can only simulate and calculate rocks with a single layered structure or a crystal structure, but there is little research on the interaction mechanism and mechanical properties of the two; (2) The research on rock specimens focuses on the two-dimensional plane, ignoring the three-dimensional ductility and joint crack tendency; (3) The discrete element calibration method is mainly based on trial and error, which has a large amount of calculation and relatively low accuracy; (4) The grain size is relatively uniform and the distribution pattern is quite different from the distribution of mineral crystals in rocks in nature. (5) There is little analysis of the crack generation and propagation process and the proportion of each type of crack, which cannot reflect the development of the rock failure process. Summary of the Invention

[0007] In order to solve the problems existing in the above-mentioned prior art, the present invention provides a modeling and calculation method for phaneritic rocks with bedding structure based on Neper-DEM, which solves the problem that the prior art lacks simulation and calculation methods for phaneritic rocks with transverse isotropy such as gneiss and granite gneiss.

[0008] Neper-DEM-based modeling and calculation methods for phaneritic rocks with bedding structures include:

[0009] Construct a polycrystalline structure with crystal characteristics in Neper, and then generate the polycrystalline structure into the geometric parameter text of the grains;

[0010] A discrete element rock specimen model is constructed in the DEM, including constructing the discrete element specimen and importing geometric parameter text as geometric information, grouping the discrete element specimens according to the geometric information, constructing multiple mutually parallel bedding planes and inter-group joint planes whose directions are consistent with the grain intersection plane in the discrete element specimens, and then assigning contact models and parameters of three types of contact modes to the discrete element specimens, including intra-group contact of discrete element specimens and indirect contact between discrete element specimens and joint planes;

[0011] Numerical simulation tests are performed on discrete element rock sample models to obtain stress-strain curves and sample failure modes. The numerical simulations are implemented using Brazilian splitting tests and triaxial tests.

[0012] The calculation parameters of the discrete element sample are calibrated according to the numerical simulation results to verify the effectiveness of the modeling and calculation method of the discrete element rock sample model.

[0013] Furthermore, the polycrystalline structure creation process includes:

[0014] The lens is constructed using the 3D Thiessen polygon tessellation method based on the crystal characteristics, namely the number of grains, grain shape, and grain distribution. The number of grains is specified in the Neper command, the grain shape is reflected by the normal distribution of equivalent diameter and sphericity, and the grain distribution is controlled by random position and seed generation weight. Under the influence of different seed weights, the boundaries between grains become non-equidistant between seeds. The greater the seed weight, the larger the associated grain. After the Thiessen polygon tessellation is completed, the various parameters of the lens are output as an STL file.

[0015] Furthermore, grouping discrete element samples according to geometric information includes:

[0016] Importing geometric information, including the normal vectors of each face of the convex grain and the node coordinates of each face, wherein multiple faces together constitute a grain entity; using the DEM to read the geometric information, generating the same node domain surface in the DEM, constructing a single Geometry based on the convex grain composed of multiple faces and naming it according to the import order;

[0017] Group the existing discrete element samples according to the named Geometry. When multiple particles fall within the Geometry polyhedron, they are identified as a separate group (Group), and their naming is consistent with the Geometry.

[0018] Repeat the above steps until all the grains constituting the polycrystal are imported and grouped.

[0019] Furthermore, the step of constructing a plurality of mutually parallel bedding planes and inter-group joint planes whose directions are consistent with the grain intersection plane includes:

[0020] A certain number of parallel joint planes are added to the polycrystalline discrete element sample, and the bedding direction is controlled by the dip and inclination angle to be consistent with the original sample and the spacing between each bedding is kept consistent;

[0021] Construct inter-group joint surfaces according to the directions of the imported grain faces, and determine whether the two contacting grains are in the same group. If they are in the same group, no processing is done; if they are not in the same group, the group names are read separately, and the group names correspond to the same-named Geometry.

[0022] Find the two corresponding Geometries based on the name, determine the common face of the Geometries, obtain the direction of the face and assign the same direction to the joint surfaces of the two contacting particles.

[0023] Furthermore, the contact models and parameters for assigning the three types of contact modes include:

[0024] Assign contact models and parameters to the discrete element model based on the established bedding planes and inter-group joint planes;

[0025] The bedding plane is assigned a smooth joint model, which has the same direction as the bedding plane and a small contact parameter to simulate a weak structural plane.

[0026] The inter-group joint surface is also assigned a smooth joint model to simulate the interface contact between crystals, and the contact parameter value is larger than that of the layer;

[0027] The particles that are not covered by the inter-group joint planes and bedding planes are assigned a flat joint model.

[0028] Furthermore, the experimental process of the numerical simulation experiment is as follows:

[0029] Create wall units at the top and bottom of the discrete element rock specimen model, use the wall units to simulate the indenter of the Brazilian splitting test, set the modulus and width of the wall units, and fix the wall loading speed.

[0030] When the wall begins to be loaded, the displacement of the wall compression, the force applied during the compression process, and the crack distribution and number of the polycrystalline model are recorded. The data and images are automatically exported through the DEM command to obtain the stress-strain curve, elastic modulus, failure mode and failure strength.

[0031] Furthermore, obtaining the stress-strain curve includes:

[0032] The axial force, axial strain, and displacement of the wall are monitored in the discrete element simulation of the Brazilian splitting test. The maximum peak force is continuously updated during the monitoring process. If the newly generated force is less than the peak force threshold, the specimen is considered to have failed.

[0033] The calculation is now completed, and the axial force and axial strain of the wall are output. The axial force is converted into stress using the formula, and the stress-strain curve is drawn. The stress σ t The calculation formula is as follows:

[0034]

[0035] where σ t is the Brazilian splitting stress, P is the wall axial force, D is the specimen diameter, and t is the specimen thickness.

[0036] Furthermore, obtaining the stress-strain curve further includes: monitoring the axial force, axial strain, and displacement of the wall in a discrete element simulation of the triaxial test, continuously updating the maximum peak force during the monitoring process, and determining that the specimen is damaged when the current force is less than a threshold value of the peak force;

[0037] The calculation is now completed, and the axial force and axial strain of the wall are output. The axial force is converted into stress using the formula, and the stress-strain curve is drawn. The stress σ c The calculation formula is as follows:

[0038]

[0039] where σ c is the axial stress, P1 is the bottom wall pressure, P2 is the top wall pressure, and D is the diameter of the cylindrical specimen.

[0040] Furthermore, the calibrating of discrete element sample calculation parameters based on the numerical simulation results includes:

[0041] The deformation parameters of the discrete element model are calibrated based on the triaxial test results. When the bedding direction is perpendicular to the loading direction, the bedding plane has the least effect on the elastic modulus of the specimen. At this time, the stiffness parameters of the smooth joint model and the flat joint model between grains are calibrated.

[0042] When the bedding direction coincides with the loading direction, the bedding plane stiffness has the greatest influence, and the stiffness parameters of the smooth joint model of the bedding plane are calibrated;

[0043] The strength parameters of the discrete element model are calibrated based on the results of the Brazilian splitting test. When the bedding plane direction is perpendicular to the loading direction, the strength parameters of the smooth joint model and the flat joint model are calibrated.

[0044] When the bedding direction coincides with the loading direction, the bedding plane stiffness has the greatest influence, and the strength parameters of the smooth joint model of the bedding plane are calibrated.

[0045] Furthermore, the process of obtaining the crack distribution of the polycrystalline model includes:

[0046] According to the different characteristics and failure modes of the three types of contact modes, the crack types are divided into 6 categories, namely, flat joint model tensile crack, flat joint shear crack, intergranular smooth joint model tensile crack, intergranular smooth joint model shear crack, bedding plane smooth joint model tensile crack, bedding plane smooth joint model shear crack;

[0047] During the test loading process, the value of the increase in the number of cracks as the strain increases is recorded, and a crack change picture is drawn to determine which type of damage accounts for the largest proportion in the test.

[0048] The beneficial effects of the present invention include:

[0049] The present invention improves the existing numerical simulation modeling methods for rocks such as gneiss, granite gneiss and schist, while taking into account the influence of the heterogeneity of the crystal structure and the transverse isotropy of the bedding structure. In the modeling method, the construction and assignment of the crystal structure and the bedding structure are simultaneously realized, making it closer to the real rock structure.

[0050] The present invention improves the polycrystal generation method in the existing discrete element software, introduces Neper software to construct a polycrystal model with three-dimensional Thiessen polygon mosaic, generates a more complex numerical model that conforms to the real crystal distribution, and realizes the import and grouping of the crystal model generated in Neper in the DEM software, effectively solving the problem of difficult crystal structure generation and single and uniform crystals in the discrete element model.

[0051] The present invention improves the problem of parameter calibration in existing discrete element software, and calibrates different parameters by considering the influence of bedding plane on rock strength respectively, thus solving the problem of time-consuming, low-efficiency and inaccurate calibration of the "trial and error method".

[0052] The present invention improves the problems of crack generation and statistics, accurately groups the cracks in the discrete element software, effectively distinguishes the tensile cracks and shear cracks generated in the crystal, crystal interface and bedding plane, and monitors the development and proportion of the number of each crack, thereby better realizing the analysis and research of rock failure modes. BRIEF DESCRIPTION OF THE DRAWINGS

[0053] Figure 1 This is a flow chart of the modeling and calculation method of phaneritic rocks with bedding structure based on Neper-DEM involved in the embodiments of this application.

[0054] Figure 2 Create grains for the Neper software Thiessen polygon tessellation method involved in this application.

[0055] Figure 3 This is the initial sample model constructed by the DEM software involved in this application.

[0056] Figure 4 The polycrystalline model of the polycrystalline structure in the Neper software involved in this application is imported into the DEM software.

[0057] Figure 5 A model with completed contact assignment for the discrete element model of the layered polycrystal involved in this application.

[0058] Figure 6 This is a flow chart of discrete element layered polycrystal parameter calibration involved in this application.

[0059] Figure 7 (a) is a comparison chart of the discrete element calibration strength and the test results when the inclination is 30° involved in this application.

[0060] Figure 7 (b) is a comparison chart of the discrete element calibration strength and the test results when the inclination is 60° involved in this application.

[0061] Figure 8 This is a comparison chart of the discrete element calibration failure mode and test results involved in this application.

[0062] Figure 9 (a) is a graph showing the variation of the failure strength with the inclination and the inclination angle in the sensitivity analysis involved in this application.

[0063] Figure 9 (b) is a diagram showing the change of failure mode with inclination and dip angle in the sensitivity analysis involved in this application.

[0064] Figure 9 (c)-(i) are diagrams showing the variation of crack proportion with inclination and dip angle in the sensitivity analysis involved in this application. DETAILED DESCRIPTION

[0065] In order to make the purpose, technical solutions and advantages of the embodiments of the present application clearer, the technical solutions in the embodiments of the present application will be clearly and completely described below in conjunction with the drawings in the embodiments of the present application. Obviously, the described embodiments are only part of the embodiments of the present application, not all of the embodiments. Therefore, the following detailed description of the embodiments of the present application provided in the drawings is not intended to limit the scope of the application for protection, but merely represents the selected embodiments of the present application. Based on the embodiments of the present application, all other embodiments obtained by those skilled in the art without making creative work are within the scope of protection of this application.

[0066] The following is combined with Figure 1-9 The specific embodiments of the present invention are described in detail;

[0067] Modeling and calculation methods of crystalline rocks with bedding structure based on Neper-DEM, such as Figure 1 Shown, including:

[0068] Construct a polycrystalline structure with crystal characteristics in Neper, and then generate the polycrystalline structure into the geometric parameter text of the grains;

[0069] A discrete element rock specimen model is constructed in the DEM, including constructing the discrete element specimen and importing geometric parameter text as geometric information, grouping the discrete element specimens according to the geometric information, constructing multiple mutually parallel bedding planes and inter-group joint planes whose directions are consistent with the grain intersection plane in the discrete element specimens, and then assigning contact models and parameters of three types of contact modes to the discrete element specimens, including intra-group contact of discrete element specimens and indirect contact between discrete element specimens and joint planes;

[0070] Numerical simulation is performed on a discrete element rock sample model to obtain stress-strain curves and sample failure modes, wherein the numerical simulation is performed using a Brazilian splitting test and a triaxial test;

[0071] The calculation parameters of the discrete element sample were calibrated according to the numerical simulation results to verify the effectiveness of the modeling and calculation method of the discrete element rock sample model, and this method was used to carry out sensitivity analysis.

[0072] In another embodiment, the polycrystalline structure creation process includes:

[0073] The lens is constructed using the 3D Thiessen polygon tessellation method based on the crystal characteristics, namely the number of grains, grain shape, and grain distribution. The number of grains is specified in the Neper command, the grain shape is reflected by the normal distribution of equivalent diameter and sphericity, and the grain distribution is controlled by random position and seed generation weight. Under the influence of different seed weights, the boundaries between grains become non-equidistant between seeds. The greater the seed weight, the larger the associated grain. After the Thiessen polygon tessellation is completed, the various parameters of the lens are output as an STL file.

[0074] Specifically in this embodiment, the size and crystal characteristics of the model are determined, and a Thiessen polygon mesh is constructed in Neper to create a polycrystalline structure with the required crystal characteristics. Neper generates specific graphics such as Figure 2 The constructed polycrystalline structure is generated as a text of the geometric parameters of the grains; a discrete element rock sample model is constructed in the DEM software, and the initial parameters of the discrete element model are adjusted to achieve the initial balance. The initial discrete element model is as follows Figure 3 As shown;

[0075] The crystal characteristics include three parameters: generation quantity, equivalent diameter and sphericity; the discrete element initial parameters include particle effective modulus, stiffness ratio, friction coefficient, particle size ratio, minimum particle size and particle density.

[0076] Some parameters of the Neper modeling and DEM modeling are as follows: (The parameters in the table are actual values ​​in the selected embodiment)

[0077] Table 1 Parameter definition table

[0078]

[0079]

[0080] In another embodiment, grouping discrete element samples according to geometric information includes:

[0081] Importing geometric information, including the normal vectors of each face of the convex grain and the node coordinates of each face, wherein multiple faces together constitute a grain entity; using the DEM to read the geometric information, generating the same node domain surface in the DEM, constructing a single Geometry based on the convex grain composed of multiple faces and naming it according to the import order;

[0082] Existing discrete element samples are grouped according to the named Geometry. When multiple particles fall within the range of the Geometry polyhedron, they are identified as a separate group (Group), and its naming is consistent with the Geometry.

[0083] Repeat the above steps until all the grains constituting the polycrystal are imported and grouped.

[0084] Specifically in this embodiment, the STL file exported from Neper contains the normal vectors of each face of the convex grain and the node coordinates of each face, and multiple faces together constitute a grain entity, with a total of 2442 grains; the DEM software can read the command of the STL file, regenerate the same nodes and faces in the DEM, and construct the convex grains composed of multiple faces into a separate Geometry and name it according to the import order; the existing discrete element samples are grouped according to the named Geometry, and when multiple particles fall within the range of the Geometry polyhedron, they are identified as a separate group (Group), and their naming is consistent with the Geometry. Repeat the above steps until all 2442 grains that constitute the polycrystal are imported and divided into 2442 groups. The imported polycrystal discrete elements are as follows: Figure 4 shown.

[0085] In another embodiment, constructing a plurality of mutually parallel bedding planes and inter-group joint planes whose directions are consistent with the grain intersection plane includes:

[0086] A certain number of parallel joint planes are added to the polycrystalline discrete element sample, and the bedding direction is controlled by the dip and inclination angle to be consistent with the original sample and the spacing between each bedding is kept consistent;

[0087] Construct inter-group joint surfaces according to the directions of the imported grain faces, and determine whether the two contacting grains are in the same group. If they are in the same group, no processing is done; if they are not in the same group, the group names are read separately, and the group names correspond to the same-named Geometry.

[0088] Find the two corresponding Geometries based on the name, determine the common face of the Geometries, obtain the direction of the face and assign the same direction to the joint surfaces of the two contacting particles.

[0089] Specifically in this embodiment, five parallel bedding planes are added to the polycrystalline discrete element sample. The bedding direction is controlled by the inclination and dip angle to be consistent with the original sample and the spacing between each bedding is consistent, and the spacing is controlled at 0.75 mm. The inter-group joint surface is constructed according to the direction of the imported grain surface, and it is determined whether the two contacting particles are the same group. If they are the same group, they are not processed. If they are not the same group, the group names are read separately. It can be seen from the above that the group name corresponds to the Geometry with the same name. The corresponding two Geometries are found according to the name, the common face of the Geometry is determined, the direction of the face is obtained, and the joint faces of the two contacting particles are given the same direction.

[0090] In another embodiment, the contact models and parameters for assigning three types of contact modes include:

[0091] Assign contact models and parameters to the discrete element model based on the established bedding planes and inter-group joint planes;

[0092] The bedding plane, i.e. the joint surface contact, is assigned a smooth joint model, which has the same direction as the bedding plane and a small contact parameter to simulate a weak structural surface.

[0093] The inter-group joint surface, i.e., the contact between discrete element specimens, is also given a smooth joint model to simulate the interface contact between crystals, and the contact parameter value is larger than that of the layer.

[0094] A flat joint model is assigned to particles that are not covered by inter-group joint planes and bedding planes, i.e., particles that are not in contact within discrete element sample groups.

[0095] Specifically, in this embodiment, the discrete element model after the three types of contact assignment are completed is as follows: Figure 5 shown.

[0096] In another embodiment, the experimental process of the numerical simulation test is:

[0097] Create wall units at the top and bottom of the discrete element rock specimen model, use the wall units to simulate the indenter of the Brazilian splitting test, set the modulus and width of the wall units, and fix the wall loading speed.

[0098] When the wall begins to be loaded, the displacement of the wall compression, the force applied during the compression process, and the crack distribution and number of the polycrystalline model are recorded. The data and images are automatically exported through the DEM command to obtain the stress-strain curve, elastic modulus, failure mode and failure strength.

[0099] In another embodiment, obtaining the stress-strain curve includes:

[0100] The axial force, axial strain, and displacement of the wall are monitored in the discrete element simulation of the Brazilian splitting test. The maximum peak force is continuously updated during the monitoring process. If the newly generated force is less than the peak force threshold, the specimen is considered to have failed.

[0101] The calculation is now completed, and the axial force and axial strain of the wall are output. The axial force is converted into stress using the formula, and the stress-strain curve is drawn. The stress σ t The calculation formula is as follows:

[0102]

[0103] where σ t is the Brazilian splitting stress, P is the wall axial force, D is the specimen diameter, and t is the specimen thickness.

[0104] The axial force, axial strain, and displacement of the wall are monitored in the discrete element simulation of the triaxial test. The maximum peak force is continuously updated during the monitoring process. When the current force is less than the peak force threshold, the specimen is judged to have failed.

[0105] The calculation is now completed, and the axial force and axial strain of the wall are output. The axial force is converted into stress using the formula, and the stress-strain curve is drawn. The stress σ c The calculation formula is as follows:

[0106]

[0107] where σ c is the axial stress, P1 is the bottom wall pressure, P2 is the top wall pressure, and D is the diameter of the cylindrical specimen.

[0108] In this embodiment, the threshold is 80%.

[0109] In another embodiment, calibrating discrete element sample calculation parameters based on numerical simulation results includes:

[0110] The deformation parameters of the discrete element model are calibrated based on the triaxial test results. When the bedding direction is perpendicular to the loading direction, the bedding plane has the least effect on the elastic modulus of the specimen. At this time, the stiffness parameters of the smooth joint model and the flat joint model between grains are calibrated.

[0111] When the bedding direction coincides with the loading direction, the bedding plane stiffness has the greatest influence, and the stiffness parameters of the smooth joint model of the bedding plane are calibrated;

[0112] The strength parameters of the discrete element model are calibrated based on the results of the Brazilian splitting test. When the bedding plane direction is perpendicular to the loading direction, the strength parameters of the smooth joint model and the flat joint model are calibrated.

[0113] When the bedding direction coincides with the loading direction, the bedding plane stiffness has the greatest influence, and the strength parameters of the smooth joint model of the bedding plane are calibrated.

[0114] Specifically in this embodiment, the main calibrated parameters are the flat joint model parameters of the grains inside the rock and the smooth joint model parameters between groups (between grains) and weak bedding planes; the deformation parameters of the discrete element model are calibrated by the triaxial test results. When the bedding direction is perpendicular to the loading direction, the bedding plane has the least effect on the elastic modulus of the specimen. At this time, the stiffness parameters of the smooth joint model and the flat joint model between the grains are calibrated; when the bedding direction coincides with the loading direction, the bedding plane stiffness has the greatest effect. The stiffness parameters of the smooth joint model of the bedding plane are calibrated; the strength parameters of the discrete element model are calibrated according to the Brazilian splitting test results. When the bedding plane direction is perpendicular to the loading direction, the strength parameters of the smooth joint model and the flat joint model are calibrated; when the bedding direction coincides with the loading direction, the bedding plane stiffness has the greatest effect. The strength parameters of the smooth joint model of the bedding plane are calibrated. The calibration process is as follows: Figure 6 As shown, the deformation parameters include the effective modulus of the flat joint model Stiffness ratio of flat joint model Intercrystalline

[0115] Normal stiffness of smooth joint model Tangential stiffness of intergranular smooth joint model Normal stiffness of smooth joint model of bedding plane Shear stiffness of smooth joint model of bedding plane The strength parameters include the tensile strength of the flat joint model Flat joint model cohesion Friction angle of flat joint model Tensile strength of intergranular smooth joint model Cohesion of intergranular smooth joint model Friction angle of intergranular smooth joint model Tensile strength of smooth joint model on bedding plane Cohesion of smooth joint model of bedding plane Friction angle of smooth joint model of bedding plane The calibration results are shown in Table 2, and compared with the test intensity Figure 7 As shown, the failure mode comparison is as follows Figure 8 shown.

[0116] Table 2 Calibration parameters

[0117]

[0118] In another embodiment, the process of obtaining the crack distribution of the polycrystalline model includes:

[0119] According to the different characteristics and failure modes of the three types of contact modes, the crack types are divided into 6 categories, namely, flat joint model tensile crack, flat joint shear crack, intergranular smooth joint model tensile crack, intergranular smooth joint model shear crack, bedding plane smooth joint model tensile crack, bedding plane smooth joint model shear crack;

[0120] During the test loading process, the value of the increase in the number of cracks as the strain increases is recorded, and a crack change picture is drawn to determine which type of damage accounts for the largest proportion in the test.

[0121] In another embodiment, a sensitivity analysis is performed on the model to analyze the effect of the change in the bedding direction on the strength, failure mode and crack propagation of the polycrystalline model. The tilt direction changes from 15° to 90°, and the inclination angle changes from 0° to 90°. The strength change results are as follows: Figure 9 shown. Figure 9 a represents the failure strength of the specimens at various inclinations and angles under the Brazilian splitting test. Figure 9 b Distribution of the six types of cracks after the specimen loading is completed: flat joint model tensile crack, flat joint shear crack, intergranular smooth joint model tensile crack, intergranular smooth joint model shear crack, bedding plane smooth joint model tensile crack, bedding plane smooth joint model shear crack. Figure 9 c- Figure 9 i is the graph showing the crack ratio changing with the inclination and inclination angle in the sensitivity analysis. Specifically, Figure 9 c- Figure 9 e are the specific crack distribution diagrams at inclinations of 15°30°45°60°75°90°, Figure 9 i represents Figure 9 c- Figure 9 The diagram of the six types of cracks involved in e, GBSC refers to intra-group shear cracks, GBTC refers to intra-group tensile cracks, BPSC refers to joint shear cracks, BPTC refers to intra-group tensile cracks, IGSC refers to inter-group shear cracks, and IGTC refers to inter-group tensile cracks.

[0122] The above-described embodiments merely represent specific implementation methods of the present application. While the descriptions are relatively specific and detailed, they should not be construed as limiting the scope of protection of the present application. It should be noted that a person skilled in the art would be able to make numerous variations and improvements without departing from the technical concept of the present application, and all such variations and improvements fall within the scope of protection of the present application.

Claims

1. Modeling and calculation method of phaneritic rocks with bedding structure based on Neper-DEM, characterized by: include: Construct a polycrystalline structure with crystal characteristics in Neper, and then generate the polycrystalline structure into the geometric parameter text of the grains; A discrete element rock specimen model is constructed in the DEM, including constructing the discrete element specimen and importing geometric parameter text as geometric information, grouping the discrete element specimens according to the geometric information, constructing multiple mutually parallel bedding planes and inter-group joint planes whose directions are consistent with the grain intersection plane in the discrete element specimens, and then assigning contact models and parameters of three types of contact modes to the discrete element specimens, including intra-group contact of discrete element specimens and indirect contact between discrete element specimens and joint planes; Numerical simulation tests are performed on discrete element rock sample models to obtain stress-strain curves and sample failure modes. The numerical simulations are implemented using Brazilian splitting tests and triaxial tests. The calculation parameters of the discrete element sample were calibrated based on the numerical simulation results to verify the effectiveness of the modeling and calculation method of the discrete element rock sample model; The polycrystalline structure creation process includes: The lens is constructed using a 3D Thiessen polygon tessellation method based on crystal characteristics, namely, the number of grains, grain shape, and grain distribution. The number of grains is specified in the Neper command, the grain shape is represented by the normal distribution of equivalent diameter and sphericity, and the grain distribution is controlled by random positions and seed generation weights. Under the influence of different seed weights, the boundaries between grains become non-equidistant between seeds. The larger the seed weight, the larger the associated grain. After the Thiessen polygon tessellation is completed, the various parameters of the lens are exported as an STL file. The contact models and parameters assigned to the three types of contact modes include: Assign contact models and parameters to the discrete element model based on the established bedding planes and inter-group joint planes; The bedding plane is assigned a smooth joint model, which has the same direction as the bedding plane and a small contact parameter to simulate a weak structural plane. The inter-group joint surface is also assigned a smooth joint model to simulate the interface contact between crystals, and the contact parameter value is larger than that of the layer; The particles that are not covered by the intergroup joint planes and bedding planes are assigned a flat joint model; The calibration of discrete element sample calculation parameters according to the numerical simulation results includes: The deformation parameters of the discrete element model are calibrated based on the triaxial test results. When the bedding direction is perpendicular to the loading direction, the bedding plane has the least effect on the elastic modulus of the specimen. At this time, the stiffness parameters of the smooth joint model and the flat joint model between grains are calibrated. When the bedding direction coincides with the loading direction, the bedding plane stiffness has the greatest influence, and the stiffness parameters of the smooth joint model of the bedding plane are calibrated; The strength parameters of the discrete element model are calibrated based on the results of the Brazilian splitting test. When the bedding plane direction is perpendicular to the loading direction, the strength parameters of the smooth joint model and the flat joint model are calibrated. When the bedding direction coincides with the loading direction, the bedding plane stiffness has the greatest influence, and the strength parameters of the smooth joint model of the bedding plane are calibrated.

2. The Neper-DEM-based phaneritic rock modeling and calculation method with bedding structure according to claim 1, characterized in that: The grouping of discrete element samples according to geometric information includes: Importing geometric information, including the normal vectors of each face of the convex grain and the node coordinates of each face, wherein multiple faces together constitute a grain entity; using the DEM to read the geometric information, generating the same node domain surface in the DEM, constructing a single Geometry based on the convex grain composed of multiple faces and naming it according to the import order; Group the existing discrete element samples according to the named Geometry. When multiple particles fall within the Geometry polyhedron, they are identified as a separate group, and their naming is consistent with the Geometry. Repeat the above steps until all the grains constituting the polycrystal are imported and grouped.

3. The Neper-DEM-based phaneritic rock modeling and calculation method with bedding structure according to claim 1, characterized in that: The method of constructing a plurality of mutually parallel bedding planes and inter-group joint planes whose directions are consistent with the grain intersection planes comprises: A certain number of parallel joint planes are added to the polycrystalline discrete element sample, and the bedding direction is controlled by the dip and inclination angle to be consistent with the original sample and the spacing between each bedding is kept consistent; Construct inter-group joint surfaces according to the directions of the imported grain faces, and determine whether the two contacting grains are in the same group. If they are in the same group, no processing is done; if they are not in the same group, the group names are read separately, and the group names correspond to the same-named Geometry. Find the two corresponding Geometries based on the name, determine the common face of the Geometries, obtain the direction of the face and assign the same direction to the joint surfaces of the two contacting particles.

4. The Neper-DEM-based phaneritic rock modeling and calculation method with bedding structure according to claim 1, characterized in that: The experimental process of the numerical simulation test is as follows: Create wall units at the top and bottom of the discrete element rock specimen model, use the wall units to simulate the indenter of the Brazilian splitting test, set the modulus and width of the wall units, and fix the wall loading speed. When the wall begins to be loaded, the displacement of the wall compression, the force applied during the compression process, and the crack distribution and number of the polycrystalline model are recorded. The data and images are automatically exported through the DEM command to obtain the stress-strain curve, elastic modulus, failure mode and failure strength.

5. The Neper-DEM-based phaneritic rock modeling and calculation method with bedding structure according to claim 1, characterized in that: The obtaining of the stress-strain curve comprises: The axial force, axial strain, and displacement of the wall are monitored in the discrete element simulation of the Brazilian splitting test. The maximum peak force is continuously updated during the monitoring process. If the newly generated force is less than the peak force threshold, the specimen is considered to have failed. The calculation is now completed, and the axial force and axial strain of the wall are output. The axial force is converted into stress using the formula, and the stress-strain curve is drawn. The calculation formula is as follows: ； in is the Brazilian splitting stress, is the wall axial force, is the sample diameter, is the sample thickness.

6. The Neper-DEM-based phaneritic rock modeling and calculation method with bedding structure according to claim 1, characterized in that: The obtaining of the stress-strain curve further comprises: The axial force, axial strain, and displacement of the wall are monitored in the discrete element simulation of the triaxial test. The maximum peak force is continuously updated during the monitoring process. When the current force is less than the peak force threshold, the specimen is judged to have failed. The calculation is now completed, and the axial force and axial strain of the wall are output. The axial force is converted into stress using the formula, and the stress-strain curve is drawn. The calculation formula is as follows: ； in is the axial stress, is the bottom wall pressure, is the top wall pressure, is the diameter of the cylindrical specimen.

7. The Neper-DEM-based phaneritic rock modeling and calculation method with bedding structure according to claim 4, characterized in that: The process of obtaining the crack distribution of the polycrystalline model includes: According to the different characteristics and failure modes of the three types of contact modes, the crack types are divided into 6 categories, namely, flat joint model tensile crack, flat joint shear crack, intergranular smooth joint model tensile crack, intergranular smooth joint model shear crack, bedding plane smooth joint model tensile crack, bedding plane smooth joint model shear crack; During the test loading process, the value of the increase in the number of cracks as the strain increases is recorded, and a crack change picture is drawn to determine which type of damage accounts for the largest proportion in the test.

Images