Bronius Baksys

Nomeda Puodziuniene
email: nomeda.puodziuniene@ktu.lt

Modeling of vibrational impact motion of mobile-based body

Keywords: vibrations, impact motion, automated assembly

The paper deals with vibrational impact motion of a mobile-based body on an inclined plane, which is characteristic for vibrational alignment of components subjected to an automated assembly. The alignment model is presented by a twodegree-of-freedom vibro-impact system in which the moving body is impacting the plane obliquely. Properties of the vibrational impact motion under kinematical excitement of the moving body and under the excitement of the plane in two perpendicular directions were investigated by numerical methods. It was determined that impact displacement occurs under transient regimes of motion from system static until dynamic state of equilibrium. Dependencies of maximum displacement of the body and average displacement velocity on force of elastic resistance, amplitude of plane vibration in the direction of the displacement, and phase angle between the perpendicular members of vibration components were established. Zones of system and excitement parameters combination when the vibrational impact motion occurs were defined.

Introduction

Operational processes of many machines and technological equipments are based on systematic and repeating impact interaction of their components or parts. Sledgehammers, vibro-impact tools, impact dampers, riveting machines, screw driving devices, equipment for vibrational drilling and cutting, etc. operate under vibro-impact regimes. For the systems and equipment of this class, periodical regimes of motion that can be established according to forced or self-excited vibration schemes are characteristic. Components or end effectors of these systems move along the one line

A c c e p t e d m a n u s c r i p t

and straight focused impact on to a fixed point are characteristic for it. Systems where straight impact takes place and the velocities of moving bodies are directed to the impact point perpendicularly with respect to the contacting surfaces are called onedimensional impact systems.

Operation of some technological machines and equipment is based on impact motion of components or bodies for which non-unidimensional impact interaction occurs. Thus, non-unidimensional dynamic models are required to describe the main characteristics of these technological systems. Processes of vibrational transportation and vibrational alignment of components during automated assembly are possible to describe by a two-dimensional dynamic model. If the bodies under transportation do not lose contact with the surface of the transporter, or if the components undergoing alignment do not lose contact with each other, then the vibrational displacement occurs within the non-impact regime. However, under the certain conditions, the impact regime of motion occurs. In this case, such systems can be described by a two-dimensional impact dynamic model.

The dynamic system for the vibrational alignment of components during automated assembly has certain specific characteristics. The system consists of two interacting components; one of them is mobile-based. By exciting one of the components using vibration, the mobile-based component moves toward the mating of joining surfaces. Because of this motion, components undergoing assembly are aligned and conditions of automated assembly are established. Vibrational non-impact motion of mobile-based body was explored in the previous paper [START_REF] Bakšys | Modeling of vibrational non-impact motion of mobile-based body[END_REF]. Under the determined excitement parameters, the mobile-based body loses contact with the surface of the immobilized body at certain instant in time. Vibrational impact motion regime occurs then, and the mobile component repeatedly and systematically impacts the other component. Such an impact system has certain specific features. First, the impact point is not fixed because the moving body is impacting different points. Therefore, the impact interaction requires a two-dimensional impact model. Also, the fixing elastic elements always oppose body displacement. Because of this property, the body may move from the position of static to the position of dynamic equilibrium.

To investigate the motion of the impact system, not only the motion of the body after lift-off from the contacting plane needs to be mathematically expressed, but also the interaction of the impacting bodies. During the impact, the kinetic energy of the contacting bodies' normalized motion is transformed into deformation energy. Because of the deformation, stresses emerge at the contact zone and stress waves from this zone

A c c e p t e d m a n u s c r i p t

propagate through the interacting bodies. Investigating the impact by means of continuous media theory, which can account for the altering contact stresses and strains during impact, the duration of the impact, and the wave phenomenon caused by impact, is possible. However, it is a complicated task to investigate such processes analytically finding a closed form solution, and it is usually successful only for simple configurations of the impacting bodies. Therefore, the more complicated tasks can be solved only numerically employing finite element method with explicit integration technique. Among the most popular commercial computer codes of this type is LS-DYNA 3D.

However, for the investigation of impact systems, more often the simpler stereomechanical impact theory is used [START_REF] Goldsmith | Impact: the theory and physical behaviour of colliding solids[END_REF]. It is based on the assumption that the duration of the impact is very short. The impact effect is described by the coefficient of velocity restitution. In this case, a detailed investigation of the impact process to determine kinematical characteristics of the impacting bodies is not needed, and the required dependencies are derived according to general propositions of rigid body mechanics. Impact is described by the impulse integral

( )dt t F I ∫ = τ 0
where F(t) is the impact interaction force. All other forces affecting the system have insignificant influence because of short impact duration. Therefore, the influence of sustained low-strength forces on the impact process can be neglected. Integral evaluation of the impact forces characteristic for the stereomechanical theory does not allow a detailed investigation of the stresses and strains emerging during the impact.

However, on the basis of this theory, it is possible to investigate the motion of the bodies subjected to the impact impulse if the external forces change slowly compared to the impact duration. The before-impact and after-impact body motions can be distinguished and related to each other by kinematical equations of impact. The essence of it is that centers of mass during impact do not change but the velocities do change stepwise. The change in velocity is determined by impulse and kinetic energy conservation laws.

To determine velocities after impact, the coefficient of velocity restitution must be known. Its value can be obtained by employing energy conservation and determining the ratio of the body kinetic energy after impact to kinetic energy dissipated during the

A c c e p t e d m a n u s c r i p t

impact. It is assumed that Newton's hypothesis is valid. According to it, the ratio of body velocity after impact (v + ) to velocity before impact (v -) is constant and does not depend on impact velocity or geometry of bodies but only on material properties of the interacting bodies, i.e.

- + = v v R ,
here R is the coefficient of velocity restitution. The velocity entirely recovers after the impact if R=1, when the impacting bodies are perfectly elastic. If the impact is perfectly plastic, then the velocity after impact is equal to zero and R=0. For real bodies,

1 0 < < R always.
Until now, vibro-impact systems with straight central impacts have been mostly investigated. For example, Babitsky [START_REF] Babitsky | Theory of vibro-impact systems and applications[END_REF] and Kobrinskii [START_REF] Kobrinskii | Dynamics of mechanisms with elastic connections and impact systems[END_REF] provided detailed analyses assuming that the dynamics of vibro-impact systems can be reduced to only oscillatory motion. A mathematical model of a vibro-impact system accounting for both oscillatory and progressive motion and capable of transferring high-frequency low-amplitude excitation into low-frequency high-amplitude response is developed in work [START_REF] Pavlovskaia | Modeling of vibro-impact system driven by beat frequency[END_REF].

Nonlinear dynamics of impacting oscillators has received considerable theoretical and experimental attention in the past [START_REF] Peterka | Dynamics of double impact oscillators[END_REF][START_REF] Oestreich | Analytical and experimental investigation of impact oscillator[END_REF][START_REF] Imamura | Dynamic behaviour in a vibro-impact mechanical system[END_REF].

A body moving under impact regime has two degrees of freedom, so its impacts to the surface of contact will be not perpendicular but oblique. In order to describe the oblique impact, it is assumed that normal components of impact velocity do not depend on tangential components, and velocities after impact v y + are related to velocities before impact v y -by the equation

- + - = y y Rv v .
Tangential components of the velocity can be calculated by two hypotheses [START_REF] Babitsky | Theory of vibro-impact systems and applications[END_REF].

According to the first hypothesis, tangential component of the velocity after impact v x + is related to the one before impact v x -by the equation

( ) - + - - = x x v f v 1 ,

A c c e p t e d m a n u s c r i p t

here f is coefficient of instant friction ( 1 0 ≤ ≤ f

).

According to the second hypothesis, tangential component of the impact momentum is proportional to the normal component. The coefficient of proportionality is the coefficient of dynamic friction. Then the tangential component of the velocity after impact is expressed by the equation

( ) - - + + - = y x x v R f v v 1 .
The second hypothesis relates normal and tangential components of impact impulse similarly as the dry friction (Coulomb) law relate normal pressure and friction force. In other words, the essence of the hypothesis is an assumption that impacting surfaces are interacting according to the dry friction law, and that this interaction under non-instantaneous and instantaneous forces is described by the same relationship. Both hypotheses are used in the paper to describe impact interaction.

Dynamic model, equations of body motion and impact

Vibrational motion under impact mode is characterized by a two-dimensional model that is similar to the non-impact motion model utilized in [START_REF] Bakšys | Modeling of vibrational non-impact motion of mobile-based body[END_REF] (Fig. 1). The only difference is that the body of mass m in certain time instances lifts-off the plane and moves along a certain trajectory above it and then impacts the plane. Impact motion can occur under kinematical t sin y y ω 0 = excitement of the body or plane in two perpendicular directions

