Method for acquiring equivalence parameters of high-porosity reef limestone

Abstract

The invention discloses a method for acquiring equivalence parameters of high-porosity reef limestone, and relates to the field of geotechnical engineering.The method comprises the steps that S1, a three-dimensional matrix-pore binary matrix is constructed based on CT scanning data; s2, based on the constructed three-dimensional matrix-pore binary matrix, performing three-dimensional pore direction vectorization so as to describe geometric features of pores in a three-dimensional space; s3, based on the constructed three-dimensional matrix-pore binary matrix, extracting pore morphological parameters; s4, based on the constructed three-dimensional matrix-pore binary matrix, calculating a pore fractal dimension; s5, based on the constructed three-dimensional matrix-pore binary matrix, analyzing pore size distribution and connectivity; and S6, calculating an equivalent elastic parameter based on the constructed three-dimensional matrix-pore binary matrix. The equivalent mechanical parameters of the high-porosity reef limestone are obtained through CT scanning, pore characteristic statistical analysis and finite element homogenization methods, and a foundation is laid for ocean engineering foundation design, geological disaster assessment and the like.

Description

Technical Field

[0001] The present invention relates to the field of geotechnical engineering, and in particular to a method for obtaining equivalent parameters of high-porosity reef limestone. Background Art

[0002] Reef limestone is the remains of reef-building coral colonies, a unique rock and soil formed through geological processes. Its unique diagenetic mechanism results in a pore structure distinct from that of ordinary rocks, resulting in a complex structure. Numerous studies have shown that the porosity of reef limestone varies greatly. The book "Sedimentary Petrology" mentions that its porosity can range from 20% to 40%. Shallow reef limestone (SRL) exhibits good pore connectivity and reaches a high porosity of 55.3 ± 3.2%, while deep reef limestone (DRL) has a dense structure with a porosity of only 4.9 ± 1.6%. Due to the diverse porosity and large proportion of reef limestone, its pore structure has a significant impact on its structure and mechanical properties.

[0003] Traditional methods rely on laboratory mechanical testing to obtain macroscopic pore parameters of reef limestone. However, the resulting test results are highly discrete and fail to reflect the impact of microscopic pores on the structure. This presents significant limitations and cannot meet engineering design requirements. Furthermore, existing techniques often ignore scale effects. Laboratory sample sizes are typically centimeters, which cannot reflect the multi-scale pore coupling effects of meter-scale rock masses encountered in actual projects.

[0004] In practical applications, CT scanning technology combined with commercial software (often using global threshold segmentation methods, such as adaptive threshold segmentation) can capture pore structure, but generally lacks quantitative statistical analysis of pore anisotropy. While this can be achieved (e.g., using adaptive threshold segmentation), the lack of quantitative statistical analysis of pore anisotropy leads to inaccurate modeling. Specifically, this is manifested in the lack of anisotropy characterization and the inability to distinguish horizontal and vertical pore morphological differences. Furthermore, the morphological parameter simplification commonly employed in existing research often reduces pores to equivalent spheres, ignoring the complex topological characteristics of real pores. This results in significant errors in the stress concentration factor caused by the pore aspect ratio.

[0005] In summary, existing technologies have systematic deficiencies in obtaining equivalent mechanical parameters of reef limestone and accurately correlating the microscopic pore structure and macroscopic mechanical behavior of reef limestone, which will directly lead to significant errors in the prediction of mechanical parameters relied upon in engineering design, threatening the safety of the structure. Summary of the Invention

[0006] In response to the above-mentioned technical problems existing in the prior art, the purpose of the present invention is to provide a method for obtaining the equivalent mechanical parameters of high-porosity reef limestone. Through CT scanning, pore characteristic statistical analysis and finite element homogenization method, the equivalent mechanical parameters of high-porosity reef limestone are obtained, which is suitable for application scenarios such as marine engineering foundation design and geological disaster assessment.

