Method for evaluating anisotropic roughness of rock fracture surface

Abstract

The invention belongs to the technical field of rock mechanics and engineering geology, and particularly relates to a method for evaluating anisotropic roughness of a rock fracture surface. A fracture surface is divided into dynamic infinitesimal grids matched with the rock grain size, the effective projection area of each infinitesimal in the specified analysis direction is calculated, the ratio of the sum of the projection areas to the area of a reference plane is counted, and the direction dependence of roughness is quantified. According to the method, by means of the double-cosine correction model, the relation between the local morphology and the mechanical response is directly established, the steep surface features of the damage sensitive area can be accurately captured, and the spatial distribution rule of the roughness can be visually displayed through the polar coordinate anisotropy rose diagram.

Description

Technical Field

[0001] The invention belongs to the technical field of rock mechanics and engineering geology, and in particular relates to a method for evaluating the anisotropic roughness of a rock fracture surface. Background Art

[0002] In the field of rock mechanics and engineering geology, the roughness of rock fracture surfaces is a core parameter that affects the mechanical behavior of rock masses (such as shear strength and seepage characteristics) and engineering stability. However, traditional evaluation methods (such as the JRC empirical method and the root mean square height method) are usually based on two-dimensional contour line analysis, and can only obtain a single scalar parameter through statistical height fluctuations or empirical comparisons, making it difficult to fully characterize the spatial heterogeneity of three-dimensional fracture surfaces.

[0003] Laboratory studies have shown that damage evolution during rock fracture exhibits significant locality and directional selectivity: only about 10%-30% of the surface area is involved in damage, and the damaged area is mainly concentrated on the rough surface with a local inclination direction opposite to the crack propagation direction, especially in steep areas with an inclination angle greater than 30°. This phenomenon shows that the directional distribution characteristics of roughness have a decisive influence on the mechanical behavior of rock mass. However, due to the lack of directional sensitivity and local morphology analysis capabilities of traditional methods, it is difficult to achieve a quantitative correlation between roughness parameters and engineering responses, thereby limiting the accuracy of rock stability prediction.

[0004] In recent years, with the breakthrough of 3D laser scanning technology, it has become possible to obtain 3D topographic data of rock fracture surfaces with sub-millimeter accuracy. Based on point cloud data, researchers have proposed 3D parameterization methods such as Fourier spectrum analysis and surface fractal dimension, but these methods have the following two limitations: On the one hand, the computational complexity is high and requires reliance on high-performance computing equipment, which makes it difficult to meet the needs of real-time engineering analysis; on the other hand, these methods are mainly based on global statistics and lack the ability to effectively analyze the anisotropic characteristics of roughness. Summary of the invention

[0005] The purpose of the present invention is to provide a method for evaluating the anisotropic roughness of rock fracture surfaces, which can not only accurately capture the steep surface characteristics of damage-sensitive areas, but also intuitively display the spatial distribution law of roughness through polar coordinate anisotropic rose diagrams.

[0006] The technical solution for achieving the purpose of the present invention is: a method for evaluating the anisotropic roughness of a rock fracture surface, dividing the fracture surface into a dynamic micro-element grid that matches the rock grain scale, and quantifying the directional dependence of the roughness by calculating the effective projected area of ​​each micro-element in a specified analysis direction and calculating the ratio of the total projected area to the area of ​​the reference plane.

[0007] Further, the specific steps include:

[0008] Step (1): pretreatment of the fracture surface of the rock specimen;

[0009] Step (2): Obtain a point cloud dataset P of the fracture surface of the specimen, and establish the point cloud dataset in a Cartesian coordinate system;

[0010] Step (3): De-noise and interpolate the point cloud data;

[0011] Step (4): Divide the rock fracture surface into a continuous rectangular grid structure consisting of a plurality of rectangular micro-elements with a side length of L;

[0012] Step (5): Perform optimal plane fitting on the continuous rectangular grid structure to establish a unified global reference plane Π;

[0013] Step (6): Establish the local plane of a single rectangular microelement, namely the rectangular microelement plane Π mn ;

[0014] Step (7): Calculate the plane Π of each rectangular microelement mn Azimuth angle relative to the specified analysis direction α ;

[0015] Step (8): Calculate the plane Π of each rectangular microelement mn The inclination angle θ relative to the global reference plane Π mn ;

[0016] Step (9): For each rectangular microelement plane Π mn , calculate its effective projected area on the global reference plane Π , and add up all the projected areas to get the sum of the effective areas of the rectangular microelement ;

[0017] Step (10): Define the total effective area Area of ​​the global reference plane The ratio is the evaluation index R(α) of the roughness of the fracture surface of the specimen along the specified analysis direction α;

