Camille Schreck

Damien Rohmer

Doug L James

Stefanie Hahmann

Marie-Paule Cani

Real-time sound synthesis for paper material based on geometric analysis

Keywords: (according to ACM CCS): I, 3, 7 [Computer Graphics]: Three-Dimensional Graphics and Realism-Animation; H, 5, 5 [Information Systems]: Information Interfaces and Presentation-Sound and Music Computing

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Thanks to the advances in the Computer Graphics field, it is now possible to create visually compelling virtual worlds. Although very common in daily life, the characteristic appearance and sound produced by a crumpled sheet of paper, or the universally recognized friction of dollar bills, are missing in virtual environments. Contrary to rigid objects whose shapes do not vary in time, or to cloth garments whose sounds are not strongly dependent on shape, the sound of paper depends heavily on its shape and also changes dramatically during crumpling phenomena. Therefore, in order to achieve real-time sound generation for interactive virtual paper models, a highly efficient analysis of the shape of paper, coupled to fast sound synthesis, is therefore required.

Although very familiar to humans, paper sound is a complex combination of different sound styles. It is a mix between continuous noisy sounds produced by frictional sliding, and discrete events produced by geometric bending and crumpling processes. These discrete sounds may also vary between long "flap" sounds when the sheet is still smooth, and more short "clac" sounds in more crumpled cases. Real-time synthesis of this variety of sounds, with strong dependence on the paper's shape and its environment, is therefore a challenging problem.

Our sound synthesis approach consists of four steps performed at run-time during the visual simulation of the sheet of paper. First, sound-producing events are detected by analyzing the paper surface animation. Second, each event is characterized using the geometry and the motion of the surface to parameterize our sound models. Third, the sound is synthesized using the estimated model parameters using computed characteristics and two pre-recorded databases of primitive sounds: one for the discrete crumpling sounds, and another for the continuous friction sounds. Fourth, the resulting sound is spatialized in 3D to account for the influence of nearby planar surfaces, such as a wall or table.

Our three main contributions are as follows. First, we develop a new shape-dependent model for friction and crumpling sounds. Our model is based on local surface regions we call resonators, which characterize the region where significant sound-related vibrations occur. Second, we propose a two-sided detection and characterization method for both friction and crumpling sounds, that incorporates their dependence on nearby free (unconstrained) boundaries of the sheet, leading to four different sound types: constrained and unconstrained crumpling sounds depending on the presence of a free border in the producing region; and similarly, constrained and unconstrained friction sounds. Third, we propose a new efficient approximation method for sound spatialization, adapted to model the sound of thin-sheet surfaces near flat surfaces, such as tables or walls.

Related Work

The generation of sound has attracted an increasing interest in computer animation over the last few years. Perceptual studies and some synthesis methods are gathered in [START_REF] Giordano | Perception and synthesis of sound-generating materials[END_REF].

Data-driven synthesis: In the traditional approach, still widely used in in motion pictures and video games, real sound effects are recorded and edited to match the visual display by Foley artists. The approach produces high quality results but is time consuming and cannot be used for interactive environments. The process can be automated for interactive applications by triggering by computer-generated events and rendering pre-recorded sound samples [START_REF] Takala | Sound rendering[END_REF]. Unfortunately, this approach lacks of the fine-grained variability and synchronization needed to interactively render complex interactions occurring with paper manipulation. * 99] propose to use wavelet trees to model sound texture and synthesize audio backgrounds. The "Sound-by-Numbers" method [START_REF] Brooks | Sound-by-numbers: Motion-driven sound synthesis[END_REF] uses a granular synthesis model driven by a low-dimensional motion signal. Concatenative sound synthesis (CSS) [S * 00] methods consist in selecting in an large database of sound units, the units that match the best the sound to be synthesized.

An et al. [START_REF] James | Motion-driven concatenative synthesis of cloth sounds[END_REF] built on these methods to automatically synthesize the sound of a physically based cloth animation. First they analyze the deformation to find crumpling and friction events that drive the synthesis of a low-quality target signal. A CSS process is then used to select the best units from a database of recorded cloth sounds. Our method is inspired by the first part of this approach, and adapts it to paper material. In the case of paper, the sound is also produced by crumpling and frictional sliding mechanisms, but in addition the sound is highly dependent on the shape of the paper.

Physically based synthesis: Another approach for simulating sound is the physically based modeling sound techniques founded on the description of the physical phenomenon involved in the sound generation. In computer graphics, these techniques usually synchronize physically based sound synthesis to physically based animation. O'Brien et al. [START_REF] Cook | Synthesizing sounds from physically based motion[END_REF] uses an FEM (Finite Element Method) to synthesize both the animation and the sound. A wide range of complex scenes (including sheet-like objects) can be simulated but the high computation time required prevents it from being efficient enough for interactive application.

Modal vibration models [START_REF] Pentland | Good vibrations: Modal dynamics for graphics and animation[END_REF][START_REF] James | Dyrt: dynamic response textures for real time deformation simulation with graphics hardware[END_REF] have been introduced to digital sound generation by Adrien [START_REF] Cani | The missing link: Modal synthesis[END_REF]. The linear modal synthesis methods have become widely used for rigid-body sound [DP98, vdDKP01, Coo02, OSG02, JBP06, ZJ10]. Moreover it is possible to model a large number of sounding objects at interactive rates using acceleration techniques [RL06, BDT * 08], and also memory compression [START_REF] Langlois | Eigenmode compression for modal sound models[END_REF]. Ren et al. [START_REF] Ren | Example-guided physically based modal sound synthesis[END_REF] presents a method to estimate modal parameters, from recorded audio clips. The linear modal synthesis has been also considered for non-linear thin-shell sounds by Chadwick et al. [START_REF] Chadwick | Harmonic shells: A practical nonlinear sound model for near-rigid thin shells[END_REF], but such models only support small deformations.

Power-law models of the acoustic emissions from crumpling elastic sheets [START_REF] Kramer | Universal power law in the noise from a crumpled elastic sheet[END_REF] and paper [START_REF] Houle | Acoustic emission from crumpling paper[END_REF] have been studied by the physic community, and have led to geometry-independent stochastic sound models of crumpling [START_REF] Fontana | Physics-based sound synthesis and control: crushing, walking and running by crumpling sounds[END_REF]. Procedural stochastic methods have also been used to create a simple tearing sound model for paper [LFD * 15], but this model cannot capture the rich sounds produced when crumpling paper.

Overview

Two kinds of sounds are usually identified for thin-sheet material:

1. Friction sounds caused by regions of the surface sliding along either another object of the scene or another part of the same surface. The sound is produced by the sliding part bumping randomly against small irregularities of the other surface and can be approximated by a colored noise (see Figure 3).

2. Crumpling sounds generated in regions where the bending direction changes suddenly, such as when the sign of surface curvature changes. The transition between two equilibrium states can produce a sudden burst of energy which is quickly dampened (see Figure 4).

Paper material is extremely rigid compared to other thin materials, such as fabric. Consequently, resonant vibrations can have a strong impact on the sound. The size and the shape of theses resonating regions influence the resulting sound. A good experiment to convince oneself of the importance of this phenomenon is to listen to the difference between the sounds produced while sliding a flat piece of paper over a table, compared to sliding a curved piece of paper enabling the sound to resonate (Figure 3). We call free borders natural edges of the paper that are curved enough to be able to vibrate (see Figure 5 (left)). As they have more degrees of freedom than interior edges or borders constrained to be flat, they have an notable influence on the sound. When the resonating region has no free border, we observe that the friction sound is close to a white noise (Figure 3 left) and is less loud than friction sounds produced by a bent surface containing therefore a free border (Figure 3 right). Similarly for the crumpling sounds, we register sharp "clac" sounds (Figure 4 left) -which you can hear while crumpling a paper into a ball -that are produced by constrained regions dampening the sound in a few milliseconds. On the opposite, when the vibrations can reach a free border a longer "flap" sound lasting few hundreds of milliseconds (Figure 4 right) -which you can hear while wiggling a sheet of paper held on one border -is produced. Based on these observations, we therefore make a difference between sounds produced by regions containing one or more free border(s), that we call unconstrained regions, and those produced by regions that contains none, i.e. constrained regions (see Figure 5 (right)).

Our key idea is to reproduce these resonance phenomena thanks to a procedural detection of representative regions in which the vibrations caused by a crumpling or friction event occur. This enables us produce, for each of these regions we call resonators, an appropriate shape-dependent sound. Figure 5: Free borders are natural edges of the paper that are curved enough to allow enough degrees of freedom to be able to oscillate (in red). A region of the surface is considered as constrained if it has no free border, and unconstrained if it contains at least one.

In order to output sounds as rich as those of real paper in real time, we use pre-recorded sound units stored in Crumpling and Friction sound databases. Due to their different nature, we need the whole sound to replay a crumpling sound, while we just use the time-averaged magnitude spectrum to reproduce a friction sound (see §5). Therefore our Friction sound database only stores spectra.

Our processing pipeline is summarized in Figure 2. Given an animated mesh modeling the deformation of the paper surface through time, we extract two types of sound source events at each time step:

(1) we compute the parts of the surface in contact with other objects, since each of these parts defines a friction event; and (2) we detect surface regions where the curvature changes sign, and use these features to identify crumpling events as in [START_REF] James | Motion-driven concatenative synthesis of cloth sounds[END_REF]. For each sound source event, we then compute a set of resonators, i.e., the regions in which vibrations caused by this event occur (see §4). Resonators are then parameterized by their shape and size (see §5) and subdivided into two classes: unconstrained resonators or constrained resonators, depending on whether or not they have a free border that is able to oscillate. For each resonator, a sound is ex- tracted from the appropriate database depending on the type of resonator and the type of sound event. This leads to the four categories of sounds described above: clac and flap sounds for respectively constrained and unconstrained resonators in case of crumpling events, and different friction sounds for constrained and unconstrained resonators in case of friction events. A specific sound is extracted based on the geometric parameters of the resonator, and then reduced or amplified depending on the magnitude of the sound-source event. Finally, the assembled sound is spatially embedded within a 3D environment, by taking into account the relative positioning of the paper with respect to surrounding obstacles such as tables or walls, and the position of the listener (see §6).

Detecting resonators

Geometric Model: At each animation step, the first stage of our method consists in detecting the sound-source events on a 3D virtual paper surface provided by an external animation or simulation engine. As we aim at providing sound synthesis in real-time, choosing an interactive animation system for paper crumpling is mandatory. In this work, we use the model proposed by Schreck et al. [SRH * 15], which is the only existing model capable of animating paper crumpling at interactive rates. It has the advantage of being built on high level geometric descriptors, which we can reuse. In particular, the approach combines coarse physical simulation with developable geometry processing steps. The geometry generation part is based on the hypothesis that the creased curves, which appear where microscopic fiber structure has been broken causing the surface not to be C 2 anymore, can be approximated by singular points and thus the developable surface of the paper can be approximated by a set of generalized cones. Although our system could work on top of other thin-sheet animation methods, this approach has the advantage that it directly provides a partition Rpaper of the sheet into developable regions defined by singular points (see Figure 6). We can distinguish between singular points respectively located on the border (yellow dots), inside the surface where fibers are broken (red dots), and anchored points (black dots) where the user's fingers guide the paper deformation. The partition Rpaper is composed of:

1. Curved regions (in purple in Figure 6) that segment sections of generalized cones, and are represented by the rulings of these cones (purple lines in Figure 6). Each ruling is defined by two points, with at least one of them belonging to the border of the paper.

2. Flat regions (in green in Figure 6) that are represented by a set of coarse triangles whose vertices are either singular points or anchored points. The dihedral angle between every pair of adjacent triangles of a same flat region is inferior to a threshold, εa.

We use this surface partition to help represent surface resonators by considering each region of Rpaper as a possible resonator. As curved regions have necessarily a curved border they are considered in the following as unconstrained regions, whereas flat regions are considered as constrained regions (see Figure 6 (left)).

Event Detection: We detect the sound producing events using a heuristic similar to the one proposed in An et al [START_REF] James | Motion-driven concatenative synthesis of cloth sounds[END_REF]. Friction events are associated with vertices in contact with an obstacle (based on proximity computations), that are spatially connected. Each connected component represents the locus of a friction event.

The mean speed of its vertices represent the amplitude of the event.

The location and time of crumpling events are estimated by analyzing the mean curvature over the surface; crumpling events are represented by connected components of vertices whose mean curvature has undergone the same change of sign between the current and the previous animation step. The associated amplitude is computed as the sum of the curvature changes for all the vertices contributing to the event [START_REF] James | Motion-driven concatenative synthesis of cloth sounds[END_REF].

As already stressed, contrary to cloth material where crumpling and vibrations only occur locally within the material, paper is very stiff (because of strong isometry constraints), and can support spatially coherent vibrations and buckling events over large regions. The sound produced by friction and crumpling events is thus much more dependent on the current geometrical state of the surface than it is for softer materials, such as cloth. Therefore, instead of directly associating a sound to each event as in [START_REF] James | Motion-driven concatenative synthesis of cloth sounds[END_REF], we first associate a set of geometrically defined resonators to each event, and then synthesize the sound accordingly. While this is a crude approximation of the complex underlying phenomena, it enables the sound to have geometric dependence which can improve realism.

Resonator Estimation: In order to compute the set of resonators for an event, defined as a collection of connected surface pieces, we need to identify the boundary of the surface region that the event causes to vibrate. We choose to define this boundary in an intuitive way as follows. The borders of the sheet naturally belong to this boundary. Moreover, paper material tends to get very stiff in the direction of highly curved folds, or "ridges" (which are aligned to some of the rulings; see Figure 6), but the sheet is more flexible in the main curvature direction, Consequently, transverse vibrations tend to be trapped in these interior resonator regions by ridges [START_REF] Gopinathan | Trapping of vibrational energy in crumpled sheets[END_REF]. Therefore, we assume that vibrations do not propagate across ridges, which we define as a ruling whose dihedral angle is above a threshold, amax; we use amax = 0.5rad in our implementation. This is done by stamping as boundaries of the propagation region all the edges with high dihedral angle, in addition to free borders and to edges with an anchored point (position of virtual fingers manipulating the sheet) or in contact with another object, since this also prevents further propagation.

Defining these limits enables us to estimate the set of resonators where a sound event propagates as follows (see Figure 7): We initialize the set with the region(s) of Rpaper where the event took place, and then recursively collect the neighboring regions until reaching an edge stamped as propagation border. Note that the edge that stops propagation may be located within a conical region of the shape. In this case, we only add the relevant sub-part of this region to the set of resonators. Each of the regions of Rpaper collected is considered as a resonator and will produce a sound as explained in the next section.

Geometry-based sound synthesis

Resonator Parameters: As already stated, our method to produce rich and natural sounds in real time is to re-play pre-recorded sound units, selected from a database. We therefore need to parameterize resonators in order to match them to a specific sound unit, easy to query in a database, and we need to register sounds in the database accordingly. The number of parameters to be kept small in order to keep the database to a reasonable size. They should also discriminate well enough the different possible sounds.

In this work, we propose to parameterize each resonator by the following two variables: lc the sum of lengths of the all free curved borders (red in Figure 8), and lr the mean length of the rulings (blue in Figure 8). As already discussed, free borders of the surface are of major importance to characterize the sound. In selecting their length, we aim at extracting the parameters of the dominant sound. By associating it with the mean length of the rulings, we sample uncorrelated variables related to the area of the resonator. We therefore claim that the couple (lr, lc) is a good choice of parameters to correlate significantly a geometrical measurement of a region to the sound it produces. Note that we further validated this choice by performing real measurements, as described in §7.1.

Sound Databases: Let us detail how we pre-record the databases of sound units to be played at run time. Firstly, such databases are specific to a given type of paper material and to a given type of obstacle for the friction sound. In our experiments we used printer paper (80g/m 2 ) and a wooden table. We also recorded databases for tracing paper and paper of a bank note.

We actually created two different databases, respectively for fric-tion sounds and for crumpling sounds. As friction sounds can easily be continuously synthesized as noise, without taking into account their temporal phase, we only store the magintude spectrum in the frequency domain of the sounds recorded in the sound database. In reverse, crumpling sounds have a specific time duration with a specific beginning and end, and are therefore stored in the temporal domain.

For these two databases, we sampled the (lr, lc) space, and for each value, we cut a rectangular sheet of paper with (lr, lc/2) edge length, bent its edge of length lc/2, recorded the sound produced by its friction, and by a change of sign of curvature (for the crumpling database). The unconstrained friction sound and the flap sounds are related to non-zero free border lengths and are the most dependent on the shape of the surface. The two other constrained types of sound, namely constrained friction and clac do not depend much on the shape of the surface and were therefore considered as special cases. To ease the recording process, we stored these two specific sounds in the (lr = 0, lc = 0) entry of their respective database. Several sound samples were stored per entry to increase variety.

Sound Synthesis: Given the databases and the resonator, sounds are synthesized at run-time in the following way: For each resonator, we select the closest recorded (lr, lc) parameter in the relevant database, considering also the special case (0, 0) when the resonator does not have any free border, and randomly select one of the corresponding sound samples. Finally, all sounds are assembled together. Friction sounds are synthesized using the inverse Fast Fourier Transform method [RDRD92, MAKMV10]: friction spectra are summed up in the frequency domain. We add random phase to the resulting magnitude spectrum by multiplying each value of the spectrum vector by the corresponding value of a vector of complex numbers {c i } = {e j ωi } where ω i is a random angle chosen uniformely in the range [-π, π]. A buffer of the global friction sound is then obtain by applying an iFFT to the resulting complex vector. In order to smooth the transition between them, two consecutive buffers are overlapped. Moreover the amplitude of the sound is linearly scaled with respect to the mean speed of the associated relative frictional motion. The crumpling sounds are played until the end of the recorded unit. Their amplitude are also scaled linearly with respect to the magnitude of the change of curvature related to the event.

Spatialization

In this section, we describe how we embed the synthesized sound in the 3D space, by taking into account both the listener's position and sound reflection on the surrounding table or walls. Such planar surfaces interact with the sound vibrations of paper material which can be modeled as a dipole.

The dipole-like radiation behavior of paper surfaces leads to significant frequency-related interference effects above planar surfaces. For example when sliding a sheet of paper on a table to bend it more, as in Figure 1 (left), one may notice a pitch-shift. This is partly caused by the fact that the friction sound is reflecting on the table and the reflected source interferes with the original source. A pitch-shift is also noticeable when one moves a sheet of paper closer to a wall while wiggling it. As the paper gets closer to the wall, the "flap" sounds change. Note that in this work, we do not consider the interference caused by self-reflections on the surface by cavities formed by the paper.

We develop below a simplified model for these phenomena for both time-domain (crumpling) and frequency-domain (friction) sounds. Transverse vibrations of the paper sheet give rise to positive and negative pressure fluctuations on either side of the sheet. If we approximate each resonator a small vibrating disk sound sources, then we can reason about each disk's sound emission and any planar reflections.

Let us consider the sound field due to a vibrating rigid disk of radius a and zero thickness. If we assume that the normal velocity is U(t), and that the disk is acoustically compact, i.e., the important wavelengths associated with the frequencies in U(t) are much larger than a, then one can derive the acoustic pressure field p generated by the motion (see [START_REF] Howe | Theory of Vortex Sound[END_REF] page 81, problem 9) as: where the listener is located at x, the center of the disk is located at y, n is the unit normal vector, and c is the speed of sound in the air.

p(R,t) ≈ 2ρ 3πc a 3 cosθ R ∂ 2 U ∂t 2 (t - R c ), (1)
The relative location is given by:

R = x -y , R = ||R|| , R = R/R , cosθ = R • n (2)
as shown on Figure 10. This vibrating disk produces a dipole pressure field, with characteristic 1/R distance scaling, and a cosine lobe about direction n. Without loss of generality, we can discard the scalar factor 2ρ 3πc (which is common to all our sound sources) and use:

p(R,t) ≈ a 3 cosθ R A(t - R c ), (3)
where A = ∂ 2 U ∂t 2 is the acceleration of the patch along n. Consider the dipole source from (3) at height d above an infinite rigid plane, as shown in Figure 10. Since the rigid plane imposes a zero boundary condition ∂p ∂n = 0 on the normal acoustic particle velocity on the surface, a fictitious image dipole source p is introduced in the sound field. The global pressure field P is then the sum of the pressure fields of the dipole source and its image source: P(x,t) = p(x,t) + p(x,t).

(4)

The image source has reflected position y and normal n, and corresponding relative coordinates R and cos θ = R • n. Using this simplified model, we approximate sound reflections for crumpling and friction sounds using time-and frequency-domain reflection models, respectively, as follows.

Time-domain crumpling sounds: we scale each current timedomain crumpling sound pressure event A(t) selected from the database to obtain its dipole-spatialized version and generate the following delayed sound source for the image dipole:

P(x,t) = a 3 cosθ R A(t) + a 3 cos θ R A(t - R -R c ). ( 5
)
This change is done to each crumpling sound before synthesis.

Frequency-domain friction sounds: we obtain the dipolespatialized version of each friction sound by placing Equation (5) in the frequency-domain (effectively replacing A(t) by A(ω)e -iωt ):

P(x, ω) = a 3 cosθ R e iω R c + a 3 cos θ R e iω R c A(ω). ( 6
)
This multiplication factor is applied to the friction spectra A(ω) prior to region-based summation and synthesis of the final sound.

Results

Please see the accompanying video for the results described in this section, and listen to the video using good headphones in a quiet environment.

