Physics-informed neural network for friction-involved nonsmooth dynamics
  problems
Friction-induced vibration (FIV) is very common in engineering areas.
Analysing the dynamic behaviour of systems containing a multiple-contact point
frictional interface is an important topic. However, accurately simulating
nonsmooth/discontinuous dynamic behaviour due to friction is challenging. This
paper presents a new physics-informed neural network approach for solving
nonsmooth friction-induced vibration or friction-involved vibration problems.
Compared with schemes of the conventional time-stepping methodology, in this
new computational framework, the theoretical formulations of nonsmooth
multibody dynamics are transformed and embedded in the training process of the
neural network. Major findings include that the new framework not only can
perform accurate simulation of nonsmooth dynamic behaviour, but also eliminate
the need for extremely small time steps typically associated with the
conventional time-stepping methodology for multibody systems, thus saving much
computation work while maintaining high accuracy. Specifically, four kinds of
high-accuracy PINN-based methods are proposed: (1) single PINN; (2) dual PINN;
(3) advanced single PINN; (4) advanced dual PINN. Two typical dynamics problems
with nonsmooth contact are tested: one is a 1-dimensional contact problem with
stick-slip, and the other is a 2-dimensional contact problem considering
separation-reattachment and stick-slip oscillation. Both single and dual PINN
methods show their advantages in dealing with the 1-dimensional stick-slip
problem, which outperforms conventional methods across friction models that are
difficult to simulate by the conventional time-stepping method. For the
2-dimensional problem, the capability of the advanced single and advanced dual
PINN on accuracy improvement is shown, and they provide good results even in
the cases when conventional methods fail.
Zilin Li, Jinshuai Bai, Huajing Ouyang, Saulo Martelli, Jun Zhao, Ming Tang, Yang Yang, Hongtao Wei, Pan Liu, Wei-Ron Han, Yuantong Gu
# Physics-informed neural network for friction-involved nonsmooth dynamics problems

###### Abstract

Nonsmooth, discontinuous dynamic behaviour of systems containing multiple contact points is an open research problem. Accurately simulating the constantly nonsmooth/discontinuous dynamic behaviour is challenging. This study presents novel physics-informed neural network (PINN) frameworks for solving nonsmooth friction-induced, or friction-involved, vibration (FIV) problems where the training of the neural network training is guided by the relevant physics laws. Four kinds of high-accuracy PINN-based methods are proposed: (1) single PINN; (2) dual PINN; (3) advanced single PINN; (4) advanced dual PINN. The single PINN method uses the nonsmooth dynamic formulation for multibody systems, as physical constraints, substituting the time iteration loop in the conventional dynamic analysis. For the dual PINN method, in conjunction with the single PINN, a new PINN is proposed as an alternative to traditional numerical methods for solving linear complementary problems (LCP). The advanced single/dual PINN methods improve the corresponding single/dual PINN methods by incorporating an interpolation technique during the training process. The success of the proposed PINN-based methods in predicting the nonlinear vibration of two typical nonsmooth dynamics problems is verified: one is a 1-dimensional system with direct rigid contact considering stick-slip oscillation, and the other is a 2-dimensional system with spring contact considering separation-reattachment and stick-slip oscillation. Overall, the
developed PINN-based methods are accurate in predicting the nonsmooth dynamic behaviour, which does not require a tiny time step as the conventional time-stepping methods.

Key words :Physics-informed neural network, Nonsmooth dynamics, Stick-slip, Contact loss, Friction-induce vibration;

## 1 Introduction

Friction is ubiquitous in daily life as well as in engineering. For centuries, humans have exploited friction and attempted to control friction. When friction is involved in a dynamic environment, the dissipating energy property of friction is generally used to reduce motion or vibration. However, in certain conditions, the stability of a system counterintuitively degrades due to friction [1-3]. Friction-induced vibration (FIV) [4-6] is one kind of self-excited vibration that associates with a variety of engineering issues, including squealing brakes[4], squeaky joints [7], chattering of cutting tools [8], inaccurate positioning of the robotics [9], etc. The dynamics of the systems with friction are protean and complex. There is a lack of comprehensive understanding of the dynamic characteristics when friction is involved. The manifestations of nonsmoothness that could be directly attributed to friction are: (1) the direction-changing feature of friction force with respect to the relative velocity; (2) stick-slip vibration, which is known as one of the main mechanisms of FIV [10]; (3) separation-reattachment events [7] resulting from unstable FIV. As the contact states change during vibration in a dynamic environment, even systems that are linear within each motion regime (e.g., stick or slip) become nonlinear and time-varying.

Research works on FIV choose to focus on two kinds of systems: (1) low-degrees-of-freedom theoretical models and (2) complicated structures. As for the former kind, models are usually nonlinear and often even nonsmooth. Stick-slip is an essential mechanism of friction-induced vibration and exists in many applications [11, 12]. Popp and Stelter [13] pointed out that even in a one-degree-of-freedom slider-belt model, dry friction-induced stick-slip vibration could be periodic, quasi-periodic and even chaotic. Leine et al. [14] proposed a switching method that is conducive for stick-slip analysis of the slider-belt model with a single contact point. With the development of unstable vibration, loss of contact may occur. Li et al. [15] used a slider-belt system with an inclined spring to reveal that separation and reattachment events could happen when unstable vibration was aroused by mode-coupling instability. In their follow-up work [16], a slider moving on an elastic disc also exhibited stick-slip vibration and the separation and reattachment phenomenon. The
conventional approach for determining FIV of a theoretical model that considers nonsmooth features such as stick-slip or separation-reattachment involves tracking the motion states (e.g., stick or slip, contact or separation), which is tedious and time-consuming. For example, Pascal [17] studied the stick-slip motion of a two-degrees-of-freedom lumped-mass-belt model when two sliders were involved, and the stick-slip transition was monitored and captured to ensure the equations of motion with respect to the right motion states were used.

Detailed finite element models of complicated structures, such as automotive brakes, are often linearised [18] for stability analysis. Such a finite element model can have many contact points. Chen et al. [19] found that in a finite element model with multiple contact points, the separation took place between various contact nodes in different regions of the contact interface. However, stick-slip oscillation could not be included in their simulation as modelling strategies suitable for problems with very few contact points are no longer valid for many contact points [14]. Dealing with the contact and friction for large contact areas is difficult [20, 21, 22, 23]. There are extensive research concerning the discontinouty in multi-body systems, such as robotics or rotor systems [24]. Moreau [25] proposed the differential measurement and scanning process theories as foundations for nonsmooth dynamic analysis. The unilateral contact problem of multibody systems can be transformed into a linear complementary problem (LCP). Several methods based on the LCP formulations for systems with rigid or elastic unilateral contact have been developed [24, 26, 27]. Two general approaches in multibody dynamics are the event-driven method and the time-stepping method [21, 22, 23]. Li et al. [28] found that for the friction-induced vibration problem with multiple contact points, the system displays repeated separation-reattachment and stick-slip events. The complex nonsmooth dynamics are hard to capture unless a tiny time step is used.

In other reseachlines, several novel approaches were introduced. One such approach is the moving least squares material point method (MLS-MPM) developed by Hu et al. [29], which considers displacement discontinuity and frictionless boundaries. Based on this method, a real-time simulator for soft robotics was also developed [3]. Another approach was proposed by Peng et al. [30], who suggested using a method based on the symplectic discrete format to tackle nonsmooth problems with friction and impact. This approach enabled accurate simulations with a larger time step. Fazeli et al. [31], on the other hand, developed a hierarchical learning approach to teach manipulation strategies that consider friction. The development of simulators that possess accuracy, efficiency, and robustness remains a significant area of interest.
Physics-informed neural network (PINN), which is proposed by Raissi et al. [32] in 2019, is considered a promising method for finding the solution of the partial differential equations (PDEs) of complex systems. PDEs, initial and boundary conditions can be applied to guide the training of deep learning neural networks. Karniadakis et al. [33] and Bai et al. [34] reviewed the development of PINN in various fields. It has been demonstrated that, combining the advantages of neural networks and physics laws, PINN is effective for dealing with problems with insufficient data and thus has the prospect of a wide range of applications. Samaniego [35] presented an energy-based PINN framework, which initially proofed the effectiveness of PINN for mechanics applications in terms of hyperelasticity, crack propagation and piezoelectric beams and plates problems. Gu's team [36, 37] developed a series of computing and modelling methods based on PINN for hydrodynamics and topology optimisation. Pfrommer et al. [38] solved the contact-induced dynamic discontinuous problem of a rigid body by the neural network. In that work, a comprehensive loss function involving governing equation, unilateral constraint law and maximum energy dissipation law was defined in their work. The new method could represent realistic impact and stick motion even when spare training data is used.

Accurately predicting nonsmooth dynamic phenomena from contact is a fundamental problem in many applications. The primary goal of this research is to propose a new methodology for the solution in the time domain of nonsmooth dynamics problems by combining neural network and the theories of nonsmooth dynamics into physics-informed neural network frameworks. The paper is organised as: (1) Section 2 presents the mix-level time-stepping modelling method; (2) Section 3 introduce the frameworks of PINN for solving LCP, and the single and dual PINN for the dynamic simulation of the complex systems with friction and contact loss; (3) Section 4 illustrates the application of proposed single PINN and dual PINN frameworks in the FIV problem with stick-slip. Eventually, the accuracy of the single/dual PINN framework in solving the FIV problem with separation-reattachment and stick-slip is determined and presented.

