Three-dimensional discrete fracture network modeling and simulation method based on fracture statistical characteristics

Abstract

The invention discloses a three-dimensional discrete fracture network modeling and simulation method based on fracture statistical characteristics, and relates to the technical field of water conservancy and hydropower engineering rock soil. According to the method, geometric and distribution characteristics of rock mass fractures are collected on site, and a random statistical model containing parameter relevance is constructed; generating a three-dimensional DFN model based on a random statistical model, and performing intersection processing, invalid fracture elimination and boundary correction optimization; the three-dimensional DFN model is embedded into a numerical platform, a fracture-rock system is constructed, Mohr-Coulomb constitutive parameters are given, crustal stress inversion, excavation supporting and operation load simulation are completed in sequence, and mechanical response data are obtained; and finally analyzing surrounding rock failure characteristics and optimizing a support scheme. The problem that a traditional model is difficult to represent a discontinuous structure and local instability is solved, the prediction precision is improved, the method is adaptive to various complex fractured rock mass projects, and reliable support is provided for stability evaluation and support design.

Description

Technical Field

[0001] This invention relates to the field of geotechnical engineering technology in water conservancy and hydropower engineering construction, and in particular to a digital rock mass modeling technology for characterizing complex geological conditions. Specifically, it proposes a three-dimensional discrete fracture network modeling and simulation method based on the statistical characteristics of fractures. Background Technology

[0002] In the construction of water conservancy and hydropower projects, the long-term stability of key structures such as dam foundations, slopes, powerhouses, and long-distance water diversion tunnels is a fundamental prerequisite for ensuring their effectiveness. Water conservancy and hydropower construction areas are often located in mountainous regions with complex geological conditions, where faults, joints, and fissures are widely developed within the rock mass. The existence of these discontinuous structural planes causes the rock mass to exhibit significant heterogeneity and anisotropy in its mechanical behavior.

[0003] Traditional rock mass stability analysis often employs mature numerical software, such as FLAC3D, to establish an equivalent continuous medium model. The influence of fractures is typically assessed simply by adjusting rock mass mechanical parameters. While this method is relatively simple, it has inherent limitations. For instance, when dealing with local instability controlled by structural planes, it struggles to accurately reflect the shear slip and tensile cracking that occur along dominant fracture surfaces under stress, leading to significant discrepancies between predicted results and actual field monitoring. Summary of the Invention

[0004] The purpose of this invention is to provide a three-dimensional discrete fracture network modeling and simulation method based on the statistical characteristics of fractures, which aims to solve the problem that traditional equivalent continuous medium models are difficult to realistically represent the discontinuous structure and local instability mechanism of complex rock masses, thereby improving the accuracy and reliability of predicting the mechanical response of rock masses under excavation and operation loads.

[0005] To achieve the above objectives, the technical method of the present invention includes: A method for modeling and simulating three-dimensional discrete fracture networks based on the statistical properties of fractures includes the following steps: S1: Statistical characterization of the geometric parameters of the rock mass fracture system, including collecting the geometric parameters and distribution characteristics of rock mass fractures in the study area, constructing a multi-dimensional parameter system, and establishing a stochastic statistical model of rock mass fractures; the multi-dimensional parameter system includes statistical analysis of fracture geometric parameters, statistical analysis of fracture spatial distribution parameters, and parameter correlation analysis; the stochastic statistical model of rock mass fractures uses field measured data as input to quantify the stochastic characteristics and intrinsic correlation of key fracture parameters; The intrinsic correlation refers to the quantitative correspondence between different parameters of rock mass fractures, including: The relationships between geometric parameters, such as the positive correlation between crack length and opening (the distance between the two walls of the crack), indicate that longer cracks are usually subjected to more stress when they form and tend to have a larger opening; and the relationship between roughness and shear strength, where a rough crack surface has greater frictional resistance and usually higher shear strength. The relationship between spatial distribution and geometric parameters, such as the correspondence between fracture dip and density, is as follows: if the study area is dominated by tectonic stress in a certain direction, the fractures in that direction, i.e. the dominant dip, will not only be more numerous and denser, but their dip angles may also be more concentrated. The relationship between size and distribution pattern, such as the relationship between fracture radius or length and spacing, is that when a large fracture forms, it induces the formation of secondary small fractures in the surrounding rock mass. Around a large fracture, the distribution of small fractures may be more dense. S2: Construction of the three-dimensional discrete fracture network (DFN) model. Based on the rock mass fracture stochastic statistical model established in step S1, the computational domain is delineated and a three-dimensional model of the engineering entity is established as a spatial constraint for fracture generation and trimming. The on-site statistical parameters are transformed into a geometric model through a Monte Carlo stochastic simulation program. After fracture intersection logic processing, invalid fracture removal and boundary truncation effect correction, the three-dimensional DFN model is completed. S3: Numerical simulation of rock mechanics response based on DFN model. The three-dimensional DFN model generated in step S2 is embedded into the numerical calculation platform, mapped to joint or structural surface units, and a "fracture-rock block" discontinuous medium system is constructed. The surrounding rock is given Mohr-Coulomb constitutive parameters, and the initial geostress inversion, engineering excavation and support simulation, and operation period load simulation are performed in sequence to obtain rock mechanics response data. Assigning Mohr-Coulomb constitutive parameters to the surrounding rock includes assigning Mohr-Coulomb constitutive parameters to rock blocks and to joints or structural planes (fractures). Assigning Mohr-Coulomb constitutive parameters to rock blocks refers to assigning Mohr-Coulomb constitutive parameters to the rock blocks themselves, which are independent rock mass units formed by fractures—the solid parts surrounded by the fracture network in the 3D DFN model. The assigned parameters are the rock blocks' own mechanical strength indices, such as compressive strength, tensile strength, internal friction angle, and cohesion, to define the failure rules of the rock blocks. For example, when the pressure on the rock block exceeds its compressive strength, it will undergo compressive fracture; when the tensile force exceeds its tensile strength, it will undergo tensile cracking, ensuring that the mechanical behavior of the rock blocks in the simulation closely matches that of real rock mass. Assigning Mohr-Coulomb constitutive parameters to joints or structural planes (fractures) refers to assigning Mohr-Coulomb constitutive parameters to the joints / structural planes, which are virtual fracture units mapped from the DFN model, corresponding to fracture surfaces in real rock mass. The parameters assigned are the shear mechanical properties of the fracture surface, such as the friction coefficient, cohesion, and shear strength of the fracture surface, in order to define the sliding rules of the fracture. For example, when the shear force on the fracture surface exceeds the shear strength, the rock blocks on both sides will undergo shear slip along the fracture surface, thus restoring the core mechanism of "slippage instability along the dominant fracture surface" in real rock masses. S4: Verification and engineering application of simulation results. Based on the mechanical response data in step S3, analyze the failure morphology and distribution characteristics of the surrounding rock, optimize the support scheme in combination with actual engineering needs, and provide a basis for stability evaluation of complex fractured rock mass engineering.

