Research on establishing numerical model of geo material based on CT image analysis

Geotechnical engineering material is a kind of heterogeneous composite material. The stone, mineral, gravel, hole, and crack contained in it have different physical and mechanical properties, and their responses to external forces are very different. For example, stress distribution, crack propagation, and failure mode are closely related to material heterogeneity and microstructure. This paper presents a finite element model method for processing CT images of geotechnical materials by using digital image technology. This paper presents a finite element model method for processing CT images of geotechnical materials by using digital image technology. The theory of digital image processing is applied to geotechnical CT image processing to realize pseudo-color enhancement of CT image, and the histograms of different geotechnical CT numbers are obtained. Canny operator is used to detect the edge of geotechnical CT image, and the binary image of geotechnical microstructure is obtained. According to the color change of pseudo-color enhanced image and CT histogram, the meso-structure and variation of rock and soil can be quantitatively analyzed. The finite element model established by this method can fully consider the heterogeneity of geotechnical materials, especially the effect of void distribution on the mechanical properties of geotechnical materials.

With the in-depth study of geotechnical mechanical properties, scholars pay more and more attention to the study of its meso-mechanical properties, resulting in a variety of meso-structural finite element models, such as lattice model, beam-particle model, random aggregate model, and random mechanical properties model [1]. Most of the above models are based on random sampling method and statistical knowledge. Although the aggregate distribution and shape of the model has been more and more close to the real structure of rock and soil, but inevitably there are differences with the real meso-structure of rock and soil. Studying the rupture process of geotechnical materials has important theoretical significance and practical value for evaluating the safety state of geotechnical engineering, understanding the stability of geotechnical engineering structure, taking reasonable supporting measures, and improving the design level of geotechnical and underground structural engineering. The nonlinearity of macroscopic mechanical behavior of geotechnical materials in the process of fracture is due to the heterogeneity of their meso-structure [2,3,4]. Under the same experimental and loading conditions, the true microstructure and the mechanical properties of various meso-media play a decisive role in the stress and deformation responses of the specimens. Therefore, further research on the effective methods to measure the spatial distribution of the meso-media in geotechnical materials determines the constitutive relations and descriptive parameters of each medium and introduces them into the numerical model, so as to study and observe the mechanical properties and failure laws of geotechnical materials more comprehensively and more in accordance with the natural laws [5, 6].

With the development of digital image processing technology, digital image processing technology is introduced to study the internal structure of rock and soil. Digital image processing technology provides convenience for computer reconstruction, processing, and simulation of image data. The application of digital image processing technology to study the damage characteristics of rock and soil, especially the damage characteristics of rock and soil, can be the development of related research fields. The effective testing technology [7,8,9,10] is provided, and the freeze-thaw rock mass can be described and described more accurately by means of numerical simulation, which reflects the actual stress of rock and soil. The finite element method based on digital image processing provides a new way to study the mechanical properties of geotechnical materials under freeze-thaw environment. Digital image processing (DIP) coupled with finite element method (FEM) is another effective way to correlate the fine microstructure of geotechnical materials with the macroscopic mechanical response [11, 12]. Yuezhongqi proposed the digital image finite element analysis method (DIP-based FEM) for geotechnical engineering materials. This method integrates digital image processing theory, geometric vector conversion technology, and finite element mesh automatic generation principle [13,14,15,16]. Firstly, the real meso-structure of material is obtained by digital image algorithm, and then vectorized. The finite element mesh of material meso-structure is generated by automatic mesh generation technology. The effect of meso-structure on stress distribution is studied by finite element method. But in this method, the digital image processing technology is only aimed at grayscale image [17, 18].

