A linear complementarity problem formulation for periodic solutions to unilateral contact problems

Presented is an approach for finding periodic responses of structural systems subject to unilateral contact conditions. No other non-linear terms, e.g. large displacements or strains, hyper-elasticity, plasticity, etc. are considered. The excitation period due to various forcing conditions—from harmonic external or contact forcing due to a moving contact interface—is discretized in time, such that the quantities of interest—displacement, velocity, acceleration as well as contact force—can be approximated through time-domain schemes such as backward difference, Galerkin, and Fourier. The solution is assumed to exist and is defined on a circle with circumference T to directly enforce its periodicity. The strategy for approximating time derivative terms within the discretized period, i.e. velocity and acceleration, is hence circulant in nature. This results in a global circulant algebraic system of equations with inequalities that can be translated into a unique linear complementarity problem (LCP). The LCP is then solved by dedicated and established methods such as Lemke’s algorithm. This allows for the computation of approximate periodic solutions exactly satisfying unilateral contact constraints on a discrete time set. The implementation efficiency and accuracy are discussed in comparison to classical time marching techniques for initial value problems combined with a Lagrange multiplier contact treatment. The LCP algorithm is validated using a simple academic problem. The extension to large-scale systems is made possible through the implementation of a Craig–Bampton type modal component synthesis. The method shows applicability to industrial rotor/casing contact set-ups as shown by studying a compressor blade. A good agreement to time marching simulations is found, suggesting a viable alternative to time marching or Fourier-based harmonic balance simulations.


Introduction
The investigation of periodic solutions of harmonically excited systems subject to unilateral contact constraints is of interest for many applications, such as predicting brake disk vibrations [1], stability analysis of delayed systems subject to material removal [2], study of limit cycles in Lur's feedback control systems for non-smooth mechanical systems [3,4] or mesh stiffness variation in gears-pairs [5,6,7]. Popularly, methods such as harmonic balance based [8] or shooting methods [9] are employed to find solutions of smooth nonlinear systems in an efficient manner. Yet, both of these approaches face some inherent downfalls when nonsmooth nonlinearities-systems exhibiting non-differentiability or discontinuities in the unknowns-are encountered. On the one hand, the harmonic balance method (HBM) is known to produce poor approximations of non-smooth functions with a finite number of harmonics, producing artifacts such as the Gibbs phenomenon [10]. Hence penalty-like approximations of the contact inequalities are introduced [1] to effectively smoothen the nonsmoothness. On the other hand, shooting methods in a contact framework can face ill-conditioned gradients, making convergence difficult and increasing computational times. Furthermore, reaching a purely steady state solution within a heavy-duty timemarching simulation may be impossible to achieve due to the high sensitivity of the solution with respect to system parameters, such as stiffness, frequency, gap.
The so-called linear complementarity problem (LCP) deals with a linear algebraic system of equations subject to inequality constraints. Originally presented in [11,12,13] as an alternative formulation for quadratic programming, a few notable solution methodologies have been introduced: the classical pivoting, also generally referred to Lemke's algorithm [11], and iterative methods such as the Gauss-Seidel method [14,15] and Newton-like methods [16] are among the most popular. With applications in many scientific fields, the LCP has found a strong foothold especially in economic engineering, mathematical programming, games theory and recently switched electronic systems [3,17,18]. The idea is that LCP solvers can find solutions of underlying linear algebraic systems subject to nonlinear Kuhn-Tucker-like conditions. The LCP method has already been applied to piecewise linear mass-spring systems in a time-marching framework, solving an LCP for every time step [19]. In a similar fashion, the application of LCPs to transient analysis frictional problems is reported in [20]. This paper presents an approach for finding periodic solutions to systems of ordinary differential equations (ODE) involving inequalitities within a unique LCP. Firstly, the general LCP formulation is outlined. Next, different time-derivative approximations are detailed for constructing the LCP system. An academic application is presented, specifically looking at a comparison of the LCP results to classic time-marching simulations. Finally, the application of the LCP method to an industrial compressor-blade geometry is explored with a detailed focus on frequency domain analysis of the responses.