( ) ε ω + = t sin A x 1 1 , t sin B y ω 1 1 = .
Impact motion is investigated under one or another case of excitement; therefore, both cases of excitement are depicted in the model scheme. The body will liftoff the inclined plane and the reaction force disappears when the amplitude of the kinematical or plane excitement in normal direction reaches a certain value. The impact motion mode then begins. The body motion above the plane is two-dimensional and directed along the x and y axes. Elastic (c x , c y ) and damping (h x , h y ) constraints are acting in these directions. Such expression of these constraints is conditional. The twodimensionality of automated assembly device constraints is determined by one or

A c c e p t e d m a n u s c r i p t

several elastic elements located along one axis (for example, the axis of the assembly device). For investigating two-dimensional motion of the body, assume that motion along the x axis has no influence on body interaction with elastic constraint c y , and similarly for motion along the y axis with elastic constraint c x . Only under these conditions will the elastic constraints function linearly.

The body motion above the plane under kinematical excitement during the interval between two impacts is described by equations:

( ) α β sin mg sin t F x c x h x m x x + = + + & & & , (1)
( )

y st y y y y c cos mg cos t t F y c h y m - - = + + α β & & & ; (2)
where ( )

t sin y c t F y ω 0 =
is the force of kinematical excitement inclined by angle β with respect to axis y; y st is a coordinate of static equilibrium.

Under excitement of the plane in two perpendicular directions by harmonic vibration with frequency ω, the equations of body motion on the plane have the following expression:

( ) α ε ω ω sin mg t sin mA x c x h x m x x + + = + + 2 1 & & & , (3)
( )

st y y y y c cos mg t sin mB y c y h y m - - = + + α ω ω 2 1 & & & . ( 4
)
Impact interaction of the body and the plane is expressed by impact equations.

Impact to the plane will be oblique because two-degree-of-freedom body motion is being investigated. For analyzing dynamic oblique impact systems, two hypotheses are used to evaluate impact. According to both hypotheses, before-impact ( ) For this impulse equation, an equivalent kinematical equation exists

- + - = y R y & & ;
where R is a coefficient of velocity restitution; + y & and - y & are velocity magnitudes.

A c c e p t e d m a n u s c r i p t

According to the first hypothesis, the tangential component of the impact impulse does not depend on the normal component and is characterized by a physical constant dependent on properties and condition of the impacting surfaces. Using this hypothesis, ( )

f x x - = - + 1 & & ; ( 5
)
where f is a coefficient of Coulomb friction ( )

1 0 ≤ ≤ f .
The second hypothesis is based on assumption that the tangential component of impact impulse is proportional to the normal component, i.e.

n t fI I =
. Then the equation for the tangential impact velocities is ( )

R y f x x + - = - - + 1 & & & . ( 6
)
The physical meaning of the second hypothesis is easier to understand, but for the approximate calculations the first hypothesis can also be used.

Using dimensionless variables and parameters, the equations for the body motion under kinematical excitement and impact equations as follows:

α µ β τ ξ γ ξ ξ sin sin sin h + = + + 2 1 & & & ; (7)
ν α µ β τ η λ η η - - = + + cos sin sin h 2 2 & & & ; ( 8
) ; R - + - = η η & & (9) ( ) - + - = ξ ξ & & f 1
(according to the first hypothesis); (10)

( ) R f + - = - - + 1 η ξ ξ & & &
(according to the second hypothesis). ( 11)

When the plane is excited in two perpendicular directions, the body motion on the plane is described by following the dimensionless equations: ( )

2 1 ω µ A g = .
Impact interaction of the body and the plane is described by the same equations ( 9),( 10),(11), as in the case of the kinematical excitement.

Impact motion of the kinematically excited body

First of all, the solution for the body impact motion under periodic regime is determined by investigating a simplified system. It is assumed that this system has no damping, and that there is no elastic element in the direction of the ξ coordinate. The motion of such a system between the impacts is described by the equations:

( ) α µ β ϕ τ ξ sin sin sin + + = & & ; (14) ( ) , cos cos sin ν α µ β ϕ τ η λ η - - + = + 2 & & ( 15
)
and interaction of the body and the plane is described by equations ( 9) and (10). The angle of impact phase ϕ is also incorporated in the equations of motion.

Parameters of periodic motion are calculated when impact velocities are the same during each impact and only one impact occurs in one period. The instant of impact is selected as the initial time point. Then the initial and boundary conditions are: Using initial conditions of motion from ( 16), the integration constants can be calculated. Expressions for ξ and η then have the following form:

, , 0, , , c + + = = = = = η η ξ ξ η ξ ξ τ 0 ( 16
( )

β
λ ν α µ ϕ τ λ β λτ ϕ λ β η λ λτ ϕ λ β λ ν α µ η + - + - + + ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ - - + ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ - - + = - 1 1 1 1 2 2 2 & (19)
Differentiating equations ( 18) and (19), and using boundary conditions from (17), the relation between velocities before impact

- -η ξ & & ,
and velocities after impact

+ + η ξ & & , is determined as . cos cos cos cos cos sin sin cos cos ; sin 2 ϕ λ β λ π ϕ λ β η λ π λ ϕ λ β λ ν α µ η ξ α πµ ξ 1 2 1 2 1 2 2 2 2 - + ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ - - + + ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ - - + - = + = + - + - & & & &
Employment of impact equations ( 9) and (10) leads to: (20)

in s f α µ π ξ 2 = - & ; ( )( ) [ ] ( )( ) πλ λ λ πλ λ ϕ β λ ν α µ λ ϕ β πλ λ η 2 1 1 2 1 2 2 2
In this expression, the impact phase ϕ is unknown. It can be determined from equation ( 19) and the boundary condition ( )

0 2 = π η
. Then, from this condition and impact equation (9), the trigonometric equation for impact phase calculation is obtained.

C cos B sin A = + ϕ ϕ ; here ( ) ( ); R 1 cos cos B ; sin cos R A + = - = πλ β λ λ π β λ 1 2 ( ) ( )( ) R sin cos C 1 1 2 - + + = λ λ π ν α µ . From the above equation C B B C A A arctg + + - ± = 2 2 2 2 ϕ .
To have periodical motion of the body, the condition 0

2 2 2 ≥ + - B C A must be
satisfied.

An impact phase depends on certain stages of body motion in respect of the plane, on the sequence, and duration of these stages. During the excitation period, body can remain in the state of rest for a while; it can also slide forward and backward, and move above the plane. Under the steady state regime these stages have cyclic sequence.

Stages consisting the motion cycle and the character of body motion depend on system and excitation parameters. Besides, the period of body motion can be equal or multiple to excitation force period. These possibilities create multiple variations of body motion in vibro-impact systems. Therefore, the body motion after impact is ambiguous and the impact phase equation yields several solutions, each of which corresponds to a different regime of motion.

Let the average velocity of body motion be defined as

( ) β ϕ α πµ τ ξ sin cos sin 2 1 / + - = = f f L T ave & ;

A c c e p t e d m a n u s c r i p t

here T τ is period of the system motion.

In order to determine the characteristics of the impact motion of a body restrained by elastic and damping restraints in x and y directions, the motion equations ( 7) and ( 8) together with impact equations ( 9) and (11), obtained using the second hypothesis, are solved by numerical methods. The peculiarity of the investigated impact system is that because of a body motion the location of the impact point is changing.

Furthermore, the elastic resistant force is acting along the body motion direction.

Because of this situation for the body transition from static to dynamic equilibrium position, the transient regimes of motion are critically important. To investigate such regimes of motion analytically is quite a complicated task. The numerical solution was performed for the motion regime where one excitation period corresponds to one period of body motion.

Examples of how the body coordinates ξ and η change during the motion are presented in Figures 2(a) and 2(b). Elastic force resisting to body motion on the inclined plate is directly proportional to parameter γ. When the stiffness coefficient of the elastic constraint in the direction of motion is increasing, the coordinate of the system dynamic equilibrium location is decreasing. Therefore, the maximum displacement of the body max ξ is also decreasing (Fig 3). Initial force pressing body to the plane is described by the parameter ν and has the greatest influence on the max ξ when the values of the γ are less then 0. The body displacement is also increasing when the angle between the direction of excitement and the normal of plane β is increasing.

Characteristics of body impact motion have a very strong dependency on initial force pressing body to the plane. The ratio of this force to amplitude of excitement determines if body motion will be with or without lift-offs. When the initial pressing force is increasing, the amplitude of body lift-off and the velocity perpendicular to the sliding plane before the lift-off is decreasing (Fig. 5). If there is no pressing force, the body will not slide on the inclined plane though the velocities before the impact will not be equal to zero. Optimal conditions of body sliding occur for 0.1 ≤ ν ≤ 0.2, where ν is the ratio of initial pressing force to the exciting force.

A c c e p t e d m a n u s c r i p t

Transient regimes of impact motion are well described by impact phase characteristics relating velocities before impact with impact phase. Examples of such characteristics are presented in Figure 6, where dependencies of impact phase on body velocities before impact in ξ direction are depicted.

When the body reaches the position of dynamic equilibrium, the regime of periodical motion about this position begins, i.e. the velocities before impact and the impact phase become constant. For the investigated system, it is characteristic that velocities before impact - ξ & in many cases are negative, i.e. velocity vectors are directed in the opposite direction to the body motion. Explanation of such phenomenon is based on the action of the elastic resisting force. This force is increasing while the body moves forward.

Regimes of body impact motion also are characterized by phase trajectory plots that relate velocity before impact with coordinate of body location on the plane (Fig. 7).

Velocities before impact are decreasing when the body moves away from the position of static equilibrium. Direction of the body velocity vector before impact changes near the position of dynamic equilibrium. Motion parameters of the impact system near the position of dynamic equilibrium have a strong dependency on the parameter γ , which expresses ratio of the body frequency of natural vibration in direction of the motion to the frequency of excitement.

Impact motion of the body on the plane excited in two perpendicular directions

The motion equations ( 12) and (13) describing body displacement above the plane and impact equations ( 9) and ( 11) are also solved by numerical methods. Using plots that show how generalized coordinates ξ and η are changing (Fig. 8), the character of the impact motion is determined. Using dependency ( ) τ η , it was determined that in the period of excitement the body impacts the plane one time.

Analyzing ( ) τ ξ plots, the character of the body motion on the plane is determined.

Usually, the body moves until the generalized coordinate ξ reaches a maximum value, then the body slides slightly backward and forward until the periodic vibrations settle about the position of dynamic equilibrium. Motion from the static equilibrium to the dynamic equilibrium is not a periodic process. During this motion, the coordinate of the impact point is changing.

A c c e p t e d m a n u s c r i p t

The parameter γ , which defines the frequency of the body free vibrations in the x axis direction, has a strong influence on the impact motion regimes. Because this frequency depends on the stiffness coefficient of the elastic constraint, it is possible to maintain that the parameter γ describes the elastic force resisting to body displacement.

Roughly, this force varies inversely with the maximum body displacement max ξ (Fig. 9).

For the impact motion regime and γ < 0.25, the values of max ξ are almost twice as large as for the non-impact regime. The influence of the parameter µ , which is inversely proportional to the amplitude of the plane vibrations along the direction of displacement, on displacement distance is significant in the mentioned γ range. Thus, when this amplitude is decreasing, the body displacement on the vibrating plane is increasing.

It is determined that, when the parameter 6 0. > γ , the body is not moving in one direction but vibrates on the plane about the position of the static equilibrium (Fig. 10).

Then, as in the case of non-impact motion, it is possible to perform automatic searching for the mating surfaces of the assembling components in the local zone. show, that when > λ 0.9, the body is not moving. Upon increasing the coefficient of velocity restitution R, the range of λ values, when the body still can move on the plane in one direction, is increasing.

Parameter σ is inversely proportional to the amplitude of the plane vibrations in the direction of the y axis. Upon increasing the amplitude, the maximum displacement of the body almost remains constant (Fig. 12). When σ >0.5, the body is not moving.

When the phase angle between perpendicular components of the plane vibrations is zero ( 0 = ε ), the body displacement is maximal (Fig. 13). Upon increasing the phase angle ε , the conditions for body displacement worsen. When

2 π ε = or 2 π ε - = and
the axes of the elliptical trajectory of excitement coincide with the axes of the coordinates, the body displacement on the plane begins to increase again.

The average velocity of the body impact motion determines the mating time of the connective surfaces of the components under automated assembly. Therefore, the influence of the system parameters to this velocity should be examined. Influence of the elastic constraint stiffness coefficient λ along the y direction upon the body displacement average velocity ave ξ & in the x direction is not significant (Fig. 15). With a positive phase angle between perpendicular components of vibrations, the average body displacement velocity is higher.

Upon increasing the σ parameter, the amplitude of excitement vibrations along the y axis is decreased and the average velocity of the body displacement ave ξ & remains almost constant (Fig. 16). Upon increasing the plane inclination angle α , the average velocity of the body displacement is also increased.

The dependence of ave ξ & upon ε is complicated but appears to have a maximum.

This maximum is achieved when the phase angle ε between perpendicular components of vibrations is close to zero (Fig. 17). Upon increasing the angle ε , independently of the plane vibrations trajectory direction orientation, ave ξ & is decreased.

Investigations show that body displacement under impact regime is possible only with certain combinations of system parameters. When designing equipment for vibrational automatic assembly, the constructional and excitement parameters should be chosen in a way to ensure proper mating of the surfaces of the assembling components.

Therefore, the zones of the system main parameters are composed where the vibrational impact motion of the body takes place (Figs. 1819).

According to the dimensionless coefficients of the body motion equations, an increase of σ means that the excitement amplitude decreases. Upon increasing the stiffness coefficient of the elastic constraint that is described by parameter λ in the η direction, the body can move for smaller values of σ , i.e. under higher amplitudes in excitement direction η (Fig. 18). The parameters γ and σ are similarly related.

The coefficient of velocity restitution R , which models material properties of the contacting bodies, has a significant influence on the body impact motion zone area (Fig. 19).

For assembling components made from stiffer materials, the value of coefficient R will be higher. As a result, the vibrational motion of the body occurs in a wider range of parameters. If the values of parameters are not directly indicated in the presented plots, they are as follows: = α 0.2; f=0.1; h 1 =h 2 =0.1; R=0.5; = µ 0.5; = ν 0.5; = λ 0.1 .

A c c e p t e d

m a n u s c r i p t

Conclusions

The impact motion of a mobile-based body on an inclined plane characteristic for a vibrational assembly is investigated under kinematical excitement of a body or excitement of a plane in two perpendicular directions. Vibro-impact system is specific because the body impact to the plane is oblique, the coordinate of the point of impact is changing and the elastic and damping forces are resisting to the body displacement.

However, such a system is a typical system applied for the vibrational assembly of components with impact motion of a mobile-based component seeking to align the mating surfaces and automatically compensate the interposition error.

Vibro-impact motion of the body proceeds from the state of static equilibrium to the state of dynamic equilibrium of the system under transient regimes in both cases of excitement. Distance between these equilibrium positions determines maximum error of the mutual interposition of assembling components under which the mating surfaces still can be matched and components can be assembled. Maximum displacement of the body which determines the coordinate of the dynamic equilibrium position mostly depends on force of elastic resistance, phase angle between perpendicular components of vibrations, and amplitude of the plane vibrations in the direction of displacement. If stiffness of the elastic element along the displacement direction is small, then maximum body displacement under vibro-impact body motion regime is almost twice as large as in case of non-impact motion regime.

Average velocity of impact displacement from static until dynamic equilibrium position determines alignment duration of assembling components. This velocity is related to elastic resistance force by a complicated dependency. Average displacement velocity not always is highest under small resistance force. Additionally, it magnitude is determined by amplitude of the vibrational excitement along the displacement direction and phase angle between components of the excitement vibration.

Impact body motion occurs only for certain combinations of system and excitement parameters. Parameter combination regions are constructed where the body can move under the impact regime. The main influences within these combination regions are stiffness of the elastic constraints, excitement amplitude in direction normal to the displacement, and coefficient of velocity restitution.

I

  normal components of impulse are related to each other by the equation -

,

  the jump length of the body moving on the plane is calculated as

3 .

 3 Increasing the angle of the plane inclination α, which is simulating the inclination of the mobile based component under automated assembly with respect to the other component, causes max ξ to increase almost according to a linear law (Fig 4).

  Upon increasing the stiffness coefficient λ of the elastic constraint in the direction of the y axis, values of max ξ have almost no change Fig. (11). Modeling results

Figures 14 -A c c e p t e d m a n u s c r i p t

 14 Figures 14 -17 present dependency of aveξ & on the system main parameters and excitement. The decisive influence on the body motion mode and average velocity is

A c c e p t e d m a n u s c r i p t A c c e p t e d m a n u s c r i p t List of Figures 1 . 8 .

 18 Model of the impact motion when body is excited kinematically or the plane is excited in two perpendicular directions. 2. Plots of body displacement: (a) -Examples of vibration of generalized body coordinates: (a) -plots of ( ) b) -plots of ( ) τ ξ : (1) -γ=0.1 µ=0.5 ν=0.5 α=0.2; (2) -γ=0.2 µ=0.5 ν=0.5 α=0.3; (3) -γ=0.2 µ=2.0 ν=0.5 α=0.2. of the parameters combination in λ-σ coordinates under vibrational displacement (hatched area). 19. Region of the parameters combination in γ-σ coordinates under vibrational displacement (hatched area).

  Fig. 1

A c c e p t e d m a n u s c r i p tFig

  Fig.2(a)