Later, Chen Sha improved Yue’s method [19,20,21], used the variation of values in the HSI (hue, saturation, intensity) color space to obtain the meso-structure of the material, and combined with FDM to analyze the mechanics of heterogeneous materials; the failure modes of granite under uniaxial compression and uniaxial direct tensile loading are simulated. Then, some scholars apply digital image processing technology to seepage analysis, Sheng et al. [22]. By processing the digital image of rock mass, the relationship between the value and the distribution of permeability coefficients of different structures of rock mass is established based on the different brightness values of the image corresponding to different structures of rock mass, and the distribution of physical parameters is input into the developed numerical model of multi-physical hazard of rock mass based on FEMLAB. The unsteady seepage process in fractured rock mass is analyzed. Compared with the statistical sampling method, this method has obvious advantages, and the results can better reflect the real physical phenomena. J. Su and Z. Liu et al. [23, 24] used the method of digital image processing to characterize the heterogeneity of rock and soil by using the value in the color space of the moon sign. The whole process of the heterogeneous rock and soil rupture under the combined action of external load and seepage disaster was studied. The strength and overall permeability of the specimen showed obvious anisotropy in macroscopic. Niknejad et al. [25] used image processing methods such as image segmentation and edge detection to process small-scale crystalline rock images, obtained the transgranular and intergranular micro-cracks, characterized the network of micro-cracks in rock and soil, and solved the permeability coefficient of samples; the effects of transgranular crack network and intergranular fracture network on the permeability of samples were studied. Y. Liu et al. [26] established the discrete element calculation model by using the CT image of the subgrade material. Sanchez-Romero et al. [27] established the concept model of the soil-rock mixture microstructure by using the digital image processing technology and established the numerical model by the vector method. L. Tan et al. [28] studied the three-dimensional microstructure of the rock and soil based on the CT image, combined with ANSYS, and a simplified numerical calculation model is established. L.I. Meng-Yi and  J. Sheng et al [29,30,31] applied CT images to establish the finite element plane model of asphalt mixture.

According to past experience, this paper takes self-made porous rock and soil materials and aggregate rock and soil materials, for example, obtains its internal structure image by CT technology, studies the method of image feature information extraction and bitmap vectorization, and establishes its numerical calculation model. The theory of digital image processing is applied to geotechnical CT image processing to realize pseudo-color enhancement of CT image, and the histograms of different geotechnical CT numbers are obtained. Canny operator is used to detect the edge of geotechnical CT image, and the binary image of geotechnical microstructure is obtained. Pseudo-color enhancement of geotechnical CT image can improve the resolution of CT image and reduce the error of visual judgment. According to the color change of pseudo-color enhancement image and CT number histogram, the meso-structure and variation of geotechnical can be quantitatively analyzed.

2.1 Digital image processing

The methods of digital image processing are divided into two categories: the airspace method and the transform domain method. Space-domain method is to treat the image as a set of pixels in the plane and directly operate the pixels in the space where the image is located, that is, directly process the two-dimensional function. There are two kinds of spatial method. One is point processing based on pixels. Each processing of the image is carried out on each pixel. The result of processing is independent of other pixels. The other one is neighborhood processing method based on the template, in which each processing of pixels is a small sub-image of the template, that is, the adjacent set of pixels to operate.

The transform domain method is to transform the image orthogonally. The image is processed indirectly in the transform domain of the image and then inversely transformed into the spatial domain. Because the transform domain algorithm has higher requirement for system memory, this paper only uses the airspace method. In the following processes of image enhancement and image segmentation, spatial filtering is used. Spatial filtering algorithm is essentially a convolution operator for image processing and is a small area of local image processing operations. Convolution operation is a transform domain method, which modifies the Fourier transform of the image to realize the image processing. Setting the original image F, the processed image is H, G is a convolution operator, and m × n is an array; the array G and digital images F of the convolution processing is defined as:

H=F*G(* represents convolution)

In order to make the calculation intuitive and simple, the convolution operation is generally converted from the operation in the transform domain to the operation in the spatial domain. The concrete method is to transform the convolution operation into the form of multiplication of two matrices. The edge of an operator is odd, and in G, the elements are symmetrically located in the center of the operator; in practice, the filter matrix with smaller order convolution operator G is the template used in detection, also known as convolution kernel.

The characteristics of digital image processing are as follows:

1. 1.

   Good reproducibility

As long as the input image and processing method remain unchanged, the digital image processing results will not change, can be well reproduced, and have good repeatability. There is no randomness in visual processing.

1. 2.

   High accuracy

