Numerical Simulation of Mesodamage Behavior of Concrete Based on Material Point Method

Concrete consists of coarse aggregates, mortar matrix, and interfacial transition zone (ITZ) between them at the mesoscale. Considering these three phases, many numerical tests have been conducted to study the mesodamage behavior of concrete, in which a variety of numerical methods have also been adopted. *ese methods are mainly based on the finite element method (FEM); however, some other methods have been proven to be helpful as well. For example, the material point method (MPM) has the advantage of building a numerical model based on pixel or voxel of the image and is capable of solving large deformation problems. In view of this, MPM is introduced in this paper. Firstly, a method for establishing the numerical specimen is put forward, considering the original sample of its mesoscopic geometric character. *en, a stochastic damage constitutive model considering the heterogeneity of the concrete is proposed. Next, the numerical model and the constitutive model are incorporated into an MPM code to conduct numerical tests. *e uniaxial tension and compression tests of a random-aggregate model and a double-aggregate specimen under uniaxial tension are then simulated numerically to validate the proposed method. Results show that the proposed method can well capture the main macroscopic mechanical behavior of concrete and the mesoscopic damage initiation and propagation. It is also found that MPM can save the time of model establishing and improve calculation efficiency. *e influences of different parameters of the proposed constitutive model are also clarified through a parametric study. *e proposed method can provide a useful tool for concrete numerical testing and for studying the mechanical behavior of concrete at mesoscale.


Introduction
Concrete is the most widely used construction material all over the world, and the study of its mechanical behaviors still remains an important issue in civil engineering and material science [1][2][3][4][5]. From the viewpoint at mesoscale, concrete is a typical composite material, which consists of three main phases: coarse aggregates, mortar matrix, and interfacial transition zone (ITZ) between them [5][6][7]. erefore, concrete is highly heterogeneous at the mesoscopic level, and thus the macroscale mechanical behaviors of concrete, which draw the attention of engineering practice, are greatly determined by its mesoscale structure.
Since the mesoscale study of concrete mechanical behaviors can obtain a more fundamental understanding of concrete damage and failure mechanisms and provide a deep insight into macroscopic mechanical behaviors, it has become a hot research area of concrete material and achieved fruitful results. Up to date, the finite element method (FEM), which helps many researchers to simulate the effects of multiphase distribution on the mechanical response and damage properties of concrete [5][6][7][8][9][10][11][12][13][14][15][16][17][18], is still being a dominant method to model different behaviors of concrete at macro-and mesoscales. However, with the development of computing technology and numerical methods, more and more approaches are being employed to explain the damage mechanism of concrete structure [19][20][21][22][23][24][25]. For example, on the basis of considering the energy dissipation during crack opening, Sageresan and Drathi [19] simulated the crack propagation in concrete with the help of meshless method. Shen et al. [21] has applied the peridynamic model to study the crack propagation path of the mesoscopic concrete, and the numerical results fit well with the experimental outcomes. Farahani et al. [22] used a radial point interpolation meshless method (RPIM) to attain the nonlinear damage solution, and the numerical results are in agreement with the experimental results. But one of the disadvantages that exist among all the numerical simulations is the time consumption of a numerical specimen. erefore, a new numerical method, the material point method, which may help to improve this shortcoming, is introduced in this study.
Material point method (MPM) is a numerical method evolved from the Particle-in-Cell method, which was first proposed by Sulsky and Chen in 1994 [26][27][28]. MPM uses Lagrange particles to discretize the material region and the Euler background grid to calculate the spatial derivative and solve field equations. With the help of the Lagrangian description, MPM can easily form an establishment of a numerical model from the material body by discretizing the continuum. MPM has already been used effectively in modeling the microstructural evolution of the polymethacrylimide (PMI) closed-cell plastic foam [29], using the microstructures reconstructed from X-ray tomography. In terms of the mesoscopic problem, due to the beneficial characteristic properties of MPM, it can easily build the relationship between the material point (Lagrange particles) and the pixel or voxel of the image or tomography [30][31][32] and further establish the numerical specimen, which helps to improve the efficiency of numerical models construction. Since MPM has a beneficial effect on solving large deformation problems, so far, MPM has helped researchers to study bulk metallic glass and nanofoams from the mesoscopic level [31,32]. Yet, it is seldom adopted for concrete numerical simulation at mesoscale. e main purpose of this study is to try to adopt MPM to study concrete structures' mesodamage behavior at mesoscale. erefore, a numerical method based on MPM is developed. e uniaxial tension and compression tests of a random-aggregate model and double-aggregate specimen under uniaxial tension are then simulated numerically to validate the applicability of MPM. In addition, a parametric study is carried out to further study the performance of the proposed method.