[0007] To achieve the above-mentioned object, the present invention adopts the following technical solution: A method for obtaining equivalent parameters of high-porosity reef limestone, comprising the following steps: S1, constructing a three-dimensional matrix-pore binary matrix based on CT scan data; S2. Based on the constructed three-dimensional matrix-pore binary matrix, three-dimensional pore direction vectorization is performed to describe the geometric characteristics of the pores in three-dimensional space; S3, extracting pore morphology parameters based on the constructed three-dimensional matrix-pore binary matrix; S4. Calculate the pore fractal dimension based on the constructed three-dimensional matrix-pore binary matrix; S5. Analyze pore size distribution and connectivity based on the constructed three-dimensional matrix-pore binary matrix; S6. Calculate the equivalent elastic parameters based on the constructed three-dimensional matrix-pore binary matrix.

[0008] Furthermore, step S2 includes the following sub-steps: S21, calculating the three-dimensional spatial gradient field in the pore area; S22. constructing a structure tensor based on the obtained three-dimensional spatial gradient field of the pore region; S23, performing eigenvalue decomposition on the constructed structure tensor, and extracting the eigenvector corresponding to the maximum eigenvalue; S24. Calculate the pore extension direction based on the eigenvector corresponding to the extracted maximum eigenvalue.

[0009] Furthermore, in step S21, the three-dimensional spatial gradient field of the pore area is calculated according to the following formula: in, is the gradient operator, They represent the partial derivatives along the x, y, and z directions, respectively, and are used to calculate the gradient field in the pore region and capture the changes in the pore boundary.

[0010] Furthermore, in step S22, a structure tensor is constructed based on the obtained three-dimensional spatial gradient field of the pore region according to the following formula: in Represents a Gaussian filtering operation.

[0011] Furthermore, step S3 includes the following sub-steps: S31, using 26-neighborhood connected domain analysis to mark pore areas; S32, calculating the sphericity of each connected pore; S33, calculating the aspect ratio of connected pores based on principal axis analysis; S34. Calculate the surface-to-volume ratio of connected pores.

[0012] Furthermore, step S4 includes the following sub-steps: S41, determining an adaptive box size sequence; S42, for each adaptive box size S k Divide the 3D image into The box network uses zero padding to handle non-divisible boundaries and counts the number of boxes containing at least one pore voxel N(S K ); S43, establishing a double logarithmic coordinate relationship; S44. Calculate the fractal dimension by linear regression.

[0013] Furthermore, step S5 includes the following sub-steps: S51, calculating the pore three-dimensional Euclidean distance field; S52, extracting the distance value of the pore area and converting it into physical size; S53, generating a pore size distribution density histogram and a cumulative curve based on the distance values ​​of the extracted pore areas; S54, using the bwconncomp function to perform 26-domain connectivity analysis; S55. Calculate the maximum proportion of connected pores.

[0014] Furthermore, step S6 includes the following sub-steps: S61. Establish material property interpolation model; S62, assembling the global stiffness matrix based on the element stiffness matrix; S63, applying fixed boundary conditions; S64, solve the displacement field; S65. Extract equivalent elastic parameters.

[0015] Furthermore, in step S61, a material property interpolation model is established using the following formula: in is a binary matrix value, and p=3 is a penalty factor.

[0016] Furthermore, in step S63, the fixed boundary conditions imposed include the displacement of the specific degree of freedom of the constrained model, Fix the first three degrees of freedom (ux, uy, uz) and release the remaining degrees of freedom.

[0017] The technical solution employed in this invention provides a method for obtaining equivalent mechanical parameters for high-porosity reef limestone. This method addresses the technical flaw in existing rock mass pore characterization techniques, which largely ignore the spatial orientation characteristics of pore groups and simplify pores into isotropic spherical structures, leading to a severe underestimation of material mechanical anisotropy. By analyzing structural tensor characteristics and using a three-dimensional directional quantization model (azimuth angle φ, inclination angle θ), combined with rose diagram visualization technology, this method achieves a fully parametric characterization of the principal extension directions of pore groups, addressing the failure of traditional homogenization methods to predict deformation in rock masses containing oriented pores.

