Method for obtaining equivalent parameters of high porosity reef limestone

Abstract

The application discloses a method for obtaining equivalent parameters of high porosity reef limestone, and relates to the field of geotechnical engineering, S1, a three-dimensional matrix of matrix-pore is constructed based on CT scanning data; S2, three-dimensional pore direction vectorization is carried out based on the constructed three-dimensional matrix of matrix-pore, so that the geometric characteristics of the pore are described in three-dimensional space; S3, pore morphological parameters are extracted based on the constructed three-dimensional matrix of matrix-pore; S4, pore fractal dimension is calculated based on the constructed three-dimensional matrix of matrix-pore; S5, pore size distribution and connectivity are analyzed based on the constructed three-dimensional matrix of matrix-pore; and S6, equivalent elastic parameters are calculated based on the constructed three-dimensional matrix of matrix-pore. Through CT scanning, pore characteristic statistical analysis and finite element homogenization method, equivalent mechanical parameters of high porosity reef limestone are obtained, which lays a foundation for marine engineering foundation design, geological disaster assessment and the like.

Description

TECHNICAL FIELD

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

[0002] Reef limestone is the remains of reef-building coral communities, and is a special rock and soil formed through geological processes. Its unique diagenetic mechanism makes its pore structure different from that of ordinary rocks, and its structure is relatively complex. A large number of studies have shown that the porosity of reef limestone fluctuates greatly. The porosity of reef limestone can reach 20%-40% as mentioned in the book Sedimentary Petrology. The pore connectivity of shallow reef limestone (SRL) is good and as high as 55.3±3.2%, while the porosity of deep reef limestone (DRL) is only 4.9±1.6% due to its dense structure. Due to the difference and large proportion of the pores of reef limestone, the pores of reef limestone have a huge impact on the structure and mechanical properties of the rock.

[0003] The traditional method relies on laboratory mechanical tests to obtain the pore parameters of macroscopic reef limestone, and the test results obtained have large dispersion, cannot reflect the influence of micro-pores on the structure, have obvious limitations, and cannot meet the requirements of engineering design. In addition, the scale effect is often ignored in the prior art, and the sample size in the laboratory is usually centimeter level, which cannot reflect the multi-scale pore coupling effect of meter-level rock mass in actual engineering.

[0004] In terms of practical application, although CT scanning technology combined with commercial software (commonly using global threshold segmentation method such as adaptive threshold segmentation) can obtain the pore structure, it generally lacks quantitative statistics of the anisotropy characteristics of the pores. Although the pore structure can be obtained (such as the adaptive threshold segmentation method), there is a lack of quantitative statistics of the anisotropy characteristics of the pores, resulting in insufficient modeling accuracy, which is specifically manifested in the lack of representation of anisotropy, and the inability to distinguish the difference in pore shape in the horizontal / vertical direction. Furthermore, the common morphological parameters in existing research simplify the pores as equivalent spheres, ignoring the complex topological characteristics of the real pores, and the error of the stress concentration coefficient caused by the aspect ratio of the pores is large.

[0005] In summary, the existing technology has systematic deficiencies in obtaining equivalent parameters of reef limestone and accurately correlating the micro-pore structure of reef limestone with the macro-mechanical behavior, which will directly lead to significant errors in the prediction of mechanical parameters relied on in engineering design, threatening the safety of the structure. SUMMARY

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

[0007] To achieve the above object, the technical scheme adopted by the present application is as follows: A method for obtaining equivalent parameters of high porosity reef limestone, comprising the following steps:

[0008] S1, constructing a three-dimensional matrix of matrix-pore based on CT scanning data;

[0009] S2, based on the constructed three-dimensional matrix of matrix-pore, performing three-dimensional pore direction vectorization to describe the geometric characteristics of the pores in three-dimensional space;

[0010] S3, based on the constructed three-dimensional matrix of matrix-pore, extracting pore morphological parameters;

[0011] S4, based on the constructed three-dimensional matrix of matrix-pore, calculating pore fractal dimension;

[0012] S5, based on the constructed three-dimensional matrix of matrix-pore, analyzing pore size distribution and connectivity;

[0013] S6, based on the constructed three-dimensional matrix of matrix-pore, calculating equivalent elastic parameters.

[0014] Further, step S2 includes the following sub-steps:

[0015] S21, calculating the three-dimensional spatial gradient field of the pore region;

[0016] S22, constructing a structure tensor based on the obtained three-dimensional spatial gradient field of the pore region;

[0017] S23, performing eigenvalue decomposition on the constructed structure tensor to extract the eigenvector corresponding to the maximum eigenvalue;

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