[0006] Furthermore, the method for constructing the three-dimensional discrete fracture network (DFN) model in step S2 includes: Determine the computational domain and spatial constraints: Based on the actual needs of the project, define the scope of the computational domain, establish a three-dimensional model of the engineering entity, clarify the spatial boundary of crack generation, and ensure that subsequent cracks only develop within the rock mass, and that cracks passing through the engineering entity are trimmed and corrected according to the entity outline. Write a Monte Carlo random simulation program: set the maximum number of iterations for random sampling and the convergence condition, wherein the convergence condition is that the relative error between the simulated fracture volume density and the measured value is ≤5%; Parameter transformation and fracture generation: Based on the volume density parameter in the random statistical model of rock mass fractures, the total number of fractures generated is determined by Poisson distribution. Fisher distribution is used to randomly sample the fracture orientation and dip angle. The equivalent radius of the fracture is extracted based on log-normal distribution. The fracture is simplified into a disk element, and the centroid coordinates of the disk strictly follow a three-dimensional uniform distribution. The generation range is set within the space 20% outside the boundary of the computational domain. Crack network optimization: Automatically identify the intersection lines between disk elements, retain valid cracks that form a connected network and mark isolated cracks, and remove invalid cracks that are completely outside the computational domain or have no interaction with engineering entities; Final model calibration: Boundary truncation effect calibration is applied to the screened cracks to complete the construction of the three-dimensional discrete crack network (DFN) model.

[0007] Furthermore, the numerical simulation method for rock mass mechanical response based on the DFN model in step S3 includes: Model embedding and system construction: The three-dimensional DFN model generated in step S2 is embedded into numerical computing platforms such as FLAC3D, and the DFN model is mapped into joint or structural surface units to construct a discontinuous medium system of fractures and rock blocks, and to clarify the spatial boundaries and interaction relationships between rock blocks and fractures. Assign Mohr-Coulomb constitutive parameters: Assign corresponding Mohr-Coulomb constitutive parameters to the rock blocks and joints / structures in the system respectively. Assign mechanical strength indicators such as compressive strength, tensile strength, internal friction angle, and cohesion to the rock blocks, and shear mechanical indicators such as friction coefficient, cohesion, and shear strength to the joints / structures. Initial geostress inversion: The multivariate regression analysis method is adopted, and the hydraulic fracturing measurement point data is selected as the fitting target. The normal displacement of the model boundary is set as the constraint condition. The initial geostress field in the calculation domain is expressed as a linear superposition of the self-weight stress field and the tectonic stress field. The regression coefficients and correction constants are solved by the least squares method to ensure that the residual between the inverted stress field and the measured value meets the engineering accuracy. Engineering excavation and support simulation: Based on the actual construction process of the project, the excavation process of tunnels and other structures is simulated, and the application of support measures is simulated simultaneously to restore the impact of excavation unloading and support reinforcement on the mechanical state of the rock mass. Operational load simulation: Apply corresponding loads according to the load requirements of the project operation phase, including equivalent concentrated force and uniformly distributed load applied according to the two-lane standard; Mechanical response data acquisition: The numerical calculation platform outputs mechanical response data such as displacement field, stress field, and plastic zone distribution of the rock mass to complete the numerical simulation of the rock mass mechanical response.

[0008] Furthermore, in step S2, the number of cracks N generated is determined using a Poisson distribution, with the estimation relationship being N≈(V·P) 32 ) / E[A], where V is the volume of the generating domain, P 32 E[A] represents the volume density of each dominant fracture group and the background fracture, and E[A] represents the average area of ​​the fracture disk.

[0009] Furthermore, in step S2, a Fisher distribution is used to randomly sample the fracture orientation and dip angle. The probability density function of the Fisher distribution... for: in Fisher's constant, The angle between the crack normal vector and the average vector is given.

[0010] Furthermore, in step S2, the equivalent radius of the crack is extracted based on a log-normal distribution, and the crack is simplified into a disk element. The generation range of the disk centroid is set within a space 20% outside the boundary of the computational domain, and the centroid coordinates of the simplified disk element are... It follows a three-dimensional uniform distribution.

[0011] Furthermore, in step S3, the initial geostress inversion will retrieve the initial geostress field within the computational domain. This can be expressed as a linear superposition of the self-weight stress field and the tectonic stress field. In the formula Indicates self-weight stress. Indicates structural stress in the X direction. Indicates structural stress in the Y direction. It is a contribution coefficient that adjusts the influence weight of self-weight and X or Y-direction tectonic stress. It is a correction constant.

[0012] The present invention discloses an apparatus comprising: at least one processor; and a memory communicatively connected to the at least one processor; wherein the memory stores instructions executable by the at least one processor, the instructions being executed by the at least one processor to enable the at least one processor to perform a method for modeling and simulating three-dimensional discrete fracture networks based on fracture statistical properties.

[0013] The present invention also discloses a computer-readable storage medium storing a computer program, which, when executed by a processor, implements a method for modeling and simulating three-dimensional discrete fracture networks based on fracture statistical characteristics.

[0014] The beneficial effects of this invention include: 1. Improve the accuracy of rock mass mechanical response prediction: This invention constructs a random statistical model of rock mass fractures based on field measured data, quantifies the random characteristics and intrinsic correlation of key fracture parameters, and combines it with a three-dimensional DFN model to realistically characterize the discontinuous structure of rock mass. This solves the problem that traditional equivalent continuous medium models cannot reflect local instability mechanisms, making the prediction of rock mass displacement, stress, and plastic zone distribution under excavation and operation loads more closely resemble engineering reality, and significantly reducing the deviation between prediction results and field monitoring.

