Fractured rock mass random field implementation method considering spatial variability under three-dimensional geological structure

Abstract

The invention discloses a fractured rock mass random field implementation method considering spatial variability under a three-dimensional geological structure. The method comprises the following steps: determining the size of a three-dimensional background grid region and the number of background nodes in each direction according to a rock mass region and the size of a rock mass unit; setting mathematical statistical characteristic variable values in the three-dimensional random field; a covariance matrix is constructed, a Karhunen-Loeve series expansion method is adopted for decomposition, and random field approximate values are generated at background grid points; based on the geological structure type of the rock mass area, performing spatial transformation operation and movement on the background grid; background grid points are quickly selected and judged, and random field values of rock mass units are obtained through shape function interpolation; the random field value is adjusted through mean value regression and standard deviation regression, and rock mass mechanical parameters are given; and when the random field ratio of the rock mass units on the two sides of the joint unit is judged to be larger than a set threshold value, weakening the mechanical property of the joint as a fracture. According to the invention, rock mass spatial variability under a three-dimensional geological structure can be simulated, and spatial distribution of fractures can be realized.

Description

Technical Field

[0001] The invention belongs to the field of rock mass geological structure simulation research, and specifically relates to a method for realizing a random field of a fractured rock mass taking into account the spatial variability under a three-dimensional geological structure. Background Art

[0002] Natural rock mass is affected by geological movement and stress history during its formation, and exhibits more complex spatial variability than soil mass. This variability is not only reflected in the mechanical behavior of the rock mass, but is also closely related to the various fracture surfaces formed inside it. The spatial variability and fracture surfaces of the rock mass often lead to higher instability and uncertainty in engineering construction, thus bringing many challenges to underground engineering and tunnel construction. As pointed out in the paper "He Manchao, Research Progress and Challenges of Deep Soft Rock Engineering, Journal of China Coal Society 39(8)(2014)1409-1417", underground excavation projects often encounter rock mass with naturally formed fractures, and many engineering examples show that rock mass failure is more dependent on fractures than intact rocks.

[0003] The main geological structure types of rock mass include parallel dip strata, synclines and anticlines, folds, faults, etc. These structural types show a high correlation with the spatial variability of rock mass. Random fields, as the main means to achieve spatial variability of rock and soil mass, have attracted the attention and application of a large number of scientific researchers. However, as a three-dimensional spatial body, most studies of rock mass use simplified two-dimensional solid models. For example, the Chinese invention patent, title: A rotational anisotropic non-stationary random field implementation method for rock and soil parameters, patent application number: 202110362390.7, reflects this point. As for the three-dimensional random field implementation method, many researchers have ignored the influence of geological structure, such as the paper "SUN Z, HUANG G, HU Y, et al. Reliability analysis of pile-reinforced slopes in width-limited failure mode considering three-dimensional spatial variation of soil strength. Computers and Geotechnics, 2023, 161: 105528.".

[0004] In the analysis of rock mass stability, cracks are also a factor that cannot be ignored. At present, the mainstream research method usually introduces cracks in the modeling stage. For example, the Chinese patent application, title: A discrete crack network generation method based on iterative function system, patent application number: CN202211146540.1, provides an effective method for crack modeling. In the study of the stability of anti-dip rock slopes, Yun Zheng et al. adopted a two-dimensional crack inclination model in the paper "Yun Zheng, Runfu Wu, Chengzeng Yan, Runqing Wang & Bin Ma. Numerical study on flexural toppling failure of rock slopes using the finite discrete element method. 2024, 83, 111.", and modeled different crack inclinations to analyze the influence of cracks on slope stability. However, this modeling method is not only time-consuming, but also simplifies the spatial distribution of cracks, and there is still a gap with the actual three-dimensional crack distribution. The crack modeling method not only consumes time, but also simplifies the cracks to a certain extent during modeling, which is different from the three-dimensional crack spatial distribution.

[0005] Therefore, further research and improvement are urgently needed to address the limitations of existing random field implementation methods. In particular, more accurate and applicable three-dimensional random field models need to be developed to effectively simulate the spatial variability and fracture distribution under complex geological structures, so as to better describe and predict the mechanical behavior of rock masses. Summary of the invention

[0006] The purpose of the present invention is to provide a method for realizing a random field of a fractured rock mass taking into account the spatial variability of a three-dimensional geological structure, so as to solve the technical difficulties such as the inability to consider the influence of the geological structure on the spatial variability in three-dimensional space and the distribution of fractures in three-dimensional space in the existing random field technology.

[0007] In order to solve the above technical problems, the present invention adopts the following technical solutions:

[0008] A method for realizing a random field of a fractured rock mass considering spatial variability under a three-dimensional geological structure comprises the following steps:

[0009] S1: Establish a three-dimensional rectangular coordinate system, set the length, width and height of the background area according to the size of the rock mass area, ensure that the background area can completely cover the rock mass area, and determine the number of grid points of the background area in three dimensions to generate a three-dimensional background grid volume;

[0010] S2: Set the mathematical statistical characteristic variables in the random field, including the mean, standard deviation, autocorrelation function and correlation length of the rock mass mechanical parameters;

[0011] S3: Construct a covariance matrix, decompose it using the Karhunen-Loeve series expansion method, approximate the eigenvalues ​​and eigenvectors, and then obtain the random field values ​​at each node in the 3D background grid;

[0012] S4: Based on the geological structure type of the rock mass area, a spatial transformation operation is performed on the background grid, and the transformed grid is moved to the center of the rock mass area to ensure that the rock mass area is located inside the transformed background grid;

[0013] S5: For any rock mass unit point, find the node of the background grid body corresponding to it, and calculate the random field value of the rock mass unit point by using the three-dimensional eight-node shape function interpolation method according to the random field values ​​of the eight background grid body nodes around the background grid body node;

