Modelling Friction-Induced Dynamic Instability Dedicated for Isogeometric Formulation

,


Introduction
Friction naturally occurs in sliding contact.Typically, as a result of friction, the kinetic energy is largely dissipated as heat, while some of the energy is conserved as vibration which as a consequence can lead to acoustic noise [1].Te resulting noise generated by friction is more complex to model mathematically and cannot be typically explained by a single mechanism [2][3][4].Te complexity can be largely owed to modelling the contact interface properties and the mechanism of friction at various scales.Te sound produced from friction can be purely as a result of unsteady excitation or by reaching a steady state through self-excitation [5].Squeal noise is largely identifed as the latter case, where as a result of perturbations around equilibrium, the system reaches a new steady-state characterizing self-excitation behavior.Tis can be understood as coalescence of modes where two modes exist at a same frequency, leading to selfexcitation between the modes under favourable conditions in the presence of friction.Even though friction is commonly associated with damping, in the case of sliding contact, it is understood that the interaction as a result of friction between the normal forces and the tangential forces relative to a contact surface gives rise to self-excitation behaviour which characterizes futter instability.Tis is evident in the literature even from early analyses based on lumped models, which greatly improved the understanding of futter-type instabilities induced by friction [6,7].
Flutter instability induced by friction can be inferred through transient analysis with the presence of a limit cycle [8] which is nevertheless expensive to model given the nonsmooth nature of friction and to capture the high frequency characteristics over human audible range [9].One of the most common ways to circumvent the expensive computation of such problems is through analyzing the perturbations around the equilibrium state of the dynamic problem.Te equilibrium state of sliding contact problem with friction can be typically characterized as steady-sliding or quasistatic equilibrium depending on the nature of the external forces present [10,11].With some assumptions made over the nature of the perturbations such that the dynamics of the system for a perturbation can be expressed linearly, the system can be analysed through eigenvalue analysis.It should be noted that such linear assumptions can only characterize instabilities very close to the equilibrium.Te presence of the limit cycle is inferred through complex conjugate pairs of eigenvalues R(λ) ± I(λ) or its equivalent form ± R(λ) + I(λ) depending on the assumed solution, and hence, such analyses are termed as complex eigenvalue analysis.A complex conjugate pair of eigenvalues only characterizes divergence from the equilibrium to reach a limit cycle but does not give any information on the actual characteristics of the limit cycle [12].Nevertheless, CEA has been widely accepted in industrial applications such as for designing of automotive and aircraft brakes to mainly predict the presence of limit cycle which can in turn indicate the presence of squeal phenomenon and also whirl phenomenon in systems such as aircraft brakes.
Even though fnite element methods are used in discretizing a complex domain, the contact interface is typically modelled through a set of linear or nonlinear springs which can be regarded as the collocation method.Te sensitivity of such an approximation made over the continuum contact interface in the estimation of steady-sliding equilibrium or the prediction of instabilities through CEA is largely not studied.In this paper, we give a simple way for approximating the integrals defned over the contact interface, the variational form of contact and friction, which is similar to Gauss point to surface (GPTS) approach [13,14] regarded as a mortar based method [15] without domain decomposition.It is intuitive that such an approximation without domain decomposition for dissimilar meshes at the contact interface is not empirically accurate.But the accuracy required depends on the quantity of interest, where in our case the sensitivity of steady-sliding equilibrium and CEA to the approximation.We focus on the variational approach without domain decomposition since it provides the best balance between computational cost and accuracy, in comparison to domain decomposition approaches which are also more complex to implement.Te variational formulation is also satisfed through collocation for comparison to show the efcacy of the GPTS approach for CEA.We do not focus on the linearization part typically associated with large sliding contact problems since to compute steady-sliding equilibrium and the linearization to be made around the equilibrium for CEA requires no defnition of velocity terms.We expand the defnitions based on the normal compliance approach [7,16] where the unknown contact boundary conditions of normal stress are defned to be displacement dependent.Tis can be regarded as the classical mesh-tying constraint enforced through the penalty approach, where the term mesh-tying expresses the coupling of the displacement feld between the domains at the contact interface.Further, we consider the isogeometric approach [17] for discretization owing to the accuracy provided in structural dynamics application [18] and contact problems [19] where the higherorder continuity of spline basis functions could be exploited [20][21][22].
Te meaning of penalty term in classical sense is essentially to satisfy the contact constraints, where an ideal penalization is of order infnity.But in a numerical context, penalization is taken to be a high value to enforce the constraints, considering the condition number of the matrices.Tis could be unrealistic for a dynamical system, where the dynamical response is sensitive to stifness at the interface which is essentially the normal compliance approach models.Tis is in contrast to the zero compliance model of Signorini conditions which the penalty approach essentially tries to satisfy.Tis means that the discretization scheme should also be accurate at lower values of penalization which could refect the contact stifness to be modelled with the normal compliance approach.Tis is because for some formulations, the accuracy of the normal stress could be improved considerably with a large value of penalization, given that the condition number of the matrices does not make it numerically unstable [19,23].Hence, we also focus on the accuracy of normal stress at low values of contact stifness which is classically ignored.
We also make some inference on numerical characteristics of CEA in relation to contact stifness and the accuracy of modelling normal stress at the contact interface, such that the sensitivity of CEA could also be implied.Tis is important in approximating the integrals defned over a contact domain such that it characterizes the contact interface in the interest of CEA independent or with little inference from experiments.Tis could be useful in applications such as optimization of systems for friction-induced dynamic instability through numerical simulations.Tis is mainly evident in the shape optimization of the contact interface itself [24,25] since the contact interface has to be modelled numerically independent of experiments for a given shape variation in optimization, given the necessary experimental inference initially, parameters inferred to model normal compliance for example.