[0015] 2. Enhanced engineering adaptability and dynamic optimization capability of the model: The random statistical model of rock mass fracture in this invention takes the specific fracture system in the study area as the object, and the statistical parameters are all from the field system survey and measurement. It can accurately match the geological environment of specific structures such as dam foundation, slope, and water diversion tunnel. At the same time, it supports dynamic supplementation and optimization of parameters. As the measured data accumulates, the parameter distribution law can be continuously corrected, and the model accuracy can be continuously improved, avoiding the coarseness of generalized models.

[0016] 3. Optimize the engineering support design and safety evaluation system: This invention accurately reveals the failure morphology and distribution characteristics of the surrounding rock through numerical simulation, and clarifies the occurrence area and mechanism of different failure types such as shallow tensile failure and deep shear failure, providing a precise basis for the optimization of support schemes and realizing "support on demand"; compared with traditional methods, it can quantitatively separate the influence of specific loads on rock mass stability, providing highly reliable technical support for the stability evaluation of complex fractured rock mass engineering, and reducing engineering safety risks.

[0017] 4. Balancing Modeling Efficiency and Computational Reliability: In the construction of the 3D DFN model, this invention utilizes optimization techniques such as convergence condition control of the Monte Carlo stochastic simulation program, invalid crack removal, and boundary effect correction. While ensuring the statistical consistency of the model (the relative error between simulated crack density and measured values ​​is ≤5%), it also reasonably controls the computational scale, avoiding the problem of low computational efficiency caused by excessive data volume, thus achieving a balance between modeling accuracy and computational efficiency.

[0018] 5. Expanding Engineering Application Scenarios and Practicality: The method of this invention is not only applicable to water conservancy and hydropower projects, but can also be extended to engineering fields involving complex fractured rock masses such as mining and transportation tunnels. Through the integrated process of "on-site statistics - three-dimensional modeling - numerical simulation - support optimization", it can provide full-process technical guidance for engineering design, construction and operation without relying on complex test equipment, and has strong engineering practicality and operability. Attached Figure Description

[0019] Figure 1 This is a flowchart illustrating the specific implementation steps of the method of the present invention.

[0020] Figure 2 This is a statistical feature diagram of joint orientation used in the present invention for modeling discrete fracture network (DFN).

[0021] Figure 3 This is a surface view of the three-dimensional discrete fracture network (DFN) model used for numerical simulation in this invention.

[0022] Figure 4 This is a statistical feature diagram of crack parameters in the three-dimensional digital rock mass model generated by this invention.

[0023] Figure 4 (a) in the figure is a histogram of the distribution of crack inclination angle in the statistical sample of the present invention.

[0024] Figure 4 (b) in the figure is the crack tendency rose diagram of the statistical sample of the present invention.

[0025] Figure 4 (c) in the figure is the histogram of the distribution of crack radius in the statistical sample of the present invention.

[0026] Figure 5 This is a characteristic diagram of displacement distribution along the axial profile of the tailrace tunnel according to an embodiment of the present invention.

[0027] Figure 5 In the example, (a) represents the total displacement and vector diagram of the tailrace tunnel after excavation in this embodiment of the invention.

[0028] Figure 5 (b) in the figure shows the incremental displacement and vector diagram of the tailrace tunnel caused by the application of highway load in the embodiment of the present invention.

[0029] Figure 6 This is a diagram showing the distribution of the plastic zone in the axial cross-section of the tailrace tunnel according to an embodiment of the present invention. Detailed Implementation

[0030] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only some, not all, of the embodiments of this invention, and are not intended to limit the invention. All other embodiments obtained by those skilled in the art based on the embodiments of this invention without inventive effort are within the scope of protection of this invention.

[0031] This invention introduces a Discrete Fracture Network (DFN) model that can directly characterize the geometric and statistical distribution characteristics of fractures. Using this method, a three-dimensional DFN model is efficiently generated by measuring fracture distribution and statistical parameters in the field. This model is then applied to numerical simulation calculations, significantly improving the understanding of the response characteristics of complex fractured rock masses during engineering construction. This effectively solves the technical challenge of traditional continuous medium models accurately characterizing complex geological structures and local instability mechanisms.

[0032] This invention proposes a three-dimensional discrete fracture network (DFN) modeling and simulation method based on the statistical characteristics of fractures. This method aims to realistically reflect the discontinuous structural characteristics of complex rock masses and improve the reliability of predicting rock mass engineering responses. The overall technical process is as follows: Figure 1 As shown, the process mainly includes four steps: obtaining statistical features of cracks, three-dimensional DFN modeling, numerical simulation based on DFN, and result verification and engineering application.