Linear complementarity problem
The formulation of an LCP lends itself explicitly to treat linear mechanical systems subject to Kuhn-Tucker conditions such as unilateral contact. Forced and possibly large-scale mechanical systems undergoing unilateral contact conditions that are Tperiodic are targeted, see Figure 1. A periodic displacement of the model discretized in space x.t/ of period T in time is assumed to exist. Accordingly, it is a solution to the following combined equations and inequations: where d.t/ is the reference position of a potentially time-dependent moving rigid wall. The impenetrability condition is expressed by inequality constraints on the gap function g.x.t// and contact force .t/ as well as the complementarity criteria in Eq. (1b). It should be noted here that the gap function is considered to be purely linear in x.t/ as expressed by Eq. (2). Systems for which the constraint matrix varies within the period of interest, e.g. large tangential displacements or strongly deformation dependent contact interfaces, cannot be accounted for within this approach. In order to simplify the computation, contact between a flexible body and a rigid wall is assumed, although flexible multibody contact problems do not pose a mathematical limitation and the equations explicitly remain the same To transform Eq. (1a) and (1b) into an equivalent LCP able to capture periodic solutions, for a given period T , a few steps are necessary. Firstly, a time-discretization is needed. Assuming an excitation period of T and hence a frequency of ! D 2=T , the support S T is discretized into k-time steps by a linear spacing of t D T =k, such that x.t/ fulfil periodicity 8t 2 S T , see Figure 2. The displacement along S T is then approximated The corresponding LCP, detailed in the next sections, is expressed as follows: find complementary vectors w and z such that: w D Az C q; w 0; z 0 and z > w D 0 where the inequalities shall be read coordinate-wise.

Formulations
The following section deals with different approaches on how to approximate the time derivative terms P x and R x within EOM (1a). Firstly the backward difference scheme (BDS) is recalled. Alternatively, a more complex approximation is developed, based on a Galerkin-like approach in time. Finally, the use of a Fourier basis is discussed, which represents a high dimension HBM (HDHBM) as proposed in [21] and presented in a contact framework in [22], specifically considering non-regularized contact conditions.

Backward difference scheme approximation
The simplest and most straightforward approximation of the time derivative terms in Eq. (1a) is based on a first order BDS in velocities and accelerations: together with the periodicity conditions in velocity and acceleration, respectively: where x i is the time discretized displacement over a time step t . This scheme represents the simplest approximation of time derivative terms, able to capture non-smoothnesses easily. Although, these approximations have a relatively high error of O.t / as compared to e.g. a central difference scheme (CDS), a more accurate CDS or higher order method would be counterproductive, as these methods only have a high accuracy for smooth functions and tend to produce oscillations within solutions of non-smooth systems [23,24]. For a single DOF, the linear differentiation operator for the BDS approximation is a circulant matrix that reflects timeperiodicity in the following: The operator L 1;t can be applied to the time discretized j -th DOF in the shape of O The hat symbol denotes that the variable is the time-discretized representation of its time-domain counterpart. It should be noted that the top right matrix entry in L 1;t is due to the imposed periodic boundary condition in time. A second order operator (7) can be applied to acceleration terms analogously: Expanding (1a) to the time discretized system using Eqs. (6) and (7), the respective matrices need to be extended with the use of the Kronecker-product: Equation (1a) can then be written in the following form: with: O D OE 1 1 ; 1 2 ; : : : ; 1 k ; 2 1 ; : : : in which the vectors O x, O and O f are not time-dependent anymore. Solving Eq. (9) for the displacement results in: Combining (11) and (2) results in an algebraic set of equations in the form of (4) where: For the mechanical systems considered within this documentation, w 2 R mk corresponds to the contact node displacements, z 2 R mk to the contact forces. The vector q 2 R mk represents the solution of the underlying linear system due to external forcing without contact constraints. A 2 R mkmk corresponds to the coefficient matrix, mapping the contact forces onto the constrained response.

Finite element time discretization
A more flexible approach than the BDS is found by using a finite element in time discretization (FETD) of the quantities of interest. In this framework, various approximations for velocity in terms of displacement and acceleration in terms of velocity can be derived in an integral sense as follows:

8‰.t/;
such that the equation of motion becomes: For example, the displacement x.t/ and weighting functionˆ.t/ will be expressed as follows: where x i are the time-discretized vectors of x.t/. Proper combinations of well-chosen shape functions to capture potentially non-smooth displacement, velocity and acceleration fields yield accurate solutions. Figure 3 shows the two possible shape functions applied in this investigation. It was found that the as- Figure 3.  Table 1 produce a good agreement with converged time-marching simulations. Further developments and improvements may yield better agreement. The assumptions Table 1. FETD shape functions made for (13) correspond to a higher order averaging on the side of velocity and a simple central difference representation of the displacement terms for a single time step: Equation (14) including the above assumptions is analogously reflected in the following: In a similar fashion, solving the equation of motion (15) using the assumptions from Table 1 results in: The linear operators L a , L b and L c , similar to (6) and (7) in spirit, can be constructed through Eqs. L a;t represents a simple central difference operator; L b;t as well as L c;t stem from an averaging due to the above assumptions. Effectively, the averaging terms provide an increased coupling of the discrete time-responses and hence smoothen out displacement and contact force oscillations which allows for a better approximation of the contact constrained response. Analogous to (8) these operators need to be expanded to the full spacial and temporal system using the Kronecker-product: The system of equations is in the form of (9) which can be translated into (4) analogously, where:

High dimension harmonic balance method
The applicability of the LCP formulation is also tested by applying a HDHBM approach on the variables of EOM (1a). The displacement is assumed to be a function of the form (16) of a truncated Fourier series with p harmonics: In a fashion similar to Eq. (6), linear operators can be constructed in the frequency domain to represent velocity and acceleration approximations. For a single DOF j in the frequency domain of the form Q x j D OEa j 0 ; a j 1 ; b j 1 ; : : : a j p ; b j p > the first order derivative operator would be the following: The second order derivative operator is computed by Following Eqs. (8) and (21), the full expanded system takes the form of (9) in the frequency domain, where: and As the inequality constraints are explicitly dealt with in time rather than in the frequency domain, a direct Fourier transform (DFT) and inverse direct Fourier transform (IDFT) need to be applied. By assuming k D 2pC1 time steps, square DFT (F) and IDFT matrices (F 1 ) can be constructed translating the frequency domain coefficients into the time domain DOFs without loss of information, thus resulting in (4), where

Implementation
Since the computational effort increases exponentially with the size of the problem, a few numerical reductions are introduced. The inversion of O Q in (11) is computationally highly inefficient. To circumvent this operation it is beneficial to numerically solve Utilizing the sparsity of O Q and the fact that generally m n a sufficiently fast implementation can be obtained.

Craig-Bampton reduction
For larger systems, e.g. the compressor blade finite element model discussed in later sections, the number of DOF can easily reach the order of 100,000 or more. The computational effort to reduce this problem to the contact DOF as mentioned before is prohibitive, as the reduction needs to be performed for every excitation frequency of interest. By incorporating a Craig-Bampton type modal component synthesis (MCS), the linear model dynamics can be reduced to a synthesis of free static modes and a truncated number of component modes [25]. The following paragraphs outline the MCS procedure. Damping is not addressed before the MCS is applied, but rather a modal damping is used on top of the MCS reduced system.
To efficiently handle boundary x b DOF-DOF that lie on the contact interface-and internal x i DOF separately, it is beneficial to rearrange the system matrices into these groups. The boundary DOF x b will be kept explicitly in the reduced model, while internal DOF x i are represented implicitly by a set of component modes: At this point generally the contact constraint matrix B is created. It selects the DOF on which contact conditions are applied, depending on the local contact interface definition, i.e. contact normal, wall function, etc. The rearranged system in terms of internal and boundary DOF mass and stiffness matrices can be written as follows: Next the internal free vibration eigenvalue system is solved, whereˆi are the internal mode shapes and ƒ i contains the internal eigenvalues on its diagonal: For the boundary DOF that are to be kept explicitly in the reduced model, static modes are computed as follows: The final MCS reduction matrixˆc b can then be written as: The mass, stiffness and contact constraint matrices are then reduced as M cb Dˆ> cb Mˆc b ; K cb Dˆ> cb Kˆc b ; and B cb Dˆ> cb B The DOF vector of the reduced system contains explicitly x b as well as component mode amplitudes a i : x cb D OEa i x b > . The resulting system can be multiple orders of magnitude smaller than the original finite element model, by considering a truncated set of modes from the internal eigenvalue set ƒ i < max . Convergence must generally be tested depending on the model, to see if the lowest system eigenmodes are represented accurately. The reduced system can now be translated into the aforementioned LCP system, while significantly reducing the computational effort in solving the linear systems. The acquired results can be translated back into the full finite element space utilizing the MCS reduction matrix as x Dˆc b x cb .

