## 3 Numerical Simulation
### Simulation Methods
#### 3.2.3 Thermal Coupling Algorithm  
Thermal coupling refers to a series of changes caused by the interaction of two physical fields: the temperature field and the stress field. In the numerical simulation, the damage caused by the freeze-thaw cycle is simulated using the discrete element thermal coupling method for subsequent research on the slope stability of dump sites.  
In thermodynamic calculations, the heat conduction equation for a continuum medium is as follows (assuming that strain changes are negligible in thermal calculations) [43]:  
\[-\frac{\partial q_{i}}{\partial x_{i}}+q_{v}=\rho C_{v}\frac{\partial\kappa}{ \partial t} \tag{1}\]  
\[q_{i}=-k_{ij}\frac{\partial\kappa}{\partial x_{j}} \tag{2}\]  
In this formula, q\({}_{i}\) is the heat flux; q\({}_{v}\) is the volume heat intensity; \(\rho\) is the density; C\({}_{v}\) is the specific heat; T is the temperature; and k\({}_{ij}\) is thermal conductivity.  
In the particle flow program, the model is composed of many particles, and the continuum heat conduction equation changes to the heat conduction between two particles, and its equation is as follows:  
\[-\sum_{p=1}^{N}Q_{p}+Q_{v}=mC_{v}\frac{\partial\kappa}{\partial t} \tag{3}\]  
\[Q=-\frac{\Delta\kappa}{\eta L} \tag{4}\]  
N is the number of contacts between two particles; Q\({}_{p}\) is the heat energy transferred in the p-th contact; Q\({}_{v}\) is the heat source intensity; m is the mass; C\({}_{v}\) is the specific heat; \(\Delta\)T is the temperature difference between the two particles; retains the heat value per unit length of contact; and L is the contact length.  
The indoor test results and the parameter calibration process lead to an exploration of how damage changes with the number of freeze-thaw cycles T. The changes in macroscopic strength and elastic modulus with the number of freeze-thaw cycles is shown in Figure 8.  
Based on the pre-pressed and cemented model, the thermodynamic calculation module is opened. The initial default temperature is 0 \({}^{\circ}\)C to simulate the freezing and thawing process of the lab test. For the SRM model, the freezing temperature of the surrounding particles is \(-\)30 \({}^{\circ}\)C, and the melting temperature is 25 \({}^{\circ}\)C. The process of heat conduction in the freeze-thaw cycle of the SRM model is shown in Figure 9.  
Based on the lab test results and parameter calibration process, the changing pattern of damage with the number of freeze-thaw cycles \(T\) is explored. The changes in the compressive strength and elastic modulus with the number of freeze-thaw cycles are shown in Figure 4.  
\begin{table}
\begin{tabular}{c c c c} \hline \hline
**Parametric** & **Value** & **Parametric** & **Value** \\ \hline emod & 40 MPa & pb\_ten & 300 kPa \\ kratio & 5.0 & pb\_coh & 60 kPa \\ pb\_emod & 40 MPa & fric & 0.8 \\ \hline \hline \end{tabular}
\end{table}
Table 2: Calibration results of the parallel bonding model parameters.  
According to the parameter calibration results, the peak stress has a linear relationship with the cohesion force pb_coh, and the elastic modulus has a linear relationship with the effective modulus mod. Therefore, as a result of the integrated trial calculation analysis, the effective modulus damage \(D_{\epsilon}(T)\) and cohesion damage \(D_{\epsilon}(T)\) can be obtained. The evolution function with the number of freeze-thaw cycles \(T\) is as follows:  
\[D_{\epsilon}(T)=1-\frac{4.37^{-\frac{T}{776}}-0.441}{3.928} \tag{5}\]  
\[D_{\epsilon}(T)=1-\frac{2.95-0.14T+0.01T^{2}-0.0002T^{3}}{2.95} \tag{6}\]  
By introducing the above formula into the model, we can obtain the changes in the physical and mechanical properties of the SRM with the number of freeze-thaw cycles during the entire thermodynamic calculation process and examine its microscopic changes.  
Figure 8: (a) Relationship between modulus of elasticity and number of freeze-thaw cycles. (b) Relationship between compressive strength and number of freeze-thaw cycles.  
Figure 9: Freeze–thaw cycle process of SRM model.