An Isogeometric Analysis Approach to Kinematics of Spatial Rigid Multibody Systems with Imperfect Joints

This paper proposes a novel generic methodology for kinematics of spatial rigid-multibody systems with and without lubricated joints. In this method CAD surface representations in the form of non-uniform rational B-splines (NURBS) are used to address the interface kinematics. This eliminates the time and effort needed to manually parameterize the interface geometry, by enabling a direct use of the engineering designs encapsulated in CAD systems. Furthermore, the use of NURBS for surface representation allows integration of tribodynamics into an isogeometric analysis (IGA) setting. The kinematic formulation is based on a new implicit matrix approach for implicitization of CAD surfaces in three-dimensional space. The construction of such implicit matrices and their properties are explained, and explicit expressions for the gap height distance, velocity and relative velocities in a general clearance joint are provided.


Introduction
The presence of clearances in mechanical joints is practically inevitable due to manufacturing tolerances, wear and material deformation (Flores et al., 2008). Friction and impact forces present in a joint clearance, mainly in absence of lubrication, contribute to increasing vibration amplitude and reduce system reliability, stability, life, and precision (Erkaya and Uzmay, 2010). Suitable mathematical models for dynamic analysis of multibody systems with clearance joints are therefore of utmost importance, perhaps even more so for emerging technologies. One example is the novel power transmission technologies such as digital fluid power where challenges arise from the non-smooth dynamical behavior, which increases pressure oscillations Pedersen, 2015, 2016a,b) and thereby vibration amplitudes. Moreover, pressure pulsations stemming from variations in the digital displacement machines' switching patterns require careful assessment of micromechanical inertial effects in lubrication films Johansen et al. (2015).
In the past decade a great number of studies about multibody dynamics with clearance joints have been published. A literature review concerning dynamics of clearance joints in multibody systems has been presented in Tian et al. (2018). A thorough methodology for kinematics and dynamics of imperfect joints in multibody mechanical systems, which covers clearances in planar revolute joints, translational joints, spherical joints, and three-dimensional revolute joints has been presented in Flores et al. (2008). Applications of this modeling approach in lubricated multibody sys-tems can be found in (Tian et al., 2013;Flores and Lankarani, 2010;xin Xu and gang Li, 2012;Tian et al., 2009Tian et al., , 2011Machado et al., 2012;Zhao et al., 2016).

Geometrical and kinematical aspects of tribological joints
In order to characterize a tribological joint, it is necessary to develop a mathematical model that describes the geometrical clearance. In such type of joints the bodies are free to move relative to each other in an unconstrained manner, in contrast with kinematically constrained joints in traditional multibody models. As an example, consider a journal bearing consisting of a shaft rotating at an angular velocity, ω, within a stationary bushing. Its basic geometry is shown in Fig Because of the eccentricity given by the distance O b O j between the axial centers of the bearing and journal the clearance gap is not uniform around the circumference of the bearing. The gap width is measured along a bearing radius. It can be identified in Figure 1 as the distance h between the corresponding points P b and P j of the bearing and journal surface, respectively. The interface geometry is defined by the gap height distribution h, which from examination of the triangle O b O j P b is found to be h = C(1 + cos θ).
Where C (C = r b − r j ) is the radial clearance, i.e., the difference between the radii of the bearing and the journal, and ( = O b O j /C) is the eccentricity ratio. The boundary wall velocities of the corresponding points P b and P j from rigid body motion is also of concern, as these can affect the pressure load capacity, friction force and fluid film thickness (Askari and Andersen, 2019). These velocities include both rotational and translational motion of the joint components. The wall boundary velocities are also important when studying impact dynamics. Here the normal relative velocity determines whether contacting bodies are approaching or separating, and the tangential relative velocity determines whether the contacting bodies are sliding or sticking (Flores et al., 2008). A further characteristic is the gap height velocity ∂h/∂t, which describes relative velocity at which two surfaces approach each other in the perpendicular direction. This is known as the squeeze film action in hydrodynamically lubricated bearings and can provide temporarily increased load capacity when a bearing is subjected to abnormally high load (Stachowiak and Batchelor, 2005).
Conclusively, in the analysis of all tribological joints, interface geometry and velocities must be well-defined, since these are directly connected to pressure distribution, impact and performance characteristics such as load capacity, friction force and oil leakage. A full description of the tribological kinematics must include the gap geometry, wall boundary velocities and gap height velocity.