Stability and bifurcations
A particularity of non-linear systems is the possibility of bifurcations of solutions. Furthermore, the existence and uniqueness of periodic solutions is not guaranteed. The found solution hence depends on the solution methodology used within the LCP algorithm itself. Throughout this documentation the LCP solver is based on Lemke's classical pivoting algorithm. The determination of the solution the algorithm will find is beyond the scope of the presented study.

Comparison
The aforementioned formulations for approximating the time derivatives within the LCP framework are compared to classic time-marching results. An explicit time-marching scheme (ETM) is utilized, where contact forces are computed using a one-step Lagrange multiplier methodology [26]. Three distinct excitation frequencies have been chosen and applied to all methodologies, in order to compare resonant and non-resonant conditions.

Rod model
To efficiently compare accuracy and performance of the different methodologies, a small scale academic non-dimensionalized problem is considered. A 10-DOF mass-spring system (simplified rod model) is used considering axial deformation only, see  The tip of the rod is subject to contact constraints and is excited due to a harmonically moving wall d.t/. External forcing is neglected, but rather the system responds solely due to contact interactions. The elementary matrices are as follows:

Simulation cases
The first case considers the excitation frequency of the moving wall to be 40% of the first natural frequency of the rod model. At this frequency participation of multiple modes is expected.
To study the response just below linear resonance, the second case studies the system response at a frequency of 90% of the first rod eigenfrequency. The final case explores the responses at an excitation that is above the first eigenfrequency at 130%.

Computation
The equation of motion is normalized with respect to the first eigenvalue ! 2 0 . The time-marching computations are executed using 40,000 time steps per period, over 100 periods. For the three chosen excitation frequencies, all solutions reach steady state within 100 cycles. A time-step convergence study has been performed prior, but is left out for the sake of brevity. The timemarching solutions are assumed to be converged and shall be used as reference solutions. The LCP computations, for both the BDS and FETD schemes are computed over 500 time steps per period. Close to resonant conditions, e.g. the ! D 1:3! 0 case, the LCP solution for both FETD and BDS using Lemke's algorithm takes notably longer (up to two order of magnitude). For these conditions, the LCP implementation used is on par with the time-marching computations to reach steady state in terms of computational time. Otherwise LCP computations run one to two order of magnitude faster than the classical time-marching strategy.

Solutions and error analysis
The contact node displacement x c over one period is plotted in Figs. 5(a), 6(a) and 7(a) for the three speeds respectively. The corresponding contact forces for the same speeds are plotted in Figs. 5(c), 6(b) and 7(b) respectively, focusing more detailed on the time interval where contact occurs, marked in gray.
For the sub-resonant case, a good agreement between the ETM and FETD solutions can be seen in the displacements in Figure 5(a), the velocities in Figure 5(b) and contact forces in Figure 5(c). The BDS solution, displays a strongly damped solu-  tion, especially in high frequency content. The HBM performs well in the free vibration interval, but has strong problems capturing the switching conditions on contact forces accurately. Rather than a smooth attachment of the contact node to the moving wall, the HBM predicts two impacts without attachment. Refining the resolution in time steps and plotting the contact forces 5(d) shows that even though the LCP satisfies the inequality constraints at the user-defined k control-points, the Gibbs phenomenon appears around the contact interactions. Hence Fourier-functions provide an ill-chosen basis within this framework.  where the integral can be exactly calculated through the timedependent shape functions. As expected the Fourier basis approach does not show a converging behaviour with more harmonics. The BDS and FETD approaches both show a similar convergence rate for small time steps. The convergence rate suggests that refining the time-step even further, would yield more computationally accurate results for both BDS and FETD approximations, yet these are subject to computational restrictions as discussed before. For the other excitation frequencies (! D 0:9! 0 and ! D 1:3! 0 ), the LCP solutions differ increasingly from the ETM result. Although the BDS approach is generally outperformed by its FETD counterpart, its efficiency and simple implementation are beneficial. Both the BDS and FETD implementations are able to capture dominant interaction dynamics, which is essential in determining critical excitation frequencies. As the applicability of the LCP approaches has been shown for a simple academic model, it is of interest to see how efficient the aforementioned methodologies are when applied to an industrial problem.

