A Research on the Macroscopic and Mesoscopic Parameters of Concrete Based on an Experimental Design Method

Concrete is a composite material that has complex mechanical properties. The mechanical properties of each of its components are different at the mesoscopic scale. Studying the relationship between the macroscopic and mesoscopic parameters of concrete can help better understand its mechanical properties at these levels. When using the discrete element method to model the macro-mesoscopic parameters of concrete, their calibration is the first challenge. This paper proposes a numerical model of concrete using the particle discrete element software particle flow code (PFC). The mesoscopic parameters required by the model need to be set within a certain range for an orthogonal experimental design. We used the proposed model to perform numerical simulations as well as response surface design and analysis. This involved fitting a set of mapping relationships between the macro–micro parameters of concrete. An optimization model was established in the MATLAB environment. The program used to calibrate the mesoscopic parameters of concrete was written using the genetic algorithm, and its macro-micro parameters were inverted. The following three conclusions can be drawn from the orthogonal test: First, the tensile strength and shear strength of the parallel bond between the particles of mortar had a significant influence on the peak compressive strength of concrete, whereas the influence of the other parameters was not significant. Second, the elastic modulus of the parallel bonding between particles of mortar, their stiffness ratio and friction coefficient, and the elastic modulus and stiffness ratio of contact bonding in the interfacial transition zone had a significant influence on the elastic modulus, whereas the influence of the other parameters was not significant. Third, the elastic modulus, stiffness ratio, and friction coefficient of the particles of mortar as well as the ratio of the contact adhesive stiffness in their interfacial transition zone had a significant influence on Poisson’s ratio, whereas the influence of the other parameters was not significant. The fitting effect of the response surface design was good.


