A random field realization method for fractured rock mass considering spatial variability under three-dimensional geological structure

Abstract

The present invention discloses a method for realizing a random field of a fractured rock mass taking into account spatial variability under three-dimensional geological structure. The method comprises the following steps: determining the size of a three-dimensional background grid area and the number of background nodes in each direction according to the rock mass area and the size of the rock mass unit; setting the value of a mathematical statistical characteristic variable in the three-dimensional random field; constructing a covariance matrix, decomposing it using the Karhunen‑Loeve series expansion method, and generating a random field approximation at the background grid point; performing spatial transformation operations and moving the background grid based on the geological structure type of the rock mass area; quickly selecting and judging the background grid point, and obtaining the random field value of the rock mass unit using shape function interpolation; adjusting the random field value by mean regression and standard deviation regression, and assigning it to the mechanical parameters of the rock mass; judging that when the random field ratio of the rock mass units on both sides of the joint unit is greater than a set threshold, weakening the mechanical properties of the joint as a fracture. The present invention can simulate the spatial variability of the rock mass under three-dimensional geological structure and can also realize the spatial distribution of fractures.

Description

Technical Field

[0001] The present 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 three-dimensional geological structure. Background Art

[0002] Natural rock masses are affected by geological movements and stress history during their formation, and exhibit more complex spatial variability than soil masses. 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 within it. The spatial variability and fracture surfaces of rock masses 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 masses with naturally formed fractures, and many engineering examples have shown that rock failure is more dependent on fractures than intact rocks.

[0003] The main geological structural types of rock masses include parallel-dipping strata, synclines and anticlines, folds, and faults. These structural types exhibit a high correlation with the spatial variability of rock masses. Random fields, as the primary means of measuring the spatial variability of rock and soil masses, have attracted considerable attention and application from researchers. However, as three-dimensional spatial entities, most rock mass studies employ simplified two-dimensional solid models. This is exemplified by the Chinese invention patent entitled "A Method for Implementing Rotationally Anisotropic Non-Stationary Random Fields for Rock and Soil Parameters," patent application number: 202110362390.7. However, many researchers have overlooked the influence of geological structure when implementing three-dimensional random field methods, as exemplified by 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 rock mass stability analysis, cracks are also a factor that cannot be ignored. Current mainstream research methods typically introduce cracks during the modeling phase. For example, a Chinese patent application titled "A discrete crack network generation method based on an iterative function system" (patent application number: CN202211146540.1) provides an effective approach for crack modeling. In their study of the stability of counter-dip rock slopes, Yun Zheng et al. (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) employed a two-dimensional crack inclination model, modeling different crack inclination angles to analyze their impact on slope stability. However, this modeling approach is not only time-consuming but also simplifies the spatial distribution of cracks, resulting in discrepancies with the actual three-dimensional crack distribution. The crack modeling approach is not only time-consuming but also involves certain simplifications in the modeling process, resulting in discrepancies with 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 of complex geological structures, thereby better describing and predicting 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 under three-dimensional geological structure, so as to solve the technical difficulties such as the inability to consider the influence of geological structure on spatial variability under three-dimensional space and the distribution of fractures in three-dimensional space in 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 three-dimensional geological structure includes 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 in 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 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 value at each node in the 3D background grid;

[0012] S4: Based on the geological structure type of the rock mass area, perform a spatial transformation 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;

[0013] S5: For any rock mass unit point, find the node of the corresponding background grid body, and calculate the random field value of the rock mass unit point using the three-dimensional eight-node shape function interpolation method based on 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 to eliminate the errors of the Karhunen-Loeve series expansion method and 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 are both larger than twice the corresponding size of the rock mass area, ensuring that after the generated 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 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 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 at points A and B are 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 ) are 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.

[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 This is derived by detecting correlations in all directions across the rock mass region, a mature technology. Refer to the Chinese patent titled "A Non-Stationary Random Field Modeling Method for Cross-Correlation of Rock and Soil Parameters," application number: CN202211404802.X. Step 2 of that patent is described in detail and will not be repeated here.

[0027] Further optimization is performed in step S3, where the covariance matrix is ​​constructed and the random field values ​​are calculated using the Karhunen-Loeve series expansion method. This method is chosen mainly because the number of 3D background grid points is large, making the traditional Cholesky decomposition method unsuitable for large matrix decomposition. Compared to 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. A value too large will reduce the efficiency of the series expansion method. Refer to Section 3.6.3 of the doctoral thesis "Analysis of the risk of instability of layered rock caverns based on random finite differences". χ is the length of the background grid body, which is the total number of nodes N = N x *N y *N z , and the random numbers satisfying the standard normal distribution; ρ 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] ∫ D ρ(X i,j,k ,X l,m,n )φ(X l,m,n )d(X l,m,n )=λ ijk φ(X i,j,k ) (3.2)

[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]

[0035] Further optimization, in step S4, the spatial transformation of the background mesh includes spatial rotation, spatial stretching or spatial distortion;

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

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

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

[0039] (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.

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

[0041] S5.1: For any rock mass unit point Yt (x' t ,y't,z' t ), find the only 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.

[0042] 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 mother 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):

[0043]

[0044] 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 ) function, and 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):