Implementation Details: In our results, we use sounds recorded at a sampling rate of 44100 Hz and play the generated sounds at the same rate. The spectra of the friction sounds contains 2048 frequencies which means that a buffer contains 2048 samples (≈ 50 milliseconds) and we use an overlap of 1024 samples between two buffers. Our implementation makes use of the C++ library STK [CS99, SC05] to read and play sounds.

Recording of the database: Our two databases contain samples of sizes lc and lr, using the values {2.5, 5, 7.5, 10, 12.5, 15, 20, 25, 30} cm. We recorded three crumpling sounds and one friction sound for each sheet size. The sounds were recorded manually under the same conditions. The microphone was placed one meter away from the sound source. For the friction sounds, the speed was measured by the time needed to go over the friction distance (2s for 20cm). The shape was retained by keeping the same amount of compression of the curved edges (30%). We recorded the unconstrained ("flap") sounds by holding both ends of a sheet (slightly curved by compression curved edges of 10 %) of the required size and twisting them to change the sign of the curvature. The constraint ("clac") sounds were isolated in a recording of crumpling paper.

Validation of the resonators' parameters

We used two parameters (lc, lr) to characterize the sound of a given resonator. To validate this choice in the case of the friction sound, we show that modifying these two parameters has a larger influence on the sound produced by friction than modifying other geometrical shape features with constant (lc, lr).

We performed the following experiment: let us consider a reference square paper sheet of size 10cm × 10cm. On one hand, we prepared rectangular paper sheets with smaller or longer edge length. Each rectangular sheet was then bent in the configuration shown in Figure 11 ((b) and (c)) such that their edge length respectively corresponds to lc/2 and lr. On the other hand, we used the reference paper sheet -fixing therefore (lc, lr) -and modified only its shape by compressing the two extremities of the sheet as illustrated in Figure 11 (a). In all these scenarios we recorded the sound produced by the friction of the sheet of paper on the same table, at same speed. The relative variation of the parameter is in the range of 50 to 75%.

To compare the different recordings, we computed the histogram corresponding to the magnitude spectrum of each sound and then compared them using the quadratic-form distance [START_REF] Pele | The quadratic-chi histogram distance family[END_REF]. Table 1 gathers our results. The first four columns show that modifying the lc or lr parameters respectively by a factor of 50 to 75% leads to a modification of their spectrum, quantified between 0.07 and 0.16. The next two columns show that keeping constant lc and lr parameters but changing the curvature of the sheet leads to a change of spectrum quantified between 0.03 to 0.05, which stays always below the influence of tuning lc or lr.

To further study the influence of shape, we also compared the sound produced by three different sheets of shape S 1 , S 2 and S 3 , shown in Figure 11 (d). Each piece of paper is compressed by a factor of 30% and the results on their respective histogram distances computed using the same approach as above are gathered in the last two columns of Table 1. We obtained distances comparable to those obtained for change of curvature, and once again lower than those obtained for changes of lr and lc. Table 1: Distances between friction spectra of rectangular sheets of paper while varying the parameters as shown in Figure 11. lr and lc/2 are the considered edge length parameters in cm, dr and dc are the respective histograms distances. comp is the relative compression ratio, dcomp the associated histogram distance. ds are the distances between friction spectra of the shapes s described in Figure 11(d) compressed of 30%. s refers to the shape pattern show in Figure 11 (right), and ds corresponds to the histograms distances. For example, the first line of the first table indicates that the histogram distance is 0.10 between the friction magnitude spectrum produced by a resonator with lr = 5 cm and the one produced by a resonator with lr = 10 cm.

Comparison with real paper material

To validate our resonator model, we compared the results obtained using our method with sounds produced by real sheets of paper of different sizes. First, we focused on friction sounds produced by sheets of paper slid over a table. We used rectangular sheets of different sizes -210x297 millimeters (A4 format), 210x148 (half A4),105x297 (A5), 105x148 (A6)-and folded them as shown in Figure 12 (left). We also recorded the sound of an A4-sized flat sheet. We compared the resulting five spectra of virtual and real sounds in log scaled space (respectively top-right and bottom-right of Figure 12). The results show that our model follows a similar evolution of the magnitude spectrum with respect to the sheet size, i.e. a general translation of the spectrum towards higher frequencies when the size decreases. This evolution -and particularly the clear difference between flat paper and folded paper being slid along the same table-can also be heard in the video. Note that our synthesized sound is still less rich in terms of frequencies variations than the real one, in particular we do not reproduce phenomena such as the stick-slip sound that can be heard when real paper begins to slide. Some more advanced interactions can also have an influence -see for instance the case of musical instruments like violon [START_REF] Pitteroff R | Mechanics of the contact area between a violin bow and a string. part i: Reflection and transmission behaviour[END_REF].

Next, we compared real and virtual crumpling flap sounds for different sizes of sheets. Considering the same sizes than in the previous experiment, we folded each sheet of paper as shown in Figure 13 and revert the sheet while recording its sound. The attached video plays the resulting sounds. In both (real and virtual) cases, we get distinctive sounds depending on the size of the sheet.

Spatialization

Here, we conducted two experiments. Figure 14 (left andmiddle) and the associated video present the result of our frequency-domain spatialization method for the friction sound. The friction sound of our virtual paper is recorded for an A4 sheet compressed by a factor of 30% first and 10% then, leading to a change of height of the sound source. One can hear in the video that the friction sound depends on the height of the sound source, as it is the case for real paper sound. We can also hear the slight pitch shift as the sheet is bent. We also show in the video the modification of sound depending on the position of the listener in front and on the side of the sheet of paper (Figure 14). Figure 15 corresponds to a sheet of paper being wiggled while changing its distance with respect to the wall. The time-domain spatialization of the crumpling sounds causes the sound shift we can hear in the video.

Computational Performance

Memory requirements: The crumpling database contains 249 sound units whose lengths are going from 200 to 800 milliseconds. The friction database contains 63 spectra sampling 2048 frequencies -between 0 and 22kHz. The memory footprint of the database is around 10MB.

Computation times:

The detection of the events and the computation of their resonators are done once per frame. Table 2 compares the mean time per frame for these audio operations to the computation time of the animation. We computed the sound of the animation represented Figure 18 (right) with three different resolutions. The computation time of the curvature depends on the resolution, so real-time performances get harder to reach when the resolution increases. Still, the time required for the audio-related operations always remains negligible compared to the time required for the animation. Playing sounds is done using different threads than computing the animation To get real time, it is necessary to generate at least 44100 audio samples per second. A thread is created for each crumpling sound. We only applied scaling operations to the recorded samples before sending them to the audio output. A single thread is used for In this section, we present some more complex examples, which are also present in the attached video. Figures 1 and 16 show, respectively in the virtual case and in the real case, the spectrogram of the sound generated by while sliding one edge of a sheet along the table thereby curving the paper, then pinching the front border causing the back to stand up. This example mixes all the different categories of sounds we considered. We note that the friction sounds are much more regular than the synthetic ones. This is mainly due to the fact the it very difficult to reproduce a smooth movement "by hand." Nonetheless, we obtain the same type of event at the same moment, notably the "flap" sound when the sheet stands up and different kinds of friction sound for the different shapes that the sheet takes. Figure 17 shows a sheet (real and virtual) held by one border and wiggled and the corresponding spectrograms. Figure 18 (left)

shows the example of a sheet being crumpled. We also simulated the sound of different categories of paper material by recording different databases, such as for tracing paper and paper of a bank note (Figure 18 (middle) and (right)).

Discussion and conclusion

We presented the first method that automatically generates plausible sounds fitting an interactive animation of paper material. We achieved real-time sound synthesis through a trade-off between speed and accuracy, leading to several noticeable limitations. First, the range of possible sounds is limited when using pre-recorded sound databases. Moreover, each database is associated with a specific type of paper and obstacle materials (for friction sounds), and changing them requires recording a new database.

Another limitation is due to the small set of physical parameters taken into account during the synthesis process. By requiring only animated mesh sequences, our method can handle a large variety of animation inputs. However, since it does not have access to physical contact force information, it can not, for instance, include dependence of friction sound volume on contact force. Moreover, our resonator model approximates paper vibrations in a resolutiondependent manner, and does not actually model vibrations directly. Challenging but familiar scenarios, such as crumpling paper into a ball, remain outside the scope of this work, since they would require detailed simulation of self contact, as well as more sophisticated modeling of sound inter-reflections and scattering.

The visual simulation of paper is also a bottleneck, and currently there are no methods that are both interactive and able to model very crumpled paper, especially with self contact. We used an interactive method to demonstrate our fast sound synthesis, and also because it enabled efficient estimation of resonators. Also although we did not obtain very crumpled paper, our results depict a large range of scenarios.