[0018] Step (11): Calculate the anisotropy index AR of the fracture surface roughness.

[0019] Furthermore, the pretreatment in step (1) includes: spraying the image acquisition enhancer uniformly on the fracture surface of the rock after failure.

[0020] Furthermore, in step (2), a three-dimensional point cloud dataset P of the fracture surface of the specimen is obtained using a three-dimensional laser acquisition system, and the point cloud dataset is established in a Cartesian coordinate system with the crack propagation direction as the x-axis; wherein the three-dimensional point cloud dataset P of the fracture surface is expressed as:

[0021] ,

[0022] Where R represents a real number set, and i represents the number of each point cloud.

[0023] Furthermore, step (3) is as follows: calculate each point cloud data z i Coordinate mean μ z and standard deviation σ z , denoise the point cloud data; the denoising condition obeys the normal distribution, so that 99.7% of the data points fall within the range of mean μ±3σ:

[0024] ,

[0025] The moving least squares method is used to interpolate and complete the point cloud data to ensure the continuity of the point cloud dataset.

[0026] Furthermore, step (4) is specifically as follows: the value of the side length L of the rectangular infinitesimal element is determined according to the following formula:

[0027] L=k*d grain

[0028] where d grain is the average grain diameter of the specimen, k is the magnification factor, and is a constant ranging from 0 to 1.

[0029] Furthermore, step (5) is specifically as follows: the global reference plane Π is expressed as:

[0030] ,

[0031] In the formula, the parameters and is the plane equation parameter, which determines the tilt direction of the plane; is the constant term of the plane equation, which represents the intercept of the plane on the z-axis; the three parameters satisfy the condition of minimizing the sum of squared errors:

[0032] ,

[0033] The normalized normal vector n0 of the global reference plane Π is:

[0034] ;

[0035] Step (6) is as follows: use the least squares method to calculate the point cloud data set P in the rectangular microelement mn Fitting is performed to obtain the local plane Π mn , where m and n represent the number of point cloud data in the rectangular microelement; local plane Π m It is expressed as:

[0036] ,

[0037] In the formula, parameter a mn With b mn is the plane equation control parameter, c mn is the constant term of the plane equation;

[0038] Normalized normal vector n mn for:

[0039] .

[0040] Furthermore, step (7) transforms the analysis direction α into a three-dimensional unit vector α(cosα, sinα, 0), and the normal vector n mn The angle between the vector α and ∈[0,π / 2] The calculation formula is:

[0041]

[0042] Azimuth The calculation formula is:

[0043] ;

[0044] In step (8), the inclination angle θ mn The calculation formula is as follows:

[0045] .

[0046] Furthermore, step (9) is specifically as follows: the effective area calculation formula corresponding to a single rectangular microelement is:

[0047] ,

[0048] In the formula, Represents the area of ​​a single microelement rectangle;

[0049] The sum of the effective areas of the entire rectangular element is:

[0050] .

[0051] Furthermore, the area calculation formula of the global reference plane Π in step (10) is as follows:

[0052] ,

[0053] In the formula, M and U represent the length and width of the fracture surface of the specimen;

[0054] The calculation formula of the evaluation index R(α) of the specimen fracture surface roughness is as follows:

[0055] ,

[0056] When R(α)∈(0,1), the larger the value, the lower the roughness of the specimen fracture surface; when R(α)=1, it indicates that the specimen is completely smooth in the analysis direction α;

[0057] The calculation formula of the anisotropy index AR in step (11) is:

[0058] ,

[0059] AR=1 means that the roughness of the fracture surface is completely isotropic; AR>1: there is anisotropy, and the larger the value, the more significant the directional difference.

[0060] Compared with the prior art, the present invention has the following significant advantages:

[0061] The present invention divides the fracture surface into a dynamic micro-element grid that matches the rock grain scale. The effective projection area of ​​each micro-element in the specified analysis direction is calculated, and the ratio of the total projection area to the reference plane area is statistically analyzed to quantify the directional dependence of the roughness. Compared with traditional methods, this technology uses a double cosine correction model (taking into account the inclination compression effect and directional matching) to directly establish the connection between local morphology and mechanical response. It can not only accurately capture the steep surface characteristics of damage-sensitive areas, but also intuitively display the spatial distribution law of roughness through the polar coordinate anisotropic rose diagram. BRIEF DESCRIPTION OF THE DRAWINGS

[0062] Figure 1 The present invention is a flow chart of the method for evaluating the anisotropic roughness of rock fracture surfaces.