[0033] Step S1: Statistical characterization of geometric parameters of the rock mass fracture system The fundamental work of this method is to systematically and reliably characterize the statistical properties of the rock fracture system in the study area and establish a stochastic statistical model of fractures that can reflect the actual engineering situation. The statistical model adopted is the stochastic statistical model of rock fractures. The stochastic statistical model of rock fractures is a mathematical model constructed based on the actual statistical characteristics of the rock fracture system in the study area and used to characterize the geometric morphology and spatial distribution of fractures. It is the basis for three-dimensional discrete fracture network modeling and simulation. By quantifying the stochastic characteristics and intrinsic correlation of key fracture parameters, the stochastic statistical model of rock fractures restores the true state of rock fractures under complex geological conditions, providing geological model support for subsequent engineering mechanical response analysis. The stochastic statistical model for rock mass fractures uses field-measured data as input. It systematically collects the geometric parameters and distribution characteristics of rock mass fractures in the study area, constructing a multi-dimensional parameter system. This system mainly includes the following elements: 1. Geometric parameter statistics of fractures: This covers fracture size parameters such as length, width, and trace length; morphological parameters such as dip direction, dip angle, aperture, and roughness. Through statistical analysis of a large number of measured samples, the probability distribution type of each parameter is determined, such as normal distribution, log-normal distribution, and power-law distribution, along with corresponding statistical parameters such as mean, variance, and standard deviation. The statistical analysis of these large numbers of measured samples yields the geometric attribute patterns of individual fractures. 2. Spatial distribution parameter statistics of fractures: This statistically analyzes the distribution patterns of fractures in the rock mass space, including fracture density (e.g., the number of fractures per unit volume or unit area), spacing (e.g., the average distance between adjacent fractures), and grouping characteristics (e.g., the existence of fracture groups with dominant development directions). Statistical analysis reveals the aggregation or dispersion patterns of fractures in three-dimensional space, reflecting the overall spatial layout characteristics of the fracture system. 3. Parameter Correlation Analysis: Considering the intrinsic correlations between different fracture parameters, such as the correlation between fracture length and aperture, and the correspondence between dip and density, irrelevant parameter interference is eliminated through statistical testing methods. A quantitative correlation model between parameters is established to avoid the problem of disconnect from engineering reality caused by independent modeling of single parameters. The stochastic statistical model of rock mass fractures is constructed based on field measured data. The parameter distribution law completely matches the actual rock mass conditions in the study area. At the same time, the random fluctuation characteristics of fracture parameters are characterized by probabilistic statistical methods, which avoids the limitations of deterministic models and ensures the objectivity and reliability of the model. It can reflect the uncertainty in the natural formation process of rock mass fractures. The stochastic statistical model of rock mass fractures takes the specific fracture system in the study area as the research object. The statistical parameters are all from the systematic survey and measurement of the area. It abandons the coarseness of generalized models and can accurately match the geological environment of specific structures such as dam foundations, slopes, and water diversion tunnels in water conservancy and hydropower projects, ensuring the model's ability to reflect engineering reality. The stochastic statistical model of rock mass fractures supports dynamic supplementation and optimization of parameters. As the accumulated field measured data increases, the parameter distribution law can be continuously corrected to improve the model accuracy. At the same time, the model can be directly coupled with the three-dimensional discrete fracture network (DFN) model to transform statistical laws into a visualized three-dimensional fracture network, providing direct input for subsequent numerical simulation calculations.

[0034] Taking a water conveyance tunnel of a large hydropower station as an example, during the tunnel excavation process, all exposed fractures in the granodiorite, schist, and slate strata were systematically recorded on-site, with a total sample size of [number missing]. The main contents of the catalog include the attitude of each fracture, such as strike, dip, and dip angle, as well as basic geometric characteristics such as length, width, and trace length. Based on this, a statistical analysis of the fracture attitude is performed using stereographic projection, such as... Figure 2 As shown in the isodensity map, two distinct high-density concentration zones can be identified, indicating that the rock mass in the study area is mainly controlled by two sets of tectonic fractures. The frequency distribution statistics of the dip angle are as follows: Figure 2 As shown in the bottom right corner, Figure 2 The frequency distribution of the dip angle in the lower right corner indicates that the dip angle is... The cracks within the range account for the total number of The rock mass is dominated by steeply dipping fractures. This result is consistent with the concentration of poles towards the periphery in the stereographic diagram, indicating that the rock mass in the study area is dominated by high-angle structural planes. The statistical characteristics of the above-mentioned attitude, dip angle, density, and size distribution will serve as input parameters and constraints for subsequent DFN modeling.

[0035] The basis for drawing the pole isodensity map is to obtain a sufficient number of measured fracture attitude data covering the study area to ensure the reliability of the statistical results. In this embodiment of the invention, the sample size reaches 11,348, which meets the statistical requirements for the number of samples. The steps to obtain the pole isodensity map are as follows: Full-section recording of fractures exposed in the rock mass of the study area, including structural fractures and joints, is performed. The strike of each fracture is recorded as the azimuth angle of the intersection line between the fracture surface and the horizontal plane (0°-360°) and the dip angle is recorded as the angle between the fracture surface and the horizontal plane (0°-90°). For some scenarios, dip data (data perpendicular to the strike, indicating the dip direction of the fracture) can be added. Oversampling in densely populated fracture areas should be avoided to ensure that the sampling points evenly cover the area affected by the project, such as the area around the tunnel and the slope rock mass. Obvious abnormal data, such as dip angles exceeding the reasonable range and incorrect strike recordings, should be removed. A pole isodensity map transforms the orientation and dip of fractures in three-dimensional space into poles on a two-dimensional plane. Specifically, it uses stereographic projection: using an imaginary unit sphere as an intermediary, the intersection of the normal vector of each fracture (the vector perpendicular to the fracture surface) and the sphere is the pole. Through stereographic projection, starting from the north or south pole of the sphere, the pole is projected onto the equatorial plane, converting it into two-dimensional coordinates, typically polar or rectangular coordinates. The fracture orientation determines the azimuth of the pole on the plane; the azimuth may be the same as the orientation or differ by 180°, requiring consistency according to projection rules. The dip angle determines the distance from the pole to the center of the plane; the larger the dip angle, the closer the pole is to the edge of the plane; the smaller the dip angle, the closer the pole is to the center. Using professional geological software such as Dips, Surfer, ArcGIS, or programming tools like Python+Matplotlib, after inputting strike and dip data, the software directly generates two-dimensional coordinates of the poles. After projecting all fracture poles onto a two-dimensional plane, isodensity lines are drawn by statistically analyzing the spatial clustering of the poles, forming a pole isodensity map. The generation of isodensity lines involves classifying the statistically obtained density values ​​into categories such as low density, medium density, and high density. Typically, the number of poles per unit area or the percentage of relative density (e.g., the number of poles per unit area to the total number of poles) is used as an indicator. Points with the same density are connected to form closed isodensity lines. The higher the density of the area, the denser the isodensity lines, ultimately presenting a high-density concentration zone on the map. In this embodiment of the invention, two concentration zones are identified, corresponding to two sets of dominant structural fractures, clarifying the occurrence characteristics of the dominant fracture sets.