Introduction
No unified method is available to describe the macroscopic parameters of particles of geotechnical materials by using their microscopic parameters. Although laboratory experiments are a feasible approach, they have such problems as low efficiency, high cost, and a large dispersion in the results. Cundall established the discrete element method (DEM) to solve this problem [1]. The DEM is a numerical method to examine the mechanics of discontinuous media based on Newton's second law of motion. It represents the given geotechnical material as a rigid particle model in which its mesoscopic and macroscopic parameters can be related. It is directly related to the geometrical characteristics of the assignment particles and particle contact between mesoscopic mechanics parameters of the rather than that at a point. This can better reflect the capability of concrete mortar to resist torque as well as the interfacial transition zone between particles. Therefore, in the concrete model in this paper, different contact models are set between aggregate and aggregate, and between mortar and mortar, so as to more truly reflect the actual performance of concrete. The parallel bonding model was used to simulate the mechanical behaviors of mortar and the aggregate, and the linear contact bonding model was used to simulate the contact-related behavior of coarse aggregate.
( Figure 1b) form the foundation of the PFC2D model. The stick mo tact force, and the parallel bond model in 2D/3D represents the line than that at a point. This can better reflect the capability of concrete as well as the interfacial transition zone between particles. The model in this paper, different contact models are set between ag and between mortar and mortar, so as to more truly reflect the actu crete. The parallel bonding model was used to simulate the mecha tar and the aggregate, and the linear contact bonding model wa contact-related behavior of coarse aggregate. The strength and deformation characteristics of concrete wer terized by using the three macroscopic parameters of peak strengt Poisson's ratio, and were measured using uniaxial compression t σu, elastic modulus E and Poisson's ratio ν were then used as test in influence of the test factors (mesoscopic parameters of concrete) o In the process of uniaxial compression in the particle discrete elem as the stress corresponding to the apex of the stress-strain curve, E the curve when the stress reaches half the peak stress for the first lute value of the ratio of lateral strain to axial strain at this time.
In simulations of the uniaxial compression test, the geomet Rmax/Rmin), load parameters ( , / , , ), parameters of t μ) and key binding parameters ( ̅ , , / , , ̅ ) affect the result ters of the particles (including particles of cement mortar (ball) a (clump)) and bonds are the focus of this study, the geometric and controlled in advance. The compression of members in engineeri because of which prismatic specimens can better reflect the comp crete than cubic specimens. A model size of 200 × 100 mm 2 was th set to 0.12, the radius ratio of the particles of cement mortar (Rmax density (ρball) was 2000 kg/m 3 , and aggregate density was 2650 kg/ age radius of particles of cement mortar ( ) to 0.5 mm. The bond normally set to one, and remains constant during the calculation; sidered. To reduce the number of unknown parameters, the stiff particles of cement mortar, their bonds, and the upper and lower The strength and deformation characteristics of concrete were preliminarily characterized by using the three macroscopic parameters of peak strength, elastic modulus, and Poisson's ratio, and were measured using uniaxial compression tests. The peak strength σ u , elastic modulus E and Poisson's ratio ν were then used as test indices (objects), and the influence of the test factors (mesoscopic parameters of concrete) on them was examined. In the process of uniaxial compression in the particle discrete element model, σ u is taken as the stress corresponding to the apex of the stress-strain curve, E is the tangent slope of the curve when the stress reaches half the peak stress for the first time, and ν is the absolute value of the ratio of lateral strain to axial strain at this time.
In simulations of the uniaxial compression test, the geometric parameters (L,W,n, R max /R min ), load parameters (E w c , k w n /k w s , ν w , µ w ), parameters of the particles (ρ, E c , k n /k s , µ) and key binding parameters (λ, E c , k n /k s , σ c , τ c ) affect the results. Because the parameters of the particles (including particles of cement mortar (ball) and aggregate particles (clump)) and bonds are the focus of this study, the geometric and load parameters were controlled in advance. The compression of members in engineering is mostly prismatic, because of which prismatic specimens can better reflect the compressive capacity of concrete than cubic specimens. A model size of 200 × 100 mm 2 was thus set. Porosity (n) was set to 0.12, the radius ratio of the particles of cement mortar (R max /R min ) was 2.5, particle density (ρ ball ) was 2000 kg/m 3 , and aggregate density was 2650 kg/m 3 . We preset the average radius of particles of cement mortar ( R ) to 0.5 mm. The bonding radius factor ( λ ) is normally set to one, and remains constant during the calculation; thus, it was not be considered. To reduce the number of unknown parameters, the stiffness parameters of the particles of cement mortar, their bonds, and the upper and lower loading plates were set to be consistent, that is, k n /k s = k n /k s = k w n /k w s , E c = E c = E w c . Parameters of the particles of aggregate were set to be consistent with those of contact stiffness, i.e., (k n /k s ) g = (k n /k s ) g and (E C ) g = (E c ) g given a wall loading rate ( v w ) of 0.02 m/s, the friction coefficient of the wall (µ w ) was zero.
Finally, the three macroscopic parameters σ u , E and ν were determined as test indicators. Parameters of the particles of cement mortar and their bonds (E c , k n /k s , σ c , τ c , and µ), the aggregate ((E c ) g , (k n /k s ) g , µ g ), and the interfacial transition zone and its bonds ((E C ) j , (k n /k s ) j , (σ c ) j , (τ c ) j , µ j ), a total of 13 mesoscopic parameters, were unknown. They were chosen as the experimental factors. y = E c , k n /k s , σ c , τ c , µ, (E c ) g , (k n /k s ) g , µ g , (E c ) j , (k n /k s ) j , (σ c ) j , (τ c ) j , µ j . The relationship is as follows: As concrete is a complex three-phase composite material with obvious anisotropy, the microscopic parameters of the aggregate, cement mortar, and the weak interfacial transition zone all need to be considered. The interfacial transition zone features both the aggregate particles and the cement paste, where the distribution of the cement particles is affected by the aggregate surface at the parts in contact. Due to the presence of many original cracks in the interior, most of the damage under external load starts from the interfacial transition zone, because of which it is the weakest part of concrete. The strength of the zone affects the overall mechanical properties and failure characteristics of concrete [9] as well as the fluctuation in its stress-strain curve. Therefore, it must be considered to enable the model to represent empirical scenarios involving concrete. The aggregate is the "skeleton" of concrete that plays an important load-bearing role. In the PFC2D software, the thickness of the interfacial transition zone cannot be expressed because the bond between particles is set through contact. In this paper, the influence of the thickness of the interfacial transition zone on the mechanical properties of concrete was ignored.
The parameters of the model are shown in Table 1.

Particle Formation
The literature has shown that when the number of particles in the sample is greater than 2000, the peak axial stress of the given sample is relatively stable, and the precision of the parameters of concrete obtained using a simulation is acceptably close to that obtained through laboratory experiments [10], The particles in the sample were filled using the radius expansion method. The number N of particles to be placed was estimated according to Equation (4). The number of particles was calculated to be 22,000, which is acceptable: where L is the height of the model (mm), W is its width (mm), n is its porosity (-), R is the average particle radius (mm), and is the rounded down result of the calculation. To restrict the magnitude of the overlap between the randomly generated particles, they were first produced at half the final particle size, and the porosity n 0 and average radius R 0 of the model were calculated. Then, the particle radius was amplified to obtain the desired particle size: The amplification coefficient m is given by:

Aggregate Generation
Past work has used the CT scanning technique to simulate the aggregate, and the scanned images are subjected to threshold processing using third-party software to obtain the boundary information to establish the model. However, this often reflects only one or more conditions while ignoring the random aggregate model. This paper uses a random algorithm for generating a convex polygonal aggregate.
In this algorithm, the circle is used as the base to generate the aggregate to ensure that the model is convex polygonal. The number of points on the circle determines the number of sides of a convex polygon. To ensure the randomness of the number of edges of the generated aggregate, a minimum value and a maximum value were set for the number of points on the circle. The range of the number of edges of the convex polygon was determined by calculating the difference between them, which was then multiplied by a random number. To round up the results, the resulting integer was set as the number of points on the circle, which was also the number of sides of the convex polygon: where N is the number of points on the boundary of the circle, N min is the minimum number of points generated. N max is the maximum number of points generated, range is a pseudo-random number uniformly distributed between zero and one, and round is a rounding of the entire function.
Assuming that one of the circles has radius R and coordinates (X,Y), the point coordinates (X 1 ,Y 1 ) on each circle are calculated by Equation (8): where θ is the anticlockwise angle between the line of each point and the center of the circle and has a positive value along the X-axis.
Given the randomness to ensure the positional angle, the normal vector of the convex polygon contains information on internal damage that is used to determine the scope of the convex polygons generated after the outline. The aim is to build a set of algorithms to calculate the convex polygons of each edge point with respect to its internal normal vector.
Based on the proposed algorithm, multiple aggregate models were built. The fewer edges the aggregate model had, the sharper was its shape. The model was used to simulate an elongated aggregate. The more edges the aggregate model had, the smoother the shape of the aggregate was, and the more round it tended to be. This model was used to simulate the circular aggregate.
To simulate aggregate of any shape for the numerical analysis of concrete, multiple particles can be combined in the particle flow program to form a "block." The block formed by a clump participates in the cyclic calculation as a rigid body. During the calculation process, the distance and contact force between the internal particles do not change with the steps of calculation, and the deformation of the block occurs only along the boundary. Previous studies have shown that the compressive strength of the aggregate is twice as high as that of cement mortar, and it is not damaged in the process of concrete compression. Therefore, the unbreakable unit "clump" is used as the object of generation of aggregate in the simulation. In the calculation cycle, contact between the particles inside the "block" is ignored, which reduces calculation time.
The range of sizes of the coarse aggregate used was 5~15 mm, and the Walraven Equation [11] (Equation (9)) was used to calculate its gradation. The percentage (P k ) of the volume of coarse aggregate in the concrete specimen was 75%. The average gradation curve obtained by using a slice of the non-circular aggregate concrete specimen and subjecting it to image processing using the "equivalent particle size" (D equivalent particle size = 2(S arregate area /π) 0.5 ) conformed to a 3D gradation curve and a 2D conversion curve. Therefore, the circular area of the aggregate was used as the equivalent area of the convex polygonal aggregate. Aggregate sizes of 5~7 mm, 7~9 mm, 9~11 mm, 11~13 mm, and 13~15 mm was used, and the results are shown in Table 2: where P k is the percentage by volume of coarse aggregate in the concrete specimen, P C (D < D 0 ) is the probability that the particle size D in the section is less than the size of the sieve D 0 , and D max is the maximum aggregate particle size. The number of particles within the range of each particle size is calculated according to the following equation: where n i is the number of particles in a certain range, p i is the probability that the particle size in the section is smaller than the sieve in the Walraven Equation, and A and A i are the sectional area and aggregate area, respectively. We set-up four to six polygons sides in MATLAB to generate a random aggregate model and imported it into the PFC2D model for use as a model of the concrete aggregate. The generation model is shown in Figure 2.
where Pk is the percentage by volume of coarse aggregate < D0) is the probability that the particle size D in the sectio D0, and Dmax is the maximum aggregate particle size. The number of particles within the range of each par to the following equation: where ni is the number of particles in a certain range, pi is size in the section is smaller than the sieve in the Walrave sectional area and aggregate area, respectively. We set-up four to six polygons sides in MATLAB t model and imported it into the PFC2D model for use as a The generation model is shown in Figure 2.