A c c e p t e d m a n u s c r i p t Fig. 4

 4 Fig.4

A c c e p t e d m a n u s c r i p t Fig. 5

 5 Fig.5

A c c e p t e d m a n u s c r i p t Fig. 6

 6 Fig. 6

A c c e p t e d m a n u s c r i p t Fig. 7

 7 Fig. 7

  Fig.8(a)

  Fig.11

  Fig.12

  Fig.13

  Fig.14

  Fig.15

  Fig.16

  Fig.17

  Fig.18

  Fig.19

& & & ,(12)

& & & &

& & & & 0 2 (17)here c ξ is the coordinate of the impact point; L is the length of the body jump.Integration of equations (14) and (15) yields:


Z K Peng
email: z.peng@sheffield.ac.uk

Z Q Lang
email: z.lang@sheffield.ac.uk

Detecting the position of nonlinear component in periodic structures from the system responses to dual sinusoidal excitations

Based on the Nonlinear Output Frequency Response Functions (NOFRFs), a novel method is developed to detect the position of nonlinear components in periodic structures. The detection procedure requires exciting the nonlinear systems twice using two sinusoidal inputs separately. The frequencies of the two inputs are different; one frequency is twice as high as the other one. The validity of this method is demonstrated by numerical studies. Since the position of a nonlinear component often corresponds to the location of defect in periodic structures, this new method is of great practical significance in fault diagnosis for mechanical and structural systems.

Periodic structures are defined as structures consisting of identical substructures connected to each other in identical manner. The real life systems which can be modelled as either finite or infinite, one-dimension or multi-dimension periodic structures range from the simple structures like periodically supported beams [START_REF] Duhame | Finite element analysis of the vibrations of waveguides and periodic structures[END_REF]~ [START_REF] Yuan | On the determination of phase constants for the study of the free vibration of periodic structures[END_REF], plates [START_REF] Mukherjee | Normal modes and acoustic response of finite one-and two-dimensional multi-bay periodic panels[END_REF][6] and building block [START_REF] Fassois | Stochastic output error vibration-based damage detection and assessment in structures under earthquake excitation[END_REF]. Analysis of the free and forced vibration and the mode analysis of linear periodic structures are of particular interests [START_REF] Duhame | Finite element analysis of the vibrations of waveguides and periodic structures[END_REF]- [START_REF] Yuan | On the determination of phase constants for the study of the free vibration of periodic structures[END_REF]. Mead [START_REF] Mead | Wave propagation in continuous periodic structures: research contributions from Southampton, 1964-1995[END_REF] provides an excellent review about periodic structure studies.

Attentions also have been paid to the study of nonlinear periodic structures [START_REF] Chi | On damped Non-linear Dynamic systems with many degrees of freedom[END_REF]- [START_REF] Marathe | Wave attenuation in nonlinear periodic structures using harmonic balance and multiple scales[END_REF]. Chi and Rosenberg [START_REF] Chi | On damped Non-linear Dynamic systems with many degrees of freedom[END_REF] have studied the existence of classical normal mode motion in onedimension non-linear mass-spring-damper systems with many degrees of freedom where all springs and / or all dampers may be strongly non-linear. Using exact and asymptotic techniques, Vakakis et al [START_REF] Vakakis | Forced Localization in a periodic chain of non-linear oscillators[END_REF] have studied the localized responses of a nonlinear periodic oscillator chain subjected to harmonic excitations with general spatial distributions. Royston and Singh [START_REF] Royston | Periodic response of mechanical systems with local nonlinearities using an enhanced Galerkin technique[END_REF] have studied the periodic response of mechanical systems with local nonlinearity using an enhanced Galerkin technique where a new semi-

A c c e p t e d m a n u s c r i p t

2 analytical framework for the study of mechanical systems with local nonlinearities has been presented. Chakraborty and Mallik [START_REF] Chakraborty | Dynamics of a weakly non-linear periodic chain[END_REF] have investigated the harmonic vibration propagation in an infinite, non-linear periodic chain using a perturbation approach where the case of cyclic one-dimension nonlinear chain has also been discussed. Marathe and Chatterjee [START_REF] Marathe | Wave attenuation in nonlinear periodic structures using harmonic balance and multiple scales[END_REF] have studied the Wave attenuation in onedimension nonlinear periodic structures using harmonic balance and multiple scale method. In engineering practices, there are considerable periodic structures that behave nonlinearly just because one or a few components have nonlinear properties, and the nonlinear component is often the component where a fault or abnormal condition occurs.

One of the well known examples is beam structures [START_REF] Neild | A discrete model of vibrating beam using time-stepping approach[END_REF] with breathing cracks, the global nonlinear behaviors of which are caused only by a few cracked elements. Therefore it is of great significance to effectively detect the position of the nonlinear component in a periodic structure. The detection of damage in large periodic structures has been studied by Zhu and Wu [START_REF] Zhu | The Characteristic receptance method for damage detection in large Mono-coupled Periodic structures[END_REF]. In their studies, the periodic structure with damage is still considered to be linear and the locations and magnitude of damage in large monocoupled periodic systems have been estimated using measured changes in the natural frequencies. Based on a one-dimensional periodic structure model, Sakellariou and Fassois [START_REF] Fassois | Stochastic output error vibration-based damage detection and assessment in structures under earthquake excitation[END_REF][16] have used a stochastic output error vibration-based methodology to detect the damage in structures where the damage elements were modeled as components of cubic stiffness.

The Volterra series approach [START_REF] Worden | A harmonic probing algorithm for the multi-input Volterra series[END_REF] is a powerful method for the analysis of nonlinear systems, and extends the familiar concept of the convolution integral for linear systems to a series of multi-dimensional convolution integrals. The Fourier transforms of the Volterra kernels are known as the kernel transforms, Higher-order Frequency Response Functions (HFRFs) [START_REF] Lang | Output frequency characteristics of nonlinear system[END_REF], or more usually Generalised Frequency Response Functions (GFRFs). These provide a convenient concept for analyzing nonlinear systems in the frequency domain. If a differential equation or discrete-time model is available for a nonlinear system, the GFRFs can be determined using the algorithm in [START_REF] Billings | Spectral analysis for nonlinear system, part I: parametric non-linear spectral analysis[END_REF]~ [START_REF] Peyton Jones | A recursive algorithm for the computing the frequency response of a class of nonlinear difference equation models[END_REF]. The GFRFs can be regarded as the extension of the classical frequency response function (FRF) for linear systems to the nonlinear case. The concept of Nonlinear Output Frequency Response Functions (NOFRFs) [START_REF] Lang | Energy transfer properties of nonlinear systems in the frequency domain[END_REF] is an alternative extension of the FRF to the nonlinear case. NOFRFs are one dimensional functions of frequency, which allow the analysis of nonlinear systems to be implemented in a manner similar to the analysis of linear systems, and which provides great insight into the mechanisms which dominate many important nonlinear behaviors. The paper is organized as follows. Section 2 gives a brief introduction to the new concept of NOFRFs. Some important properties of the NOFRFs for locally nonlinear MDOF systems, which were first revealed in the authors' recent studies [START_REF] Peng | Nonlinear Output Frequency Response Functions of MDOF Systems with Multiple Nonlinear Components[END_REF], are introduced in Section 3. The novel method for the nonlinear component position detection is presented in Section 4. In Section 5, three numerical case studies are used to verify the effectiveness of the proposed method. Finally conclusions are given in Section 6.

Output Frequency Response Functions of Nonlinear Systems

Output Frequency Response Functions under General Input

The definition of NOFRFs is based on the Volterra series theory of nonlinear systems. The Volterra series extends the well-known convolution integral description for linear systems to a series of multi-dimensional convolution integrals, which can be used to represent a wide class of nonlinear systems [START_REF] Lang | Output frequency characteristics of nonlinear system[END_REF].

Consider the class of nonlinear systems which are stable at zero equilibrium and which can be described in the neighbourhood of the equilibrium by the Volterra series

i n i i n N n n d t u h t x τ τ τ τ ) ( ) ,..., ( ) ( 1 1 1 ∏ ∑ ∫ ∫ = = ∞ ∞ - ∞ ∞ - - = (1)
where x(t) and u(t) are the output and input of the system, ) ,..., (

n 1 n h τ τ
is the nth order Volterra kernel, and N denotes the order of the Volterra series representation. Lang and Billings [START_REF] Lang | Output frequency characteristics of nonlinear system[END_REF] derived an expression for the output frequency response of this class of nonlinear systems to a general input. The result is

⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ = ∀ = ∫ ∏ ∑ = + + = - = ω ω ω ω σ ω ω ω π ω ω ω ω n n n i i n n n n N n n d j U j j H n j X j X j X ,..., 1 1 1 1 1 ) ( ) ,..., ( ) 2 ( 1 ) ( for ) ( ) ( (2)

A c c e p t e d m a n u s c r i p t 4

This expression reveals how nonlinear mechanisms operate on the input spectra to produce the system output frequency response. In (2), )