[0045]

[0046] Further optimization, 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:

[0047] S6.1: In mean regression, calculate the actual mean μ of the rock mass regional random field RF(Y) 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 can get the random field of rock mass unit point after regression:

[0048]

[0049] 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, that is, the rock mass regional random field after mean regression processing; μ X is the mean value of rock mass parameters, which is obtained through actual detection and calculation of rock mass.

[0050] S6.2: In standard deviation regression, calculate RF(Y) according to formula (6.3) new The actual standard deviation σ real ; Then, calculate RF(Y according to formula (6.4) 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;

[0051]

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

[0053]

[0054] 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 rock mass.

[0055] 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.

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

[0057] S7.1: For any two rock mass elements with a common surface, generate joint elements between them that are identical to the common surface, so that the rock mass elements become discontinuous with each other.

[0058] S7.2: Set the relevant mechanical parameters Z of the joint element g Equal to the mean value of the corresponding rock mass mechanical parameter μ 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.

[0059] By inserting joint elements between rock mass elements, the rock mass region is transformed into a discontinuous model and the relevant mechanical parameters of the joint elements are determined. This method can effectively simulate the mechanical behavior of rock masses in a fractured state and is applicable to the finite element method (FDEM), discrete element method (DEM), and discontinuous analysis method (DDA).

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

[0061] S8.1: Based on the standard deviation σ of rock mass mechanical parameters X , set the threshold Determine the joint weakening coefficient W based on the mechanical parameters of the rock mass and the mechanical parameters of the cracks;

[0062] S8.2: Calculate the ratio of the random field values ​​of the rock mass units on both sides of the joint unit and determine:

[0063] 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;

[0064] If the ratio 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 cracks in the rock mass, and is generally taken as 0.1%-10%.

[0065] Adjusting the mechanical parameters of a joint based on the ratio of the random field values ​​of the rock units on either side of the joint is more realistic. Due to differences in stress history, the mechanical properties of the two sides of the rock mass may differ significantly, and when this difference is too large, cracks are more likely to form.

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

[0067] 1. Based on the prior application (202411201166X, title: A random field modeling method for spatial variability of rock parameters under complex geological structures), the present invention continues to study and improves 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, the spatial distribution of cracks is realized by introducing joint units, setting thresholds and other operations, avoiding crack modeling operations, and solving 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 existing random field technology.

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

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

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

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

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

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

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

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

[0076] To make the objectives, technical solutions, and advantages of the embodiments of the present invention more clear, the technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the embodiments described are only part of the embodiments of the present invention, not all of them. All other embodiments obtained by ordinary technicians in this field based on the embodiments of the present invention without making any creative efforts shall fall within the scope of protection of the present invention.

[0077] Example 1:

[0078] 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 considering spatial variability under a three-dimensional geological structure specifically includes the following steps: Figure 1 As shown:

[0079] 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 in the background area in three dimensions to generate a three-dimensional background grid.

[0080] In this embodiment, the length, width and height of the rock mass area are set to 45m, 45m and 30m respectively. The shape of the rock mass unit is a tetrahedron. The engineering modeler sets its size between 1-2m and the number of units is 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 unit modeling diagram.

[0081] 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. 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.

[0082] 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:

[0083] In this embodiment, the mean tensile strength of the rock mass in this area 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 the 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.

[0084] 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:

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

[0086]

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

[0088] In this embodiment, the number of nodes in the three-dimensional background grid is 753,571. The traditional Cholesky decomposition method is not suitable for large matrix decomposition. However, the Karhunen-Loeve series expansion method is significantly more efficient than the traditional Cholesky decomposition matrix for large matrix decomposition. Its calculation formula is as follows:

[0089]

[0090] 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 about 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 number of p expansion terms in the decomposition determines the accuracy of the random field approximation. Generally speaking, a value of 20 is sufficient. A value that is too large may reduce the efficiency of the series expansion method. Refer to Section 3.6.3 of the doctoral thesis "Risk Analysis of Instability of Layered Rock Caverns Based on Random Finite Differences". χ is the length of the total number of background grid points N = N x *N y *N z =753571, and the random number 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 ρ in 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).

[0091] S4: Based on the geological structure type of the rock mass area, perform a spatial transformation 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:

[0092] In this example, 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 a parallel rock formation geological structure: rotate at any 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.

[0093] 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 using the three-dimensional eight-node shape function interpolation method based on the random field values ​​of the eight background grid body nodes surrounding the background grid body node, which specifically includes the following steps:

[0094] 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 <t j ,z k-1 <z’ t <z k )'s background mesh node (x i ,y j ,z k ), and because (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 <t j ,z k-1 <z’ t <z k By excluding any background grid node that meets the requirements, you can quickly select the (x i ,y j ,z k ); where (x' t ,y' t ,z' t ) represents the mechanical parameter Y of the rock mass unit pointt , t=1,2...s, s=85339 is the number of rock mass units in the background grid area D.

[0095] 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 mother 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).

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

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