[0019] Further, the three-dimensional spatial gradient field of the pore region is calculated according to the following formula in step S21,

[0020]

[0021] wherein, is a gradient operator, respectively represent partial derivatives along x, y, and z directions, which are used to calculate the gradient field of the pore region and capture the changes in pore boundaries.

[0022] Further, the structure tensor is constructed based on the obtained three-dimensional spatial gradient field of the pore region according to the following formula in step S22,

[0023]

[0024] wherein represents a Gaussian filter operation.

[0025] Further, the step S3 comprises the following sub-steps:

[0026] S31, using 26-neighborhood connected domain analysis to mark the pore region;

[0027] S32, calculating the sphericity of each connected pore;

[0028] S33, calculating the aspect ratio of the connected pore based on the principal axis analysis;

[0029] S34, calculating the surface-volume ratio of the connected pore.

[0030] Further, the step S4 comprises the following sub-steps:

[0031] S41, determining the adaptive box size sequence;

[0032] S42, for each adaptive box size S k dividing the three-dimensional image into a box network of , using zero padding to handle non-integer division boundaries, and counting the number of boxes N(S K ) containing at least one pore voxel;

[0033] S43, establishing a double logarithmic coordinate relationship;

[0034] S44, calculating the fractal dimension by linear regression.

[0035] Further, the step S5 comprises the following sub-steps:

[0036] S51, calculating the three-dimensional Euclidean distance field of the pore;

[0037] S52, extracting the distance value of the pore region and converting the physical size;

[0038] S53, generating the pore size distribution density histogram and cumulative curve based on the extracted distance value of the pore region;

[0039] S54, performing 26-field connected domain analysis using the bwconncomp function;

[0040] S55, finding the maximum connected pore occupancy ratio.

[0041] Further, the step S6 comprises the following sub-steps:

[0042] S61, establishing a material property interpolation model;

[0043] S62, assembling a global stiffness matrix based on the element stiffness matrix;

[0044] S63, applying fixed boundary conditions;

[0045] S64, solving displacement field;

[0046] S65, extracting equivalent elastic parameters.

[0047] Further, in the step S61, the material property interpolation model is established by the following formula:

[0048]

[0049] wherein is a binary matrix value, and p=3 is a penalty factor.

[0050] Further, in the step S63, the fixed boundary conditions applied include restricting the displacement of the specific degrees of freedom of the model,

[0051] the first three degrees of freedom (ux, uy, uz) are fixed, and the remaining degrees of freedom are released.

[0052] The technical scheme adopted by the present application has the beneficial effects that the method for obtaining equivalent parameters of high-porosity reef limestone can solve the technical defects that most of the existing rock mass pore characterization techniques ignore the spatial orientation characteristics of pore groups, simplify pores as isotropic spherical structures, and seriously underestimate the anisotropy of material mechanics. Through structural tensor feature analysis and three-dimensional directional quantization model (azimuth angle φ, inclination angle θ), combined with rose diagram visualization technology, the full parameterization characterization of the main extension direction of the pore group is realized, and the deformation prediction failure problem of the traditional homogenization method in the rock mass containing directional pores is solved.

[0053] The conventional analysis of the existing commercial image processing software (such as Avizo) lacks systematic quantitative statistics of specific pore shape parameters such as sphericity, aspect ratio, and surface-to-volume ratio. This leads to the technical defects that non-spherical pores, especially high-aspect-ratio tubular pores, are equivalent to spherical pores, and the stress concentration effect of the pore periphery is seriously underestimated or incorrectly calculated. By establishing a three-dimensional index coordination analysis system of sphericity (Sphericity)-aspect ratio (Aspect Ratio)-surface-to-volume ratio (Surface-to-Volume Ratio), the pore morphology spectrum characteristics are revealed through kernel density distribution statistics, and the misjudgment problem of the traditional method for complex pore periphery stress concentration is solved.

[0054] The traditional fractal dimension calculation method such as fixed step box counting method has a sharp decrease in accuracy under the condition of high porosity (>30%), which leads to the technical defect that the correlation between micro-roughness and macro-mechanical response is not exhibited. Through the adaptive power scaling box counting algorithm (i.e. the box size is dynamically adjusted according to the power of 2), the fractal dimension of the pore surface is accurately calculated, which provides a key parameter basis for establishing the quantitative correlation between fractal dimension and macro modulus.

[0055] The traditional homogenization algorithm regards the pores as blank areas without geometric characteristics, ignores the complex influence of pore directionality, morphological topology and surface fractal characteristics on the elastic modulus, and the damage effect caused by micro-pores is seriously underestimated. BRIEF DESCRIPTION OF DRAWINGS