Modelling
We give a general overview of continuum formulations in this section which is useful to defne the variational formulation of contact and friction, for which the discretization will be expanded in the following section.We start with the dynamical problem of contact and friction considering contact between an elastic domain and a rigid domain for the sake of simplicity §.From this, we defne the steady-sliding equilibrium problem § and the problem of the dynamics of perturbation around the equilibrium §.Followed by, we 2 Shock and Vibration defne the problem for the dynamics of perturbation explicitly §, where we defne explicitly the quantities of interest required to be calculated from steady-sliding equilibrium.Tis is to understand the approximation of the integrals defned over a contact domain for given quantities of interest which are to be defned explicitly in the interest of studying contact formulations given in §.Further, for the dynamics of the perturbation problem defned explicitly, contact is considered between two elastic domains so as to give clarity in the expansion of the contact formulation in §.Finally, we express the discretized form of the problems, also in matrix form for all the preceding continuum problems defned §.

Initial-Boundary Value Problem.
Te initial-boundary value problem for a domain in contact with a rigid body in the presence of friction can be given as where u: with zΩ defning the boundary of Ω, and  v n defning the unit normal vector on zΩ.Under isotropic material consideration, the constitutive equations can be defned as ) are 3D Lamé parameters expressed in terms of Young's modulus E and Poisson's ratio ].Te kinematic relation for the strain tensor ε under infnitesimal displacement is defned as ε � 1/2(∇u + ∇u T ).For a system with domains Ω (a) and , where X represents the material coordinates.⃖ X (b) on Ω (b) is defned as the outward normal projection  v n from zΩ (a) .
Signorini conditions (1b) and Coulomb's law (1c) defne elegant mathematical models for contact and friction, respectively, in macroscopic view.But such conditions have been found to lack the characteristics of contact interface in reality [7].Tis leads to the view of the normal compliance approach [16] to defne an approximation of the characteristics at a contact interface, also through which regularization for the multivalued mapping of (1b) can be achieved.We do not focus on the regularization of (1c) which is not useful in the scope of steady-sliding equilibrium and with the nature of the perturbations to be considered.Normal compliance is defned by σ n as the function of relative normal displacement (b) � 0 in this case) at the interface in relation to the initial gap g n , given as where (.) + allows only a positive value.Tis can be extended to friction as where the parameters c n , m n , c t , and m t characterize interface properties and are to be determined from experiments.We consider the friction model purely based on Coulomb's law to focus on numerical modelling, and hence, c t � μc n and m t � m n in this case.With Γ C known as the active set which satisfes the normal compliance constraint (.) + at a given time, the problem can be expressed as variational inequality [26], given as where δu is the directional derivative.Te solution to the above dynamical problem is often discussed in the context of nonsmooth mechanics [27,28] which we do not focus here.But through regularization of the nondiferentiable friction functional, mainly as ‖ _ u t ‖ ⟶ ψ(‖ _ u t ‖), the abovementioned inequality can be expressed as variational equation [7,11].In this case, ψ corresponds to a continuous function over velocity for transition from stick to slip state.Given the initial conditions u 0 and _ u 0 which satisfy the inequality at the initial time, the system can be numerically integrated with time.
We focus on modelling futter-type dynamic instability through classical theories of linear analysis, where the frstorder efect of perturbations around a fxed point is analyzed.Hence, the stability of the dynamical system with frictional contact can be characterized by determining the Shock and Vibration fxed points which are typically of quasistatic or steadysliding equilibrium depending on the external forces and defning the dynamics for the perturbation around such fxed points.We only focus on steady-sliding equilibrium for the following discussions.

Steady-Sliding
Equilibrium.Te dynamical system can be expressed as steady-sliding equilibrium when no net acceleration is present in frictional contact.Tis can be seen as a series of unchanging equilibrium states with respect to time, where the equilibrium characteristics remain the same except for the relative velocity at the contact interface.Nevertheless, the problem can be expressed purely by static forces.Te steady-sliding equilibrium explicitly characterizes slip condition at constant sliding velocity, and hence, at the equilibrium, σ t � μσ n on Γ C .Further, empirically, the knowledge of sliding direction  v k at any point on Γ C is defned to be known a priori, and hence, σ t � μσ n  v k .Tis means that the nondiferential friction term of (4) as a result of the inequality (3) can be expressed as variational equation at pure slip state (Te normal compliance terms in modelling σ n make it nonlinear.But eventually through linearization of the normal compliance model, bilinearity could be achieved which will be discussed in the following).Since the equilibrium characteristics remain the same for all the time, the equilibrium state could be expressed with no time-dependent forces but purely by static forces.Te problem of steady-sliding equilibrium can hence be expressed as Te nonlinear problem can be solved through the Newton-Rhapson method to obtain the steady-sliding solution u.It should be noted that for steady-sliding case with friction, nonunique solutions can exist [10,29].