[0063] Figure 2 This is a diagram of spraying an image acquisition enhancer on a rock fracture surface according to the present invention.

[0064] Figure 3 This is the distribution diagram of the coordinate system of the rock fracture surface of the present invention.

[0065] Figure 4 This is a graph of the missing region completed using the moving least squares method of the present invention.

[0066] Figure 5 This is a diagram of the rectangular microelement network structure of the present invention.

[0067] Figure 6 It is a global reference plane map established based on the rock fracture surface in the present invention.

[0068] Figure 7 A schematic diagram of calculating the azimuth angle of a micro-element rectangle relative to a specified analysis direction in the present invention.

[0069] Figure 8 It is a schematic diagram of calculating the inclination angle of the micro-element rectangle relative to the global reference plane diagram according to the present invention.

[0070] Fig. 9 Schematic diagram of the spatial distribution law of fracture surface roughness of the present invention. DETAILED DESCRIPTION

[0071] The present invention is further described in detail below in conjunction with the accompanying drawings.

[0072] A method for evaluating the anisotropic roughness of a rock fracture surface comprises the following steps:

[0073] Step (1): Pre-treat the fracture surface of the rock specimen by evenly spraying an image acquisition enhancer on the fracture surface of the rock after failure.

[0074] Step (2): Use a three-dimensional laser acquisition system to obtain a three-dimensional point cloud dataset P of the fracture surface of the specimen, and establish the point cloud dataset in a Cartesian coordinate system with the crack propagation direction as the x-axis.

[0075] The fracture surface 3D point cloud dataset P is expressed as:

[0076] ,

[0077] Where R represents a real number set, and i represents the number of each point cloud.

[0078] Step (3): Calculate each point cloud data z i Coordinate mean μ z and standard deviation σ z , denoise the point cloud data.

[0079] The noise reduction condition follows a normal distribution, so that 99.7% of the data points fall within the range of the mean μ±3σ:

[0080] ,

[0081] Step (4): Use the moving least squares method to interpolate and complete the point cloud data to ensure the continuity of the point cloud data set.

[0082] Step (5): Divide the rock fracture surface into a continuous network structure consisting of rectangular elements with a side length of L.

[0083] The side length of the rectangular microelement is determined by the average grain diameter d grain Determined by the amplification factor k (k is a constant related to the number of rectangular differentials on the fracture surface, ranging from 0 to 1), it is expressed as:

[0084] L=k*d grain ,

[0085] Step (6): Perform the best plane fitting on the continuous rectangular grid structure by the least square method to establish a unified global reference plane Π.

[0086] The global reference plane Π is expressed as:

[0087] Z=a0x+b0y+c0,

[0088] In the formula, parameters a0 and b0 are the plane equation parameters, which determine the tilt direction of the plane; c0 is the plane equation constant term, which represents the intercept of the plane on the z-axis. These three parameters satisfy the condition of minimizing the sum of squared errors:

[0089] ,

[0090] The normalized normal vector n0 of the global reference plane Π is:

[0091] ,

[0092] Step (7): Use the least squares method to calculate the point cloud data set P within the rectangular element. mn (m and n represent the numbers of the point cloud data in the rectangular microelement) to fit and obtain the local plane equation Π mn , calculate the normalized normal vector n of each local plane mn .

[0093] The local plane fitting equation is:

[0094] ,

[0095] In the formula, parameter a mn With b mn is the plane equation control parameter, c mn is the constant term of the plane equation.

[0096] Normalized normal vector n mn for:

[0097] ,

[0098] Step (8): Calculate the plane Π of each rectangular microelement mn Azimuth relative to the specified analysis direction α∈[0,2π] ∈[0,π / 2].

[0099] Convert the analysis direction α into a three-dimensional unit vector α (cosα, sinα, 0), the normal vector n mn The angle between the normal and the vector α ∈[0,π / 2] The calculation equation is:

[0100] ,

[0101] Azimuth The calculation equation is:

[0102] ,

[0103] Step (9): Calculate the plane Π of each rectangular microelement mn The inclination angle θ relative to the global reference plane Π mn ∈[0,π / 2].

[0104] The inclination angle calculation equation is:

[0105] ,

[0106] Step (10): For each rectangular microelement plane Π mn , calculate its effective projected area on the global reference plane Π , and add up all the projected areas to get the sum of the effective areas of the rectangular microelement .

[0107] Global reference plane Π area The calculation equation is:

[0108] ,

[0109] Where M and U represent the length and width of the fracture surface of the specimen, and L is the side length of the rectangular element.

[0110] The calculation equation for the effective area corresponding to a single microelement is:

[0111] ,

