X-ray Computed Tomography (X-CT) technique is one of the most advanced non-destructive testing techniques, which can obtain the internal structure of materials. The density difference between the MKP and the pores leads to the difference of X-ray absorption coefficient and gray value in the X-CT images. Therefore, the two materials can be distinguished according to the gray value division (Sun et al., 2014). The X-ray projections were obtained with an exposure time of 0.32 s at an accelerating voltage of 150 kV and 140 μA beam current using a tungsten target.

In nanoindentation test, the elastic modulus of the specimens was measured by nanoindentation instrument manufactured by Agilent Company, USA with the Berkovich indenter (a positive triangular pyramid with an angle of 65.35° between the center line and the side). The appropriate loading depth and indentation spacing should be selected for nanoindentation test (Al-Amoudi, 2002; Němeček et al., 2009; Wang et al., 2009; Chen et al., 2010). Since in this study the nanoindentation test was simulated, a larger indentation depth was required to match the X-CT resolution and obtain a sufficient number of finite element meshes. The maximum indentation depth of the nanoindentation instrument used in this study is 20 μm. Hence, the indentation depth was set to 20 μm and the indentation point spacing was set to 200 μm.

Figure 1A displays a typical Scanning Electron Microscope (SEM) image of the nanoindentation. However, the indentation area includes not only the matrix MKP but also the pores, shown in Figure 1B. Therefore, the material properties, volume fractions and distribution characteristics of the matrix and pores both affect the test results of the nanoindentation test. The nanoindentation process can be divided into three stages: first, the elastic deformation of the specimen occurred at the loading stage, followed by the plastic deformation with the increase of the load. Then the constant loading stage appeared when the maximum indentation depth was reached. Finally, the unloading stage took place, reflecting the elastic recovery of the indentation points. The typical load-displacement curve of the nanoindentation is shown in Figure 1C. The curve consists of three parts: the loading stage, the constant loading and the unloading stage. The loading time of the loading stage was 1,000 s. When the indentation depth reached 20 μm, the constant loading kept for 100 s, followed by the unloading stage for 300 s. The contact stiffness S is fitted by the upper elastic part of the unloading curve. According to the Oliver-Pharr principle (Oliver and Pharr, 2011), the elastic modulus of the indentation point can be calculated by Equation (1).

where E and ν indicate the elastic modulus and Poisson's ratio of tested material, γ is a correction factor (γ = 1.034), S is the contact stiffness, A is the contact area, ν_i and E_i denote the parameters of indenter (ν_i = 0.07, E_i = 1,141 GPa).

In the mesoscopic numerical analysis of cement-based materials, it is important to study the numerical morphology and gradation of aggregates, pores and other phases, which directly affect the mechanical properties of the materials (Wang et al., 2016).

In order to improve the accuracy of numerical simulation in material mesomechanics, a 3D random concave-convex aggregate modeling method with grid pre-generated was proposed (Wang B. et al., 2018). Before placing the aggregate, grid partition was needed to record all node information and unit information in the model. In the three-dimensional polar coordinate system, the position of any point in space could be determined by three parameters r, θ, φ. It is impossible to use infinite number of spherical aggregate surface nodes in the numerical modeling. It is found through trial calculation that the spherical aggregate established by the following method could meet the calculation requirements: take r as the aggregate radius and take a point, respectively, at every 45° in the direction of θ and φ, forming an approximate spherical area composed of 26 points in the space. The three-dimensional spherical aggregate had a total of 26 nodes and the surface was divided into 48 triangular regions. The generation of the concave-convex aggregate could be achieved when the apex of each aggregate fluctuates randomly with respect to its initial position. After the single concave-convex aggregate was generated, the aggregate library could be generated according to the required gradation.

