NONSMOOTH MODAL ANALYSIS WITH BOUNDARY ELEMENT METHOD

. Numerical schemes based on the Boundary Element Method are proposed to perform Nonsmooth Modal Analysis. The latter aims at ﬁnding continuous families of periodic orbits of mechanical systems featuring unilateral contact constraints. In this contribution, a simple one-dimensional rod system is targeted. The frequency response, in the form of backbone-curve diagrams, and displacement ﬁeld are presented. The proposed results compare well with existing studies on this topic.


INTRODUCTION
Within the framework of structural dynamics, linear modal analysis is a daily used tool in industry, aiming at predicting vibratory resonances of periodically forced mechanical systems by searching for continuous families of periodic solutions exhibited by the underlying autonomous (i.e. unforced) system. However, various challenges arise when possibly large-scale nonlinear dynamical systems are targeted and for which nonlinear modal analysis is needed instead [7]. Nonlinearity has commonly two forms: smooth or nonsmooth function of the state. The nonlinearity is said to be smooth when the governing equation involves a function of the state which is differentiable. Such functions are commonly polynomial. The nonlinearity is said to be nonsmooth otherwise, as for instance, unilateral contact conditions. Nonsmooth Modal Analysis (NSM) is a version of modal analysis dedicated to such systems. A few numerical schemes exist to perform NSM, with inherent limitations [11,12]. In the present work, it is proposed to explore the numerical capabilities of two numerical schemes grounded on the Boundary Element Method (BEM) to perform NSM.

SYSTEM OF INTEREST
The system of interest, shown in Figure 1, is a one-dimensional bar of length L with a unilateral contact condition at its right tip. The unknown displacement field is u.x; t /, where u.x; t/ L Figure 1: System of interest x is the space coordinate and t is time. Wave speed c is space-independent and the classical one-dimensional wave equation governs the dynamics, that is Unilateral contact on the right tip is expressed as a Signorini boundary condition g.t/ 0; u ;x .L; t / Ä 0; g.L; t /u ;x .L; t / D 0: where g.t/ D g 0 u.L; t/. Also, since the bar is clamped on the left, u.0; t / D 0 applies. In context of nonlinear modal analysis [7], families of periodic solution forming modal manifolds are computed through a numerical scheme which assumes that no impact law is required as the contact interface but instead, a switch between free-flight and sticking phases occurs.

METHODOLOGY
In this paper, two numerical schemes based on the Boundary Element Method are implemented. The BEM forms a family of methods for which the boundary of the domain of interest is the main ingredient of the formulation and full domain discretization is not required when initial condition and body force vanish [6]. It has the notable benefit of reducing the dimension of the formulation. Among the various incarnations of BEM, the Time Domain BEM (TD-BEM) and Frequency Domain BEM (FD-BEM) are selected to perform nonsmooth modal analysis. The first strategy combines TD-BEM and a simplified shooting method, while the second one combines FD-BEM to Harmonic Balanced Method (HBM).