[0018] Conventional analysis using existing commercial image processing software (such as Avizo) lacks systematic quantitative statistical analysis of specific pore morphological parameters such as sphericity, aspect ratio, and surface-to-volume ratio. This results in the misinterpretation of non-spherical pores, especially high-aspect-ratio tubular pores, as spherical, leading to a serious underestimation or miscalculation of stress concentration effects around pores. By establishing a three-dimensional collaborative analysis system combining sphericity, aspect ratio, and surface-to-volume ratio, and revealing pore morphology spectrum characteristics through kernel density distribution statistics, this approach addresses the problem of traditional methods misjudging stress concentration around complex pores.

[0019] Traditional fractal dimension calculation methods, such as fixed-step box counting, suffer from a technical drawback: their accuracy plummets at high porosity levels (>30%), leading to a lack of correlation between microroughness and macroscopic mechanical response. This paper addresses the technical limitations of this approach by employing an adaptive power-scaling box counting algorithm (where the box size is dynamically adjusted by powers of 2) to accurately calculate the fractal dimension of the pore surface. This provides a key parameter foundation for establishing a quantitative correlation between fractal dimension and macroscopic modulus.

[0020] To address the technical flaw of traditional homogenization algorithms, which treat pores as featureless blank areas and ignore the combined influence of pore directionality, morphological topology, and surface fractal features on the elastic modulus, and thus severely underestimate the damage effect of microscopic pores, a mapping model from "pore characteristic matrix → dynamic penalty factor → equivalent elastic tensor" was established to achieve the transfer of microscopic pore parameters to macroscopic mechanical responses. BRIEF DESCRIPTION OF THE DRAWINGS

[0021] Figure 1 This is a flow chart of a method for obtaining equivalent parameters of high-porosity reef limestone shown in an embodiment of the present invention; Figure 2 1 is a distribution histogram of the angle between the pore extension direction and the z-axis generated by a method for obtaining equivalent parameters of high-porosity reef limestone shown in an embodiment of the present invention; Figure 3 This is a horizontal projection rose diagram of pore extension azimuth distribution generated by a method for obtaining equivalent parameters of high-porosity reef limestone shown in an embodiment of the present invention; Figure 4 This is a connected pore sphericity distribution diagram generated by a method for obtaining equivalent parameters of high-porosity reef limestone shown in an embodiment of the present invention; Figure 5 This is a connected pore aspect ratio distribution diagram generated by a method for obtaining equivalent parameters of high-porosity reef limestone shown in an embodiment of the present invention; Figure 6 This is a connected pore surface-to-volume ratio distribution diagram generated by a method for obtaining equivalent parameters of high-porosity reef limestone shown in an embodiment of the present invention; Figure 7 This is a fractal dimension fitting diagram and a double logarithmic coordinate diagram of the box counting method generated by a method for obtaining equivalent parameters of high-porosity reef limestone shown in an embodiment of the present invention; Figure 8 is a pore size distribution density histogram generated by a method for obtaining equivalent parameters of high-porosity reef limestone shown in an embodiment of the present invention; Figure 9 This is a cumulative pore size distribution curve generated by a method for obtaining equivalent parameters of high-porosity reef limestone shown in an embodiment of the present invention; Figure 10 This is a thermodynamic diagram of equivalent elastic modulus generated by a method for obtaining equivalent parameters of high-porosity reef limestone shown in an embodiment of the present invention. DETAILED DESCRIPTION

[0022] The present invention will be described in detail below with reference to the accompanying drawings and embodiments.