( ω j X
is the spectrum of the system output, )

( ω j X n represents the nth order output frequency response of the system,

+ + - ∞ ∞ - ∞ ∞ - ∫ ∫ = (3)
is the nth order Generalised Frequency Response Function (GFRF) [START_REF] Lang | Output frequency characteristics of nonlinear system[END_REF], and . Equation ( 2) is a natural extension of the well-known linear relationship

∫ ∏ = + + = ω ω ω ω σ ω ω ω n n n i i n n d j U j j H ,..., 1 1 1 )
) ( ) ( ) ( ω ω ω j U j H j X = , where ) ( ω j H
is the frequency response function, to the nonlinear case.

For linear systems, the possible output frequencies are the same as the frequencies in the input. For nonlinear systems described by equation (1), however, the relationship between the input and output frequencies is more complicated. Given the frequency range of an input, the output frequencies of system (1) can be determined using the explicit expression derived by Lang and Billings in [START_REF] Lang | Output frequency characteristics of nonlinear system[END_REF] [START_REF] Lang | Output frequencies of nonlinear systems[END_REF].

Based on the above results for the output frequency response of nonlinear systems, a new concept known as the Nonlinear Output Frequency Response Function (NOFRF) was recently introduced by Lang and Billings [START_REF] Lang | Energy transfer properties of nonlinear systems in the frequency domain[END_REF]. The NOFRF is defined as

∫ ∏ ∫ ∏ = + + = = + + = = ω ω ω ω ω ω ω ω σ ω σ ω ω ω ω n n n n i i n n i i n n n d j U d j U j j H j G ,..., 1 ,..., 1
1 1 1 ) ( ) ( ) ,..., ( ) ( (4)
under the condition that 0

) ( ) ( ,..., 1 1 ≠ = ∫ ∏ = + + = ω ω ω ω σ ω ω n n n i i n d j U j U (5) Notice that ) ( ω j G n
is valid over the frequency range of )

( ω j U n
, which can be determined using the algorithm in [START_REF] Lang | Output frequency characteristics of nonlinear system[END_REF].

By introducing the NOFRFs

)

( ω j G n , N n , 1 =
, equation (2) can be written as

∑ ∑ = = = = N n n n N n n j U j G j X j X 1 1 ) ( ) ( ) ( ) ( ω ω ω ω (6)
which is similar to the description of the output frequency response for linear systems. The NOFRFs reflect a combined contribution of the system and the input to the system is that it is one dimensional, and thus allows the analysis of nonlinear systems to be implemented in a convenient manner similar to the analysis of linear systems. Moreover, there is an effective algorithm [START_REF] Lang | Energy transfer properties of nonlinear systems in the frequency domain[END_REF] available which allows the estimation of the NOFRFs to be implemented directly using system input output data.

Output Frequency Response Functions under Harmonic Inputs

When system (1) is subject to a harmonic input

) cos( ) ( β ω + = t A t u F (7)
Lang and Billings [START_REF] Lang | Output frequency characteristics of nonlinear system[END_REF] showed that equation (2) can be expressed as

∑ ∑ ∑ = = + + = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = = N n k k k k n n N n n kn k n n j A j A j j H j X j X 1 1 1 1 1 ) ( ) ( ) , , ( 2 1 )
( ) ( ω ω ω ω ω ω ω ω ω (8)
where

⎩ ⎨ ⎧ = 0 | | ) ( ) ( sign β ω k j k e A j A i if { } otherwise , , 1 , 1 , n i k k F k i = ± = ∈ ω ω (9)
Define the frequency components of nth order output of the system as n Ω . Then according to equation ( 8), the frequency components in the system output can be expressed as

∪ N n n 1 = Ω = Ω (10)
and n Ω is determined by the set of frequencies

{ } n i F k k k i n , , 1 , | 1 = ± = + + = ω ω ω ω ω (11) From equation (11), it is known that if all n k k ω ω , , 1 are taken as F ω -, then F nω ω - = . If k of these are taken as F ω , then F k n ω ω ) 2 ( + - =
. The maximal k is n. Therefore the possible frequency components of )

( ω j X n are n Ω ={ } n k k n F , , 1 , 0 , ) 2 ( = + - ω (12) Moreover, it is easy to deduce that } , , 1 , 0 , 1 , , , { 1 N N k k F N n n ∪ - - = = Ω = Ω = ω (13)

A c c e p t e d m a n u s c r i p t

6 Equation ( 13) explains why superharmonic components are generated when a nonlinear system is subjected to a harmonic excitation. In the following, only those components with positive frequencies will be considered.

The NOFRFs defined in equation ( 4) can be extended to the case of harmonic inputs [START_REF] Peng | Resonances and resonant frequencies for a class of nonlinear systems[END_REF] as

∑ ∑ = + + = + + = ω ω ω ω ω ω ω ω ω ω ω ω ω kn k n kn k n n k k n k k k k n n H n j A j A j A j A j j H j G 1 1 1 1 1 ) ( ) ( 2 1 ) ( ) ( ) , , ( 2 1 )
( (n = 1,…, N ) (14) under the condition that 0 ) ( ) ( 2 1 ) ( 1 1 ≠ = ∑ = + + ω ω ω ω ω ω kn k n k k n n j A j A j A ( 15
)
Obviously, )

( ω j G H n is only valid over n
Ω defined by equation [START_REF] Chakraborty | Dynamics of a weakly non-linear periodic chain[END_REF]. Consequently, the output spectrum )