The flow path of placing aggregate was as follows (Wang B. et al., 2017): firstly, the aggregates to be placed were sorted according to the radius from large to small. Secondly, the central point of the non-throwing unit was selected to place aggregate based on whether there was a geometric boundary point inside the aggregate. If there were any, the placing failed. Thirdly, if the interior of the aggregate did not contain the geometric boundary point and the outer boundary element point of the delivered aggregate, the delivery was successful. Fourthly, when the delivery failed, the aggregate was rotated to continue judging and the position was reselected after a certain number of rotations. If all the delivery failed in the center position of all the non-throwing areas, the aggregate failed to be delivered. In order to prevent the delivery process from entering an infinite loop, the failed aggregate was stored in the specified set. Fifthly, the aggregate information was updated after delivering successfully. In this study, when the random aggregate placement method was adopted to establish the MKP nanoindentation model, the pores were regarded as the aggregates and put into the MKP matrix.

Since the elastic modulus of the diamond indenter was much larger than the elastic modulus of the MKP, the deformation of the indenter during the indentation process was negligible. Consequently, the indenter in the nanoindentation model was set as a rigid body without deformation and was meshed with the mesh size of 5 μm. During the contact between the indenter and the MKP, the contact surface was assumed to be smooth. Material properties were the extremely critical input parameters for the finite element model. Material parameters were not required for the indenter as a rigid body. The material properties of the MKP were verified in the previous paper: the density was 1.864 g/cm³, the elastic modulus was 37.3 GPa, the Poisson's ratio was 0.2, and the compressive strength and tensile strength were 75.6 and 11.3 MPa, respectively (Li et al., 2018a). The above MKP material parameters were input into the nanoindentation model.

In the X-CT nanoindentation model, the six degrees of freedom of all nodes on the bottom surface of the MKP substrate were limited to zero, indicating that it was completely fixed, while the nodes of the other five surfaces were not limited. Five degrees of freedom of the indenter was restricted except Z direction, so they could only move in the directions of loading and unloading. Besides, the tip of the rigid indenter was regarded as a reference point applied with a displacement of about 20 μm. The model was solved by ABAQUS dynamic analysis with the total running time of 1.4 × 10⁻⁶. The loading time-amplitude conformed to the following rules. When the time was between 0 s and 10⁻⁶, the indenter was in the loading stage with the displacement from 0 to 20 μm. When the time was in the range of 10⁻⁶ to 1.1 × 10⁻⁶, the indenter was in the constant loading without movement. During the period of 1.1 × 10⁻⁶ to 1.4 × 10⁻⁶, the indenter was in the unloading phase with the displacement reducing from 20 to 16.2 μm. The relationship between amplitude and time is as follows:

when t ≤ 10⁻⁶,

when 10⁻⁶ < t < 1.1 × 10⁻⁶,

when 1.1 × 10⁻⁶ < t < 1.4 × 10⁻⁶,

where A is the amplitude, t represents time.

The MKP constitutive relation was plastic damage model and determined in previous studies (Li et al., 2018a) as follows:

when ε < 0.4ε_c,

when 0.4ε_c ≤ ε ≤ ε_c,

when ε > ε_c,

where subscript c indicates “compressive,” ε_c is the corresponding peak compressive strain, f_c denotes the compressive strength of the MKP, namely 75.6 MPa.

when ε ≤ ε_t,

when ε > ε_t,

where subscript t indicates “tensile,” ε_t is the corresponding peak tensile strain, f_t denotes the peak tensile strength of the MKP, namely 11.3 MPa.

According to the literature (Birtel and Mark, 2006), the classical damage factor can be calculated from Equations (12)–(14):

where d indicates the damage factor of the MKP, σ is the stress, E represents the elastic modulus, εⁱⁿ and ε^pl denote the inelastic strain and the plastic strain, b is the constant (b_c = 0.7, b_t = 0.1), subscript c and t indicate “compressive” and “tensile,” respectively.

Since the elastic modulus of MKP is 37,300 MPa, the compressive and tensile damage factor of the MKP are expressed as:

where the relationship between σ and ε can be calculated from Equations (6)–(11) in section The Constitutive Relation of the MKP. Thus, the curve of the classical damage factor is shown as the blue curves in Figure 4.

The above boundary conditions, MKP material parameters, constitutive relations and damage factors were input into the X-CT nanoindentation model of the P₁₄ indentation point. Then the simulated load-displacement curve based on the classical damage factor were expressed as the red curve in Figure 5, in which the simulated peak load is 2.25 N. The blue curve in Figure 5 indicates the test load-displacement curve with a peak load of 2.04 N. The relative difference was used to compare the simulated peak loads with the experimental peak loads, defined as