Initial-Boundary
Value Problem for Perturbation.Te perturbation of the displacement feld relative to the steadysliding equilibrium can be given as where  u(t) and t correspond to perturbed displacement feld and time, respectively.Te dynamics of the perturbation can be deduced by substituting u +  u in the variational equation form of (4) (Te variational equation can be deduced by regularization of the nondiferentiable friction functional of the variational inequality (4).For the sake of simplicity, we do not focus on regularization since at the steady-sliding equilibrium and for perturbations around the equilibrium, the system is considered to operate beyond the discontinuity of Coulomb's friction law, and hence, regularization has no efect) and subtracting with (5) to obtain Te idea is to analyze the dynamics of the system for a small perturbation such that the stability of the dynamical system can be characterized through the evolution of the dynamics for perturbation close to the steady-sliding equilibrium, where we hypothesize that the evolution of the dynamics can be sufciently characterized linearly.Tis brings the question of linearizing normal compliance.
We start with the discussion of linearizing the operator (.) + under necessary conditions.Te perturbation  u n on Γ C may not result in the same active set Γ C , which introduces nonlinearities defned by (.) + (Separation on Γ C defnes u n − g n < 0 ⟹ p (u n − g n ) m n + ⟶ 0 which can only be taken into account nonlinearly.Tis is also true for the case contrary to separation).Tis can be linearized with the hypothesis that the onset of instability to model occurs close to the equilibrium u such that the perturbation  u is small to defne Γ C as stationary.But this may not be true since Γ C is stationary only when u n +  u n − g n > 0 (in the view of normal compliance, the perturbation  u n for no variation of Γ C defnes the variation of σ n .Tis would have been impossible with Signorini conditions which can be regarded as the zero compliance model, where  u n > 0 leads to strict contact separation on Γ C ), implying that σ n ≪ 0 must be true on Γ C at steady-sliding equilibrium, conceivably as a result of external forces.Tis can lead to a new active set Γ C in the current confguration X + u(X) that satisfes u n (X) +  u n (X + u(X)) − g n (X) > 0 on zΩ, which depends on the nature of the perturbation  u n (X + u(X)) and σ n (X) at steady-sliding equilibrium.Perhaps, this could reveal if the linearization of the operator (.) + is realistic for a given system.To simplify, we consider Γ C as stationary for the perturbation  u, which evades 4 Shock and Vibration the operator (.) + .With further linearization for the polynomial terms of the normal compliance, the frst-order efect of the perturbation  u close to u can be expressed as An illustration of the linearization is shown in Figure 1.It is also possible that the nonlinearity defned by the normal compliance parameters could mean that the linearization at very low values of u n − g n could yield very low m n c n (u n − g n ) m n − 1 .Hence, the assumption of stationary Γ C could be justifed in this case.
We provide discussions on the presumed characteristics of friction for the perturbation  u to express the dynamics of the perturbation to be linear.Even though steady equilibrium is expressed by static forces, for the perturbation  u(t), the state of friction can change from slip state of steady-sliding equilibrium, where without regularization for friction, the dynamics of the perturbation resembles the variational inequality of equation ( 4).Hence, this requires further assumptions on the nature of the perturbation  u such that the system is always in slip state, which is reasonable for small perturbations and also preserves the variational equation form.For the given example, we consider contact between an elastic domain and a rigid domain, and in this case, the assumption of stationary Γ C at sliding is fairly easy to consider on the elastic domain.With the dynamics considering two elastic domains in contact, for steady-sliding equilibrium, the forces at the contact interface are stationary, and hence, Γ C can also be taken to be stationary even in the presence of relative motion.But for the perturbation  u(t) of such problems, the variation of forces on Γ C at slip state means that the assumption of stationary Γ C may not be realistic.Tis is true unless _  u(t). v k for the dynamical system (4) is negligibly relative to the natural frequencies of the system, which is taken to be the case here.Tis is evident in applications such as brake system, where the rotational frequency of brake discs is much lower than the natural frequencies of the brake system.
We defne p(X) which can be interpreted as contact stifness for perturbations close to steady-sliding equilibrium.For the following discussions, we do not relate to the experimental determination of the parameters c n and m n .We mainly focus on the fnite element discretization of the integrals defned over Γ C for the dynamics of the perturbation and the consequence of p(X) on the discretization schemes.Hence, the focus on the infuence of discretization schemes in the estimation of steady-sliding equilibrium is not considered in this study.u efectively defnes Γ C at steady-sliding equilibrium, also at which normal compliance is linearized to defne p(X).But Γ C and p(X) can be defned explicitly, for which a new initial-boundary value problem ( 9) can be considered.Tis is to simplify the comparison of discretization schemes for (17) with constant Γ C and p(X), rather than to deduce from steady-sliding equilibrium where Γ C and p(X) could indeed vary depending on the discretization scheme.

Initial-Boundary Value
Problem for Perturbation: Explicitly Defned.Te preceding problem (1a) was expressed for contact between an elastic domain and a rigid body to simplify the expressions in focus on the nature of the problem.But to generalize the following problem for contact between elastic bodies, the subscript k is introduced to distinguish the domains in contact.With Γ C and p(X) known a priori, either from steady-sliding equilibrium or defned explicitly, the strong form of the dynamics for  u can be expressed as Te parameters p(X) and μ are defned to be constant on Γ C between all the domains in contact.Te variational formulation of the problem can hence be expressed as Possible separation at low contact stress Figure 1: Linearization of contact defned by the normal compliance approach for the perturbation  u at the steady-sliding equilibrium u.

Shock and Vibration
Te traction force t C can be decomposed as Hence, with the expansion of the traction forces for n k domains in contact, equation (10) can be expressed as To simplify, we consider contact between two domains Ω (a) and Ω (b) , with the derivation of traction forces on Ω (a) .
Hence, the inner products C can be expressed as Te traction forces on Ω (b) can similarly be defned from the conservation of momentum as σ (a)  n � −σ (b)  n and σ (a) t � −σ (b) t .

Finite Element Approximation. We introduce the function space
in which we seek the solution u, where the Sobolev space 3 , where (H − 1/2 (Γ C )) 3  is the dual of the space (H 1/2 (Γ C )) 3 .
Te solution u ∈ V can be expressed as With the fnite element approach, we defne a fnite-dimensional function space h V ⊂ V, and hence, there is some bound on v i ∈ h V. Te approximation of u in h V can be expressed as u ≈ h u ∈ h V. We focus on the defnition of the space h V with the isogeometric approach in the next section.But in general, we give the matrix form of the variational formulations derived with h V.
Te fnite element approximation of u ∈ h V for equation ( 5) considering the two domains Ω (a) and Ω (b) in contact can be expressed as . v n .∼ k defnes the domains which are not k, i.e., in the case of contact between Ω (a) and Ω (b) , k � a implies ∼ k � b and vice-versa.h g n implies the realization of g h with the parameterization of domains.Te abovementioned expression can be expressed in the matrix form as 6 Shock and Vibration where K is the stifness matrix, F C and F F are the nonlinear traction force vectors of contact and friction, respectively, and with F being the classical force vector implying external forces.Te linearization of the integrals defned over Γ C in solving the steady-sliding equilibrium (5) through the Newton-Rhapson method takes the similar form of. Tis means that the tangent matrices of and , respectively, to be discussed.Similarly, with the substitution of equation in equation (7), for  u ∈ h V, results in the following form: which can be expressed in the matrix form as Te properties of the matrices can be given as follows, where M and K are symmetric and positive defnite matrices.K C is also symmetric which essentially defnes the conservation of momentum at the interface for normal contact.While K F is nonsymmetric which defnes the nonconservative nature of friction at slip state.Hence, all the contact forces are expressed as displacement-dependent forces linearly.Considering a solution of the form Θe λt for  U, the characteristics of the eigenvalue problem which we refer to here as CEA can be expressed as It is evident that λ determines the state of the perturbed solution  U. Given an eigenvalue of the form λ � R(λ) + I(λ), the I(λ) part models the oscillatory behaviour while the stability of the given mode is defned by R(λ).R(λ) > 0 can be understood to characterize an unstable mode and defnes the divergent nature of the solution.Similarly, R(λ) < 0 can be understood to characterize a stable mode which can be interpreted as damping.Te occurrence of a pair of eigenvalues ± R(λ) + I(λ) defnes the presence of a limit cycle which is of interest for applications concerning futter-type instability induced by friction.Damping can also be factored through the defnition of a matrix C, typically through modal or Rayleigh damping.In the scope of modelling rotational inertia efects, gyroscopic matrix G can also be defned, where C and G are both velocity dependent.We factor out damping and gyroscopic efects in our model to mainly focus on contact and friction modelling.

Isogeometric Approach for Discretization of the Contact
Integrals.We defne the space h V through the isogeometric approach to take advantage of the properties of the spline basis functions.Te higher-order continuity of the spline basis functions between the knots provides superior approximation properties.Te infuence of the continuity in contact problems has also been largely studied for various formulations in the context of penalization or mixed formulation approach to enforce the contact constraint [19].Tis is mainly true for dissimilar meshes between the domains in contact at the contact interface, where simple discretization schemes through collocation were shown to provide improved accuracy.Te higher-order solution continuity is a consequence of the geometric continuity provided by spline parameterization.Te geometric continuity greater than c 0 is also useful to defne the actual normal of the contact surface and also avoids numerical discrepancies in large displacement contact problems, as a result of the faceted c 0 continuity of classical fnite element meshes.

NURBS for Approximation. B-spline basis functions can be defned by Cox de Boor's formula as follows:
where p is defned recursively for p > 0 to obtain a curve of degree p, which starts with a piecewise constant at p � 0. Naturally, a uniform knot vector can be defned as ξ � ξ 1 , ξ 2 , . . ., ξ n+p+1  , where any ξ i − ξ i+1 is uniformly spaced.Te knot vector need not be equidistant, and the multiplicity of a knot ξ i by M in the knot vector decreases the continuity by c p− M across the knot ξ i , which defnes a nonuniform knot vector.Trough B-spline basis functions and a knot vector ξ � ξ 1 , . . ., ξ n+p+1  , a B-spline curve can be defned with the basis functions and its coefcients as follows: Te coefcients P i ∈ R d are the control points, with d being the dimension of the space.Te defnition of a weighing parameter w i > 0 associated with a respective basis Shock and Vibration function N i , normalized defnes rational B-splines where it respects the partition of unity, given as follows: Te parameter w i provides a new dimension for controlling the geometry through projective transformation, while the afne transformation is achieved by P i .Hence, the combination of nonuniform knot vectors and rational basis functions defnes NURBS.Te higher dimensional NURBS are a natural extension of its 1-dimensional precursor through tensor product defnition where the order of the tensor is the same as the dimension of the geometry.Hence, a tensor product NURBS surface can be defned as with n × m net of control points P i,j .Similarly, a tensor product NURBS volume can be defned as where the knot vectors are given as Te abovementioned expression can be simply expressed in the matrix form as X v (Ξ) � R(Ξ)P.We also give the following notations, where to express the same degree for all the knot vectors defning a NURBS parameterization in R 3 is simply defned as D x , i.e., (p � q � r � x): � D x .Similarly, for the same continuity y between the internal knots in all the knot directions is expressed as C y .
Te parameterization of a domain Ω ∈ R 3 as an initial geometric description through NURBS can be expressed as where X defnes the mapping from the parametric domain  Ω to the physical domain Ω, and for simplicity, we consider the parameterization of the domain Ω through only a single patch: Te analysis-suitable parameterization X 6 (for simplicity of the notation, we defne X to be the default notation for analysis-suitable parameterization of a domain Ω ∈ R 3 ) can be achieved through the refnement of � X ⟶ X with one or several of the refnement methods: h, p , and k, where X can be defned as X v (Ξ) � R(Ξ)P to take into account of the modifed knot vectors and control points.Te isogeometric approach for the approximation of a solution is achieved through the bases R i,j,k , where for the vector-valued function space h V, the vectorial defnition of the bases R i,j,k ∈ R 3 can be defned as Tis can be expressed in the matrix form as R i,j,k (Ξ) � R i,j,k (Ξ)I, with I being the identity matrix.Hence, , In an abstract sense, the bases R i,j,k (Ξ) in parametric space are transformed into the bases ϕ i,j,k (x, y, z) in physical space using the push-forward operator °, where the bases ϕ i are defned with the property Hence, the approximation of a feld variable on Ω is defned through all the bases ϕ i spanning the fnite-dimensional function space Φ.For the isogeometric approach, we express the fnite-dimensional space h V as Φ and its associated bases v i as ϕ i .Te approximation of  u∈ Φ can be defned as h  u �  ∀i∈Ω ϕ i  u i , expressed in the matrix form as h  u � N(X)  U, where Te variational form of contact and friction given in equations ( 13) and ( 14) expressed in the function space Φ can be defned as Te expansion of the abovementioned expressions also generalizes for the defnition of the tangent matrices for the nonlinear terms of contact and friction in equation (15).Te parameterization of the domains Ω (a) and Ω (b) with NURBS can be expressed as X (a) � R (a) (Ξ (a) ) and X (b) � R (b) (Ξ (b) ).

Te active sets Γ (a)
C and Γ (b) C is defned with the property g n � 0: X (a) .v n � ⃖ X (b) .v n , where  v n in this case is taken to be the outward normal projection from Shock and Vibration projection exists that maps X (a) on Γ (b) C as ⃖ X (b) .For the following explanations, we detail the derivation of traction forces on Γ (a) C , which can be extended to Γ (b) C depending on the considered discretization scheme.
were  v: Essentially, the abovementioned equations are a construct of the matrix form as For k � l � a, the defnition of the integrals is straight forward.In contrast, for k � b and l � a, we have terms of the form  Γ (a) r dΓ (a)  C which can be expanded as where the integrals are simultaneously defned over the bases of the two domains in contact since C ).Even though the defnition of integral for ϕ (b)  j ( ⃖ X (b) ) on Γ (a) C exists through the projection ⃖ X (b) (X (a) ), for dissimilar meshes at the contact interface, the defnition of numerical quadrature scheme for the integrals of such form requires domain decomposition to fnd the common span: Tis means that the integral exists between ϕ (a) i (X (a) ) and ϕ (b) j ( ⃖ X (b) ) only on the span where the projection ⃖ X (b) (X (a) ) exists and hence requires a quadrature scheme specifc on the span.
Inferring from the variational form of equations ( 29) and (30) with the expansion provided in equation (30), the following relation should hold: where it verifes the conservation of linear momentum at the contact interface.Since we deal with contact between fat surfaces, for the normal compliance approach which can be viewed as coupling of displacement feld between the contact surfaces, conservation of angular momentum is implicitly satisfed.
We give some intuition through the collocation approach to satisfy the conservation of momentum at the Shock and Vibration contact interface and extend it to the defnition of quadrature scheme without domain decomposition.For collocation, the integral  Γ C ϕ i dΓ C on one of the domains, in this case, Γ (a) C -can be defned as i , where i (a) is the set of points on Γ (a)  C which depends on the collocation scheme [21,30].It should be noted that this is in contrast to the collocation of the strong form.Hence, the relation ( 32) satisfed discretely at the collocation points takes the form as follows: where i ϕ :� i ϕ I, i ϕ (a) :� ϕ (a) ( i X (a) ), and i ϕ (b) :� ϕ (b) ( ⃖ X (b)i X (a) )).Tis satisfes conservation of momentum at the contact interface, even though the integral  Γ (a) C cannot be defned accurately.For any i, the following relation holds: Tis means that any quantity defned at a collocation point i over i ϕ (a)   i is projected equally over the bases in Tis is as a result of the properties of the NURBS basis functions similar to classical fnite element basis functions which satisfy partition of unity.Since the integral is satisfed only at discrete points, the solution may not be to the necessary accuracy.Te collocation strategy can be replaced by a numerical quadrature scheme as i where i (a) , in this case, corresponds to the quadrature points with i w being the quadrature weights.But the notion of i w on C ). Tis is where the higher-order continuity of the spline basis functions and increasing the number of quadrature points can be useful, which will be discussed in the next section.Hence, equations ( 29) and ( 30) approximated through a quadrature scheme can be expressed as i w p i N (a) .v n   T i N (a) .v n   i N (a) .v n i w μp i N (a) .v t   T i N (a) .v n   i N (a) .v k where and i N: Te integral defned over the parametric space  Ω for the isogeometric approach takes the form.
where i R: � R( i Ξ) and i Ξ being the collocation point in the parametric space and J is the Jacobian matrix for the mapping Ξ ⟶ X. i Ξ (b) is the corresponding map of ⃖ X (b) in the parametric space, which can be determined through the Newton-Rhapson method in solving for X (b) )).Hence, there exists a mapping i Ξ (a) ⟶ i Ξ (b) for which X (a) ( i Ξ (a) ) � X (b) ( i Ξ (b) ).
From the conservation of momentum at the interface, the following relation holds σ (a)  n � −σ (b) n and σ (a) t � −σ (b) t .Hence, the traction stresses on Γ (b)  C can similarly be defned as 10 Shock and Vibration It should be noted that for the abovementioned equations, even though the traction forces are defned over Γ (b)  C as 〈 h σ (b) , ϕ (b)  i 〉 Γ (b) , the quadrature points' set i (a) is determined only on Γ (a) C , where its corresponding projection on Γ (b) C is defned through the projection ⃖ X (b) .Tis also includes Jacobian | i J (a) | which is evaluated only on Γ (a) C , similar to equations ( 36) and (37).Tis greatly simplifes the fnite element discretization of the integrals which may otherwise require the domain decomposition approach to satisfy the integrals exactly, where this as a consequence can have an efect on accuracy which is discussed in §.With the abovementioned defnitions, the matrices K (a−b) C and K (a−b) F for the system take the form.
Since the quadrature rule or collocation is defned with respect to only one surface, it is commonly termed as halfpass.Te role of Γ (a)  C and Γ (b) C can be switched, and the average over the two half-pass can be taken into account to defne the so called full-pass [31].A variation of the full-pass approach is also given in [32] based on surface potentials.
Considering equation (30) given in the general form, the submatrices of the form K (.)  C and K (.) F , in contrast to the form K (.,.) C and K (.,.)  F , of the system matrices K (a−b)