( ω j X of nonlinear systems under a harmonic input can be expressed as

∑ ∑ = = = = N n n H n N n n j A j G j X j X 1 1 ) ( ) ( ) ( ) ( ω ω ω ω (16) When k of the n frequencies of n k k ω ω , ,
k n j k n n n F n e C A k n j A + - = + - (17)
Moreover

β β ω ω ω ω ω ) 2 ( ) 2 ( | | 2 1 | | ) , , , , , ( 2 1 )
) 2 ( ( k n j k n n n k n j k n n k n F F k F F n n F H n e C A e C A j j j j H k n j G + - + - - - - = + - ) , , , , , ( k n F F k F F n j j j j H - - - = ω ω ω ω (18) where ) ,..., ( 1 n n j j H ω ω
is assumed to be a symmetric function. Therefore, in this case, )

( ω j G H n over the nth order output frequency range n Ω = { } n k k n F , , 1 , 0 , ) 2 ( = + - ω is equal to the GFRF ) ,..., ( 1 n n j j H ω ω evaluated at , 1 F k ω ω ω = = = F n k ω ω ω - = = = + 1 , n k , , 0 = .
Substituting equations ( 17) and ( 18) into ( 16), it can be derived that the output spectrum )

( ω j Y of nonlinear systems subjected to a harmonic input can be expressed as

[ ] ∑ + - = - + - + = 2 / ) 1 ( 1 ) 1 ( 2 ) 1 ( 2 ) ( ) ( ) ( k N n F n k F H n k F jk A jk G jk X ω ω ω ( N k , , 1 , 0 = ) (19)
where [•] denotes the operator of taking the integral.

NOFRFs of Nonlinear Periodic Structures

Consider the one-dimension nonlinear periodic structures where the Lth component is nonlinear, which have be used in [START_REF] Fassois | Stochastic output error vibration-based damage detection and assessment in structures under earthquake excitation[END_REF], [START_REF] Chakraborty | Dynamics of a weakly non-linear periodic chain[END_REF], [START_REF] Marathe | Wave attenuation in nonlinear periodic structures using harmonic balance and multiple scales[END_REF] and [START_REF] Zhu | The Characteristic receptance method for damage detection in large Mono-coupled Periodic structures[END_REF], shown in Figure 1.

∑ = Δ = Δ P i i i LD w S 1 ) ( ( 20
)
where P is the degree of the polynomial. Without loss of generality, further assume

n L , 1 ≠ . Denote ∑ ∑ = - = - - + - = P i i L L i P i i L L i x x r x x w NonF 2 1 2 1 ) ( ) ( (21)
' 2 0 0 0 0 ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - = - - L n L NonF NonF NF ( 22
)
Then the motion of the nonlinear oscillator in Figure 1 is described in a matrix form as.

) (t F NF Kx x C x M + - = + + ( 23
)
where M is the system mass matrix,

⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ = n m m m M 0 0 0 0 0 0 2 1 and ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ - - + - - + - - + = - - n n n n n n c c c c c c c c c c c c c C 0 0 0 0 0 0 1 1 3 3 2 2 2 2 1 ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ - - + - - + - - + = - - n n n n n n k k k k k k k k k k k k k K 0 0 0 0 0 0 1 1 3 3 2 2 2 2 1
are the system mass, damping and stiffness matrix respectively. ' 1 ) , , (

n x x x =
is the displacement vector, and

m n k n m n-1 m L-1 m 1 k L-1 k 1 x 1 x L-1 x n-1 x n u(t) m L+1 k L+1 m L k L x L x L+1
A c c e p t e d m a n u s c r i p t

8 ' 1 ) 0 , , 0 ), ( , 0 , , 0 ( ) ( J n J t u t F - - =
(24) is the external force vector acting on the J th mass of the oscillator.

The system described by equation ( 23) is a typical locally nonlinear periodic structure. The Lth nonlinear component can lead the whole system to behave nonlinearly. In this case, the Volterra series can be used to describe the relationships between the displacements ) (t

x i ( n i , , 1 =
) and the input force

) (t u as below Z j Z Z N j j j i i d t u h t x τ τ τ τ ) ( ) ,..., ( ) ( 1 1 1 ) , ( ∏ ∑ ∫ ∫ = = ∞ ∞ - ∞ ∞ - - = (25)
where ) ,..., ( 1) , (

j j i h τ τ
is the jth order Volterra kernel associated to the ith mass. In the frequency domain, the relationship ( 16) can be expressed as

∑ ∑ = = = = N l l l i N l l i i j U j G j X j X 1 ) , ( 1
) , ( ) ( ) ( ) ( ) ( ω ω ω ω ( n i , , 1 = ) (26) where ) ( ) , ( ω j G l i
is the lth order NOFRF associated to the ith mass.

Without loss of generality, assume L<J, as revealed in [START_REF] Peng | Nonlinear Output Frequency Response Functions of MDOF Systems with Multiple Nonlinear Components[END_REF], for any two consecutive masses, the NOFRFs of system (23) satisfy following relationships

) ( ) ( ) ( ) ( ) ( ) ( ) , 1 ( ) , ( 1 , ) 2 , 1 ( ) 2 , ( 1 , 2 ω
ω ω λ ω ω ω λ j G j G j j G j G j N i N i i i N i i i i + + + + = = = = ( 1 1 - ≤ ≤ n i ) (27) ) ( ) ( ) ( ) ( ) ( ) ( 1 , ) , 1 ( ) , ( ) 1 , 1 ( ) 1 , ( 1 , 1 ω λ ω ω ω ω ω λ j j G j G j G j G j i i Z Z i Z i i i i i + + + + = = = ( 2 1 - ≤ ≤ L i or 1 - ≤ ≤ n i J ; N Z ≤ ≤ 2 ) (28) ) ( ) ( ) ( ) ( ) ( ) ( 1 , ) , 1 ( ) , ( ) 1 , 1 ( ) 1 , ( 1 , 1
ω λ ω ω ω ω ω λ j j G j G j G j G j i i Z Z i Z i i i i i + + + + = ≠ = ( 1 1 - ≤ ≤ - J i L , N Z ≤ ≤ 2 ) (29)
Based on these relationships of the NOFRFs, a novel method can be developed to determine the position of the nonlinear element in system [START_REF] Peng | Nonlinear Output Frequency Response Functions of MDOF Systems with Multiple Nonlinear Components[END_REF].

As equations ( 27)~(29) are the basis of this study, but the derivation and justification of these equations can take quite a lot of contents and, therefore, the elaboration of them are not given in this paper. However, for the sake of completion, a numerical case is given here to justify equations (27)~(29). The numerical case study is conducted on a damped 8-DOF oscillator whose fourth spring (L = 4) was nonlinear. As widely used in modal analysis, the damping was assumed to be proportional to the stiffness, e.g., K C μ = . The values of the system parameters are taken as 1 16) and ( 17), the frequency components of the outputs of the 8 masses can be written as 30), it can be seen that, using the method in [START_REF] Lang | Energy transfer properties of nonlinear systems in the frequency domain[END_REF], two different inputs with the same waveform but different strengths are sufficient to estimate the NOFRFs up to 4P th P order. Therefore, in this numerical study, two different inputs are used with A=0.8 and A=1.0 respectively. The simulation studies were conducted using a fourth-order Runge-Kutta method to obtain the forced response of the system.