## 2 Theoretical Approach and Numerical Implementation

### The linear complementary problem in unilateral contact

Mathematically, unilateral contact can be described by linear complementary conditions. Specifically, Fig.1(a) shows the complementary relations for single-point rigid contact in the normal
direction. The distance between the object of interest and the other contact surface is denoted by \(g_{\rm Ni}\), and the contact force by \(\lambda_{\rm Ni}\).

These two quantities satisfy the relationships that if \(g_{\rm Ni}\)\(\geq\)0, \(\lambda_{\rm Ni}\)=0; otherwise if \(g_{\rm N}\)=0, \(\lambda_{\rm Ni}\)\(\geq\)0. For systems with \(m\) contact points, \(g_{\rm Ni}\) (\(i\)=1, 2,..., \(m\)) is the \(i\)th element of vector \(\mathbf{g}_{\rm N}\) and \(\lambda_{\rm Ni}\) (i=1, 2,..., \(m\)) is the \(i\)th element of vector \(\mathbf{\lambda}_{\rm N}\). The complementary relationship is given in the vector form:

\[0\mbox{$\leq$}\underline{g}_{\rm N}\,\mbox{$\perp$}\,\mbox{$\lambda_{\rm N}$} \mbox{$\geq$}0 \tag{1}\]

Additionally, for spring contact shown in Fig.1(b), \(\Omega_{\rm Ni}\)=\(k_{\rm c}g_{\rm Ni}\)+\(\lambda_{\rm Ni}\) is introduced. The following revised relationship now applies: if \(\lambda_{\rm Ni}\)=0, \(\Omega_{\rm Ni}\)\(\geq\)0, otherwise if \(\lambda_{\rm Ni}\)\(\geq\)0, \(\Omega_{\rm Ni}\)=0. For multiple spring contact, \(\Omega_{\rm Ni}\) (\(i\)=1, 2,..., \(m\)) is the \(i\)th element of vector \(\mathbf{\Omega}_{\rm N}\). The complementary relationship is written as:

\[0\mbox{$\leq$}\mathbf{\Omega}_{\rm N}\,\mbox{$\perp$}\,\mbox{$\lambda_{\rm N }$}\mbox{$\geq$}0 \tag{2}\]

When friction is considered in the tangential direction, a general friction law can be described by a function of friction force, namely \(\lambda_{\rm Ti}\), with respect to the relative velocity \(\gamma_{\rm Ti}\), which is shown in Figure 1 (b). It is generally assumed that sliding friction \(\lambda_{\rm Ti}\) follows Coulomb's law of friction, which is expressed as \(\lambda_{\rm Ti}\)=\(\mu\)\(\lambda_{\rm Ni}\). Like the notation rules stated in the normal contact, \(\lambda_{\rm Ti}\) and \(\gamma_{\rm Ti}\) (\(i\)=1, 2,..., \(m\)) are the elements of vector \(\lambda_{\rm T}\) and \(\gamma_{\rm T}\) for describing the friction and relative velocity of the multi-point contact problem. Through mathematical transformation, the friction law can be decomposed into two complementary relations:

\[0\mbox{$\leq$}\,\mbox{$\lambda_{\rm R}^{\top}$}\,\mbox{$\perp$}\,\mbox{$\gamma _{\rm R}$}\mbox{$\geq$}\,0\,,\,\,\,\,\,0\mbox{$\leq$}\,\mbox{$\lambda_{\rm L}^{ \top}$}\,\mbox{$\perp$}\,\mbox{$\gamma_{\rm L}$}\mbox{$\geq$}\,0 \tag{3}\]

where \(\gamma_{\rm R}\) = \(\frac{1}{2}\)(\(\mid\)\(\gamma_{\rm T}\)\(\mid\)+\(\gamma_{\rm T}\)), \(\gamma_{\rm L}\) = \(\frac{1}{2}\)(\(\mid\)\(\gamma_{\rm T}\)\(\mid\)-\(\gamma_{\rm T}\)), \(\lambda_{\rm R}\)=\(\mu\)\(\lambda_{\rm N}\)+\(\lambda_{\rm T}\) and \(\lambda_{\rm L}\)=\(\mu\)\(\lambda_{\rm N}\) - \(\lambda_{\rm T}\),

with the relationship