Finally, in future work, we would also like to take into account the specific properties of developable surfaces to allow for more accurate vibration modeling without sacrificing computational efficiency.

  Other data-driven methods are inspired by granular synthesis techniques [Roa04]. Dubnov et al. [DBJEY * 02] and Bar-Joseph et al. [BJLW

Figure 2 :

 2 Figure2: Overview: After computing the sound-producing events, each of them is associated with a set of resonators. Each of those resonators leads to a sound unit after querying it on a pre-computed sound data-base. Each sound unit is spatialized to adapt to our scene and the final resulting sound is obtained by mixing all the contributing sounds.

Figure 3 :

 3 Figure 3: Spectrogram of recorded friction sound reveals slowly varying spectra in time, which is well approximated by a colored noise model. (left) Constrained friction sound obtained for a flat sheet of paper. (right) Unconstrained friction sound produced by a curved sheet.

Figure 4 :

 4 Figure 4: Spectrogram of recorded crumpling sound: (left) "clac" sound produced by a constrained resonating region. (right) "flap" sound produced by an unconstrained region.

Figure 6 :

 6 Figure 6: Geometric structure of surface animation model: The surface is approximated by a set of generalized cones (purple) and flat region (green). There are two kinds of singular points -interior singular points (red), where the fibers are broken and which are apices on cones, and border singular points (yellow) delimiting curved regions. Image courtesy of [Schreck et al. 2015].

Figure 7 :

 7 Figure 7: Resonator estimation (for a crumpling event): (left) curved and flat regions, (middle-left) detecting sound event, (middle-right) finding region of event and recursively finding the regions in which the vibrations can propagate, (right) resonators.

Figure 8 :

 8 Figure 8: Resonator parameterization using the mean length of its (blue) rulings (lr) and the sum of the lengths of the (red) free borders (lc).

Figure 9 :

 9 Figure 9: Sound Databases: Ordered in a matrix form according to the dimensions lr and lc. The crumpling database contains short sound clips, whereas the friction database contains sound spectra.

Figure 10 :

 10 Figure 10: Scattering: Geometric configuration of a vibrating paper patch above a planar surface, along with its image source.

Figure 11 :

 11 Figure 11: Validation experiment: We compare sounds obtained while varying three parameters: amount of compression (a), lr (b) and lc (c). We also measured the sounds obtained for 3 sheet shapes with the same compression, lr and lc: S1 is a square 10 cm x 10cm (d top); S2 is a trapezoid whose bases, of respective sizes 4cm and 16cm, are the curved borders (d middle); S3 is a trapezoid whose bases of sizes 4cm and 16cm, are the flat borders (d bottom).

Figure 12 :

 12 Figure 12: Real-virtual comparison of friction sound spectra for different sheet sizes.

Figure 13 :

 13 Figure 13: Real-virtual comparison of "flap" sounds for different sheet sizes, obtained by changing the curvature of bent paper.

Figure 14 :

 14 Figure 14: Spatialization: The video shows the result of spatialization for friction sound from two different points of view.

Figure 15 :

 15 Figure15: Spatialization: The sound may be modified by the proximity of a wall. It is noticeable for example when one wriggles a sheet of paper near a wall. A change a pitch can be heard while moving the sheet away from the wall. The accompanying video shows similar effect obtains thanks to the spatialization step.

  (a) (Left) Virtual paper. (Right) Real paper. (b) Spectrogram obtained for (Top) virtual paper, and (Bottom) real paper. Both examples show 1.5s of the sound

Figure 17 :

 17 Figure 17: An A4 sheet of paper held by two corners is lightly shaken.

Figure 18 :

 18 Figure 18: Different paper materials being crumpled can be synthesized in recording different databases: (Left) printer paper, (Middle) tracing paper, (Right) a bank note.

Acknowledgements: This research was partially funded by the ERC advanced grant no. 291184 EXPRESSIVE.


L Zuo

A Curnier

NON-LINEAR REAL AND COMPLEX MODES OF CONEWISE LINEAR SYSTEMS

Conewise (piecewise) linear, non-gyroscopic and gyroscopic autonomous systems are studied in detail First, the general properties of such systems are examined. Then the non-linear real and complex modes are obtained for both non-gyroscopic and gyroscopic cases, by means of analytical or numerical tools. Several interesting characteristics of non-linear modes are found and compared to those of linear modes. Accurate analytical estimates of non-linear mode periods are formulated for both non-gyroscopic and gyroscopic systems. The validity of such estimates is discussed in relation to a non-linear mode bifurcation phenomenon, which is located by using both characteristic multipliers and Poincare mapping.

I. INTRODUCTION

Actual engineering systems often involve components which are in unilateral contact, operate with clearances, are made of bimodular elastic materials, etc. Such systems are governed by equations of motion containing piecewise linear terms. The example which has motivated this study is that of rotating shafts with a transverse crack. Indeed, when a cracked shaft revolves, the crack alternatively opens and closes and the shaft stiffness can be roughly modelled by a half-spacewise constant function [1][START_REF] Gasch | Institution of Mechanical Engineers Conference[END_REF][START_REF] Grabowski | [END_REF][4][START_REF] Zuo | 6tection des fissures dans les arbres de machines tournantes[END_REF]. Natural frequencies of rotating shafts are important in system design in order to avoid contingent resonances [START_REF] Tondl | Some Problems of Rotor Dynamics[END_REF]. Crack detection is even more crucial for obvious safety reasons. The reduction of the "natural frequencies" with the crack depth could serve as an indicator of crack initiation and propagation [1,[START_REF] Mayes | Institution of Mechanical Engineering Conference[END_REF][START_REF] Ismail | [END_REF].

Owing to its practical and theoretical importance, the dynamic behavior of piecewise linear (PL) systems has been the subject of many investigations. The free responses of single-degree-of-freedom (1-DOF) PL oscillators were mentioned by Timoshenko and Young in their book on vibrations [START_REF] Timoshenko | Vibration Problems in Engineering[END_REF]. In subsequent studies, most authors considered periodical excitations and used classical analytical tools for their investigations. When dealing with a beam in intermittent contact with a spring, Masri et a/. [START_REF] Masri | [END_REF] sought an analytical solution separately in the contact and gap regions and then matched these solutions through a condition of continuity at the switching point. A similar approach was followed by Natsiavas [I I, 12], who examined some particular oscillators with bilinear and trilinear restoring forces, and studied the existence and stability of their motions. For other PL models, Maezawa eta/. [13][14][15] used Fourier series for determining the superharmonic resonances up to second order; Comparin and Singh [16) employed the harmonic balance method to study the main resonance; Choi and Noah [17] combined the fast Fourier transform (FFT) with a Newton-Raphson scheme for approximating the arbitrarily forced periodic vibrations. Recently, Kim and Noah [18] have modified this method in order to enhance its reliability and efficiency.

During the past ten years, further studies on the dynamics of forced 1-DOF PL systems have been carried out with the help of more elaborate non-linear analysis tools: Poincare section, averaging method, local bifurcation theory, center manifolds, normal forms, etc. [START_REF] Guckenheimer | Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields[END_REF][START_REF] Lichetenberg | Regular and Stochastic Motion[END_REF]. It has been found that, in certain circumstances, the periodic responses of such systems under excitations may go through a bifurcation and exhibit complicated and unpredictable behavior known as chaos. A detailed analysis of a periodically forced PL oscillator was presented by Shaw and Holmes [21] in 1983, while extensive experimental and numerical results on the characteristics of such systems, in particular the period-doubling bifurcation, were supplied by Thompson et a/. [START_REF] Thompson | [END_REF][23][START_REF] Thompson | Nonlinear Dynamics and Chaos[END_REF][START_REF] Thompson | [END_REF] and Mahfouz and Badrakhan [26]. Similar observations in either mechanical or electrical systems were also reported by Schulman [27], Natsiavas [II], Chua and Hasler [28], Chua and Matsumoto [29,30], and Kim and Noah [18]. In addition, some methods for eliminating chaos in PL systems were proposed by Elvey [31], Thompson and Elvey [32] and Shaw [START_REF] Shaw | The Supression of chaos in periodically forced oscillators[END_REF].

Up to now, few results on the dynamics of multi-degree-of-freedom (MDOF) PL systems have been published. Recently, the frequency response characteristics of a MDOF system with clearances were investigated by Comparin and Singh [34] with the help of the harmonic balance method. Although their study provides a general method for the development of simplified approximate solutions, it is limited to harmonically forced systems.

In this work, we study autonomous conewise linear multi-degree-of-freedom gyroscopic and non-gyroscopic systems which, to our knowledge, have not been treated before. Such systems are of both theoretical and practical importance because they are the simplest models of cracked beams (non-gyroscopic) and cracked rotating shafts (gyroscopic) as well as of many other mechanical systems subject to unilateral contacts [1][START_REF] Gasch | Institution of Mechanical Engineers Conference[END_REF][START_REF] Grabowski | [END_REF][4][START_REF] Zuo | 6tection des fissures dans les arbres de machines tournantes[END_REF].

Our approach consists of combining analytical methods with numerical ones for investigating the fundamental dynamic characteristics of such systems; namely, their non-linear real and complex modes. In the light of the analytical results obtained for a 1-DOF model and a particular 2-DOF model, the concept of non-linear normal mode introduced by Rosenberg for smooth and symmetric non-linear systems [35][36][37][START_REF] Rosenberg | Advances in Applied Mechanics[END_REF][START_REF] Whin | [END_REF] is extended to conewise linear systems. An efficient numerical algorithm is developed to extract all the non-linear modes of gyroscopic and non-gyroscopic systems and compare them with linear modes. Then, an analytical estimate of the natural frequencies of conewise linear systems is established in terms of the system parameters. This estimate not only provides additional insight into system behavior, but also a check for the numerical results obtained. Finally, the stability of the non-linear modes is discussed by using modern analysis tools: characteristics multipliers and Poincare mapping. The bifurcation of 2-DOF non-linear modes is revealed.

CONEWISE LINEAR SYSTEMS

DEFINITION

The dynamic equation of the n-DOF piecewise linear gyroscopic systems considered here can be stated, together with its initial conditions, as follows. Find x(t):

[0, T]-+ ~' such that Mx + Gx + f(x) = o, x(O) = X 0 , x(O) = v 0 .

(I)

In this equation, M = Mr is a symmetric posiltve definite matrix representing the system mass, G = -Gr is a skew matrix describing the gyroscopic effect and f(x) is a half spacewise linear elastic force: i.e., a force defined, by means of two constant symmetric matrices S+ = s: and S_ = S': over two complementary half-spaces H;l' :, and H;l"_ respectively, as f(x)=Sx)x={S+x, VxEH;l':,={xlxrn>O}}•

( 2 ) ( S_x, VxEH;l" ={xlxrn~O} Thus, the two regions of linearity are half-spaces separated by a hyperplane xrn = 0, oriented by its unit normal vector n pointing inward R':, . Strictly speaking, models defined by equations (I) and (2) should hence be called half-spacewise linear systems. In such systems, the non-linearity is concentrated at the interface which includes the origin. It occurs when two system components are in unilateral contact in the rest position and one of the components has a negligible mass in comparison to the other. Systems involving several such unilateral contacts can give rise to more complicated subdivisions of H;l" into polyhedral convex cones centered at the origin [START_REF] Rockafellar | Convex Analysis[END_REF], by means of several hyperplanes. These systems would be more appropriately called conewise linear than piecewise linear. All the results presented in this paper on half-spacewise linear systems for didactic reasons can be directly extended to conewise linear systems.

A simple 2-DOF "gyroscopic" example is shown in Figure I. It consists of a mass point attached to a rotating table by two springs and in unilateral contact with a third (massless) one. For this example, equations (I) and ( 2) are more explicitly defined by M=mi,

G=( 2 ~Q -2 0 mQ)- n=(~)- S(x) = K(x)-m!J 2 1 2 , -{K+ =K:, xEH;l 2 + ={xlx1 >0}} K(x) - r ' -{ 0} . K_ = K_, X E H;l_ = X I x, ~ (3)
In equations (3), K+ and K_ represent the stiffnesses of the springs with and without contact, respectively, and 1 2 is the two-dimensional identity matrix. It should be emphasized that the symmetric matrix S(x) is not necessarily positive definite, due to the gyroscopic term -m!J 2

(5)

and 0 ~] E ~'"', 0 2.2. PROPERTIES
The main properties of MDOF half-spacewise linear systems are now deduced from their definition, as a preliminary analysis. All of these properties remain valid for conewise linear systems with minor adjustments.

Continuity off(x): existence and uniqueness of the solution x(t)

For mechanical and mathematical consistency, the internal force f(x) is assumed to be continuous across the interface. The following characterization is crucial.

Proposition. The half-spacewise linear force f(x) defined in equation ( 2) is continuous across the interface xTn ~ 0 if and only if the half-spacewise constant stiffness S(x) satisfies the jump condition

S+ = S_ + ODDT, ( 6
)
where o is a real coefficient. As a consequence, the solution x(t) of equations (l) exists and is unique.

Proof Continuity across the hyperplane~".. n W. requires that Vtjt'n = 0, S+ t = S_ t. It follows that the difference matrix s+ -s_ must be singular and in fact of rank one. It can thUS be expreSSed aS the dyadic prOdUCt Of a vectorb by the normal Vector D: S+-S_ =bOT. Owing to the symmetry of S+ and S_, it is further deduced that b =on. Conversely, post-multiplying equation (6) by t yields (S+-S_ )t =buTt= 0, which implies continuity. Now, it can be shown that the elastic force f(x) is Lipschitz continuous [START_REF] Zuo | 6tection des fissures dans les arbres de machines tournantes[END_REF], with constant max (liS+ II, IIS-11); so is the function b(y) with constant max (IIH+ II, IIH-11). It follows from the theory of differential equations that the solution y(t) of equations (4) exists and is unique. Therefore, the solution x(t) of equations (l) exists and is unique. 0 Equation [START_REF] Tondl | Some Problems of Rotor Dynamics[END_REF] expresses that the jump in the stiffness matrix must be normal to the plane interface xTn = 0. There is no tangential jump. The coefficient o measures the magnitude of the jump and hence the intensity of the non-linearity. Foro = 0, the linear case is recovered. The non-differentiability of f(x) across an interface is characteristic of the piecewise linear systems (l) and [START_REF] Gasch | Institution of Mechanical Engineers Conference[END_REF]. This means that all the non-linearity of such systems is concentrated at the interface, and implies that the equation cannot be linearized about the origin for the purpose of finding an approximate solution.

Positive homogeneity •off(x): similarity of solutions x(t)

The function f(x) is neither additive nor homogeneous, in contrast to a linear function. However, it remains half-spacewise (more generally conewise) additive and positively homogeneous. This has a major consequence for the solutions x(t).

Proposition. The half-spacewise linear elastic force f(x) is positively homogeneous of degree l:

f(..\x) = J.f(x), \f). >0, lfxE~'. (7)
As a consequence, ifx(t) is the solution of equations (I) with x(O) = x,, x(O) = v 0 , then J.x(t)

is the solution of equations (I) with x(O) = J.x,, x(O) = J.v 0 ,

M(J.i) + G(H) + f(h) = 0, \f). >0. ( 8
)
Proof Since W,_ and H;l' C are half-spaces, for any coefficient ,\ > 0 and any x e W,_ or x e H;l'C, one has ,\ x e H;l':-or ,\ x e H;l" • Therefore, in view of its definition [START_REF] Gasch | Institution of Mechanical Engineers Conference[END_REF], the function f(x) is positively homogeneous of degree !, as defined in equation [START_REF] Mayes | Institution of Mechanical Engineering Conference[END_REF]. Now, let x(t) be the solution of equations (l) for