Current technology can transform a digital image into a two-dimensional array of arbitrary size. For a computer, no matter the size of the array or the number of bits of pixels, its processing program is the same. That is to say, as long as the array parameters in the program are changed, the processing of any number of high-precision images can be realized in principle.

1. 3.

   Quantification

Digital image processing is easy to obtain quantitative results, which is an incomparable advantage of other processing methods.

1. 4.

   Flexibility

Digital image processing can not only complete linear operations but also complete nonlinear operations. All operations that can be expressed by mathematical formulas or logical expressions can be implemented by digital processing.

2.2 CT scanning principle

CT recognition technology is an advanced nondestructive testing technology to reconstruct CT images by transmitting X-rays to rotate around objects and collecting attenuation information of X-rays. The working principle of CT identification technology can be expressed as fixing X-ray source and detection receiver on the same rack, synchronously scanning them with the detected object, and scanning the scanner rack at each angle of rotation; after each scan, the scanner rack rotates to the next angle for the next scan, so repeating the process repeatedly, you can collect many sets of scan data. If the scanner frame is shifted every time and scanned once, 256 data can be obtained. Then, the scanner frame is scanned once for every rotation of 1°. If the scanner frame is rotated 180°, 256 × 180 = 46,080 data can be obtained. Finally, the scanner information is processed, and a real digital image of a scanned object can be obtained.

The mathematical basis and principle of CT nondestructive identification technology can be described as CT scanning equipment mainly consists of radioactive sources and detectors. The X-ray emitted by the source can penetrate any non-metallic material, and the penetration ability of the X-ray is different when the wavelength of the X-ray is different. The absorptive capacity of different non-metallic materials to X-ray wavelength at the same wavelength is also different. Generally, the higher the density of matter, the higher the atomic number of matter, the stronger the absorption of X-ray. When the X-ray emitted by the radiation source in the CT machine penetrates the detected object, the intensity of the X-ray attenuates because it is absorbed by the object. The attenuation follows the equation:

In the form: I₀ is X radiation which penetrates the light intensity of the detected substance (ev/m²·s), I is X which is the intensity of light passing through the detected substance (ev/m²·s), μ_m is the unit mass absorption coefficient for detected objects (cm²/g), and ρ is the density of measured substances (g/cm³); x is the X ray penetration length (cm).

Generally speaking, the absorption coefficient per unit mass of the object to be detected is related to the wavelength of the radiation emitted by the source. If the radiation X is constant, the wavelength is also fixed. Therefore, the absorption coefficient per unit mass and the density of the substance to be measured can be combined together, that is, the absorption coefficient per unit volume:

Formula: for the measured object to X-ray absorption coefficient per unit volume; for water 0.1, its absorption coefficient is M. During the scanning process, projection value P is used to record and express X radial initial ray intensity I₀ and X, the relationship between the attenuation I value of the rays passing through the measured object. The computer Xattenuation information of the ray projection value P follows the pattern:

In the form: X_i is the X of every interval of the ray path. μ_i is the local attenuation coefficient, and the total attenuation of the X rays depends on the local attenuation coefficient. μ is the line integral for scanning path. When the interval of the scanning path is accumulated, the small thickness X_i should be taken into account. CT scanning is essentially a process of accurately measuring many line integrals. Therefore, in order to obtain satisfactory scanning image quality, it is necessary to record enough projection values and measure them within 180°.

2.3 Gray histogram

Gray histogram is the probability density distribution function (PDF) of gray level, and its integral function (or area function) is the cumulative distribution function (CDF), which reflects the relationship between gray level and frequency of gray level in an image. In the gray histogram, the abscissa is the gray level, and the ordinate is the frequency (pixel number) of the gray level appearing in the image.

The mathematical expression of grayscale histogram can be written as:

The upper level S_k is the k gray level of the image f(x, y), and the number of pixels n_k in the middle of f(x, y) gray value k. n is the total number of pixels in the image. From the definition, an estimate of the frequency of S_k occurrence is given by P(S_k), so the histogram shows the distribution of the gray value in the image.