Industrial application
This section discusses the application of the LCP formulations to a modern compressor blade. The goal is to simulate a single compressor blade rotating over a given speed range, while subject to contact interactions with a rigid casing [27]. Any interaction phenomena with neighbouring blades as well as centrifugal effects are neglected. This does not explicitly present a limitation of the LCP methodology. Similar to the previous case presented in section 4, the equation of motion is normalized with respect to the first eigenfrequency ! 2 0 . The studied blade geometry represents the size and shape of a modern compressor blade geometry. The mesh consists of 11,000 quadratic tetrahedral elements, resulting in 69,000 DOF. The typical mesh associated with the blade geometry is shown in Figure 9. In the MCS the first 24 modes are preserved, and in addition to the 3 DOF per contact node, resulting in 33 DOF for the reduced model. The contact of the blade with its surrounding casing is modeled by a radially moving wall-similar the the rod model set up-periodically contacting three chosen interface contact nodes at the blade tip: at the leading edge (LE), mid-chord (MC) and trailing edge (TE), see Figure 9. The speed-range, and hence excitation frequency, is chosen to be ! 2 OE0:1; 0:8, since rotational speeds in high-pressure compressors generally do not exceed the first blade eigenfrequency. Analogous to the rod model, the system responds purely due to contact interactions, rather than linear external forcing. Post-processing is performed on the frequency domain signal of each simulation, obtained through the FFT of the contact node displacement.
The excitation frequency range is discretized into 70 equally spaced frequencies i 2 OE0:1; 0:8. After simulating each discrete excitation frequency in ETM and both BDS and FETD LCP approaches, and performing an FFT of the time signals considering ten periods, Figs. 10(a), 10(b) and 10(c) are plotted. Displayed are the harmonic response coefficient amplitudes of the LE contact node displacement x c in color corresponding to their response frequency versus the excitation frequency. The inclined dashed lines represent integer multiple harmonics with respect to the excitation frequency. In turbo-machinery applications these are generally referred to as Engine Orders (EO). EO D 1 corresponds to a response occurring at the excitation frequency, i.e. the engine speed. Higher EO, eg. EO D 2; 3; 4; : : : correspond to responses occurring at integer multiples of the excitation. The horizontal solid line across all excitation frequencies at f =! 0 D 1 corresponds to the first natural frequency of the underlying linear structure. and ! e D 0:7. The response at ! a , even though it occurs at a EO crossing, shows high amplitudes for all response frequencies. This is due to the fact that the time marching algorithm does not reach a periodic steady state. At ! e a peak response is visible on EO-lines as well as half EO-lines, see e , suggesting a subharmonic response, e.g. requiring two excitation periods to complete one response cycle. Both of these phenomena are exclusive to the time marching computations, where no explicit enforcing or periodicity is given. Due to the single-period periodic nature of the results stemming from the LCP implementations, sub-harmonic or nonsteady state solutions are not reflected in both BDS and FETD results. This is indicated by the zero magnitude response, in between EO-lines. Hence the response peaks lie purely on EOlines. Figure 10(b) shows two distinct peaks that are also visible in the ETM frequency map, i.e. ! c D 0:35 and ! d D 0:53, see c and d . The peak at ! b that is visible for time-marching is about two orders of magnitude smaller in the BDS computations. The FETD computations is able to capture the interaction at ! b more accurately.
Plotting the Fourier coefficient amplitudes for the range of discret rotational speeds/excitation frequencies !, considering a single period for the first 10 harmonics, results in Figure 11. Figure 11(a) shows the harmonic content of the time-marching solutions. Figures 11(b) and (c) show the harmonic content of the LCP solutions for BDS and FETD discretizations. The non-steady state reaching time-marching solution at ! a D 0:21 ( a ) is clearly visible in 11(a), yet it has only little qualitative information as it is not periodic. The same is true for the subharmonic response e at ! e D 0:7. Appart from these, the peaks b , c and d show distinct responses on the 4th, 3rd and 2nd harmonics respectively, corresponding to their respective EO crossings.
The overly damped nature of the BDS implementation comes to light in Figure 11(b). While showing minimal response amplitudes for most excitation frequencies, only the periodic re-  There are no notable contributions above the 4th harmonic, as it seems that higher harmonic content is damped out. The FETD implementation of the LCP on the other hand is able to capture the three periodic responses clearly. Not only do higher harmonics still contribute to the response, but also the main response harmonics respond on a far higher amplitude. Figure 10(d) shows the frequency domain results considering two full excitation periods within the LCP formulation. This consideration is able to capture sub-harmonic results of order EO D 1=2. The crossing e appears within Figure 10(d) in a similar fashion to Figure 10(a) at ! e D 0:7. As two periods are considered, even though the time-step size is to be maintained, the computational effort is significantly larger when trying to find sub-harmonic responses, especially for high order sub-harmonics.