Outline of MPM
In MPM, deformations of a body are simulated by the material points, which carry all the information of the continuum, including density, velocity, acceleration, stress and strain, and other parameters. Despite being a meshless method, MPM is often regarded as an extension of FEM.
is is because the underlying Eulerian grid in MPM is used to solve the equilibrium equations for each-time increment just like a computational finite element mesh used in FEM. e main difference between MPM and FEM is that all the information is mapped from grid nodes to material points at the end of a time step while at the following step the mapping changes from material points to grid nodes, thus preventing MPM from mesh distortion that happened in traditional FEM. MPM is specifically characterized as two algorithms, the implicit and explicit algorithms. In this study, the explicit algorithm is used, which is more popular than the implicit algorithm in the MPM society. Please note that quasi-static problems are considered here; however, the explicit algorithm, as a true dynamic procedure, has proven valuable in solving quasi-static problems as well.
Generally, the explicit computation scheme of MPM is proceeding as follows: first, state variables are interpolated to the grid nodes at a certain time step, and the governing equations of motion are thus solved at the grid nodes. en, the state parameters, including density, velocity, and stress, are updated by mapping grid nodes to material points at the current time step. Finally, the used Eulerian grid will be replaced by a new fixed grid while keeping the material points still at the next time step. Since each material point contains a fixed amount of mass for all time, mass conservation is automatically satisfied. Assume a one-phase dynamic problem governed by the discretized momentum equation: where M is the mass matrix, a is the acceleration, F ext is the external force, and F int is the internal force. Equation (1) is derived from the weak form of the momentum balance equation by applying the Galerkin procedure [26]. At the start of each step, information mapped from material points to grid nodes is initialized by the interpolation functions. en the Gauss quadrature used in FEM is adopted in MPM. Information regarding the state of a certain material volume Ω p is stored at the material points; integrals can therefore be calculated based on material points and the internal force term in equation (1) can be written as where B q contains the derivatives of the shape functions computed at the material point location, σ p is the total stress of a single material point, and n p is the number of material points contained in the considered element. Nodal accelerations a are then solved by the governing equations of motion in equation (1). anks to the shape functions, the incremental velocities can be obtained. Incremental strains can be calculated from nodal velocities of material points, and the stress, displacements, locations, and other variables can thus be updated in material points.
At the end of a time step, the Eulerian grid can be abandoned as all the updated information has been mapped from grid nodes back to material points. A detailed solution description can be found in [27][28][29][30]. e numerical program applied in this study is updated according to the outline of MPM and compiled by MATLAB.