where L_E indicates the experimental peak load, L_S is the simulated peak load. It can be seen that the relative error of the peak load of the classical simulation is 10.3% and the classical simulated load-displacement curve is above the experimental curve. The poor simulation result was because the classical simulation did not consider the damage contributed by the pores. Since the matrix in the nanoindentation model contained two phases: the MKP and pores, the pores inevitably affected the simulation results, so the influence coefficient of the pores on the damage factor was introduced as χ. Then the modified damage factor considering the influence of the pores was expressed as

The influence coefficient of P₁₄ indentation point was debugged. It was found that the simulation result with the influence coefficient χ 1.262 was the closest to the experimental result. The modified damage factor with the influence coefficient of 1.262 is shown as the purple curve in Figure 4 and the yellow curve of Figure 5 displays the modified simulated load-displacement curve. The modified simulated peak load was 2.11 N with the relative error of only 3.4%. The modified simulation accuracy was greatly improved, indicating that it was necessary to consider the influence of pores on the damage. The modified simulated nephogram is shown in Figure 6.

Figure 6A shows the simulated stress nephogram of the model at the time 10⁻⁶. It can be seen that the shape of the indentation area is a triangular pyramid and the mesh deformity in contact with the three edges of the indenter is severe due to the stress concentration. The stress at the bottom of the indentation area is larger and the upward stress decreases gradually. As a result of the extrusion of the indenter, the upper end of the indentation region produces a crowding effect. The stress distribution of the middle section in the X direction is displayed in Figure 6B at time 0 and 10⁻⁶, respectively. It can be seen that the stress in the lower part of the indenter is larger and gradually decreases toward the far side. The stress nephogram shows an elliptic shape and extends outward in layers. During the loading process, the pores in the lower region of the indenter were squeezed and seriously deformed, such as the red marked pores being extruded from the initial quasi-circular into a long strip. In addition, the existence of pores may result in stress concentrations, just as the sudden increase of stress around the pore marked with yellow. Moreover, pores may affect the path of stress transfer. For instance the pink marked pore hindered the development of stress. Figure 6C is the distribution of the compressive damage of the model when the time is 10⁻⁶. It is observed that the element damage in the larger area around the indenter is serious, that is to say, the pressing of the indenter greatly affects the stress state of the surrounding area. Figure 6D displays the distribution of the compressive damage in the middle section of the Y direction at time 0 and 10⁻⁶, respectively. It can also be seen that the loading of the indenter leads to severe distortion of the pores, such as the green marked pores being significantly flattened. In addition, the development of damage may be affected by the pores, such as the damage transmission path in pink marked area was affected by the pores.

It was the same as the process of debugging the influence coefficient of P14 indentation point with the influence coefficient of 1.262 in the fourth paragraph of this subsection. The influence coefficient of the porosity of the other indentation point model on the damage factor was debugged, thereby making the simulation results most consistent with the experimental results. Thus, the optimal influence coefficient values and the optimal simulated peak loads of each indentation point were obtained as shown in Table 3. There are two data at each indentation point, the optimal influence coefficient values and the optimal simulated peak loads, respectively. For example, the data “1.248, 1.84” for the indentation point P₁₁ indicates that the optimal influence coefficient value and the optimal simulated peak load at the indentation point P₁₁ are 1.248 and 1.84 N, respectively. The average relative error of the peak load simulated by the X-CT nanoindentation model was 3.8%, which indicates that the accuracy of the model was greatly improved after debugging the optimal influence coefficient.

The strength of the material was affected not only by the porosity, but also by the pore distribution. The pore structure of the material was complex and the pore size reflected the dominant characteristics of the pore structure (Hou et al., 2019). Therefore, the pore size was selected as an important factor affecting the damage factor. Considering the distribution and contribution of pore size, the effect of pore distribution on the damage factor was defined as λ.