Verifying Model Feasibility
Because the mesoscopic parameters cannot be direc numerical simulations were used, and the results were ments to verify the feasibility of the numerical model. Ac the elastic modulus of concrete increases with Ec, and its and ̅ . Ec and kn/ks have a significant influence on the ela Poisson's ratio. The compressive strength and elastic mod

Verifying Model Feasibility
Because the mesoscopic parameters cannot be directly obtained from experiments, numerical simulations were used, and the results were compared with those of experiments to verify the feasibility of the numerical model. According to the literature review, the elastic modulus of concrete increases with E c , and its peak strength increases with σ c , and τ c . E c and k n /k s have a significant influence on the elastic modulus of concrete and its Poisson's ratio. The compressive strength and elastic modulus of the interfacial transition zone should be lower than those of cement mortar by the same proportion [12,13]. Based on this, a trial-and-error method was used to model the mesoscopic parameters following a set of mesoscopic parameters (Table 3), with the macroscopic parameters of C30 concrete as reference (30 MPa, compressive strength; elastic modulus, 30 GPa; Poisson's ratio, 0.2) to obtain the numerical macroscopic parameters σ u = 32.23 MPa, E = 28.37 GPa, and ν = 0.197. The error was very small, which shows that the numerical model is feasible. Table 3. Values of microparameters of model.