Generation of Geometric Models of Concrete at Mesoscale.
Before proceeding to validate the MPM program by conducting uniaxial tension and compression tests of a random-aggregate specimen, a numerical specimen consisting of random coarse aggregates, mortar, and ITZ needs to be generated. In order to obtain a numerical 2 Advances in Civil Engineering specimen, a geometric model is needed first.
us, a method for generating a geometric model is introduced in this study. e dimensions of the specimen are 150 mm by 150 mm. For ordinary concrete, coarse aggregates generally occupy 40%∼50% volume of concrete and form the material skeleton, which plays an important role in mechanical properties in concrete. In this study, we choose 45% as the target volume fraction of coarse aggregates. e coarse aggregates may be gravels or crushed stones. Gravels have a round shape and are generally modeled by circular or elliptical particles, while crushed stones have an angular or irregular shape and are generally modeled by polygonal particles [5,6]. is study chose gravel aggregates as a representative without losing generality, and the elliptical particles with an aspect ratio between 0.5 and 1 are adopted to model coarse aggregates. A similar procedure with [6] is utilized to generate and arrange coarse aggregates in a repeated manner until the target volume fraction of coarse aggregates is achieved.
e concrete mesoscopic geometric model is established by programming by means of MATLAB. e aggregates are placed in the following procedure: (1) Program initialization: the size of a concrete geometric model is initiated as 150 by 150 mm, which means the axis length is 150 by 150. Meanwhile, the minimum distance between each generated aggregate and the existing aggregate is 0.5 mm in addition to the same limitation with it and the boundary of the model as well.
(2) Aggregate generation: according to the size range of different stones, the position coordinates of each ellipse, dropping angle, and the radius of the major and minor axis of aggregate are produced randomly and then placed into the material point distribution model. (3) Placing determination: firstly, generated aggregate is checked whether it is inside the sample or not. By determining the distance between the elliptical center coordinates and the four endpoints of the long and short axes, the maximum and minimum x and y coordinates are found, and the distance between each coordinate and the boundary of the model can be compared. If the maximum x coordinate is smaller than the difference between the maximum boundary and the minimum distance of the model, which is 149.5 mm, and the minimum x is also greater than the sum of the minimum boundary and the minimum distance, namely, 0.5 mm, then the coarse aggregate is regarded as effectively placing. en move on to the next step. Otherwise, the generated aggregate needs to be given up and return to Step (2).
(4) Judging intersecting relation: the intersecting relationships between the current ellipse and existing ellipse are judged by the distance between the center of current and each existing ellipse. If the distance is greater than the sum of two respective long axis radii of the current ellipse and existing ellipse, it can be assumed that there is no intersection between them, and then go to Step (5). However, if the distance between these two centers is determined to be less than the sum of their short axis radius, then they definitely have a cross section, meaning this aggregate should be abandoned and move back to Step (2). (5) Judging tangent relation: firstly, the coordinates of current ellipse and each existing ellipse are rotated and translated into the matrix form, and then the eigenvalues and determinants of the two matrices are solved, respectively, to define whether there are two different positive roots; if they exist, meaning they are not tangent, then enter Step (6); otherwise, this random placement is aborted and repeat Step (2). (6) When an aggregate is successfully placed in the model, it is numbered and the current volume fraction area of aggregate and the total volume are calculated. (7) Repeat Steps (2)-(6) until the total volume of aggregate reaches the target volume fraction and end the generation process. e generation process of one of the geometric image models is shown in Figure 1.