[0056] Figure 1 is a flow chart of a high-porosity reef limestone equivalent parameter acquisition method according to an embodiment of the present application;

[0057] Figure 2 is a pore extension direction and z-axis angle distribution histogram generated by the high-porosity reef limestone equivalent parameter acquisition method according to an embodiment of the present application;

[0058] Figure 3 is a pore extension azimuth angle distribution horizontal projection rose diagram generated by the high-porosity reef limestone equivalent parameter acquisition method according to an embodiment of the present application;

[0059] Figure 4 is a connected pore sphericity distribution diagram generated by the high-porosity reef limestone equivalent parameter acquisition method according to an embodiment of the present application;

[0060] Figure 5 is a connected pore aspect ratio distribution diagram generated by the high-porosity reef limestone equivalent parameter acquisition method according to an embodiment of the present application;

[0061] Figure 6 is a connected pore surface-volume ratio distribution diagram generated by the high-porosity reef limestone equivalent parameter acquisition method according to an embodiment of the present application;

[0062] Figure 7 is a fractal dimension fitting diagram and a double logarithmic coordinate diagram of the box counting method generated by the high-porosity reef limestone equivalent parameter acquisition method according to an embodiment of the present application;

[0063] Figure 8 is a pore diameter distribution density histogram generated by the high-porosity reef limestone equivalent parameter acquisition method according to an embodiment of the present application;

[0064] Figure 9 is a pore diameter distribution cumulative curve generated by the high-porosity reef limestone equivalent parameter acquisition method according to an embodiment of the present application;

[0065] Figure 10 is an equivalent elastic modulus thermal diagram generated by the high-porosity reef limestone equivalent parameter acquisition method according to an embodiment of the present application. DETAILED DESCRIPTION

[0066] The application will be described in detail below with reference to the accompanying drawings and embodiments.

[0067] Embodiment one

[0068] Referring to the accompanying drawings Figure 1 , the application provides a method for obtaining equivalent parameters of high porosity reef limestone, comprising the following steps:

[0069] S1, constructing a three-dimensional matrix of matrix-pore based on CT scanning data;

[0070] S2, based on the constructed three-dimensional matrix of matrix-pore, performing three-dimensional pore direction vectorization to describe the geometric characteristics of the pores in three-dimensional space;

[0071] S3, based on the constructed three-dimensional matrix of matrix-pore, extracting pore morphological parameters;

[0072] S4, based on the constructed three-dimensional matrix of matrix-pore, calculating pore fractal dimension;

[0073] S5, based on the constructed three-dimensional matrix of matrix-pore, analyzing pore size distribution and connectivity;

[0074] S6, based on the constructed three-dimensional matrix of matrix-pore, calculating equivalent elastic parameters.

[0075] Step S2 includes the following sub-steps:

[0076] S21, calculating the three-dimensional space gradient field of the pore region according to the following formula.

[0077]

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

[0079] S22, based on the obtained three-dimensional space gradient field of the pore region, constructing a structure tensor according to the following formula.

[0080]

[0081] wherein represents a Gaussian filtering operation.

[0082] S23, performing eigenvalue decomposition on the constructed structure tensor , extracting the eigenvector corresponding to the maximum eigenvalue .

[0083] where is the eigenvector matrix, which stores the pore orientation information, is the eigenvalue diagonal matrix, which represents the importance weight of the orientation, represents the eigenvector corresponding to the largest eigenvalue.

[0084] S24, based on the eigenvector corresponding to the largest eigenvalue extracted calculate the pore extension direction.

[0085] The eigenvector corresponding to the largest eigenvalue of each pore voxel point is taken as the extension direction of the pore at that point

[0086] Azimuth angle:

[0087] Dip angle:

[0088] where, is the eigenvector corresponding to the largest eigenvalue three-dimensional component, used to calculate the directional angle, is the azimuth angle, is the dip angle.

[0089] The pore extension direction output according to the calculation result and the angle distribution diagram with the z-axis are shown in Fig. 6, and the output pore extension azimuth angle distribution rose diagram is shown in Fig. 7. Figure 2 The pore extension direction output according to the calculation result and the angle distribution diagram with the z-axis are shown in Fig. 6, and the output pore extension azimuth angle distribution rose diagram is shown in Fig. 7. Figure 3 The pore extension direction output according to the calculation result and the angle distribution diagram with the z-axis are shown in Fig. 6, and the output pore extension azimuth angle distribution rose diagram is shown in Fig. 7.