[0023] Example 1 Refer to the attached Figure 1 The embodiment of the present invention provides a method for obtaining equivalent parameters of high-porosity reef limestone, comprising the following steps: S1, constructing a three-dimensional matrix-pore binary matrix based on CT scan data; S2. Based on the constructed three-dimensional matrix-pore binary matrix, three-dimensional pore direction vectorization is performed to describe the geometric characteristics of the pores in three-dimensional space; S3, extracting pore morphology parameters based on the constructed three-dimensional matrix-pore binary matrix; S4. Calculate the pore fractal dimension based on the constructed three-dimensional matrix-pore binary matrix; S5. Analyze pore size distribution and connectivity based on the constructed three-dimensional matrix-pore binary matrix; S6. Calculate the equivalent elastic parameters based on the constructed three-dimensional matrix-pore binary matrix.

[0024] Step S2 includes the following sub-steps: S21. Calculate the three-dimensional spatial gradient field in the pore area according to the following formula.

[0025] in, is the gradient operator, They represent the partial derivatives along the x, y, and z directions, respectively, and are used to calculate the gradient field in the pore region and capture the changes in the pore boundary.

[0026] S22. Construct a structure tensor based on the obtained three-dimensional spatial gradient field of the pore area according to the following formula.

[0027] in Represents a Gaussian filtering operation.

[0028] S23. Perform eigenvalue decomposition on the constructed structure tensor , extract the eigenvector corresponding to the maximum eigenvalue .

[0029] in is the eigenvector matrix, storing the pore direction information, is the eigenvalue diagonal matrix, representing the importance weight of the direction, Indicates the eigenvector corresponding to the maximum eigenvalue.

[0030] S24, based on the eigenvector corresponding to the extracted maximum eigenvalue Calculate the pore extension direction.

[0031] The eigenvector corresponding to the maximum eigenvalue of each pore voxel point is used as the extension direction of the pore at that point Azimuth: inclination: in, is the eigenvector corresponding to the largest eigenvalue The three-dimensional components of are used to calculate the direction angle, is the azimuth, is the inclination angle.

[0032] Refer to the attached diagram for the distribution of the angle between the pore extension direction and the z-axis output according to the calculation results. Figure 2 As shown in the figure, the output pore extension azimuth distribution rose diagram is shown in the attached figure. Figure 3 shown.

[0033] Step S3 includes the following sub-steps: S31. Use 26-neighborhood connected domain analysis to mark the pore areas.

[0034] S32. Calculate the sphericity of each connected pore.

[0035] in Indicates the degree to which the pores are close to spherical, with a value range of 0 to 1, 1 being a perfect sphere, V being the pore volume, and S being the pore surface area.

[0036] According to the calculation results, the sphericity of the connected pores is output as shown in the attached figure. Figure 4 shown.

[0037] S33. Calculate the aspect ratio of the connected pores based on principal axis analysis.

[0038] in represents the aspect ratio of connected pores, and are the longest and shortest axis lengths of the connected pores, respectively.

[0039] Refer to the attached figure for the aspect ratio of connected pores output according to the calculation results. Figure 5 shown.

[0040] S34. Calculate the surface-to-volume ratio of connected pores.

[0041] in Expresses the ratio of surface area to volume.

[0042] According to the calculation results, the surface ratio of connected pores is output as shown in the attached figure. Figure 6 shown.

[0043] Step S4 includes the following sub-steps: S41. Determine an adaptive box size sequence.

[0044] in A dynamic sequence representing adaptive box sizes, (Image size); the image size refers to the size of the CT scan three-dimensional image of the reef limestone.

[0045] S42, for each adaptive box size S k The CT scan 3D image of reef limestone is divided into The box network uses zero padding to handle non-divisible boundaries and counts the number of boxes containing at least one pore voxel N(S K ).

[0046] S43. Establish a double logarithmic coordinate relationship: in represents the minimum number of boxes required to cover the holes, Indicates the reciprocal of the box size.