Variety Variable Value
Mortar-mortar interface

Research Methods
A preliminary analysis showed that to get a set of mapping relations between the mesoscopic parameters of concrete, 13 factors need to be considered on three test indices. Considering the nonlinear relationship between the macroscopic and microscopic parameters, the number of levels should be three or more. A long time and considerable effort are required to study the multi-level and multi-factor problems of a comprehensive test. The DOE is an efficient procedure for planning experiments so that the results obtained can be analyzed to yield valid and objective conclusions. It begins by determining the objectives of a given experiment and selecting the process variables for it. We used the DOE to arrange our test to reduce the time needed and applied it to determine the bonding between microparticles of the model during contact in a sensitivity analysis.
The test points selected by using an orthogonal test design are more representative and can give reliable research conclusions by significantly reducing the number of tests. This design can be used to screen out the microscopic parameters of concrete that have a significant influence on its macroscopic parameters, and to estimate the quantitative relationship between the two types of parameters. Although a screening test can achieve the same goal with fewer trials, it can select only two horizontal values for each factor and cannot consider nonlinear influences. Incorrect conclusions are thus easily drawn when the scope being considered is wide. We used the orthogonal experimental design for parameter sensitivity analysis.
Response surface design represents the nonlinear relationship between the test indices and factors. We used the second-order response surface equation. The experimental design idea is shown in Figure 3.

Determining the Research Interval
An orthogonal experiment was carried out to determine the first sites to arrange orthogonal test table and generate the test scheme. The mesoscopic parameters of the concrete used as reference, shown in in Table 3, were used.