[0036] Stereographic projection transforms the orientation and dip of fractures in three-dimensional space into poles on a two-dimensional plane. The position of the pole on the plane directly corresponds to the dip angle of the fracture, which is the angle between the fracture surface and the horizontal plane, ranging from 0° to 90°. The larger the dip angle, the closer the pole is to the periphery of the plane; the larger the dip angle, the farther the pole is from the center of the plane and the closer it is to the periphery. The smaller the dip angle, the closer the pole is to the center. The radial distance of the pole on the stereographic map is positively correlated with the fracture dip angle. Through field measurements of 11,348 fractures, it was found that high-angle fractures with dip angles of 61°-90° accounted for 72.1%. The distribution of poles on the stereographic map is a visual mapping of the fracture dip angle distribution. The concentrated poles on the periphery represent the projected set of high-angle fractures; the sparse poles in the center represent the projected set of low-angle fractures. Therefore, the concentration of poles towards the periphery directly confirms that the rock mass in the study area is dominated by high-angle structural planes.

[0037] Step S2: Construction of the 3D Discrete Fragment Network (DFN) Model Based on the fracture statistical characteristic model established in step S1, this step aims to construct a three-dimensional DFN model that can realistically reflect the spatial distribution law of rock mass structural surfaces and meet the requirements of numerical calculation scale control. Taking a water conveyance tunnel project as an example, the model is first delineated... The computational domain is defined, and a 3D tunnel solid is established within it as a spatial constraint for fracture generation and trimming. A Monte Carlo stochastic simulation program is written in Python, with the maximum number of iterations for random sampling set to [value missing]. The program then uses the relative error between the simulated fracture volume density and the measured value within 5% as the convergence condition for terminating the loop. The program converts the field statistical parameters into a geometric model based on the volume density of each dominant fracture group and the background fracture. The number of fractures generated was determined using the Poisson distribution. The estimation relationship is as follows: in To generate the domain volume, Let be the average area of ​​the fractured disk. Fisher distribution is used to randomly sample the fracture strike and dip angle to reproduce the anisotropy of the attitude; its probability density function is expressed as: In the formula, It is the probability density function of the Fisher distribution, used to describe the random distribution of fracture attitude, including strike and dip angle. In step S2, the construction of the three-dimensional DFN model is used. Based on the statistical data of fracture attitude measured in the field, such as the strike and dip angle of the dominant fracture group, fracture orientation that conforms to the actual distribution is generated by random sampling. Combined with the Fisher distribution, the strike and dip angle of the fracture are randomly selected to avoid homogenization of fracture attitude and restore the distribution characteristics of the real rock mass, where the dominant fracture group dominates and a small number of scattered fractures supplement. The sampling results directly determine the spatial orientation of each fracture in the DFN model, providing a realistic geometric basis for subsequent mechanical response simulation. Fisher's constant, The angle between the crack normal vector and the average vector is given. Based on the log-normal distribution, the equivalent radius is extracted, and the crack is simplified to a disk element. The centroid of the disk is generated within a space extending 20% ​​beyond the boundary of the computational domain. The coordinates of the centroid of the simplified disk element are... Strictly adhering to a three-dimensional uniform distribution to eliminate boundary void effects, the disks are geometrically constructed as zero-thickness patches, but are given a statistically consistent hydraulic gap width attribute. When handling fracture intersection logic, the program automatically identifies intersection lines between disks, retaining valid fractures forming a connected network and marking isolated fractures. Finally, invalid fractures completely outside the computational domain or with no interaction with engineering entities are removed, and boundary truncation effect correction is implemented to complete the construction of the DFN model.

[0038] Model validation results are as follows Figure 3 and Figure 4 As shown: The DFN model intuitively reproduces the steeply inclined ( , percentage The simulated fracture spatial morphology, dominated by [missing information], clearly shows its geometric distribution and spatial relationship with the water conveyance tunnel. Quantitative statistics further indicate that the dip angle (mean) of the simulated fractures... ), Predominant orientation and size distribution (mean) The results are in high agreement with the measured data. This indicates that the constructed DFN model accurately represents the actual fracture system in a statistical sense and can provide high-precision geometric input for subsequent numerical simulations such as FLAC3D.

[0039] In this invention, establishing a three-dimensional tunnel entity as a spatial constraint for crack generation and trimming clarifies the scope of the tunnel and the surrounding rock mass. Subsequent cracks can only appear within this rock mass. Furthermore, if a crack crosses the tunnel, it must be trimmed to match the tunnel's shape; for example, if the tunnel is circular, the crack is cut off at the circular boundary to avoid discrepancies between the generated cracks and the actual engineering scope. A random generator is written in Python. This generator operates on a structured randomness; since most cracks in reality have steep angles, the program will not randomly generate flat cracks. The program repeatedly attempts to generate cracks; in this embodiment, it is set to attempt a maximum of 100,000 times. After each batch of cracks is generated, the crack density is calculated (e.g., how many cracks are in 1 cubic meter of virtual rock mass). The simulated density is compared with the actual density measured on-site. If the error is ≤5% (e.g., the actual density is 10 cracks / cubic meter, and the simulated result is 9.5-10.5 cracks / cubic meter), the attempt is stopped. The process of converting on-site statistical parameters into a geometric model involves first going to the construction site and actually measuring many cracks. For example, in this embodiment of the invention, a total of 11,348 cracks were measured. Key information recorded includes: ① main crack groups (e.g., two groups of cracks are the main ones); ② crack volume density (e.g., the number of cracks per unit volume). The program uses the measured data as rules, combined with the computational domain and the scale control requirements of the numerical simulation, to reasonably optimize the number of cracks generated through Poisson distribution. This avoids reduced computational efficiency due to excessive data volume. Poisson distribution is a classic probability distribution in statistics used to describe the number of random events occurring within a unit of time or space. Using Poisson distribution to calculate the random number of cracks generated in a certain volume of rock mass in stochastic modeling in fields such as geotechnical engineering and hydraulic engineering is an existing technology and will not be elaborated here. Taking the computational domain of 100m×100m×120m (volume 1.2 million cubic meters) in this embodiment as an example, although the initial estimate of the total theoretical number of cracks based on the measured density is relatively large, considering the balance between computational efficiency and accuracy in numerical simulation, effective cracks are selected, and isolated cracks, cracks outside the computational domain, and invalid cracks that do not interact with the engineering entity are eliminated. Finally, 800 effective cracks that conform to statistical laws are generated. According to the statistical distribution of measured data, 70% belong to the first group of cracks, 30% belong to the second group and scattered cracks. This invention draws out the scope of the tunnel and the surrounding rock mass, and defines the generation rules. The crack data measured on site, such as density and main direction, are input into the program. The program generates cracks one by one until the number of cracks and the error between the actual number and reality are ≤5%. Finally, a 3D model that is almost consistent with the distribution of cracks in the real rock mass is obtained, so that the virtual cracks can reproduce the cracks in reality as much as possible. In subsequent simulated construction, it can be determined where collapse is likely and how many anchor bolts are needed for reinforcement.