[0112] In the formula, Represents the area of ​​a single microelement rectangle.

[0113] The sum of the effective areas of the entire rectangular element is:

[0114] ,

[0115] Step (11): Define the total effective area Area of ​​the global reference plane The ratio is the evaluation index R(α) of the roughness of the fracture surface of the specimen along the specified direction α.

[0116] The calculation equation for the evaluation index of fracture surface roughness is:

[0117] ,

[0118] When R(α)∈(0,1), the larger its value is, the lower the roughness of the fracture surface of the specimen is; when R(α)=1, it indicates that the specimen is completely smooth in the analysis direction α.

[0119] Step (12): Calculate the roughness index R(α) of the fracture surface of the specimen at each analysis angle with a calculation interval of 10°. The ratio AR of the maximum value to the minimum value of R(α) is defined as the anisotropy index of the fracture surface of the specimen.

[0120] The calculation equation for the anisotropy index is:

[0121] ,

[0122] AR=1 means that the roughness of the fracture surface is completely isotropic; AR>1: there is anisotropy, and the larger the value, the more significant the directional difference.

[0123] Example 1

[0124] The present invention is further described in detail by taking a granite specimen as an example in conjunction with the accompanying drawings.

[0125] like Figure 1 As shown, the present invention proposes a method for evaluating the anisotropic roughness of a rock fracture surface, comprising the following steps:

[0126] Step 1: Place a granite specimen with a diameter of 50 mm and a thickness of 25 mm in a hydraulic testing machine and generate a fracture surface by the Brazilian splitting method. Spray the nano-alumina reinforcement evenly on the fracture surface to form a uniform reflective layer (see Figure 2 ). After spraying, let it dry for 12 hours.

[0127] Step 2: Use a 3D laser scanning system to collect the fracture surface point cloud dataset P = {(x i ,y i ,z i )}, the scanning resolution is 0.05mm, and the scanning rate is 976000 points / second. The scanned point cloud data set is established in the Cartesian coordinate system with the crack propagation direction as the x-axis (see Figure 3 ).

[0128] Step 3: Calculate the average z coordinate μ in the point cloud dataset z =1.8mm and standard deviation σ z =0.08mm, remove the satisfied =0.24mm abnormal point.

[0129] Step 4: Use the moving least squares method to fill in the missing areas (see Figure 4 The red area in the figure). The kernel function bandwidth h=4d grain=2mm, the polynomial order is 2. The point cloud data density after interpolation is 30000 points / mm 2 .

[0130] Step 5: Measure the average grain diameter d of the granite specimen using a thin section microscope grain =0.5, and the magnification factor k=0.4. The fracture surface is divided into 12,000 rectangular micro-elements with a side length of 0.2 mm (see Figure 5 ), which assumes that the granite grain size distribution is uniform.

[0131] Step 6: Fit the divided rectangular micro-element network structure by the least squares method to obtain the plane equation of the global reference plane Π z = 0.008x − 0.012y + 1.8 (see Figure 6 ), calculate the normalized normal vector n0 of plane Π.

[0132] Step 7: Fit the point cloud data set within a single microelement rectangle using the least squares method to obtain the local plane fitting equation Π mn , and calculate its normalized normal vector n mn .

[0133] Step 8: Calculate the azimuth of the micro-rectangle in the analysis direction α (See Figure 7 ).

[0134] Step 9: Calculate the rectangular microelement plane Π mn The inclination angle θ relative to the global reference plane Π mn (See Figure 8 ).

[0135] Step 10: Calculate the area of ​​the global reference plane Π 480mm 2 , according to the azimuth of the rectangular element With the inclination angle θ mn Calculate the effective projected area of ​​a single microelement , and add them together to get the total effective projected area .

[0136] Step 11: Total Effective Area With global reference plane area The ratio is the evaluation index R(α) of the fracture surface roughness of the specimen along the specified direction α.

[0137] Step 12: Fig. 9The roughness evaluation index R(α) corresponding to each analysis angle of the granite fracture surface is shown. It can be found that the roughness along the crack extension direction (α=0°) is the lowest, with a value of 0.552; the roughness of the fracture surface perpendicular to the crack extension direction (α=90°) is the highest, with a value of 0.788. The anisotropy index of the roughness of the fracture surface is 1.428, which has strong roughness anisotropy.

Claims

1. A method for evaluating the anisotropic roughness of a rock fracture surface, characterized in that: The fracture surface is divided into a dynamic micro-element grid that matches the rock grain size. The directional dependence of the roughness is quantified by calculating the effective projected area of ​​each micro-element in the specified analysis direction and calculating the ratio of the total projected area to the area of ​​the reference plane.