[0014] S6: Perform mean regression and variance regression on the random field of rock mass unit points, eliminate the errors of the Karhunen-Loeve series expansion method and the shape function interpolation method on the mean and variance, and assign them to the mechanical parameters of the rock mass unit;

[0015] S7: inserting joint units between rock mass units and determining relevant mechanical parameters of the joint units;

[0016] S8: Set the threshold and joint weakening coefficient, adjust the mechanical parameters of the joint according to the ratio of the random field values ​​of the rock units on both sides of the joint unit, and effectively distinguish between joints and cracks.

[0017] Further optimization, in step S1, the width D of the constructed three-dimensional background area D x Length D y and height D z The sizes of the generated background mesh body are both larger than twice the corresponding sizes of the rock mass area, ensuring that after the background mesh body is deformed and translated in step S4, the background mesh body can still completely cover the entire rock mass area.

[0018] According to the number of rock mass units in the three-dimensional direction of the rock mass area, the width D of the three-dimensional background area D is set. x The number of nodes N x Length D y The number of nodes N y , height D z The number of nodes N z , ensuring that the spacing between nodes in the background grid is smaller than the rock mass unit size, the total number of grid nodes N in the background area D is expressed as N = N x *N y *N z .

[0019] Further optimization, in step S2, set any two points A(x i ,y j ,z k )、B(x l ,y m ,z n ) are X i,j,k and X l,m,n , 1≤i≤N x ,1≤j≤N y ,1≤k≤N z , 1≤l≤N x ,1≤m≤N y ,1≤n≤N z ; The rock mass mechanical parameters include the elastic modulus, tensile strength, etc. of the rock.

[0020] Mean μ X is the average value of the mechanical parameters of all rock mass units in the background area D, and the standard deviation σ X represents the mechanical parameters X at each grid node in the background area D i,j,k Distance from mean μ X degree of deviation.

[0021] The autocorrelation function is a function of the relative distance (Δx, Δy, Δz) between any two points A and B in the background area D. The rock mass mechanical parameters X i,j,k and X l,m,n The correlation will only decrease with the increase of relative distance and will approach a non-negative constant c, 0≤c<1.

[0022] The three-dimensional autocorrelation function ρ(X i,j,k ,X l,m,n ) is exponential and square exponential, and their expressions are shown in formulas (2.1) and (2.2) respectively:

[0023]

[0024]

[0025] Where Δx i,l =|x i -x l |,Δy j,m =|y j -y m |,Δz k,n =|z k -z n |, L x , L y and L z They represent the relevant lengths of rock mechanical parameters in the three dimensions respectively.

[0026] The correlation length is also called the fluctuation range, which means that when the distance between two points in each direction is less than the correlation length in the corresponding direction, the mechanical parameters will show a stronger correlation in that direction. x , L y and L z It is obtained by detecting the correlation in all directions of the rock mass area, which is an existing mature technology. Refer to the Chinese patent, title: A non-stationary random field modeling method for the mutual correlation of rock and soil parameters, application number: CN202211404802.X, step 2 of the patent is described in detail, which will not be repeated here.

[0027] Further optimization, in step S3, the covariance matrix is ​​constructed, and the random field value is calculated using the Karhunen-Loeve series expansion method. The main reason for choosing this method is that the number of three-dimensional background grid points is huge, making the traditional Cholesky decomposition method unsuitable for large matrix decomposition. Compared with the traditional Cholesky decomposition, the Karhunen-Loeve series expansion method has higher computational efficiency, and its calculation formula is shown in (3.1):

[0028]

[0029] Where X is the mechanical parameter X of the grid point in the background area D i,j,k The collection of X = {X i,j,k |1≤i≤N x ,1≤j≤N y ,1≤k≤N z}; RF(X) is the set of random field values ​​generated by the set X; μ X is the mean value, σ X is the standard deviation, p is the covariance matrix ρ NN The pth expansion term of the decomposition. The number of expansion terms determines the accuracy of the random field approximation. Generally speaking, a value of 20 can meet the requirements. Too large a value will reduce the efficiency of the series expansion method. Refer to Chapter 3.6.3 of the doctoral thesis "Risk Analysis of Instability of Layered Rock Caverns Based on Random Finite Difference". χ is the length of the total number of background grid nodes N = N x *N y *N z , and the random numbers satisfying the standard normal distribution; p is NN The eigenvalue set λ ijk After sorting in descending order, the eigenvalue of p is numbered, and φ p (X) corresponds to λ p The eigenvector of ijk and φ(X i,j,k ) satisfies formula (3.2):

[0030]

[0031] This formula uses the covariance matrix ρ NN As a kernel function, it is defined as ρ(X i,j,k ,X l,m,n ),N=N x *N y *N z ρ NN is the autocorrelation function value between any two nodes of the background grid, and its specific form is shown in formula (3.3):

[0032]

[0033] Among them, ρ(X i,j,k ,X l,m,n The specific expression of ) depends on the type of autocorrelation function, refer to formula (2.1) or (2.2); In addition, φ(X i,j,k ) are mutually orthogonal, and their integral with themselves is 1, as shown in formulas (3.4) and (3.5):

[0034] ∫ D φ(X i，j，k )φ(X l，m，n )dV=δ(X i，j，k , X l，m，n ) (3.4)

[0035]

[0036] Further optimizing, in step S4, the spatial transformation of the background mesh body includes spatial rotation, spatial stretching or spatial distortion;

[0037] The geological structure morphology is set by simulating the deformation of the background grid. The specific operations are as follows:

[0038] (a) The background grid is rotated in any direction to simulate the geological structure of the parallel rock layer. The rotation angle θ can be freely selected in the range of 0°-360°.