[0090] Step S3 includes the following sub-steps:

[0091] S31, using 26-neighborhood connected component analysis to mark the pore region.

[0092] S32, calculate the sphericity of each connected pore.

[0093]

[0094] where represents the degree of pore approaching to a sphere, the value range is 0 to 1, 1 is perfect sphere, V is the pore volume, and S is the pore surface area.

[0095] The connected pore sphericity output according to the calculation result is shown in Fig. 8. Figure 4 The connected pore sphericity output according to the calculation result is shown in Fig. 8.

[0096] S33, based on the principal axis analysis, calculate the aspect ratio of the connected pore.

[0097]

[0098] where represents the aspect ratio of the connected pore, and respectively the longest and shortest axis length of the connected pores.

[0099] The aspect ratio of the connected pores output according to the calculation result is shown in FIG. 4. Figure 5

[0100] S34, calculate the surface-volume ratio of the connected pores.

[0101]

[0102] wherein represents the ratio of the surface area to the volume.

[0103] The surface-volume ratio of the connected pores output according to the calculation result is shown in FIG. 5. Figure 6

[0104] Step S4 includes the following sub-steps:

[0105] S41, determine the adaptive box size sequence.

[0106]

[0107] wherein represents the dynamic sequence of the adaptive box size, (image size); wherein the image size is the size of the CT scan three-dimensional image of the reef limestone.

[0108] S42, for each adaptive box size S k divide the CT scan three-dimensional image of the reef limestone into a box network of , use zero padding to handle non-integer boundary, and count the number of boxes N(S K ) containing at least one pore voxel.

[0109] S43, establish a double logarithmic coordinate relationship:

[0110]

[0111] wherein represents the minimum number of boxes required to cover the pores, represents the reciprocal of the box size.

[0112] S44, calculate the fractal dimension by linear regression:

[0113]

[0114] wherein, represents the logarithm of the number of boxes required for pore coverage, represents the logarithm of the reciprocal of the box size.

[0115] ​​The fractal dimension fitting figure output according to the calculation result and the double logarithmic coordinate figure of the box counting method are shown in Figs. 9 and 10. Figure 7

[0116] Step S5 includes the following sub-steps:

[0117] S51, calculate the three-dimensional Euclidean distance field of the pores:

[0118]

[0119] S52, extract the distance value of the pore region and convert the physical size;

[0120] S53, generate the pore size distribution density histogram and the cumulative curve based on the extracted distance value of the pore region;

[0121] The pore size distribution density histogram output based on the extracted distance value of the pore region is shown in Fig. 6, and the pore size distribution cumulative curve output based on the extracted distance value of the pore region is shown in Fig. 7. Figure 8 Figure 9

[0122] S54, perform 26-field connected component analysis by using the bwconncomp function;

[0123] S55, calculate the maximum connected pore proportion.

[0124] Step S6 includes the following sub-steps:

[0125] S61, establish a material property interpolation model:

[0126]

[0127] wherein E e is the equivalent elastic modulus, E matrix is the matrix elastic modulus, E pore is the pore elastic modulus, is a binary matrix value, and p=3 is a penalty factor;

[0128] S62, assemble the global stiffness matrix based on the element stiffness matrix KE;

[0129]

[0130] wherein K is the global stiffness matrix, e is the element after the finite element mesh is discretized, is the equivalent elastic modulus of the element,

[0131] ​​​For pre-computed unit stiffness matrix, the assembly process adopts sparse matrix technique.

[0132] S63, apply fixed boundary conditions:

[0133] Constrain displacement of model-specific DOFs;

[0134] Fix the first three DOFs (ux, uy, uz)

[0135] Release the remaining DOFs

[0136] S64, solve the displacement field:

[0137]

[0138] where U is the nodal displacement vector, K is the global stiffness matrix, F is the load vector.

[0139] S65, extract equivalent elastic parameters:

[0140]

[0141] where Q is the fourth-order equivalent elastic tensor, and σ and ε represent macroscopic stress and macroscopic strain, respectively.

[0142] Generate three-dimensional elastic parameter visualization: display the spatial distribution of Q11, Q22, Q33; use a heat map to show the anisotropy characteristics Referring to the drawings Figure 10 shown.