2. The method according to claim 1, characterized in that The specific steps include: Step (1): pretreatment of the fracture surface of the rock specimen; Step (2): Obtain a point cloud dataset P of the fracture surface of the specimen, and establish the point cloud dataset in a Cartesian coordinate system; Step (3): De-noise and interpolate the point cloud data; Step (4): Divide the rock fracture surface into a continuous rectangular grid structure consisting of a plurality of rectangular micro-elements with a side length of L; Step (5): Perform optimal plane fitting on the continuous rectangular grid structure to establish a unified global reference plane Π; Step (6): Establish the local plane of a single rectangular microelement, namely the rectangular microelement plane Π mn ; Step (7): Calculate the plane Π of each rectangular microelement mn Azimuth angle relative to the specified analysis direction α ; Step (8): Calculate the plane Π of each rectangular microelement mn The inclination angle θ relative to the global reference plane Π mn ; Step (9): For each rectangular microelement plane Π mn , calculate its effective projected area on the global reference plane Π , and add up all the projected areas to get the sum of the effective areas of the rectangular microelement ; Step (10): Define the total effective area Area of ​​the global reference plane The ratio is the evaluation index R(α) of the roughness of the fracture surface of the specimen along the specified analysis direction α; Step (11): Calculate the anisotropy index AR of the fracture surface roughness.

3. The method according to claim 2, characterized in that Step (1) pretreatment includes: spraying the image acquisition enhancer uniformly on the fracture surface of the rock after failure.

4. The method according to claim 2, characterized in that: In step (2), a three-dimensional point cloud data set P of the fracture surface of the specimen is obtained using a three-dimensional laser acquisition system, and the point cloud data set is established in a Cartesian coordinate system with the crack propagation direction as the x-axis; The fracture surface 3D point cloud dataset P is expressed as: ， Where R represents a real number set, and i represents the number of each point cloud.

5. The method according to claim 1, characterized in that Step (3) is as follows: Calculate each point cloud data z i Coordinate mean μ z and standard deviation σ z , denoise the point cloud data; the denoising condition obeys the normal distribution, so that 99.7% of the data points fall within the range of mean μ±3σ: ， The moving least squares method is used to interpolate and complete the point cloud data to ensure the continuity of the point cloud dataset.

6. The method according to claim 1, characterized in that Step (4) is as follows: the value of the side length L of the rectangular element is determined according to the following formula: L=k*d grain where d grain is the average grain diameter of the specimen, k is the magnification factor, and is a constant ranging from 0 to 1.

7. The method according to claim 6, characterized in that Step (5) is as follows: The global reference plane Π is expressed as: ， In the formula, the parameters and is the plane equation parameter, which determines the tilt direction of the plane; is the constant term of the plane equation, which represents the intercept of the plane on the z-axis; the three parameters satisfy the condition of minimizing the sum of squared errors: ， The normalized normal vector n0 of the global reference plane Π is: ； Step (6) is as follows: use the least squares method to calculate the point cloud data set P in the rectangular microelement mn Fitting is performed to obtain the local plane Π mn , where m and n represent the number of point cloud data in the rectangular microelement; local plane Π m It is expressed as: ， In the formula, parameter a mn With b mn is the plane equation control parameter, c mn is the constant term of the plane equation; Normalized normal vector n mn for: 。 8. The method according to claim 7, characterized in that Step (7) converts the analysis direction α into a three-dimensional unit vector α (cosα, sinα, 0), and the normal vector n mn The angle between the vector α and ∈[0,π / 2] The calculation formula is: ， Azimuth The calculation formula is: ； In step (8), the inclination angle θ mn The calculation formula is as follows: 。 9. The method according to claim 8, characterized in that Step (9) is as follows: The effective area calculation formula corresponding to a single rectangular microelement is: ， In the formula, Represents the area of ​​a single microelement rectangle; The sum of the effective areas of the entire rectangular element is: 。 10. The method according to claim 9, characterized in that The area of ​​the global reference plane Π in step (10) The calculation formula is as follows: ， In the formula, M and U represent the length and width of the fracture surface of the specimen; The calculation formula of the evaluation index R(α) of the specimen fracture surface roughness is as follows: ， When R(α)∈(0,1), the larger the value, the lower the roughness of the specimen fracture surface; when R(α)=1, it indicates that the specimen is completely smooth in the analysis direction α; The calculation formula of the anisotropy index AR in step (11) is: ， AR=1 means that the roughness of the fracture surface is completely isotropic; AR>1: there is anisotropy, and the larger the value, the more significant the directional difference.

Images