Conclusion
A methodology is presented for finding periodic solutions of externally forced mechanical systems subject to unilateral contact constraints. By formulating the contact constrained equation of motion in a linear complementarity problem form and solving the system for an entire time-discretized period in a unique set of algebraic equations, an efficient strategy is developed. Three approaches of dealing with time derivative discretization are presented: a backward difference, a finite element in time and a Fourier based approach. The backward difference approximation shows strong artificial damping for a small number of time steps per period, unable to capture but the most dominant contact interactions. Using the finite element approach with the same number of time steps, these interactions are captured more accurately. Both the backward difference and finite element in time approaches show an asymptotic convergence with decreasing time step size, that perform similar for both but are subject to computational constraints. The Fourier based approach produces undesirable oscillations around switching conditions, generally known as the Gibbs phenomenon, thereby rendering the harmonic basis ill-chosen. A strong limitation exists on the size of the problem, and the number of time steps used per period. For small problems with DOF of order 10, computational limitations of about 1,000 time steps per period are experienced, making time step convergence difficult. . The performance and numerical limitations of the proposed method are highly dependent on the solver, as well as hardware and software used to perform the simulations. Great improvements on the problem size and time-step resolution are expected to be possible by optimizing the solver as well as the set-up of the LCP system.
An application of this methodology to an industrial case is documented. The contact interactions of a modern compressor blade-like geometry with its surrounding casing are considered. Over the predetermined speed range, multiple interactions are encountered using all three methodologies. Within the timemarching framework, solutions that do not reach a periodic steady-state as well as sub-harmonic solutions are found, that do not appear in the linear complementarity problem solutions, due to the single-period periodic assumptions. Solving a linear complementarity problem that imposes periodicity over multiple periods is able to capture sub-harmonic solutions, yet is even further constrained in the number of time-steps required.
Due to the strict conditions on periodicity within the LCP strategy, not all phenomenon during contact interactions can be accounted for. At this point time-marching algorithms may be more physically accurate. Yet, the LCP is able to capture the dominant dynamics at critical excitation frequencies. More sophisticated approaches in the description of time derivative terms may increase performance of the algorithm, allowing for a fine time-step size and hence more accurate solutions. Overall, the presented work shows the potential for efficient and fast computations of non-smooth periodic contact-constrained solutions with no artificial implementation of common penalty-based strategies.  Figure 11. Harmonic amplitudes for ETM, BDS and FETD methods