[0143] As can be seen from the above embodiments, the beneficial effects of the present application are that the disclosed method for obtaining equivalent parameters of high-porosity reef limestone realizes full-parameterization characterization of the main extension direction of the pore group by structural tensor feature analysis and a three-dimensional directional quantization model (azimuth angle φ, inclination angle θ), combined with rose diagram visualization technology, solves the deformation prediction failure problem of traditional homogenization methods in rock bodies containing directional pores. By establishing a three-dimensional index coordination analysis system of sphericity (Sphericity)-aspect ratio (AspectRatio)-surface-to-volume ratio (Surface-to-Volume Ratio), the pore morphology spectrum characteristics are revealed by kernel density distribution statistics, solving the misjudgment problem of traditional methods for stress concentration around complex pores. By using the adaptive power scaling box counting algorithm (i.e., the box size is dynamically adjusted according to the power of 2), the fractal dimension of the pore surface is accurately calculated, providing a key parameter basis for establishing the quantitative correlation between fractal dimension and macroscopic modulus. By establishing a "pore feature matrix → dynamic penalty factor → equivalent elastic tensor" mapping model, the transmission from micro-pore parameters to macro-mechanical response is realized.

[0144] It will be apparent to those skilled in the art that various modifications and variations can be made to the present application without departing from the spirit or scope of the application. Thus, it is intended that the present application cover modifications and variations of this application provided they come within the scope of the appended claims and their equivalents.

Claims

1. A method for obtaining equivalent chemical parameters of high-porosity reef limestone, comprising the following steps: S1. Construct a three-dimensional matrix-pore binary matrix based on CT scan data; S2. Based on the constructed three-dimensional matrix-pore binary matrix, the three-dimensional pore direction is vectorized to describe the geometric features of the pores in three-dimensional space. S3. Based on the constructed three-dimensional matrix-pore binary matrix, extract pore morphology parameters; S4. Calculate the fractal dimension of the pores based on the constructed three-dimensional matrix-pore binary matrix; S5. Based on the constructed three-dimensional matrix-pore binary matrix, analyze the pore size distribution and connectivity; S6. Based on the constructed three-dimensional matrix-pore binary matrix, calculate the equivalent elastic parameters; Step S2 includes the following sub-steps: S21. Calculate the three-dimensional spatial gradient field of the pore region; S22. Construct the structural tensor based on the obtained three-dimensional spatial gradient field of the pore region; S23. Perform eigenvalue decomposition on the constructed structural tensor and extract the eigenvector corresponding to the largest eigenvalue. S24. Calculate the pore extension direction based on the feature vector corresponding to the extracted maximum feature value; In step S21, the three-dimensional spatial gradient field of the pore region is calculated according to the following formula. in, For gradient operators, These represent the partial derivatives along the x, y, and z directions, respectively, used to calculate the gradient field in the pore region and capture changes in the pore boundary; In step S22, the structural tensor is constructed based on the obtained three-dimensional spatial gradient field of the pore region according to the following formula. Where <·> represents Gaussian filtering operation.

2. The method for obtaining equivalent chemical parameters of high-porosity reef limestone according to claim 1, characterized in that, Step S3 includes the following sub-steps: S31. Use 26-neighbor connected component analysis to mark the pore region; S32. Calculate the sphericity of each connected pore; S33. Calculate the aspect ratio of the connected pores based on the main axis analysis; S34. Calculate the surface-to-volume ratio of connected pores.

3. The method for obtaining equivalent chemical parameters of high-porosity reef limestone according to claim 1, characterized in that, Step S4 Includes the following sub-steps: S41. Determine the adaptive box size sequence; S42, For each adaptive box size S k Divide the 3D image into S K ×S K ×S K The box network is constructed using zero-filling to handle non-divisible boundaries, and the number of boxes N(S) containing at least one pore voxel is counted. K ); S43. Establish a double logarithmic coordinate relationship; S44. Calculate the fractal dimension using linear regression.

4. The method for obtaining equivalent chemical parameters of high-porosity reef limestone according to claim 1, characterized in that, Step S5 includes the following sub-steps: S51. Calculate the three-dimensional Euclidean distance field of the pores; S52. Extract the distance values ​​of the pore region and convert them into physical dimensions; S53. Based on the distance values ​​of the extracted pore regions, generate a histogram of pore size distribution density and a cumulative curve; S54. Use the bwconncomp function to perform neighborhood connectivity analysis (26). S55. Calculate the percentage of the largest connected pores.

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

6. The method for obtaining equivalent chemical parameters of high-porosity reef limestone according to claim 5, characterized in that: In step S61, a material property interpolation model is established using the following formula: E e =E matrix r p +E pore (1-p p ) Where E e E is the equivalent elastic modulus. matrix E represents the matrix elastic modulus. pore Let ρ be the pore elastic modulus, ρ be the binary matrix value, and p = 3 be the penalty factor.

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

Images