C
and K (a−b) F , respectively, correspond to the case k � l.With the relation (34), K (.)  C and K (.)  F can be lumped.Hence, with equation ( 30) for l � k, lumping can be achieved by defning ϕ (k)  j ≡ δ i,j ,δ i,j being the Kronecker delta product, where the inner product ϕ (k) j , ϕ (l�k) i   is expressed with the property.
l ≠ k defnes the matrices of the form K (.,.) C and K (.,.) F , where the relation (34) holds.Tis preserves the conservation of momentum at the interface.
Collocation can be achieved with the weak form considering i w � 1, | i J| � 1, and i ∈ I being a suitable set of collocation points for equations ( 36)-(39).It should be noted that for collocation schemes based on the variational form, the efect of the area for a collocation point can be taken into account through i w and | i J| to improve the accuracy [30].It should be noted that as with the classical collocation method, collocation with the strong form can also be considered, where the strong form of contact and friction boundary conditions are considered, as in equation ( 9), which is shown to provide more accuracy for the basis functions with higher-order continuity.For the isogeometric approach, Greville and Botella points are of interest which exploit the properties of NURBS.Also, for a given NURBS parameterization, the number of Greville and Botella points is fxed which is in contrast to the Gauss quadrature scheme which can be defned independent of the NURBS parameterization.Nevertheless, the main interest with the collocation scheme is the reduction of computational cost associated with the reduced number of evaluations, where the number of Greville or Botella points is typically very low in comparison to the number of quadrature points for a reasonable accuracy.We only consider Greville points in the following section.Greville abscissae in a direction are the average of the knots in the knot vector for the direction.Further, the Greville points lie on the surface which makes it useful to defne the projection ⃖ X, and the number of Greville points is the same as the number of control points.