[0040] Where is Fisher's constant, the The values ​​are explained as follows: Fisher's constant in this invention The value is an index describing the degree of concentration of fracture orientation, such as strike and dip angle. The larger the value, the more concentrated the fracture attitude, such as dip angle and strike, that is, the more obvious the dominant fracture group; The smaller the value, the more dispersed the fracture orientation is, with no obvious dominant direction. Its specific value is not fixed but is determined entirely based on fracture data measured in the study area. The fracture distribution varies from site to site. The values ​​will inevitably be different. (Determined) First, sufficient field fracture data must be collected. For example, technicians visit sites such as tunnels and slopes to systematically measure and record the key attitudes of the fractures, primarily their strike (the direction of the intersection of the fracture surface and the horizontal plane) and dip angle (the angle between the fracture surface and the horizontal plane). In this embodiment, a total of 11,348 samples were collected. After obtaining the field data, professional software such as geological analysis software Dips or Surfer, or programming tools like Python + SciPy, were used to transform the attitudes of all fractures into direction vectors in three-dimensional space, similar to turning the dip direction of each fracture into an arrow. A Fisher distribution model was then used to fit the distribution of these direction vectors, and adjustments were made. The value is adjusted until the directional distribution calculated by the model most closely matches the field-measured fracture orientation distribution. For example, if the model predicts that 70% of the fractures have a dip angle of 60-90°, which is almost identical to the field-measured 72.1%, then the model is considered to be the closest match to the measured data. The value is the Fisher constant at the project site.

