Research on back analysis of meso - parameters of hydraulic cemented sand and gravel based on Box - Behnken design response surface

: Cemented sand and gravel ( CSG ) is a new type of dam - building material. Aiming at the cumbersome process and long calculation time of traditional methods to invert the meso - parameters, a mesophase parameter inversion method based on Box - Behnken Design response surface was proposed. By constructing a response surface simulation test scheme with di ﬀ erent inversion parameters ( elastic modulus of aggregates, mortars and interfaces, and interfacial tensile strength ) , the stochastic aggregate model is used to complete the numerical simulation of the damage process, and di ﬀ erent results are obtained. The equation between the response variable ( stress at di ﬀ erent loading times ) and the independent variable ( inversion parameter ) was veri ﬁ ed, and the rationality of the response model was veri ﬁ ed; the action mechanism of mesophase parameters at di ﬀ erent stages on the mechanical properties of the specimen was analysed. The test results are brought into the response surface model, and the meso - parameters are obtained by inverse analysis. The stress – strain curve obtained by numerical simulation with this parameter has an error of 1.1% at the peak stress and 3.27% at the peak strain. The accuracy is high, the number of test groups is much smaller than other conventional inversion methods, and has feasibility of application in CSG.


Introduction
Cemented sand and gravel (CSG) is a composite material consisting of water, non-screened aggregates, and a small amount of cementitious material [1], and its mechanical properties determine the usability of CSG [2][3][4][5][6][7]. Both CSG materials and concrete are three-phase composites [8][9][10], consisting of aggregate phase, mortar phase, and interface phase. However, the non-screening and non-washing characteristics of CSG aggregate make its mechanical properties different from traditional concrete materials. And the aggregates are selected locally, so the properties cannot be generalised. It is necessary to study the difference between CSG and ordinary concrete materials. The meso-analysis method can accurately grasp the difference and interaction between CSG and its components, especially the mechanical properties of the aggregate and mortar, and can accurately determine the force characteristics and failure conditions of CSG. It is a method commonly used by most researchers. The acquisition of meso-component parameters is particularly important [11].
CSG has little research on its parameters and parameter acquisition as a new building material. Scholars [12] have used many uniaxial compression and loading and unloading experiments of CSG to use the linear elastic constitutive model as their constitutive model. The parameters are obtained by experimental means [13]. Batmaz [14] simulated the nonlinear mechanical characteristics of CSG by piecewise linear method and applied the method and model to actual engineering. More research on the mechanical properties of CSG is still in macroscopic experiments [15][16][17], and there are fewer meso-chromatographic analyses.
Due to the scale effect and the fractal [18], the mesocomponents of the materials are often obtained by back-analysis methods. Genetic algorithms [19][20][21][22][23], for example, are prone to "prematureness" and require independent coding, which does not guarantee both accuracy and simplicity. At the same time, support vector machine models [24] and particle swarm algorithms [25,26] are also more complex. The response surface method is an optimisation method that combines the experimental design method and mathematical statistics proposed by Box and Wilson [27], and has good robustness [28,29]. Pei et al. [30] combined the response surface model and the genetic algorithm to perform inversion analysis of the thermodynamic parameters; Li et al. [31] verified that the response surface method is an optimal method of mathematical modelling, which can effectively reduce the number of experiments and can investigate the interaction between influencing factors. Li [32] and Wang [33] established response surfaces of different orders and analysed the concrete meso-parameters in inversion analysis, indicating that the inversion method of geotechnical materials has high prediction accuracy.
The documents mentioned above confirm the feasibility of a response surface in engineering back analysis. Still, it has not been studied to use it to obtain the mesocomponent parameters of CSG materials. This article solves the non-simultaneous response surface equations obtained from the BBD simulation test scheme by numerical simulation-response surface method experiment to obtain the fine-scale parameters of the test blocks, make theoretical pavement, and case study for the experimental study of CSG meso-parameters inversion, and propose new ideas for engineering applications.

Parametric inverse analysis 2.1 Principles of the response surface method
Response surface methodology is an experimental improvement and model optimisation method. It works by finding representative points in the global range, testing the representative points, obtaining the resultant function, and then explicitly fitting the global range of input variables (parameters) to the response variables through regression. A large number of applications are available [34][35][36].
According to Taylor's polynomial theory, the response variable (dependent variable) can be fitted with a polynomial as a function of the input variable (independent variable). A second-order response surface model is used on synthesis, represented by equation (1).
where equation y represents the response, m represents the number of independent variables, that is, the number of controlled factors, x i represents the i independent variable, and β 0 represents the constant term coefficient, β i represents the coefficient of the first term, β ij represents the coefficient of the quadratic term and the quadratic cross term, and ε represents the random error.