Orthogonal Experimental Design
Orthogonal experiment design refers to the examination of multiple factors using experimental design method. According to orthogonality, representative points are lected from the overall test that have the characteristics of uniform dispersion and n comparison. Orthogonal experimental design is the main method of factorial des When three or more factors are involved in the experiment that may interact, the workload becomes large and difficult to implement. The orthogonal design of an exp ment is a good choice in such cases. The main tool of orthogonal experimental desig the orthogonal table. According to the requirements of the number of factors, the leve factors, and whether interaction among them occurs, the participants looked up the responding orthogonal table, and selected representative points from a comprehen test based on the orthogonality of the table. By doing so, it became possible to ach results equivalent to those of a large number of full tests using the minimum numbe tests, Therefore, the orthogonal table design is an efficient and economical multi-fa test design method. The distribution of orthogonal test points is shown in Figure 4.
According to the research interval shown in Table 4, the three factor levels were termined as shown in Table 5. The SPSS experimental design software was used to c struct the orthogonal experiment design scheme. In the analysis of variance of the orth onal experiment, leaving an appropriate blank column improves the accuracy of the sults. Thus, L643 14 was used to arrange the orthogonal experiment table and the 14th umn was left out. The test scheme and results are shown in Table 6. Table 4. Research interval of microparameters.

Variety Variable Value Range
Mortar-mortar interface Figure 3. Idea of the experimental design.

Determining the Research Interval
An orthogonal experiment was carried out to determine the first sites to arrange an orthogonal test table and generate the test scheme. The mesoscopic parameters of the C30 concrete used as reference, shown in in Table 3, were used.

Orthogonal Experimental Design
Orthogonal experiment design refers to the examination of multiple factors using the experimental design method. According to orthogonality, representative points are selected from the overall test that have the characteristics of uniform dispersion and neat comparison. Orthogonal experimental design is the main method of factorial design. When three or more factors are involved in the experiment that may interact, the test workload becomes large and difficult to implement. The orthogonal design of an experiment is a good choice in such cases. The main tool of orthogonal experimental design is the orthogonal table. According to the requirements of the number of factors, the level of factors, and whether interaction among them occurs, the participants looked up the corresponding orthogonal table, and selected representative points from a comprehensive test based on the orthogonality of the table. By doing so, it became possible to achieve results equivalent to those of a large number of full tests using the minimum number of tests, Therefore, the orthogonal table design is an efficient and economical multi-factor test design method. The distribution of orthogonal test points is shown in Figure 4.   According to the research interval shown in Table 4, the three factor levels were determined as shown in Table 5. The SPSS experimental design software was used to construct the orthogonal experiment design scheme. In the analysis of variance of the orthogonal experiment, leaving an appropriate blank column improves the accuracy of the results. Thus, L 64 3 14 was used to arrange the orthogonal experiment table and the 14th column was left out. The test scheme and results are shown in Table 6. Table 4. Research interval of microparameters.

Variety Variable Value Range
Mortar-mortar interface Aggregate-mortar interface  Table 6. Orthogonal test scheme and results.

Range Analysis
Range analysis is also called the direct analysis method. It has the advantages of simple calculation, a visual representation, and the ease of understanding of the results. It is the most commonly used method for analyzing the results of orthogonal experiments. The trends of changes in the test index average effect are shown in Figures 5-7. A simple calculation can be used to obtain the influence of the factors on the test index. Range analysis was carried out using Minitab, and the results were plotted (Figures 8-10).

Range Analysis
Range analysis is also called the direct analysis method. It has the advantages of simple calculation, a visual representation, and the ease of understanding of the results. It is the most commonly used method for analyzing the results of orthogonal experiments. The trends of changes in the test index average effect are shown in Figures 5-7. A simple calculation can be used to obtain the influence of the factors on the test index. Range analysis was carried out using Minitab, and the results were plotted (Figures 8-10).              The results show that σ c and τ c had a significant influence on σ u because their sizes determined the difficulty of relative sliding between particles after bond failure. E c and k n /k s , which are directly related to the bond stiffness coefficient, had a significant influence on E, k n /k s and µ influenced ν.

Analysis of Variance
Range analysis of the orthogonal experimental design has the advantages of simplicity and requiring a small number of calculations, but it cannot estimate the size of the error or the importance of the influence each factor on the results. To compensate for the deficiencies of range analysis, the analysis of variance was used to analyze the results of the orthogonal experiment (Table 7). The results of the analysis of variance were as follows: (1) σ c and τ c had a significant influence on σ u while the influence of other parameters was not significant. (2) E c , k n /k s , µ, (E c ) j and (k n /k s ) j had a significant influence on E while the influence of other parameters was not significant. (3) E c , k n /k s , µ and (k n /k s ) j had a significant influence on ν in the interfacial transition zone of the mortar particles.
Combined with the above analysis, the relationship can be simplified as follows: The degree of fitting R 2 of each fitting equation was between 0.7598 and 0.9777, indicating that the fitting effect was relatively good. However, due to limitations on the number and form of distribution of test points in the orthogonal experimental design, the Equation in Table 8 cannot accurately predict the macro-micro parameter relationship beyond the test points. Therefore, to increase the accuracy of the fitting Equation, another numerical test was needed.