[0041] In this invention, extracting the equivalent radius refers to randomly selecting a value as the equivalent radius of the fracture using a log-normal distribution probability statistical method during the construction of the 3D DFN model. Each fracture is then simplified into a disk element, thereby quantifying the size of the fracture and standardizing its geometric shape. This approach balances the randomness of fracture size with the operability of modeling. Fractures in nature are irregular in shape, with no fixed pattern in length, width, or extension, making them unsuitable for direct geometric modeling. This invention uses the equivalent radius to represent the radius of a circular disk, allowing the calculation of the disk area (equivalent to the fracture's equivalent area). This standardizes the spatial size of the fracture, enabling computer models to recognize and calculate irregular fractures. Extensive engineering measurement data shows that the distribution of fracture size, length, and radius in rock masses conforms to a log-normal distribution, meaning that most fracture sizes are concentrated near the mean, while a few fracture sizes are larger or smaller. The radius is extracted using a log-normal distribution, not by arbitrary assignment. Instead, values ​​are randomly selected based on the actual statistical laws governing crack sizes in nature. This ensures that the simulated crack size distribution closely matches the actual engineering situation, avoiding model distortion due to deviations in size settings from reality. After extracting the equivalent radius, each crack is abstracted as a disk with zero thickness and the equivalent radius as its radius. The planar shape of the disk corresponds to the actual extension range of the crack, the equivalent radius determines the size of the crack, and the spatial orientation and inclination angle of the disk are determined by Fisher distribution sampling.

[0042] Step S3: Numerical simulation of rock mass mechanical response based on DFN model This step uses the high-fidelity DFN model generated in step S2 as the geometric core, embedding it into a numerical computing platform such as FLAC3D, aiming to analyze the rock mass mechanical behavior controlled by real fractures. First, the DFN model is mapped to joint or structural surface elements to construct a discontinuous "fracture-rock block" medium system, and Mohr-Coulomb constitutive parameters of the surrounding rock are assigned based on engineering geological surveys. Taking a tailrace tunnel of a hydropower station crossing a highway section as an example, a... Computational model, As the three-dimensional dimension of the computational domain during numerical simulation, the simulation process encompasses initial in-situ stress inversion, tunnel excavation and support, and operational load loading. For the critical initial in-situ stress inversion stage, a multiple regression analysis method based on field measured data is employed. Hydraulic fracturing method measurement data is selected as the fitting target, and the normal displacement of the model boundary is set as the constraint condition to determine the initial in-situ stress field within the computational domain. This can be expressed as a linear superposition of the self-weight stress field and the tectonic stress field. In the formula The regression coefficients are obtained by solving the least squares method, thereby ensuring that the residuals between the inverted stress field and the measured values ​​meet the engineering accuracy requirements. This represents the self-weight stress, which is the pressure exerted by the mountain above on that point. Self-weight stress = weight of the rock mass above ÷ area under stress. It represents tectonic stress in the X direction, such as the compressive or tensile force generated by crustal movement in the horizontal X direction; This represents tectonic stress in the Y direction, such as the compressive or tensile force generated by crustal movement in the horizontal Y direction. Tectonic stress is measured through field tests, often using the hydraulic fracturing method. A small hole is drilled in the rock mass, and pressure is applied into the hole until the rock cracks. By measuring the pressure and direction of the cracks, the tectonic stress in the X and Y directions of the crust can be inferred. It is a contribution coefficient that adjusts the influence weights of self-weight and X / Y tectonic stress, and determines the contribution ratio of the three types of stress to the total initial geostress. This is a correction constant that corrects for minor deviations between the model assumptions and actual field conditions. Multiple measuring points, such as 10, are selected at the engineering site. The actual initial geostress at each measuring point is directly measured using the hydraulic fracturing method, and then substituted into the aforementioned initial geostress field. The formula establishes a system of equations: For each measurement point, an equation can be listed; 10 measurement points correspond to 10 equations. The least squares method is used for fitting, i.e., solving this system of equations using professional software such as Excel, Matlab, or FLAC3D, to find a suitable system. , To minimize the error between the calculated values ​​and the actual measured values ​​at all measuring points, i.e., to ensure the residuals meet engineering accuracy, the following formula is determined: , The value of .

[0043] The tunnel excavation and support were then simulated sequentially, followed by the application of equivalent concentrated forces during the operation period according to a two-lane standard. and uniformly distributed load The highway load effect.

[0044] Numerical simulation results are as follows Figure 5 As shown: Under the condition of excavation only, the surrounding rock exhibits uniform convergence, with a total displacement of approximately After the superimposed highway load, the additional deformation of the tunnel roof and the densely cracked area increased significantly, with the maximum incremental displacement exceeding [missing information]. Furthermore, the deformation is primarily transmitted along the dominant fracture surface. Through comparative analysis, this method not only achieves continuous simulation of the entire process from initial equilibrium to operational load, but also uses the "incremental displacement" index to quantitatively separate the influence of specific loads on rock mass stability, thus providing a reliable basis for support optimization and safety evaluation in complex fractured rock mass engineering.

[0045] In this invention, the DFN model is mapped to joint or structural surface elements. The location, shape, and direction of cracks in the DFN model data are consistent with the real rock mass. The joint or structural surface elements are virtual crack elements that the computer program can recognize; that is, each real crack corresponds to a virtual crack in the computer program. The program stores the location, size, and shape of each crack. In a real rock mass, a crack divides a large rock mass into two rock blocks, A and B. The computer model clearly defines that blocks A and B are separate, and the crack in the middle is their boundary. Subsequent calculations consider whether block A will slide along the crack to block B. Based on engineering geological surveys, the Mohr-Coulomb constitutive parameters of the surrounding rock are assigned to the rock blocks and cracks in the model, and rules for stability are set respectively. For example, a maximum pressure value that the rock block can withstand is set, exceeding which will cause it to break; a maximum tensile force value that the rock block can withstand is set, exceeding which will cause it to fracture; a friction coefficient of 0.6 is set for the crack surface, and a maximum shear force value that the crack can withstand is set, exceeding which the two rocks will slide along the crack. Traditional models treat rock masses as solid blocks, which can lead to errors in calculating localized collapses caused by cracks. For example, in reality, a single rock might slide down a crack, but traditional models would calculate it as a uniform deformation of the entire rock mass. In contrast, the crack-rock block model of this invention restores the structure of the real rock mass. Combined with the stability rules determined by the Mohr-Coulomb constitutive parameters of the surrounding rock, this invention can calculate results that are closer to reality. This allows for the determination of reinforcement methods in actual engineering projects, such as adding anchor bolts in areas with dense cracks to prevent rocks from sliding.

[0046] The initial geostress field within the computational domain The initial geostress field is expressed as a linear superposition of the self-weight stress field and the tectonic stress field, aiming to recreate the actual initial stress environment at the engineering site. This ensures that the initial conditions in the computer program simulation are completely consistent with the stress state of the actual rock mass, guaranteeing the accuracy of simulation results for subsequent tunnel excavation and operational loads. The initial geostress field represents the stress borne by the rock mass before excavation and construction. For example, the weight of the mountain above pressing on the rock mass below is the self-weight stress; the compressive and tensile forces generated by crustal movement are tectonic stresses, such as the compression of the rock mass by plate tectonics. The superposition of these two stresses constitutes the initial geostress field of the rock mass. The initial geostress field is the fundamental stress state affected by all subsequent engineering projects. For example, tunnel excavation will disrupt this balance, leading to rock mass deformation and cracking. The calculation sets the starting point for subsequent simulations. The program simulates the rock mass response, such as whether it will collapse after excavation. It is necessary to obtain the rock mass stress value under the initial state first.

[0047] Step S4: Verification of simulation results and demonstration of technological advantages Building upon the aforementioned modeling and numerical analysis, this step further clarifies the rationality and engineering value of the DFN-based fracture modeling method through results verification. Figure 6Taking the tailrace tunnel axial profile as an example, the simulation results reveal the non-uniform distribution characteristics of the plastic zone in the surrounding rock: near the free face of the tunnel, tensile failure is the main manifestation of unloading during excavation, while the deep surrounding rock mainly forms shear failure along the dominant fracture zone. Unlike the symmetrical smooth "butterfly-shaped" plastic zone commonly found in traditional equivalent continuous medium models, the failure morphology predicted by this model is irregular blocky or toothed, and is highly coupled with the spatial geometry of the three-dimensional DFN fracture network. This indicates that fractures substantially dominate the stress redistribution and failure process in the model, and can realistically reproduce discontinuous deformation phenomena such as block slippage and wedge instability.

[0048] Based on the above analysis of the failure mechanisms, the integrated "fracture numerical analysis" method proposed in this invention can provide precise basis for support optimization: for shallow areas where tensile failure is the main cause, it is recommended to use systematic anchor bolts and shotcrete to strengthen the integrity of the rock mass; for shear failure areas where shear failure is concentrated along the fractures, long anchor cables or shear piles should be arranged to effectively "stitch" the potential slip surface. In summary, this invention effectively overcomes the shortcomings of traditional continuous medium methods in reflecting the influence of real joint networks by organically linking on-site fracture statistics, three-dimensional DFN modeling, and discontinuous numerical simulation, providing a practical and highly reliable technical solution for the stability evaluation and support design of complex fractured rock mass engineering.

[0049] The above embodiments are merely illustrative examples for clear explanation and are not intended to limit the implementation. Those skilled in the art will recognize that other variations or modifications can be made based on the above description. It is neither necessary nor possible to exhaustively list all possible implementations. However, obvious variations or modifications derived therefrom are still within the scope of protection of this invention.

Claims

1. A method for modeling and simulating three-dimensional discrete fracture networks based on the statistical properties of fractures, characterized in that, The method includes: S1: Statistical characterization of the geometric parameters of the rock mass fracture system, including collecting the geometric parameters and distribution characteristics of rock mass fractures in the study area, constructing a multi-dimensional parameter system, and establishing a stochastic statistical model of rock mass fractures; the multi-dimensional parameter system includes statistical analysis of fracture geometric parameters, statistical analysis of fracture spatial distribution parameters, and parameter correlation analysis; the stochastic statistical model of rock mass fractures uses field measured data as input to quantify the stochastic characteristics and intrinsic correlation of key fracture parameters; S2: Construction of the three-dimensional discrete fracture network model. Based on the rock mass fracture stochastic statistical model established in step S1, the computational domain is delineated and a three-dimensional model of the engineering entity is established as a spatial constraint for fracture generation and trimming. The on-site statistical parameters are transformed into a geometric model through a Monte Carlo stochastic simulation program. After fracture intersection logic processing, invalid fracture removal and boundary truncation effect correction, the three-dimensional DFN model is completed. S3: Numerical simulation of rock mechanics response based on DFN model. The three-dimensional DFN model generated in step S2 is embedded into the numerical calculation platform, mapped to joint or structural surface units, and a discontinuous medium system of fractures and rock blocks is constructed. The surrounding rock is given Mohr-Coulomb constitutive parameters, and the initial in-situ stress inversion, engineering excavation and support simulation, and operation period load loading simulation are performed in sequence to obtain rock mechanics response data. S4: Verification and engineering application of simulation results. Based on the mechanical response data in step S3, analyze the failure morphology and distribution characteristics of the surrounding rock, optimize the support scheme in combination with actual engineering needs, and provide a basis for stability evaluation of complex fractured rock mass engineering.

2. The method for modeling and simulating three-dimensional discrete fracture networks based on fracture statistical properties according to claim 1, characterized in that: The method for constructing the three-dimensional discrete fracture network (DFN) model in step S2 includes: Determine the computational domain and spatial constraints: Based on the actual needs of the project, define the scope of the computational domain, establish a three-dimensional model of the engineering entity, clarify the spatial boundary of crack generation, and ensure that subsequent cracks only form within the computational domain of the rock mass. Write a Monte Carlo random simulation program: set the maximum number of iterations for random sampling and the convergence criteria; Parameter transformation and fracture generation: Based on the volume density parameter in the random statistical model of rock mass fractures, the total number of fractures generated is determined by Poisson distribution. Fisher distribution is used to randomly sample the fracture orientation and dip angle. The equivalent radius of the fracture is extracted based on the log-normal distribution, and the fracture is simplified into a disk element. Crack network optimization: Automatically identify the intersection lines between disk elements, retain valid cracks that form a connected network and mark isolated cracks, and remove invalid cracks that are completely outside the computational domain or have no interaction with engineering entities; Model calibration: Boundary truncation effect calibration is applied to the screened fractures to complete the construction of the three-dimensional discrete fracture network (DFN) model.

3. The method for modeling and simulating three-dimensional discrete fracture networks based on fracture statistical properties according to claim 1, characterized in that: The numerical simulation method for rock mass mechanical response based on the DFN model in step S3 includes: Model embedding and system construction: The three-dimensional DFN model generated in step S2 is embedded into numerical computing platforms such as FLAC3D, and the DFN model is mapped to joint or structural surface units to construct a discontinuous medium system of fractures and rock blocks; Assign Mohr-Coulomb constitutive parameters: Assign corresponding Mohr-Coulomb constitutive parameters to rock blocks and joints / structural planes in the system, respectively; Initial geostress inversion: The multivariate regression analysis method was adopted, and the hydraulic fracturing measurement data was selected as the fitting target. The normal displacement of the model boundary was set as the constraint condition, and the initial geostress field in the calculation domain was expressed as a linear superposition of the self-weight stress field and the tectonic stress field. Engineering excavation and support simulation: Based on the actual construction process of the project, the excavation process of tunnels and other structures is simulated, and the application of support measures is simulated simultaneously to restore the impact of excavation unloading and support reinforcement on the mechanical state of the rock mass. Operational load simulation: Apply corresponding loads according to the load requirements of the project operation phase, including equivalent concentrated force and uniformly distributed load applied according to the two-lane standard; Mechanical response data acquisition: The mechanical response data of the rock mass is output through the numerical calculation platform to complete the numerical simulation of the rock mass mechanical response. The mechanical response data includes, but is not limited to, the displacement field, stress field, and plastic zone distribution of the rock mass.

4. The method for modeling and simulating three-dimensional discrete fracture networks based on fracture statistical properties according to claim 1, characterized in that: In step S2, the number of fractures N generated is determined using a Poisson distribution, and the estimation relationship is N≈(V·P) 32 ) / E[A], where V is the volume of the generating domain, P 32 E[A] represents the volume density of each dominant fracture group and the background fracture, and E[A] represents the average area of ​​the fracture disk.