Time Domain Boundary Element Method and Shooting Method
TD-BEM is a version of BEM which requires the discretization of both time and space. By pre-multiplying the governing equation (1) with its fundamental solution [4] and performing integration in time and space and then integration by parts, the formulation can be transformed into Boundary Integral Equation (BIE), which here reads [3] u.x; t/ D where p represents traction on the boundary ; u .x; t; ; / and p .x; t; ; / are the displacement and stress fundamental solutions. Any internal state for x 2 0 I LOE at time t is dependent on the boundary displacement and traction at past times as well as initial states. In order to use TD-BEM, space and time are first discretized. Space integration of the initial states (last two terms in (3)) is performed by discretizing space into n elements and n C 1 nodes with x D L=n such that x k D kx for k D 0; : : : ; n. Time step is set to t D x=c to guarantee stability of TD-BEM, such that t i D it for i D 0; : : : ; n. Equation (3) is discretized by reading it on the boundary, that is x D 0 and then x D L for all t i for all i D 0; : : : ; n; it now becomes where u i OEu.0; t i /; u.L; t i / > and p i OEp.0; t i /; p.L; t i / > are displacement and traction vectors on boundary, with size of 2 1. The quantities H i i , G i i , H ij , and G ij are 2 2 coefficients matrices, evaluated from time-domain integration in (3) of fundamental solutions u and p . Vectors u 0 D OEu.x k ; 0/ kD0:::;n and v 0 D OEv.x k ; 0/ kD0:::;n are the discretized initial displacement and velocity fields, with size n C 1. The corresponding coefficient matrices of size 2 .n C 1/ are U i and V i stemming from space integration of fundamental solutions in (3). For systems without Signorini boundary conditions, Equation (4) can be solved at every time step starting at i D 1 where half of the quantities in u i and p i are unknown, that is .u.0; t i /; u.L; t i // or .u.0; t i /; p.L; t i // or .u.L; t i /; p.0; t i // or .p.0; t i /; p.L; t i //.
In order to account for unilateral contact conditions, TD-BEM is coupled with the floating boundary method [8]. It sees Signorini boundary condition as a switch between nonhomogeneous Dirichlet and homogeneous Neumann boundary conditions. This technique does not require an impact law in the BEM format, thus avoiding the chattering or energy dissipation that might exist in other methods [1].
To perform NSM, a shooting method is used. It calculates, here in the discretized framework, the difference between the initial state q 0 D .u 0 ; v 0 / and the corresponding final state q.q 0 ; T / after a yet unknown period T assumed to exist. For a periodic motion, the equality q.q 0 ; T / D q 0 holds. In this work, v 0 D 0 is assumed [11], and the problem to be solved thus reduces to: Find vector u 0 and strictly positive integer m, where T D mt, such that where O G mj and O H mj are .n C 1/ 2 coefficient matrices and O U m is a .n C 1/ .n C 1/ coefficient matrix. Also, u m D OEu.x k ; mt/ kD0:::;n . Since some p j and u j in (5) are unknown, Equation (4) is invoked, that is: with the notations H D : Combining linear systems (5) and (6) yields where I is a .n C 1/ .n C 1/ identity matrix. In the floating boundary method, the Signorini condition is considered as a switch between Dirichlet and Neumann conditions. As soon as the contact duration T c is set, the time step at which the switch happens is known, so does the boundary condition (Dirichlet or Neumann) at each time step at x D L. Thus Equation (7) can be solved by reorganizing the known and unknown entries in u and p, according to the boundary condition at x D L at each time step. Accordingly, by first assuming only one contact switch per period, such a strategy systematically skim through given intervals of values on period T and contact duration T c (for a chosen t). However, physically unacceptable solutions can be found. Accordingly, admissibility is systematically checked to ensure complementary conditions (2) are satisfied. Non-admissible solutions with non-admissible penetration or contact forces are discarded. Such way of Signorini boundary condition enforcement is similar in spirit with a precedent scheme [12], but also shows differences in detail.

