Experimental-Computational Approach to Investigate Nanoindentation of Magnesium Potassium Phosphate Hexahydrate (MKP) With X-CT Technique and Finite Element Analysis

The magnesium phosphate cement (MPC) is a carbon-free cementitious material, widely used in solidification of nuclear waste, heavy metals, and repair and reinforcement. The magnesium potassium phosphate hexahydrate (MgKPO4·6H2O, MKP) is the main hydration product of MPC, seriously affecting the mechanical properties of the MPC. Therefore, this paper presented an experimental-computational approach to investigate the mechanical properties of the MKP through nanoindentation with X-ray Computed Tomography (X-CT) technique and finite element analysis. Firstly, the micro-mechanical properties and structural distribution characteristics of the MKP were tested based on the nanoindentation and the X-CT technique, respectively. Then, the 3D structure grid model of the MKP was obtained based on X-CT data, imported into the ABAQUS software for the finite element simulation. Besides, considering the effect of porosity and pore distribution on the damage, the modified MKP constitutive relation was proposed and input into the X-CT nanoindentation model and the RAP nanoindentation model, respectively. It was found that those two models can effectively describe the mechanical and deformation characteristics of the MKP, which verified the correctness of the modified constitutive relationship of MKP. Finally, the influence of pore distribution on the nanoindentation results was predicted based on the RAP nanoindentation model.