[0098] S6.1: In mean regression, calculate the actual mean μ of the rock mass regional random field RF(Y) 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 can get the random field of rock mass unit point after regression:

[0099]

[0100] 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 is the mean value of rock mass parameters, which is obtained through actual detection and calculation of rock mass.

[0101] S6.2: In standard deviation regression, calculate RF(Y) according to formula (6.3) new The actual standard deviation σ real ; Then, calculate RF(Y according to formula (6.4) t ) newWith μ 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;

[0102]

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

[0104]

[0105] 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 , satisfying the set mean and standard deviation. Furthermore, 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, which now has the spatial variability of the tensile strength of the rock mass unit under the geological structure.

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

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

[0108] S7.1: For any two rock mass units with a common surface, generate joint units between them that are identical to the common surface, so that the rock mass units become discontinuous with each other, e.g. Figure 3 As shown, Figure 3 (a) in the equation is a continuous model. Figure 3 (b) in the figure is a discontinuous model.

[0109] 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. g Equal to the mean value of the corresponding rock mass mechanical parameter μ X , that is, Z={Z g=100|g=1,2,…337514}, where Z g Represents the mechanical parameter value of the g-th joint unit.

[0110] 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.

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

[0112] S8.1: Based on the standard deviation of rock mass mechanical parameters of 10.0, which is 0.1 times the mean, and considering the density of rock mass cracks, set the threshold R of the ratio of random field values ​​of rock mass units. Y =1-0.1 2 =0.99; In addition, the strength of 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;

[0113] S8.2: Calculate the ratio of the random field values ​​of the rock mass units on both sides of the joint unit and determine:

[0114] 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.

[0115] After the above processing, the spatial variability of fractured rock mass under three-dimensional parallel stratigraphic structure is realized, as shown in Figure 2. Figure 5 As shown in (a)-(c) in the figure. Figure 5 (a) shows the spatial variability of rock units parallel to the rock layer structure; Figure 5 (b) is the spatial variability of the joint units containing fractures under 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 the real one.

[0116] Example 2:

[0117] 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.

[0118] Using the case model of Example 1, Figure 4As 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.

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

[0120] In this example, 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 point 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 mesh, and (c) is a schematic diagram of the initial background mesh after spatial stretching. The model formation process, i.e., other steps, can be referred to in Example 1 and will not be described in detail.

[0121] 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 of the geological structure;

[0122] 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.

[0123] Example 3:

[0124] 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.

[0125] 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. 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.

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

[0127] 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 2In (a) and (d), (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 be referred to in Example 1 and will not be described in detail.

[0128] After processing in steps S1-S8, the effect of spatial variability of fractured rock mass under three-dimensional fold geological structure is shown in Figure 2. 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 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.

[0129] 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's description 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 in 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 rock mass mechanical parameters; 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 value at each node in the 3D background grid; S4: Based on the geological structure type of the rock mass area, perform a spatial transformation 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; S5: For any rock mass unit point, find the node of the corresponding background grid body, and calculate the random field value of the rock mass unit point using the three-dimensional eight-node shape function interpolation method based on 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 to eliminate the errors of the Karhunen-Loeve series expansion method and 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 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 is x Length D y and height D z are larger than the corresponding size of the rock mass area, ensuring that the generated background mesh body can still completely cover the entire rock mass area after deformation and translation in step S4; 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 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 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 at points A and B are 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 ) are 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 Represents the relevant lengths of rock mass mechanical parameters in 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 three-dimensional geological structure according to claim 3, characterized in that: In step S3, a covariance matrix is ​​constructed and decomposed using the Karhunen-Loeve series expansion method to approximately obtain the eigenvalues ​​and eigenvectors, and the random field values ​​at each node in the background grid are obtained. 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 ), B=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 three-dimensional geological structure according to claim 4, characterized in that: In step S4, the spatial transformation of the background mesh includes spatial rotation, spatial stretching or spatial distortion; The geological structure morphology is set by simulating the deformation of the background mesh. The specific operations are as follows: (a) Rotate the background grid in any direction to simulate the geological structure of parallel rock layers; (b) Select any two nodes on the boundary surface of the background grid to determine a straight line, and stretch the background grid 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 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 only 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 mother 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 ) function, and 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 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, calculate the actual mean μ of the rock mass regional random field RF(Y) 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 can get the random field of 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' r ) at the updated value; RF(Y) new ={RF(Y t ) new |t=1,2...s},is RF(Y t ) new The set of, that 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, calculate RF(Y) according to formula (6.3) new The actual standard deviation σ real ; Then, calculate RF(Y according to formula (6.4) 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_RF(Y t )=RF(Y t ) new -μ X (6.4) 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 rock mass.

8. The method for realizing a random field of a fractured rock mass considering spatial variability under 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 identical to 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 parameter μ 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 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 Determine the joint weakening coefficient W based on the mechanical parameters of the rock mass and the mechanical parameters of the cracks; S8.2: Calculate the ratio of the random field values ​​of the rock mass units on both sides of the joint unit and determine: 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 the ratio 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 units on both sides, and then multiply them by the joint weakening coefficient to weaken the joint mechanical parameters to the crack level.

Images