Basic steps of inversion of the response surface method
BBD is a commonly used experimental design method for response surface optimisation, which is suitable for exploring the relationship between 2 and 5 factors. Therefore, this article adopts the BBD experimental design, and selects three levels of factors, namely, the maximum value, the minimum value, and the median value. The steps of parameter inversion based on the response surface method are [37]: (1) Determine the inversion parameters and other basic parameters. (2) Determine the parameter levels of the independent variables, and design the simulation experiment scheme. (3) Complete the numerical simulation calculation process, construct the response surface function, and verify the rationality of the response surface equation. (4) The laboratory test results are used as response variables, and are brought into the response surface model function for optimisation and solution, and the parameters to be inverted are obtained. (5) Perform mesoscopic numerical simulation using the obtained parameters to obtain the simulation results. The degree of agreement with the laboratory test results is compared and analysed to verify the rationality of the inversion parameters.

Analysis on the reasonableness of response surface model
The response surface model reflects the implicit functional relationship between the independent and response variables. Therefore, whether there is a statistical relationship in the mathematical sense between the independent variable and the response variable requires a significance test, analysis of variance, and a precision test.
Whether the constructed response surface model is reasonable or not needs to be tested. The control parameters for model testing are: R 2 and adjusted R-Squared (R adj 2 ), F-test method, P-value test method, residual distribution test method, etc. R 2 is used to evaluate the goodness of fit of the model and describes how well the input variables explain the output variables; the closer R 2 is to 1, the better the model is. It is calculated as follows: where T is the number of check points in the design space, x t rs is the target value of the response surface calculation of the t-th parameter, and x t p is the target value of the finite element model calculation of the t-th parameter, x t is the average value of the finite element model and the calculation results of the finite element model. The adjusted R 2 (R adj ) imposes penalties on the number of items that increase but do not improve the accuracy of the model, to increase the simplicity of the model while ensuring the reliability of the model. It is usually analysed together with the predictive determination coefficient (R pred ). R pred reflects the predictive ability of the model and requires that the difference between the two does not exceed 0.2 as shown in the following equations (3) and (4): The adjusted R 2 (R adj ) increases model simplicity while ensuring model reliability by imposing a penalty on the number of terms that increase but do not improve the accuracy of the model. It is usually analysed in conjunction with the predicted R 2 (R pred ), which reflects the predictive power of the model and requires a difference of no more than 0.2 as shown in equations (5) and (6).
where P is the number of features, R adj lies between 0 and 1, and the larger the better, but if it is really close to 1, the over-fitting situation may occur.
Use the F-test method for hypothesis testing: where SSR is the regression sum of squares of R, SSE is the sum of squares of the errors. p is the number of nonconstant terms in the RSM function, and p and ( ) − − m p 1 are the degrees of freedom of SSR and SSE, respectively. For a given confidence level α, the RSM model is consid- The signal-to-noise ratio, or Adeq Precision, is used to analyse the model's interference items further when the model's regression items are obtained. As long as the accuracy of the model is greater than 4, it is generally considered that the accuracy of the model is reasonable, and the calculation formula is shown in equation (8): where S represents the output value of the independent variable, and N represents the error value caused by the independent variable in the calculation process.

Analysis of calculation examples
3.1 Numerical simulation tests

Stochastic aggregate model
With the support of ANSYS's APDL for secondary development, the model was created from a meso-level perspective treating the CSG material as a three-phase composite model consisting of natural sand and gravel aggregate, mortar matrix, and the interface between the mortar matrix and the aggregate, using 100 × 100 × 100 test blocks with a calculated aggregate volume rate of 70% and a secondary mix of aggregates with Fuller's gradation, using the Walraven formula. As in equation (7), the number of aggregates in different intervals is generated by deriving the probability of occurrence of an inner truncated circle P C (D < D 0 ) with an aggregate diameter of D < D 0 at any point, through the threedimensional Fullerian gradation of the aggregate, from a twodimensional cross-section. where D 0 is the diameter of the sieve hole, mm, D max is the maximum aggregate particle size, mm, D min is the minimum size of the aggregate, and P K is the percentage of aggregate volume, which is taken as 0.7. Considering the non-sieving characteristics of the aggregate of CSG, the finite element model (FEM) 1 established by the random aggregate model is shown in Figure 1.

Realisation of the uniaxial compression process
The constitutive numerical simulation model adopts the linear elastic model, and the failure theory adopts the first strength theory. The uniaxial compression simulation test uses the displacement loading method, the same as the indoor test using the loading method. The analysis type is static analysis. The automatic time step is used to control the load substeps and kill the elements whose tensile stress exceeds the tensile strength of each phase material to simulate the failure process.