[0047] S44. Calculate the fractal dimension by linear regression: in, represents the number of boxes required for the pore coverage, Represents the logarithm of the reciprocal of the box size.

[0048] Refer to the attached diagram for the fractal dimension fitting diagram and the double logarithmic coordinate diagram of the box counting method based on the calculation results. Figure 7 shown.

[0049] Step S5 includes the following sub-steps: S51. Calculate the 3D Euclidean distance field of pores: S52, extracting the distance value of the pore area and converting it into physical size; S53, generating a pore size distribution density histogram and a cumulative curve based on the distance values ​​of the extracted pore areas; The pore size distribution density histogram output based on the extracted pore area distance value is shown in the attached figure. Figure 8 As shown in the figure, the cumulative pore size distribution curve based on the extracted pore area distance value is shown in the attached figure. Figure 9 shown.

[0050] S54, using the bwconncomp function to perform 26-domain connectivity analysis; S55. Calculate the maximum proportion of connected pores.

[0051] Step S6 includes the following sub-steps: S61. Establish material property interpolation model: in E e is the equivalent elastic modulus, E matrix is the matrix elastic modulus, E pore is the poroelastic modulus, is a binary matrix value, and p=3 is a penalty factor; S62, assembling the global stiffness matrix based on the element stiffness matrix KE; in K is the global stiffness matrix, e is the unit after the finite element mesh is discretized, is the unit equivalent elastic modulus, To precompute the element stiffness matrix, the assembly process uses sparse matrix technology.

[0052] S63. Apply fixed boundary conditions: Constrain the displacement of specific degrees of freedom of the model; Fix the first three degrees of freedom (ux, uy, uz) Release the remaining degrees of freedom S64. Solve the displacement field: in U is the node displacement vector, K is the global stiffness matrix, F is the load vector.

[0053] S65. Extract equivalent elastic parameters: where Q is the fourth-order equivalent elastic tensor, σ and ε represent macroscopic stress and macroscopic strain, respectively.

[0054] Generate 3D elastic parameter visualization: Display the spatial distribution of Q11, Q22, and Q33; Use heat maps to display anisotropic characteristics. Figure 10 shown.

[0055] The above examples demonstrate the beneficial effects of the present invention. The disclosed method for obtaining equivalent mechanical parameters for high-porosity reef limestone utilizes structural tensor feature analysis and a three-dimensional directional quantization model (azimuth angle φ, inclination angle θ), combined with rose diagram visualization technology, to achieve a fully parametric characterization of the principal extension direction of pore groups. This addresses the failure of traditional homogenization methods to predict deformation in rock masses with oriented pores. By establishing a collaborative analysis system of three-dimensional indices (sphericity, aspect ratio, and surface-to-volume ratio), kernel density distribution statistics reveal pore morphology spectrum characteristics, addressing the problem of traditional methods misjudging stress concentration around complex pores. An adaptive power-scaling box counting algorithm (i.e., dynamically adjusting the box size by powers of 2) accurately calculates the fractal dimension of the pore surface, providing a key parameter foundation for establishing a quantitative correlation between fractal dimension and macroscopic modulus. By establishing a "pore characteristic matrix → dynamic penalty factor → equivalent elastic tensor" mapping model, the transfer of microscopic pore parameters to macroscopic mechanical response is achieved.

[0056] Obviously, those skilled in the art may make various changes and modifications to the present invention without departing from the spirit and scope of the present invention. Thus, if such changes and modifications fall within the scope of the claims and their equivalents, the present invention is intended to include such changes and modifications.

Claims