Generation of Concrete Numerical Specimen.
Based on the geometric model of random aggregates, a method for the generation procedure of the numerical specimen is proposed. Firstly, a geometric image is obtained and the coarse aggregates and mortar matrix are separated clearly. en, the ITZ is ignored at first, and thus part of the material point mesh is occupied by mortar, and the remaining zone is filled with coarse aggregates. Finally, the ITZ is inserted into the material point mesh. e specific steps are shown as follows: (1) Material point mesh initiation: the size of the numerical model is initialized to 150 by 150 units, which contains 4 material points in each unit, so the final mesh initialization is set to 150 by 150 by 4 material points. After that, each material point is assigned a coordinate. In addition, it has been concluded that the volume ratio, when ranging from 0.1% to 10%, between material points and numerical model, can have a better effect on the computation accuracy and efficiency [22]. In order to obtain more accurate results, we choose 90,000 material points to conduct numerical experiments without affecting the computational efficiency. (2) Aggregates coordinate determination: according to the number, coordinate, placing angle, and radius of the long and short axis of each ellipse, all the internal material point coordinates of each aggregate are defined one by one. By determining the elliptic focus of current aggregate, the distance between all material points in the coordinates and current two focus points is determined one after the other. If the distance is no more than the diameter of the long axis  (2) and (3), followed by the geometric model being discretized by material point mesh. According to the geometric image, different types of pixels build a relationship with material points, which indicates that the concrete sample has been discretized by different material points. Particles carrying different information are discretely mapped to the initialized material point mesh so as to characterize the three-phase mesoscopic structure of the material point distribution model, as shown in Figure 2. (5) Material point distribution model establishment: three different kinds of material points, carrying mortar, aggregate, and ITZ geometric information, respectively, are assigned the corresponding parameter information, including elastic mechanics parameters, plastic mechanics parameters, and survival probability. e establishment of a material point distribution model of different mesoscopic structure is then completed. Take Figure 1 as an example; the final material point model for the numerical specimen consists of 92628 material points. Among all these material points, there are 49386, 40614, and 2628 particles for aggregate, mortar, and ITZ, respectively.
By combining the geometric image model with the material point method, a numerical specimen with a size of 150 by 150 of concrete mesostructure with random elliptical aggregates, volume fraction of which is 45%, is established. e average time of the whole modeling process is approximately 10 minutes, which is 66.7% quicker than that of using traditional FEM [5].

Stochastic Constitutive Model Used in MPM
e stress renewal of a material point is obtained by calculating the updated stress increment from the updated strain, the damage variable, and the elastic modulus of the material point and by adding the stress retained in the previous step: where D is the damage variable corresponding to a specific material point of mortar or ITZ at a certain time step n, E p is Young's modulus of the material point, and ε p is the updated strain.
Considering the primary characteristic of concrete is its stochastic heterogeneity at mesoscale, the classical Weibull's failure theory is resorted to [33]. According to Weibull's failure theory, the material strength complies with the Weibull distribution. e probability density function of a Weibull random variable is where m > 0 is the shape parameter and λ > 0 is the scale parameter of the distribution. e cumulative distribution function for the Weibull distribution is e probability density function and cumulative distribution function for the Weibull distribution are shown in Figure 3. erefore, for concrete of a representative material point, the survival probability P S (V ref ) under a type of maximum stress σ can be calculated as  In the present study, two failure modes, the tensile failure mode and shear failure mode, are used to characterize the material points in mortar and ITZ. For the material points representing aggregate, their failure modes need not be considered as they are considered linear elastic. Weibull's theory [33] based on the maximum principal stress is adopted here. According to Weibull's theory, the survival probability from tensile damage P st of material samples of a material point under a maximum principal stress σ 1 can be calculated as

Advances in Civil Engineering
in which f t is the mean uniaxial tensile strength at which the fraction of 1/e ≈ 36.8% of a material point survives. m, as the shape parameter, is reflecting the homogeneity of the material and how rapidly the survival probability reduces as σ 1 approaches f t . e survival probability from shear damage P ss is given by In order to implement the constitutive model in the material point program, at the beginning of the simulation test, a random number r between 0 and 1 will be assigned to all the material points of mortar and ITZ. At each time step,  Advances in Civil Engineering the survival probability P s of each material point is calculated according to equations (7) and (8) and compared with the random value r assigned in the beginning. If the survival probability P s of the material point is bigger than r, the material point is considered to be undamaged, and the damage variable D is taken as 0. If P s is less than r, the material point is taken to be damaged and invalid, and damage variable D is set to (1 − ξ). Since the information of all material points in the previous step continue to be used in the next step, it is necessary to adjust ξ into a minor number so as to maintain numerical stability. In this paper, ξ is taken to be 1 × 10 − 6 . e failed material point is referred to as a "ghost material point" and no longer contributes to the volume integration. Using k as an index for the ghost material point, equation (2) can be written as in which n k is the total number of ghost material points. Because of the piecewise property of the shape function, F int ghost reaches its maximum at node k, which is the nearest node on the computational mesh to material point k. In the node-duplicating scheme of the traditional FEM [5,34], the node k has to be duplicated and separated in order to generate a crack. e degrees of freedom have to be renumbered and a new element connection has to be formed, which requires major computer programming efforts. In MPM, on the other hand, it is the discrete material points that are separated instead of the background mesh. e separation of material upon failure takes place automatically and does not require any further treatment in computer programming. e "ghost points" do not carry any stress but still contribute to the mass matrix as it is sensible to consider their effect of inertia on the damage behavior.
Please note that, in most of the existing studies [1,3,[5][6][7]34], the heterogeneity at mesoscale is considered by generating numerical specimens with randomly distributed aggregates and ITZs. In this study, with the incorporation of Weibull distribution, the heterogeneity at a smaller scale is considered; that is, the mortar itself is also heterogeneous and so is the ITZ. is makes the numerical specimen more alike to a real one.