[0039] (b) Select any two nodes on the boundary surface of the background grid body to determine a straight line, and stretch the background grid body along the normal direction outside the boundary of the straight line to simulate the geological structure of the anticline or syncline.

[0040] (c) On any boundary surface of the background grid, distortion is performed along the normal direction of the boundary to simulate the folded geological structure.

[0041] Further optimization, the step S5 specifically includes the following steps:

[0042] S5.1: For any rock mass unit point Y t (x' t ,y' t ,z' t ), find the unique corresponding background mesh node (x i ,y j ,z k ), so that it satisfies: x i-1 <x' t <x i ,y j-1 <y' t <y j ,z k-1 <z' t <z k , t=1,2...s, s represents the number of rock mass unit points in the rock mass area.

[0043] S5.2: Replace {(x a ,y b ,z c )|i-1≤a≤i,j-1≤b≤j,k-1≤c≤k} are used as the parent grid unit points, and the random field values ​​of the corresponding rock mass unit points are calculated using the three-dimensional eight-node shape function interpolation method. The calculation formula is as follows (5.1):

[0044]

[0045] Among them, the shape function T abc (ε a ε,η b η,ζ c ζ)=(1 / 8)*(1+ε a ε)(1+η b η)(1+ζ c ζ), ε i-1 =η j-1 =ζ k-1 =-1,ε i =η j =ζ k =1; (ε,η,ζ) is about (x' t ,y' t ,z' t ), the relationship between the two is (ε,η,ζ)=[J] -1 (x' t ,y' t ,z' t ), J is the Jacobian determinant, which is calculated by formula (5.2):

[0046]

[0047] Further optimization, in step S6, the error caused by the Karhunen-Loeve series expansion method is eliminated by mean regression and standard deviation regression, and the specific processing steps are as follows:

[0048] S6.1: In mean regression, the actual mean μ of the rock mass regional random field RF(Y) is calculated according to formula (6.1): real ; Then, according to formula (6.2), each element of RF(Y) t ) multiplied by the set mean μ X With μ real The ratio of , we get the random field of the rock mass unit point after regression:

[0049]