Tribological kinematics with NURBS surfaces
The ambition to apply computational tools for tribological design optimization is severely constrained by joint specific analytical geometric and kinematic models. An approach to address this issue is to merge the design model and the analysis model. This idea is based on the recent concept of Isogeometric Analysis (IGA) introduced by Hughes et al. (Hughes et al., 2005). The IGA framework allows to compute the analysis solution on the exact geometry defined by NURBS (the de facto standard CAD representation), instead of a discretized geometry leading to gains in both efficiency and accuracy. In addition, the potential ability of IGA to integrate CAD geometry, FEA and design optimization is considered a key advantage, as a strong interaction between the models of design, analysis and optimization is essential to obtain an optimal design (Wall et al., 2008). Consequently, the use of NURBS in the tribodynamic model formulations provides new possibilities for the geometrical optimization of the tribological interfaces. The topic of this paper belongs to the field of dynamic analysis of mechanisms with clearance joints. A thorough search of the relevant literature yielded only few works concerning dynamic analysis of clearance joints in the IGA framework. In Liu et al. (2016) IGA was successfully applied to solve the Reynolds equation for lubrication of a piston-cylinder interface of an internal combustion engine. The piston inertial dynamics was coupled with the lubricated hydrodynamics of the piston-cylinder interface to simulate a complete working cycle of the piston-cylinder system. In this study IGA showed to be more efficient for solving the Reynolds equation compared to classical FEM. The authors pointed out that IGA might play an important role in promoting accuracy and efficiency in tribodynamic simulation and should in the future be extended with aspects such as thermal effects, elastic effects and mixed friction. An IGA framework was developed by Pi and Zhang (2019) to study contact interactions in dry revolute clearance joints of planar multibody systems. In this approach a NURBS multipatch model was used to describe the interface geometry and the contact-detection search is performed using the closest point procedure (CPP), which is commonly used in contact mechanics problems. The authors studied sticking and sliding friction and concluded that the proposed modeling and simulation within the framework of IGA has a clear advantage in terms of its ability to include e.g. exact geometry of the contact interface, different constitutive laws, deformations and strains into the contacts mechanics formulation. Ćerimagić et al. (2018) proposed an isogeometric approach to the coupling of fluid film mechanics with the interface dynamics of a hydraulic radial piston motor. The modeling methodology was exemplified for the pistoncylinder interface in the motor. This was the first attempt to develop a generic tribodynamic framework useful for design, analysis and optimization of different motor topologies. This paper consolidates the work in (Ćerimagić et al., 2018) by giving a rigorous mathematical definition of the tribological kinematics for a generic joint.
This paper shows that when the tribological surfaces are given in a NURBS form, possibly extracted from a CAD model, the clearance geometry can be constructed automatically by using an implicitization technique. The proposed kinematical methodology in this paper allows developing tribodynamic simulation tools using NURBS representations for all design and analysis tasks.
Thus, the objective of this paper is to present an isogeometric analysis approach to kinematics of a multibody system with spatial clearance joints. In the proposed methodology, the interface geometry is considered in NURBS forms and the kinematics is formulated exclusively on that basis. This is considered the first necessary step in developing a generic multibody tribodynamic simulation tool that can take into account various machine design topologies.

Outline
The rest of the paper is organized as follows. Section 2 briefly introduces the NURBS surface maps used to describe the interface geometry and introduces the mathematical model of a generic joint with clearance in a multibody system. Section 3 provides an in depth description of the algorithm to find the corresponding surface points and ultimately the film or gap geometry. The gap height velocity and relative velocities are derived in Section 4. Conclusions are drawn in Section 5. It needs to be highlighted that the present paper aims to present a general methodology for kinematic modeling of real joints in multibody mechanical systems. The focus is on providing rigorous mathematical definitions, therefore numerical examples are not in the scope of this paper.
2 Joint kinematics using NURBS

NURBS interface geometry
In this section the modeling of the interface geometry of a general tribological joint is described. The surfaces of the interacting bodies, see Figure 2, are considered to be represented by NURBS (Non-Uniform Rational B-Splines). NURBS are chosen here due to their common use in computer graphics packages and their ability to represent a broad range of geometries. A NURBS surface is defined by the geometrical map where the rational basis functions are given by This mapping is defined by two knot vectors Ξ 1 and Ξ 2 of degree d = (d 1 , d 2 ), respectively, as well as the weights w i,j and the control points For example, the surface map in Equation (1) is used to define the surface geometries in Figure 3. The B-spline basis functions N i,d (ξ) of degree d in the parameter ξ used in the construction of NURBS surfaces are defined recursively (Piegl and Tiller, 1995), starting with the step functions for degree zero (d = 0) For degrees higher than zero (d > 0) they are defined by a linear combination of (d−1)-degree basis functions