Results and Discussion
Te variational formulation satisfed through collocation with Greville points is simply termed as the collocation method, unless otherwise specifed.For the following discussions, we use the term variational formulation to specify the use of quadrature scheme in an approximate sense as discussed in the preceding section, as it preserves the integral sense of the variational formulation.Te results are shown with only half-pass formulation.It should be noted that the results implied with half-pass for a given formulation cannot be extended to full-pass.Since even if the patch test could be passed, it does not guarantee conditions for LBB stable [33] or the implications in CEA.Hence, for a given formulation, the full pass must be tested to imply its stability in CEA which is not considered.D x and C y are considered to be the same for all the domains in a given analysis.For the variational formulation, we defne the number of Gauss abscissae between the knot spans in a direction to be the same as the order of the NURBS in that direction, unless otherwise specifed.Te results are mainly focused on the accuracy of the formulations inferred through a contact patch test, and the inferences are related to the sensitivity of CEA for a given formulation.Even though we only focus on the accuracy of σ n through the considered contact patch test, it should be Shock and Vibration noted that σ t for steady-sliding equilibrium is given by the relation |σ t | � μ|σ n |.Hence, the accuracy of σ n will have a subsequent efect on σ t , though the sensitivity of σ n and σ t on steady-sliding equilibrium is not considered in this study.Further, all the results of CEA are obtained with μ � 0.7.To proceed with the following discussions, we introduce the global contact stifness P at the interface Γ C as P �  Γ C p(X) dΓ C .Hence, a given P is constant at the interface, and the local contact stifness p(X) varies depending on the accuracy of the formulations.
For the contact patch test, we consider p(X) to be constant on Γ C , and hence, an ideal discretization should yield p(X) � P/Area(Γ C ).But for a given P, the local contact stifness p(X) can vary depending on the approximation in obtaining a formulation.Te discretization of p(X) has direct infuence on contact normal stress σ n through which the accuracy of a discretization scheme can be inferred to some extent since for some formulations, the accuracy depends on the condition number of the matrices where higher accuracy can be attained with a large value of P. For this reason, we focus on considerably low values of P in the interest of the normal compliance approach.Te setting for the contact patch test is shown in Figure 2, where we consider two cuboid domains in contact with two diferent meshes visible through the knot vectors of the NURBS parameterization.Te bottom face of Ω (b) is fxed, while the top face of Ω (a) is applied a uniform pressure of 1000 N/m 2 and constrained to move normal to Γ (b)  C .Te material parameters E � 10 5 and ] � 0 are defned to be the same for both the domains.Analytically, one can infer that |σ n | � 1000 N/m 2 on Γ C , which provides good benchmark to compare σ n calculated with diferent formulations.Te contact patch test can be defned by considering equation (15) for μ � 0 and  v k � 0. With further simplifcation of c n � p(X) and m n � 1, the matrix form of the problem can be solved for K To begin with, we show the comparison of variational formulation with lumping and without lumping in Figure 3. Surprisingly, the approach with lumping has reduced oscillation of σ n for a given P, where for the higher value of P � 10 14 , the approach even passes the contact patch test to machine precision.For CEA 19, we consider a simple discpad system which can give rise to dynamic instability induced by friction.Similar to the contact patch test, p(X) is defned explicitly to be constant on Γ C , as discussed in §.Tis is because p(X) determined through steady-sliding equilibrium can vary depending on the discretization scheme to defne an unbiased comparison of CEA results.Even though the lumped approach is observed to be superior in the patch test, it proves to be numerically unstable for CEA, shown in Figures 4 and 5.We defne the following notations u (d) and u (p) to express the displacement feld of the disc and the pad, respectively.Even though the scale of the displacement feld has no sense for the computed eigenvectors, comparison can be made relatively with the displacement feld for a given eigenvector.It should be noted that we deal with complex eigenvectors where the imaginary part of an eigenvector characterizes phase diference of the displacement feld, which will not be considered here.Te mode shapes in Figure 4 considered for comparison are chosen to be close in frequency and has the same characteristics in relation to the mode shape of the disc.For the lumped approach, it can be observed that the relative displacement at the contact interface of the disc is much lower in relation to the displacement feld elsewhere.Te arrow on which u (d) .v (d)  n is plotted also traverses the boundary of the contact interface zΓ (d)  C (Boundary of the contact interface can be defned as zΓ C for the contact interface Γ C ) where strong solution gradient for u (d) .v (d)  n can be observed.Tis is in contrast to the approach without lumping where u (d) .v (d)   n is smooth across zΓ (d)  C .In this case, the dynamic response of the approach with lumping and without lumping is observed to be diferent as shown with the results in Figure 6 where the natural frequencies I(λ) and its corresponding instabilities R(λ) largely do not coincide.Empirically, the behaviour observed at the contact interface for the lumped approach seems unrealistic since Γ (d)  C and Γ (p) C are observed to be locked with each other relative to the rest of the domains in contact.Tis could also be a consequence of the properties of NURBS basis functions and hence can be verifed with classical fnite element basis functions, which is not considered here.It should be noted that the role of mesh refnement at zΓ (d)  C plays an important role in capturing the solution gradient since the higher-order continuity of the spline basis functions will have strong propensity to defne a smooth solution across zΓ (d)  C .In our case, multipatch parameterization of disc is considered where one of the patches consist of Γ (d)  C .Hence, across zΓ (d) C , h− refnement was achieved in the patch that contains Γ (d)  C such that the knots can trace along side zΓ C , strategies based on RBQ [34], or in the scope of local refnement T-splines [35] or THB splines [36] have to be considered.For the lack of understanding of the numerical implication with lumping, we only focus on the approach without lumping for the following discussions.Nevertheless, the above discussions reveal that the contact patch test cannot be purely accounted for the accuracy in CEA.
As it can be inferred from the contact patch test in Figure 3, the accuracy of σ n can be improved with increase in P. For the contact problems defned by static analysis, when the contact constraint is satisfed with the penalty approach, it is necessary that p is large where in an ideal case p ⟶ ∞.But for dynamical systems, contact stifness models the interface properties which have infuence on the resulting dynamical response.It has been observed that the value of contact stifness in this case is not typically higher than the penalty value typically used in satisfying the contact constraint [37].In this case, the variational formulation with higher-order continuity and increased number of Gauss quadrature points can yield better results.Te necessary accuracy indeed depends on the sensitivity of the approximation in estimating the steady-sliding equilibrium, the efect it has on linearization for CEA, and the evaluation of R(λ) and I(λ) with CEA.In the interest of computational cost, as the order of continuity C x− 1 is increased (k Shock and Vibration refnement) with increase in degree D x (p refnement), the increase in the number of control points is of very low order.Even though the increase in the degrees of freedom of the fnite element matrices is small for simultaneous p and k refnement, this can reduce the sparsity of the matrices and hence can increase the computational cost in terms of both memory and time.Hence, C x− 1 can be considered only up to a reasonable value of D x depending on the required accuracy.
With the following results, we analyze the signifcance of the higher-order continuity of the spline basis functions and number of Gauss quadrature points on the accuracy of σ n .
Te infuence of the higher-order continuity is shown in Figure 7 for various D x and C y considering P � 10 6 .It can be observed that the oscillation in σ n decreases with increase in continuity.Tis can be clearly seen with comparing the cases of D 2 C 0 and D 2 C 1 defned with equal number of Gauss points between the knot spans ∈ R 2 (span ∈ R 2 : [ξ i , ξ i+1 ] × [η j , η j+1 ], considering that Γ C is defned on the surface parameterized by knot vectors ξ and η), where D 2 C 1 shows reduced oscillation of σ n .While for D 2 C 0 , the solution is observed to be only c 0 continuous at the knots.Te results of D 2 C 0 and D 2 C 1 are shown in comparison to D 4 C 3 for which the oscillation of σ n is negligible.Te accuracy of σ n can be 14 Shock and Vibration improved with increasing the number of Gauss quadrature points, which is shown in Figure 8.For the case of D 4 C 3 in Figure 9, the number of Gauss points is also increased with the degree by default as specifed before, and hence, the addition of Gauss points also contributed to increase in accuracy along with the contribution from higher-order continuity.Nevertheless, even D 4 C 3 does not satisfy the contact patch test at P � 10 , shown in Figure 9.It can be inferred that even though the accuracy of σ n increases profoundly with initial refnement of continuity and increasing the number of Gauss quadrature points, after a certain order of continuity and number of Gauss points, the accuracy can only get up to a certain digit.Hence, from the asymptotic convergence of σ n for increase in continuity and number of Gauss points, it can be inferred that the contact patch test cannot be passed to machine precision.It can also be observed that the oscillation of σ (b)  n is greater than σ (a)  n for D x C y .Tis can be hypothesized as resulting from the quadrature scheme defned only over C , i.e., the quadrature points, Jacobian, and the quadrature weights are evaluated only on Γ (a) C , while only their corresponding projection is defned over Γ (b)  C .Finally, we provide some intuition for possibly why the accuracy improves with increase in continuity and number of quadrature points for the variational formulation.As discussed in the preceding section, the approximation of the integrals defned over Γ C is simplifed with considering the quadrature scheme on only one of the domains in contact, while its projection is defned on the other.Tis is not very intuitive since the quadrature quantities such as quadrature weights i w and the determinant of Jacobian | i J| of the quadrature points i ∈ I are evaluated in relation to only one of the contact domains, in this case Γ (a)  C .Nevertheless, the projection on to Γ (b)  C will preserve conservation of momentum at the interface.When the domain to be projected on Γ (b)  C is C 0 , i w and | i J| evaluated over Γ (a)  C and projected on to a knot span ∈ R 2 parameterizing Γ (b)  C are confned within the C .But for higher-order continuity, i w and | i J| projected on a knot span ∈ R 2 are higher-order continuous across the knots and hence are distributed on the knot spans ∈ R 2 of the patch rather than be confned within a knot span ∈ R 2 when the continuity is C 0 on the patch.Hence, the accuracy can be inferred to increase with increase in continuity.Te accuracy associated with increase in the number of quadrature points can be explained as follows.As the number of quadrature points increase, the quadrature weight associated with a quadrature point decreases.Hence, for dissimilar meshes on Γ (a) C and Γ (b) C , lower the quadrature weights means lower the error in distribution of quadrature quantities across the knots in a given direction.Te efect of P on the accuracy of σ n for the collocation method and the variational formulation is shown in Figure 10, considering D 2 C 1 parameterization.In addition with Greville points, the collocation scheme is also considered with equally spaced grid points in the parametric space over the patch.Te total number of grid points is defned to be 81, and the total number of Greville points is 25 resulting from the NURBS parameterization.Meanwhile, the variational formulation was considered with 12 × 12 number of Gauss points in a knot span ∈ R 2 .It is evident that the variational formulation yields less oscillation of σ n for 16 Shock and Vibration both the extremes of P, where for P � 10 6 and P � 10 14 , the diferences in oscillation are ≈ 34 N/m 2 and ≈ 0.01 N/m 2 , respectively.While for the collocation method with the grid points, the diference in oscillation for P � 10 6 is ≈ 1700 N/m 2 which reduced drastically for P � 10 14 to be ≈ 0.06 N/m 2 .Collocation with the Greville points has the worst oscillation for the extremes of P, where for P � 10 6 and P � 10 14 , the diferences in oscillation are ≈ 3200 N/m 2 and ≈ 1800 N/m 2 , respectively.Tis is likely since the accuracy of the collocation scheme with Greville points is known to be highly biased depending on the choice of surface considered for collocation [21].To illustrate this,  18 Shock and Vibration we consider a new parameterization as shown in Figure 11 for contact patch test 2. Te results are shown in Figure 12 for the collocation method with Greville points and the variational formulation with the same setting as the preceding analysis.Considerable improvement for P � 10 14 can be observed with collocation of the Greville points, where the diference in oscillation is ≈ 0.16 N/m 2 in comparison to ≈ 0.05 N/m 2 for the variational formulation.Nevertheless, still for P � 10 6 , the diference in oscillation is ≈ 2300 N/m 2 in comparison to ≈ 0.055 N/m 2 with the variational formulation.For the case of a static analysis with contact, for which the contact constraint is satisfed with the penalty approach, given the proper choice of surface considered for collocation, collocation scheme with Greville points could be a good choice and more computationally efcient.But for accuracy at low values of P from the view point of the normal compliance approach, the variational formulation is more robust and the accuracy of σ n is less sensitive to P and the choice of surface.Te bias in the choice of surface for collocation could be overcome by the two-pass formulation though its numerical properties mainly in the context of CEA are not considered in this study.Te oscillation of σ n can indeed have an efect on steady-sliding equilibrium and the estimation of p(X) and Γ C for linearization around the equilibrium, for which the variational formulation could be highly efcient in comparison with the collocation method.
With various formulations in approximating the integrals defned over Γ C , the only diference is the relative variation of p(X) on Γ C for a given P. Te relative variation of p(X) infers the variation of σ n owing to the inaccuracies in the approximation of the integrals defned over Γ C , as observed with the contact patch test.Hence, the sensitivity of CEA for oscillations in σ n can be inferred through comparing the collocation method and the variational formulation, which show diferent characteristics of oscillation for σ n depending on P. Figure 13 shows comparison of CEA results for various values of P with the collocation method achieved through Greville points and the variational formulation.For collocation with Greville points, proper choice of contact surface is considered as inferred from the preceding discussion.It can be observed that for P � 10 9 , the natural frequencies I(λ) and their corresponding instabilities R(λ) do not coincide between the formulations.Te collocation method relatively seems to overestimate the instabilities R(λ) for I(λ) < 10KHz.For P � 10 12 , I(λ) and R(λ) are identical.Surprisingly, for P � 10 14 , the results are largely identical but still not as identical as P � 10 12 .We hypothesize this to be as a result of the increasing condition number of the matrices which shows the extremum of P considered.Tis indeed resembles the static case of the contact patch test in Figure 12 where for a large value of P, the accuracy of σ n is nearly identical for the collocation