Frequency Domain Boundary Element Method (FD-BEM) is a frequency-domain form of BEM which targets periodic solutions. Via a Fourier transform in time at frequency
the wave equation (1) becomes the well-known one-dimensional Helmholtz equation where Ä D !=c is the frequency number. Similar to TD-BEM, FD-BEM is based on a weighted residual form of (9), where the weight function is the fundamental solution to the Helmholtz equation [2] with the corresponding BIE where Q p is the Fourier Transform of p. Reading (10) on the boundary leads to Through (11), Q u at frequency ! can be solved for a combination of Dirichlet and/or Neumann boundary conditions. The time domain solution is then recovered through an inverse Fourier Transform. However, such a formulation is not well adapted when u.L; t / and p.L; t/ are both constrained by the complementary condition (2) expressed in time domain. This complementarity condition can be equivalently recast into the equality [9] r.t/ Á p.L; t/ C max.0; .u.L; t / g 0 / p.L; t// D 0; 8t (12) where 2 R C . In the remainder, the residual r.t/ will be made to vanish in an integral sense only. To insert (12) into FD-BEM with the aim of finding periodic solutions, Fourier expansions are considered, that is where ! 0 D 2 =T is the base frequency. The corresponding base frequency number is Ä 0 D ! 0 =c. Knowing that u.0; t/ D 0, it follows from (11) that tan nÄ 0 L nÄ 0 .a n cos n! 0 t C b n sin n! 0 t/ Since initial velocity is assumed to be zero [11], all sin terms vanishes. Only constant and cos terms remain in equation (13) and equation (14). The truncated Fourier series of u.L; t / and p.L; t/ with m harmonics only is then substituted into (12). The Harmonic Balance Method (HBM) [5] is used to solve for the unknown coefficients .a 0 ; a 1 ; a 2 ; : : : ; a m /. As a special case of Galerkin techniques, HBM reads where the following family of functions is selected: (1, cos ! 0 t, cos 2! 0 t, : : :, cos m! 0 t).
With m C 1 equations, m C 1 unknown coefficients can be solved for by the (Non-Smooth) Newton method. The function p.L; t/ is then reconstructed through (13) and so are other needed functions used to build the Nonsmooth Modes of Vibration. Basic continuation is implemented [7] in order to construct the desired backbone curves, by increasing ! 0 by a small value and solve again, with previous solution as initial guess.

RESULTS
In this section, the following is considered: L D 1, c D 1 and g 0 D 0:5. For TD-BEM, the bar is discretized into n D 50 elements in space with element length x D L=n and time step is set accordingly. frequency-dependence of the total energy, and compares well with existing solutions [11], except for backbone curve near linear mode. Since contact period length decreases as frequency ! decreases [11], short period of contact cannot be evaluated precisely with m D 15 harmonics. Due to inaccurate Fourier expansion and residual function, HBM struggles converging to useful periodic solutions. Thus this part of solution is skipped.
It can be observed that TD-BEM captures more than one possible solution for a given frequency of vibration, which shows similar pattern of nonsmooth modes distribution with existing results [10]. Two such solutions feature the same frequency but distinct contact durations, as shown in figure 3 (see points a and b). This is induced by the full internal resonance in the investigated system [12]. Although solution b shows a low-frequency pattern similar to solution a, extra energy is introduced due to higher frequency mode. Comparatively FD-BEM/HBM currently cannot guarantee convergence to internal resonance solutions. More advanced continuation technique might be required to find internal resonance branches.
Displacement solved by FD-BEM/HBM at point c on backbone curve is depicted in Figure 4. Overall, the displacement field is similar to solution at point a in Figure 3. It is worth noting that Signorini conditions in FD-BEM/HBM is enforced in integral sense through (15), unlike in TD-BEM the complementary condition is met pointwise, ie for every time step. Accordingly, residual pointwise penetration occurs as shown in figure 5.
In terms of computational efficiency, FD-BEM/HBM is much superior to TD-BEM/Shooting. With the suggested discretization settings, TD-BEM is about 700 times slower than FD-BEM/HBM to find the whole backbone curve. When n increases, such difference is expected to increase, since the total number of combination between T and T c to be skimmed increases in O.n 2 /, while for FD-BEM/HBM the step size of continuation can be arbitrarily controlled.

CONCLUSIONS
In this paper, two methods have been introduced to perform nonsmooth modal analysis on a simple academic system. Both methods have shown capability to converge to admissible solutions.
The TD-BEM/shooting methodology has shown good capability of capturing highly-detailed step-wise admissible solutions. However, its challenging implementation leads to high computational costs. Although it has shown no numerical energy dissipation in this one-dimensional case, this might be untrue for higher dimensional systems, in which case it would become ineligible for such modal investigations.
Comparatively, FD-BEM/HBM features a much lower computational cost, at the cost of solving Signorini conditions in an integral sense only. This formulation is also energy preserving by construction and could be extended to problems in two or three dimensions. More advanced continuation techniques, such as pseudo-arclength [7], are expected to be coupled into existing FD-BEM/HBM scheme.