## 4 Sliding Joint  
In this work, the sliding joint is implemented using the Lagrange multiplier method. The constraint of sliding joint is composed of three position constraint equations and one constraint equation, which is expressed according to three Lagrange multipliers and one non-generalized degree of freedom that represents the arc length parameter \(\zeta\). The complete constraint equations can be written as:  
\[\left\{\begin{array}{c}\mathbf{\Phi}=\mathbf{x}_{\text{A}}-\mathbf{x}_{\text {B}}(\zeta)=\mathbf{0}\\ \frac{\partial\mathbf{x}_{\text{B}}}{\partial\zeta}\cdot\boldsymbol{\lambda}= 0\end{array}\right. \tag{50}\]  
where the vectors \(\mathbf{x}_{\text{A}}\) and \(\mathbf{x}_{\text{B}}\) denote the position of joint point on slider and carrier beam, respectively. The constrain \(\mathbf{\Phi}\) in Equation (50) actually leads to a spherical joint, and the Lagrange multiplier \(\boldsymbol{\lambda}=\left[\lambda_{1},\ \lambda_{2},\ \lambda_{3}\right]^{ \text{T}}\) is employed to connect with the constraint forces applied to both bodies. When the vector \(\boldsymbol{\lambda}\) is perpendicular to the sliding direction, the constraint force in this direction is released so that the joint point can move along the midline of carrier beam under the frictionless condition.  
In the simulation, denoting \(\mathbf{q}^{\boldsymbol{\wedge}\text{T}}=\left[\mathbf{q}^{\text{T}},\ \zeta\right]\) as the generalized coordinates, the dynamics equilibrium equations in Equation (48) would be further written as:  
\[\left\{\begin{array}{c}\mathbf{F}^{\text{int}}(\mathbf{q})+\mathbf{F}^{ \text{ine}}(\mathbf{q})-\mathbf{F}^{\text{ext}}(\mathbf{q})+\mathbf{\Phi}^{ \text{T}}_{\boldsymbol{\lambda}}\boldsymbol{\lambda}=\mathbf{0}\\ \frac{\partial\mathbf{x}_{\text{B}}}{\partial\zeta}\cdot\boldsymbol{\lambda}=0 \\ \mathbf{\Phi}(\mathbf{\dot{\alpha}},t)=\mathbf{0}\end{array}\right. \tag{51}\]  
The Jacobian matrix of the equilibrium equation can be evaluated as:  
\[\boldsymbol{\Delta_{\text{Jac}}}=\left[\begin{array}{cc}\mathbf{K}^{\text{ int}}+\mathbf{K}^{\text{ine}}-\mathbf{K}^{\text{ext}}+(\mathbf{\Phi}^{\text{T}}_{ \boldsymbol{\lambda}})&\mathbf{\lambda}&\frac{\partial\mathbf{\Phi}^{\text{T}}_ {\boldsymbol{\lambda}}}{\partial\zeta}\boldsymbol{\lambda}\\ (\frac{\partial\mathbf{x}_{\text{B}}}{\partial\zeta})_{\mathbf{q}}\cdot \boldsymbol{\lambda}&\frac{\partial^{2}\mathbf{x}_{\text{B}}}{\partial\zeta} \cdot\boldsymbol{\lambda}\end{array}\right]&\mathbf{\Phi}^{\text{T}}_{ \boldsymbol{\lambda}}\\ \mathbf{\Phi}_{\boldsymbol{\lambda}}&\mathbf{0}\end{array}\right] \tag{52}\]  
in which the matrices \(\mathbf{K}^{\text{int}}\), \(\mathbf{K}^{\text{ine}}\), and \(\mathbf{K}^{\text{ext}}\) are interpreted as the Jacobian matrix of generalized internal force, generalized inertia force, and generalized external force. The notation \((\bullet)_{\text{q}}\) denotes the derivative relative to generalized coordinates \(\mathbf{q}\). Note that the position vector \(\mathbf{x}_{\text{B}}\) is calculated from interpolation, giving:  
\[\mathbf{x}_{\text{B}}=N_{\text{I}}x_{\text{B}I} \tag{53}\]  
The subscript \(I\) is determined by the interpolation function. It is remarked that for the case of sliding along very flexible beam, the IGA method takes more advantage on accuracy and efficiency when using the \(C^{k}\)-continuous NURBS basic function (k \(\geqq 1\)).