## 2 Wheel-Rail Contact Model
### Wheel-Rail Contact Force Calculation  
The wheel-rail normal contact force based on the Hertz elasticity theory is given as:  
\[F_{N}=K_{nr}\delta_{c}^{3/2} \tag{6}\]  
where \(K_{nr}\) is the wheel-rail normal contact stiffness, and \(\delta_{c}\) is the wheel-rail normal compression amount. According to an empirical formula, \(K_{nr}\) is expressed as:  
\[K_{nr}=(1/G)^{3/2} \tag{7}\]  
where \(G\) is the wheel-rail contact constant \(\left(\text{m}/\text{N}^{2/3}\right)\). In contrast, the general expression for \(K_{nr}\) derived from the Hertz theory is given as [24]:  
\[K_{nr}=\frac{2\sqrt{2}}{3}\frac{q_{k}}{(\delta_{1}+\delta_{2})}\left(\frac{ \rho}{R_{nr}}\right)^{1/2}\sqrt{R_{nr}} \tag{8}\]  
where \(\delta_{1}\) and \(\delta_{2}\) are the calculated parameters related to the Poisson's ratio and elastic modulus of the wheel and rail; \(\rho\) is the parameter related to \(R_{ny}\) and \(R_{ry}\) at the contact points; \(R_{ny}\) and \(R_{ry}\) are the curvature radii of the wheel and rail profiles; and \(R_{ux}\) is the rolling circle radius of the wheel. \(\rho/R_{wx}\) is obtained from the corresponding calculation table. The intermediate calculation parameter \(q_{k}\) is expressed as:  
\[q_{k}=[m/2K(e)]^{3/2} \tag{9}\]  
where \(m\) and \(K(e)\) are the intermediate parameters in the solution obtained using the Hertz theory. To clarify the difference between the empirical formula and the theoretical formula, the cases of wheel-tread contact and wheel-flange contact were analyzed. The different colored line segments on the wheel and rail profiles represent areas with different curvature radii in Figure 4. For wheel-tread contact, the contact points on the wheel and rail are located on lines AB and CD, respectively, as shown in Figure 4. In this case, the applicability of the empirical formula has been verified [24]. In contrast, for wheel-flange contact, the contact points are located on the lines EF and GH, respectively, as shown in Figure 4. In this case, \(R_{ny}\) and \(R_{ry}\) at the contact points are significantly different from that in wheel-tread contact. When using the theoretical formula, the calculation result of \(K_{nr}\) must be significantly different from that in wheel-tread contact according to Equation (8). However, \(K_{nr}\) will not change compared to that in wheel-tread contact when using the empirical formula according to Equation (7). Hence, using the empirical formula tends to cause large calculation errors of normal contact stiffness in this case. After obtaining the normal contact force within each contact patch, the wheel-rail tangential force is first calculated using [PERSON] linear theory [29]. Considering the non-linear relationship between wheel-rail creepages and creep forces for large creepages, the  
Figure 4: Distribution range of wheel-rail contact point positions in wheel-tread contact and flange contact. Shen-Hedrick-Elkins theory [5] is applied for the nonlinear correction of the tangential contact force. Thus, the wheel-rail contact forces in each contact patch can be calculated, represented by blue line segments with arrows in Figure 5. Notably, the contact forces should be converted to an absolute coordinate system in dynamic analysis, which is given as:  
\[F_{Y}=\sum F_{Yi},\ F_{Z}=\sum F_{Zi},\ M_{X}=\sum M_{Xi},\ M_{Z}=\sum M_{Zi}(i=1,2,\ldots n) \tag{10}\]  
where \(F_{Y}\), \(F_{Z}\), \(M_{X}\) and \(M_{Z}\) are composed by the wheel rail contact force and spin moment; \(i\) is the indicator of each wheel-rail contact point, and \(n\) is the number of contact points; and \(F_{Yi}\), \(F_{Zi}\), \(M_{Xi}\) and \(M_{Zi}\) can be expressed as:  
\[ \begin{array}{l}F_{Xi}=F_{Ni}\cdot\mathbf{i}+F_{xi}\cdot\mathbf{i}+F_{yi}\cdot\mathbf{i},F_{Yi}=F_{Ni}\cdot\mathbf{j}+F_{xi}\cdot\mathbf{j}+F_{yi}\cdot\mathbf{j}\\ F_{Zi}=\mathbf{N}_{i}\cdot\mathbf{k}+F_{xi}\cdot\mathbf{k}+F_{yi}\cdot\mathbf{k},M_{Xi}=M_{i} \cdot\mathbf{i},M_{Zi}=M_{i}\cdot\mathbf{k}\end{array}