Establishment of response surface model and reliability analysis
Response surface design is based on design-export software, and image processing. The initial parameters were selected concerning the results of the sensitivity analysis of the relevant literature [38]. The linear elastic model was adopted according to the fine view, and the modulus of elasticity of the aggregate (E G ), the modulus of elasticity of the mortar (E m ), the modulus of elasticity of the interface (E i ), and the tensile strength of the interface (f i ) were selected as input variables. The parameters were taken according to three levels and the four factors and three levels were taken in Table 1. Since the indoor test process used displacement loading to obtain uniaxial compression stress-strain curves, and the parameters to be sought were four, the stress values at four points (corresponding after different loading steps) were used as the response quantities, and the strain values corresponding to the above response quantities are shown in Table 2.
where ε 1 , ε 2 , and ε 3 are located before the peak strain and ε max is the peak strain. Other parameters include aggregate Poisson's ratio (μ G ), mortar Poisson's ratio (μ m ), interface Poisson's ratio (μ i ), aggregate tensile strength (f g ), and mortar tensile strength (f m ). These four quantities are less sensitive to the overall intensity of CSG. Referring to the relevant literature and based on the parameters to be inverted, the values of other parameters are shown in Table 3.
The response surface simulation test design is shown in Table 4, and 29 sets of simulation test programs are obtained.
Performing uniaxial compression numerical simulation experiments based on the parameters in the table above, 29 sets of peak stress and stress-strain curves are obtained, and a quadratic polynomial with cross terms as the response surface function is used. For convenience, the response surface function based on the coded factor is expressed as follows: where f f f f , , , 1 2 3 max are the corresponding stress values at ε 1 , ε ε , , 2 3 and ε max , respectively. A is the modulus of elasticity of the aggregate, B is the modulus of elasticity of the mortar, C is the modulus of elasticity of the interface, and D is the tensile strength of the interface.
In the process of obtaining the four response surface equations, the regression equation coefficients (R squared ), regression equation prediction coefficients (R pred ), regression equation correction coefficients (R pred ) and p-values (p-test) and F-values (F-test) were analysed to remove unreasonable terms, adjust out-of-fit terms, compare R adj with R pred , etc. Polynomial processing was carried out to correct the response surface model and determine the model's overall reliability and accuracy, as shown in Table 5.
It can be seen from the table that the p-value and Fvalue are reasonable, and the four response surface models are overall significant, so the parameter selection is reliable. In addition, the difference between R adj and R pred is less than 0.2, indicating that the response surface model has good predictive ability, and all the prediction accuracy is greater than 4. Therefore, it is considered that the model accuracy meets the expected requirements, and the four response surface models are overall reliable. The residual normal distribution probability diagram of the response surface model is shown in Figure 2.
The residual normal distribution plot shows that the data fall on a straight line, which means that the errors are distributed regularly. In addition, the distribution of the residuals shows that the residuals are normally random, which shows the accuracy of the model [37].     The level of the parameters reflects the sensitivity of the meso-parameters, which is consistent with the initial assumptions and the verification conclusions of the cited literature. It further illustrates the rationality of the response surface model and judges that the control items are reasonable.

Response surface model analysis
According to the response surface model (equations (10)- (14)), a 1D graph is drawn corresponding to the intermediate level value of each parameter and the influence of the primary term of each variable on the overall stress value is analysed; a 3D response surface diagram of the interaction of factors is drawn, and the interaction of various parameters is observed visually (Figures 3 and 4).
As can be seen from the figure, σ σ and 1 2 have similar single-factor action plots and interaction term action plots, indicating that the stress-strain relationship is stable in the early stage of the test, and it can be observed   that the stress magnitude of the model is basically determined only by the relationship with E G , E m , and E i at small strains, and the degree of influence is: E m > E G > E i ; that is, at small strains, the stress magnitude of CSG mainly related to E G , E m , and E i is positively correlated, and the relationship with f i is not significant, indicating that the relative deformation between different components is small at this stage, and the interface is subject to little tension, and the material stress is mainly determined by the modulus of elasticity of each component, while under small strains, although the different modulus of elasticity of each component will cause the deformation difference, and thus produce tensile stress, the effect of tensile stress is not obvious because the tensile stress between materials does not exceed the tensile strength of each component. However, because the parameter taking interval is relatively small and the coefficients of each variable are also small, the increase in σ is not obvious as E G , E m , and E i increase ( Figure 5).
During the transition from σ 1 and σ 2 to σ 3 , the degree of σ affected by E m > E G > E i remains unchanged, but σ 3 is more affected by f i ; a 3D response surface diagram of the  interaction of factors is drawn, and the interaction of various parameters are visually observe as shown in Figure 6. During the transition from σ 1 and σ 2 to σ 3 , the development trends of the interaction terms AB and BC remain unchanged, and the overall trends of the interaction terms AD and BD are the same. σ 3 increases with the increase in f i and E i . σ increases with the increase in the four parameters as a whole, that is, as a whole, the elastic modulus of the model increases with the increase in E G , E m , E i , and f i .
It can be seen visually on Figures 7 and 8 that the two terms E i and f i have the greatest influence on the model σ max and are positively correlated, with f i having a greater influence than E i ; this is inversely proportional to E g and E m because the damage criterion uses the maximum tensile stress damage criterion and material damage is caused by insufficient interface deformation and small tensile strength leading to interface damage first. In contrast, the aggregate and mortar damage is minimal (due to the general literature) [39,40].

Indoor test verification analysis
According to the test procedure, three 100 × 100 × 100 test blocks were made. The cement used is ordinary Portland cement (P·O42.5), the aggregate is sand and gravel from a river bed, and the fly ash is purchased from a power plant class I ash. Standard grade II aggregates are used. The test coordination ratio is shown in Table 6. The specimens were cured for 28 days for the uniaxial compression test. In the test, the universal testing machine and displacement control method were used, the rate was 0.6 mm/min. Three groups of specimens were made for uniaxial compression test, and the test results were in accordance with the specification [41], so the test results were considered reliable. According to the specification, a stress-strain curve is selected as the inversion curve, as shown in Figure 10.
The measured values are as follows in Table 7. Bringing the measured values in the above table as the inversion target into the response surface model, which requires a guaranteed rate of more than 97%, the     Table 8.
To verify the accuracy of parameter inversion, two groups of finite element models with the same gradation as in Figure 1 are re-established, as shown in Figure 9. The abovementioned inversion mesoscopic parameters are used to bring into three finite element models to simulate the failure process of the indoor uniaxial compression test.
The stress-strain relationship curve is obtained and compared with the laboratory test results, as shown in Figure 10.
It can be seen from Figure 10 that the stress-strain curves of the three numerical models are in good agreement with the indoor test curves, and the results are reliable.
The error calculation results are shown in Table 9. The analysis found that the strength of CSG was significantly lower than that of ordinary concrete materials [42]. The maximum error of peak strain is 6.182%, and the maximum error of peak stress is 3.471%, which appear in different finite element models, and the overall error is small. The optimal experimental group is FEM 1, the inversion model has high accuracy, and the inversion results can be considered reliable. Compared with other inversion methods [43,44], the above four-parameter inversion analysis method of response surface only uses 29 sets of numerical simulation schemes, which greatly reduces the number of calculation sets and improves the calculation efficiency.

Conclusion
This article uses the indoor test-response surface methodfinite element test to complete the material parameter analysis of the mesophase of cementitious sand and gravel materials. The reliability of the inversion parameters is verified, and the specific conclusions are as follows: (1) Based on the response surface theory, the mesoscopic phase parameter inversion method proposed in this article is feasible. Compared with the traditional inversion method, the experimental design of this method is simple, and the number of inversion groups is small. Provides a reference for the interdisciplinary application of response surface methods. (2) It can be seen from the response surface model that in the process of reaching the peak stress, the relationship of each component to the mechanical properties of the material is also changing. The initial stress is mainly positively correlated with E g , E m , and E i , and f i is not significant. The peak stress is most affected by f i and E i , which is a positive correlation, and has a negative correlation with E g and E m as a whole. With the increase in strain, each phase has different elastic modulus, so different deformation occurs, and tensile failure occurs. (3) The parameter inversion results are brought into the finite element model. Compared with the laboratory test, the stress-strain curve is generally consistent.
The optimal group has an error of 1.1% at the peak stress and 3.27% at the peak strain, with high accuracy and reliable inversion results. The rationality of the response surface method inversion is verified, which lays a foundation for the subsequent theoretical research based on this method.

Acknowledgments:
We would like to submit the enclosed manuscript entitled "Research on Back Analysis of Mesoparameters of Hydraulic Cemented Sand and Gravel Based on Box-Behnken Design Response Surface," which we wish to be considered for publication in "Science and Engineering of Composite Materials." No conflict of interest exists in the submission of this manuscript, and the manuscript is approved by all authors for publication. I would like to declare on behalf of my co-authors that the work described was original research that has not been published previously, and not under consideration for publication elsewhere, in whole or in part. The data is authentic and reliable, and the code is usable. All the authors listed have approved the manuscript that is enclosed.  Conflict of interest: The authors declare no conflict of interest.
Data availability statement: All data, models, and code generated or used during the study appear in the published article.