8 1 = = = m m , , 10 5531 . 3
× = = = = k k r μ =0.01
) ( ) ( ) ( ) ( ) ( 3 ) 3 , ( 1 ) 1 , ( F F H i F F H i F i j U j G j U j G j X ω ω ω ω ω + = ) 2 ( ) 2 ( ) 2 ( ) 2 ( ) 2 ( 4 ) 4 , ( 2 ) 2 , ( F F H i F F H i F i j U j G j U j G j X ω ω ω ω ω + = ) 3 ( ) 3 ( ) 3 ( 3 ) 3 , ( F F H i F i j U j G j X ω ω ω = ) 4 ( ) 4 ( ) 4 ( 4 ) 4 , ( F F H i F i j U j G j X ω ω ω = ) 8 , , 1 ( = i (30) From equation (
The evaluated results of ) (

1 F H j G ω , ) ( 3 F H j G ω , ) 2 ( 2 F H j G ω and ) 2 ( 4 F

H j G

ω for all masses are given in Table 1. According to relationships (27)~(29), it is known that the following relationships should be tenable. ) From the NOFRFs in Table 1,

i i F H i F H i F H i F H i F i i j j G j G j G j G j ω λ ω ω ω ω ω λ + + + + = = = for 7 , 6 , 2 , 1 = i ) ( ) ( ) ( ) ( ) ( ) ( 1 , 3 ) 3 , 1 ( ) 3 , ( ) 1 , 1 ( ) 1 , ( 1 , 1 F i i F H i F H i F H i F H i F i i j j G j G j G j G j ω λ ω ω ω ω ω λ + + + + = ≠ = for 5 , 4 , 3 = i ) 2 ( ) 2 ( ) 2 ( ) 2 ( ) 2 ( ) 2 ( 1 , 4 ) 4 , 1 ( ) 4 , ( ) 2 , 1 ( ) 2 , ( 1 , 2 F i i F H i F H i F H i F H i F i i j j G j G j G j G j ω λ ω ω ω ω ω λ + + + + = = = for 7 , , 1 = i (31) Table 1. The evaluated results of ) ( 1 F H j G ω , ) ( 3 F H j G ω ) 2 ( 2 F H j G ω and ) 2 ( 4 F H j G ω ) ( 1 F H j G ω (×10P -6 P) ) ( 3 F H j G ω (×10P -9 P) ) 2 ( 2 F H j G ω (×10P -9 P) ) 2 ( 4 F H j G ω (×10P -10 P) Mass 1 -1.9442+2.8776i ( ) ( ) ( ) ( ) ( ) ( 1 , 3 ) 3 , 1 ( ) 3 , ( ) 1 , 1 ( ) 1 , ( 1 , 1 F
F i i jω λ + , ) ( 1 , 3 F i i jω λ + ) 2 ( 1 , 2 F i i j ω λ + and ) 2 ( 1 , 4 F i i j ω λ + ) ( 1 , 1 F i i jω λ + ) ( 1 , 3 F i i jω λ + ) 2 ( 1 , 2 F i i j ω λ + ) 2 ( 1 , 4 F i i j ω λ + i=1 0.
) ( 1 , 1 F i i jω λ + , ) ( 1 , 3 F i i jω λ + , ) 2 ( 1 , 2 F i i j ω λ + and ) 2 ( 1 , 4 F i i j ω λ + ( 7 , , 1 = i
) can be evaluated, and the results are given in Tables 2. The results shown in Tables 2 have a strict accordance with the relationships in (31). Therefore, the numerical study verifies relationships (27)~(29).

The Nonlinear Component Position Detection Method

) ( ) ( ) ( ) ( ) ( ) ( ) ( 1 , ) 1 2 , 1 ( ) 1 2 , ( ) 3 , 1 ( ) 3 , ( ) 1 , 1 ( ) 1 , ( F i i F D i F D i F i F i F i F i j j G j G j G j G j G j G ω λ ω ω ω ω ω ω + + + + + + = = = = = ( , 2 , 1 = D ) ( 2 1 - ≤ ≤ L i or 1 - ≤ ≤ n i J ) (32) ) ( ) ( ) ( ) ( ) ( ) ( ) ( 1 , ) 1 2 , 1 ( ) 1 2 , ( ) 3 , 1 ( ) 3 , ( ) 1 , 1 ( ) 1 ,
( F i i F D i F D i F i F i F i F i j j G j G j G j G j G j G ω λ ω ω ω ω ω ω + + + + + + = = = = ≠ ( , 2 , 1 = D ) ( 1 1 - ≤ ≤ - J i L ) (33) ) 2 ( ) 2 ( ) 2 ( ) 2 ( ) 2 ( ) 2 ( ) 2 ( 1 , ) 2 , 1 ( ) 2 , ( ) 4 , 1 ( ) 4 , ( ) 2 , 1 ( ) 2 , ( F i i D i D i F i F i F i F i j j G j G j G j G j G j G ω λ ω ω ω ω ω ω + + + + = = = = = ) 3 ( ) 3 ( ) 3 ( ) 3 ( ) 3 ( ) 3 ( ) 3 ( 1 , ) 1 2 , 1 ( ) 1 2 , ( ) 5 , 1 ( ) 5 , ( ) 3 , 1 ( ) 3 , ( F i i F D i F D i F i F i F i F i j j G j G j G j G j G j G ω λ ω ω ω ω ω ω + + + + + + = = = = = ) , 2 , 1 ( = D ) 1 , , 1 ( - = n i
(34) According to equation [START_REF] Billings | Spectral analysis for nonlinear system, part I: parametric non-linear spectral analysis[END_REF], the first harmonic components of

) (t x i ( n i , , 1
=
) can be written as

[ ] ∑ = - - = 2 / 1 1 2 ) 1 2 , ( ) ( ) ( ) ( N k F k F H k i F i j A j G j X ω ω ω ( n i , , 1
) ( ) ( ) ( ) ( ) ( ) ( 1 1 , 2 / 1 1 2 ) 1 2 , 1 ( 1 , F i F i i N k F k F H k i F i i F i j X j j A j G j j X ω ω λ ω ω ω λ ω + + = - - + + = = ∑ ( 2 1 - ≤ ≤ L i or 1 - ≤ ≤ n i J ) (36) Therefore, ) ( ) ( ) ( 1 1 , F i F i F i i j X j X j ω ω ω λ + + = ( 2 1 - ≤ ≤ L i or 1 - ≤ ≤ n i J ) (37)
For the masses located between the nonlinear spring and the input force, substituting equation ( 33) into (35) yields,

[ ] [ ] ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( 1 1 , 2 / 1 1 2 ) 1 2 , 1 ( 1 , 1 ) 1 , 1 ( 1 , 2 / 1 1 2 ) 1 2 , 1 ( 1 , 1 ) 1 , 1 ( 1 , 1 F i F i i N k F k F H k i F i i F F H i F i i N k F k F H k i F i i F F H i F i i F i j X j j A j G j j A j G j j A j G j j A j G j j X ω ω λ ω ω ω λ ω ω ω λ ω ω ω λ ω ω ω λ ω + + = + + + + + + = + + + + + + = + ≠ + = ∑ ∑ ( 1 1 - ≤ ≤ - J i L ) (38) Obviously, ) ( ) ( ) ( 1 1
,

F i F i F i i j X j X j ω ω ω λ + + ≠ ( 1 1 - ≤ ≤ - J i L ) (39)
According to equation [START_REF] Billings | Spectral analysis for nonlinear system, part I: parametric non-linear spectral analysis[END_REF], the second harmonic components of

) (t x i ( n i , , 1
=
) can be written as

[ ] ∑ - = = 2 / ) 1 ( 1 2 ) 2 , ( ) 2 ( ) 2 ( ) 2 ( N k F k F H k i F i j A j G j X ω ω ω ( n i , , 1 = ) (40)
Substituting equation ( 34) into (40) yields

[ ] [ ] ) 2 ( ) 2 ( ) 2 ( ) 2 ( ) 2 ( ) 2 ( ) 2 ( ) 2 ( 1 1 , 2 / ) 1 ( 1 2 ) 2 , 1 ( 1 , 2 / ) 1 ( 1 2 ) 2 , ( F i F i i N k F k F H k i F i i N k F k F H k i F i j X j j A j G j j A j G j X ω ω λ ω ω ω λ ω ω ω + + - = + + - = = = = ∑ ∑ ( 1 , , 1 - = n i ) (41) Consequently ) 2 ( ) 2 ( ) 2 ( 1 1 , F i F i F i i j X j X j ω ω ω λ + + = ) 1 , , 1 ( - = n i (42)
Similarly, it can be deduced that

) ( ) ( ) ( 1 1 , F i F i F i i jD X jD X jD ω ω ω λ + + = ( 2 ≥ D , 1 , , 1 - = n i ) (43)
When using two sinusoidal inputs, the frequencies of which are 1 F ω and 2 F ω , to excite the nonlinear system under study, it is easy to extract the harmonic components from the

A c c e p t e d m a n u s c r i p t

12 FFT spectra of the system responses, denote them as

) ( 1 ) , 1 ( F i F jD X ω and ) ( 2 ) , 2 ( F i F jD X ω ( n i , , 1 = ; , 2 , 1 = D
) respectively. It can be known from (37), ( 39) and (43) that

) ( ) ( ) ( 2 ) 1 , 2 ( 2 ) , 2 ( 2 1 ,
F i F F i F F i i j X j X j ω ω ω λ + + = ( 2 1 - ≤ ≤ L i or 1 - ≤ ≤ n i J ) (44) ) ( ) ( ) ( 2 ) 1 , 2 ( 2 ) , 2 ( 2 1 ,
F i F F i F F i i j X j X j ω ω ω λ + + ≠ ( 1 1 - ≤ ≤ - J i L ) (45) ) ( ) ( ) ( 1 ) 1 , 1 ( 1 ) , 1 ( 1 1 ,
F i F F i F F i i jD X jD X jD ω ω ω λ + + = ) 1 , , 1 ( - = n i (46) If P F F 2 1 ω ω =
and P is an integral, then for the P th harmonic components of

) , 1 ( i F X ) , , 1 ( n i =
, obviously, equation ( 46) can be rewritten as

) ( ) ( ) ( ) ( 2 1 , 1 ) 1 , 1 ( 1 ) , 1 ( 1 1
, F i i F i F F i F F i i j jP X jP X jP ω λ ω ω ω λ + + + = = ) 1 , , 1 ( - = n i (47)
From ( 44), ( 45) and (47) it can be known that ) (

F i F F i F F i F F i F j X j X jP X jP X ω ω ω ω + + = ( 2 1 - ≤ ≤ L i or 1 - ≤ ≤ n i J ) ) ( ) ( ) ( 2 ) 1 , 2 ( 2 ) , 2 ( 1 ) 1 , 1 ( 1 ) , 1 (
F i F F i F F i F F i F j X j X jP X jP X ω ω ω ω + + ≠ ( 1 1 - ≤ ≤ - J i L ) (48) and ) ( ) ( ) ( ) ( 2 ) 1 , 2 ( 2 ) , 2 ( 1 ) 1 , 1 ( 1 ) , 1 (
The relationships given in (48) and (49) provide a simple way to detect the position of nonlinear components in the MDOF nonlinear systems using dual sinusoidal excitations. Usually, for simplicity and convenience, 1 F ω can be chosen as

F F ω ω = .
2) Calculate the FFT spectra of the two sets of system responses, denote them as

) , 1 ( i F X and ) , 2 ( i F X , ) , , 1 ( n i = . 3) Extract the first harmonics ) ( 2 ) , 2 ( F i F j X ω from ) , 2 ( i F X
, and the second harmonics

) 2 ( 1 ) , 1 ( F i F j X ω from ) , 1 ( i F X , ) , , 1 ( n i = . 4) Calculate ) ( / ) ( 2 ) 1 , 2 ( 2 ) , 2 ( F i F F i F j X j X ω ω + and ) 2 ( / ) 2 ( 1 ) 1 , 1 ( 1 ) , 1 ( F i F F i F j X j X ω ω + , denote them as ) ( 2 1 , 2 F i i F j R ω + and ) 2 ( 1 1 , 1 F i i F j R ω + , ) 1 , , 1 ( - = n i . 5) Find out the masses where ≠ + ) ( 2 1 , 2 F i i F j R ω ) 2 ( 1 1 , 1 F i i F j R ω +
, then the component on the right side of the furthest left mass is the nonlinear one.

It is worth noting here that, if the force position is located on the left side of the nonlinear component, that is, L J ≤ , then, at the 5 th step of the above detection procedure, the nonlinear component is the one on the right side of the furthest right mass which satisfies

+ ) ( 2 1 , 2 F i i F j R ω ) 2 ( 1 1 , 1 F i i F j R ω +
. In addition, if there are no nonlinear components in the system, then there is no superharmonic component in the system response spectrum. Therefore there is no need to use the above detection procedure.

The novel nonlinear component position detection method requires only two tests where the MDOF system excited by two different sinusoidal forces. This is obviously very easy to carry out in practices. In the following section, the effectiveness of this method will be demonstrated using numerical studies.

Numerical Studies

In order to verify the nonlinear component position detection method, the same oscillator used in Section 3 was adopted.

Case Study 1 (L < J): π respectively, and were imposed on the 6 th mass of this system, that is J = 6. The responses of the system were obtained using a fourth-order Runge-Kutta method to integrate equation (31). The second super-harmonics were used to calculate ) 2 ( The results given in Table 6 clearly show that

Table 3, ) 2 ( 1 ) , 1 ( F i F j X ω and ) 2 ( 1 1 , 1 F i i F j R ω + (i=1,…,8) ) 2 ( 1 ) , 1 ( F i F j X ω ( 8 10 - × ) ) 2 ( 1 1 , 1 F i i F j R ω + i=1 -0.
, 1 F i i F j R ω + and ) ( 2 1 , 2 F i i F j R ω + (i=1,…,7) ) 2 ( 1 1 , 1 F i i F j R ω + ) ( 2 1 , 2 F i i F j R ω + i=1 0. 2 ( 1 1
) 2 ( 1 1 , 1 F i i F j R ω + ≠ ) ( 2 1 , 2 F i i F j R ω + at i = 3,4,5.
According to the detection method in Section 4, it can be known that the component on the right side of the 3 rd mass is the nonlinear one, that is, the 4 th spring component.

Case Study 2 (L = J):

In this case study, the 4 th spring are nonlinear, that is L=4. The sinusoidal forces used in this case are the same as the ones used in Case 1, but were imposed on the 4 th mass of this system, that is J=4, so L=J. The results of

) 2 ( 1 ) , 1 ( F i F j X ω and ) 2 ( 1 1 , 1 F i i F j R ω + (i=1,…,8) are given in
, ) 2 ( 1 ) , 1 ( F i F j X ω and ) 2 ( 1 1 , 1 F i i F j R ω + (i=1,…,8) ) 2 ( 1 ) , 1 ( F i F j X ω ( 8 10 - × ) ) 2 ( 1 1 , 1 F i i F j R ω + i=1
) , 2 ( F i F j X ω and ) ( 2 1 , 2 F i i F j R ω + (i=1,…,8) ) ( 2 ) , 2 ( F i F j X ω ( 5 10 - × ) ) ( 2 1 , 2 F i i F j R ω + i=1 0. ( 2
F i i F j R ω + and ) ( 2 1 , 2 F i i F j R ω + (i=1,…,7) ) 2 ( 1 1 , 1 F i i F j R ω + ) ( 2 1 , 2 F i i F j R ω + i=1
i i F j R ω + ≠ ) ( 2 1 , 2 F i i F j R ω + only at i = 3.
According to the detection method in Section 4, it can be known that the component on the right side of the 3 rd mass is the nonlinear one, that is, the 4 th spring component.

Case Study 3 (L>J):

In this case study, the 7 th spring are nonlinear, that is L=7. The sinusoidal forces used in this case are the same as the ones used above cases, but were imposed at the 4 th mass of this system, that is J=4, so L>J. The results of According to the detection method in Section 4, it can be known that the component on the right side of the 6 th mass is the nonlinear one, that is, the 7 th spring component.

) 2 ( 1 ) , 1 ( F i F j X ω , ) 2 ( 1 1 , 1 F i i F j R ω + , ) ( 2 )

Conclusions and Remarks

Based on the properties of NOFRFs, a novel method is developed to detect the position of the nonlinear component in a periodic structure. The detection procedure requires exciting the nonlinear system under study twice using two sinusoidal inputs separately. The frequencies of the two inputs are different; one frequency is twice as high as the other one. Three numerical studies have been used to demonstrate the effectiveness of this method. The distinct advantage of this method is that it only needs the test data under two sinusoidal input forces, which can be readily carried out in practices. Since the positions of the nonlinear components in periodic structures often correspond to the locations of faults, the nonlinear component position detection method is of practical significance in the fault diagnosis for mechanical and structural systems. It is worthy to note here that in real periodic structures there are always minor non-linearities between each component and, such structures therefore should be better to model as weakly nonlinear chains which have been by Chakraborty and Mallik [START_REF] Chakraborty | Dynamics of a weakly non-linear periodic chain[END_REF]. For such structures, the problem is then changed to detect the components with strong nonlinear property in the weakly nonlinear chains. It is a more complicated problem we are planning to investigate in future studies.

A c c e p t e d m a n u s c r i p t 1 1

 1 Introduction

A c c e p t e d m a n u s c r i p t 3 In

 3 this paper, a novel method is derived based on the NOFRF concept to detect the position of the nonlinear component in a periodic structure. The detection procedure requires exciting the nonlinear systems twice using two sinusoidal inputs separately. The frequencies of the two inputs are required to be different; one frequency is twice as high as the other one. Numerical studies verify the effectiveness of the method. The new method is of great practical significance in fault diagnosis for mechanical and structural systems.

A c c e p t e d m a n u s c r i p t

1 are

 1 taken as F ω and the remainders are as F ω -, substituting equation (9) into equation (15) yields,

A c c e p t e d m a n u s c r i p t 7

 7

Figure 1 ,

 1 Figure 1, a locally nonlinear multi-degree freedom oscillator Assume the restoring forces ) (Δ LS S and ) (Δ LD S of the Lth spring and damper are the polynomial functions of the deformation Δ and Δ respectively, e.g., ∑ = Δ = Δ P

  NOFRFs up to the 4P th P order is considered, according to equations (

  type force of frequency F ω , according to the definition of NOFRFs under a harmonic input in Section 2.2, it is known from equations (27)~(29), that

A c c e p t e d m a n u s c r i p t 11 For

 11 the masses which are on the left of the nonlinear spring or on the right of the input force, substituting equation (32) into (35) yields [ ]

1 ) 1 F ω and 2 F

 112 Excite the nonlinear system by two sinusoidal inputs separately whose frequencies are

A c c e p t e d m a n u s c r i p t 15 Table 6

 156

Table 5 .Table 4

 54

	R	1 + 1 , i F i	j	ω	F	1	, which were extracted from the FFT spectra of the
	responses of the system subjected to	) u . The results of the second super-harmonics are ( 1 t
	given in Table 3 together with the calculated values of	R	1 + 1 , i F i	(	j	2 ω	F	1	)	(i=1,…,7). Table 4
	gives the first harmonics of the FFT spectra of the responses of the system subjected to

Table 5 ,

 5 Moduli)

Table 6

 6

	, and the results of	X	(	F	2	,	i	)	(	j	ω	F	2	)		and	R	, F i	1 + 2 i	(	j	ω	F	2	)	(i=1,…,8) are given
	in Table 7. Table 8 gives the moduli of	R	1 + 1 , i F i	(	j	2 ω	F	1	)	and	R	, F i	1 + 2 i	(	j	ω	F	2	)	(i=1,…,7).

Table 7 ,

 7 )