Reference body
Corresponding body Figure 2: Generic joint with clearance in a multibody system.

Kinematic chain
In the general setting shown in Figure 2, the tribological interface is divided into a reference body and a corresponding body, which both have general threedimensional motion. In order to fully describe the joint geometry, the gap height vector denoted by H must be established. The first step towards a mathematical description of the gap geometry is to define coordinate systems in which it will be described. The coordinate systems defined by their origin and orthogonal coordinates, Ox, are shown in Figure 2. Here the frame O G x G is the fixed (inertial) frame, hereafter referred as the global frame. The frames O r x r and O c x c are rigidly attached to the center of mass of the reference and corresponding body, respectively. The surface description in Equation (1) is coordinate-independent. This means, the properties of the surface itself are not affected by the choice of coordinate system. Thus, in this paper the surface points S r (ξ ξ ξ) and S c (ξ ξ ξ) are considered in some arbitrary surface frame O sr x sr and O sc x sc , respectively, which are rigidly attached to their bodies. To simplify notation in the position analysis, the configuration of the system in Figure 2 is managed in the configuration space SE(3) (special Euclidean group), which is defined as Each element of the group SE(3) is a rigid body transformation that maps one coordinate system into another. The rigid body transformation is described by a homogeneous transformation matrix, or simply transformation matrix for short, where R corresponds to the orientation of the rigid body and p the translation. The inverse of the transformation matrix is given by Therefore, the transformation matrix representing the coordinate frame O r x r relative to the global frame O G x G , has the form The configuration of the surface frame O sr x sr relative to the global frame O G x G is then defined by the composition rule as Figure 3: Examples of surface geometries represented by NURBS with their control points (filled dots) and control net (dashed line).
In this work, the gap height in the tribological joint is defined with respect to the boundary geometry of the reference body, as illustrated in Figure 2. More precisely, the gap height is the Euclidean distance between the corresponding points S r (ξ ξ ξ) and S c (ξ ξ ξ) measured along the normal n Ωr . Where the normal n Ωr of the reference body surface is a normalized vector that is normal to the surface at a given point. It is computed from the cross product of the partials ∂/∂ξ 1 and ∂/∂ξ 2 at that point t Ωr,ξ1 = ∂ ∂ξ 1 S r (ξ ξ ξ) and t Ωr,ξ2 = ∂ ∂ξ 2 S r (ξ ξ ξ) (9) A programmable formula for implementation of NURBS surface derivatives can be found in Tiller, 1998, 1995). The basis vectors {t Ωr,ξ1 , t Ωr,ξ2 } which span the tangent space at the relevant point S r (ξ ξ ξ) can be used to define a local Cartesian basis. This new vector basis denoted by {ê ξ ,ê ζ ,ê η } is defined such that the first basis vectorê ξ is parallel to t Ωr,ξ1 , andê η is orthogonal to it. Moreover, it is oriented such that the second axis, that isê ζ is collinear with the normal n Ωr : e η = t Ωr,ξ2 − (t Ωr,ξ2 ·ê ξ )ê ξ t Ωr,ξ2 − (t Ωr,ξ2 ·ê ξ )ê ξ This local Cartesian basis, see Figure 4, facilitates coupling between the multibody-system formulation and the conservation laws defined on the fluid domain. The transformation matrix representing this new frame with respect to the surface frame O sr x sr can be expressed as where the correlation factor κ is utilized to ensure the normal n Ωr points towards the corresponding body, such that κ = −1, ifê ζ points into the reference body 1, ifê ζ points out from the reference body By expressing the coordinates of point S c (ξ ξ ξ) in the frame O sc x sc in homogeneous coordinates, the gap height vector H ∈ R 3 connecting the points S r (ξ ξ ξ) and S c (ξ ξ ξ) via the normal n Ωr is defined by where 0 is just the zero point (0, 0, 0) ∈ R 3 . Furthermore, note that the fourth component is zero in the gap height vector [H 0] T ∈ R 4 , because this is a vector and not a point. By post-multiplication of the transformation matrices in Equation (15), the gap height vector reads Conclusively, in order to find the gap height vector H in this generic formalism one must know the corresponding surface point S c (ξ ξ ξ). This is not a trivial task since the surface of the corresponding body designated Ω c is defined by a NURBS parameterization.