Uniaxial Tension and Compression Tests of a Random-Aggregate Specimen.
In order to validate the proposed method, both the uniaxial tension test and uniaxial compression test are conducted numerically. e loading schemes of the numerical simulation are depicted in Figure 4. e material parameters are listed in Table 1. As for Young's modulus E and Poisson's ratio ], the parameters are generic and suitable for usual concrete [5]. As for the parameter m and the strength parameters f t and c, a calibration process is carried out to determine their values, for there has been no reference available yet. As for the strength parameter φ, we take the reference values of usual concrete to reduce the calibration effort, because the variation of strength parameter φ is relatively smaller when compared with that of strength parameter c. e strength of ITZ is taken to be approximately 60% weaker than mortar according to FEM numerical testings [5]. e macroscale results for uniaxial tension and uniaxial compression tests obtained by numerical simulation are compared with several laboratory experimental results. Relations between normalized stress σ/σ p and normalized strain ε/ε p are shown in Figure 5, in which σ p and ε p are peak stress and peak strain, respectively. It can be seen that although a certain variation in experimental results is exhibited, the stress-strain curves for uniaxial tension and compression tests can be well captured by the MPM program.
e mesoscale damage evolution processes at mesoscale for uniaxial tension test and uniaxial compression test are shown in Figures 6 and 7 respectively. In order to demonstrate the damage evolution process, the damage status of all material points is depicted corresponding to 25%, 50%, 75%, 100%, 125%, and 150% of peak strain. It is clear that the damage initiation and propagation process can be well captured by the MPM program. Meanwhile, localized damage can also be well reproduced. Note that, even for the uniaxial compression case, the main damage type at mesoscale is the tensile type, which coincides with the findings reported by Zhu and Tang [35].

Uniaxial Tension Test of a Double-Aggregate Specimen.
To further validate the proposed method, the laboratory uniaxial tension test of a double-aggregate specimen conducted by Corr et al. is numerically simulated. In the laboratory test [36], a total of 7 two-aggregate specimens were tested by affixing steel plates to the specimen for mounting in a direct tension grip. e studied specimen dimensions are 75 mm by 75 mm. e peak stresses for all two-aggregate specimens range from 1.40 MPa ∼ 2.09 MPa, while the peak strains range from 45.92 µm/m to 103.21 µm/m.
A numerical specimen corresponding to the physical specimens is generated, as shown in Figure 8. e generation procedure of the numerical specimen is similar to the random elliptical aggregate concrete specimen. e loading scheme of the numerical simulation is depicted in Figure 9, along with the stress-strain response obtained by the numerical test. e values of peak stress and strain coincide well with those obtained by the laboratory tests. e damage evolution processes at the mesoscale obtained by the numerical test are shown in Figure 10. An obvious localized damage zone can be found near the top of the two aggregates, which agrees well with the laboratory test (see Figures 11(a) and 11(c)). Note that Xu and Chen also simulated this laboratory test [5], but they obtained two localized damage zones. One was near the top of the two aggregates, while the other was near the bottom of them (see Figure 11(b)). ey attribute this discrepancy to the fact that the constitutive model they employed lacks the capability of describing strain localization in the real specimen, and the bond strength at the tops of aggregate may be lower than that at the bottoms. On the contrary, this discrepancy can just be overcome thanks to the stochastic properties of the constitutive model used in this study; the bond strength at the top of aggregates is indeed lower than that at the bottom of them in some stochastic cases.   Advances in Civil Engineering