INTRODUCTION
The production of widely used Portland cement consumes a lot of energy and emits a large amount of carbon dioxide, causing serious environmental pollution. Green building materials are urgently needed in today's world to avoid the deterioration of ecosystems and the intensification of global warming. Magnesium phosphate cement (MPC) is a carbon-free cement that does not emit carbon dioxide during the production process (Walling and Provis, 2016). MPC is considered as a new type of environmentally friendly cement (Haque and Chen, 2019). Additionally, it is believed as a cementitious material formed by acid-base chemical reaction between water, magnesium oxide, and phosphate. The main chemical reaction equation of potassium MPC is shown in Equation (1): MgO + KH 2 PO 4 + 5H 2 O→MgKPO 4 ·6H 2 O where potassium hexahydrate hexahydrate (MKP) is the main hydration product of the MPC (Vinokurov et al., 2018a).
Furthermore, MPC is widely used in solidification of nuclear waste, heavy metals and repair and reinforcement due to its advantages of fast hardening, early high strength, and high viscosity He et al., 2019;Mestres et al., 2019;Zhenghua et al., 2019), which attracts many scholars to do enormous research on hydration mechanism, durability, and mechanical properties (Ma et al., 2014;Li et al., 2015Li et al., , 2016.
Liquid radioactive waste (LRW) containing elements such as uranium and thorium is the product of nuclear industry activities which has great harm to the environment and human health. Long-term controlled storage or disposal of LRW is one of the key links in the safe management of radioactive waste (Stefanovsky et al., 2017). Solidification/stabilization (S/S) is the mixing of LRW and binder, thus fixing the LRW in the binder for long-term safe disposal. The S/S is a very effective technique for handling large amounts of the LRW and the MPC is one of the most promising materials for solidifying the radioactive waste (Vinokurov et al., 2019). Four different leaching tests were conducted to determine the effect of the MPC on solidifying heavy metal-containing waste liquid. The results showed that the MPC could successfully solidify the solution containing Cd, Cr, Cu, N, Pb, or Zn (Buj et al., 2010). The radioactive waste containing metal uranium is incompatible with conventional ordinary Portland cement (OPC)-based encapsulation matrices. The reason is that the high alkaline environment and high free water content in the OPC result in volume change of the system and hydrogen generation. Since the MPC has lower pH value and less free water content, it can encapsulate the radioactive waste containing active metals such as uranium (Covill et al., 2011). In addition, the researchers found that the MPC can solidify radioactive waste within actinide and rare earth elements. The method of the MPC curing the LRW has the advantages of simple technology and high physical and chemical stability (Vinokurov et al., 2009(Vinokurov et al., , 2018b. Heavy metals such as Zn, Pb, Cd, As contaminate soils, which are also increasingly serious environmental problems. Heavy metals not only pose a threat to the human health and the environment, but also deteriorate the mechanical properties of the soil, limit the reuse of contaminated sites, and even pose a safety threat to engineering in the polluted areas (Du et al., 2014). As mentioned above, the S/S technique is an effective and economical remediation technology. It mixes the binder MPC with contaminated soil to physically fix or chemically bind harmful contaminants, thus preventing heavy metals from migrating to the environment and enhancing the strength of the soil (Wang Y.S. et al., 2018;Xu et al., 2018). Aiming at the environmental problem of high lead content in soil polluted by lead-acid batteries, researchers found that the MPC could effectively fix lead in soil and convert lead in polluted areas into less mobile contaminant binding forms at an acceptable cost, effectively increasing the strength of polluted soil . The factors affecting the solidification effect of the MPC are metal concentration in soil and water-binder ratio of the MPC, curing age and dosage (Wang P. et al., 2017(Wang P. et al., , 2018. In addition, municipal solid waste incineration (MSWI) fly ash pollution is highly toxic, threatening human living environment and health. The MPC can also effectively reduce the toxic pollution of the MSWI fly ash (Fan et al., 2018).
Energy demand and resource extraction activities are one of the major environmental concerns of modern society. The MPC has a good repair and reinforcement effect on concrete structures due to its characteristics of fast hardening and early high strength. The MPC can be injected into cracks from old buildings or bonded with CFRP to reinforce damaged structures, thus delaying the demolition of old buildings and the construction of new buildings, reducing energy, and resource consumption . In addition, the MPC can be produced from industrial by-products, which has a positive impact on the environment and sustainability (Maldonado-Alameda et al., 2017). For example, boron-containing magnesium oxide (B-MgO) was a byproduct of the production of Li 2 CO 3 from salt lakes. B-MgO can be used as raw materials to produce MPC (Tan et al., 2014;Formosa et al., 2015). The MPC even can upcycle construction wood waste into rapidshaping cementbonded particleboards, reducing environmental burden and increasing economic value (Wang L. et al., 2017).
In summary, the mechanical properties of the MPC have an important influence on the solidification of heavy metal contaminated soil and the repair and reinforcement of engineering structures. The main hydration product MKP determines the mechanical properties of the MPC. However, most researches only focus on the application and macromechanical properties of the MPC, which do not involve the mechanical properties of the MKP. Therefore, in this paper, the micro-mechanical properties and structural characteristics of the MKP are studied in detail.
Nanoindentation technique and X-ray Computed Tomography (X-CT) technique are the state of the art technologies, which can measure the micro-mechanical properties and 3D structure distribution characteristics of materials. The micro-mechanical properties and creep behavior of cement paste were attained by the nanoindentation technique. The factors affecting the test results of nanoindentation and the contact creep function of various phases in cement paste were determined (Liang et al., 2017a,b;Wei et al., 2017). CT technology was used for in-situ monitoring of water and ion intrusion into cement paste and the erosion process was visualized in three dimensions (Yang et al., 2015(Yang et al., , 2018. Based on the CT images, pore-scale modeling and micromechanical modeling of cement paste could be obtained Zhang, 2017). The combination of X-CT and finite element method can more accurately simulate the properties of the materials (Skarzynski and Tejchman, 2016). Therefore, based on nanoindentation and X-CT techniques, the micro-mechanical properties and structural distribution characteristics of the MKP were investigated. The nanoindentation model of the MKP was established based on the X-CT and the random aggregate placement method. From the perspective of damage factor, the modified constitutive relation of the MKP considering porosity and pore distribution was proposed. Firstly, the 3D microstructure distribution characteristics of the MKP were obtained based on the X-CT technique. The MKP was tested by the nanoindentation to attain the elastic modulus and the indentation load-displacement curves. Then based on the X-CT images, the MKP 3D structure grid model was obtained and imported into ABAQUS software as the nanoindentation model. In addition, the influence equation of pores on the damage factor in MKP constitutive was proposed. The effect of porosity and pore distribution on the MKP damage was studied in detail and the modified MKP constitutive relation considering porosity and pore distribution was verified. Finally, based on the random aggregate placement method (RAP), the RAP nanoindentation model of the MKP was established and the effect of pore distribution on the nanoindentation results was predicted.

Materials and Specimen Preparation
The experimental material in this research was magnesium potassium phosphate hexahydrate (MgKPO 4 ·6H 2 O, MKP). The MKP granules were obtained from a pharmaceutical reagent factory in Xuzhou City, Jiangsu Province, China, with a purity of over 99%.
The elastic modulus and internal structure of MKP were tested by nanoindentation and X-CT techniques, respectively. Before the nanoindentation test, the specimens needed to be pretreated: the mass ratio of MKP particles to epoxy resin is 1/10, and these two materials were stirred for 5 min, then the epoxy resin hardened after 6 h as the specimens for the nanoindentation test. Since the nanoindentation test requires high flatness on the surface of specimens, it is necessary to perform smooth pretreatment on the test piece (Zhao et al., 2005;Zheng et al., 2008;Han et al., 2012). The surfaces of the specimens were polished in the order of 400-4,000 mesh sandpaper, canvas, and silk. The surface roughness of the polished specimens was tested by an atomic force microscope to ensure that the roughness was <100 nm. After smoothing the MKP, the testing region was marked and then the structure distribution characteristics of the region was obtained by X-CT.

X-CT Test
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 Xray 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.

Nanoindentation Test
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).

Random Aggregate Placement Method
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 . 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.

THE RESULTS OF EXPERIMENTS AND RANDOM AGGREGATE PLACEMENT
The Results of X-CT A total of 15 indentation area on the MKP specimens was scanned by X-CT. Three hundred and twenty two-dimensional CT slices with the resolution of 1,000 × 1,000 pixels were obtained by each X-CT scan and the spatial resolution was 0.5 µm voxel. Then the images were processed by the AVIZO software to gain the structural distribution characteristics and volume fraction of each phase in each indentation area. Subsequently, the tetrahedral mesh generated by AVIZO software was imported into ABAQUS software as the X-CT nanoindentation model of the MKP.
The X-CT images of the MKP was imported into the AVIZO software. Based on the difference of gray value between the MKP and pores, the slices were divided into two phases. Determining the threshold value of the phases was the premise to correctly distinguish MKP from pore. Therefore, a MKP particle with the size of about 3 mm was tested by mercury intrusion porosimeter (MIP) and X-CT, respectively. The porosity of the MKP particle was 28.43% from the MIP. Then the threshold of the X-CT slices of the particle was debugged. It was found that when the parts with the gray range of 0-149 and 149-207 were adopted to divide the pores and the MKP, the voxel number of the MKP and the pores were 5.33E8 and 2.12E8. So, that is to say, the porosity of the MKP particle is 2.12E8 ÷ (2.12E8 + 5.33E8) = 28.46%, which was almost consistent with the test result. Therefore, the gray threshold of distinguishing the MKP from pore was 149. For example, the gray value 149 was used as the threshold to divide the MKP and pores, when the X-CT data of the P 14 indentation point region was processed. The 160th slice of the P 14 indentation point is shown in Figure 2A. The dark gray area was the pores and the light gray area was the MKP. Figure 2B displays the result of the threshold segmentation of the 160th slice, where the red part indicates pores and the blue part represents the MKP. The distribution characteristics of the pores in the P 14 indentation point are shown in Figure 2C and the green curve in Figure 2D denotes the statistical pore size distribution. The yellow curve in Figure 2D displays the pore distribution function curve, which will be discussed in detail in section The Effect of Porosity on the Damage Factor. Additionally, the voxel number of the two phases was 2.9412E8 and 2.528E7, respectively. In other words, the volume fractions of the MKP and pores after the threshold segmentation were 2.9412E8 ÷ (2.9412E8 + 2.528E7) = 92.1% and 2.528E7 ÷ (2.9412E8 + 2.528E7) = 7.9%, respectively. Similarly, the porosity of each indentation point region could be obtained. There were 15 indentation points numbered P 11 -P 35 , as shown in Table 1.
Finally, the tetrahedral meshes were generated based on the X-CT model, as shown in Figure 2E, in which the gray portion indicates the MKP and the red portion represents the pores. The tetrahedral mesh model was imported into the finite element software ABAQUS and then the pore set was deleted. The remaining MKP set was used as the compressive matrix of the nanoindentation model. A triangular pyramid model with a height of 30 µm was established by CAD and imported into the finite element software ABAQUS as the indenter of the nanoindentation model. The X-CT nanoindentation model of the MKP was obtained by assembling the MKP matrix with the triangular pyramid indenter, as shown in Figure 2F, where the green part represents the MKP matrix and the yellow part denotes the triangular pyramid indenter. When the indenter was assembled with the MKP substrate, the upper surface of the triangular pyramid indenter was parallel to the upper surface of the MKP substrate, in which the central axes of the indenter and the MKP substrate were coincident. Additional details on the X-CT nanoindentation model were discussed in section Simulation Based on X-CT Results.

The Results of Nanoindentation
In this study, 15 marked areas were selected on the MKP specimens for nanoindentation test and the indentation points were numbered P 11 -P 35 with the elastic modulus and the peak loads shown in Table 2. There are two data at each indentation point, the elastic modulus and the peak loads, respectively. For example, the data "22.52, 1.91" for the indentation point P 11 indicates that the elastic modulus and the peak load at the indentation point P 11 are 22.52 GPa and 1.91 N, respectively. It can be seen that the elastic modulus of all the indentation points are in the range of 17.66-26.14 GPa. The average elastic modulus of the indentation points is 21.91 GPa. Besides, the peak loads fluctuate in the range of 1.57-2.23 N and the average peak load is 1.88 N. The fluctuations of the elastic modulus and the peak loads are mainly caused by the difference in volume fraction and distribution characteristics of the matrix and pores in each indentation area.

The Results of Random Aggregate Placement
First, a matrix part of 500 × 500 × 160 µm was built in ABAQUS and was meshed with the grid unit size of 1.5 µm and the element type of linear hexahedral elements C3D8R. Then, based on the random aggregate placement method (RAP), the RAP nanoindentation model of the P 14 indentation point was established. The pores were randomly placed in the matrix part and the porosity was set to 7.9%. Thus, the RAP nanoindentation model corresponding to the X-CT nanoindentation model was obtained. As shown in Figure 3A, the blue part indicates the MKP set and the white part represents the pores set. Subsequently, the pore set was deleted and the remaining MKP matrix was assembled with the triangular pyramid indenter. Hence, the RAP nanoindentation model of the P 14 indentation point was displayed in Figure 3B, where the blue part denotes the MKP matrix and the yellow part indicates the triangular pyramid indenter. When the indenter was assembled with the MKP substrate, the upper surface of the triangular pyramid indenter was parallel to the upper surface of the MKP substrate, in which the central axes of the indenter and the MKP substrate were coincident. Additional details on the model were discussed in section Simulation Based on the Random Aggregate Placement Method (RAP). Similarly, the RAP nanoindentation models of other indentation points were obtained.

SIMULATION BASED ON X-CT RESULTS
In section Simulation Based on X-CT Results, the nanoindentation test of the MKP was simulated based on the X-CT nanoindentation model. Firstly, the input parameters and boundary conditions of the X-CT model were determined. Secondly, considering the effect of pores on damage, the modified MKP constitutive relation was proposed and input into the X-CT nanoindentation model. Then, the equations of influence of porosity and pore distribution on the damage factor were proposed. Finally, the correctness of the modified damage factor equation was verified with the average relative error of 3.2% by comparing the numerical results with the experimental results.

The Input Parameters and Boundary Condition
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 3 , 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 −6 . The loading time-amplitude conformed to the following rules. When the time was between 0 s and 10 −6 , the indenter was in the loading stage with the displacement from 0 to 20 µm. When the time was in the range of 10 −6 to 1.1 × 10 −6 , the indenter was in the constant loading without movement. During the period of 1.1 × 10 −6 to 1.4 × 10 −6 , 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 −6 , A = t × 10 6 (2) when 10 −6 < t < 1.1 × 10 −6 , when 1.1 × 10 −6 < t < 1.4 × 10 −6 , A = −2.1 t − 10 −6 2 × 10 12 + 1 where A is the amplitude, t represents time.

The Constitutive Relation of the MKP
The MKP constitutive relation was plastic damage model and determined in previous studies (Li et al., 2018a) as follows: 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.

Damage Factor of the MKP and Simulation Results
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, ε in 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 14 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 loaddisplacement 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 14 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 −6 . 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 −6 , 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 −6 . 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 −6 , 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 11 indicates that the optimal influence coefficient value and the optimal simulated peak load at the indentation point P 11 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 Effect of Porosity and Pore Distribution on the Damage Factor The influence of pores on the damage factor came from two aspects: on the one hand, the larger the porosity was, the greater the damage factor was. On the other hand, under the same porosity, different pore distributions led to differences in damage factors. Therefore, the total influence coefficient χ of pores on the damage factor was a comprehensive reflection of the influence of porosity and pore distribution. Then, the influence coefficient of porosity and pore distribution on the damage factor were defined as β and λ, respectively. Thus, the relationship of the three coefficients was assumed to be: Considering the influence of porosity and pore size on the damage factor, the modified MKP damage factor is expressed as follows: After obtaining the total pore influence coefficient of 12 indentation points, the effects of porosity and pore distribution on the damage factors were further discussed.

The Effect of Pore Size Distribution on the Damage Factor
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, The value of λ i is the influence coefficient of the pore median of the i-th interval. c indicates the total porosity, namely, c = n i = 1 c i . The smaller the interval is, the closer the Equation (26) is based on the Equation (25).

The Effect of Porosity on the Damage Factor
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 14 in Figure 2D, the pore distribution function of the P 14 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 λ 1 , λ 2 , λ 3 , λ 4 , λ 5 , λ 6 , λ 7 , λ 8 , and λ 9 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 14 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 2 represents the correlation coefficient of the fitted equation. The R 2 -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:

Verification of the Modified Damage Factor Equation
The damage factors calculated by Equations (29) and (30) were input into the X-CT nanoindentation model and the three indentation points of P 15 , P 25 , and P 35 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.

Verification of the RAP Nanoindentation Model
The nanoindentation test of P 14 was simulated based on the random aggregate placement method (RAP). The relevant parameter settings for the RAP nanoindentation model were identical to those of the X-CT nanoindentation model in section The X-CT Nanoindentation Model. The MKP material properties were the same as those of the X-CT model: the MKP density was 1.864 g/cm 3 , 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. Besides, the same boundary conditions were set: the bottom of the matrix was fixed and the indenter could only move in the Z direction. Additionally, the modified damage factor considering the pores were input into the RAP model. Finally, the simulated peak load of P 14 by the RAP model was 2.17 N with the relative error of 6.4%. The simulated nephogram of the RAP model is displayed in Figure 8.
It can be seen that the RAP nanoindentation model also exhibited stress concentration, severe mesh deformity and crowding effect. In addition, the pores in the lower region of the indenter were seriously deformed and the pores affected the transmission of stress and damage. That is to say, the simulation results of the RAP nanoindentation model were consistent with those of the X-CT model. According to the above steps, the simulation results of P 15 , P 25 , P 35 were obtained by the RAP model, listed in Table 5. It was calculated that the average relative error of the peak load of was 5.9%, which verified the validity of the RAP nanoindentation model.

Prediction of Nanoindentation Test
The pore distribution certainly affected the results of the nanoindentation test, but it was difficult to find the indentation points with the same porosity in the X-CT model to study the effect of the pore distribution on the results alone. Therefore, the RAP model can be used to study the effect of pore distribution by controlling both porosity and pore distribution characteristics. The effect of the pore distribution on the total influence coefficient and peak load was predicted with the same porosity of 9% and the influence function k (D) = 2.45D 0.18 . Firstly, the pore distribution was in the range of 0.0005-0.03 mm. Then, the pore distribution dominated by small pores, medium pores and large pores were, respectively, set, listed in Table 6. Subsequently, according to Equations (29) and (30), the total influence coefficient and peak load prediction results are exhibited in Table 6 and Figure 9 displays the simulated displacement-load curve. It can be seen that at the same porosity, with the increase of the distribution of the large pore size, the total influence coefficient of MKP increased, the peak load decreased and the load-displacement curve as a whole was lower. It indicated that the effect of pore size on peak load should not be neglected.

CONCLUSIONS
In this study, the micro-mechanical properties and structural distribution characteristics of the MKP were investigated by nanoindentation and X-CT techniques. The MKP nanoindentation model was established based on X-CT and random aggregate placement method. The constitutive relationship of the MKP was modified considering the effects of porosity and pore distribution. The influence of pore distribution on the results was predicted by RAP model. The conclusions are as follows: 1. The elastic modulus of the indentation points of the MKP are in the range of 17.66-26.14 GPa with the average elastic modulus of 21.91 GPa while the peak loads of the indentation points fluctuate in the range of 1.57-2.23 N and the average peak load is 1.88 N. 2. The nanoindentation plastic damage model based on X-CT method can effectively describe the mechanical and deformation characteristics of the MKP.
3. The average relative error of simulated peak load is only 3.2% based on the X-CT model, which indicates that the hypothesis of the relationship between the influence coefficients β and λ, the MKP pores conforming to the exponential distribution and the influence function of pore size distribution are correct. That is, from the perspective of the damage factor, the modified mkp constitutive relation considering porosity and pore distribution is effective. 4. The mechanical and deformation characteristics of the MKP can also be described effectively by the RAP model, which verifies the valid of the modified constitutive relation of the MKP again.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation, to any qualified researcher.

AUTHOR CONTRIBUTIONS
YL provided the research direction and funding equipment support. GZ was responsible for testing, simulating, and writing papers. ZW provided the detailed research ideas and suggestions for the research.