\[\mbox{$\lambda_{\rm R}$}\mbox{$\geq$}\,\mbox{$
\[\gamma_{\tau}=-\gamma_{\rm R}+\gamma_{\rm L}\ \ \ {\rm and}\ \ \ \lambda_{\rm R}=2\mu\dot{\lambda}_{\rm N}-\dot{\lambda}_{\rm L} \tag{4}\]

or between the velocity and impulse of the friction force:

\[0\leq\Lambda_{\rm R}^{\tau}\perp\gamma_{\rm R}\geq 0\,,\ \ \ 0\leq\Lambda_{\rm L}^{ \tau}\perp\gamma_{\rm L}\geq 0 \tag{5}\]

where \(\Lambda_{\rm R}=\int\lambda_{\rm R}{\rm d}t\) and \(\Lambda_{\rm L}=\int\lambda_{\rm L}{\rm d}t\)

### Equation of motion of the multibody system with unilateral contact

For the dynamic modelling of the multibody systems with unilateral contact, the equation of motion in matrix form is generally written as Eq.(6) :

\[{\bf M}\ddot{\bf q}-{\bf h}-{\bf W}_{\rm N}\dot{\lambda}_{\rm N}-{\bf W}_{\rm T }\dot{\lambda}_{\rm T}={\bf 0} \tag{6}\]

in which \({\bf M}\) is the mass matrix, \({\bf h}\) is the force vector for the system regardless of the unilateral conditions, \(\dot{\lambda}_{\rm N}\) and \(\dot{\lambda}_{\rm T}\), are the normal and tangential force vectors, respectively. \({\bf W}_{\rm N}\) is the displacement vector defined in the global coordinates. The constraint functions in the normal direction and tangential direction are expressed as \({\bf g}_{\rm N}={\bf W}_{\rm N}^{\tau}{\bf q}\), and \({\bf g}_{\rm T}={\bf W}_{\rm T}^{\tau}{\bf q}\).\({\bf W}_{\rm N}\) and \({\bf W}_{\rm T}\) can be obtained by \({\bf W}_{\rm N}=\left(\frac{\partial\xi_{\rm N}}{\partial{\bf q}}\right)^{\rm T}\) and \({\bf W}_{\rm T}=\left(\frac{\partial\gamma_{\rm T}}{\partial{\bf u}}\right)^{ \rm T}\).

Then, take the first derivative of \({\bf g}_{\rm N}\) and \(\gamma_{\rm T}\) with respect to time, the normal velocity and tangential acceleration can be expressed as:

\[\gamma_{\rm N}={\bf W}_{\rm N}^{\tau}\dot{\bf q}+\bar{\omega}_{\rm N}\,,\ \ \dot{\gamma}_{\rm T}={\bf W}_{\rm T}^{\tau}\ddot{\bf q}+\frac{\partial{\bf W}_{ \rm T}^{\tau}}{\partial t}\dot{\bf q}+\dot{\bar{\omega}}_{\rm T} \tag{7}\]

\[{\rm in\ which}\ \ \ \bar{\omega}_{\rm N}=\frac{\partial{\bf W}_{\rm N}^{ \tau}}{\partial t}\ \ \ {\rm and}\ \ \ \dot{\bar{\omega}}_{\rm T}=\frac{\partial^{2}{\bf W}_{\rm T}^{ \tau}}{\partial t^{2}}\.\]

By discretising the motion states over time interval \(\Delta t=t_{\rm E}-t_{\rm A}\), one obtained the time-domain discretised equation of motion and the state vector:

\[\begin{cases}{\bf M}^{\rm(A)}\Delta{\bf u}-{\bf h}^{\rm(A)}\Delta t-{\bf W}_{ \rm N}^{\rm(A)}\Lambda_{\rm N}-{\bf W}_{\rm T}^{\rm(A)}\Lambda_{\rm T}={\bf 0} \\ \Delta{\bf q}=({\bf u}^{\rm(A)}+\Delta{\bf u})\Delta t\end{cases} \tag{8}\]

in which \(\Lambda_{\rm N}=\lambda_{\rm N}\Delta t\), \(\Lambda_{\rm T}=\lambda_{\rm T}\Delta t\) and superscript (A) represents the corresponding physical quantity at instant \(t_{\rm A}\).
The increment of constraint function \(\Delta\mathbf{g}_{{}_{\mathrm{N}}}\) and \(\Delta\mathbf{\gamma}_{{}_{\mathrm{T}}}\) over \(\Delta t\) can be obtained according to Eq.(7):

\[\Delta\mathbf{g}_{{}_{\mathrm{N}}}=\mathbf{W}_{{}_{\mathrm{N}}}^{{\mathrm{(AT)} }}\Delta\mathbf{q}+\hat{\mathbf{o}}_{{}_{\mathrm{N}}}\Delta t\,\ \ \Delta\mathbf{\gamma}_{{}_{\mathrm{T}}}=\mathbf{W}_{{}_{\mathrm{T}}}^{{ \mathrm{T}}}\Delta\mathbf{u}+\frac{\hat{\mathbf{o}}\mathbf{W}_{{}_{\mathrm{T}} }^{{\mathrm{T}}}}{\hat{\mathbf{o}}\mathbf{q}}\,\Delta\mathbf{q}+\hat{\hat{ \mathbf{o}}}_{{}_{\mathrm{T}}}\Delta t \tag{9}\]

In corporation with Eq.(8), motion states and the constraint function at the next time instance \(t_{\mathrm{E}}\) are derived:

\[\mathbf{g}_{{}_{\mathrm{N}}}^{{\mathrm{(E)}}}=\mathbf{G}_{{}_{\mathrm{NN}}} \mathbf{A}_{{}_{\mathrm{N}}}^{{\mathrm{(E)}}}\Delta t+\mathbf{G}_{{}_{\mathrm{NT }}}\mathbf{A}_{{}_{\mathrm{T}}}^{{\mathrm{(E)}}}\Delta t+\mathbf{g}_{{}_{ \mathrm{N}}}^{{\mathrm{(A)}}}+\mathbf{W}_{{}_{\mathrm{N}}}^{{\mathrm{(AT)}}} \mathbf{u}^{{\mathrm{(A)}}}\Delta t+\mathbf{G}_{{}_{\mathrm{N}}}\Delta t+\hat{ \mathbf{o}}_{{}_{\mathrm{N}}}\Delta t \tag{10}\]

\[\mathbf{\gamma}_{{}_{\mathrm{T}}}^{{\mathrm{(E)}}}=\mathbf{G}_{{}_{\mathrm{ TN}}}\mathbf{A}_{{}_{\mathrm{N}}}^{{\mathrm{(E)}}}+\mathbf{G}_{{}_{\mathrm{TT}}} \mathbf{A}_{{}_{\mathrm{T}}}^{{\mathrm{(E)}}}+\mathbf{G}_{{}_{\mathrm{T}}} \mathbf{h}^{{\mathrm{(A)}}}\Delta t+\frac{\hat{\mathbf{o}}\mathbf{W}_{{}_{ \mathrm{T}}}^{{\mathrm{(AT)}}}}{\hat{\mathbf{o}}\mathbf{q}}\,\mathbf{u}^{{ \mathrm{(A)}}}\Delta t+\hat{\hat{\mathbf{o}}}_{{}_{\mathrm{T}}}\Delta t \tag{11}\]

in which \(\mathbf{G}_{{}_{\mathrm{N}}}=\mathbf{W}_{{}_{\mathrm{N}}}^{{\mathrm{(A)}}} \mathbf{M}^{{}^{\mathrm{-1(A)}}}\), \(\mathbf{G}_{{}_{\mathrm{NN}}}=\mathbf{W}_{{}_{\mathrm{N}}}^{{\mathrm{(A)}}} \mathbf{M}^{{}^{\mathrm{-1(A)}}}\mathbf{W}_{{}_{\mathrm{N}}}^{{\mathrm{(A)}}}\), \(\mathbf{G}_{{}_{\mathrm{NT}}}=\mathbf{W}_{{}_{\mathrm{N}}}^{{\mathrm{(A)}}} \mathbf{M}^{{}^{\mathrm{-1(A)}}}\mathbf{W}_{{}_{\mathrm{T}}}^{{\mathrm{(A)}}}\), \(\mathbf{G}_{{}_{\mathrm{T}}}=\mathbf{W}_{{}_{\mathrm{N}}}^{{\mathrm{T(A)}}} \mathbf{M}^{{}^{\mathrm{-1(A)}}}\mathbf{W}_{{}_{\mathrm{T}}}^{{\mathrm{(A)}}}\),

\[\mathbf{G}_{{}_{\mathrm{T}}}=\mathbf{W}_{{}_{\mathrm{T}}}^{{\mathrm{(A)}}} \mathbf{M}^{{}^{\mathrm{-1(A)}}}+\frac{\hat{\mathbf{o}}\mathbf{W}_{{}_{\mathrm{ T}}}^{{\mathrm{T}}}}{\hat{\mathbf{o}}\mathbf{q}}\,\mathbf{M}^{{}^{\mathrm{-1(A)}}} \Delta t\,\ \ \mathbf{G}_{{}_{\mathrm{TN}}}=\mathbf{G}_{{}_{\mathrm{T}}}\mathbf{W}_{{}_{ \mathrm{N}}}^{{\mathrm{T(A)}}}\,\ \ \mathbf{G}_{{}_{\mathrm{TT}}}=\mathbf{G}_{{}_{\mathrm{T}}} \mathbf{W}_{{}_{\mathrm{T}}}^{{\mathrm{T(A)}}}\,.\]

There are six unknown vectors in the four equations in Eqs. (8), (10) and (11). To solve the problem, another two complementary equations need to be introduced. In the following, the theoretical formulations for two contact conditions are presented.

#### 2.2.1 Time-stepping formulations for the unilateral rigid contact at displacement-velocity level

Unilateral rigid contact can be enforced using a complementary relation between \(\mathbf{g}_{\mathrm{N}}\) and \(\lambda_{\mathrm{N}}\). The friction is independent of the state solutions in the normal direction, which follow the complementary relation given in Eq.(1). According to Eq.(4), \(\mathbf{\gamma}_{{}_{\mathrm{T}}}^{{\mathrm{(E)}}}\) and \(\mathbf{\Lambda}_{{}_{\mathrm{T}}}^{{\mathrm{(E)}}}\) in Eq.(11) can be replaced by \(\mathbf{\gamma}_{{}_{\mathrm{L}}}^{{\mathrm{(E)}}}\) and \(\mathbf{\Lambda}_{{}_{\mathrm{R}}}^{{\mathrm{(E)}}}\). And then, by introducing \(\mathbf{\Lambda}_{{}_{\mathrm{R}}}{=}2\mathbf{\mu}\mathbf{\Lambda}_{{}_{ \mathrm{N}}}-\mathbf{\Lambda}_{{}_{\mathrm{L}}}\) according to Eq.(4), the LCP equations in the matrix form can be obtained:

\[\underbrace{\left[\begin{array}{cc}\mathbf{g}_{{}_{\mathrm{N}}}^{{\mathrm{(E) }}}\\ \mathbf{\gamma}_{{}_{\mathrm{L}}}^{{\mathrm{(E)}}}\Delta t\\ \mathbf{\Lambda}_{{}_{\mathrm{R}}}^{{\mathrm{(E)}}}\Delta t\end{array}\right]}_{ \mathbf{\gamma}}\underbrace{\left[\begin{array}{cc}\mathbf{G}_{{}_{\mathrm{ NN}}}+\mathbf{G}_{{}_{\mathrm{NT}}}\mathbf{H}&-\mathbf{G}_{{}_{\mathrm{NT}}}&\mathbf{0}\\ -\mathbf{(G}_{{}_{\mathrm{TN}}}+\mathbf{G}_{{}_{\mathrm{TT}}}\mathbf{H})&\mathbf{G }_{{}_{\mathrm{TT}}}&\mathbf{I}\\ 2\mathbf{\mu}&-\mathbf{I}&\mathbf{0}\end{array}\right]}_{\mathbf{\Lambda}} \underbrace{\left[\begin{array}{cc}\mathbf{G}_{{}_{\mathrm{N}}}^{{\mathrm{(A)}}} \Delta t\\ \mathbf{\Lambda}_{{}_{\mathrm{R}}}^{{\mathrm{(E)}}}\Delta t\\ \mathbf{\gamma}_{{}_{\mathrm{R}}}^{{\mathrm{(E)}}}\Delta t\end{array}\right]}_{ \mathbf{\gamma}}+\underbrace{\left[\begin{array}{cc}\mathbf{g}_{{}_{\mathrm{N}}}^{{ \mathrm{(A)}}}+\mathbf{W}_{{}_{\mathrm{N}}}^{{\mathrm{(AT)}}}\Delta t\mathbf{u}^{ {\mathrm{(A)}}}+\mathbf{G}_{{}_{\mathrm{N}}}\mathbf{h}^{{\mathrm{(A)}}} \Delta t+\hat{\mathbf{o}}_{{}_{\mathrm{N}}}\Delta t\\ \mathbf{\Lambda}_{{}_{\mathrm{R}}}^{{\mathrm{(E)}}}\Delta t\end{array}\right]}_{ \mathbf{\Lambda}} \tag{12}\]

with

\[0\leq\mathbf{g}_{{}_{\mathrm{N}}}^{{\mathrm{(E)}}}\perp\mathbf{\Lambda}_{{}_{ \mathrm{N}}}^{{\mathrm{(E)}}}\geq 0,\ 0\leq\mathbf{\gamma}_{{}_{\mathrm{L}}}^{{\mathrm{(E)}}}\perp\mathbf{\Lambda}_{{}_{ \mathrm{L}}}^{{\mathrm{(E)}}}\geq 0\ \text{and}\ 0\leq\mathbf{\gamma}_{{}_{\mathrm{R}}}^{{\mathrm{(E)}}}\perp\mathbf{\Lambda}_{{}_{ \mathrm{R}}}^{{\mathrm{(E)}}}\geq 0\]
.2.2 Time-stepping formulations for the unilateral spring contact at the displacement-velocity level

In contrast with direct contact, the complementary condition of the spring contact is between \(\mathbf{\Omega}_{\text{N}}\) and \(\mathbf{\dot{\omega}}_{\text{N}}\). Meanwhile, the friction force depends on the state solutions in the normal direction when contact is maintained. By substituting Eq. (9) into \(\mathbf{\Omega}_{\text{N}}\), one gets:

\[\mathbf{\Omega}_{\text{\tiny NE}}=\left(\mathbf{I}\,/\,\Delta t^{2}+\mathbf{k} _{\text{c}}\mathbf{G}_{\text{\tiny NN}}+\mathbf{k}_{\text{c}}\mathbf{G}_{\text {\tiny NT}}\mathbf{\mu}\right)\mathbf{\Lambda}_{\text{N}}\Delta t-\mathbf{k}_{ \text{c}}\mathbf{G}_{\text{\tiny NT}}\mathbf{\Lambda}_{\text{L}}\Delta t+ \mathbf{k}_{\text{c}}\left(\mathbf{g}_{\text{\tiny NA}}+\mathbf{W}_{\text{ \tiny NA}}^{\text{T}}\Delta\mathbf{\mu}_{\text{A}}+\mathbf{G}_{\text{\tiny N}} \mathbf{h}_{\text{\tiny A}}\Delta t+\widehat{\mathbf{\omega}}_{\text{\tiny N}} \Delta t\right) \tag{13}\]

Collecting Eqs. (10), (13) and \(\mathbf{\Lambda}_{\text{\tiny R}}\!\simeq\!2\mathbf{\mu}\mathbf{\Lambda}_{\text {\tiny N}}\!-\!\mathbf{\Lambda}_{\text{\tiny L}}\) in matrix form gives the LCP formulations of the whole system:

\[\underbrace{\begin{bmatrix}\mathbf{\Omega}_{\text{\tiny N}}^{\text{(E)}}\\ \mathbf{\gamma}_{\text{\tiny L}}^{\text{(E)}}\Delta t\end{bmatrix}}_{\text{ \tiny Y}}=\underbrace{\begin{bmatrix}\mathbf{I}\,/\,\Delta t^{2}+\mathbf{k}_{ \text{c}}\mathbf{G}_{\text{\tiny NN}}+\mathbf{k}_{\text{c}}\mathbf{G}_{\text {\tiny NT}}\mathbf{\mu}&-\mathbf{k}_{\text{c}}\mathbf{G}_{\text{\tiny NT}}& \mathbf{0}\\ -\left(\mathbf{G}_{\text{\tiny TN}}+\mathbf{G}_{\text{\tiny TT}}\mathbf{\mu} \right)&\mathbf{G}_{\text{\tiny TT}}&\mathbf{I}\\ 2\mathbf{\mu}&-\mathbf{I}&\mathbf{0}\end{bmatrix}\underbrace{\begin{bmatrix} \mathbf{\Lambda}_{\text{\tiny N}}^{\text{(E)}}\Delta t\\ \mathbf{\Lambda}_{\text{\tiny N}}^{\text{(E)}}\Delta t\end{bmatrix}}_{\text{ \tiny x}}}+\underbrace{\begin{bmatrix}\mathbf{k}_{\text{c}}\left(\mathbf{g}_{ \text{\tiny N}}^{\text{(A)}}+\mathbf{W}_{\text{\tiny N}}^{\text{(AT)}} \mathbf{u}^{\text{(A)}}\Delta t+\mathbf{G}_{\text{\tiny N}}\mathbf{h}^{\text{ (A)}}\Delta t+\widehat{\mathbf{\omega}}_{\text{\tiny N}}\Delta t\right)\\ -\mathbf{\gamma}_{\text{\tiny T}}^{\text{(A)}}\Delta t-\left(\mathbf{G}_{\text{ \tiny T}}\mathbf{h}^{\text{(A)}}+\frac{\partial\mathbf{W}_{\text{\tiny T}}^{ \text{T}}}{\partial\mathbf{q}}\mathbf{u}^{\text{(A)}}+\dot{\mathbf{\omega}}_{ \text{\tiny T}}\right)\Delta t^{2}\\ \mathbf{0}\end{bmatrix}}_{\text{\tiny R}} \tag{14}\]

Eventually, Eqs. (8) and (14) together form the equations for the nonsmooth dynamic analysis of the complex frictional system with multi-contact points. The mix-level time-stepping method is called the conventional method in the following part.

### Physics-informed neural network for nonsmooth dynamics

Compared to the event-detection strategy, the time-stepping strategy is more suitable for problems with a considerable amount of contact points as it avoids the repeated detection of the state transition. However, as precision cannot be compromised, a tiny time step is required for solving complex nonlinear vibrations at the cost of a large computation time.

To achieve high-accuracy calculation of nonsmooth dynamics we propose a new approach by utilizing the physics-informed neural network (PINN). A PINN consists of the feedforward neural network (FNN) and a physics-informed loss function [34], as shown in Figure 2. In a PINN, the FNN
is used to capture the relationship between the input data and the output data, while the physics-informed loss function embeds the physics laws and quantifies the performance of the FNN [39].

As observed from Figure 2, an FNN comprises three parts, namely the input layer, the hidden layers and the output layer. When using an FNN, input data are fed into the input layer and transported forward to the next adjacent layer. The final predictions of an FNN are output from the output layer. An L-layer FNN can be mathematically expressed as

\[{}^{0}\textbf{a}={}^{0}\textbf{a}, \tag{15}\] \[{}^{l}\textbf{a}=\sigma\big{(}{}^{l-1}\textbf{w}\cdot{}^{l-1} \textbf{a}+{}^{l-1}\textbf{b}\big{)},\] \[{}^{L}\textbf{a}={}^{L-1}\textbf{w}\cdot{}^{L-1}\textbf{a},\]

Where \({}^{0}\textbf{a}\) and \({}^{l}\textbf{a}\) denote the input layer (the 0\({}^{\text{th}}\) layer) and the output layer (the \(L^{\text{th}}\) layer), respectively. \(l=1\), \(2\), \(\ldots\), \(L{-}1\) denotes the number of hidden layers. **w** and **b** are the weights and biases in the FNN. \(\varSigma\) denotes the activation function that adds nonlinearity to the FNN [40]. More details regarding the activation function are given in the following context. On the right-hand side of Figure 2, the outputs of the FNN are used to formulate the physics-informed loss function. It is worth noting that, when solving (partial differential equations) PDEs via FNNs, the partial differential terms in the target PDEs can be analytically obtained through automatic differentiation [41]. Through training algorithms [42], the parameters of the FNN are iteratively modified to decrease the value of the physics-informed loss function. Finally, the training process converges when the given criteria are satisfied.

Figure 2: An example of a physics-informed neural network (PINN).
To calculate the nonsmooth dynamic problem, two fundamental solution problems are involved, namely, the solution of the LCP equation and the calculation of the dynamic response. In what follows, different neural network structures for the LCP equation and dynamics response are introduced in detail, respectively.

#### 2.3.1 PINN for solving LCP equations

This part presents the idea of combining the neural network and the theoretical description of LCP. The general LCP equation is in the form of Eq. (16):

\[\begin{cases}\mathbf{y}=\mathbf{Ax}+\mathbf{b}\\ \mathbf{y}^{\top}\mathbf{x}=\mathbf{0}\\ \mathbf{x}\geq\mathbf{0,y}\geq\mathbf{0}\end{cases} \tag{16}\]

in which \(\mathbf{x}\) and \(\mathbf{y}\) are vectors.

Instead of the conventional methods (pivoting or the quadratic programming method (QP method)) to get the solutions, we define two functions:

\[\mathbf{f}=\mathbf{y}-\mathbf{Ax}+\mathbf{b} \tag{17}\]

\[\mathbf{r}=\mathbf{y}^{\top}\mathbf{x} \tag{18}\]

The loss function in the following form is defined:

\[\mathcal{L}=\frac{\sum\limits_{i=1}^{N}f_{i}^{2}+\sum\limits_{i=1}^{N}r_{i}^{2} }{N} \tag{19}\]

in which \(f_{i}\) and \(r_{i}\) are respectively the element in vector \(\mathbf{f}\) and \(\mathbf{r}\), and \(N\) is the element number of the vector. Physics-informed loss function \(\mathcal{L}\) serves as the optimisation object in the training process of the neural network.

The architecture of the PINN for solving the LCP equation (Eq. (17)) is shown in Figure 3. By starting with \(a_{0}\), \(a_{1}\),..., \(a_{n}\) as the input, the neural network is used for predicting the output data \(\mathbf{x}\) and \(\mathbf{y}\) data. In this procedure, several variables need to be determined. One critical choice is the activation function between the layers. Several nonlinear activation functions can be used, such as the tanh, mish and ReLU, shown in Figure 4.
As **x** and **y** in LCP need to be either positive or zero values, the ReLU function given in Figure 3 can provide such an output satisfying the LCP condition. However, the output of the ReLU function has a flat area when inputs are located at the left side of the origin point. In this case, the training of neural networks may trap by the flat area, resulting in more computational time to convergence. To alleviate this problem, a new activation function is designed based on the ReLU function, which is in the following form:

\[g=\max\left(0,a+c_{1}\right)+c_{2} \tag{20}\]

where \(c_{1}\) and \(c_{2}\) are the correction factors. The modified ReLU function could adjust the data near the origin point, shown in Figure 4 (d). By using the modified ReLU function, the output of initialised neural networks is more likely fall out of the flat area, which largely saves the training time.

After obtaining the outputs \(\mathbf{x}\) and \(\mathbf{y}\), the optimisation of the physics-informed loss function is accomplished by L-BFGS algorithm. As shown in the overall procedure of the LCP PINN (Figure 3), the training ends until the error is smaller than the tolerance.

Figure 3: The process of the PINN framework for LCP.
Several examples are calculated by the LCP PINN framework as well as two conventional methods (Pivoting and QP method) for comparison, with \(A=\begin{bmatrix}1&-1\\ -1&0\end{bmatrix}\) and \(b=\begin{bmatrix}-0.009&0.02\end{bmatrix}^{T}\). For the sake of conciseness, the solutions are listed in Table 1. This example uses a 3-layer neural network with 5 neurons at each layer. It is proved that the proposed LCP PINN framework can give accurate solutions.

\begin{table}
\begin{tabular}{c c c c c c} \hline  & Accurate result & PINN framework & Pivoting method & QP method \\ \hline
**x** & [0.009 0] T & [0.009 0] T & [0.009 0] T & [0.0090004 4.8204583e-09] T \\ \hline
**y** & [0 0.011] T & [0 0.011] T & [0 0.011] T & [0 0.011] T & [4.0007060e-07 0.0109996] T \\ \hline \end{tabular}
\end{table}
Table 1: The solutions by different methods

Figure 4: The nonlinear activation functions: (a) Tanh function; (b) Mish function; (c) ReLu function; (d) the modified ReLu function.
#### 2.3.2 PINN for vibration of multibody systems with unilateral contact

For the initial condition problem, the second-order differential equation is general in the form \(\ddot{q}=f\left(t,q,\dot{q}\right)\), and the initial conditions are \(q\left(t_{0}\right)=q_{0}\) and \(\dot{q}\left(t_{0}\right)=\dot{q}_{0}\). Implicit Runge-Kutta method is an high-order integral method to solve the initial condition problem. The \(R^{\text{th}}\)-order implicit Runge-Kutta method can be expressed as a transformation form, given in Eqs. (21) and (22):

\[\dot{q}_{k}^{\left(i\right)}=\dot{q}_{k}^{\left(i+\kappa_{k}\right)}-\Delta t \sum_{r=1}^{R}a_{\lambda_{r},f}\left(t^{\left(i\right)}+c_{r}\Delta t,\dot{q}_ {r}^{\left(i+\kappa_{r}\right)}\right)\text{, }k=1,2,\ldots R \tag{21}\]

\[\dot{q}^{\left(i\right)}=\dot{q}^{\left(i+1\right)}-\Delta t\sum_{k=1}^{R}b_{ \lambda_{r}}f\left(t^{\left(i\right)}+c_{r}\Delta t,\dot{q}_{k}^{\left(i+ \kappa_{k}\right)}\right) \tag{22}\]

in which \(\Delta t\) is one time increment over discrete time, \(a_{\lambda_{r}}\), \(b_{\lambda_{k}}\) and \(c_{r}\) are the coefficients. In Eqs. (21) and (22), the known term at time step \(I\) are expressed by the unknown terms at time step \(i\)+1. When a neural network is introduced, using \(q\) as the input data, \([\begin{array}{cccc}\dot{q}_{k}^{\left(i+\kappa_{0}\right)},\dot{q}_{k}^{ \left(i+\kappa_{0}\right)},\ldots,&\dot{q}_{k}^{\left(i+1\right)}\end{array}]\) can be predicted which is designed as the output data. Subsequently, the neural network is trained until \(\dot{q}_{r}^{\left(i\right)}\) and \(\dot{q}^{\left(i\right)}\) obtained by Eqs. (21) and (22) converge to the corresponding values that has been known.

As to multi-degree-of-freedom systems, the general form of the equation of motion is expressed as:

\[\mathbf{M}\ddot{\mathbf{q}}+\mathbf{C}\dot{\mathbf{q}}+\mathbf{K}\mathbf{q}= \mathbf{f}_{\text{e}} \tag{23}\]

in which \(\mathbf{M}\), \(\mathbf{K}\), \(\mathbf{C}\) are the mass, stiffness and damping matrix, \(\mathbf{q}\) is the displacement vector and \(\mathbf{f}_{\text{e}}\) is the external force. \(\mathbf{v}\) refers to as \(\dot{\mathbf{q}}\). The conventional way to solve this second-order differential equation of dynamical systems is using numerical integral methods, such as the explicit/implicit Runge-Kutta (RK) methods. Here, Figure 5 shows the framework of a PINN for dynamic simulation of multi-degree-of-freedom systems. Although the time integral idea of the current PINN is based on the implicit R-K method but as introduced above, it is actually different from conventional implicit R-K method.

In the framework, the input data are the elements of \(\mathbf{q}\)= [ \(q_{1}\),... \(q_{l}\),..., \(q_{N}\)]T (\(n\)=1,2,..., \(N\)) and the network is designed to output \(N\times R\) number of data. \(N\) is the number of the system degrees of freedom, and \(r\) is the number of orders used for the time integral. Then the output data is reshaped to \(N\) vectors that each has \(R\) elements. Each vector is assigned to the velocity \(\mathbf{v}_{r}^{\left(i+1\right)}=[\begin{array}{cccc}\dot{q}_{r}^{\left(i+ \kappa_{0}\right)},\dot{q}_{r}^{\left(i+\kappa_{0}\right)},\ldots,&\dot{q}_{r }^{\left(i+1\right)}\end{array}]T\) corresponding to \(r\) number of quantities of each degree-of-freedom. After obtaining the physics
quantities \(\mathbf{v}_{{}_{N}}^{{}^{(i+1)}}\), the displacement and the acceleration at time step \(i\)+1 can be obtained through the dynamic equations (Eqs. (i)-(iii)) shown in Figure 5. **R** is the coefficients matrix about \(a_{{}_{x},{}_{x}}\) and \(b_{{}_{r}}\). Thus, according to the idea of Eqs. (21) and (22), which is written as a matrix form Eq. (iv) in Figure 5, the velocity vectors \(\mathbf{v}_{{}_{n}}^{{}^{(i)}}=[\mathbf{v}_{{}_{1,n}},...,\mathbf{v}_{{}_{r,n} },...,\mathbf{v}_{{}_{R+1,N}}]\) can be obtained. The displacement and velocity can be learned by optimizing the loss function in the following form:

\[\mathcal{L}_{\text{D}}=\sum_{n=1}^{N}\sum_{r=1}^{R+1}\left(\mathbf{v}_{{}_{r,n }}-\mathbf{v}^{{}^{(i)}}\right) \tag{24}\]

The coefficient tables from the 2\({}^{\text{th}}\)-order to 200\({}^{\text{th}}\)-order implicit Runge-Kutta method can be obtained in [https://github.com/maziarraissi/PINNs](https://github.com/maziarraissi/PINNs)[32].

#### 2.3.3 Single PINN for the transient dynamic analysis of nonsmooth vibration

When unilateral contact is considered, nonsmooth characteristics are involved in the dynamic system. Figure 6 shows the framework of the specific PINN, which is referred to as the single PINN framework in the following. Compared with the PINN for dynamic simulation given in Figure 5, the equations of motion Eqs. (b) and (c) used in the neural network are:

\[\mathbf{h}_{{}_{n}}^{{}^{(i+1)}}=\mathbf{K}_{{}_{n}}\mathbf{q}_{{}_{n}}^{{}^{( i+1)}}+\mathbf{C}_{{}_{N}}\mathbf{v}_{{}_{n}}^{{}^{(i+1)}}+\mathbf{F}_{{}_{e}} \tag{25}\]

\[\mathbf{a}_{{}_{n}}^{{}^{(i+1)}}=\mathbf{M}^{-1}\mathbf{h}_{{}_{n}}^{{}^{(i+1) }}+\mathbf{W}_{{}_{N}}\mathbf{\lambda}_{{}_{N}}+\mathbf{W}_{{}_{T}}\mathbf{ \lambda}_{{}_{T}} \tag{26}\]

Figure 5: The framework of PINN for dynamic simulation.
in which \(\mathbf{K}\), and \(\mathbf{C}_{\mathrm{s}}\) are the stiffness and damping matrix when the unilateral constraints are removed, and the contact force \(\lambda_{\mathrm{s}}\) and friction force \(\lambda_{\mathrm{\tau}}\) are time-varying variables, which can be determined by the LCP equations introduced in Section 2.2.1 or 2.2.2.

The calculation process of the single PINN framework is: (1) use the system information, such as the system parameters and the initial conditions to calculate \(\lambda_{\mathrm{s}}\) and \(\lambda_{\mathrm{\tau}}\), based on the time-stepping method according to the contact type; (2) transfer \(\lambda_{\mathrm{s}}\) and \(\lambda_{\mathrm{\tau}}\) to the dynamic equation part in the PINN for dynamic calculation. The operation mechanism of the PINN for dynamic calculation has been explained in previous section. In terms of the neural network design, the activation function for the dynamic calculation is tanh. The dynamic response can be learned by optimising the loss function \(\mathcal{L}_{\mathrm{D}}\) given in Eq. (24).

#### 2.3.4 Dual PINN for the transient dynamic analysis of nonsmooth vibration

One the basis of the LCP PINN for linear complementary problem proposed in Section 2.3.1 and the PINN for multibody systems with unilateral contact, a new physics-informed neural network is proposed, which is shown in Figure 7. In this method, two physics-informed networks are involved. Thus, this method is referred to as the dual PINN framework.

As the normal contact and tangential contact can be enforced through LCP, the first neural network is basically the LCP PINN. Compared with the PINN for LCP in Section 2.3.1, the outputs of the

Figure 6. The framework of the single PINN framework.
FNN for LCP are assigned to physical variables. It should be noticed that for different contact types, the physical variables are different: (1) for direct rigid contact (Section 2.2.1), outputs are assigned to \(\Lambda_{\mbox{\tiny N}}\), \(\Lambda_{\mbox{\tiny RE}}\), \(\Lambda_{\mbox{\tiny LE}}\), \(\Lambda_{\mbox{\tiny LE}}\), \(\mathbf{g}_{\varepsilon}\), \(\mathbf{r}_{\mbox{\tiny LE}}\) and \(\mathbf{r}_{\mbox{\tiny RE}}\); (2) for rigid contact with springs (Section 2.2.2), outputs are assigned to \(\Lambda_{\mbox{\tiny N}},\Lambda_{\mbox{\tiny RE}}\), \(\Lambda_{\mbox{\tiny LE}}\), \(\mathbf{\Omega}_{\mbox{\tiny E}}\), \(\mathbf{r}_{\mbox{\tiny LE}}\) and \(\mathbf{r}_{\mbox{\tiny RE}}\). Then \(\mathbf{x}\!=\!\left[\Lambda_{\mbox{\tiny N}}\)\(\Lambda_{\mbox{\tiny LE}}\)\(\mathbf{r}_{\mbox{\tiny RE}}\right]^{\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!
When contact force \(\lambda_{{}_{\rm N}}\) and friction force \(\lambda_{{}_{\rm T}}\) are outputted, they are transferred to the physics part of the second PINN for the dynamic analysis of the multibody system with unilateral contact. The basic idea of this PINN has been introduced in Section 2.3.4. Eventually, the dynamic response during an arbitrary time duration can be learned through training.

## 3 One-dimensional stick-slip vibration with direct contact

This section presents the application of PINN frameworks to solve the stick-slip vibration of Figure 8. It is assumed that the contact in the normal direction is maintained during vibration. The vertical position of the slider mass is a constant value of \(0.\lambda_{{}_{\rm N}}\)is known in advance, which is the compression force. The equations of motion for systems with friction contact can be reduced to:

\[\mathbf{M}\ddot{\mathbf{x}}+\mathbf{h}+\mathbf{W}_{\rm T}\dot{\lambda}_{{}_{ \rm T}}=\mathbf{0} \tag{27}\]

The LCP equation given in Eq.(12) is degenerated to:

\[\underbrace{\left\{\begin{matrix}\tau_{{}_{\rm L}}^{\rm(E)}\\ \mathbf{\Lambda}_{\rm R}^{\rm(E)}\Delta t\end{matrix}\right\}}_{\rm y}= \underbrace{\left[\begin{matrix}\mathbf{G}_{{}_{\rm TT}}&\mathbf{I}\\ -\mathbf{I}&\mathbf{0}\end{matrix}\right]}_{\rm A}\underbrace{\left[ \begin{matrix}\mathbf{\Lambda}_{\rm L}^{\rm(E)}\\ \mathbf{\gamma}_{\rm R}^{\rm(E)}\end{matrix}\right]}_{\rm A}+\underbrace{ \left[\begin{matrix}-\left(\mathbf{G}_{{}_{\rm TN}}+\mathbf{G}_{{}_{\rm TT}} \mathbf{\mu}\right)\mathbf{\Lambda}_{\rm N}^{\rm(E)}-\gamma_{{}_{\rm T}}^{\rm(A )}-\left(\mathbf{G}_{{}_{\rm T}}\mathbf{h}^{\rm(A)}+\frac{\partial\mathbf{W}_{ \rm T}^{\rm T}}{\partial\mathbf{q}}\mathbf{u}^{\rm(A)}+\dot{\mathbf{\phi}}_{{}_{ \rm T}}\right)\Delta t\\ \mathbf{2}\mathbf{\mu}\mathbf{\Lambda}_{\rm N}^{\rm(E)}\end{matrix}\right]}_{ \rm\dot{\mathbf{h}}} \tag{28}\]

with \(0\leq\gamma_{{}_{\rm L}}^{\rm(E)}\perp\Lambda_{\rm L}^{\rm(E)}\geq 0\) and \(0\leq\gamma_{{}_{\rm R}}^{\rm(E)}\perp\Lambda_{\rm R}^{\rm(E)}\geq 0\).

### Mechanical model

The novel PINN approach proposed here is demonstrated by simulating a common 1-dimensional stick-slip friction and vibration problem representing a slider-belt system (Figure 8). The system is modelled as a point mass \(m\) held by a linear spring and viscous damper, supported by a rigid belt moving at a constant velocity \(v_{0}\) whereas \(k_{0}\) is the spring stiffness, \(c_{0}\) is the damping, \(F_{\rm n}\)is the normal compression force, and \(\mu\) is the coefficient of friction, assumed constant. The coordinate \(x\) describes the horizontal displacement. The vertical position of this model is a constant value assumed to be 0. \(\lambda_{{}_{\rm N}}=F_{{}_{\rm n}}\) is known in advance. The constraint function in the tangential direction are expressed as \(g_{{}_{\rm T}}=\dot{x}-v_{o}t\). According to the definitions, one gets: \(W_{{}_{\rm T}}=1\), \(\widetilde{w}_{{}_{\rm N}}=\widetilde{w}_{{}_{\rm N}}=\widetilde{w}_{{}_{\rm T }}=0\), \(\widetilde{w}_{{}_{\rm T}}=-v_{o}\).
The two kinds of Coulomb-Stribeck laws of friction are (Figure 9):

\[\mu=\frac{\mu_{\ast}}{1+\delta\big{|}v_{\ast}\big{|}}\,,\,\,\,\,\mu=\mu_{\ast}+ \big{(}\mu_{\ast}-\mu_{\ast}\big{)}e^{-\alpha\mathbf{k}\cdot\mathbf{l}} \tag{29}\]

### PINN scheme

When the system is nonlinear/nonsmooth, the requirement for accuracy is stricter as the results of the next step are sensitive to the results of the previous step. To explore the applicability and accuracy of the new PINN strategy, two specific PINN frameworks for the single/dual PINN strategy are used for stick-slip simulation. Only the friction force is calculated from the conventional LCP algorithm involved in the single PINN framework or the LCP PINN part in the dual PINN framework. The LCP formulations are now Eq. (28).

Figure 8: Model I: the single mass-belt system with friction.
### Numerical Simulation

For the single-point frictional oscillator with stick-slip motion, the numerical methods, for example, the smoothing method or the bisection method, can provide very accurate results. Here, the switching model method gives ground truth results. In this section, the mix-level time-stepping method presented in Section 2 is used to simulate stick-slip, which is called the conventional method. In the numerical simulations, system parameters are \(m\)=1, \(k_{0}\)=1, c\({}_{0}\)=0.

#### 3.3.1 Friction law with a slowly decreasing rate

Firstly, the PINN frameworks are examined in a situation when the friction force gently decays with the relative velocity. The nonlinear friction law that is employed has been given in Eq. (29) and depicted in Figure 9. Figure 10 shows the results of the phase plot of the vibration calculated by two PINN strategies. The black line with the triangle markers is about the ground truth results calculated by the switching method. The green line with the square markers is about the result of the single PINN framework, and the blue dashed line with a circle mark is about the result of the dual PINN framework. The straight line in the phase plot indicates stick motion, and the curve represents the slip motion. The vibration of the system alternates between the stick phase and the slip phase.

The comparison in Figure 10 shows that the results from both the single PINN framework and dual PINN framework are identical to the ground truth results. This indicates that the idea of employing nonsmooth dynamics theory for the multi-point contact problem in a neural network works very well

Figure 10: The phase plot of stick-slip vibration.
in dealing with nonlinear stick-slip vibration. The 4th-order PINN and 10th-order PINN give very close results. Thus, only one group of results are provided.

The dynamic behaviour of nonlinear systems with stick-slip or loss of contact (that would appear in 2-DoF models) depends on the initial condition. Hence, capturing the transition point, which provides the initial conditions for the next phase, is important. Time-stepping method does not require the transition check or event check, while the price is a sufficiently small time step required for high accuracy. Figures 11-12 illustrate the comparisons between conventional time-stepping methods and two PINN frameworks when the time step \(\Delta t\)=0.01. The calculation of the conventional LCP method and the conventional RK4 LCP method is based on a numerical integration method. The latter is a high-order numerical method. Figure 11 shows that when \(\Delta t\)=0.01, the conventional LCP method could still predict the transitions between stick and slip motion. However, higher-order RK4 LCP method, in this case, is not stable enough to predict the correct stick-slip. Notably, both the single and Dual PINN frameworks are high-order methods and could produce the correct stick-slip limit cycle. Furthermore, from Figure 12, it is known that the accuracy of PINN frameworks in both the displacement and the velocity level is quite good. The conventional RK4 LCP method could initially predict the correct results until the failure of detecting the second transition from slip to stick around 18s (illustrated in Figure 12 (b)). As high-order algorithms, the PINN frameworks are more reliable.

Figure 11: The phase plot of stick-slip vibration when \(\Delta t\)=0.01.
The frequency spectrum shown in Figure 13 indicates that the system with stick-slip motion has one primary vibration frequency, one doubling frequency and one tripling frequency. From the frequency point of view, when a method (Conventional RK4 LCP method) is incapable of simulating correct stick-slip transitions, it not only predicts distorted time history responses but also provides incorrect vibration frequencies. Stick-slip vibration has a higher requirement for numerical methods.

When the time step increases to \(\Delta t\)=0.08, Figure 14 shows that the two PINN frameworks could maintain the correct prediction of stick-slip vibration. However, the two conventional methods could

Figure 12. Time history of the displacement and the velocity \(\Delta t\)=0.01.

Figure 13. The frequency of stick-slip vibration \(\Delta t\)=0.01.
not give accurate results. Figure 14 (b) is the time history of the velocity, in which the constant velocity means that the system is in the stick phase. Otherwise, the slider is slipping. Figure 14 (b) explains the reason for the failure of the conventional methods. With the large time step, the conventional method cannot identify the transition instants between stick and slip motion.

#### 3.3.2 Friction law with a sharply decreasing rate

In this section, the simulations using the friction law shown in Figure 9 (b) are conducted. The sharp decrease of the friction force with a velocity near zero velocity reflects the friction behaviour of many real frictional pairs. This kind of friction law causes more difficulties in numerical calculations. On the one hand, using a smaller time step could ensure the detection of stick-slip transitions which has been known as a key factor for the accuracy of nonsmooth vibration. On the other hand, the time step could not be too small, considering the expenses of numerical calculation. In the following examples, the neural network that is used in the PINN for dynamic simulation is a 7-layers feedforward neural network in which each layer contains 20 neurons. The network used in the PINN for LCP is a 3-layers feedforward neural network in which each layer contains 5 neurons.

Figure 15 shows the results of different methods at the time step that could predict correct stick-slip vibration. As illustrated in Figure 15 (b), the conventional high-order method needs a very small time step compared with the PINN frameworks. The dual PINN framework needs a time step that is 50 times larger than the Conventional RK4 LCP method. Although the advantage of the Single PINN

Figure 14: The time response of stick-slip vibration when \(\Delta t\)=0.08.
framework in terms of the time step is not as outstanding as the dual PINN framework, it could predict stick-slip vibration accurately when the time step is not too small.

Figure 16 compares the time history results of the different methods based on LCP strategy for multi-contact systems when \(\Delta t\)=0.005. Figure 16 (b) shows that the Single PINN framework and the two conventional methods at the same time step fail to switch to stick motion after the slider starts to slide on the belt. Consequently, the wrong displacement track begins, as shown in Figure 16 (a). The above results indicate the importance of the numerical method in correct state detection and the excellent performance of Dual PINN.

Figure 15: The correct phase plot of stick-slip by different methods: (a) conventional LCP \(\Delta t\)=0.0001; (b) Conventional RK4 LCP when \(\Delta t\)=10\({}^{\text{-6}}\); (c) Single PINN when \(\Delta t\)=0.0001; (d) Dual PINN when \(\Delta t\)=0.005.
## 4 Two-dimensional nonsmoothness FIV with spring contact

### Mechanical model

Figure 17 is one classic and minimal model showing stability bifurcation, known as mode-coupling instability. It is a 2-DoF model consisting of a mass-spring-damper sliding on a rigid moving belt. In the horizontal direction, the mass \(m\) is held by a spring \(k_{1}\) and a damper \(c_{1}\). In the vertical direction, \(m\) is linked with a ground damper \(c_{2}\). Additionally, \(m\) is constrained by an inclined spring \(k_{3}\), which would cause an asymmetry stiffness that accounts for mode-coupling instability. Between the mass and the belt, the contact stiffness is \(k_{2}\). An assumed massless slider is in contact with the moving belt. A preload \(F_{\text{p}}\) is applied to bring the slider into contact with the belt before the belt starts to move. The friction at the interface is assumed to follow the Coulomb-Stribeck friction law. The sliding friction force \(F_{\text{T}}\) is proportional to the normal contact force.

The coefficient of friction \(\mu\) affects the natural frequencies of the 2-DoF system [35, 36]. As \(\mu\) increases, the two natural frequencies of the system move closer and eventually merge causing instability of the system. As horizontal vibration grows, the interface between the slider and the belt opens, due to vertical motion, before contant is restored at a later time. The transition of the states in the normal direction is called the separation-reattachment phenomenon. Additionally, when the Coulomb-Stribeck friction law is considered, the slider would undergo horizontal stick-slip vibration

Figure 16: The time response of stick-slip vibration with friction law when \(\Delta t\)=0.005: (a) the time history of the displacement; (b) the time history of the velocity.
at the interface, potentially promoting further instability of the system. Both types of nonsmooth phenomena can be present in this model, which is unlike the one-degree-of-freedom model in Section 3.

Figure 17. Model II: 2-DoF slider-belt model with dry friction.
dual PINN framework, an interpolation technique is used. According to the order of the numerical integration method, the force vector \(\lambda_{\mathrm{r}}^{\left(i\right)}\) and \(\lambda_{\mathrm{n}}^{\left(i\right)}\) in conventional methods are now transformed to \(\lambda_{\mathrm{r}\cdot\mathbf{x}\mathrm{x}}^{\left(i\right)}\) and \(\lambda_{\mathrm{n}\cdot\mathbf{x}\mathrm{x}}^{\left(i\right)}\) matrices based on the following equations:

\[\lambda_{\mathrm{r}\cdot\mathbf{x}\mathrm{x}}^{\left(i\right)}=\lambda_{ \mathrm{r}}^{\left(i-1\right)}+\mathbf{s}\left(\lambda_{\mathrm{r}}^{\left(i \right)}-\lambda_{\mathrm{r}}^{\left(i-1\right)}\right) \tag{30}\]

\[\lambda_{\mathrm{n}\cdot\mathbf{x}\mathrm{x}\mathrm{x}}^{\left(i\right)}= \lambda_{\mathrm{n}}^{\left(i\right)}+\mathbf{s}\left(\lambda_{\mathrm{n}}^{ \left(i\right)}-\lambda_{\mathrm{n}}^{\left(i-1\right)}\right) \tag{31}\]

in which \(\mathbf{s}\) is the vector formed by the time segment calculated based on the coefficient c\({}_{k}\) (\(k\)=1,2,...,\(q\)) of the \(q^{\mathrm{th}}\)-order RK method. Then \(\lambda_{\mathrm{r}\cdot\mathbf{x}}^{\left(i\right)}\) and \(\lambda_{\mathrm{n}\cdot\mathbf{x}}^{\left(i\right)}\) are passed into the physics part of the PINN for dynamic response \(\lambda_{\mathrm{r}}^{\left(i\right)}\) and \(\lambda_{\mathrm{n}}^{\left(i\right)}\) in Eq. (c) in Figures 5 and 6. The new methods are referred to as the advanced Single PINN and the advanced Dual PINN in the following. The solution procedure of the two advanced PINN frameworks is shown in Figure 18.

### Numerical simulation

Although stick-slip and separation-reattachment are allowed in this problem, the occurrence of these phenomena depends on the system parameters. As mode-coupling instability happen, the critical coefficient of friction \(\mu_{\mathrm{c}}\) at bifurcation can be obtained via eigenvalue analysis. Figure 19 shows the change of the real and imaginary parts of the eigenvalues during sliding contact with the coefficient

Figure 18: Framework of the advanced Single/Dual PINN for the 2D nonsmooth problem.
of friction for Example 1 (see section 4.3.1). Instability occurs at \(\mu_{\epsilon}=\)0.83 (Figure 19). Two examples with parameters \(\mu=\)0.4 and \(\mu=\)1.2, are investigated.

#### 4.3.1 Example 1

By using the single PINN and dual PINN schemes, numerical simulations are completed by setting the initial condition to \(x(t=0)=y(t=0)=\)-10, and \(v_{x}(t=0)=v_{y}(t=0)=v_{0}\).The size of the neural network is identical to the problem in Section 3 (Figure 20). The results of the root-shooting method are used as the ground truth results. In Figures 20 and 21, the solid and hollow squares are for the advanced single PINN method, when a 4\({}^{\text{th}}\)-order and a 10\({}^{\text{th}}\)-order integration formulations are used in the PINN for dynamic simulation. The solid triangles and hollow circles are for the advanced dual PINN method when a 4\({}^{\text{th}}\)-order and a 10\({}^{\text{th}}\)-order integration formulations are involved.

Figure 20 shows the time history of the contact force. In this example, the contact force alternates between positive values, which represent contact between the slider and the belt, and zero value, which represents the loss of contact. Thus, separation and reattachment happen during the vibration, which brings vertical nonsmoothness. From the results, one can see that both the advanced single PINN and the advanced dual PINN framework can produce very good results for a problem when repeated separation-reattachment events happen. As illustrated in the enlarged plot of the contact force (Figure 20 (b)), when the 4\({}^{\text{th}}\)-order PINN is used, the result calculated by the dual PINN framework is not as good as that of the single PINN. When the 10\({}^{\text{th}}\)-order integration method is used, the accuracy of the dual PINN framework (red circles) can be improved from 4\({}^{\text{th}}\)-order dual PINN

Figure 19: The evolution of the real part (left) and imaginary part (right) of the eigenvalue with the friction coefficient.
(blue triangle). Similarly, the accuracy of the single PINN framework can be also improved by using a high-order integration method.

Figure 21 illustrates that the displacement and velocity responses calculated by the two advanced PINN strategies are very accurate compared with the root-shooting results. As shown by the time history of the velocity given in Figure 21 (c)-(d), stick-slip vibration does not appear in this case. These results enable comparisons between the PINN strategies that employ the dynamic equations with unilateral contact and the ground truth results, which are calculated by the rooting-searching method for the single-point contact problem. In the following, the comparisons between the conventional numerical methods for the multi-point contact problem with the PINN strategies are made.

Figure 20: The time history of the contact force. (a) \(t\)=[0\(\sim\)10]s; (b) \(t\)=[9\(\sim\)10]s.
Table 2 shows the root mean squares (RMS) of the contact force, the displacement and the velocity calculated by different methods. The difference between the conventional LCP method and the RK4 LCP method is that the former uses Morue's algorithm, and the latter adopts the idea of 4th-order RK. Single PINN and Dual PINN are the methods without the interpolation that are introduced in Sections 2.2.3 and 2.2.4. To give a clearer illustration of the data in Table 2, Figure 22 shows the RMS error of the horizontal and vertical displacement. By using the results of the root-shooting method as the reference value, the absolute errors of the corresponding RMS results of the six methods for nonsmooth simulation of multi-point contact problems are calculated. One can find that: (1) with the interpolation technology, the advanced PINN frameworks improve the accuracy of the PINN framework. (2) compared with the conventional methods, the advanced PINN frameworks have the advantage in accuracy when the same time step is used ( \(\Delta\)\(t\)=0.001). The accuracy of the conventional methods could be improved when the time step is smaller ( \(\Delta\)\(t\)=0.0005), which means lower efficiency. As shown by the velocity data listed in Table 2, the proposed advanced single and dual PINN methods can produce accurate results.

\begin{table}
\begin{tabular}{c c c c c c} \hline \hline MethodRMS & Contact force & Horizontal displacement & Vertical displacement & Horizontal velocity & Vertical velocity \\ \hline \hline \end{tabular}
\end{table}
Table 2: RMS results of different numerical methods.

Figure 21: The time history of displacement and velocity: (a) horizontal displacement vs time; (b) vertical displacement vs time; (c) horizontal velocity vs time; (d) vertical velocity vs time.
Figure 22: The RMS error of the displacement response:(a) horizontal displacement; (b) vertical displacement.
#### 4.3.2 Example 2

Figures 23 and 24 show the response of the contact force and the horizontal velocity when \(\mu\)=1.2 and the initial conditions are \(x(t=0)=y(t=0)=-1\), and \(v_{x}(t=0)=v_{y}(t=0)=v_{0}\). The first observation is that separation and reattachment happen in this example, as the contact force repeatedly alternants between positive values and zero representing contact and loss of contact, respectively. Then from Figure 24, it can be seen that stick-slip vibration in the horizontal direction occurs as well. More clearly, Figure 24 (b) shows that when both normal and tangential nonsmoothness are involved in a problem, the state change is unpredictable. Figure 23 shows whether the advanced single PINN or advanced dual PINN framework can accurately detect the separation-reattachment events. The zoomed-in view is presented in Figure 23 (b). Compared with the results of the 10th-order advanced single and dual PINN cases, adding one PINN could lower the accuracy, which also happens in Example 1. However, when a higher-order numerical integration scheme is used, the accuracy of both single and dual PINN frameworks is improved, and it is not a big concern. As shown in Figure 24, stick-slip vibration occurs from around 4.6s. Consistently with the contact force results, when the order of the numerical method used is ten, accurate stick-slip vibration can be produced by both the advanced Single and dual PINN frameworks. Even though the dynamic response is highly nonlinear, PINN frameworks work quite well.
In the following, the results of the new PINN frameworks and the conventional methods that all nonsmooth calculations of multi-point contact systems are compared. In Table 3, numerical calculations of the two conventional methods are carried out with different time steps (\(\Delta\,t=\)0.001, \(\Delta\,t=\)0.0005, and \(\Delta\,t=\)0.0001). \(\times\) denotes the failure of detection of stick-slip transitions, which in the view of nonsmooth simulations is unacceptable. When \(\Delta\,t=\)0.0001, conventional methods could give satisfied results. With respect to the PINN frameworks, when the 4th-order methods are used, a smaller time step is needed as when \(\Delta\,t=\)0.001 would not correctly predict stick-slip transitions at 8.88 second (see Figure 24(b)). Furthermore, the 10th-order single and dual PINN frameworks can make an accurate transition detection at \(\Delta\,t=\)0.001, which is 10 times larger than the time step of the conventional methods, and the accuracy of the overall time- domain response is high.

\begin{table}
\begin{tabular}{l c c c c} \hline \hline Method & RMS & Contact force & Horizontal displacement & Vertical displacement \\ \hline Root-Shooting Method & & 224.449 & 0.3997 & 0.7629 \\ \hline Conventional LCP method & \(\Delta\,t=\)5\(\times\)10\({}^{-4}\) & \(\times\) & \(\times\) & \(\times\) \\  & \(\Delta\,t=\)10\({}^{-4}\) & 222.331 (0.94\%) & 0.3980 (0.43\%) & 0.7530 (1.30\%) \\ \hline \hline \end{tabular}
\end{table}
Table 3: RMS results of different numerical methods. Cross (\(\times\)) denotes results that fail to present correct stick-slip transitions

Figure 24: The time history of the velocity. (a) \(t\)=[0\(\sim\)10] s; (b) \(t\)=[8.7\(\sim\)9.7]s, in which the dotted zone represents the sticking state, the grey zone represents the slipping state, the white zone represents the separation state.
Through the above numerical analysis of examples 1 and 2, one can conclude that when determining the vibration of a system considering separation-reattachment and stick-slip, the numerical methods and the time step must be chosen very carefully. The error in predicting nonsmooth vibration comes not only from the discernible error accumulation due to the numerical integration but also from the unpredictable failure of detecting nonsmooth transitions, which is a fatal issue. The new PINN frameworks proposed in this work are shown to be well capable of dealing with nonsmoothness in two dimensions (the normal and tangential directions). Moreover, compared with conventional methods, the advanced single and deal PINN frameworks are less dependent on the time step length to give good predictions of the stick-slip and separation-reattachment transitions. Because PINN frameworks rely on the functions from conventional physical equations, in which the numerical error from the traditional numerical methods can be offset through the training of the neural network. The accuracy of PINN frameworks can be improved in two ways: decreasing the time step and increasing the order of the numerical integration formulations in PINN frameworks.

## 5 Conclusions

This work proposes novel physics-informed neural network (PINN) frameworks for simulating the vibration of frictional systems with multiple contact points. Firstly, a tailored PINN based on the mathematical expressions of the linear complementary problem (LCP) is proposed. Then, the PINN for the dynamic simulation of complex systems based on numerical integration formulations of the dynamic equations is established. Furthermore, based on the two PINN frameworks, four PINN-based
methods to replace conventional numerical methods for solving the vibration of frictional systems with multiple contact points are proposed. They are the single PINN framework, the dual PINN framework, and the advanced single/dual PINN framework, respectively.

The applications of these frameworks in solving the direct contact problem with only the stick-slip vibration (1-D problem) and spring contact problem with both separation-reattachment and stick-slip (2-D problem) lead to the following findings: (1) all four PINN frameworks can accurately detect stick-slip and separation-reattachment transition events in both 1-D and 2-D nonsmooth problems; (2) compared with the conventional numerical methods, the PINN frameworks do not necessitate a tiny time step. High-order PINN frameworks allow a larger time step with better accuracy and greater stability than the conventional methods; (3) the advanced single/dual PINN frameworks outperform the single/dual PINN frameworks in identifying state transitions, which is important for achieving high accuracy; (4) the accuracy of PINN strategies can be improved by simply increasing the numerical integration order. Overall, the newly developed methods have illustrated the potential to accurately predict non-smooth dynamic behaviour, indicating encouraging prospects for future research.

## Acknowledgements

The authors are grateful for the financial support from Postdoctoral Research Foundation of China (No. 2019M652564). Support from the Australian Research Council to SM and YG (FT180100338; IC190100020), and the Songshan Laboratory Project (No: 221100211000-01) are also gratefully acknowledged.