CT scanning image is a gray image composed of a large number of projection data to solve the density distribution function f(x, y) of the detected material. Essentially, it is a gray image composed of CT numbers reflecting the change of the density inside the rock mass. The distribution of different materials inside the rock mass can be reproduced in the image through the change of the gray value of the image. The evolution of meso-damage of rock and soil during loading process will cause the change of internal structure of rock and soil, and the change of gray value appears on CT image. Therefore, we can study the evolution process of rock and soil damage based on the gray histogram of image.

The gray scale of the image, the average brightness, and contrast of the whole image are the important basis for further processing. At the initial stage of loading, the gray distribution of CT images is relatively narrow, causing the details of the image to be unclear. Gray histogram correction technology can make the gray distance of the image open or make the gray distribution uniform, thereby increasing the contrast and achieving the purpose of image enhancement. The principle is to normalize the histogram first and normalize the original gray scale [Z₁, Z_k]; it is [0, 1]. Set any one of the grayscale Z normalization to r, any gray level Z^' of the transformed image is normalized to S. r, s satisfy 0 ≤ r ≤ 1, 0 ≤ s ≤ 1.

Histogram correction is the calculation of the following formula:

(6)

Formula: T(r) is a transformation function; the following conditions must be satisfied:

1. 1.

   T(r) is in the interval of 0 ≤ r ≤ 1; T(r) is a single valued function and monotonically increases.

2. 2.

   T(r) in 0 ≤ r ≤ 1 satisfies 0 ≤ T(r) ≤ 1.

The first condition guarantees that the order of the gray level of the image from white to black remains unchanged. The second condition guarantees that the gray value of the pixel after mapping is within the allowable range.

Gray level equalization is the abbreviation of gray level histogram equalization. It is a kind of histogram correction method. By calculating every pixel in the image, the input image is transformed into the output image with the same number of pixels at each gray level, that is, the histogram of the original image is transformed into the form of uniform distribution. This enlarges the dynamic range of pixel gray value, highlights the details that are not clearly seen in the original image, and thus achieves the effect of enhancing the overall contrast of the image. Probability theory has proved that any cumulative distribution function (CDF) obeys normalized uniform distribution. This conclusion makes the gray equalization algorithm evolve into integral form.

Formula: h(x) is the gray frequency function A for the image area, that is, the total pixel L is the maximum gray level, generally take 255. For discrete images, the upper form evolves into the following summation form:

The general steps of grayscale equalization are:

1. 1.

   The gray level S_k(k = 0, 1,  … , L − 1) of the original image is counted, and the number of pixels n_k of each gray level is calculated.

2. 2.

   The gray histogram of the original image is calculated by \( P\left({S}_k\right)=\frac{n_k}{n}\kern1.25em k=0,1,\cdots L-1 \)

3. 3.

   Calculate the cumulative histogram of the original image.

\( {t}_k=\sum \limits_{i=0}^k\frac{n_i}{n}=\sum \limits_{i=0}^k{P}_s\left({S}_i\right) \) 0 ≤ S_k ≤ 1 k = 0, 1, … , L − 1

1. 4.

   Rounding calculation: \( {t}_k=\operatorname{int}\left[\left(L-1\right){t}_k+\frac{k}{L}\right] \)

2. 5.

   Mapping relations: S_k → t_k

3. 6.

   Count the number of pixels n_k of the new histogram gray level and calculate the new histogram \( P\left({t}_k\right)=\frac{n_k}{n} \)

2.4 Image edge processing

The edge of an image is the reflection of discontinuous gray change and texture structure change of local characteristics. It marks the end of one region and the beginning of another region. It is the collection of pixels with step change or roof change in the gray level of surrounding pixels. Edge is the place where the gray level changes most dramatically. Traditional edge detection is to use this feature to differentiate or calculate the second-order differential of each pixel to determine the edge pixels. The peak value of the first-order differential image corresponds to the edge point of the image, and the zero-crossing point of the second-order differential image corresponds to the edge point of the image. For simple first-order derivative operation of image, it can only detect edges in specific directions because of its fixed orientation, so it is not universal. In order to overcome the disadvantage of first derivative, image gradient operator is used to calculate the image. The gradient corresponds to the information of the first derivative, and the gradient operator is the first derivative operator. For a continuous image f(x, y), the gradient at the point (x, y) can be expressed as a vector. Suppose G_x, G_y, it is used to represent f(x, y), along the direction x and direction y gradient, and the gradient vector can be expressed as:

The gradient direction is the fastest change direction of the image gray value. θ_g indicates gradient direction:

The rate of change in direction θ_g (the magnitude of the gradient) is:

g is a gradient operator, and the amplitude calculation is modular to 2 (corresponding to Euclidean distance). The equivalent norm can also be used.

Or the infinite norm is used:

According to the actual situation in the field of digital image processing to use discrete form, the above differential difference is used instead of:

The second derivative information is a sign of the change of the first derivative. In some cases, if the gray level changes uniformly, the boundary may not be found by using the first derivative only, then the second derivative can provide useful information. The two derivative information is generally based on formulas in the field of digital image processing.

Based on the first derivative operator, this kind of operator mainly detects the change of the edge gray level feature. There are several operators which use the first derivative extremum detection as follows:

(a) Roberts operator: A template as shown in Fig. 1. Edge detection based on this operator is generally very sensitive to noise.

Roberts operator template

Full size image

(b) Sobel operator: A template as shown in Fig. 2 is convoluted vertically and horizontally to obtain luminance difference. Different pixel positions are weighted. Generally, edge detection based on this operator is not accurate.

Sobel operator template

Full size image

(c) Prewitt operator: Two templates, as shown in Fig. 3, are used to convolute the two templates, the upper and lower neighbors in the vertical direction and the left and right neighbors in the horizontal direction.

Prewitt operator template

Full size image

Canny algorithm is used in many operations in this paper. Canny operator is a multi-scale spatial edge detection operator. It determines the image edge by finding the local maximum of the image gradient. The Canny operator can suppress noise, and there will not be many isolated edge pixels. The specific steps are as follows:

1. 1.

   First, two dimensional Gauss filter templates are used to filter the image to eliminate noise.

2. 2.

   Find the sum of the partial derivatives G_x and G_y along the two directions of the image gray level by using the differential operator and get the direction of the gradient.

3. 3.

   The gradient direction of the edge is roughly divided into four kinds (horizontal, vertical, 45°, and 135° directions). Each direction is compared with different neighboring pixels to determine local maxima. If the gray value of a pixel is not the largest compared with the gray value of the two pixels before and after the gradient direction, then the pixel value is set to 0, that is, not the edge.

4. 4.

   Use histogram to calculate the two thresholds. If the detection result is larger than the low threshold but smaller than the high threshold, it is necessary to see if there are edge pixels in the neighborhood pixels of this pixel that exceed the high threshold: if there is, it is the edge; otherwise, it is not the edge.

Comparing the edge detection results of several existing operators, Roberts operator can only roughly detect the outline, because this operator is very sensitive to gradient, so it produces some isolated points, which make the edge connection discontinuous and more discontinuous; Sobel operator and Prewitt operator do not detect the edge in some directions, and there is a considerable length of the break. This is because the two operators are not isotropic; Laplacian operator is a linear, isotropic second-order differential operator; the detection effect is relatively clear, but the detection results of the boundary are too wide, and the positioning accuracy of the edge is not high. Canny operator suppresses noise, and there are no isolated edge pixels. The location is accurate and the detection effect is the best. Therefore, the Canny operator is used to realize the edge detection of frozen soil CT images.