Shock and Vibration 19
method and the variational formulation.One can hypothesize from the patch test that the approximation error which is profound at lower values of P is implied through the oscillation of σ n and consequently also has an efect on CEA.
In contrast, at high values of P, CEA is less sensitive to the approximation error similar to the patch test, even with the collocation method which only models concentrated points of force.Tis means that for low values of p, the dynamics of the system are much sensitive to the relative variation of p(X) on Γ C , which implies that characterizing the contact interface with global contact stifness P may not be realistic with smaller values of local contact stifness p(X).Despite the lack of the numerical benchmark to compare the results of CEA, the variational formulation can be considered to be the most accurate notably at low values of P, inferring from the results of the contact patch test.For the static analysis of the contact patch test, considerable improvement in the accuracy of σ n was achieved with increasing the continuity and the number of Gauss quadrature points.Tis property could be useful if CEA is highly sensitive to the error in the discretization of p(X).Tis is verifed in Figure 14 where no considerable variation could be observed except for small shift in frequencies at high frequencies, I(λ) > 13KHz, which could be as a result of convergence.Hence, it can be inferred that at low values of P, CEA is not sensitive to the scale of improvement in the accuracy of approximation provided by increasing the continuity or the number of Gauss quadrature points for the variational formulation, but it is certainly sensitive to the scale of improvement in the accuracy of approximation between the collocation and the variational formulation.To further show the sensitivity of collocation, we compare the collocation method between nearly identical number of grid points and Greville points.Te number of equally spaced grid points in the parametric space was defned to be 100 for the 104 Greville points inferred from the NURBS parameterization.In this case, the only diference is the positions in which the collocation is defned.Te comparison is shown in Figure 15, where the results are observed to be diferent.