5. The method for modeling and simulating three-dimensional discrete fracture networks based on fracture statistical properties according to claim 1, characterized in that: In step S2, Fisher distribution is used to randomly sample the fracture orientation and dip angle. The probability density function of Fisher distribution... for: in Fisher's constant, The angle between the crack normal vector and the average vector is given.

6. The method for modeling and simulating three-dimensional discrete fracture networks based on fracture statistical properties according to claim 1, characterized in that: In step S2, the equivalent radius of the crack is extracted based on a log-normal distribution, and the crack is simplified into a disk element. The generation range of the disk centroid is set within a space 20% outside the boundary of the computational domain, and the centroid coordinates of the simplified disk element are... It follows a three-dimensional uniform distribution.

7. The method for modeling and simulating three-dimensional discrete fracture networks based on fracture statistical properties according to claim 1, characterized in that: In step S3, the initial geostress inversion will retrieve the initial geostress field within the computational domain. This can be expressed as a linear superposition of the self-weight stress field and the tectonic stress field. In the formula Indicates self-weight stress. Indicates structural stress in the X direction. Indicates structural stress in the Y direction. It is a contribution coefficient that adjusts the influence weight of self-weight and X or Y-direction tectonic stress. It is a correction constant.

8. A device, characterized in that, include: At least one processor; And a memory communicatively connected to the at least one processor; wherein the memory stores instructions executable by the at least one processor to enable the at least one processor to perform the method as described in any one of claims 1 to 7.

9. A computer-readable storage medium storing a computer program, characterized in that: When the computer program is executed by a processor, it implements the method of any one of claims 1 to 7.

Images