Response Surface Design
Response surface design is a statistical method used to solve multi-variable problems by using reasonable experimental design methods, obtaining certain data through experiments, and using multiple quadratic regression equations to fit the functional relationships between the factors and the response values. and requires two to seven factors for a range of tests to identify nonlinear relationships (second order). It can be used to carry out tests with two to seven factors, and to estimate the nonlinear relationship (second order) between the test index and the factors. The fitting equation is as shown in Equation (14): where β is an undetermined parameter, x k is the k th independent variable, p is the number of factors, and ε is the error term.
To calculate the test factors and indicators in general, the quadratic regression equation between the levels of each factor number should be greater than or equal to three. Using a comprehensive design to estimate the above undetermined parameters leads to the problem of too much residual freedom and too high a cost. Response surface design, by testing only reasonable arrangements, avoids these problems. The Box-Behnken test is a common response surface design.

Central Composite Test
As shown in Figure 11, the test point designed for the two-factor center composite test consisted of three parts: cubic points, axial points, and center points. Cubic points represent the two-level factorial design, center points represent the nonlinear test and error estimates, and axial points consider conditions of the nonlinear effect.
To enable the design to meet the requirements of rotation and sequencing, the factor level of the axial point is usually calculated according to Equation (15) [13]: where α is the number of factors.
To enable the design to meet the requirements of rotation and sequencing, the level of the axial point is usually calculated according to Equation (15) [13]: where α is the number of factors.

Box-Behnken Experimental Design
Compared with the central composite test mentioned above, the Box-Behnk has fewer tests designed with the same number of factors. There is no axial point test design, because of which its horizontal setting does not exceed the range of safe ation. Axial points generated by the central composite test with axial points may the error threshold for safe operation area. As is shown in Figure 12, the Box-Behnk takes each test point as the middle point of each edge of the cube and uses the cente to estimate the test error.

Box-Behnken Experimental Design
Compared with the central composite test mentioned above, the Box-Behnken test has fewer tests designed with the same number of factors. There is no axial point in the test design, because of which its horizontal setting does not exceed the range of safe operation. Axial points generated by the central composite test with axial points may exceed the error threshold for safe operation area. As is shown in Figure 12  To enable the design to meet the requirements of rotation and sequencing, the fa level of the axial point is usually calculated according to Equation (15) [13]: where α is the number of factors.

Box-Behnken Experimental Design
Compared with the central composite test mentioned above, the Box-Behnken has fewer tests designed with the same number of factors. There is no axial point in test design, because of which its horizontal setting does not exceed the range of safe op ation. Axial points generated by the central composite test with axial points may exc the error threshold for safe operation area. As is shown in Figure 12, the Box-Behnken takes each test point as the middle point of each edge of the cube and uses the center p to estimate the test error.

Test Plan
According to Equation (11), two factors that had a significant impact on the peak intensity corresponding to its response surface. Therefore, a two-factor CCD test table with 13 test sites (Table 9) was used for the experiment. According to Equations (12) and (13), the elastic modulus and Poisson's ratio had significant impact factors of five and four. For the response surface corresponding to the elastic modulus and Poisson's ratio, five and four variables were considered respectively. Considering that the axial point generated by CCD exceeded the actual value represented by the factors, a five-factor BBD test table (Table 6) with 46 test points and a four-factor BBD test table (Table 7) with 27 test points were selected for the experiment. In Table 9, point type "0" represents the center, "−1" represents the axial point, and "1" represents the cubic point. In Tables 10 and 11, point type "0" represents the center point and type "2" represents the middle point. The other factors that had no significant influence on σ u , E and ν are not discussed. Their values were (E c ) g = 30 GPa, (k n /k s ) g = 2, µ g = 1, (σ c ) j = 15 MPa, (τ c ) j = 15 MPa, and µ j = 1.

