# Efficient Wheel-Rail Stick-Slip Numerical Modeling for Railway Traction Vehicles
## 2 Materials and Methods
### Zero-Crossings Detection  
In the slip mode, transitions to stick may occur when the relative velocity at a contact point vanishes. However, confirming this mode change requires zero relative velocity and a contact force within the friction cone (\(\left\lvert{F}\right\rvert\leq F_{sMax}\)). Finding these events necessitates a root-finding algorithm, which can be computationally expensive, as discussed in [35]. An alternative approach proposed by the Karnopp model [41] (Figure 2) defines a \"stick band\" around zero velocity. Stick occurs when the relative velocity falls within this band (\(\left\lvert{v}\right\rvert<\eta\)), and slip occurs when the relative velocity falls outside the band (\(\left\lvert{v}\right\rvert\geq\eta\)). While this introduces a discontinuity in the system, it eliminates the need for computationally expensive zero-velocity detection and allows for efficient simulations. The friction law defined in Equation (1) will change as in Equation (2):  
\[F_{f}\in \begin{cases}\quad F_{s}(t,v),\text{ if }\left\lvert{v}\right\rvert> \eta,\text{ (slip)}\\ \quad F_{a}\in\left[0\,F_{sMax}\right]\text{ if }\left\lvert{v}\right\rvert \leq\eta,\text{ (stick)}\end{cases}  \tag{2}\]  
The Karnopp model offers the following advantages:  
1. Eliminates difficulties with zero-velocity detection. 2. Enables efficient simulations. 3. Allows for distinguishing stick and slip periods even with a wide stick band. However, a drawback exists: setting the relative acceleration to zero during the stick phase can cause a \"drift-off effect\" for large intervals. This effect leads to numerical instability in the ODE solver [17]. Despite this limitation, the Karnopp model's benefits make it the method of choice for all simulations presented in this paper. Numerical instability arises during simulations due to the \"drift-off effect\" mentioned earlier. A common technique to address this is the Baumgarte stabilization method [31]. This method essentially adds corrective forces to counter drift by introducing a linear combination of acceleration, velocity, and position terms into the constraint acceleration equation. The physical interpretation is that these forces act to enforce the velocity and position constraints [42]. A detailed analysis of Baumgarte stabilization and its impact on numerical efficiency can be found in [43]. However, the Baumgarte method is not always necessary. In the case of static friction, position information is not critical for maintaining numerical stability during the stick phase. Here, the acceleration constraint depends solely on the relative velocity. To improve upon the Karnopp model, ref. [17] introduced a linear term in the static friction force that guides the relative velocity towards zero. This linear term helps eliminate potential chattering that might occur at the boundary of the zero-velocity interval. Alternative friction models have also been proposed. For example, ref. [44] suggests a model that uses the width of the stick band to smoothly transition between the kinetic and static friction values. This approach eliminates the discontinuity caused by the switch from slip to stick. Unlike previous methods, this model introduces a non-linear term derived from the continuity condition imposed on the friction force. In practice, the width of the stick band (\(\eta\)) should be significantly smaller than the typical amplitude of the velocity to avoid affecting the solution and solver tolerance. Similarly, the solver tolerance (Tol) should be much smaller than \(\eta\), following the relationship \(\max(\eta)>>\eta>>\) Tol. Setting overly strict tolerances (\(\eta\) and Tol) can significantly increase computation time without substantial benefit. The method presented in [36] has been shown to produce good results even with large \(\eta\) values, while also offering considerable improvements in numerical efficiency. It is important to note that there is no single \"one-size-fits-all\" solution for numerical stabilization in constrained systems integrated with standard ODE solvers.