Due to the limitation of the experimental conditions, the size of the specimen is 2550 mm. During the processing, some large aggregate will be cut off. In order to understand the actual size distribution of the rock and soil inside the specimen, the rock and soil in the CT image of the specimen in the initial state are identified and analyzed statistically. The purpose of determining the basic dimensions of research objects is to understand the nature of the research questions. Fracture of geotechnical materials is an evolutionary problem, and its mechanism cannot be understood thoroughly only by macroscopic phenomenological method. At the same time, the microscopic or microscopic mechanisms of many complex phenomena are much simpler, and their combination and interaction can show complex behavior. Through the CT scanning experiment, the CT scanning images (two-dimensional images) of a certain scanning plane of rock and soil can be obtained. A set of parallel two-dimensional scanning CT images can be obtained by multiple parallel scanning, and the cross-sectional images of meso-structure of frozen rock and soil in different positions can be obtained. These two-dimensional scanning CT images can be observed. Three-dimensional reconstruction of fault sequence images is realized by MATLAB. MATLAB software provides a variety of matrix operations, image manipulation, image display functions, and image processing tools and has been widely used in image processing. MATLAB has many functions that can be used for image processing. MATLAB’s powerful image processing function makes image processing more convenient. This paper uses the regionprops function in MATLAB to obtain the pore structure parameters, such as equal area circle radius, shape center coordinates, equal area ellipse long and short axis, and long axis angle. Because the maximum value of the pore shape coefficient is 0.87, the minimum value is 0.64, the average value is 0.72, and the pore shape is more regular, using equal area circle respectively. The two models were meshed with the same accuracy, and the number of grids was 1207 and 1202 respectively.

Heterogeneous distribution of aggregates, cement slurry, and pores in rock and soil forms different speckle structures. The non-linearity of macroscopic mechanical behavior of rock and soil under load is due to the heterogeneity of its meso-structure. In the process of non-linear evolution, the disorder of the internal structure of rock and soil may be strongly amplified and ultimately show significant macroscopic differences. Under basically the same experimental and loading conditions, the real meso-structure of rock samples and the mechanical properties of various meso-media play a decisive role in the stress and deformation response of the specimen. The CT image of porous materials in Fig. 4 has a typical structural diagram, which can be directly set threshold value for segmentation and extraction of pores, the original gray scale image is transformed into binary image, and the image is scanned by CT. Different geotechnical materials are shown in Fig. 4, and the CT scanning image is obtained.

Different geotechnical materials and CT scanning images

Full size image

The size of each voxel in the CT image is 40 × 40 × 40 μm³, and there are still micro-pores in the voxel. The gray value of each voxel reflects its internal structure. When the fracture or pore size is smaller than the spatial resolution, that is, the micro-crack or pore exists in the voxel, this part of the pore cannot be distinguished on the pixel scale, but these micro-cracks or pore changes will lead to the change of voxel density.

It can be seen that there is a large number of micropores smaller than CT resolution in voxels. In the above statistical analysis of gray value, before the appearance of cracks larger than CT resolution scale, the change of gray value is mainly caused by micro-cracks or pores in voxels. It can be seen that the interface between aggregate and matrix is easy to be destroyed with the increase of load when there is a large stress value at the interface between aggregate and matrix. For the above porous materials, vectorization method can also be used to model, but if the extracted pore edge is not smoothed, it will lead to excessive number of elements and increase the amount of calculation.

The edge image is vectorized by software. As shown in Fig. 5, the thresholds are selected to segment the aggregate and pore regions. Fracture or pore size extracted from CT images should be larger than the spatial resolution of CT, so the extraction of fracture or pore is more limited by the spatial resolution of the system.

Image processing of different geotechnical images

Full size image

Because the gray value between the crack or pore and the matrix has a certain transition, it is not a sharp boundary, and the boundary is not just right for the whole pixel, as well as part of the volume effect; when the threshold is selected for image segmentation, there will always be some voxels, their porosity is high, resulting in the gray value to reach the threshold value. It is divided into pores larger than CT resolution scale. If a single threshold is chosen, some crack details may be lost. In order to extract cracks more accurately, the segmentation of images with different sizes or in different regions can be carried out by using multi-threshold method. However, extracting cracks from a series of three-dimensional CT images will be a time-consuming task. The larger the fracture size, the easier it is to extract. Based on this, the method of fracture extraction is more suitable for analysis of cracks larger than CT resolution scale. In addition, it should be noted that the porosity calculated from CT images is less than that measured by mercury injection and helium gas injection due to the limitation of the resolution of CT images.

At this time, the numerical model can be established by vectorization of the feature area edge. The grayscale histogram is shown in Fig. 5.