First, the actual pore size distribution of the material was complex. Previous studies indicated that the pore distribution functions of cement-based materials and rock masses were similar, which all had an exponential function distribution (Ju et al., 2008). It was also assumed that the pore size distribution of MKP conformed to the exponential function distribution, as shown in Equation (23):

where D indicates the pore size (mm), E and F represent the parameter of the exponential function. The values of the parameters E and F can be solved by fitting the hole structure results based on the X-CT statistics.

Secondly, gray correlation analysis was a good method to reflect the correlation between various factors, which could be used to study the influence of pore size on strength. The correlation value had a power function relation with pore size. The larger the pore size was, the lower the strength of the material was Zhang and Zhang (2007); Li et al. (2018b). In addition, the influence coefficient was zero when the pore size was zero. Therefore, the influence function of pore size was assumed to be:

where λ(D) is the influence coefficient associated with the pore size, G and H denote the parameter of the influence function. By fitting the relationship between the strength and the test results of pore structure with gray correlation analysis, the values of the parameters G and H can be determined.

Combining the MKP pore distribution function with the pore size influence function, the influence coefficient of the pore size distribution on the damage factor was obtained as follows:

where λ is the influence coefficient depending on the pore size distribution, D_min and D_max indicate the minimum diameter and the maximum diameter (mm) of the pores, respectively.

Due to the complexity of the integration process, the Equation (25) was substituted by the Equation (26) for the sake of simplicity:

where c_i and λ_i represent the porosity and the influence coefficient of the i-th interval when the pore size range is divided into n intervals. c_i is related to the pore size distribution, namely, ci=∑DiDi+1p(D)dD. The value of λ_i is the influence coefficient of the pore median of the i-th interval. c indicates the total porosity, namely, c=∑i = 1nci. The smaller the interval is, the closer the Equation (26) is based on the Equation (25).

According to the porosity results based on X-CT statistics, the porosity of the MKP nanoindentation area was between 3 and 23%. Therefore, the maximum porosity in this study was considered as 30%. According to the pore distribution curve of P₁₄ in Figure 2D, the pore distribution function of the P₁₄ indentation point was fitted as shown in the Equation (27). The yellow curve in Figure 2D displays the pore distribution function curve.

According to the Equation (27), the pore size range was divided into nine intervals, as listed in Table 4. The grading porosity indicates the percentage of the pore volume fraction in this interval. The regression values of the parameters G and H in Equation (21) were 2.4 and 0.18, respectively, by fitting the relationship between the load and the pore distribution of each indentation point. The parameters λ₁, λ₂, λ₃, λ₄, λ₅, λ₆, λ₇, λ₈, and λ₉ were determined as 0.745, 0.878, 0.976, 1.052, 1.102, 1.137, 1.170, 1.207, and 1.253, respectively, according to the Equation (26). Thus, the influence coefficient λ of the pore distribution was 1.06. Because the total pore influence coefficient χ of P₁₄ was 1.262, the influence coefficient β of porosity on damage was 1.191 according to the Equation (21). Similarly, the influence coefficient of the porosity of each indentation point could be obtained as Table 3. Figure 7 displays the scatter plot about the coefficient β and the porosity. Then the relationship between the porosity and the influence coefficient β was gained by fitting the scatter plot:

where c denotes the porosity, R² represents the correlation coefficient of the fitted equation. The R²-value is 0.98, which means that the fitted curve has high accuracy.

Considering the influence of porosity and pore size on the damage factor, the modified MKP damage factor is expressed as follows:

The damage factors calculated by Equations (29) and (30) were input into the X-CT nanoindentation model and the three indentation points of P₁₅, P₂₅, and P₃₅ were simulated, respectively, to verify the accuracy of the modified damage factor equation. Table 4 shows some input parameters and simulation results. The simulated peak loads of the three indentation points were 1.90, 1.89, and 1.79 N with the relative errors of 3.1, 2.4, and 4.2%, respectively. The average relative error was 3.2%. Thus, the assumptions that the relationship between the influence coefficient β and λ, the MKP pores conforms to the exponential distribution and the influence function of the pore size distribution were correct. That is, the modified damage factor equation considering the pores influence was valid.