[0050] Among them, RF(Y t ) new RF(Y k ) After mean regression processing, at the rock mass unit point Y t (x' t ,y' t ,z' t ) at the updated value; RF(Y) new = {RF(Y t ) new |t=1,2...s}, is RF(Y t ) new The set of is the rock mass regional random field after mean regression processing; μ X It is the mean value of rock mass parameters, which is obtained through actual detection and calculation of rock mass.

[0051] S6.2: In standard deviation regression, RF(Y) is calculated according to formula (6.3) new The actual standard deviation σ real ; Then, according to formula (6.4), calculate RF(Y t ) new With μ X The difference Dis_RF(Y t ); Finally, according to formula (6.5),

[0052] Calculate the random field value of each rock mass unit point after standard deviation regression processing;

[0053]

[0054] Dis_RF(Y t )=RF(Y t ) new -μ X (6.4)

[0055]

[0056] Among them, RF(Y t ) final RF(Y t ) new After standard deviation regression processing, at the rock mass unit point Y t (x' t ,y' t ,z't) update value; RF(Y) final = {RF(Y t ) final |t=1,2...s}, is RF(Y t ) final The set of is the final value of the random field in the rock mass area; σ X is the standard deviation of rock mass parameters, which is obtained through actual detection and calculation of the rock mass.

[0057] When performing mean regression and standard deviation regression, the order of operations must remain unchanged: mean regression should be performed before standard deviation regression; if standard deviation regression is performed first and then mean regression, the standard deviation will change, affecting the accuracy of the final result.

[0058] Further optimization, the step S7 specifically includes the following steps:

[0059] S7.1: For any two rock mass units with a common surface, generate joint units between them that are the same as the common surface, so that the rock mass units become discontinuous with each other.

[0060] S7.2: Set the relevant mechanical parameters Z of the joint element g Equal to the mean value μ of the corresponding rock mass mechanical parameters X , that is, Z={Z g =μ X |g=1,2,…G}, where Z g Represents the mechanical parameters of the gth joint unit, and G is the number of joint units.

[0061] By inserting joint units between rock mass units, the rock mass area is transformed into a discontinuous model, and the relevant mechanical parameters of the joint units are determined. This method can effectively simulate the mechanical behavior of rock mass in a fractured state and is suitable for finite element method FDEM, discrete element method DEM and discontinuous analysis method DDA.

[0062] Further optimization, the step S8 specifically includes the following steps:

[0063] S8.1: Based on the standard deviation σ of rock mass mechanical parameters X , set the threshold According to the mechanical parameters of the rock mass and the mechanical parameters of the cracks, the joint weakening coefficient W is determined;

[0064] S8.2: Calculate the ratio of the random field values ​​of the rock mass units on both sides of the joint unit and make a judgment:

[0065] If the ratio is less than the threshold R Y , then multiply the joint mechanical parameters by the average value of the random field values ​​of the rock mass units on both sides;

[0066] If it is greater than or equal to the threshold R Y , then multiply the joint mechanical parameters by the average value of the random field values ​​of the rock mass units on both sides, and then multiply them by the joint weakening coefficient to weaken the joint mechanical parameters to the crack level; W varies according to the density and development of the cracks in the rock mass, and is generally taken as 0.1%-10%.

[0067] Adjusting the mechanical parameters of the joint according to the ratio of the random field values ​​of the rock units on both sides of the joint unit is more in line with the actual situation. Due to the difference in stress history, the mechanical properties on both sides of the rock mass may be significantly different, and when this difference is too large, cracks are more likely to form.

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

[0069] 1. Based on the prior application (202411201166X, name: A random field modeling method for spatial variability of rock mass parameters under complex geological structures), the present invention continues to study by improving the traditional random field method, including operations such as transformation of the random field background grid and grid correction, on the basis of realizing the spatial variability under a variety of common three-dimensional geological structures, and by introducing joint units, setting thresholds and other operations, the spatial distribution of cracks is realized, the crack modeling operation is avoided, and the technical problems such as the inability to consider the influence of geological structure on spatial variability in three-dimensional space and the distribution of cracks in three-dimensional space in the existing random field technology are solved.

[0070] 2. The present invention can be effectively applied to numerical simulation methods such as finite discrete element method, discrete element method, and non-continuous analysis method, providing a reliable implementation technology for simulation analysis of related geological structures. BRIEF DESCRIPTION OF THE DRAWINGS

[0071] Figure 1 A flowchart of the method for realizing the random field of fractured rock mass with spatial variability under three-dimensional geological structure;

[0072] Figure 2 Schematic diagram of the initial background grid constructed in step S4 and three spatial transformation operations;

[0073] Figure 3 A schematic diagram of inserting a joint unit in step S7;

[0074] Figure 4 It is a case map of the rock mass area and a map of the rock mass unit of the case;

[0075] Figure 5 It is the realization effect and real picture of the spatial variability of fractured rock mass under three-dimensional parallel stratigraphic structure;

[0076] Figure 6 It is the realization effect and real picture of the spatial variability of fractured rock mass under three-dimensional anticline geological structure;

[0077] Figure 7 This is the realized effect and real picture of the spatial variability of fractured rock mass under folded geological structure. DETAILED DESCRIPTION

[0078] In order to make the purpose, technical solution and advantages of the embodiments of the present invention clearer, the technical solution of the present invention will be clearly and completely described below in conjunction with the accompanying drawings. Obviously, the described embodiments are part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0079] Embodiment 1:

[0080] In this embodiment, for a fractured rock mass model with spatial variability under a parallel rock layer structure with an inclination angle of 30°, a method for realizing a fractured rock mass random field with spatial variability under a three-dimensional geological structure is specifically implemented by the following steps: Figure 1 As shown:

[0081] S1: Establish a three-dimensional rectangular coordinate system, set the length, width and height of the background area according to the size of the rock mass area, ensure that the background area can completely cover the rock mass area, and determine the number of grid points of the background area in three dimensions to generate a three-dimensional background grid.

[0082] In this embodiment, the length, width and height of the rock mass area are set to be 45m, 45m and 30m respectively, and the shape of the rock mass unit is a tetrahedron. The engineering modeling personnel set its size between 1-2m and the number of units to 85339. Figure 4 As shown in the case of Figure 4 (a) is a schematic diagram of the rock mass area. Figure 4 (b) in the figure is the rock mass unit modeling diagram.

[0083] To ensure that the background area can contain the rock area no matter how it is transformed, its area size should be twice as large as the rock area. Let the width D of the three-dimensional background area D be x =90m, length D y =90m and height Dz = 90m, the background grid unit size is set between 1m, that is, x i -x i-1 =1m,y j -y j-1 =1m,z k -z k-1 =1m, the maximum number of nodes of the background grid is N=N x *N y *N z =91×91×91=753571.

[0084] S2: Mathematical statistical characteristic variables in the statistical random field, including the mean, standard deviation, autocorrelation function and correlation length of rock mass mechanical parameters:

[0085] In this embodiment, the mean tensile strength of the rock mass in the region is detected to be 100.0 MPa, and the standard deviation is about 0.1 times the mean, i.e., 10.0 MPa. The autocorrelation function uses an exponential type. The rock mass is divided into multiple layers vertically, so the correlation is small, and is set to L z = 2m, and has a high correlation in both the horizontal and vertical directions, set as the correlation length L in the length and width directions of the model x =L y = 45m, and considering that the rock mass beyond the correlation length still has a small correlation, the lower limit of the exponential autocorrelation function is set at 0.1. The rock layers are parallel and have an angle of 30° with the horizontal ground.

[0086] Set any two points A(x i ,y j ,z k )、B(x l ,y m ,z n ), combined with the above information, we can conclude:

[0087] μ X =100.0MPa,σ X =10.0Mpa

[0088]

[0089] S3: Construct the covariance matrix, decompose it using the Karhunen-Loeve series expansion method, approximate the eigenvalues ​​and eigenvectors, and then obtain the random field values ​​at each node in the three-dimensional background grid.

[0090] In this embodiment, the number of nodes in the three-dimensional background grid is 753571, and the traditional Cholesky decomposition method is not suitable for large matrix decomposition. Compared with the traditional Cholesky decomposition matrix, the Karhunen-Loeve series expansion method is significantly more efficient in dealing with large matrix decomposition, and its calculation formula is as follows:

[0091]

[0092] Where X is the tensile strength X of the grid points in the background area D i,j,k The collection of X = {X i,j,k =X i,j (x i ,y j ,z k )|1≤i≤91,1≤j≤91,1≤k≤91}; RF(X) is the set of random field approximations generated by the set X; 100.0 is the mean, 10.0 is the standard deviation, and p represents the covariance matrix ρ in formula (3.3) NN The pth expansion term of the decomposition. The number of expansion terms determines the accuracy of the random field approximation. Generally speaking, a value of 20 can meet the requirements. Too large a value may reduce the efficiency of the series expansion method. Refer to Chapter 3.6.3 of the doctoral thesis "Risk Analysis of Instability of Layered Rock Caverns Based on Random Finite Difference". χ is the length of the total number of background grid points N = N x *N y *N z =753571, and the random number that satisfies the standard normal distribution; p is NN The eigenvalue set λ ijk After sorting in descending order, the eigenvalue of p is numbered, and φ p (X) corresponds to λ p The eigenvector of ijk and φ(X i,j,k ) satisfies formula (3.2); this formula uses the covariance matrix ρ NN ρ(X i,j,k ,X l,m,n ) as kernel function, N = 753571; ρ NN is the autocorrelation function value between any two nodes of the background grid, as shown in formula (3.3).

[0093] S4: Based on the geological structure type of the rock mass area, perform a spatial transformation operation on the background mesh and move the transformed mesh to the center of the rock mass area to ensure that the rock mass area is located inside the transformed background mesh:

[0094] In this embodiment, the center of the background area is (45,45,30), and it is moved to the center of the rock area (22.5,22.5,15); then, the background grid is rotated in any direction to achieve the parallel rock layer geological structure: rotate at any inclination angle θ = 30°. Figure 2 In (a) and (b), (a) is the initial background mesh, and (b) is a schematic diagram of the initial background mesh after spatial rotation.

[0095] S5: For any rock mass unit point, find the node of the background grid body corresponding to it, and calculate the random field value of the rock mass unit point by using the three-dimensional eight-node shape function interpolation method according to the random field values ​​of the eight background grid body nodes around the background grid body node, which specifically includes the following steps:

[0096] S5.1: Mechanical parameters Y for rock mass element points t (x' t ,y' t ,z' t ), we need to find a solution that satisfies (x i-1 <x' t <x' i y j-1 <y' t <y j ,z k-1 <z' t <z k )'s background mesh node (x i ,y j ,z k ), and due to (x' t ,y' t ,z' t ) is inside the background grid, there is only one (x i ,y j ,z k ) meets the conditions, and in the judgment process, i-1 <x' t <x i ,y j-1 <y' t <y j ,z k-1 <z' t <z k By eliminating any background grid nodes that meet the requirements, you can quickly select (x i ,y j ,z k ); where (x' t ,y' t ,z' t ) represents the mechanical parameter Y of the rock mass unit pointt The position of , t=1,2...s, s=85339 is the number of rock mass units in the background grid area D.

[0097] S5.2: Replace {(x a ,y b ,z c )|i-1≤a≤i,j-1≤b≤j,k-1≤c≤k} corresponding to the eight nodes are used as the parent grid units, and the three-dimensional eight-node shape function interpolation method is used to calculate the random field values ​​of the corresponding rock unit points. The calculation formula is as shown in (5.1).

[0098] S6: Perform mean regression and variance regression on the random field of rock mass unit points to eliminate Karhunen-Loev e The series expansion method and shape function interpolation method are used to calculate the errors in the mean and variance and assign them to the mechanical parameters of the rock mass unit.

[0099] In this embodiment, the specific steps are as follows:

[0100] S6.1: In mean regression, the actual mean μ of the rock mass regional random field RF(Y) is calculated according to formula (6.1): real ; Then, according to formula (6.2), each element of RF(Y) t ) multiplied by the set mean μ X With μ real The ratio of , we get the random field of the rock mass unit point after regression:

[0101]

[0102] Among them, RF(Y t ) new After mean regression, at the rock mass unit point Y t (x' t ,y' t ,z' t ) at the updated value; RF(Y) new = {RF(Y t ) new |t=1,2...85339},is RF(Y t ) new The set of , that is, the rock mass regional random field after mean regression processing; μ X It is the mean value of rock mass parameters, which is obtained through actual detection and calculation of rock mass.

[0103] S6.2: In standard deviation regression, RF(Y) is calculated according to formula (6.3) new The actual standard deviation σ real ; Then, according to formula (6.4), calculate RF(Y t )new With μ X The difference Dis_RF(Y t ); Finally, according to formula (6.5),

[0104] Calculate the random field value of each rock mass unit point after standard deviation regression processing;

[0105]

[0106] Dis_RF(Y t )=RF(Y t ) new -μ X , t=1,2...85339 (6.4)

[0107]

[0108] Among them, RF(Y t ) final RF(Y t ) new Rock mass unit Y after standard deviation regression t (x' t ,y' t ,z' t ) Update value; RF(Y) final = {RF(Y t ) final |t=1,2...85339},is RF(Y t ) final The set of , satisfies the set mean and standard deviation. In addition, the order of mean regression and standard deviation regression cannot be changed. The random field value is assigned to the tensile strength of the rock mass unit. At this time, the tensile strength of the rock mass unit already has the spatial variability under the geological structure.

[0109] S7: Insert joint units between rock mass units and determine relevant mechanical parameters of the joint units.

[0110] In this embodiment, the following steps are specifically included:

[0111] S7.1: For any two rock mass units with a common surface, generate joint units between them that are the same as the common surface, so that the rock mass units become discontinuous with each other, such as Figure 3 As shown, Figure 3 (a) is a continuous model. Figure 3 (b) is a discontinuous model.

[0112] S7.2: The number of joint elements inserted in this model is 337514, and the relevant mechanical parameters Z of the joint elements are set. gEqual to the mean value μ of the corresponding rock mass mechanical parameters X , that is, Z = {Z g =100|g=1,2,…337514}, where Z g Represents the mechanical parameter value of the gth joint unit.

[0113] S8: Set the threshold and joint weakening coefficient, adjust the mechanical parameters of the joint according to the ratio of the random field values ​​of the rock units on both sides of the joint unit, and effectively distinguish between joints and cracks.

[0114] In this embodiment, the following steps are specifically included:

[0115] S8.1: Based on the standard deviation of rock mass mechanical parameters of 10.0, which is 0.1 times the mean value, and considering the density of rock mass cracks, the threshold R of the ratio of the random field values ​​of rock mass units is set. Y =1-0.1 2 =0.99; In addition, the strength at the crack is assumed to be 1% of the joint strength, and the mean tensile strength of the rock mass is μ X =100MPa, so the joint weakening coefficient W is set to 0.01;

[0116] S8.2: Calculate the ratio of the random field values ​​of the rock mass units on both sides of the joint unit and make a judgment:

[0117] If the ratio is less than 0.99, the joint mechanical parameters are multiplied by the average of the two random field values; if it is greater than or equal to 0.99, the joint tensile strength Z is multiplied by the average of the two random field values ​​and then multiplied by the joint weakening coefficient 0.01 to weaken it to the tensile strength level of the crack.

[0118] After the above processing, the spatial variability of fractured rock mass under three-dimensional parallel stratigraphic structure is realized, such as Figure 5 As shown in (a)-(c) in the figure. Figure 5 (a) shows the spatial variability of rock mass units parallel to the rock layer structure; Figure 5 (b) is the spatial variability of the joint units containing fractures in the geological structure. Figure 5 (c) in the figure is the cross-section of the geological structure; 5(d) in the figure is the real geological structure map. Figure 5 As can be seen from the comparison between (a)-(c) and (d), the present invention can generate a three-dimensional parallel stratum subsurface spatial variability fractured rock model that is close to reality.

[0119] Embodiment 2:

[0120] In this embodiment, for a fractured rock mass model with spatial variability under anticline structure, a fractured rock mass random field considering spatial variability under three-dimensional geological structure is calculated.

[0121] Using the case model of Example 1, Figure 4 As shown, the geological structure of the area is changed to anticline structure. The sizes of the rock mass area are still 45m, 45m and 30m respectively. The rock mass unit is tetrahedral in shape, with a size between 1-2m and a number of 85339, while the number of joint units is 337514. To ensure that the background area can contain the rock mass area no matter how it is transformed, its area size must be twice that of the rock mass. Let the width D of the three-dimensional background area D be x =90m, length D y =90m and height D z = 90m, the background grid unit size is set between 1m, that is, x i -x i-1 =1m,y j -y j-1 =1m,z k -z k-1 =1m, the maximum number is N=N x *N y *N z =91×91×91=753571. The mean value of rock mass tensile strength is also expressed in μ X =100.0MPa, standard deviationσ X =10.0MPa, related length L x =L y =45m,L z =2m; the autocorrelation function is also exponential, with a lower limit of c=0.1.

[0122] Any two nodes on the boundary surface of the background grid are selected to determine a straight line, and the background grid is stretched along the normal direction outside the boundary of the straight line to simulate the geological structure of the anticline or syncline.

[0123] In this embodiment, the anticline structure is located in the XOZ plane, parallel to the Y direction. The dip angles on both sides of the anticline are 45°, and the turning end is at the horizontal center of the area, that is, 45 / 2=22.5m. Figure 2 In (a) and (c), (a) is the initial background grid body, and (c) is a schematic diagram of the initial background grid body after spatial stretching. The model formation process, i.e., other steps, can refer to Example 1 and will not be described in detail.

[0124] After processing in steps S1-S8, the spatial variability of fractured rock mass under the three-dimensional anticline geological structure is realized as follows: Figure 6 As shown in (a)-(c) in the figure. Figure 6 (a) shows the spatial variability of the rock units in the anticline structure; Figure 6 (b) shows the spatial variability of joint units containing fractures under the anticline structure; Figure 6(c) is a cross-section diagram of the geological structure; Figure 6 (d) in the figure is the real geological structure map. Figure 6 As can be seen from the comparison between (a)-(c) and (d), the present invention can generate a fractured rock model with spatial variability under a three-dimensional anticline structure that is close to the real one.

[0125] Embodiment 3:

[0126] In this embodiment, for a fractured rock mass model with spatial variability under a fold structure, a fractured rock mass random field considering spatial variability under a three-dimensional geological structure is calculated.

[0127] Using the case model of Example 1, Figure 4 As shown, the geological structure of the area is changed to fold structure. The sizes of the rock mass area are still 45m, 45m and 30m respectively. The rock mass unit is tetrahedral in shape, with a size between 1-2m and a number of 85339, while the number of joint units is 337514. To ensure that the background area can contain the rock mass area no matter how it is transformed, its area size must be twice that of the rock mass. Let the width D of the three-dimensional background area D be x =90m, length D y =90m and height D z = 90m, the background grid unit size is set between 1m, that is, x i -x i-1 =1m,y j -y j-1 =1m,z k -z k-1 =1m, the maximum number is N=N x *N y *N z =91×91×91=753571. The mean value of rock mass tensile strength is also expressed in μ X =100.0MPa, standard deviationσ X =10.0MPa, related length L x =L y =45m,L z =2m; the autocorrelation function is also exponential, with a lower limit of c=0.1.

[0128] On any boundary surface of the background grid, distortion is performed along the normal direction of the boundary to simulate the folded geological structure.

[0129] In this embodiment, the fold structure is located in the XOZ plane, parallel to the Y direction. The difference between the peak and the valley of the fold is 10m, and it consists of two peaks and valleys. Figure 2(a) and (d), where (a) is the initial background mesh, and (b) is a schematic diagram of the initial background mesh after spatial distortion. The model formation process, i.e., other steps, can refer to Example 1 and will not be described in detail.

[0130] After processing in steps S1-S8, the realization effect of spatial variability of fractured rock mass under three-dimensional fold geological structure is shown in Figure Figure 7 As shown in (a)-(c) in the figure. Figure 7 (a) shows the spatial variability of rock units in fold structures; Figure 7 (b) shows the spatial variability of joint units containing fractures under the fold structure; Figure 7 (c) is a cross-section diagram of the geological structure; Figure 7 (d) in the figure is the real geological structure map. Figure 7 As can be seen from the comparison between (a)-(c) and (d), the present invention can generate a fractured rock model with spatial variability under a three-dimensional fold structure that is close to the real one.

[0131] The above descriptions are merely embodiments of the present invention and are not intended to limit the patent scope of the present invention. Any equivalent structure or equivalent process transformation made using the contents of the present invention specification and drawings, or directly or indirectly applied in other related technical fields, are also included in the patent protection scope of the present invention.

Claims

1. A method for realizing a random field of a fractured rock mass considering spatial variability under three-dimensional geological structure, characterized in that: The following steps are involved: S1: Establish a three-dimensional rectangular coordinate system, set the length, width and height of the background area according to the size of the rock mass area, ensure that the background area can completely cover the rock mass area, and determine the number of grid points of the background area in three dimensions to generate a three-dimensional background grid volume; S2: Set the mathematical statistical characteristic variables in the random field, including the mean, standard deviation, autocorrelation function and correlation length of the rock mass mechanical parameters; S3: Construct the covariance matrix, decompose it using the Karhunen-Loeve series expansion method, approximate the eigenvalues ​​and eigenvectors, and then obtain the random field values ​​at each node in the three-dimensional background grid; S4: Based on the geological structure type of the rock mass area, a spatial transformation operation is performed on the background grid, and the transformed grid is moved to the center of the rock mass area to ensure that the rock mass area is located inside the transformed background grid; S5: For any rock mass unit point, find the node of the background grid body corresponding to it, and calculate the random field value of the rock mass unit point by using the three-dimensional eight-node shape function interpolation method according to the random field values ​​of the eight background grid body nodes around the background grid body node; S6: Perform mean regression and variance regression on the random field of rock mass unit points, eliminate the errors of the Karhunen-Loeve series expansion method and the shape function interpolation method on the mean and variance, and assign them to the mechanical parameters of the rock mass unit; S7: inserting joint units between rock mass units and determining relevant mechanical parameters of the joint units; S8: Set the threshold and joint weakening coefficient, adjust the mechanical parameters of the joint according to the ratio of the random field values ​​of the rock units on both sides of the joint unit, and effectively distinguish between joints and cracks.

2. The method for realizing a random field of a fractured rock mass considering spatial variability under a three-dimensional geological structure according to claim 1, characterized in that: In step S1, the width D of the constructed three-dimensional background area D x Length D y and height D z The sizes of the generated background mesh body are larger than the corresponding sizes of the rock mass area, ensuring that after the background mesh body is deformed and translated in step S4, the background mesh body can still completely cover the entire rock mass area. According to the number of rock mass units in the three-dimensional direction of the rock mass area, the width D of the three-dimensional background area D is set. x The number of nodes N x Length D y The number of nodes N y , height D z The number of nodes N z , ensuring that the spacing between nodes in the background grid is smaller than the rock mass unit size, the total number of grid nodes N in the background area D is expressed as N = N x *N y *N z .

3. The method for realizing a random field of a fractured rock mass considering spatial variability under a three-dimensional geological structure according to claim 2, characterized in that: In the step S2, Set any two points A(x i ,y j ,z k )、B(x l ,y m ,z n ) are X i,j,k and X l,m,n , 1≤i≤N x ,1≤j≤N y ,1≤k≤N z , 1≤l≤N x ,1≤m≤N y ,1≤n≤N z ; mean μ X is the average value of the mechanical parameters of all rock mass units in the background area D, and the standard deviation σ X represents the mechanical parameters X at each grid node in the background area D i,j,k Distance from mean μ X the degree of deviation; The autocorrelation function is a function of the relative distance (Δx, Δy, Δz) between any two points A and B in the background area D. The rock mass mechanical parameters X i,j,k and X l,m,n The correlation will only decrease with the increase of relative distance and will approach a non-negative constant c, 0≤c<1; The three-dimensional autocorrelation function ρ(X i,j,k ,X l,m,n ) is exponential and square exponential, and their expressions are shown in formulas (2.1) and (2.2) respectively: Where Δx i,l =|x i -x l |,Δy j,m =|y j -y m |,Δz k,n =|z k -z n |;L x , L y and L z Respectively represent the relevant lengths of the rock mass mechanical parameters in the three dimensions, L x , L y and L z It is obtained by detecting the correlation of various directions in the rock mass area, or is artificially specified based on engineering experience.

4. The method for realizing a random field of a fractured rock mass considering spatial variability under a three-dimensional geological structure according to claim 3, characterized in that: In step S3, a covariance matrix is ​​constructed, and the Karhunen-Loeve series expansion method is used to decompose it, and the eigenvalues ​​and eigenvectors are approximately obtained to obtain the random field values ​​at each node in the background grid. The calculation formula is shown in (3.1): Where X is the mechanical parameter X of the grid point in the background area D i,j,k The collection of X = {X i,j,k |1≤i≤N x ,1≤j≤N y ,1≤k≤N z }; RF(X) is the set of random field values ​​generated by the set X; μ X is the mean value, σ X is the standard deviation; p is the covariance matrix ρ NN The pth expanded term of the decomposition; χ is the total number of background mesh nodes N = N x *N y *N z , and the random numbers satisfying the standard normal distribution; p is NN The eigenvalue set λ ijk After sorting in descending order, the eigenvalue of p is numbered, and φ p (X) corresponds to λ p The eigenvector of ijk and φ(X i,j,k ) satisfies formula (3.2): ∫ D p(X i,j,k ,X l,m,n )φ(X l,m,n )d(X l,m,n )=λ ijk φ(X i,j,k (3.2) This formula uses the covariance matrix ρ NN As a kernel function, it is defined as ρ(X i,j,k ,X l,m,n ), N=N x *N y *N z ρ NN is the autocorrelation function value between any two nodes of the background grid, and its specific form is shown in formula (3.3): Among them, ρ(X i,j,k ,X l,m,n The specific expression of ) depends on the type of autocorrelation function, refer to formula (2.1) or (2.2); In addition, φ(X i,j,k ) are mutually orthogonal, and their integral with themselves is 1, as shown in formulas (3.4) and (3.5): ∫ D φ(X i,j,k )φ(X l,m,n )dV=δ(X i,j,k ,X l,m,n (3.4) 5. The method for realizing a random field of a fractured rock mass considering spatial variability under a three-dimensional geological structure according to claim 4, characterized in that: In step S4, the spatial transformation of the background grid body includes spatial rotation, spatial stretching or spatial distortion; The geological structure morphology is set by simulating the deformation of the background grid. The specific operations are as follows: (a) Rotate the background grid in any direction to simulate the geological structure of parallel rock layers; (b) selecting any two nodes on the boundary surface of the background grid body to determine a straight line, and stretching the background grid body along the normal direction outside the boundary of the straight line to simulate the geological structure of the anticline or syncline; (c) On any boundary surface of the background grid, distortion is performed along the normal direction of the boundary to simulate the fold geological structure.

6. The method for realizing a random field of a fractured rock mass considering spatial variability under a three-dimensional geological structure according to claim 5, characterized in that: The step S5 specifically includes the following steps: S5.1: For any rock mass unit point Y t (x, t ,y t ,,z, t ), find the unique corresponding background mesh node (x i ,y j ,z k ), so that it satisfies: x i-1 <x, t <x i ,y j-1 <y t , <y j ,z k-1 <z, t <z k , t=1,2...s, s represents the number of rock mass unit points in the rock mass area; S5.2: Replace {(x a ,y b ,z c )|i-1≤a≤i,j-1≤b≤j,k-1≤c≤k} are used as the parent grid unit points, and the random field values ​​of the corresponding rock mass unit points are calculated using the three-dimensional eight-node shape function interpolation method. The calculation formula is as follows (5.1): Among them, the shape function T abc (ε a ε,η b η,ζ c ζ)=(1 / 8)*(1+ε a ε)(1+η b η)(1+ζ c ζ), ε i-1 =η j-1 =ζ k-1 =-1,ε i =η j =ζ k =1; (ε,η,ζ) is about (x, t ,y t ,,z, t ), the relationship between the two is (ε,η,ζ)=[J] -1 (x, t ,y t ,,z, t ), J is the Jacobian determinant, which is calculated by formula (5.2):

7. The method for realizing a random field of a fractured rock mass considering spatial variability under a three-dimensional geological structure according to claim 6, characterized in that: In step S6, the error caused by the Karhunen-Loeve series expansion method is eliminated by mean regression and standard deviation regression. The specific processing steps are as follows: S6.1: In mean regression, the actual mean μ of the rock mass regional random field RF(Y) is calculated according to formula (6.1): real ; Then, according to formula (6.2), each element of RF(Y) t ) multiplied by the set mean μ X With μ real The ratio of , we get the random field of the rock mass unit point after regression: Among them, RF(Y t ) new RF(Y k ) After mean regression processing, at the rock mass unit point Y t (x, t ,y t ,,z, t ) at the updated value; RF(Y) new = {RF(Y t ) new |t=1,2...s}, is RF(Y t ) new The set of is the rock mass regional random field after mean regression processing; μ X is the mean value of rock mass parameters, obtained through actual detection and calculation of rock mass; S6.2: In standard deviation regression, RF(Y) is calculated according to formula (6.3) new The actual standard deviation σ real ; Then, according to formula (6.4), calculate RF(Y t ) new With μ X The difference Dis_RF(Y t ); finally, according to formula (6.5), the random field value of each rock mass unit point after standard deviation regression processing is calculated; Dis_R F(Y t )=RF(Y t ) new -μ X (6.4) Among them, RF(Y t ) final is RF(Y t ) new After standard deviation regression processing, at the rock mass unit point Y t (x, t ,y t ,,z, t ) Update value; RF(Y) final = {RF(Y t ) final |t=1,2...s}, is RF(Y t ) final The set of is the final value of the random field in the rock mass area; σ X is the standard deviation of rock mass parameters, which is obtained through actual detection and calculation of the rock mass.

8. The method for realizing a random field of a fractured rock mass considering spatial variability under a three-dimensional geological structure according to claim 7, characterized in that: The step S7 specifically includes the following steps: S7.1: For any two rock mass units with a common surface, generate joint units between them that are the same as the common surface, so that the rock mass units become discontinuous with each other; S7.2: Set the relevant mechanical parameters Z of the joint element g Equal to the mean value μ of the corresponding rock mass mechanical parameters X , that is, Z={Z g =μ X |g=1,2,…G}, where Z g Represents the mechanical parameters of the gth joint unit, and G is the number of joint units.

9. The method for realizing a random field of a fractured rock mass considering spatial variability under a three-dimensional geological structure according to claim 8, characterized in that: The step S8 specifically includes the following steps: S8.1: Based on the standard deviation σ of rock mass mechanical parameters X , set the threshold According to the mechanical parameters of the rock mass and the mechanical parameters of the cracks, the joint weakening coefficient W is determined; S8.2: Calculate the ratio of the random field values ​​of the rock mass units on both sides of the joint unit and make a judgment: If the ratio is less than the threshold R Y , then multiply the joint mechanical parameters by the average value of the random field values ​​of the rock mass units on both sides; If it is greater than or equal to the threshold R Y , then multiply the joint mechanical parameters by the average value of the random field values ​​of the rock mass units on both sides, and then multiply them by the joint weakening coefficient to weaken the joint mechanical parameters to the crack level.

Images