1. A method for obtaining equivalent parameters of high-porosity reef limestone, comprising the following steps: S1, constructing a three-dimensional matrix-pore binary matrix based on CT scan data; S2. Based on the constructed three-dimensional matrix-pore binary matrix, three-dimensional pore direction vectorization is performed to describe the geometric characteristics of the pores in three-dimensional space; S3, extracting pore morphology parameters based on the constructed three-dimensional matrix-pore binary matrix; S4. Calculate the pore fractal dimension based on the constructed three-dimensional matrix-pore binary matrix; S5. Analyze pore size distribution and connectivity based on the constructed three-dimensional matrix-pore binary matrix; S6. Calculate the equivalent elastic parameters based on the constructed three-dimensional matrix-pore binary matrix.

2. The method for obtaining equivalent parameters of high-porosity reef limestone according to claim 1, characterized in that: Step S2 includes the following sub-steps: S21, calculating the three-dimensional spatial gradient field in the pore area; S22. constructing a structure tensor based on the obtained three-dimensional spatial gradient field of the pore region; S23, performing eigenvalue decomposition on the constructed structure tensor, and extracting the eigenvector corresponding to the maximum eigenvalue; S24. Calculate the pore extension direction based on the eigenvector corresponding to the extracted maximum eigenvalue.

3. The method for obtaining equivalent parameters of high-porosity reef limestone according to claim 2, characterized in that: In step S21, the three-dimensional spatial gradient field of the pore area is calculated according to the following formula: in, is the gradient operator, They represent the partial derivatives along the x, y, and z directions, respectively, and are used to calculate the gradient field in the pore region and capture the changes in the pore boundary.

4. The method for obtaining equivalent parameters of high-porosity reef limestone according to claim 3, characterized in that: In step S22, a structure tensor is constructed based on the obtained three-dimensional spatial gradient field of the pore region according to the following formula: in Represents a Gaussian filtering operation.

5. The method for obtaining equivalent parameters of high-porosity reef limestone according to claim 1, characterized in that: Step S3 includes the following sub-steps: S31, using 26-neighborhood connected domain analysis to mark pore areas; S32, calculating the sphericity of each connected pore; S33, calculating the aspect ratio of connected pores based on principal axis analysis; S34. Calculate the surface-to-volume ratio of connected pores.

6. The method for obtaining equivalent parameters of high-porosity reef limestone according to claim 1, characterized in that: Step S4 The following sub-steps are included: S41, determining an adaptive box size sequence; S42, for each adaptive box size S k Divide the 3D image into The box network uses zero padding to handle non-divisible boundaries and counts the number of boxes containing at least one pore voxel N(S K ); S43, establishing a double logarithmic coordinate relationship; S44. Calculate the fractal dimension by linear regression.

7. The method for obtaining equivalent parameters of high-porosity reef limestone according to claim 1, characterized in that: Step S5 includes the following sub-steps: S51, calculating the pore three-dimensional Euclidean distance field; S52, extracting the distance value of the pore area and converting it into physical size; S53, generating a pore size distribution density histogram and a cumulative curve based on the distance values ​​of the extracted pore areas; S54, using the bwconncomp function to perform 26-domain connectivity analysis; S55. Calculate the maximum proportion of connected pores.

8. The method for obtaining equivalent parameters of high-porosity reef limestone according to claim 1, characterized in that: Step S6 includes the following sub-steps: S61. Establish material property interpolation model; S62, assembling the global stiffness matrix based on the element stiffness matrix; S63, applying fixed boundary conditions; S64, solve the displacement field; S65. Extract equivalent elastic parameters.

9. The method for obtaining equivalent parameters of high-porosity reef limestone according to claim 8, characterized in that: In step S61, a material property interpolation model is established using the following formula: in E e is the equivalent elastic modulus, E matrix is the matrix elastic modulus, E pore is the poroelastic modulus, is a binary matrix value, and p=3 is a penalty factor.

10. The method for obtaining equivalent parameters of high-porosity reef limestone according to claim 8, characterized in that: In step S63, the fixed boundary conditions imposed include the displacement of the specific degree of freedom of the constrained model, Fix the first three degrees of freedom (ux, uy, uz) and release the remaining degrees of freedom.

Images