A method for evaluating anisotropic roughness of rock fracture surfaces

Abstract

The present invention belongs to the field of rock mechanics and engineering geology technology, and specifically relates to a method for evaluating the anisotropic roughness of rock fracture surfaces. The fracture surface is divided into a dynamic micro-element grid that matches the rock grain scale. The effective projected area of ​​each micro-element in the specified analysis direction is calculated, and the ratio of the total projected area to the reference plane area is statistically analyzed to quantify the directional dependence of the roughness. With the help of a double cosine correction model, the present invention directly establishes 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 of roughness through polar coordinate anisotropic rose diagrams.

Description

Technical Field

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

[0002] In the fields of rock mechanics and engineering geology, the roughness of rock fracture surfaces is a key parameter that influences rock mechanical behavior (such as shear strength and seepage characteristics) and engineering stability. However, traditional assessment methods (such as the JRC empirical method and the root mean square height method) are typically based on two-dimensional contour line analysis. These methods can only derive 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 localization and directional selectivity: only approximately 10%-30% of the surface area is damaged, and the damaged area is primarily concentrated on rough surfaces with a local inclination opposite to the direction of crack propagation, particularly in steep areas with inclination angles greater than 30°. This phenomenon suggests that the directional distribution of roughness has a decisive influence on the mechanical behavior of rock masses. However, due to the lack of directional sensitivity and local morphology resolution capabilities of traditional methods, it is difficult to quantitatively correlate roughness parameters with engineering responses, thus limiting the accuracy of rock stability predictions.

[0004] In recent years, breakthroughs in 3D laser scanning technology have made it possible to obtain submillimeter-level 3D topographic data of rock fracture surfaces. Based on point cloud data, researchers have proposed 3D parameterization methods such as Fourier spectrum analysis and surface fractal dimension. However, these methods have two limitations: First, they are computationally complex and rely on high-performance computing equipment, making them difficult to meet the needs of real-time engineering analysis. Second, these methods rely primarily on global statistics and lack the ability to effectively analyze anisotropic 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 objectives of the present invention is: a method for evaluating the anisotropic roughness of rock fracture surfaces, which divides the fracture surface 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 a specified analysis direction and calculating the ratio of the total projected area to the area of ​​the reference plane.

[0007] Furthermore, the method specifically includes the following steps:

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

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

[0010] Step (3): Denoise and interpolate the point cloud data;

[0011] Step (4): Divide the rock fracture surface into a continuous rectangular grid structure consisting of multiple rectangular 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 element, namely the rectangular element plane Π mn ;

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

[0015] Step (8): Calculate the π of each rectangular element plane 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 elements ;