According to Fig. 6, the meso-structure image of rock and soil is a binary image, and the closed interface between cracks, cavities, and rock and soil is formed by the pixels of gray degree. Since there is no finite element software to accept the image as input data, if the data of geotechnical micro-structure image is input by a pixel, according to the principle of finite element mesh generation, each pixel will become a node of finite element, so that the generated mesh will be too dense and the finite element model will be generated. This is very difficult and will lead to the failure of further finite element numerical analysis. Therefore, mesoscopic structure images can not directly generate finite element meshes. In order to realize the next step of mechanical analysis, the discrete pixel data of mesoscopic structure images must be transformed into geometric vector information to realize the generation of finite element meshes.

Histogram

Full size image

Taking rock and soil as the research object, combining with micromechanics, experimental mechanics, statistical theory, industrial CT technology, and finite element numerical simulation method, the evolution process of internal damage of rock and soil under uniaxial compression is studied. The generation and development process of strain localization before failure is analyzed by strain field. The main work and understanding are as follows:

1. The damage variable and damage constitutive model of concrete are further studied. Based on the meso-test of damage evolution characteristics of concrete, different damage variables of concrete materials are defined by different methods, and corresponding damage evolution equations are established. The damage constitutive relation is established on the basis of the piecewise damage evolution equation established by the number change. The damage evolution equation and damage constitutive relation established by the change of density damage number are related to the deterioration of macroscopic mechanical properties.

2. The difference of gray mean and mean square deviation of CT images of different layers under different stress states reflects the inhomogeneity of the specimen. With the increase of load, the change of specimen volume causes the change of gray value. The gray mean value of the whole image goes through a process of increasing first and then decreasing, while the gray variance decreases first and then increases correspondingly. Both of them change nonlinearly with the deformation of specimen. The percentage distribution of cracks in the specimen under different stress states was obtained, and the change of micro-cracks in voxels at the interface between aggregate and cement slurry was obvious.

3. The larger the fracture size, the easier it is to extract. Based on this, the method of fracture extraction is more suitable for analysis of cracks larger than CT resolution scale. In addition, it should be noted that the porosity calculated from CT images is less than that measured by mercury injection and helium gas injection due to the limitation of the resolution of CT images.

CDF:

Cumulative distribution function

FEM:

Finite element method

PDF:

Probability density distribution function

Not applicable.

About the authors

Gang Luo (corresponding author, e-mail:20044592@163.com) received the Ph.D. degree in tunnel and underground engineering from Southwest Jiaotong University, Chengdu, China. He is currently a teacher in the School of Highway, Chang’an University, Xi’an, China. His research interests include submerged floating tunnel, controlled blasting technology and tunnel engineering.

Shaokang Pan received the bachelor’s degree in civil engineering from Yangtze University, Jingzhou, China. He is graduate student in the School of Highway, Chang’an University, Xi’an, China. His research direction is tunnel engineering.

Yulong Zhang received the bachelor’s degree in civil engineering from Huaiyin Institute of Technology, Jiangsu, China. He is graduate student in the School of Highway, Chang’an University, Xi’an, China. His research direction is tunnel engineering.

Hanghang Jia received the bachelor’s degree in civil engineering from Tianjin Chengjian University, Tianjin, China. He is graduate student in the School of Highway, Chang’an University, Xi’an, China. His research direction is tunnel engineering.

Liang Chen received the bachelor’s degree in civil engineering from Chengdu University of Technology, Chengdu, China. He received the master’s degree and Ph.D in bridge and tunnel engineering from Southwest Jiaotong University, Chengdu, China. He is currently a engineer in Guangxi Communications Investment Group Co., Ltd., Nanning, China. His research direction is tunnel engineering.

Funding

This work was supported by national natural science foundation of China (No.51708042, 51508037).

Availability of data and materials

Please contact author for data requests.

Authors and Affiliations

1. School of Highway, Chang’an University, Xi’an, Shaanxi, China

   Gang Luo, Shaokang Pan, Yulong Zhang & Hanghang Jia

2. Guangxi Communications Investment Group Co., Ltd., Nanning, Guangxi, China

   Liang Chen

Contributions

All authors take part in the discussion of the work described in this paper. Author GL wrote the first version of the paper. Authors SP, YZ, and HJ did part experiments of the paper, and LC proposed some modifications. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Gang Luo.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Reprints and permissions