Table 8

 8

	, Moduli	, 1	1	(	2	1	)

  (i=1,…,8) are given in Table9and Table10respectively. Table11gives the moduli of

																						X	(	F	2	,	i	j	ω	F
	and	R	, F i	1 + 2 i	(	j	ω	F	2	)
										R	1 + 1 , i F i	(	j	2 ω	F	1	)	and	R	, F i	1 + 2 i	(	j	ω	F	2	)	(i=1,…,7).
										Table 9,	X	(	F	, 1	i	)	(	j	) ω and 2 1 F	R	1 + 1 , i F i	(	j	2 ω	F	1	)	(i=1,…,8)
											X	(	F	, 1	i	)	(	j	) ω ( 2 1 F	×	10 -	7	)	R	1 + 1 , i F i	(	j	2 ω	F	1	)
								i=1				-0.0691 + 0.0115i	0.5078 -0.1765i
								i=2					-0.1284 -0.0220i
																						0.5730 -0.3611i
								i=3					-0.1430 -0.1286i
																						0.5046 -0.4142i
								i=4					-0.0444 -0.2913i
																						0.4783 -0.3913i
								i=5					0.2429 -0.4103i
																						0.4842 -0.3798i
								i=6					0.7221 -0.2811i
																						-1.1682 -0.3649i
								i=7				-0.4947 + 0.3951i
																						0.7568 + 0.6112i
								i=8				-0.1405 + 0.6356i

Table 10

 10

			,			X	(	F	2	,	i	)	(	j	ω	F	2	)	and	R	, F i	1 + 2 i	(	j	ω	F	2	)	(i=1,…,8)
	X	(	F	2	,	i	)	(	j	ω	F	2	)	(	×	10 -	5	)	R	, F i	1 + 2 i	(	j	ω	F	2	)
	i=1			0.1013 + 0.1193i			0.5078 -0.1765i
	i=2			0.1052 + 0.2714i			0.5730 -0.3611i
	i=3		-0.0822 + 0.4218i			0.5047 -0.4142i
	i=4		-0.5073 + 0.41953i			1.3192 + 0.9672i
	i=5		-0.0985 + 0.3902i			1.3736 + 1.1183i
	i=6			0.0960 + 0.2059i			0.9580 + 1.2488i
	i=7			0.1409 + 0.0313i			0.7570 + 0.6107i
	i=8			0.1330 -0.0660i

Acknowledgements

The authors gratefully acknowledge the support of the Engineering and Physical Science Research Council, UK, for this work.