x(O) ~ x.. i(O) ~ v 0 : i.e., Mx + Gi + f(x) ~ 0.
Multiplying this equation by any .l. > 0 gives, in view of equation [START_REF] Mayes | Institution of Mechanical Engineering Conference[END_REF], ,\Mx + .l.Gi + U(x) ~ M(.l.X) + G(H) + f(h) = 0. Therefore .l.x(t) is also a solution of equations (! ). Since the solution is unique, it must be the solution for the initial conditions [START_REF] Ismail | [END_REF] shows that solutions corresponding to homothetic initial conditions are homothetic (similar). This property is the most interesting characteristic of the halfspacewise linear systems under consideration and, more generally, of conewise linear systems, which is reminiscent of the superposition principle. In phase space, positive homogeneity of the solution with respect to the initial condition is expressed by y(.l.y,, t) ~ .l.y(y,, t), V.l. >0. It follows by integration that q(x, i) = txTMi + txTS(x)x ~ c, q(x, i) = 0, [START_REF] Masri | [END_REF] where q is a half-spacewise quadratic form representing the total energy of the system and c a motion constant depending on the initial conditions. Expression ( 10) implies that the system is conservative, and that the gyroscopic matrix G does not contribute to the system energy. In phase space, equations [START_REF] Masri | [END_REF] can be written as:

x(O) = h 0 , x(O) ~ !.v 0 • 0 Equation
q(y) = tyTB(y)y = c, B(y) = {B+, Vy E G;l~}• B_, Vye G;l _ q(y) = 0, 2.2.4.
Positiveness of q(y): stability of the equilibrium point y, = 0 (II)

Obviously, the only equilibrium point of system (I) and ( 2) is the origin y, = 0. From equations (!!), it follows that if S+ and S_ are positive definite, the half-spacewise quadratic form q satisfies q(y) > 0,