[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 an 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. The point cloud dataset is established in a Cartesian coordinate system with the crack propagation direction as the x-axis. The three-dimensional point cloud dataset P of the fracture surface is expressed as:

[0021] ,

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

[0023] Furthermore, step (3) is specifically 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 element. mn Fitting is performed to obtain the local plane Π mn , where m and n represent the number of point cloud data in the rectangular element; local plane Π m 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 Angle with vector α ∈[0,π / 2] is calculated as:

[0041]

[0042] Azimuth The calculation formula is:

[0043] ;

[0044] The inclination angle θ in step (8) 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 element is:

[0047] ,

[0048] Where, Represents the area of ​​a single infinitesimal rectangle;

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

[0050] .

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

[0052] ,

[0053] Where 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 indicates 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 size. By calculating the effective projected area of ​​each micro-element in the specified analysis direction and statistically analyzing the ratio of the total projected area to the area of ​​the reference plane, the directional dependence of the roughness is quantified. 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 polar coordinate anisotropic rose diagrams. BRIEF DESCRIPTION OF THE DRAWINGS

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

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

[0064] Figure 3 This is the coordinate distribution diagram 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 This is a schematic diagram of the present invention for calculating the azimuth angle of the micro-element rectangle relative to the specified analysis direction.

[0069] Figure 8 This is a schematic diagram of the present invention for calculating the inclination angle of a microelement rectangle relative to a global reference plane.

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

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

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

[0073] Step (1): Pre-treating 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 set of real numbers, 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 dataset.

[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 element is determined by the average grain diameter d of the specimen. 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 that 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 dataset P within the rectangular element. mn (m and n represent the numbers of the point cloud data in the rectangular element) and fit to 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 element mn Azimuth angle relative to the specified analysis direction α∈[0,2π] ∈[0,π / 2].

[0099] The analysis direction α is converted 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 element mn The inclination angle θ relative to the global reference plane π mn ∈[0,π / 2].

[0104] The equation for calculating the inclination angle 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 elements .

[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 infinitesimal element is:

[0111] ,

[0112] Where, Represents the area of ​​a single infinitesimal rectangle.

[0113] The sum of the effective area 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 specimen fracture surface 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 indicates 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 will be further described in detail with reference to the accompanying drawings, taking a granite specimen as an example.

[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 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 area (see Figure 4 The red area in the figure). Wherein, the kernel function bandwidth h=4d grain=2mm, the polynomial order is 2. The point cloud data density after interpolation is 30,000 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 microelements with a side length of 0.2 mm (see Figure 5 ), which assumes that the granite has a uniform particle size distribution.

[0131] Step 6: Fit the divided rectangular micro-element network structure by the least square 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 infinitesimal 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 angle of the infinitesimal rectangle in the analysis direction α (See Figure 7 ).

[0134] Step 9: Calculate the rectangular element 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 and the inclination angle θ mn Calculate the effective projected area of ​​a single infinitesimal element , and add them together to get the total effective projected area .

[0136] Step 11: Total Effective Area and the 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: Figure 9The roughness evaluation index R(α) for each analysis angle of the granite fracture surface is displayed. It is found that the roughness is lowest along the crack propagation direction (α = 0°), with a value of 0.552; the roughness is highest perpendicular to the crack propagation direction (α = 90°), with a value of 0.788. The anisotropy index of the fracture surface roughness is 1.428, indicating 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. The specific steps include: Step (1): pretreatment of the fracture surface of the rock specimen; Step (2): Obtain the point cloud dataset P of the specimen fracture surface and establish the point cloud dataset in a Cartesian coordinate system; Step (3): Denoise and interpolate the point cloud data, specifically: 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σ: |z i -m z |>3s z , The moving least squares method is used to interpolate and complete the point cloud data to ensure the continuity of the point cloud dataset; Step (4): Divide the rock fracture surface into a continuous rectangular grid structure consisting of multiple rectangular elements with a side length of L; The value of the side length L of the rectangular infinitesimal 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; 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 element, namely the rectangular element plane Π mn ; Step (7): Calculate the plane π of each rectangular element mn Azimuth angle relative to the specified analysis direction α Step (8): Calculate the plane π of each rectangular element 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 elements Specifically, the calculation formula for the effective area corresponding to a single rectangular element is: Where S (mn) (α) represents the area of ​​a single infinitesimal rectangle; The sum of the effective area of ​​the entire rectangular element is: Step (10): Define the total effective area and the global reference plane π area S base 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.

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

3. The method according to claim 1, characterized in that In step (2), a three-dimensional laser acquisition system is used to obtain a three-dimensional point cloud dataset P of the fracture surface of the specimen, and the point cloud dataset 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: P={(x i ,y i ,z i )∈R 3 ∣i=1,2,...,N}, Where R represents a set of real numbers, and i represents the number of each point cloud.

4. The method according to claim 3, characterized in that Step (5) is specifically as follows: the global reference plane Π is expressed as: Z=a0x+b0y+c0, 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: The normalized normal vector n0 of the global reference plane π is: Step (6) is as follows: the point cloud dataset P in the rectangular element is subjected to the least square method. mn Fitting is performed to obtain the local plane Π mn , where m and n represent the number of point cloud data in the rectangular element; local plane Π m Expressed as: Π mn :z=a mn x+b mn y+c mn , 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:

5. The method according to claim 4, characterized in that Step (7) converts the analysis direction α into a three-dimensional unit vector α(cosα, sinα, 0), and the normal vector n mn Angle with vector α The calculation formula is: Azimuth The calculation formula is: The inclination angle θ in step (8) mn The calculation formula is as follows:

6. The method according to claim 5, characterized in that The area S of the global reference plane Π in step (10) base The calculation formula is as follows: S base =M·U·L 2 , Where 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 indicates that the roughness of the fracture surface is completely isotropic; AR > 1: anisotropy exists, and the larger the value, the more significant the directional difference.

Images