Parametric Study
As mentioned in the previous section, there has been no reference available yet for different parameters when conducting the simulation of uniaxial tension and compression tests for the random-aggregate specimen. erefore, the influences of the homogeneity parameter m, the strength parameter for tension f t , and the strength parameter for shear c are analyzed, respectively, so as to study the influence of different parameters in the MPM program. When a parameter is studied, the other parameters remain the same as those in Table 1.
e first parameter to be studied is the homogeneity parameter m, which is a probability parameter. According to Weibull's theory, the larger the parameter m is, the more homogeneous the material is. e parameter m is taken to be 1.0, 1.5, 2.0, 3.0, and 6.0, respectively. e numerical results for both tension and compression tests are shown in Figure 12. It can be seen that, with the increase of m, the peak strength of the specimen becomes larger, while the specimen becomes more brittle. It is due to this fact, with the increase of the homogeneity of the specimen, that the mesoscale damage tends to appear at the same time.  Advances in Civil Engineering and 14. It can be said that both strength parameters have a direct influence on peak strength. e parameter f t has a more evident effect on tension peak strength than on the compression peak strength, while the parameter c mainly affects the compression peak strength and has a slight influence on the tensile peak strength.
It should be noted from Figure 13(b) that when the peak compressive strength reaches 25 MPa, the effects of the increase of parameter f t begin to wear off. Although the peak strength still increases, the increase is not obvious. is is because the strength parameter c remains unchanged; the interior shear strength of the concrete structure is relatively fixed. Due to the fact that the mesodamage type is the tensile failure, when the compressive strength reaches a certain peak, the internal structure has already met the shear failure, which makes the peak stress also stabilize at a certain point. Similarly, as can be seen from Figure 14(a), parameter c also has a small impact on tensile peak strength but rather mainly affects compressive peak strength. It is also similar to the rule of influence of parameter f t ; that is, the tensile failure is met earlier than the shear failure.
is explains why the parameter c has a more direct effect on compression tests of concrete mesostructure but less effect on tension tests.

Conclusions
Concrete is a typical multiscale heterogeneous material. A novel numerical method based on MPM is proposed to study the mesodamage behavior. According to the original sample of the geometric image of the concrete, a method for establishing the numerical specimen for MPM is presented.
In order to consider the heterogeneity of concrete at mesoscale, a stochastic damage constitutive model is also put forward. en the numerical specimen and the proposed constitutive method are incorporated into an MPM code complied by MATLAB, and the program is validated by simulating the uniaxial tension and compression tests of a random-aggregate specimen and by simulating the uniaxial tension test of a double-aggregate specimen. It is found that the proposed method can well capture the main macroscale mechanical behavior of concrete and reproduce the damage initiation and propagation at mesoscale. Meanwhile, this approach can help to save the time of the model establishing and improving calculation efficiency. e parametric study of the model has revealed some properties of the proposed constitutive model. With the increase of homogeneity parameter m, the peak strength of the specimen becomes larger, while the specimen becomes more brittle.
e strength parameters f t and c have a clear influence on tensile strength and compressive strength, respectively. e proposed method is testified to be a helpful tool for concrete numerical testing and for studying the mechanical behavior of concrete at mesoscale.

Data Availability
e data used to support the findings of this study are included within the article.

12
Advances in Civil Engineering