q(y) =0, Vy#O. ( 12
)
According to stability theory, q(y) is a Lyapunov function. Therefore the origin is a stable equilibrium point for the system. It should be pointed out that for gyroscopic systems (G # 0), the positiveness of S+ and S_ is only a sufficient condition of stability.

Having finished with the general features of conewise linear systems, we progress to a more detailed investigation of particular cases.

PIECEWISE LINEAR 1-DOF CASE

To begin with, let us consider a piecewise linear autonomous 1-DOF model:

Mx + K(x)x = 0, x(O) = x 0 , x(O) = v,, K(x)=K+>O, 'lx>O; K(x) = K_ > 0, 'lx.;:; 0.
The exact solution of equations ( 13) can be written as follows: for n =I, 2, ... ,

x(t) = A(x 0 , v 0 }p(l 0 , 1), p(t 0 ,t)= max (o.-1 -sinw+(t -t2n_ 2 )) 1€(/2n-2•12n-l) (I)+ + min (o. --1 -sin w_ (t-t2n_ 1 >). IE{I2n-J•12n) (I)_ t 0 = --tan ---, I _1 (w+x 0 ) OJ+ Vo 1t 1211-1 = 1211-2+-. w+ 1t l211 = 1211-I +-, w w_=JK_/M. ( 13
) (14a, b) (14c) (14d) (14e)
Here x 0 and v 0 are initial conditions (with x 0 > 0, without loss of generality). Figure 2(a) shows the solution as a function of time, and Figure 2(b) its trajectory in phase space.

Notice that the general solution x(l) =A (x 0 , v 0 )p(l 0 , 1) in equations ( 14) is a twice differentiable periodic piecewise sinusoidal function (the acceleration is continuous but not differentiable). Its period is given by T = (T+ /2) + (T_ /2) = (nfw+) + (nfw_) = njM{(ljjl(,.) + (1/JIC. )}. (15) An equivalent natural frequency w,, of the system can be defined as follows:

w,, = 2n/T = 2w+w_j(w+ + w_) = 2 JK:. JIC.IJM<JK:. + JlC. ). (16) It is obviously a dominant component of the spectrum of x(l) in the frequency domain. It was found numerically that if the system is excited by a sinusoidal force with frequency wrq' Mx + K(x)x =A sin (w,,l + 4> ), (17) then the maximum amplitude of the motion increases linearly to infinity (see Figure 3). Such a resonance phenomenon suggests that w,, plays a role similar to that of the natural frequency win a 1-DOF linear system. This is indeed the case when the piecewise linear system is excited by any function/(1) = Cft(1 0 , 1), where p(1 0 , 1) is defined in equations ( 14) and C and 1 0 are arbitrary coefficients:

Mx + K(x)x = Cft(1 0 , 1). ( 18
)
The forced response of the system is:

x(l) = (C/2M)Ip(1 0 , 1). (19)
1•0 1-1 0•5 " 0•0 0•0 ~ -0•5 -1•0 -1-5 -2•2 ~-----:"::------:"':----~:----::' -1•5 -1•0 -0•5 0•0 0•5 1•0 1•5 -2•2 -1-1 0•0 1-1 2•2
x, Clearly, x(t) diverges linearly as t-+ + oo. It is well known that under an excitation A sin (rot + </> ), the response of a linear system with natural frequency w is -Atcos(wt+</>)/2, which is comparable to equation [START_REF] Guckenheimer | Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields[END_REF] and to the graph in Figure 3(a).

COMMUTATIVE 2-DOF CASE

For further reference, let us recall some basic facts about linear systems. According to classical modal analysis theory, the general solution of a 2-DOF linear system is a linear combination of two fundamental solutions called linear modes. Each mode is composed of a harmonic eigenfunction e'w' and a corresponding eigenvector v ,p 0, defined by

( -ro 2 M + iroG + S)v = 0. (20)
Equation [START_REF] Lichetenberg | Regular and Stochastic Motion[END_REF] results in real eigenvalues w and real eigenvectors v if G = 0 and complex conjugate ones if G ,P 0. The system launched in the direction of an eigenvector oscillates in a sinusoidal motion on this eigenvector. In configuration space, its trajectory is a segment if G = 0 (see Figure 4(a)) and an ellipse if G ,p 0 (see Figure 4(b)).

Most of the above characteristics of linear systems can be extended to a particular category of conewise linear systems; namely, non-gyroscopic ones with commutative conewise constant stiffnesses. Indeed, Jet G = 0 (thus S = K) and M = mi 2 for simplicity (for M ,P ml 2 , it suffices to replace K+ and K_ everywhere by M-1 K+ and M-1 K_ or, better, by M-1 1 2 K+ M-1 1 2 and M-1 1 2 K_ M-1 1 2 to preserve their symmetry). Assume further that K+ and K_ are not only positive definite and continuous, but also commutative: (21) The commutativity condition (21) and the continuity condition (6) yield K+ DDT= nD'K+ and K_DDT = DDTK_, which imply that K+ D = (DTK+D)D and K_D = (DTK_ D)D, respectively. It follows that D is an eigenvector of both K+ and K_. Due to their symmetry, the other eigenvector is t, orthogonal to D. Therefore, V = (D, t) is the real orthogonal matrix which diagonalizes K+ and K_ simultaneously, K_ =VL_ ¥', [START_REF] Thompson | [END_REF] L+ and L_ being diagonal matrices containing the eigenvalues (w:', w,+') of K+ and (w;;', w,- 2 ) of K_, respectively. It can be shown that L+ and L_ are related through the continuity relation L+ = L_ + oe,ef, e, =(I, 0)\ (23) which implies that w:' = w;;' + /i and w,+' = w,-2 = w;. In this special case, the system equations (I) can be uncoupled into two independent equations by means of the linear transformation x = Vz, = (z., z,)L z, + w(z,)z, = 0, zl + w;z, = 0, [START_REF] Thompson | Nonlinear Dynamics and Chaos[END_REF] where w(z 11 ) = w;t 2 , Vzn > 0, ro(z 11 ) = w; 2 , Vzn ~ 0. The motions in the directions Z 11 and Zr are independent (see Figure 5(a)). According to the results for the 1-DOF system, the first mode z, = A,p(t 0 , t) is piecewise sinusoidal as described in equation ( 14), whereas the second mode z, =A, sin (w,t + </>) is sinusoidal (see Figure 5(b)). It follows that any solution x(t) can be represented by a linear combination of nz,(t) and tz,(t), x(t) = nA,p(t 0 , t) + tA, sin (w,t + </>,), [START_REF] Thompson | [END_REF] where A., A, t 0 and </>, are real numbers determined by initial conditions. Inspection of equation ( 25) reveals that the periodic functions p(t 0 , t) and sin (wt + 4>) play the role of two eigenfunctions associated to the two orthogonal eigenvectors n and t. These eigenfunctions and eigenvectors are fundamental, since the former determine the natural frequencies of the system and the latter uncouple its equations. It is straightforward to extend these properties to MDOF commutative half-spacewise linear systems. It can be concluded that the "rectilinear mode" concept remains meaningful for these particular non-linear systems.

In other terms, modal decomposition is not strictly restricted to linear systems, which are both additive and homogeneous. It also applies, in the specific form [START_REF] Thompson | [END_REF], to half-spacewise additive, positively homogeneous, commutative systems.

Similar motions were studied by Rosenberg in a slightly different context [35][36][37][START_REF] Rosenberg | Advances in Applied Mechanics[END_REF]. He investigated discrete, undamped, strongly non-linear, symmetric, smooth systems satisfying --'-::---L-----,L where x E IR" and U(x): IR"--+ IR is a twice continuously differentiable potential energy function satisfying the central symmetry conditions U( -x) = U(x) and U(O) = 0. He was able to characterize the periodic solutions by the following properties: all co-ordinates vary according to a periodic function (not necessarily a sinusoidal one); all displacements and velocities vanish simultaneously once every half period of the motion; all co-ordinates of the system can be expressed, at any given time, as simple-valued functions of any one co-ordinate. He called such periodic solutions normal modes. Notice that the periodic motions z.(t) and z,(t) of a half-spacewise linear commutative system satisfy all the above conditions, despite the fact that the system is not symmetric with respect to the origin and non-smooth. To distinguish these periodic motions from the normal modes of Rosenberg, we refer to them as non-linear modes for reasons that will be further explained later.

• x + iJU(x)jiJx = 0, ( 26
) 0•4 (a) 0•2 0•0 -0•2 -0•4 -0•6 C..,.--L.,

GENERAL CONEWISE LINEAR MDOF SYSTEMS

Motivated by the results of the previous sections and inspired by Rosenberg's normal mode concept, we now look for the periodic motions of more general conewise linear systems. We aim, for example, at answering the following questions: what do the periodic solutions of a general (non-commutative) half-spacewise linear 2-DOF system look like? How many are there? What are their characteristics? Since, for non-commutative matrices, s+ and s_, there is no regular matrix v simultaneously diagonalizing them, the system equations (I) cannot be uncoupled (even when G = 0), and the global free response of general systems can no longer be obtained analytically. A numerical approach must be used instead.

First of all, we characterize a periodic solution of MDOF conservative systems (I) with the following five conditions of periodicity and energy conservation: i(Xo, X,, T) = io, q(x., io) =c. (27) Here Xo and X, are the initial conditions inducing the periodic solution, Tis its period, q is the quadratic form of the system energy and c its value for Xo, i 0 • By introducing an augmented phase vector a, and a corresponding residual vector r as [ x(x0 , X,, T)-X 0 ] r= i(x 0 ,i, T)-io , q(x 0 , X,)-c (28) and by using the phase space vector y defined in equation ( 5), we can put conditions (27) in the form

r(a) = r(yo, T) = (Y(Yo(• T))-Yo)= O. q Yo -c (29)
Then, we develop an efficient numerical algorithm based on Newton's iteration method and Newmark's integration technique to extract all the periodic solutions at the predetermined energy level c. More specially, the solution of equations (29}-a system of five non-linear equations in five unknowns in a 2-DOF case-is obtained by applying the Newton's recurrence formula (30) wherej is the iteration counter, rj and (orfoa)j denote the residual vector and the Jacobian matrix evaluated at aj. At each iteration, Newmark's implicit unconditionally stable finite difference method (trapezoidal rule),

k = I, ... ,N, NLJt = Tj, ( 31
)
is used to estimate the phase vector yj = y(aj) after an approximate period P, starting from y 0 • The difficult part of this algorithm lies in the calculation of the Jacobian matrix (cf., Appendix A). Finally, the algorithm (30H31) gives us the periods T and the initial conditions Yo of the periodic solutions.

To demonstrate the capacity of the algorithm and study the periodic solutions, we use the 2-DOF rotating system shown in Figure I as a gyroscopic model. A nongyroscopic system is obtained from the same model by setting the rotation speed to zero

(Q = 0): non-gyroscopic system, Mx + K(x)x = 0; gyroscopic system, Mx + Gx + [K(x)-mQ 2 1 2 x] = 0.
M, G and K(x) are given in equations [START_REF] Grabowski | [END_REF]. The next four graphs illustrate the periodic solutions of these two systems. Figure 6 presents their trajectories in configuration space, Figure 7 shows the phase portraits of x 1 of the solution with longest period, Figure 8 shows the time history of such a solution, and Figure 9 depicts the trajectories of two periodic solutions in the configuration space for different constants of motion.

Besides assessing the algorithm capabilities, the numerical study has revealed the following interesting characteristics of the periodic solutions of conewise linear systems.

The trajectories of the periodic solutions of non-gyroscopic systems at a given energy level are open curves in the configuration space (instead of straight lines in the linear case), which neither pass through the origin nor are othogonal at their intersection in general (in contrast to the normal modes of Rosenberg) (see Figure 6(a)). The curvature comes from the non-linearity, the intersection offset from the lack of the central symmetry, and the loss of orthogonality is due to the lack of commutativity S+ S_ # S_ S+. The time history of the two components of the first periodic solution (Figure 8(a)) shows that they do not vanish simultaneously, but that they do reach their maximum and minimum at the same time. In phase space, the orbits are closed curves symmetric with respect to the velocity axis, instead of ellipses (see Figure 7(a)).

The trajectories of the periodic solutions of gyroscopic systems are closed curves (instead of ellipses) (see Figure 6(b)). The closure is due to a phase difference between the two components of the solution, as can be seen in Figure 8(b ). The deviation from ellipses, both in configuration space and in phase space, comes from the non-linearity (see Figures 6(b

) and 7(b)).

The periods of the periodic solutions are constant and independent of the energy level (i.e., of the vibration amplitude) as for linear systems. This is because the nonlinearity of the systems considered is concentrated at the interface passing through the origin.

The trajectories of the periodic solutions with the same period but at different energy levels form a set of homothetic curves (see Figures 9(a) and(b)). This property is due to the positive homogeneity of degree I of the solution with respect to the initial conditions discussed before.

For weak non-linearities s+ = s_ + Jn n', J <I! II s_ II, the number of structurally stable periodic solutions with different periods is equal to the number of DOF of the system, just as in the linear case. By "structurally stable" we mean that there exists no infinitesimal perturbation of the system parameters that alters the qualitative behavior of the system solution [START_REF] Parker | Practical Numerical Algorithms for Chaotic Systems[END_REF].

1•100 1-4 (a) (b) 1•0 0•525 0•6 0•2 " 0•050 " -1•200 L_ _ _l _ _ __L_L__L _ _ __J -1-8 L___L,--J.
,-,---L_l::-::--'::c:'-!:-c:':--:-' The characteristics of periodic solutions listed above show a close resemblance between conewise linear systems and linear systems. They suggest the following definition.

-1-400 -0•825 -0•250 0•325 0•900 -2•2 -1•8 -1•4 -1•0-0•6-0•2 0•2 0•6 1•0 x,
Definition. Any structurally stable periodic solution of a conewise linear system is called a non-linear real mode when G = 0, and a non-linear complex mode when G # 0. By analogy with linear modes, the periodic solutions with the same period and similar trajectories but different energy levels are said to belong to the same non-linear mode. Each non-linear mode can be considered to have a generalized eigenvalue that is inversely proportional to the period of the solution, and a generalized eigenfunction describing the non-linear modal curve.

1•0 (a) 0•5 0•0 -0•5 -1•0 -2•5 , 1•5 (b) 1•0 0•5 0•0 -0•5 -1•0 -1•5 -1•5 ' ' ___ .... ' ', / / / ... ____ _ -2•0 -1•5 -1•0 -0•5 -1•0 -0•5 0•0 ' ' 0•0 0•5 0•5 ' \ \ \ ' 1•0 1•0
0•5 ~ 0•0 c ~ ' ' -0•5 ' 8 / / • -a -1•0 6 -1-5 (a) 0 I 2 3 4 5 1•0 0•5 ~ c 0•0 a -0•5 ~ y • -1•0 ! -1-5 -2•0 -2•5 (b) 0 2 8 10
Time, t Figure 8. Time histories of the periodic solution with longest period: --, x 1 ; --, x 2 • (a) Non-gyroscopic systems; (b) gyroscopic system.

ANALYTICAL ESTIMATE OF NON-LINEAR MODE PERIODS

Having located and characterized the non-linear real and complex modes with constant periods, we proceed to formulate an analytical estimate of these periods in terms of system parameters, using the 1-DOF system analysis as a guide. This estimate provides a simple yet precise index for crack initiation and propagation in rotating shafts.

The idea is to find an equivalent n-DOF linear system, the linear mode periods of which are equal or close to the periods of the periodic solutions of the n-DOF conewise linear system:

{ Mi + Gx + S(x)x = 0} {T, = Znjw"''} { Mi + Gx + s .. x = 0 } (JZ) x(t) = x(t + T,) = > 1-ro;,M + iro,,G + S,,l = 0 . z = 1, n 1•5 1•5 (a) (b) 1•0 1•0 0•5 0•5 0•0 ~ 0•0 ~ -0•5 -0•5 -1•0 -1•0 -1•5 -1•5
-2•0 ~.:2:~:::::::.-:-::----c:-:!---:'-:--.,J

-2•0 -1•5 -1•0 -0•5 0•0 0•5 1•0 1•5 -3•2 -2•2 -1•2 -0•2 0•8 1•8
x, It was found in section 3 (equations ( 16)) that the equivalent natural frequency w,, of a 1-DOF PL system can be written as a function of the system parameters M, K+ and K_ as w,, = 2n/T = 2w+ w_/(w+ + w_) = (2/jM){(I/A) + (1/JfC)}. [START_REF] Shaw | The Supression of chaos in periodically forced oscillators[END_REF] This formula shows that the stiffness of the equivalent linear 1-DOF system is

K,, = Mw;, = 4{(1/A) + (1/JfC)) -2 (34)
This formula can be straightforwardly extended to a MDOF non-gyroscopic conewise linear system, for which G = 0 and S(x) = K(x) is the conewise constant stiffness matrix with K+ and K_ positive definite. Indeed, since the square root of a positive definite matrix is well defined, we can determined S,, by analogy with equation ( 34) as:

S,, = K,, = 4(jiC -I+ jiL -l)-2 (35)
The spectral decompositions of K+ and K_ into their respective orthonormal eigenvector matrices V + and V _ (V: = V:;: 1 , V~ = V: 1 ) and diagonal positive eigenvalue matrices L+ and L _ yields Consequently, K,, in equation ( 35) can be computed as

K = 4 (V ~-lyT +V ~-IVT)-'. eq +VL+ + -VL- - (36)
Finally, the solution of the eigenvalue problem IK,,-w;,MI = 0 gives the equivalent frequencies the inverses of which are approximations to the periods of the non-linear modes according to equations (32). These period estimates are exact for commutative MDOF systems. To assess their accuracy in the non-commutative case, we consider at first a particular 2-DOF conewise linear system, where the natural frequencies in the two linear regions satisfy wj = 3w ;t and w 0 = 3w;-. In this case, it can be demonstrated that one of the non-linear modes has a period equal to T = (n/wi) + (n/wl). (37) In Figure I 0 it is shown that the approximate formula (36) provides a very good estimate of the non-linear mode periods for small non-linearity, while the algorithm (30) and (31) and exact formula (37) yield identical results. Incidentally, this last fact confirms that the algorithm developed in section 5 gives accurate results. The non-linear mode periods calculated from the analytical estimate (36) are also compared with the numerical results ---,---,--~----.---, Coefficient of non-linearity, a Without losing generality, we examine the eigenvalue problem in the positive region. It can be shown that the 2n generalized eigenvalues of the system (39) obtained from IA+-.<8+1=0 are purely imaginary and conjugate and, thus, come in pairs: (jw;', -jwi), ... , (jw:, -jw,;"). They correspond to n pairs of complex conjugate eigenvectors (u;', iii'), ... , (u:, ii,;") which form a basis for the n-dimensional complex space. Moreover, the real and imaginary parts of these eigenvectors constitute on orthogonal basis u+ = (ut', uti, ... ' u:', u:'), u~ i= V~1 for the 2n-dimensional real space satisfying [START_REF] Parker | Practical Numerical Algorithms for Chaotic Systems[END_REF] In equation [START_REF] Parker | Practical Numerical Algorithms for Chaotic Systems[END_REF], I, is the 2n-dimensional identity matrix. We next perform a linear transformation on equation ( 39), [START_REF] Rudiger | From Equilibrium to Chaos-Practical Bifurcation and Stability Analysis[END_REF] and pre-multiply it by U~ to obtain U~A+ U+z+ = ~B+ U+:i+. According to the orthogonality relations [START_REF] Parker | Practical Numerical Algorithms for Chaotic Systems[END_REF], the equation is reduced to :i+ = E~ z+, which yields a set of second order differential equations, z+ +W+z+ =0. [START_REF] Huseyin | Vibrations and stability of multiple parameter systems[END_REF] where W + is a positive definite diagonal matrix defined by [START_REF] Moon | Chaotic Vibration-An Introduction for Applied Scientists and Engineers[END_REF] The same procedure can be carried out in the negative region, which results in z_ + W_z_ =0.

7-00,

(45)

Applying the inverse transformations z+ = U:;:'y and z_ = U='y to equations ( 43) and (45) respectively, yields a standard non-gyroscopic equation: y+ R(y)y = 0, R( )= {R+ "'U+ W+U:;:', y R_ ,U_ W_U='•

The system equivalent to system ( 46) is then defined as y + R,,y = 0,

lfyE ~~} ltye ~"' • (46) (47)
where R,, is a function of R+ and R_. Since W + and W _ ~ositive definite, the inverses of the square roots of R+ and R_ are equal to U+v'W+ -'u:;:' and U_~-'U='

respectively. Therefore, R,. can be obtained by analogy with equation (35) as

R =4(U f\V-Iu-'+U f\V-Iu-')-2 (48)
eq +V"+ + -V""'--

Finally, it can be shown (see Appendix B) that the following two eigenvalue problems possess the same distinct eigenvalues:

1-w'I,+R,,I=O = 1-w'M+iwG+S,,I=O.

(49)

As the dimension of the first eigenvalue problem is twice that of the second, each of its distinct eigenvalues is repeated. 'a 4

~ : f-><--X--><--X--X--X--X--~--><--><--><-~--X--l<--X--X-*-X--><-1

0•0 0•4 0•8 1•2 1•6 2•0
Coefficient of non-linearity, 8

Figure 12. Analytically estimated and numerically calculated periods of gyroscopic system. Analytical:-, first mode; ---, second mode. Numerical: Q, first mode; x, second mode.

To validate this formulation, we have re-examined numerically the previous 2-DOF gyroscopic model (Figure I) rotating at Q = 8 radfs, the other parameters being given by equation [START_REF] Rosenberg | Advances in Applied Mechanics[END_REF]. The estimates of the non-linear modes periods have been obtained for the same interval of the non-linearity coefficient J E [0, 2]. The analytically estimated and numerically calculated periods are presented in Figure 12. The coincidence of the two methods is even better than in the non-gyroscopic case. Note that the analytical estimate can be used as an initial guess for the numerical algorithm, which usually results in a drastic reduction in the number of iterations.

STABILITY AND BIFURCATION OF CONEWISE LINEAR SYSTEMS

The restriction on the magnitude of the non-linearity that we experienced in using the analytical formulation (36) for non-gyroscopic conewise linear systems requires further study of their characteristics. It is well known that a strongly non-linear system suffers important qualitative changes, such as instability and bifurcation of solution, chaotic motion, etc. In order to trace the entire evolution of the system characteristics as the non-linearity increases and to determine the extent to which our analytical estimate is adequate, we use again the 2-DOF system [START_REF] Rosenberg | Advances in Applied Mechanics[END_REF] to examine the stability of non-linear real modes in phase space for each value of b. For completeness, two methods, namely characteristic multipliers and the Poincare map, are employed, although they are known to be equivalent through Floquet theory [START_REF] Guckenheimer | Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields[END_REF].

CHARACTERISTIC MULTIPLIERS

The stability of any non-linear mode is determined by its characteristic multipliers, also called Floquet multipliers. As a start, let us briefly recall why.

The Lyapunov definition of stability of a periodic solution y(l, y~) requires that the distance between the non-linear mode and the perturbed motion, lld(l)ll = lly(y~ +do. 1)-y(y~' 1)11 (51) The 2n eigenvalues of the monodromy matrix [START_REF] Parker | Practical Numerical Algorithms for Chaotic Systems[END_REF][START_REF] Rudiger | From Equilibrium to Chaos-Practical Bifurcation and Stability Analysis[END_REF] (52) are referred to as characteristic multipliers. They play an important role in deciding whether or not the initial perturbation d 0 decays or grows as t-+ oo (k-+oo). It can be observed that if 3k E [I, 2n] for which IJ..I > I, the non-linear mode is unstable.

• (a) i--o. ., 1-4 = ~ ~ 1•3 ~ 1•2 ~ ~ -------- = 1-1 Bifurcation point r ~ •~ 1•0 ~/ :;; 0•9 ' ' 0•0 0•4 0•8 1-2 1•6 2•0 7
Notice that the monodromy matrix 8yj8y 0 associated with a non-linear mode in expression (52) is a part of the Jacobian matrix 8r/8a introduced in equation ( 30) and evaluated at the last iteration of the numerical procedure presented there (see Appendix A for a justification). Therefore, after having located a non-linear mode by numerical iteration, it is straightforward to examine its stability.

For the previously discussed system (Figure I, equations ( 38)), the variations of the maximum amplitudes and periods of the non-linear modes are plotted against the non-linearity coefficient bin Figures 13(a) and (b) respectively. For small b, two non-linear modes determined numerically are found neutrally stable. (In this particular case, the second mode period varies very little with b.) At a critical value b, one of the modes becomes unstable and a third periodic solution suddenly appears.

Im

Im Im

2l 0•2 < • ~ 0•0 0 -0•2 ! -0•4 -0•6 -0•8 0 2 3 5 Time, t 1•5 (b)
--------~-:_"':..":.."':,.":::."

1•0 / - / / ~' 0•5 / ! I \ I ' 0•0 r r ' ' > -0•5 \ I ' I ' -; ' ' --- -1•0 ' -/
,,_ ___ :.:::::.:..:...-- The instability phenomenon can be interpreted in the eigenvalue complex plane as schematized in Figure 14. For small o, all the characteristics multipliers of the mode are located on the unit circle. As o increases, a pair of multipliers relative to one mode moves on the unit circle towards -I. At a critical value b, = 0•96 this pair merges at -l and o > o, one eigenvalue moves inside and the other outside the circle, on the real axis. It follows from the stability criterion that the corresponding non-linear mode becomes unstable afterwards. Differentiating equations (54) with respect to the entering positions yields: (oz,._ J!oz,._,) = P(A+ ,t"'-') + P(A+ ,t"'-')z,.(ot"'-2 /oz,._,) = z':;-2 , (56a) (oz,.foz,._,) = P(A_ ,t2k-') + P(A_, t"' -')z,._, (ot"'-'/Bz,._,) = z"-'• (56b) Upon compiling these leaps, the Jacobian matrix DP of the Poincare map takes the form

-}.5 -0•8 -0•6 -0•4 -0•2 0•0 0•2 0•4 0•6 Displacement 1•2 (c) 0•8 / / 0•4 T. / / ;\~ / / / " 0•0 / " -0•4 Bifurca~ :~~ \ ____ .......... -0•8 -1•2 -H -1•0 -0•6 -0•2 0•2 0•6 1
DP = z~-~z~-2 ••• z'1!.-'z':-2 ••• z~.
(57) After localizing the travelling times t,._, and t2k-l and intercepting points z,._ 2 and z,._, (k = l,N), the matrices P(A+,t"-2 ), P(A_,t"-1 ), z':;-2 , z':;-1 and DP can be successively evaluated. They are explicitly given in reference [START_REF] Zuo | 6tection des fissures dans les arbres de machines tournantes[END_REF].

The same criterion of stability as for the characteristic multipliers is then applied to the eigenvalues ofDP. To return to the example [START_REF] Rosenberg | Advances in Applied Mechanics[END_REF], exactly the same instability phenomenon is observed as with multipliers. Another neutrally stable mode with double period appears. The maximum eigenvalue of the bifurcating mode is replotted against ~ in Figure 17 to show the perfect agreement between the two methods.

To recapitulate, a 2-DOF conewise linear system has two neutrally stable non-linear modes below a non-linearity threshold(~ < ~,) and it can have three or more beyond it.

CONCLUSIONS

Mechanical systems involving unilateral contacts (whether microscopic or macroscopic)

give rise to conewise linear models. In particular, half-spacewise linear models have been widely adopted, as a first order approximation, for describing the transverse vibrations of cracked beams (non-gyroscopic case) and rotating shafts (gyroscopic case) (1)[START_REF] Gasch | Institution of Mechanical Engineers Conference[END_REF][START_REF] Grabowski | [END_REF](4)[START_REF] Zuo | 6tection des fissures dans les arbres de machines tournantes[END_REF]. The present study sheds some light on the multi -degree-of-freedom models. The main results are summarized here.

(i) The existence and uniqueness of the solution has been established on the grounds of a continuity condition at the interface and the similarity of solutions with similar initial conditions (at different energy levels) has been proved on the basis of the positive homogeneity of the governing force.

(ii) Below a critical non-linearity or bifurcation threshold (~ < ~,). conewise linear systems bear a strong resemblance to linear systems, in terms of their fundamental periodic solutions herein called non-linear modes whether real (non-gyroscopic case) or complex (gyroscopic one): their number is equal to the number of DOF, their periods are constant and energy independent, their trajectories are similar curves at different energy levels, etc.

(iii) Beyond the bifurcation threshold (b > bJ, a bifurcation of the second non-linear mode occurs and the resemblance to linear systems is reduced. The threshold is much higher in the gyroscopic case (b, > 2, [START_REF] Zuo | 6tection des fissures dans les arbres de machines tournantes[END_REF]) than in the non-gyroscopic one (b, = 0•96) due to the well known stabilizing effect of rotation [START_REF] Rudiger | From Equilibrium to Chaos-Practical Bifurcation and Stability Analysis[END_REF].

(iv) The analytical estimate of the non-linear mode periods is very accurate below the bifurcation threshold and, hereby, could become a valuable tool in engineering practice. Indeed, for crack detection in rotating shafts, it gives the precise relationship between the natural frequencies and the system parameters, and thus shows how the former will vary with the crack depth [START_REF] Zuo | 6tection des fissures dans les arbres de machines tournantes[END_REF]. Note that even in the non-gyroscopic case, the bifurcation threshold is situated far enough to interfere with crack detection concerns.

1 2 •Figure I .

 2I Figure I. The simplest gyroscopic model. X, fixed reference; x, rotating reference.

3 .

 3 Conservation of energy q: motion constant c Pre-multiplying equation (I) by iT yields, upon taking the anti-symmetry of G into account (i.e., xTGx ~ 0), xTMx + xTS(x)x = o.

Figure 2 .

 2 Figure 2. Free response of l~DOF piecewise linear model. (a) Displacement (--), velocity (-0-) and acceleration (---); (b) phase portrait.

20

 20

Figure 3 .

 3 Figure 3. Resonance of 1-DOF model. (a) Displacement history; (b) phase portrait.

Figure 4 . 2 -

 42 Figure 4. 2-DOF linear model:--, first mode;--, second mode; -0-, maximum equipotential. (a) Real non-gyroscopic modes; (b) complex gyroscopic modes.

Figure 5 .

 5 Figure 5. Non-linear modes of 2-DOF commutative conewise linear model: --, z,.(t); -, z,(t); -0-, maximum equipotential. (a) Configuration space; (b) time histories.

Figure 6 .

 6 Figure 6. Periodic trajectories in configuration space: --, period T 1 ; --, period T 2 , -0-, maximum equipotential. (a) Non-gyroscopic system; (b) gyroscopic system.

Figure 7 .

 7 Figure 7. Phase portraits of the periodic solution with longest period: --, (x 1 ,v 1 }; --, (x 2 ,v 2 }. (a) Non-gyroscopic system; (b) gyroscopic system.

Figure 9 .

 9 Figure 9. Similarity of periodic solutions in conewise linear systems: --, period T 1 ; --, period T 2 ; -0-. maximum equipotential. (a) Non-gyroscopic system; (b) gyroscopic system.

Figure lO .

 lO Figure lO. Comparison of theoretical, numerical and estimation results:--, theoretical values; Q, numerical data; ---, analytical estimates. w t = w 1 +b.

~

  14r-----~----.------r-----,-----, 112 f-"'-o..,. ~ tor---~~~~~~~~~~~~~~ ~ 8 § 6 c

  to zero when 1 _, oo, for an arbitrary small initial perturbation lld(O)II <{ IIYt 11. A Taylor expansion of d(1) yields d(l) = {8y(y~, 1)/ayndo plus terms of higher order.Because of the periodicity of y, y(y~, 1 + T) = y(y~, 1) for 1 = kT, one has d(kT), {8y(y~, T)/8y~)'cJo.

Figure 13 .

 13 Figure 13. Bifurcation diagram of non-linear modes: -0-, first mode;---, second mode;--, bifurcated mode. (a) Amplitudes; (b) periods.

Figure 14 .

 14 Figure 14. Eigenvalues ofmonodromy matrix in complex plane. (a) Before bifurcation; (b) bifurcation; (c) after bifurcation.

-0 x, Figure 15 .

 15 Bifurcated non• linear modes. (a) Time history:--, x 1 ; ---, x 2 • (b) Phase portrait:--, (x 1 , vd; ---, (x 2 ,v 2 ). (c) Configuration space: ---, first mode; --, second mode; -+-, bifurcated mode; -0-, maximum equipotential.

8 Figure 17 .

 817 Figure 17. Maximum eigenvalues of equations (22) and (28): -, by characteristic multipliers; -Q, by Poincare section.

ACKNOWLEDGMENT

We wish to thank Mrs M. Schatzman for showing us how to shorten the proof of the second proposition concerning the positive homogeneity of the solutions.

for the 2-DOF conewise linear model sketched in Figure 1 when the disc is at rest. The system parameters used are as follows:

The periods calculated by the analytical estimate and the numerical iteration are plotted in Figure 11 vs. the non-linearity coefficient 0 E [0, 2]. Again, for a small amount of non-linearity (J <A= 1), a perfect agreement between these two methods can be observed. In this particular case, the period of the second mode varies very little with the non-linearity coefficient J.

GYROSCOPIC CASE

Now we consider the system with a gyroscopic anti-symmetric matrix G. In this case, the equivalent natural frequencies of the system depend not only on S(x) and M, but also on G, which is neither symmetric nor positive definite. As a consequence, S,, defined in (32) can no longer be obtained through (36). To establish a correct estimate for S,,, we put equation (1) in the form [START_REF] Whin | [END_REF] where the indexes "+" and "-" denote the two regions of linearity as usual, B± = B~ are defined in equations ( 11) and S±)--AT -G -±• [START_REF] Rockafellar | Convex Analysis[END_REF] The numerical results show that the additional periodic solution which appears after the critical point b, has a period equal to the double of that of the unstable periodic solution.

This newly born periodic solution traces a different trajectory in phase space, but keeps exactly the same characteristics as the first two non-linear modes. Its time history is shown in Figure 15(a), its phase portrait in Figure 15(b), and its trajectory in configuration space togeher with the other two non-linear modes in Figure 15(c). Such a bifurcation with period doubling as () = b, is a typical phenomenon, known as the flip bifurcation, which occurs in many non-linear dynamic systems [START_REF] Guckenheimer | Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields[END_REF][START_REF] Thompson | Nonlinear Dynamics and Chaos[END_REF][START_REF] Rudiger | From Equilibrium to Chaos-Practical Bifurcation and Stability Analysis[END_REF][START_REF] Moon | Chaotic Vibration-An Introduction for Applied Scientists and Engineers[END_REF]. Further examination of characteristic multipliers shows that for large non-linearity (() > bJ, the unstable non-linear mode continues to coexist with the other two modes.

.2. POINCARE MAP

The stability of periodic solutions can also be discussed from a more geometric point of view. Instead of examining the continuous trajectory of the system in phase space, one looks at a time-sampled sequence z(t 1 ), z(t 2 ), ••• , z(t,), .. . , abbreviated by z 1 , ••• , z,.

These points are obtained by intercepting the trajectory with a hyperplane called a Poincare section. The stability of a periodic trajectory can then be discussed in terms of the stability of the corresponding Poincare map.

Since the definition of this map relies on the knowledge of the solution of the differential equation ( 4), the sequence z, cannot be computed unless a general solution of this equation is available. In our autonomous conewise linear systems, the trajectory in any one cone can be obtained analytically, as long as the entering condition is known. Thus, for establishing the Poincare map for a half-spacewise linear MDOF system, it is natural to use the hyperplane 1: across which the system stiffness is discontinuous as a Poincare If the solution point hits the hyperplane 2N times within a period (see Figure 16), then it has travelled N times in each half-space. In view of the closed form solution in each halfspace, three successive points z,. _ 2 , z,. _ 1 and z,. on the Poincare section are related through two square 3 x 3 matrices [5], z,. = P(A_, r"-')z,. _,,

where the travelling times can be found in terms of the entering positions z,._, and z,._ 2 :

r"'-'=t"'-'(z

Figure 16. Poincare section and trajectory.

APPENDIX A: CALCULATION OF JACOBIAN MATRIX (8r/8a) AND MONODROMY MATRIX 8y(yt, T)/8yt

The definitions of r, a and q in section 5 imply the following equalities: iJy 0 fiJy 0 =I,. E R"' ""'; i!y 0 fiJT = 0 E R"'" '• The Jacobian matrix in equation ( 30) takes therefore the form or= ~(y(y0 , n-Yo)= (iJy(yo. T)fi!Yo-I,. i!y(yo. nfo"J\.

Moreover, in equation (A I) i!y(yo, T)/i!T = dy(yo, T)fdT = Y<Yo• n and from equation (II)

If we denote the monodromy matrix as i!y(yo , T)fi!Yo = J E ~;!'" x ,., the Jacobian matrix (A I) becomes:

or (J(yo. n-I,.

In order to obtain J, we first differentiate equation ( 4) with respect to y 0 to obtain

and then solve simultaneously this matrix equation and the system equation ( 4):

Here we use the fact that J(O) = iJy(y 0 , O)/ily 0 =I,. E !;!,. x ,._ H(y) is defined in equation [START_REF] Zuo | 6tection des fissures dans les arbres de machines tournantes[END_REF].

Integration of equation (A7) allows us to evaluate the square matrix (A4), the column vector (A2) and the line vector (A3) at time T and the initial condition y 0 • The Jacobian matrix (AI) follows. For this purpose, we employ the unconditionally stable trapezoidal rule (31). We do not use the closed form solution in each linear region because the number of times a periodic solution crosses the border line is unknown; therefore, the switching points are not defined a priori.

Finally, with the predicted values y 0 , and T, we obtain the Jacobian matrix (AI), as well as the residual r = y(y 0 , T)-Yo for each iteration of equation (30). Once the initial condition leading to a periodic solution yt and its period T are determined, the last evaluation of equation (A4) at y6 and T gives the monodromy matrix: The n-DOF gyroscopic system equivalent to the piecewise linear system (I) is defined as Mx + Gx + s,,x =0.

Another equilvalent 2n-DOF system represented by a non-gyroscopic equation is y+R.,y=O.

(BI) (B2)

We now demonstrate that the eigenvalues of system (B2) are double, and that its distinct eigenvalues are identical to that of the system (Bl).

Transforming where the matrix W is diagonal and positive definite: