Optimization of Ball-Milling Parameters for Enhanced Compressive Strength of Coal Gangue Cement Paste: A Response Surface Methodology Approach

The production of cement, the most widely used cementitious material, is associated with the emission of approximately 0.913 tons of carbon dioxide per ton of cement produced, significantly contributing to global greenhouse gas emissions and escalating energy consumption worldwide [1]. The search for alternative materials to cement is therefore crucial in mitigating the environmental impacts of its production. Coal gangue, a by-product of coal mining and cleaning processes, constitutes approximately one-tenth of the total coal quantity [2]. The accumulation of coal gangue, which may contain heavy metal ions, poses a substantial environmental threat, with estimates suggesting that the cumulative amount could exceed 700 million tons without appropriate control and treatment measures. Given its composition, predominantly silicon-aluminum, coal gangue presents itself as a promising alternative to cement, offering significant potential for reducing cement consumption and promoting environmental sustainability [3], [4], [5].

Research has indicated that CGP exhibits promising activity levels as a cement substitute. Guo et al. [6] reported that the activity of CGP, at 102%, surpasses that of red mud, making it a more suitable replacement for cement in solid waste applications. However, as noted by Liu et al. [7], the reactivity of low-kaolinite calcined coal gangue, which is similar to siliceous fly ashes, is lower than that of calcareous fly ashes (CFA), natural pozzolan (Po), ground granular blast furnace slags (S), and calcined clay (CC), which impedes the development of strength in blended mortars. Nevertheless, it has been demonstrated that the co-utilization of coal gangue and limestone in a 2:1 ratio can enhance later-stage strength by 5-12% compared to corresponding calcined coal gangue mortars [7].

Enhancing the physical activity of coal gangue is considered a more effective approach to improving the performance of coal gangue mortar. This improvement also facilitates chemical synergy with other calcium siliceous materials, as highlighted by Liu et al. [7]. The mechanical activation of coal gangue primarily influences its physical properties, particularly particle size distribution and morphology. Coal gangue, a high-silicon, low-calcium material, is typically composed of small particles (0.5 cm-1.5 cm) and consists mainly of stable kaolin minerals. To be effectively utilized as an auxiliary cementitious material, coal gangue must be ground into a fine powder. The physical properties of CGP, however, vary depending on the grinding parameters employed, making the selection of optimal grinding parameters critical to achieving the most cost-effective and efficient use of coal gangue as a cement replacement.

This study investigates the relationship between various grinding parameters and the compressive strength of coal gangue cement paste using RSM. RSM, a combination of mathematical, statistical, and experimental design methods [8], [9], [10], is particularly well-suited for exploring complex mathematical models between response outputs (response variables) and influencing factors (variables) in systems with unknown behavior. RSM has been widely applied in civil engineering, as demonstrated by Tian et al. [11], who utilized RSM coupled with a multi-objective particle swarm optimization algorithm to predict and optimize the variables affecting the initial fluidity and mechanical strength of road geopolymer grouting material. Their results showed that optimal solutions could be predicted and verified with only a 5% deviation from experimental results. Hou et al. [12] also employed an RSM model to establish functional relationships between variables (water-cement ratio, magnesia-to-phosphate ratio, and borax content) and response variables (setting time and early-age compressive strength) in the development of magnesium phosphate cement (MPC)-based materials for subway segment repair. Similarly, Mtarfi et al. [13] applied RSM to optimize the effects of fly ash and grinding agents on the compressive strength of cement at 7, 28, and 90 days, ultimately obtaining the optimal mix design for the paste.

Therefore, this research employs RSM to investigate the effects of grinding parameters (ball-to-material ratio, milling time, and milling speed) on the compressive strength of coal gangue cement paste at different substitution rates. The study’s findings, based on RSM analysis and micro-experiments, reveal the influence of these parameters on the properties of CGP and the mechanical performance of coal gangue cement paste. Finally, numerical optimization through RSM is utilized to identify affordable and feasible grinding parameters that promote the effective utilization of coal gangue as a sustainable construction material.

The particle size analyzer was used for material characterization of CGP to observe the changes in grinding parameters on the particle size of CGP. A scanning electron microscope (FlexSEM1000) was used for observing the surface topography of CGP. The acceleration voltage of the scanning electron microscope was set to 5.0 kV, and the magnification was set to 50 times.

The cement (P.O. 42.5) is sourced from the Tian Jia An cement plant located in Huainan City, Anhui Province, China. Its chemical composition is shown in Table 1. $\mathbf{S i O}_{\mathbf{2}}$+$\mathbf{C a O}$+$\mathbf{A l}_{\mathbf{2}} \mathbf{O}_{\mathbf{3}}$ occupies about 90.3% of the total mass. The particle size distribution of cement is shown in Figure 1, and the median particle size is 13.65$\mathrm{\mu m}$.

The coal gangue was sourced from a coal mine in Huaibei City, Anhui Province, China. The size of the primary coal gangue is about 1 cm to 1.5 cm, with a smooth surface and sharp edges and corners (Figure 2). $\mathbf{S i O}_{\mathbf{2}}$+$\mathbf{C a O}$+$\mathbf{A l}_{\mathbf{2}} \mathbf{O}_{\mathbf{3}}$ accounted for 93.85% of the total mass, which is similar to cement material. Moreover, GB/T 29162-2012 (coal gangue classification) considers that w($\mathbf{C a O}$+$\mathbf{M g O}$)$ < $10% belongs to aluminum-silicon gangue. 0.30$ < $m($\mathbf{A l}_{\mathbf{2}} \mathbf{O}_{\mathbf{3}}$)/m($\mathbf{S i O}_{\mathbf{2}}$)$\leq$0.5 belongs to the middle-grade silicon-aluminum ratio coal gangue. The coal gangue w($\mathbf{C a O}$+$\mathbf{M g O}$) in this paper is 1.19%, and m($\mathbf{A l}_{\mathbf{2}} \mathbf{O}_{\mathbf{3}}$)/m($\mathbf{S i O}_{\mathbf{2}}$) is 0.41, so the coal gangue belongs to the middle-level silicon-aluminum ratio coal gangue.

The primary factors influencing the mechanical grinding of CGP include the ball-to-material ratio, milling speed, and milling time. The ball-to-material ratio is defined as the mass ratio of zirconia beads to coal gangue particles. When the volume of the grinding cylinder remains constant, an increase in the pellet ratio enhances grinding efficiency and refines the particle size of the CGP. However, if the pellet ratio is too low, the CGP may become tightly adsorbed onto the grinding beads, thereby reducing the effectiveness of the grinding process. In this study, the sizes of the zirconium beads used are 10 mm, 5 mm, and 3 mm, with a mass ratio of 1:1:2. For instance, a pellet ratio of 1 indicates that 100 g of coal gangue particles and 25 g and 50 g of zirconium beads, with particle sizes of 10 mm, 5 mm, and 3 mm respectively, are added to each grinding cylinder.

The grinding procedure is as follows: First, the required quantities of zirconium beads and coal gangue particles were weighed according to the predetermined pellet ratio and sequentially poured into a clean grinding cylinder. The coal gangue particles were placed on top, and the zirconium beads at the bottom. If the zirconium beads are placed on top, they may idle or rotate ineffectively during grinding, causing the CGP to adhere to the surface of the grinding cylinder, which can negatively impact grinding efficiency and complicates sampling. After securing the pressure cover, the power supply was started, and the specific milling time and milling speed were set. Upon completion of grinding, the CGP was removed, further dried, and cooled to room temperature.

The test specimens, measuring 2mm×2mm×2mm, were prepared in a clean laboratory, with manual mixing time uniformly set to 120 seconds. After mixing, the material was poured into a non-deformable steel mold, then placed in a standard curing room (98% humidity, temperature 20±2°C) for 24 hours. Following demolding, the specimens were cured for an additional 90 days, after which their compressive strength was measured. Six test blocks were prepared for each sample group to ensure accurate measurement of compressive strength and to minimize errors.

RSM is a combination of mathematical method, statistical method, and experimental design used to explore the complex mathematical model between response output (response variables) and influencing factors (factors) of an unknown system. The design of the RSM model is mainly divided into two stages: The first stage is the factor experiment design. Three independent factors are ball-to-material ratio (1$\sim$3), milling time (0.67h$\sim$2h), and milling speed (100rpm-300rpm). The mix proportion of the experiment was designed by a central composite circumscribed (CCC) design. The second stage is to establish a suitable second-order or higher-order response surface regression model to fit the real model, explain the influence of single factors or interaction factors on the response surface variables (strength-20% and strength-50%), and guide the actual grinding process to achieve the lowest cost-effectiveness in the end. This is the core content of the RSM model establishment. At last, we obtain affordable and feasible grinding parameters for coal gangue via numerical optimization for a defined target or goal.

A well-designed factorial experiment is crucial for the successful development of the model. The commonly used factorial designs of RSM mainly include the central composite design and the Box-Behnken design. In fact, the Box-Behnken design is more appropriate for model factorial design with four or more factors [14], [15], and the CCC design is most suitable for designing three-factor experiments. Therefore, this paper selects the CCC in the central composite design, which is the most suitable for designing a three-factor experiment. Figure 3 shows the schematic diagram of the factorial design of the CCC.

In Figure 3, the corner points (blue points), axis points (yellow points), and center points (black points) of the cube represent different levels of factors. The corner points represent the 8 factors in the experimental design; the factor level of +1 (-1) corresponds to the maximum value (minimum value) of the original value range of the factor. The axis points express the six factors in the experimental design; the factor level of +1.682 (-1.682) corresponds to 1.682 times the maximum or minimum value of the factor range. The RSM improves the goodness of fit of the model by designing factors beyond the original range. The factor level of the center point is 0, which corresponds to the middle value of the original value range of the factors (usually set to 6 factors). The setting of six center points can basically avoid the problem of model accuracy caused by the center point measurement error. However, if the prediction variance of the model is too large, it is a feasible option to increase the number of center points to reduce the variance and improve the accuracy of the model. The authenticity and stability of the center data point (the number of centers in this article is 6) is the key to improving the accuracy of the response surface model. It can be seen in Table 2 that the 90-day compressive strength values of the six center points are all very close to the average value. The strength values change little.

The total number of tests means the sum of all corner points, axis points, and center points, for a total of 20 groups. The specific CCC design results are shown in Table 3. In addition, there are two response variables in the factorial design test, including strength-20% and strength-50%, which represent the 90d strength values corresponding to 20% and 50% of cement replaced by coal gangue.

The above content is about the data preparation stage before modeling. The authenticity and reliability of the data are the basis of the reliability of the model. Hence, the data acquisition stage should avoid data points that deviate from the actual value as much as possible.

After obtaining a true and reliable dataset, it is necessary to use a second-order response surface polynomial (formula such as Eq. (1)) to fit three factors and response variables. $y$ represents the response variable (strength-20% and strength-50%), $A$ means ball-to-material ratio, $B$ represents milling time, $C$ expresses milling speed. $A\, B\, C$ are the coding factors. $b$ is intercept. $a_{1 \sim 3}, b_{1 \sim 3}, c_{1 \sim 3}$ are coefficients.

Besides, the paper can model the complex functional relationship between the response variable and the coding factors by incorporating higher-order terms, thus enhancing the model's complexity. This approach can reduce prediction variance and enhance the model's applicability.

Then, the surface polynomial coefficients, model variance, and prediction model fit should be evaluated, and the prediction variance value of the model is reduced by continuously adjusting the structure of the response surface polynomial in order to improve the fitting goodness of the model. Finally, a reliable RSM prediction model is obtained, and this stage is considered to be the optimization model stage. In the final verification stage, the usability and credibility of the model are verified.

Before establishing the model, an in-depth analysis of the surface polynomial is necessary, encompassing the estimation of coefficients in terms of coded factors, the analysis of variance for the model, and fitting statistics. Evaluating the coding factor coefficients serves as a foundational step, ensuring the stability and independence of factors and thereby mitigating multicollinearity issues. Multicollinearity refers to potential high correlations between factors in the regression model, distorting regression estimations and potentially causing modeling failure. Model variance analysis assesses the significance of single and interaction factors in the model. This analysis enables the elimination of redundant elements—insignificant factors or interactions—to simplify the regression model structure. It serves as the foundation for enhancing significance (significant factors or interactions) to improve model accuracy. The fitting statistics gauge the model's fitting degree, which is essential for assessing the disparity between predicted values and real values. A higher degree of model fitting corresponds to smaller differences, indicating better predictability. Such a model demonstrates practical applicability in engineering scenarios.

The process of establishing the RSM model involves seeking a mathematical formula that explains the functional relationship between factors and response variables. Typically, previous publications have leaned towards utilizing second-order polynomial regression equations to capture the intricate functional connections between three factors and response variables [16], [17]. However, the second-order polynomial regression model tends to overlook significant ternary or even quaternary interaction factors. Although incorporating ternary or quaternary interaction factors creates a more complex structure in the modified model, it has the potential to more effectively represent the functional relationship between the three factors and response variables. The comparison and analysis presented below assess the reduced quadratic model against the modified model, aiming to derive the optimal regression model. The following section initiates an evaluation and analysis of the coding factor coefficients.

The coefficient estimate signifies the expected change in response per unit change in factor value, assuming all other factors remain constant. In an orthogonal design, the intercept denotes the overall average response derived from all experimental runs. These coefficients represent adjustments relative to this average, based on the specific factor settings. The coefficient estimate constitutes the factor's coefficient within the regression polynomial, while the intercept serves as a constant term.

The standard error and Variance Inflation Factor (VIF) serve as evaluative metrics for coefficients. The standard error reflects the degree of dispersion among the coding factor coefficients. A larger standard error indicates greater instability in the regression coefficients. Table 4 indicates that the regression coefficients derived from the polynomial regression equation, incorporating higher-order correction terms, exhibit enhanced stability.

The VIF quantifies the extent of multicollinearity within a multiple regression model, essentially testing factor correlations. In cases where factors are orthogonal, VIFs register as 1. VIFs exceeding 1 signify the presence of multicollinearity, with higher VIF values indicating more pronounced factor correlations. VIFs under 10 are considered tolerable.

Table 4 is the analysis table of the coefficient of the response variable (strength-20%). It can be seen that the modified model adds some higher-order terms. The standard error of the regression model coefficients is close to 0 and obviously smaller than the reduced quadratic model (the regression coefficient is more stable). Meanwhile, the VIF corresponding to the modified high-order surface polynomial is larger. The VIF is much smaller than 10 and belongs to the tolerable range, which shows that the multi-collinearity is less obvious. This provides theoretical reasons for the applicability of modifying high-order models. Table 5 is the coefficient analysis table for the response variable strength (50%). It can be seen that the results are similar to those in Table 4. The modified model exhibits better performance.

The analysis of variance (ANOVA), which includes the sum of squares ($SS$), degrees of freedom ($df$), mean square error ($MS$), F-value, and P-value, is used to evaluate the significance of each factor and their interactions within the model [16]. This analysis ensures the reliability of the model and optimizes its structure to the greatest extent possible.

When the F-value exceeds the critical value $F_{0.05}\left(d f, d f_e\right)$ it indicates that the fitted model is significant. If the F-value exceeds the critical value $F_{0.01}\left(d f, d f_e\right)$, it indicates that the fitted model is highly significant. Here, $df_e$ represents the degrees of freedom associated with the error term. The empirical critical value of $F$ for evaluating the reduced quadratic model is provided in the note to subtable (a) of Table 6. It is observed that the model’s F-value of 5.08 (which is greater than 4.88 but less than 10.45) indicates that the model is significant, though not highly significant. There is only a 0.70% chance that an F-value this large could occur due to noise [17].

Similarly, subtable (b) of Table 6 provides the empirical critical value of F for the modified model. The F-value of 410.64 is much greater than $F_{0.01}$(9.77), demonstrating that the modified model is highly significant and worthy of adoption.

T-tests were employed to identify nonsignificant factors in order to derive the most simplified response regression model [18]. P-values from the T-tests that are less than 0.0500 indicate that the model terms are significant, while P-values less than 0.0001 indicate that the terms are highly significant. The P-values for the modified model, including single factors and interaction factors shown in subtable (b) of Table 6, are notably smaller than those in subtable (a) of Table 6. Specifically, the terms B, C, AC, and A² are significant in the model (see subtable (a) of Table 6), whereas B, C, AC, A$^2$, C$^2$, ABC, A$^2$C, and A$^2$B$^2$ are highly significant in the modified model (see subtable (b) of Table 6). Additionally, P-values greater than 0.1000 indicate that the model terms are not significant. If a model contains many insignificant terms (excluding those required to maintain hierarchy), model reduction may enhance its performance. Consequently, the reduced quadratic model in subtable (a) of Table 6 and the modified model in subtable (b) of Table 6 have been adjusted to remove redundant terms and include necessary ones.

Similarly, an analysis of variance was conducted for the response variable strength-50%. The results of the three-factor ANOVA with strength-50% as the response variable are presented in Table 7. It is observed that the F-value for the reduced quadratic model is 2.59, which is less than the critical value $F_{0.05}(4.82)$ as shown in subtable (a) of Table 7. This indicates that there is a 7.31% chance that an F-value this large could occur due to noise. The $P-value$ (0.0731), being greater than 0.05, suggests that the model terms are not significant. Only A and A$^2$ are significant model terms. However, the presence of numerous non-significant terms (including B, C, AC, BC, B$^2$, and C$^2$) renders the reduced quadratic model insignificant, making it unsuitable for subsequent analysis.

In comparison, it is evident that the response variable strength-50% model (modified model) in subtable (b) of Table 7 demonstrates higher significance, with a P-value of less than 0.0001. Additionally, most of the P-values for the factors in this model indicate obvious significance, meaning these factors make significant contributions to the modified model. Therefore, the modified models presented in subtable (b) of Table 6 and subtable (b) of Table 7 are more suitable for describing the relationship between the response variables (strength-20% and strength-50%) and the three influencing factors (ball-to-material ratio, milling time, and milling speed).

The statistics used to assess the model's fit include the standard deviation (Std. Dev.), mean value, coefficient of variation (C.V.%), determination coefficient (R$^2$), adjusted determination coefficient (Adjusted R$^2$), and predicted determination coefficient (Predicted R$^2$). These fitting indices are crucial for evaluating the predictiveness and utility of the model, which is key to the establishment of a reliable RSM model. The Std. Dev. reflects the degree of deviation between the predicted and actual data points, with smaller values indicating less deviation. The R$^2$ value is used to evaluate the fitting degree of the regression model, ranging from 0 to 1, where a higher R$^2$ value indicates a higher degree of fit and precision of the regression equation. The Adjusted R$^2$ accounts for the effect of sample size on R$^2$, with higher values indicating better fit. Conversely, a negative Predicted R$^2$ suggests that the overall mean might be a better predictor of the response than the current model, and in such cases, higher-order models might offer improved predictive performance.

Table 8 demonstrates that the standard deviation (Std. Dev.) of the modified model is significantly smaller than that of the reduced quadratic model in the fitting analysis of the two response variables. Both R$^2$ and Adjusted R$^2$ exceed 0.995, approaching 1, which strongly indicates that the modified model possesses a higher degree of fit. Additionally, the negative Predicted R$^2$ further suggests that the reduced quadratic model is not suitable for describing the relationship between the factors and response variables, indicating the need for higher-order terms to enhance the model's suitability.

In conclusion, the evaluation of the coding factor coefficient (standard error and VIF), model analysis of variance (F-value and P-value), and model fitting (Std. Dev., R$^2$, Adjusted R$^2$, and Predicted R$^2$) clearly shows that the modified model is more appropriate for representing the real system than the reduced quadratic model. It better captures the functional relationship between the factors and the response variables. When the range of response variables approaches the optimal value, the second-order surface polynomial fitting may exhibit curvature effects, necessitating the inclusion of higher-order terms to improve the model. Furthermore, Figure 4 presents the scatter plots of the predicted versus actual values for the response variables strength-20% and strength-50% in the modified model. All data points align closely with the y=1 line, directly indicating that the predicted values are very close to the actual values.

The process of establishing the model involves several key steps. Firstly, it is crucial that the 20 groups of tests conducted during the factor design stage, which serve as the original data for the model, are both authentic and reliable. Next, a second-order surface polynomial, created using the Design-Expert 12.0 software, is employed to describe the relationship between the factors and the response variables. By incorporating necessary terms and eliminating redundant ones (as detailed in Sections 3.1.1 - 3.1.3), a more accurate and reasonable modified model is developed, which closely approximates the real system. Finally, the higher-order function formulas for strength-20% and strength-50% are presented in Eqs. (2a) and (2b).

Section 3.2.1 establishes the fitting function relationship between the three factors and the response variables, with the goal of analyzing the influence of different individual factors and their interactions on the response variables. For a more intuitive analysis, the interaction effects of two factors on strength-20% are illustrated in Figure 5, Figure 6 and Figure 7. These figures are used to explain the influence of various factors, including A: ball-to-material ratio, B: milling time, C: milling speed, as well as their interactions (such as BC, AC, AB), on the response variable (strength-20%).

The 3D surface model of subgraph (a) of Figure 5 reveals that when A=2 (factor level is 0), B and C jointly affect the response variable strength-20%. The red dot in the middle represents the space coordinate corresponding to the center point in the factor design. The response values of strength-20% are increasing with the increase in milling time or milling speed. That is, the maximum stength-20% requires the longest milling time (2h) and the maximum milling speed (300rpm). Subgraph (b) of Figure 5 explains the common influence of milling time and milling speed on the response variable strength-20% under the condition of A = 2 (the ball-to-material ratio is 2). The response variable strength-20% is plotted at 2 MPa intervals. With milling speed as the only parameter, the compressive strength values can be increased by about 8$\sim$10 MPa, while with milling time as the only parameter, the maximum strength value can be increased by about 4$\sim$6 MPa. This shows that the milling speed has more influence on the response variable strength-20% than the milling time in the range of the set factors in the range of the given independent variables.

The 3D surface model in subgraph (a) of Figure 6 shows that when B=1.335h, A (ball-to-material ratio) and C (milling speed) have a common influence on the response variable strength-20%. Subgraph (b) of Figure 6 is a diagram of the interaction between A and C under the condition B = 1.335h. It can be shown that when A is constant, the compressive strength value increases with the improvement of C. However, with C as a constant variable, the result differs from the above. When C$\leq$131rpm, the compressive strength value increases with the increase of A (less than 6 MPa). When C = 131rpm$\sim$191rpm, the compressive strength values slightly increase (less than 2 MPa). When C$\geq$191rpm, the contour line presents a semi-circular state, which indicates that the compressive strength values first increase and then decrease with the increase of A. Similarly, when the ball-to-material ratio is fixed, the compressive strength values basically increase with the improvement in milling speed. Furthermore, when A is near 2.5, the strength difference gradually decreases. An excessive ball-to-material ratio makes it hard to enhance the compressive strength of the specimen only by increasing the milling speed. It is concluded that subgraph (b) of Figure 6 intuitively shows that the milling speed has a greater impact on the response variable strength-20% than the ball-to-material ratio.

The 3D surface model in subgraph (a) of Figure 7 shows that when C=200rpm, A (ball-to-material ratio) and B (milling time) have a common influence on the response variable strength-20%. Obviously, increasing A (the only variable) will cause a trend of first increasing and then decreasing the value of strength-20%, just as the implication is expressed by the two semicircular rings in subgraph (b) of Figure 6. The C=200rpm, B$\geq$1.22h and A$\leq$2 are lower cost-effective parameters, corresponding to the strength-20% value (22MPa$\sim$24MPa). Furthermore, in the comparison of Subgraph (b) of Figure 5 and subgraph (b) of Figure 6, the equivalent strength line in subgraph (b) of Figure 7 is looser, which also explains that the milling speed is the most influential factor in limiting the growth of strength-20%. In summary, for a 20% replacement of CGP, the influence of factors on compressive strength follows the following order: milling speed is the most significant, followed by milling time, and then the ball-to-material ratio.

Figure 8, Figure 9 and Figure 10 declare the interaction impact of two factors (BC, AC, and AB) on strength-50%. The 3D surface model in subgraph (a) of Figure 8 clearly shows the growth of the response variable (strength-50%) along with the increase in milling time or milling speed. Subgraph (b) of Figure 8 explains the common influence of B and C on the response variable strength-50% under the condition of A=2 (the ball-to-material ratio is 2). The response variable strength-50% is plotted at 2MPa intervals, similar to the above subgraph (b) of Figure 7. It is found that the fixed C and the adjustment of B basically increase the corresponding variable strength-50% at five intervals (less than 10MPa). Moreover, when the B is set to a constant variable, it is not ideal to improve the strength-50% value by the increment of C. The above result shows that when the ball-to-material ratio is 2, milling time has more influence on the response variable strength-50% than milling speed for the response variable strength-50%.

The 3D surface model in subgraph (a) of Figure 9 shows the common influence of A and C on the response variable strength-50%. Subgraph (b) of Figure 9 shows the interaction between ball-to-material ratio and milling speed on the response variable strength-50% under the condition B=1.335h (milling time is 1.335h). It can be seen that the ball-to-material ratio has a greater influence on the response variable strength-50% than the milling speed. When C$\leq$151rpm, the response variable strength-50% increases with the increment of A (less than 10MPa). When C=151rpm$\sim$212rpm, the response variable strength-50% has no obvious improvement trend (less than 2MPa). The large half-ring line appears (under the condition of C$\geq$212rpm). It shows that when B=1.335h, A$\geq$2.78 and C$\geq$212rpm, the response variable strength-50% will decrease (see subgraph (b) of Figure 9). Furthermore, it seems difficult to increase the response variable strength-50% by increasing the milling speed by fixing the ball ratio.

In conclusion, the ball-to-material ratio has a greater influence on the response variable strength-50%. Moreover, when milling time is 1.355h, the milling speed is larger ($\geq$212rpm) and the ball material is larger ($\geq$2.78), the negative effect of strength should be avoided, resulting in unnecessary waste of resources.

The 3D surface model in subgraph (a) of Figure 10 shows the common influence of A and B on the response variable strength-50%. Subgraph (b) of Figure 10 exhibits an interaction impact between ball-to-material ratio and milling time on the response variable strength-50%. The effect of ball-to-material ratio is greater than that of milling time, and there are half-ring lines in both ball-to-material ratio and milling time. In the processing of CGP, the line of semi-rings (right data points) should be avoided to pursue higher economic costs in practical engineering.

With CGP as an admixture, different substitution rates have the corresponding optimal milling parameters. It is concluded that the influence of a single factor on the response variable strength-20% is as follows: milling speed $ > $ milling time $ > $ ball-to-material ratio. The influence of a single factor on the response variable strength-50% is as follows: ball-to-material ratio $ > $ milling time $ > $ milling speed. Then, in Sections 3.3-3.4, the theoretical reasons for the influence of physical properties (morphology and particle size distribution) and chemical properties (molecular structure) on the strength of CGP corresponding to different milling parameters are studied.

In Figure 11, the ball-to-material ratios of Run 0, Run 17, Run 1 and Run 15 are 0, 0.318, 2, and 3.68, respectively, under the conditions of milling time of 1.335h and milling speed of 200rpm. The most probable particle size indicates the particle size value corresponding to the highest point in the particle size curve; that is, there are the most particle size values in the particles. The most probable particle sizes of Run 0, Run 17, Run 1, and Run 15 were 724.43$\mathrm{\mu m}$, 275.42$\mathrm{\mu m}$, 181.97$\mathrm{\mu m}$, 7.85$\mathrm{\mu m}$, respectively. The range of particle size distribution was 0.72$\mathrm{\mu m}\sim$2187.76$\mathrm{\mu m}$, 0.72$\mathrm{\mu m}\sim$2187.76$\mathrm{\mu m}$, 0.63$\mathrm{\mu m}\sim$363.08$\mathrm{\mu m}$, 0.41$\mathrm{\mu m}\sim$15.14$\mathrm{\mu m}$.

The results show: (1) The particle size distribution of CGP tends to move to a small size distribution with the increase of the ball-to-material ratio (Figure 11). (2) Compared with Run 0, the most probable particle size of Run 17 (ball-to-material ratio = 0.318) decreases obviously, and the range of particle size has no change. Run 1 (ball-to-material ratio=2) has no obvious decrease in the most probable pore size, but the range of particle size decreases obviously. Run 15 (ball-to-material ratio=3.68) has an obvious decrease in the most probable pore size and the range of particle size. This shows that increasing the ratio of ball to material can effectively reduce the particle size of CGP, and the key point is to reduce all large particles (Figure 12).

In comparison with Run 0, Run 17 (ratio=0.318) has fewer large particles of coal gangue. There is no obvious difference between the CGP of Run 1 (ratio=2) and Run 17. But there are no large particles of coal gangue in Run 15, and at the same time, there are more fine particles and ultra-fine particles. The above conclusion is consistent with the conclusion of the particle size analysis. In a word, particle size analysis and morphology observation show that increasing the ball-to-material ratio can effectively reduce the particle size of CGP, and its aim is to minimize all large particles.

In Figure 13, the milling speeds of Run 0, Run 17, Run 1, and Run 15 are 31.8 rpm, 200 rpm, and 368 rpm, respectively, under the conditions of the ball-to-material ratio = 2 and milling time = 1.335 h. The most probable particle sizes of Run 0, Run 4, Run 1, and Run 10 were 724.43$\mathrm{\mu m}$, 363.07$\mathrm{\mu m}$, 181.97$\mathrm{\mu m}$, and 3.31$\mathrm{\mu m}$, respectively. The range of particle size distribution was 0.72$\mathrm{\mu m}\sim$2187.76$\mathrm{\mu m}$, 0.95$\mathrm{\mu m}\sim$2187.76$\mathrm{\mu m}$, 0.63$\mathrm{\mu m}\sim$363.08$\mathrm{\mu m}$, and 0.32$\mathrm{\mu m}\sim$275.42$\mathrm{\mu m}$.

The results show: (1) Increasing the speed causes the particle size distribution of coal gangue to move to a smaller particle size distribution, similar to the result of increasing the ball-to-material ratio (Figure 11). Also, increasing the milling speed to 200rpm can effectively reduce the particle size range (the ball-to-material ratio is 2 and the milling time is 2h), but when the milling speed increases to 368 rpm, the particle size range has little change. Most of the particle sizes in Run 10 are smaller (D60=5.46$\mathrm{\mu m}$), and there are some particles larger than 100$\mathrm{\mu m}$.

There is little difference between Run 0 and Run 4 (milling speed=31.8rpm). Compared with Run 0 and Run 4, Run 1 (milling speed=200rpm) has fewer large particles. There is no large particle gangue in Run 15 in Figure 14.

The above contents show that increasing the milling speed is more effective for easy grinding particles of coal gangue, but it is not obvious for reducing the overall particle size range (Figure 13). The volume ratio of CGP between 15.14$\mathrm{\mu m}$ and 275.42$\mathrm{\mu m}$ in the Run 10 interval accounts for about 15% of the total volume. Possibility: the milling speed is too high, resulting in some CGP agglomeration, thus restricting the further decrease of grain size [19]. Since in the process of crushing large CGP into powder, it disturbed the atomic arrangement on the surface of CGP, weakened the surface electrostatic repulsion, increased the specific surface activation energy, and increased the surface adsorption effect, which undoubtedly makes the ultra-fine CGP agglomerate again into hard aggregates under the action of electrostatic and adsorption in the grinding process. It is also possible that because the ball material ratio is small (ball material ratio=2), the mechanical friction chance of grinding steel balls on CGP is reduced, so that there are still a small number of large particles.

The effective grinding of coal gangue is promoted by optimizing the grinding parameters. It is conducive to enhancing the micro-aggregate effect of CGP in cement so as to improve the mechanical properties of cement paste. The grinding parameters in this paper are limited by the efficiency of the ball mill, and the low-speed grinding cannot effectively stimulate the activity of the minerals on the surface of CGP. Therefore, mechanical activation is not obvious in promoting cement hydration [20]. Researchers can improve the chemical activity of coal gangue by calcination [21], alkali excitation [22], and microwave activation [23]. These treatment methods are conducive to the conversion of quartz or kaolin minerals in coal gangue into active silica with low crystallinity [24] or metakaolin [25]. It collaborates with cement to carry out the hydration reaction, thus realizing the reinforcement of the internal structure of cement stone [26].

The aim of the model establishment can be explored to analyze the influence of various single factors and interaction factors on response variables, especially to guide actual engineering practice. So, this article performs multi-objective optimization to find affordable and feasible grinding parameters to meet a desirable defined goal. Multi-objective optimization was evaluated using the global desirability function Eq. (3).

Based on this equation, solutions have $D$ values, including factors, response variables, and their combinations, close to 1, indicating that the obtained solution is closer to the defined goal or target [27]. n represents the number of factors and response variables existing in the optimization, corresponding to ball-to-material ratio, milling time, milling speed, strength-20% and strength-50%. $r_i$ is the relative importance of each factor or response, and all of the relative importance are all set to most important. The weight factor is uniformly set to 1.

Table 9 exhibits the potential solutions after numerical optimization. Factor A represents the ball-to-material ratio, factor B represents milling time, and factor C represents milling speed. M-1 and M-2 are solutions based on the goal of maximization, while N-1$\sim$N-4 are solutions based on the principle of maximum impact factor. The most influential factor of the response variable strength-20% is the milling speed, and the most significant factor of the response variable strength-50% is the ball-to-material ratio. Therefore, according to the influencing factors of the response variable strength-20%, the milling speed target is maximized, the ball-to-material ratio target is minimized, and the target milling time with medium influencing factors is set to maximize and minimize, corresponding to N-1 and N-2, respectively (Table 9). Similarly, according to the influencing factors of the response variable strength (50%), the target of the ball-to-material ratio in Table 9 is maximize, and the target of the milling speed is minimize. The target settings of maximize and minimize for milling time with medium influencing factors correspond to N-3 and N-4, respectively.

In addition, in the original data of this paper, the range of strength-20% value is 17.42 MPa$\sim$30.85 MPa, and the range of strength-20% value is 7.90 MPa$\sim$14.27 MPa. So, the paper sets a defined goal or target, which should be limited to the original data: strength-20% is 25 MPa, strength-50% is 10 MPa. In the maximization target (M-1 and M-2), the targets of response variables strength-20% and strength-50% are set as: in range → 25 MPa$\sim$30.85 MPa, in range → 10 MPa$\sim$14.27 MPa, and in range → 17.42 MPa$\sim$30.85 MPa, in range → 7.90 MPa$\sim$14.27 MPa, respectively. It can be shown that the combination desirability of M-1 and M-2 is 0.977 and 1, respectively, and the D values of the factor and response variable are also greater than 0.93, which shows that the solutions obtained by M-1 and M-2 are credible. M-1 and M-2 have the same milling time (2h) and milling speed (300 rpm). However, the ball-to-material ratio of M-1 is smaller, and the corresponding strengths of 20% and 50% are 7.06% and 1.06% larger than those of M-2. The work on M-1 feels better, which is explained by the previous discussion results (the ball material ratio is not positively correlated with the strength). Besides, this means that if you want to get a defined goal or target, the lower limit of the response variable goal should be set to the defined goal or target, such as M-1.

Therefore, the setting of the response variable goal of N-1$\sim$N-4 is the same as that of M-1. The combination desirability values of N-1 and N-4 are 0.963, 0.871, 0.670, and 0.409, respectively. The combined desirability values of N-3 and N-4 are too low to be feasible solutions and should be discarded. After the above analysis, only M-1, N-1, and N-2 are more reliable optimization solutions.

In Table 10, the response variable values corresponding to the grinding parameters of M-1, N-1, and N-2 can reach the defined goal or target. The comparison results show that the strength-20% (M-1, N-1, and N-2) are all 25 MPa, and the strength-50% corresponding to M-1 is 0.25 MPa and 0.84 MPa (not more than 1 MPa) larger than N-1 and N-2, respectively. Compared to M-1, which has the largest grinding parameter value, it is more reasonable to choose N-1 and N-2. The ball-to-material ratios of N-1 and N-2 are 1.18 and 1.40, respectively, and the milling time is 2h (120s) and 0.843h (50.58s), respectively. The ball-to-material ratio is not much different. However, the milling time of N-2 is 57.85% less than that of N-1, and the milling time of each batch can be reduced by 69.42s. In other words, when the milling time is increased from 50.58s to 120s, the ball-to-material ratio is slightly reduced from 1.40 to 1.18, and the corresponding strength-50% value is only reduced by 0.59 MPa. Obviously, the N-2 is more reasonable and economical.

The quantitative functional relationship between different grinding parameters and the mechanical properties of coal gangue cement slurry has been successfully established using the RSM model in this study. By analyzing the impact of various grinding parameters, the model has identified the optimal test parameters, and its usability and applicability have been verified. This approach provides valuable insights for the rational utilization of coal gangue as a solid waste material. The key findings are as follows:

(1) A RSM model was successfully developed to describe the functional relationship between three grinding parameters and the response variables.

(2) The three factors examined were the ball-to-material ratio, milling time, and milling speed, while the two response variables were the compressive strengths measured at 90 days.

(3) The ranking of the influence factors for the response variable strength-20% was determined as follows: milling speed $ > $ milling time $ > $ ball-to-material ratio. For the response variable strength-50%, the ranking was: ball-to-material ratio $ > $ milling time $ > $ milling speed.

(4) Milling speed was found to be more effective in producing ultra-fine particles from coal gangue, although it did not significantly reduce the overall particle size range (15.14$\mathrm{\mu m}$ to 275.42$\mathrm{\mu m}$). In contrast, the ball-to-material ratio was more effective in reducing and narrowing the overall particle size range (0.4$\mathrm{\mu m}$ to 15.14$\mathrm{\mu m}$).

(5) The optimal grinding parameters were identified as N-2 (ball-to-material ratio: 1.40, milling time: 0.843 hours, milling speed: 300 rpm), which were found to be both cost-effective and feasible.

These findings underscore the significance of using RSM for optimizing the grinding parameters in coal gangue cement slurry, contributing to enhanced mechanical properties and promoting the sustainable use of coal gangue in construction materials.