Quantitative Equation
We used the test design scheme of Minitab to design the response surface experiment (Tables 9-11), and the results were imported for response surface analysis to obtain the second-order response surface Equation of σ u , E and ν, as shown in Table 12.

Domain of Definition
The table shows that the degree of fitting of each response surface σ u , E and ν was greater than 0.988, indicating a good fitting effect and adequate reflection of information on the test site.

Calibration Program for Microscopic Parameters
The mesoscopic parameters of the PFC2D model cannot be obtained directly. Multiple numerical simulations are needed to determine them corresponding to the macroscopic physical parameters through continual selection and trial calculations to establish the relationship between them. This process is called parameter calibration. Microscopic parameter calibration is essentially a mapping problem between macroscopic and microscopic parameters. As long as a qualitative or quantitative relationship between them can be identified, the macroscopic parameters of the model can be calculated from its microscopic parameters, and the latter can be matched based on the former. We use an optimization method to calibrate the mesoscopic parameters. Based on the macro -meso parameters obtained from the previous section, the calibration problem is transformed into an optimization problem on the basis of the quantitative relationship between the mesoscopic parameters. The characteristics of the objective function and constraints on it are used to choose the appropriate optimization algorithm.

Optimization Model and Algorithm
The purpose of mesoscopic parameter calibration is to find a set of appropriate mesoscopic parameters for the particle discrete element model for concrete. We assumed that the macroscopic parameters of the granular discrete element concrete model derived from an empirical Equation or a quantitative relation are y i * (i = 1, 2, 3 . . . . . . ), and those of concrete measured through laboratory tests are y i lab (I = 1, 2, 3 . . . ,s). Then, the goal of mesoscopic parameter calibration can be converted into finding a set of mesoscopic parameters x 1 , x 2 , x 3 , . . . , x k (k is the total number of microscopic parameters involved) such that the difference between y i * and y i lab is as small as possible. Considering the different orders of magnitude of various macroscopic parameters, the dimensionless number for relative deviation |y i * − y i lab |/y i lab is used to describe the difference between y i * and y i lab . The ultimate goal is to find a set of microscopic parameters x 1 , x 2 , x 3 , . . . , x k that minimizes the relative deviation of all macroscopic parameters.
Finally, we set the range of values of each mesoscopic parameter to [a j ,b j ] (j = 1, 2, 3, . . . ,k). The general form of the calibration model for mesoscopic parameters is then given by: Objective function: Constraint: a j ≤ x j ≤ b j (1 ≤ j ≤ k) According to the above numerical test and analysis, the fitting Equation for the macroscopic parameters σ u , E and ν in the specified domain were obtained ( Table 12). The microscopic parameters used as variables included E c , k n /k s , σ c , τ c , µ, (E c ) j and (k n /k s ) j . We set σ u , E and ν as macroscopic parameters for y1, y2, and y3, and used the parameters E c , k n /k s , σ c , τ c , µ, (E c ) j and (k n /k s ) j for x1, x2, x3, x4, x5, x6, and x7 to obtain the optimal model as follows: Objective function:

Implementation of Microscopic Parameter Calibration Program
MATLAB is among the most widely used scientific computing software, and is powerful and easy to use. It provides the user with the freedom of an interactive programming interface as well as strong functional encapsulation in the form of a toolbox for users. Our work here required the MATLAB's built-in genetic algorithm toolbox.
The algorithm used was the default one for the fitness function in the genetic algorithm toolbox. The lower the smallest value of the individual fitness value is, the close the given individual is to the optimum solution. The fitness function of the problem of microscopic parameter calibration has the same form as the objective function in Equation (17), which can be expressed as Equation (19): Fitness function: is the fitness value, σ u lab is a constant obtained by the peak intensity of an indoor uniaxial compression test of rocks, E lab is a constant representing the elastic modulus measured using the same test, ν lab is a constant representing Poisson's ratio, and σ u * , E * and ν * are dependent variables calculated by the quantitative Equation in Table 12.
The genetic algorithm toolbox allows the user to give equality, inequality, and boundary constraints. According to the optimization model given in Equation (17), only the boundary constraints were set to limit the upper and lower limits of each mesoscopic parameter, as shown in Equation (20): The evolutionary process in the genetic algorithm is one of iterative optimization completed by the evolution of the population from one generation to the next. Each generation of population consists of several individuals, and the "excellence" of individuals or the fitness of the environment is evaluated by the fitness function. "Excellent" individuals are chosen for the crossover and mutation operator to form a new generation of the population, whereas suboptimal individuals are "erased." Each generation forms a new group in this way. The optimal solution is obtained when the process of evolution of the individuals is completed.
The process of solving the optimization model shown in Equation (17) by using the genetic algorithm is described as follows ( Figure 13 (1) The initial population {x} 0 is generated, and N individuals are randomly generated (300 in this paper).  (1) The initial population {x} 0 is generated, and N individuals are randomly generated (300 in this paper). (2) Fitness evaluation is performed for each individual in the population: 1 Assume that k = 1. MATLAB's built-in graphical user interface (GUI) module can help users write interactive programs. Based on the MATLAB environment, the author wrote the mesoscopic parameter calibration program and designed its interface based on the GUI module. To enable the program to run independently of MATLAB's operational environment, the MATLAB complier module was used to encapsulate the above program, and its icon and interface were designed and added. Finally, a calibration program for the microscopic parameters that could run independently was obtained.

Verifying Effect of Calibration Program
Tables 9 and 10 are used to assess the range of calculation of the microscopic parameter calibration program, which is shown in Table 13. Table 13. Calculation range of the microscopic parameter calibration program. Five groups of macroscopic parameters were selected for parameter calibration, and the corresponding optimal solution was obtained from the above program (Table 14). When the obtained microscopic parameters were input into the DEM, the random characteristics of particle structure and aggregate distribution led to differences in the results of the simulation even if the microscopic parameters of the model remained consistent. The results were thus calculated five times to ensure the accuracy of the verification test. Table 15 shows that the coefficients of variation of the various macroscopic parameters measured by the uniaxial compression numerical test of the concrete model were between 2.0928% and 6.3849%, indicating that the calculated mechanical parameter calculations were stable. The relative error between the measured macroscopic parameters and the target macroscopic parameters is between 0.0093%~6.4246%, which has achieved a very good calibration effect. This verified the effect of calibration and precision of the program to calibrate the microscopic parameters.

Conclusions
This paper proposed and tested a discrete element model of concrete using the particle discrete element software PFC2D. Based on the ideas of an orthogonal experimental design and response surface design, a set of mapping relationships between the macro-micro parameters of concrete were obtained. The MATLAB environment was used to write a procedure to calibrate the mesoscopic parameters. The work here can serve as a reference for future research. The main conclusions of this study are as follows: (1) An orthogonal test of the analysis of variance showed that σ c and τ c had a significant influence on σ u , while the influence of the other parameters was not significant. E c , k n /k s , µ (E c ) j and (k n /k s ) j had a significant influence on E, and E c , k n /k s , µ and (k n /k s ) j had a significant influence on ν in the interfacial transition zone of mortar particles, while the influence of the other parameters was not significant. (2) The response surface design method was adopted to obtain a set of quantitative relations of the macro -meso parameters of concrete within a certain range that accurately described the nonlinear relationship between them. In addition, by increasing their range of values, the mesoscopic parameters could be calibrated more precisely and, thus, more accurately inverted. (3) The selection of an appropriate contact model of concrete, test design method, mathematical model, and optimization algorithm is key to achieving good parameter calibration.
The research results of this paper aim to provide a parameter inversion method. Due to reasons of time, test equipment and so on, no real concrete block testing was carried out in the laboratory, and only the numerical simulation part was done. In this paper, for the convenience of calculation, some parameters of the model are selected in advance, and only 13 meso-parameters within a certain range were selected as test factors. If possible, we will continue to study the relationship between other meso-parameters and macro-parameters in the future. In addition, this paper only selects concrete in the case of aggregate gradation for research, which only provides a method of mutual inversion between macro and micro parameters of concrete. If other aggregate gradations are used, this model is also applicable according to the research ideas in this paper. If we simulate materials such as reinforced concrete, more problems need to be considered. Although the modeling will be more complex, the method in this paper is also applicable. In the future, we plan to conduct further research to solve more practical problems.