## 3 Construction of MBD-DEM Coupling Model  
The interaction between the track subsystem and the particles on the road surface is shown in Figure 10. As shown in Figure 10, it is assumed that the particle and the track subsystem make contact with each other at point P. Corresponding local coordinates of the contacted particle and the track subsystem were built. It is assumed that the direction cosine matrix between the two local coordinates \(A_{k}\) could be expressed as Equation (23).  
\[A_{k}=\begin{bmatrix}\cos\varphi_{k}&-\sin\varphi_{k}\\ \sin\varphi_{k}&\cos\varphi_{k}\end{bmatrix} \tag{23}\]  
where \(\varphi_{k}\) represents the angle between the X-axis in the global coordinate system and the positive X-axis direction in the local coordinate system of the track plate.  
\begin{table}
\begin{tabular}{c c} \hline \hline
**Parameters** & **Value** \\ \hline \(M\) (kg) & 1300 \\ \(l_{1}\), \(l_{2}\) (m) & 1.1 \\ \(l_{1}\), \(l_{2}\) (kg/m) & 6.44 \\ \(k_{\mathrm{e}1}\), \(k_{\mathrm{e}2}\) & 16,000 \\ \(k_{3}\) & 48,000 \\ \(c_{\mathrm{e}1}\), \(c_{\mathrm{e}2}\) & 4000 \\ \(c_{3}\) & 5000 \\ \(\theta\) (\({}^{\circ}\)) & 170 \\ \hline \hline \end{tabular}
\end{table}
Table 7: Parameters of symmetrical suspension.  
Figure 9: The vertical accelerations of the jointed point C obtained by the MBD model and the mathematical model.  
Thus, point P in the global coordinate system (GCS) could be expressed as Equation (24) as follows:  
\[l_{k}^{p}=l_{k}+A_{k}p_{k}^{p} \tag{24}\]  
where \(l_{k}^{p}\) is the coordinate matrix of the contact point P in the global coordinate system; \(l_{k}\) is the coordinate matrix of the center of mass of a rigid body in global coordinates; and \(p_{k}^{P}\) is the coordinate matrix of point P in particle local coordinates.  
By taking the derivative of Equation (24) with respect to time, the velocity \(l_{k}^{p}\) of the contact point P can be expressed as follows:  
\[l_{k}^{\tilde{p}}=\dot{l}_{k}+IA_{k}p_{k}^{p}\varphi_{k} \tag{25}\]  
where \(\dot{l}_{k}=\big{(}\dot{x_{k}},\dot{y_{k}}\big{)}^{\text{T}}\) is the velocity of the mass center of the tracked vehicle in the GCS; \(\dot{\varphi}_{k}\) is the angular velocity of the mass center of the tracked vehicle in the GCS, \(I=\begin{pmatrix}0&-1\\ 1&0\end{pmatrix}\)[18].  
Due to the advantage of the discrete element method (DEM), which could effectively reproduce the road deformations caused by the compression of the tracked vehicle, DEM was used to build the road surface in the EDEM 2021.2 software for a comprehensive understanding of the contact mechanism between the vehicle and the road surface, and the authentic dynamic response characteristics. For improving the accuracy of the DEM model, the parameters of the particles, as shown in Table 8, were decided by experiments [19]. With these parameters of the particles, the DEM model of the road surface was built in the EDEM 2021.2 software.  
Due to the coupling interface between the Edem and Recurdyn software, the data could interact with each other. Consequently, the MDB-DEM model shown in Figure 11  
Figure 10: Schematic of contact between track and particles.  
\begin{table}
\begin{tabular}{c c} \hline \hline
**Parameters** & **Value** \\ \hline Density (kg/m\({}^{3}\)) & 2800 \\ Poisson’s ratio & 0.28 \\ Restitution coefficient & 0.5 \\ Static friction coefficient & 1.09 \\ Rolling friction coefficient & 0.05 \\ JKR Surface Energy (J\(\cdot\)m\({}^{-2}\)) & 13.05 \\ \hline \hline \end{tabular}
\end{table}
Table 8: The parameters of the particles filled in the road for the DEM model.  
was built by coupling the MDB model and the DEM model. The MBD-DEM model was solved with the parameters shown in Table 9.