Finding the Corresponding Surface Point on Ω c
In order to find the corresponding point a few steps need to be taken. The basic idea is to convert the NURBS surface of the corresponding body into defining implicit matrices called an "implicit matrix representation" of the NURBS surface. The implicit form of the surface permits a fast algorithm for computing point/surface intersection. The general algorithm is shown in Figure 5.  The algorithm starts by decomposing the NURBS surface of the corresponding body into a number of Bézier patches. The reason for this decomposition is that it simplifies the implicitization process immensely by switching from NURBS basis functions to Bernstein polynomials. The Bézier patches are then implicitized by taking the implicit matrix approach of Busé (2014), which produces an implicit matrix representation (Mreps for short) of each of the Bézier patches. Lastly, a simple intersection search algorithm is used to find the corresponding point. Each of these steps are described in the following.

Bézier decomposition
Before the NURBS surface is implicitized it is first subdivided into Bézier patches by applying knot insertion. Bézier decomposition is accomplished upon repetition of all interior knots of a knot vector up to multiplicity equal to the order of the knot vector (d + 1).
Given an open knot vector Ξ = {ξ 1 , . . . , ξ m } of degree d, with the bounds ξ 1 , . . . , ξ d+1 = 0 and ξ m−d , . . . , ξ m = 1 and let ξ ∈ [ξ k , ξ k+1 ) be a new knot that is already present in the original knot vector Ξ. Although, it could as well be a knot that is not already in the original knot vector. Then the (n + 1) new B-spline basis functions N i,d (ξ) can be computed from Equation (3) and (4), with the new knot vector Ξ = Ξ {ξ}. Using projective control points P w i = w i (x i , y i , z i , 1) (for efficient processing), the (n + 1) new projective control points Q w i are computed from the original P w i by where By computing the new control points (Bézier control points) using the formula in Equation (17), the knot insertion corresponds to a change of vector space basis, i.e. the surface is not changed. In order to get back to the three-dimensional Euclidean space each component of the Q w i is divided by the weight w i . Details on the Bézier decomposition can be found in Piegl and Tiller (1995).
and it is said to be in degree ν. Each entry Γ i,j in the M-rep is a linear polynomial in the variables x, y, z.
A M-rep as defined in Equation (23) has the following drop-of-rank property: * If ν ≥ ν 0 , then r ν ≥ m ν , such that the rank (rk) of the M-reps evaluated at any point x = (x, y, z) is rk(M ν (x)) ≤ m ν ; * If ν ≥ ν 0 , then rk(M ν (x)) < m ν if and only if x is a point in the collection of points Ω in the surface.
The above properties constitute the drop-of-rank property, which validates the use of M-reps as an implicit representation of a parametric surface. The proofs are not presented, but the interested reader is referred to Busé (2014) and the references therein. Now by letting the 4-tuple of polynomials (g 0 , g 1 , g 2 , g 3 ) ∈ P ν (ξ) 4 be expressed in the tensorproduct Bernstein basis B ν = {β i (ξ 1 )β j (ξ 2 )} ν1,ν2 i=0,j=0 , the polynomials (g k (ξ ξ ξ)) 3 k=0 can be written as To find the polynomials (g k (ξ ξ ξ)) 3 k=0 that satisfy Equation (22), the problem is converted to the component space under appropriate choice of the Bernstein basis. For this reason, the matrix S ν is introduced, which is the matrix representation of the linear map 3 i=0 g i (ξ ξ ξ)f i (ξ ξ ξ) from Equation (22), see diagram in Figure 6 for clarification. The matrix S ν must satisfy the following matrix equality relation Where B ν and B ν+d are row vectors. The matrix S ν is thus a coefficient matrix, which decomposes the right-hand side with respect to the basis B ν+d = {β i (ξ 1 )β j (ξ 2 )} ν1+d1,ν2+d2 i=0,j=0 . Consequently, the dimension of S ν yields dim(S ν ) = (ν 1 + d 1 + 1)(ν 2 + d 2 + 1) × 4m ν where m ν = (ν 1 + 1)(ν 2 + 1). To compute the matrix S ν , let f (ξ ξ ξ) be a parametric tensor product function, defined by Here c i,j is just a placeholder variable for the Bézier weights w i,j and weighted control points w i,j Q i,j . In order to derive a programmable formula for implementation, it can be useful to introduce global indicesĨ,J. These are expressed as where k and l are any pair of integers for which 0 k ν 1 and 0 l ν 2 . Then, the matrix S ν can be filled according to the formula The coefficients of the polynomials (g k (ξ ξ ξ)) 3 k=0 that satisfy Equation (22) are equivalent to the vector components of all the vectors γ γ γ that lie in the null space of the matrix S ν , that is For each free variable x i1 , . . . , x ir ν of S ν , γ γ γ j is the i j 'th column of Null(rref(S ν )) 1 plus the i j 'th unit vector. Thus, the null space of S ν is spanned by this set of linearly independent vectors γ γ γ j , Null(S ν ) = span(γ γ γ 1 , γ γ γ 2 , . . . , γ γ γ rν ) The M-rep is then given by the formula, Where each block matrix M ν,i for i = 0, . . . , 3 is built from the components of the vectors γ γ γ ∈ R 4·mν spanning the null space of S ν .

