Abstract
Flutter-type dynamic instability induced by friction is a highly nonlinear phenomenon and computationally expensive to model through transient analysis. An efficient way to make inference of such instabilities in a dynamical system is through analyzing the first-order effect of a perturbation at one of its equilibrium with eigenvalue analysis. The contact characteristics of such dynamical systems are typically modelled through the normal compliance approach with inference from experiments. In this case, the dynamical response of the system is implied to be sensitive to the contact stiffness modelled through the normal compliance approach. Typically, with the normal compliance approach, the continuum of the contact interface is approximated through a set of nonlinear springs which can be interpreted as a collocation method. Such approximations or the numerical implication of contact formulations in general for such problems is not largely studied. We focus on a variational formulation-based contact formulation without domain decomposition which is computationally efficient with small sacrifice in accuracy, where we imply that the dynamical response can be robustly modelled with the given accuracy. Further, we expose the inadequacy of the collocation method for such problems, where the dynamical system is observed to be sensitive to the extent of inaccuracy as a result of collocation for low values of contact stiffness. The inferences numerically imply the characteristics of the dynamical system for variation in contact stiffness.
1. 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]. The resulting noise generated by friction is more complex to model mathematically and cannot be typically explained by a single mechanism [2–4]. The complexity can be largely owed to modelling the contact interface properties and the mechanism of friction at various scales. The 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 identified as the latter case, where as a result of perturbations around equilibrium, the system reaches a new steady-state characterizing self-excitation behavior. This can be understood as coalescence of modes where two modes exist at a same frequency, leading to self-excitation 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 flutter instability. This is evident in the literature even from early analyses based on lumped models, which greatly improved the understanding of flutter-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. The 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. The presence of the limit cycle is inferred through complex conjugate pairs of eigenvalues or its equivalent form 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 finite 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. The 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 defined 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. The variational formulation is also satisfied through collocation for comparison to show the efficacy 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 definition of velocity terms. We expand the definitions based on the normal compliance approach [7, 16] where the unknown contact boundary conditions of normal stress are defined to be displacement dependent. This 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 field 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 higher-order continuity of spline basis functions could be exploited [20–22].
The meaning of penalty term in classical sense is essentially to satisfy the contact constraints, where an ideal penalization is of order infinity. But in a numerical context, penalization is taken to be a high value to enforce the constraints, considering the condition number of the matrices. This could be unrealistic for a dynamical system, where the dynamical response is sensitive to stiffness at the interface which is essentially the normal compliance approach models. This is in contrast to the zero compliance model of Signorini conditions which the penalty approach essentially tries to satisfy. This means that the discretization scheme should also be accurate at lower values of penalization which could reflect the contact stiffness to be modelled with the normal compliance approach. This 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 stiffness which is classically ignored.
We also make some inference on numerical characteristics of CEA in relation to contact stiffness and the accuracy of modelling normal stress at the contact interface, such that the sensitivity of CEA could also be implied. This is important in approximating the integrals defined over a contact domain such that it characterizes the contact interface in the interest of CEA independent or with little inference from experiments. This could be useful in applications such as optimization of systems for friction-induced dynamic instability through numerical simulations. This 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.
2. Modelling
We give a general overview of continuum formulations in this section which is useful to define 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 define the steady-sliding equilibrium problem § and the problem of the dynamics of perturbation around the equilibrium §. Followed by, we define the problem for the dynamics of perturbation explicitly §, where we define explicitly the quantities of interest required to be calculated from steady-sliding equilibrium. This is to understand the approximation of the integrals defined over a contact domain for given quantities of interest which are to be defined explicitly in the interest of studying contact formulations given in §. Further, for the dynamics of the perturbation problem defined 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 defined §.
2.1. Initial-Boundary Value Problem
The initial-boundary value problem for a domain in contact with a rigid body in the presence of friction can be given aswhere , , with defining the boundary of , and defining the unit normal vector on . Under isotropic material consideration, the constitutive equations can be defined as , where and are 3D Lamé parameters expressed in terms of Young’s modulus and Poisson’s ratio . The kinematic relation for the strain tensor under infinitesimal displacement is defined as . For a system with domains and in contact, , where represents the material coordinates. on is defined as the outward normal projection from .
Signorini conditions (1b) and Coulomb’s law (1c) define 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]. This leads to the view of the normal compliance approach [16] to define 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 defined by as the function of relative normal displacement (Considering is a rigid body, in this case) at the interface in relation to the initial gap , given aswhere allows only a positive value. This can be extended to friction aswhere the parameters , , , and 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, and in this case. With known as the active set which satisfies the normal compliance constraint at a given time, the problem can be expressed as variational inequality [26], given aswhere is the directional derivative. The 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 nondifferentiable friction functional, mainly as , 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 and which satisfy the inequality at the initial time, the system can be numerically integrated with time.
We focus on modelling flutter-type dynamic instability through classical theories of linear analysis, where the first-order effect of perturbations around a fixed point is analyzed. Hence, the stability of the dynamical system with frictional contact can be characterized by determining the fixed points which are typically of quasistatic or steady-sliding equilibrium depending on the external forces and defining the dynamics for the perturbation around such fixed points. We only focus on steady-sliding equilibrium for the following discussions.
2.2. Steady-Sliding Equilibrium
The dynamical system can be expressed as steady-sliding equilibrium when no net acceleration is present in frictional contact. This 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. The steady-sliding equilibrium explicitly characterizes slip condition at constant sliding velocity, and hence, at the equilibrium, on . Further, empirically, the knowledge of sliding direction at any point on is defined to be known a priori, and hence, . This means that the nondifferential friction term of (4) as a result of the inequality (3) can be expressed as variational equation at pure slip state (The normal compliance terms in modelling 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. The problem of steady-sliding equilibrium can hence be expressed as
The nonlinear problem can be solved through the Newton–Rhapson method to obtain the steady-sliding solution . It should be noted that for steady-sliding case with friction, nonunique solutions can exist [10, 29].
2.3. Initial-Boundary Value Problem for Perturbation
The perturbation of the displacement field relative to the steady-sliding equilibrium can be given aswhere and correspond to perturbed displacement field and time, respectively. The dynamics of the perturbation can be deduced by substituting in the variational equation form of (4) (The variational equation can be deduced by regularization of the nondifferentiable 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 effect) and subtracting with (5) to obtain
The 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 sufficiently characterized linearly. This brings the question of linearizing normal compliance.
We start with the discussion of linearizing the operator under necessary conditions. The perturbation on may not result in the same active set , which introduces nonlinearities defined by (Separation on defines which can only be taken into account nonlinearly. This is also true for the case contrary to separation). This can be linearized with the hypothesis that the onset of instability to model occurs close to the equilibrium such that the perturbation is small to define as stationary. But this may not be true since is stationary only when (in the view of normal compliance, the perturbation for no variation of defines the variation of . This would have been impossible with Signorini conditions which can be regarded as the zero compliance model, where leads to strict contact separation on ), implying that must be true on at steady-sliding equilibrium, conceivably as a result of external forces. This can lead to a new active set in the current configuration that satisfies on , which depends on the nature of the perturbation and at steady-sliding equilibrium. Perhaps, this could reveal if the linearization of the operator is realistic for a given system. To simplify, we consider as stationary for the perturbation , which evades the operator . With further linearization for the polynomial terms of the normal compliance, the first-order effect of the perturbation close to can be expressed as
An illustration of the linearization is shown in Figure 1. It is also possible that the nonlinearity defined by the normal compliance parameters could mean that the linearization at very low values of could yield very low . Hence, the assumption of stationary could be justified in this case.
We provide discussions on the presumed characteristics of friction for the perturbation to express the dynamics of the perturbation to be linear. Even though steady equilibrium is expressed by static forces, for the perturbation , 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 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 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, can also be taken to be stationary even in the presence of relative motion. But for the perturbation of such problems, the variation of forces on at slip state means that the assumption of stationary may not be realistic. This is true unless for the dynamical system (4) is negligibly relative to the natural frequencies of the system, which is taken to be the case here. This 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 define which can be interpreted as contact stiffness for perturbations close to steady-sliding equilibrium. For the following discussions, we do not relate to the experimental determination of the parameters and . We mainly focus on the finite element discretization of the integrals defined over for the dynamics of the perturbation and the consequence of on the discretization schemes. Hence, the focus on the influence of discretization schemes in the estimation of steady-sliding equilibrium is not considered in this study. effectively defines at steady-sliding equilibrium, also at which normal compliance is linearized to define . But and can be defined explicitly, for which a new initial-boundary value problem (9) can be considered. This is to simplify the comparison of discretization schemes for (17) with constant and , rather than to deduce from steady-sliding equilibrium where and could indeed vary depending on the discretization scheme.
2.4. Initial-Boundary Value Problem for Perturbation: Explicitly Defined
The 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 is introduced to distinguish the domains in contact. With and known a priori, either from steady-sliding equilibrium or defined explicitly, the strong form of the dynamics for can be expressed as
The parameters and are defined to be constant on between all the domains in contact. The variational formulation of the problem can hence be expressed as
The traction force can be decomposed as
Hence, with the expansion of the traction forces for domains in contact, equation (10) can be expressed as
To simplify, we consider contact between two domains and , with the derivation of traction forces on . Hence, the inner products and can be expressed as
The traction forces on can similarly be defined from the conservation of momentum as and .
2.5. Finite Element Approximation
We introduce the function space in which we seek the solution , where the Sobolev space and the norm: . A subspace can be defined as the restriction of on , i.e., , defined through the normal compliance approach. Even though is typically defined to be in , the unknown a priori conditions of on imply , where is the dual of the space .
The solution can be expressed as . With the finite element approach, we define a finite-dimensional function space , and hence, there is some bound on . The approximation of in can be expressed as . We focus on the definition of the space with the isogeometric approach in the next section. But in general, we give the matrix form of the variational formulations derived with .
The finite element approximation of for equation (5) considering the two domains and in contact can be expressed aswhere and . defines the domains which are not , i.e., in the case of contact between and , implies and vice-versa. implies the realization of with the parameterization of domains. The abovementioned expression can be expressed in the matrix form aswhere is the stiffness matrix, and are the nonlinear traction force vectors of contact and friction, respectively, and with being the classical force vector implying external forces. The linearization of the integrals defined over in solving the steady-sliding equilibrium (5) through the Newton–Rhapson method takes the similar form of. This means that the tangent matrices of and to solve for has the form and , respectively, to be discussed.
Similarly, with the substitution of equation in equation (7), for , results in the following form:which can be expressed in the matrix form as
The properties of the matrices can be given as follows, where and are symmetric and positive definite matrices. is also symmetric which essentially defines the conservation of momentum at the interface for normal contact. While is nonsymmetric which defines 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 for , 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 . Given an eigenvalue of the form , the part models the oscillatory behaviour while the stability of the given mode is defined by . can be understood to characterize an unstable mode and defines the divergent nature of the solution. Similarly, can be understood to characterize a stable mode which can be interpreted as damping. The occurrence of a pair of eigenvalues defines the presence of a limit cycle which is of interest for applications concerning flutter-type instability induced by friction. Damping can also be factored through the definition of a matrix , typically through modal or Rayleigh damping. In the scope of modelling rotational inertia effects, gyroscopic matrix can also be defined, where and are both velocity dependent. We factor out damping and gyroscopic effects in our model to mainly focus on contact and friction modelling.
2.6. Isogeometric Approach for Discretization of the Contact Integrals
We define the space through the isogeometric approach to take advantage of the properties of the spline basis functions. The higher-order continuity of the spline basis functions between the knots provides superior approximation properties. The influence 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]. This 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. The higher-order solution continuity is a consequence of the geometric continuity provided by spline parameterization. The geometric continuity greater than is also useful to define the actual normal of the contact surface and also avoids numerical discrepancies in large displacement contact problems, as a result of the faceted continuity of classical finite element meshes.
2.7. NURBS for Approximation
B-spline basis functions can be defined by Cox de Boor’s formula as follows:where is defined recursively for to obtain a curve of degree , which starts with a piecewise constant at . Naturally, a uniform knot vector can be defined as , where any is uniformly spaced. The knot vector need not be equidistant, and the multiplicity of a knot by in the knot vector decreases the continuity by across the knot , which defines a nonuniform knot vector. Through B-spline basis functions and a knot vector , a B-spline curve can be defined with the basis functions and its coefficients as follows:
The coefficients are the control points, with being the dimension of the space. The definition of a weighing parameter associated with a respective basis function , normalized defines rational B-splines where it respects the partition of unity, given as follows:
The parameter provides a new dimension for controlling the geometry through projective transformation, while the affine transformation is achieved by . Hence, the combination of nonuniform knot vectors and rational basis functions defines NURBS. The higher dimensional NURBS are a natural extension of its 1-dimensional precursor through tensor product definition where the order of the tensor is the same as the dimension of the geometry. Hence, a tensor product NURBS surface can be defined aswhich is supported by knot vectors and , for the domain , with net of control points . Similarly, a tensor product NURBS volume can be defined aswhere the knot vectors are given as and . The abovementioned expression can be simply expressed in the matrix form as . We also give the following notations, where to express the same degree for all the knot vectors defining a NURBS parameterization in is simply defined as , i.e., . Similarly, for the same continuity between the internal knots in all the knot directions is expressed as .
The parameterization of a domain as an initial geometric description through NURBS can be expressed as , , where defines 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: . The analysis-suitable parameterization (for simplicity of the notation, we define to be the default notation for analysis-suitable parameterization of a domain ) can be achieved through the refinement of with one or several of the refinement methods: , , and , where can be defined as to take into account of the modified knot vectors and control points. The isogeometric approach for the approximation of a solution is achieved through the bases , where for the vector-valued function space , the vectorial definition of the bases can be defined as
This can be expressed in the matrix form as , with being the identity matrix. Hence,
In an abstract sense, the bases in parametric space are transformed into the bases in physical space using the push-forward operator , where the bases are defined with the property . Hence, the approximation of a field variable on is defined through all the bases spanning the finite-dimensional function space . For the isogeometric approach, we express the finite-dimensional space as and its associated bases as . The approximation of can be defined as , expressed in the matrix form as , where
The variational form of contact and friction given in equations (13) and (14) expressed in the function space can be defined as
The expansion of the abovementioned expressions also generalizes for the definition of the tangent matrices for the nonlinear terms of contact and friction in equation (15). The parameterization of the domains and with NURBS can be expressed as and . The active sets and is defined with the property , where in this case is taken to be the outward normal projection from to . This means that , hence for that parametrizes , a projection exists that maps on as . For the following explanations, we detail the derivation of traction forces on , which can be extended to depending on the considered discretization scheme.were . Essentially, the abovementioned equations are a construct of the matrix form as
For , the definition of the integrals is straight forward. In contrast, for and , we have terms of the form which can be expanded aswhere the integrals are simultaneously defined over the bases of the two domains in contact since and . Even though the definition of integral for on exists through the projection , for dissimilar meshes at the contact interface, the definition of numerical quadrature scheme for the integrals of such form requires domain decomposition to find the common span: . This means that the integral exists between and only on the span where the projection exists and hence requires a quadrature scheme specific 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 verifies the conservation of linear momentum at the contact interface. Since we deal with contact between flat surfaces, for the normal compliance approach which can be viewed as coupling of displacement field between the contact surfaces, conservation of angular momentum is implicitly satisfied.
We give some intuition through the collocation approach to satisfy the conservation of momentum at the contact interface and extend it to the definition of quadrature scheme without domain decomposition. For collocation, the integral on one of the domains, in this case, – can be defined as , where is the set of points on 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) satisfied discretely at the collocation points takes the form as follows:where , , and . This satisfies conservation of momentum at the contact interface, even though the integral cannot be defined accurately. For any , the following relation holds:
This means that any quantity defined at a collocation point over is projected equally over the bases in and . This is as a result of the properties of the NURBS basis functions similar to classical finite element basis functions which satisfy partition of unity. Since the integral is satisfied only at discrete points, the solution may not be to the necessary accuracy. The collocation strategy can be replaced by a numerical quadrature scheme as where , in this case, corresponds to the quadrature points with being the quadrature weights. But the notion of on may not be realistic when is defined for . This 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 aswhere and .
The integral defined over the parametric space for the isogeometric approach takes the form.where and being the collocation point in the parametric space and is the Jacobian matrix for the mapping . is the corresponding map of in the parametric space, which can be determined through the Newton–Rhapson method in solving for . Hence, there exists a mapping for which .
From the conservation of momentum at the interface, the following relation holds and . Hence, the traction stresses on can similarly be defined as
It should be noted that for the abovementioned equations, even though the traction forces are defined over as , the quadrature points’ set is determined only on , where its corresponding projection on is defined through the projection . This also includes Jacobian which is evaluated only on , similar to equations (36) and (37). This greatly simplifies the finite 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 effect on accuracy which is discussed in §. With the abovementioned definitions, the matrices and for the system take the form.
Since the quadrature rule or collocation is defined with respect to only one surface, it is commonly termed as half-pass. The role of and can be switched, and the average over the two half-pass can be taken into account to define 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 and , in contrast to the form and , of the system matrices and , respectively, correspond to the case . With the relation inferred from (34), and can be lumped. Hence, with equation (30) for , lumping can be achieved by defining , being the Kronecker delta product, where the inner product is expressed with the property. defines the matrices of the form and , where the relation inferred from (34) holds. This preserves the conservation of momentum at the interface.
Collocation can be achieved with the weak form considering , , and 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 effect of the area for a collocation point can be taken into account through and 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 fixed which is in contrast to the Gauss quadrature scheme which can be defined 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 define the projection , and the number of Greville points is the same as the number of control points.
3. Results and Discussion
The variational formulation satisfied through collocation with Greville points is simply termed as the collocation method, unless otherwise specified. 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. The 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. and are considered to be the same for all the domains in a given analysis. For the variational formulation, we define 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 specified. The 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 through the considered contact patch test, it should be noted that for steady-sliding equilibrium is given by the relation . Hence, the accuracy of will have a subsequent effect on , though the sensitivity of and on steady-sliding equilibrium is not considered in this study. Further, all the results of CEA are obtained with . To proceed with the following discussions, we introduce the global contact stiffness at the interface as . Hence, a given is constant at the interface, and the local contact stiffness varies depending on the accuracy of the formulations.
For the contact patch test, we consider to be constant on , and hence, an ideal discretization should yield . But for a given , the local contact stiffness can vary depending on the approximation in obtaining a formulation. The discretization of has direct influence on contact normal stress 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 . For this reason, we focus on considerably low values of in the interest of the normal compliance approach. The setting for the contact patch test is shown in Figure 2, where we consider two cuboid domains in contact with two different meshes visible through the knot vectors of the NURBS parameterization. The bottom face of is fixed, while the top face of is applied a uniform pressure of and constrained to move normal to . The material parameters and are defined to be the same for both the domains. Analytically, one can infer that on , which provides good benchmark to compare calculated with different formulations. The contact patch test can be defined by considering equation (15) for and . With further simplification of and , the matrix form of the problem can be solved for . Given and known a priori, for , can be expressed as .
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 for a given , where for the higher value of , the approach even passes the contact patch test to machine precision. For CEA 19, we consider a simple disc-pad system which can give rise to dynamic instability induced by friction. Similar to the contact patch test, is defined explicitly to be constant on , as discussed in §. This is because determined through steady-sliding equilibrium can vary depending on the discretization scheme to define 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 define the following notations and to express the displacement field of the disc and the pad, respectively. Even though the scale of the displacement field has no sense for the computed eigenvectors, comparison can be made relatively with the displacement field for a given eigenvector. It should be noted that we deal with complex eigenvectors where the imaginary part of an eigenvector characterizes phase difference of the displacement field, which will not be considered here. The 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 field elsewhere. The arrow on which is plotted also traverses the boundary of the contact interface (Boundary of the contact interface can be defined as for the contact interface ) where strong solution gradient for can be observed. This is in contrast to the approach without lumping where is smooth across . In this case, the dynamic response of the approach with lumping and without lumping is observed to be different as shown with the results in Figure 6 where the natural frequencies and its corresponding instabilities largely do not coincide. Empirically, the behaviour observed at the contact interface for the lumped approach seems unrealistic since and are observed to be locked with each other relative to the rest of the domains in contact. This could also be a consequence of the properties of NURBS basis functions and hence can be verified with classical finite element basis functions, which is not considered here. It should be noted that the role of mesh refinement at plays an important role in capturing the solution gradient since the higher-order continuity of the spline basis functions will have strong propensity to define a smooth solution across . In our case, multipatch parameterization of disc is considered where one of the patches consist of . Hence, across , refinement was achieved in the patch that contains such that the knots can trace along side and on . But for more complicated shape of , strategies based on RBQ [34], or in the scope of local refinement 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.
(a)
(b)
(a)
(b)
As it can be inferred from the contact patch test in Figure 3, the accuracy of can be improved with increase in . For the contact problems defined by static analysis, when the contact constraint is satisfied with the penalty approach, it is necessary that is large where in an ideal case . But for dynamical systems, contact stiffness models the interface properties which have influence on the resulting dynamical response. It has been observed that the value of contact stiffness 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. The necessary accuracy indeed depends on the sensitivity of the approximation in estimating the steady-sliding equilibrium, the effect it has on linearization for CEA, and the evaluation of and with CEA. In the interest of computational cost, as the order of continuity is increased ( refinement) with increase in degree ( refinement), the increase in the number of control points is of very low order. Even though the increase in the degrees of freedom of the finite element matrices is small for simultaneous and refinement, this can reduce the sparsity of the matrices and hence can increase the computational cost in terms of both memory and time. Hence, can be considered only up to a reasonable value of depending on the required accuracy.
With the following results, we analyze the significance of the higher-order continuity of the spline basis functions and number of Gauss quadrature points on the accuracy of . The influence of the higher-order continuity is shown in Figure 7 for various and considering . It can be observed that the oscillation in decreases with increase in continuity. This can be clearly seen with comparing the cases of and defined with equal number of Gauss points between the knot spans ∈ (span ∈ , considering that is defined on the surface parameterized by knot vectors and ), where shows reduced oscillation of . While for , the solution is observed to be only continuous at the knots. The results of and are shown in comparison to for which the oscillation of is negligible. The accuracy of can be improved with increasing the number of Gauss quadrature points, which is shown in Figure 8. For the case of in Figure 9, the number of Gauss points is also increased with the degree by default as specified before, and hence, the addition of Gauss points also contributed to increase in accuracy along with the contribution from higher-order continuity. Nevertheless, even does not satisfy the contact patch test at , shown in Figure 9. It can be inferred that even though the accuracy of increases profoundly with initial refinement 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 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 is greater than for . This can be hypothesized as resulting from the quadrature scheme defined only over , i.e., the quadrature points, Jacobian, and the quadrature weights are evaluated only on , while only their corresponding projection is defined over . 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 defined over is simplified with considering the quadrature scheme on only one of the domains in contact, while its projection is defined on the other. This is not very intuitive since the quadrature quantities such as quadrature weights and the determinant of Jacobian of the quadrature points are evaluated in relation to only one of the contact domains, in this case . Nevertheless, the projection on to will preserve conservation of momentum at the interface. When the domain to be projected on is , and evaluated over and projected on to a knot span ∈ parameterizing are confined within the knot span ∈ as a result of , where for dissimilar meshes between and , the projection can lead to unrealistic definition of quadrature quantities on . But for higher-order continuity, and projected on a knot span ∈ are higher-order continuous across the knots and hence are distributed on the knot spans ∈ of the patch rather than to be confined within a knot span ∈ when the continuity is on the patch. Hence, the accuracy can be inferred to increase with increase in continuity. The 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 and , lower the quadrature weights means lower the error in distribution of quadrature quantities across the knots in a given direction.
The effect of on the accuracy of for the collocation method and the variational formulation is shown in Figure 10, considering parameterization. In addition with Greville points, the collocation scheme is also considered with equally spaced grid points in the parametric space over the patch. The total number of grid points is defined to be 81, and the total number of Greville points is 25 resulting from the NURBS parameterization. Meanwhile, the variational formulation was considered with number of Gauss points in a knot span ∈ . It is evident that the variational formulation yields less oscillation of for both the extremes of , where for and , the differences in oscillation are and , respectively. While for the collocation method with the grid points, the difference in oscillation for is which reduced drastically for to be . Collocation with the Greville points has the worst oscillation for both the extremes of , where for and , the differences in oscillation are and , respectively. This 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, we consider a new parameterization as shown in Figure 11 for contact patch test 2. The 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 can be observed with collocation of the Greville points, where the difference in oscillation is in comparison to for the variational formulation. Nevertheless, still for , the difference in oscillation is in comparison to with the variational formulation. For the case of a static analysis with contact, for which the contact constraint is satisfied 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 efficient. But for accuracy at low values of from the view point of the normal compliance approach, the variational formulation is more robust and the accuracy of is less sensitive to and the choice of surface. The 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. The oscillation of can indeed have an effect on steady-sliding equilibrium and the estimation of and for linearization around the equilibrium, for which the variational formulation could be highly efficient in comparison with the collocation method.
With various formulations in approximating the integrals defined over , the only difference is the relative variation of on for a given . The relative variation of infers the variation of owing to the inaccuracies in the approximation of the integrals defined over , as observed with the contact patch test. Hence, the sensitivity of CEA for oscillations in can be inferred through comparing the collocation method and the variational formulation, which show different characteristics of oscillation for depending on . Figure 13 shows comparison of CEA results for various values of 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 , the natural frequencies and their corresponding instabilities do not coincide between the formulations. The collocation method relatively seems to overestimate the instabilities for . For , and are identical. Surprisingly, for , the results are largely identical but still not as identical as . We hypothesize this to be as a result of the increasing condition number of the matrices which shows the extremum of considered. This indeed resembles the static case of the contact patch test in Figure 12 where for a large value of , the accuracy of is nearly identical for the collocation method and the variational formulation. One can hypothesize from the patch test that the approximation error which is profound at lower values of is implied through the oscillation of and consequently also has an effect on CEA. In contrast, at high values of , 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. This means that for low values of , the dynamics of the system are much sensitive to the relative variation of on , which implies that characterizing the contact interface with global contact stiffness may not be realistic with smaller values of local contact stiffness .
(a)
(b)
(c)
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 , inferring from the results of the contact patch test. For the static analysis of the contact patch test, considerable improvement in the accuracy of was achieved with increasing the continuity and the number of Gauss quadrature points. This property could be useful if CEA is highly sensitive to the error in the discretization of . This is verified in Figure 14 where no considerable variation could be observed except for small shift in frequencies at high frequencies, , which could be as a result of convergence. Hence, it can be inferred that at low values of , 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. The number of equally spaced grid points in the parametric space was defined to be 100 for the 104 Greville points inferred from the NURBS parameterization. In this case, the only difference is the positions in which the collocation is defined. The comparison is shown in Figure 15, where the results are observed to be different.
4. Conclusion
Numerically, the results imply that for low values of contact stiffness, the dynamical response is sensitive to the relative variation of contact stiffness at the interface. This factor proves to be sensitive with approximations in defining a contact formulation. We show that the defined collocation method is inaccurate at low contact stiffness and hence can be largely deemed to be inefficient with the normal compliance approach. Even though the collocation formulation was defined 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. This is also true in classical finite element analysis with methods such as node-to-node and node-to-surface formulations. While for large values of contact stiffness, the dynamics is less sensitive to the relative variation of contact stiffness and converges with increase in contact stiffness. Hence, to satisfy Signorini conditions with the penalty approach where the penalty value can be interpreted as high contact stiffness, the collocation method can be accurate with the right choice of collocation points and a collocation surface for the half-pass formulation. For computational efficiency while preserving good accuracy, collocation points based on NURBS parameterization such as points defined with Greville abscissae which in general implies few number of collocation points can be even efficient.
For the normal compliance approach, notably at low contact stiffness, preserving the integral nature of the variational formulation, even approximately without domain decomposition proves to be efficient and robust. We show that CEA could be analyzed robustly and accurately with the accuracy of such formulations. The inference from the results realized with the isogeometric approach can also be extended to classical finite element methods, except for the higher-continuity of the spline basis functions can be more accurate than with the classical finite element 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.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.