Conclusion
Numerically, the results imply that for low values of contact stifness, the dynamical response is sensitive to the relative variation of contact stifness at the interface.Tis factor proves to be sensitive with approximations in defning a contact formulation.We show that the defned collocation method is inaccurate at low contact stifness and hence can  C in parametric space.20 Shock and Vibration be largely deemed to be inefcient with the normal compliance approach.Even though the collocation formulation was defned from the variational formulation, the analogy of modelling concentrated points of force can be in general extended to other collocation-type formulations and hence also the inaccuracy associated with it.Tis is also true in classical fnite element analysis with methods such as nodeto-node and node-to-surface formulations.While for large values of contact stifness, the dynamics is less sensitive to the relative variation of contact stifness and converges with increase in contact stifness.Hence, to satisfy Signorini conditions with the penalty approach where the penalty value can be interpreted as high contact stifness, the collocation method can be accurate with the right choice of collocation points and a collocation surface for the half-pass formulation.For computational efciency while preserving good accuracy, collocation points based on NURBS parameterization such as points defned with Greville abscissae which in general implies few number of collocation points can be even efcient.
For the normal compliance approach, notably at low contact stifness, preserving the integral nature of the variational formulation, even approximately without domain decomposition proves to be efcient and robust.We show that CEA could be analyzed robustly and accurately with the accuracy of such formulations.Te inference from the results realized with the isogeometric approach can also be extended to classical fnite element methods, except for the higher-continuity of the spline basis functions can be more accurate than with the classical fnite element c 0 basis functions.Such comparisons have been largely discussed in the literature.We also showed that the contact patch test cannot be relied to imply accuracy for CEA, where we referred to a lumped approach which is superior in the patch test, proving to be numerically unstable for CEA.

Figure 2 : 5 Figure 3 :
Figure 2: Description of the contact patch test.