Intersection search
Having computed the M-reps of the Bézier patches, the points on the surface can be found by directly evaluating M ν (x) in Equation (32). The points on the surface corresponds to the points for which M ν (x) has not full rank, following from the drop-of-rank property. A more targeted approach to find the corresponding point S c (ξ ξ ξ) is to set up a minimization problem, where the real evaluation function from Busé (2014) is used as objective function. The real evaluation function is defined as: Definition 2. Let R 3 be a vector space, then the real evaluation function of M ν at x ∈ R 3 is defined by where σ i (M ν (x)) are the singular values of the Mreps M ν (x) ∈ R mν ×rν . The real evaluation function has the following properties, * ∀ x ∈ R 3 : D(x) ≥ 0; * D(x) = 0 if and only if x ∈ Ω ⊂ R 3 .
The first property of the above implies that the real evaluation function is non-negative for every point x ∈ R 3 . The second property is a direct implication of the drop-of-rank property of M ν , which states that the real evaluation function is zero if and only if the point of evaluation belongs to the point set Ω representing the surface geometry.
A point x ∈ R 3 between the corresponding points S r (ξ ξ ξ) and S c (ξ ξ ξ), see Figure 2, is given by the relation Where the gap height vector H is given with respect to the basis {ê ξ ,ê ζ ,ê η }, thus it can be written as H = (0, h, 0). This means, that the function evaluation can be performed by simply varying the one-dimensional parameter h in Equation (35). The gap height vector H can be found by solving the following minimization problem, arg min where the variable h is constrained to the domain of positive reals. The minimization problem can be solved numerically using the golden section search algorithm cf. Arora (2012). To bracket the minimum h min 2 , the function D(x) is evaluated and compared at a sequence of points x. The sequence of points x is found by incrementing the variable h in Equation (35) where δ is a small number. Having found the lower and upper limits h L and h U , respectively, spanning the interval of uncertainty I, the interval can be reduced upon calculation of the intermediate points h a and h b , which are located on the interval at a distance of 0.382I from either end, see Figure 7. More details and algorithmic implementation of the golden section search can be found in Arora (2012). Note that both of these wall velocities are expressed in the reference body frame on the basis that the lubrication surface is coincident with the reference body surface, wherein hydrodynamic conservation laws are expressed.

Conclusion
A general methodology for tribological kinematics has been presented. Explicit expressions for the gap height distance and velocity and the wall velocities in a general clearance joint are provided. The whole framework is developed in a NURBS setting, thus the suitability of the proposed method is entirely based on the assumption of a readily available CAD model. This allows an exact geometry representation of the tribological interface and ability to consider more complex joint shapes. Moreover, it reduces the time and effort needed to manually parameterize the interface geometry. The approach can be easily extended to the whole dynamics analysis of systems of rigid bodies using the Newton-Euler method, as the kinematics framework is described in terms of Cartesian coordinates. Thus, the proposed kinematics approach is considered a key component in the development of a generic multibody tribo-dynamic simulation tool. Such a simulation tool is potentially beneficial in a wide-ranging design analysis taking into account various machine topologies suitable for some targeted application system. In that regard, but also tribodynamic system modeling in general, the computational effort is a fundamental issue. The generality of the proposed method possibly adds to the computational cost and is therefore considered a drawback in relation to computational efficiency. To improve practicality one must consider code optimization and parallel computing.