D J Rodrigues
email: david.rodrigues@bristol.ac.uk

A R Champneys

M I Friswell

R E Wilson

Automatic two-plane balancing for rigid rotors

We present an analysis of a two-plane automatic balancing device for rigid rotors. Ball bearings, which are free to travel around a race, are used to eliminate imbalance due to shaft eccentricity or misalignment. The rotating frame is used to derive autonomous equations of motion and the symmetry breaking bifurcations of this system are investigated. Stability diagrams in various parameter planes show the coexistence of a stable balanced state with other less desirable dynamics.

Introduction

The primary cause of vibration in rotating machinery is mass imbalance, which occurs when the principal axis of the moment of inertia is not coincident with the axis of rotation. In the case of a rigid rotor the usual balancing procedure involves adding or subtracting correction masses in two distinct planes so that the principal axis is recentred and realigned. However, effects such as thermal deformation and material erosion can cause the rotor's mass distribution to change, and in such circumstances the balancing procedure may have to be repeated. This limitation motivates the study of self-compensating balancing devices, in which masses redistribute themselves so as to eliminate any imbalance. One such device is the automatic ball balancer (ABB), which consists of a series of balls that are free to travel around a race which is filled with a viscous fluid. The first study of an ABB was carried out by Thearle in 1932 [START_REF] Thearle | A new type of dynamic-balancing machine[END_REF], and the existence of a stable steady state at rotation speeds above the first critical frequency was demonstrated. More recently, the equations of motion for a planar Jeffcott rotor with an ABB have been derived using Lagrange's method [START_REF] Lee | Analytical and experimental analysis of a self compensating dynamic balancer in a rotating mechanism[END_REF][START_REF] Chung | Dynamic analysis of an automatic dynamic balancer for rotating mechanisms[END_REF][START_REF] Green | Bifurcation analysis of an automatic dynamic balancer for eccentric rotors[END_REF]. In particular Green et al. [START_REF] Green | Bifurcation analysis of an automatic dynamic balancer for eccentric rotors[END_REF] use rotating coordinates to obtain an autonomous system of governing equations. The stability boundaries for the fully nonlinear system are then computed using numerical continuation techniques. In 1977, Hedaya and Sharp [START_REF] Hedaya | An analysis of a new type of automatic balancer[END_REF] extended the autobalancing concept by proposing a two-plane device that would be able to compensate for both shaft eccentricity and shaft misalignment, see Fig. 1. However, the ABB models which are based on a Jeffcott rotor are unable to explain phenomena that are related to shaft misalignment as they do not include any tilting motions. These out-of-plane motions are considered in the studies of Chung and Jang [START_REF] Chung | Dynamic response and stability analysis of an automatic ball balancer for a flexible rotor[END_REF], Chao et al. [START_REF] Chao | Non-planar dynamic modeling for the optical disk drive spindles equipped with an automatic balancer[END_REF] and Sperling et al. [START_REF] Sperling | Simulation of two-plane automatic balancing of a rigid rotor[END_REF], but only linear stability analyses are provided. Here we extend the planar analysis of [START_REF] Green | Bifurcation analysis of an automatic dynamic balancer for eccentric rotors[END_REF] so that we may model the tilting motion of the ABB and analyse the two-plane automatic balancer for rigid rotors. The rest of this paper is organised as follows. In section 2 we use rotating coordinates to derive autonomous equations of motion for the ABB. The steady states of the system are considered in section 3, and we focus on using numerical bifurcation theory to investigate the effect that the physical parameters have on the stability of the balanced state. In section 4 we supplement the findings of the bifurcation analysis by providing numerical simulations of the dynamics, and investigating the effect of the initial conditions. Finally, in section 5 we draw conclusions and discuss possible directions for future work.

A c c e p t e d m a n u s c r i p t 2 Derivation of the Equations of Motion

In this section, we derive an autonomous system of governing equations for the automatic balancer by describing the position of the rotor centre with respect to rotating coordinates and applying Lagrange's method [START_REF] Chung | Dynamic analysis of an automatic dynamic balancer for rotating mechanisms[END_REF][START_REF] Green | Bifurcation analysis of an automatic dynamic balancer for eccentric rotors[END_REF][START_REF] Chung | Dynamic response and stability analysis of an automatic ball balancer for a flexible rotor[END_REF]. For a background to the rotordynamics theory see for example [START_REF] Genta | Dynamics of Rotating Systems[END_REF].

Mechanical Setup

The mechanical device that we wish to model is illustrated in Fig. 1, and is based on a rigid rotor which has been fitted with a two-plane automatic balancer [START_REF] Hedaya | An analysis of a new type of automatic balancer[END_REF][START_REF] Sperling | Simulation of two-plane automatic balancing of a rigid rotor[END_REF]. The rotor has mass M , moment of inertia tensor J, and is mounted on two compliant linear bearings which are located at S 1 and S 2 . The automatic balancer consists of a pair of races that are set normal to the shaft in two separate planes. Each race contains two balancing balls of mass m, which move through a viscous fluid and are free to travel, at a fixed distance R from the shaft axis. The position of the ith ball is specified by the axial and angular displacements z i and α i , which are written with respect to the Cξηz rotor axes.

l 2 l 1 z 1,2 -z 3,4 α 1 α 2 α 3 α 4 R R C ξ η z k 1 k 2 c 1 c 2 S 1 S 2
Fig. 1: Schematic diagram of a two-plane automatic balancer.

Co-ordinate frames

In order to describe the position and orientation of the rotor, it is helpful to consider the following frames of reference, as shown in Fig. 2. We begin with an inertial space frame OXY Z with origin at O and Z-axis oriented along the undeflected bearing centreline. A rotation about the Z axis by the spin angle θ results in the rotating frame OxyZ. This transformation can be written as

x = R 1 X,
where x and X are the column vectors of coordinates in the rotating and inertial frames respectively, and

R 1 =   cos θ sin θ 0 -sin θ cos θ 0 0 0 1   .
The torsional behaviour of the rotor lies outside the scope of the present study, and so we shall only consider the special case of constant-speed operation, in which θ = Ωt. The rotor's lateral motion can be described by introducing a frame Cx ′ y ′ z ′ with origin at the geometric shaft centre C, and axes parallel to those of the rotating frame OxyZ. We neglect any motion in the axial direction, and so, the position vector of the geometric centre r C lies in the x-y plane. The rotor may also perform an out-of-plane tilting motion that can be described as follows: firstly we define an intermediate axes Cx ′′ y ′′ z ′′ related to Cx ′ y ′ z ′ by a rotation of an angle φ y about the y ′ axis, then we rotate Cx ′′ y ′′ z ′′ about x ′′ by an angle φ x , which results in a body frame Cξηz that is fixed with respect to the rotor. These transformations can be combined to give

ξ = R 3 R 2 x ′ , (1)
where ξ and x ′ are the co-ordinates in the body and primed axes respectively and

R 2 =   cos φ y 0 -sin φ y 0 1 0 sin φ y 0 cos φ y   , R 3 =   1 0 0 0 cos φ x sin φ x 0 -sin φ x cos φ x   .
We are now able to transform between body coordinates and the space reference frame, however as shown in Fig. 3, small errors in the rotor's mass distribution will cause the body axes Cξηz to differ from the principal axes of the moment of inertia. The eccentricity ǫ, which gives rise to the static imbalance, is defined as the distance between the shaft centre C and the rotor's centre of mass G. Also the symmetry axis p 3 corresponding to the polar moment of inertia may be misaligned to the shaft axis by an angle χ, and this results in a couple imbalance. The symmetry of the rotor enables us to take the misalignment to be about the η axis without detracting from the generality of the model. The matrix expressing this Before moving on to the derivation of the Lagrangian, it is worth noting that we have chosen the order of rotations so that the spin Ω is taken about the space Z axis as opposed to the body z axis. In the standard formulation for a linear four degree of freedom rotor, both approaches lead to equivalent equations of motion [START_REF] Genta | Dynamics of Rotating Systems[END_REF]. However, for our case where we must include the geometric nonlinearities due to the balancing balls, we find that autonomous equations can only be derived when Ω is taken about the space Z axis [START_REF] Rodrigues | Automatic balancing of a rigid rotor with misaligned shaft[END_REF].

Lagrangian formulation

The nonlinear equations of motion for the system are derived from Lagrange's equations

d dt ∂T ∂ qk - ∂V ∂ qk - ∂T ∂q k - ∂V ∂q k + ∂F ∂ qk = 0, (2)
where T , and V , are the kinetic and potential energies respectively, F is Rayleigh's dissipation function and the q k are the generalised co-ordinates. In our case the generalised coordinates are given by q = (x, y, φ x , φ y , α 1 , α 2 , α 3 , α 4 ) , and so the motion is governed by a total of eight independent equations of motion. We proceed by deriving in turn the kinetic energy, the potential energy and Rayleigh's dissipation function.

Kinetic Energy

The kinetic energy of any rigid body can be decomposed as the sum of the translational energy of the centre of mass and the rotational energy of the body about its centre of mass. By treating the balls as points, the kinetic energy of the system can be written as

T = 1 2 M ṙ2 G + 1 2 Ω T JΩ + 1 2 4 i=1 m ṙ2 Bi . (3)

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

where ṙG and ṙBi are the velocities of the rotor's centre of mass and the ith ball respectively, and Ω is the angular velocity of the rotor. By using the inverse of the transformation defined in (1), we may write the position vector of the centre of mass, in the rotating frame as

r G =   x y 0   + R T 2 R T 3   ǫ cos β ǫ sin β 0   , =  
x + ǫ (cos β cos φ y + sin β sin φ x sin φ y ) y + ǫ sin β cos φ x ǫ (-cos β sin φ y + sin β sin φ x cos φ y )   .

The rotating coordinates x and y have been chosen to describe the lateral motion because their use enables an autonomous formulation in which the the spin angle Ωt does not explicitly enter into the governing equations. The position vector of the i-th ball can be similarly given by

r Bi =   x y 0   + R T 2 R T 3   R cos α i R sin α i z i   , =   x + R (cos φ y cos α i + sin φ x sin φ y sin α i ) + cos φ x sin φ y z i y + R cos φ x sin α i -sin φ x z i R (-sin φ y cos α i + sin φ x cos φ y sin α i ) + cos φ x cos φ y z i   .
In contrast to the eccentricity ǫ, the race radius R cannot be assumed to be a small quantity. As a consequence, we must take extra care in terms relating to the ball positions when making any small angle assumptions on φ x and φ y . Thus, we prefer to compute the kinetic energy using the complete expressions for the displacements and make any approximations only after the equations of motion have been derived. When calculating the translational kinetic energies, it is also necessary to express the velocities as they are seen by an observer in a space-fixed frame. Since we have written the position vectors with respect to rotating coordinates, we must use the following operator relation between derivatives to account for the change in frames

d dt s = d dt r + Ω r × .
Here the subscripts s and r denote the time derivatives as seen by an observer in the space and rotating frames respectively, and Ω r = [0, 0, Ω] T is the angular velocity of the rotating frame with respect to an observer fixed in space. Hence, the velocity of the centre of mass in the space frame is given by

ṙG = dr G dt s = dr G dt r + Ω r ×r G , =   ṙG,1 ṙG,2 ṙG,3   , (4)
where ṙG,1 = ẋ + ǫcos β sin φ y φy + sin β cos φ x φx sin φ y + sin φ x cos φ y φy -Ω (y + ǫ sin β cos φ x ) , ṙG,2 = ẏǫ sin β sin φ x φx + Ω (x + ǫ (cos β cos φ y + sin β sin φ x sin φ y )) , ṙG,3 = ǫcos β cos φ y φy + sin β cos φ x φx cos φ y -sin φ x sin φ y φy .

The velocities of the balls are computed in the same manner, so that

ṙBi =   ṙBi,1 ṙBi,2 ṙBi,3   , (5)
A c c e p t e d m a n u s c r i p t with ṙBi,1 = ẋ + Rsin φ y φy cos α i -cos φ y sin α i αi + cos φ x φx sin φ y sin α i + sin φ x cos φ y φy sin α i + sin φ x sin φ y cos α i αi + -sin φ x φx sin φ y + cos φ x cos φ y φy z i

-Ω (y + R cos φ x sin α i -sin φ x z i ) , ṙBi,2 = ẏ + Rsin φ x φx sin α i + cos φ x cos α i αicos φ x φx z i +Ω (x + R (cos φ y cos α i + sin φ x sin φ y sin α i ) + cos φ x sin φ y z i ) , ṙBi,3 = Rcos φ y φy cos α i + sin φ y sin α i αi + cos φ x φx cos φ y sin α i -sin φ x sin φ y φy sin α i + sin φ x cos φ y cos α i αi + -sin φ x φx cos φ y -cos φ x sin φ y φy z i .

Next we turn to the expression for the rotational energy, which is most easily derived by using the principal axes Gp 1 p 2 p 3 of Fig. 3. It is in this frame that the inertia tensor J, takes the diagonal form

J =   J t 0 0 0 J t 0 0 0 J p   ,
in which J t and J p are the transverse and polar moments of inertia respectively. The angular velocity Ω is composed from the following three angular velocity vectors: θ directed along the Z-axis, φy which is directed along the y ′ axis and φx which lies along the x ′′ axis. The rotation matrices of Section 2.2 can now be used to transform these vectors so that they are written with respect to the principal axes, and taking the sum we obtain the angular velocity vector Ω as

Ω = R 4     φx 0 0   + R 3   0 φy 0   + R 3 R 2 R 1   0 0 θ    , =    
cos χ φxθ sin φ y -sin χφy sin φ x + θ cos φ x cos φ y φy cos φ x + θ sin φ x cos φ y sin χ φxθ sin φ y + cos χφy sin

φ x + θ cos φ x cos φ y     .
Hence the rotational energy can now be easily computed as

1 2 Ω T JΩ = 1 2 J t cos χ φx -Ω sin φ y -sin χ -φy sin φ x + Ω cos φ x cos φ y 2 + φy cos φ x + Ω sin φ x cos φ y 2 (6) + 1 2 J p sin χ φx -Ω sin φ y + cos χ -φy sin φ x + Ω cos φ x cos φ y 2 .
Finally substituting the equations (4-6) into equation ( 3) completes the derivation of the kinetic energy for the system.

Potential Energy and Damping

The potential energy V , arises from the elastic deflection of the support bearings at S 1 and S 2 , as shown in Fig. 1. The coordinates of these supports are given in the rotor-fixed frame by [0, 0, l 1 ] T and [0, 0, -l 2 ] T , and we assume that each bearing is linear and isotropic with stiffness k 1 , k 2 and damping c 1 , c 2 respectively. By transforming coordinates to the space frame, the deflection d Si of each support may be written as

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

d Si = R T 1     x y 0   + R T 2 R T 3   0 0 li     -   0 0 li   , =     cos θ x + cos φ x sin φ y li -sin θ y -sin φ x li sin θ x + cos φ x sin φ y li + cos θ y -sin φ x li cos φ x cos φ y li -li     , where l1 = l 1 and l2 = -l 2 .
Hence the potential energy can be computed as

V = 1 2 2 i=1 k i d 2 Si , = 1 2 k 1 (x + cos φ x sin φ y l 1 ) 2 + (y -sin φ x l 1 ) 2 + (cos φ x cos φ y l 1 -l 1 ) 2 (7) + 1 2 k 2 (x -cos φ x sin φ y l 2 ) 2 + (y + sin φ x l 2 ) 2 + (-cos φ x cos φ y l 2 + l 2 ) 2 .
It is usual to assume that the vibrational coordinates (φ x , φ y ) together with the imbalance errors ǫ and χ are small and O(δ), with δ a small quantity. If we are to neglect terms of O(δ 2 ) in the equations of motion, then all terms of order O(δ 3 ) which enter into Lagrange's equation will not contribute to the system, and can be ignored. A routine calculation now gives the more familiar expression

V = 1 2 k 11 x 2 + y 2 + k 12 (xφ y -yφ x ) + 1 2 k 22 φ 2 y + φ 2 x + O(δ 3 ),
where

k 11 = k 1 + k 2 , k 12 = k 1 l 1 -k 2 l 2 , k 22 = k 1 l 2 1 + k 2 l 2 2 .
Next we turn to the dissipative effects, which arise from the damping at the supports, together with the viscous drag on the balls c b , as they pass through the fluid in the race. These terms are most easily calculated through the use of Rayleigh's dissipation function

F = 1 2 2 i=1 c i ḋ2 Si + 1 2 4 i=1 c b αi 2 ,
which after some manipulation we may write as

F = 1 2 c 11 ( ẋ -Ωy) 2 + ( ẏ + Ωx) 2 +c 12 ( ẋ -Ωy) φy + Ωφ x + (y + Ωx) -φx + Ωφ y (8) + 1 2 c 22 φy + Ωφ x 2 + -φx + Ωφ y 2 + 1 2 4 i=1 c b αi 2 + O δ 3 ,
where

c 11 = c 1 + c 2 , c 12 = c 1 l 1 -c 2 l 2 , c 22 = c 1 l 2 1 + c 2 l 2 2 .
The fully geometric nonlinear equations of motion are generated by substituting (3), ( 7) and ( 8) into Lagrange's equations. Before displaying them here, we shall first continue with non-dimensionalisation.

A c c e p t e d m a n u s c r i p t 2.3.3 Non-dimensionalisation and the Equations of Motion

We now consider the following dimensionless state variables

x = x R , ȳ = y R , φx = φ x , φy = φ y , ᾱi = α i for i = 1 . . . 4,
and the dimensionless time

t = Ω 1 t.
Here Ω 1 is the critical frequency associated with cylindrical whirl, and is given by, Ω 1 = k 11 /M . Furthermore, we introduce the dimensionless parameters

Ω = Ω Ω 1 , ǭ = ǫ R , m = m M , z = z R , J = J M R 2 , Jp = J p M R 2 , χ = χ, cb = c b mΩ 1 R 2 , k11 = k 11 M Ω 2 1 ≡ 1, k12 = k 12 M Ω 2 1 R , k22 = k 22 M Ω 2 1 R 2 , c11 = c 11 M Ω 1 , c12 = c 12 M Ω 1 R , c22 = c 22 M Ω 1 R 2 .
The equations of motion with respect to these non-dimensional quantities are computed with the aid of the computer algebra system Maple, and are given as follows up to order O(δ 2 ). The overbars have been omitted for clarity and we start with the equations for the rotor motion, which are

M 0 0 M ẍ + C -Ω (2M -G) Ω (2M -G) C ẋ + K -Ω 2 (M -G) -ΩC ΩC K -Ω 2 (M -G) x = Ω 2     ǫ cos β χ (J t -J p ) ǫ sin β 0     + m 4 i=1     (Ω + αi ) 2 0 αi 0 0 (Ω + αi ) 2 0 αi -αi 0 (Ω + αi ) 2 0 0 -αi 0 (Ω + αi ) 2         cos α i z i cos α i sin α i z i sin α i     +     0 n φy 0 n -φx     . (9)
where

M = 1 + 4m i mz i i mz i J t + i mz 2 i , G = 0 0 0 J p , C = c 11 c 12 c 12 c 22 , K = k 11 k 12 k 12 k 22 ,
are the mass, gyroscopic, damping and stiffness matrices respectively and x = [x, φ y , y, -φ x ] T is the vector of the rotor degrees of freedom. We have used -φ x instead of φ x so that (9) assumes a more regular pattern. Finally the n φy and n -φx terms are given by

n φy = m 4 i=1 ẍ -2Ω ẏ -Ω 2 x (φ y cos α i -φ x sin α i ) -cos 2 α i φy -2 αi φx + Ω 2 + 2Ω αi φ y + sin α i cos α i φx + 2 αi φ y + Ω 2 + 2Ω αi φ x + αi φ x , n -φx = m 4 i=1 ẍ -2Ω ẏ -Ω 2 x φ y sin α i + sin 2 α i φx + 2 αi φy + Ω 2 + 2Ω αi φ x -ÿ + 2Ω ẋ -Ω 2 y φ x sin α i -sin α i cos α i φy -2 αi φx + Ω 2 + 2Ω αi φ y .

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

Coupled with ( 9) are the equations of motion for the balancing balls which are

αi + c b αi = ẍ + z i φy -2Ω ẏ -z i φx -Ω 2 (x + z i φ y ) sin α i -ÿ -z i φx + 2Ω ẋ + z i φy -Ω 2 (y -z i φ x ) cos α i , for i = 1, . . . 4. (10)
We note that by taking m = 0 in (9), we recover the equations of motion for a four degree of freedom rotor in the rotating frame [START_REF] Genta | Dynamics of Rotating Systems[END_REF]. Similarly, by setting the tilt angles φ x = φ y ≡ 0, the system reduces to the equations of motion for the planar automatic balancer [START_REF] Green | Bifurcation analysis of an automatic dynamic balancer for eccentric rotors[END_REF].

It is also worth remarking that there are certain asymmetries in the n φy and n -φx terms corresponding to the O(δ) moments of the balancing balls. These arise because of the inherently asymmetrical nature of an Euler angle formulation. With reference to Fig. 2, we recall that φ x is a rotation about the body ξ axis while φ y is taken, not about the body η axis, but about an intermediate y ′ axis. With the standard rotordynamic assumptions of small tilt angles and small eccentricities these asymmetries do not survive the linearisation process. Here however, where the lever arm of the race cannot be neglected, extra terms including these asymmetries are present [START_REF] Genta | Some considerations on the basic assumptions in rotordynamics[END_REF].

Steady state bifurcation analysis

In this section we consider the stability of the steady state solutions of the system given by equations [START_REF] Genta | Dynamics of Rotating Systems[END_REF][START_REF] Rodrigues | Automatic balancing of a rigid rotor with misaligned shaft[END_REF]. In order to simplify the bifurcation analysis, we shall restrict attention to a symmetric subsystem which displays all the primary features of the autobalancing process. A situation that is common in practice occurs when the rotor's centre of mass is located at the midspan and both supports are identical, so that

k 1 = k 2 = k, c 1 = c 2 = c, l 1 = l 2 = l, which leads to k 12 = c 12 = 0.
With this set-up, the lateral and inclinational motions of the rotor remain coupled, but only through the motion of the balancing balls. In addition, we assume z 1,2 = -z 3,4 = z, so that the autobalancing planes are equidistant from the midspan. The rest of this section is organised as follows: in Section 3.1 we obtain the steady state solutions, noting various conditions governing their existence, and in Section 3.2 we use the continuation package AUTO [START_REF] Doedel | AUTO97: Continuation and bifurcation software for ordinary differential equations[END_REF] to compute bifurcation diagrams showing the boundaries of stability in various parameter planes.

Steady state solutions

Steady state solutions are obtained by setting all time derivatives in the equations of motion (9-10) to zero. Moreover, if we also set the vibrational coordinates (x, y, φ x , φ y ) = 0, we arrive at the following conditions for a balanced steady state

cos α 1 + cos α 2 + cos α 3 + cos α 4 = s x , sin α 1 + sin α 2 + sin α 3 + sin α 4 = s y , cos α 1 + cos α 2 -cos α 3 -cos α 4 = c x , sin α 1 + sin α 2 -sin α 3 -sin α 4 = c y , (11)
where

s x = -ǫ cos β m , s y = -ǫ sin β m , c x = -χ (J t -J p ) mz , c y = 0.
These are precisely the conditions for equilibrium between the centrifugal forces and moments acting on the rotor due to the imbalance and balancing balls [START_REF] Sperling | Simulation of two-plane automatic balancing of a rigid rotor[END_REF]. By rewriting the above equation with respect to

where

ᾱ12 = α 1 + α 2 2 , α12 = α 1 -α 2 2 , ᾱ34 = α 3 + α 4 2 , α34 = α 3 -α 4 2 .
Because the balls are assumed to be identical, we may without loss of generality take 0 ≤ α12 , α34 ≤ π 2 and equation ( 12) can now be solved to give the following physically unique balanced state ᾱ12 = arctan

s y + c y s x + c x + 0 if s x + c x ≥ 0, π if s x + c x < 0. , α12 = arccos 1 4 (s x + c x ) 2 + (s y + c y ) 2 , ᾱ34 = arctan s y -c y s x -c x + 0 if s x -c x ≥ 0, π if s x -c x < 0. , α34 = arccos 1 4 (s x -c x ) 2 + (s y -c y ) 2 .
(

) 13
This solution exists provided that the arguments for both the arccoses have modulus less that one, which ensures that the balls have enough mass to counteract the imbalance of the system. In what follows we refer to (13) as the balanced state B. The steady state configurations for the special case of a rotor with a static imbalance are illustrated in Fig. 4. The balanced state B is shown on the left, but there are also various other steady state solutions for which the rotor remains out of balance. In these cases we consider the equation for the balls [START_REF] Rodrigues | Automatic balancing of a rigid rotor with misaligned shaft[END_REF], which yields

(x + z i φ y ) sin α i = (y -z i φ x ) cos α i , i = 1 . . . 4, (14)

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

where x + z i φ y and yz i φ x , are recognised as the x and y deflections of the race centres. If for both of the races (x + z i φ y , yz i φ x ) = 0, then we may rewrite (14) as

tan α i = y -z i φ x x + z i φ y , i = 1 . . . 4,
with solutions

α 1 = α 2 + pπ, α 3 = α 4 + qπ, p , q ∈ Z.
Therefore the balls in each race are either coincident, or in-line with the race centre and opposite each other. For example, in state C1 all the balls have the same angle and this configuration results in the excitation of a cylindrical whirl. For state C2, the balls in each race are again coincident but there are now two distinct angular displacements. This arrangement generates a couple imbalance and leads to a conical type whirl. The inline states will not be discussed further because they are found to be always unstable and are never involved in bifurcations with the balanced state.

Next we consider the state C3, in which (x + z i φ y , yz i φ x ) = 0 is satisfied by one of the two races, so that its centre stays fixed at the undeflected position. The balls within this fixed race are now not constrained by equation ( 14), although the balls in the opposing race are again either coincident or inline.

As an example, the arrangement for C3 depicted in Fig. 4 shows coincident balls in the top race while the bottom race balls are split. The corresponding state, where the roles of the races are swapped, is physically equivalent, because the symmetry between the two races is preserved by a static imbalance. However, this will not be the case for a general imbalance, as we shall discuss in section 3.2.2. Finally, we note that if (x + z i φ y , yz i φ x ) = 0, for both of the races, then we are again led to the balanced state B, which has been covered previously.

Numerical continuation

For the remainder of this study we shall consider a rotor with the following inertial and stiffness parameters

J t = 3.25, J p = 0.5, l = 3, k 11 ≡ 1, k 12 = 0, k 22 = 9. ( 15
)
These values are based on a solid cylindrically shaped rotor with a height of six times its radius. A rotating machine on compliant bearings with this geometry would typically undergo a two-plane balancing procedure before going into service. The approximate natural frequencies for the cylindrical and conical modes occur respectively at

Ω 1 = k 11 M ≡ 1, and Ω 2 = k 22 J t -J p ≃ 1.81, (16)
and we shall now investigate the results for both the static and dynamic imbalance cases.

Static Imbalance

As mentioned above, a static imbalance occurs when the principal axis is displaced, but parallel to, the rotation axis so that there is a shaft eccentricity ǫ = 0, but no shaft misalignment χ = 0. For this special case, we are able to recover all the bifurcations and steady state solutions that have been found previously in the studies of a planar automatic balancer [START_REF] Chung | Dynamic analysis of an automatic dynamic balancer for rotating mechanisms[END_REF][START_REF] Green | Bifurcation analysis of an automatic dynamic balancer for eccentric rotors[END_REF]. However, because our model includes the extra inclinational degrees of freedom, the rotor is now able to destabilise via out of plane motions, and so, the steady states will necessarily have smaller regions of stability.

In Fig. 5 we present results for the following parameters

(i) (ii) (ii) (iii) (iii) (a) (b) (c) (d) (e) (e) (e)
(b,c) m m m m m m c m c Ω Ω Ω 1 Ω 2 α 1 L 2 L 2 PF B,C1,C3 PF B,C1,C3 PF B,C1,C3 PF B,C1,C3 PF B,C1,C3 H con H con H con H cyl

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

For this value of the eccentricity, the balanced state B exists provided that

m ≥ m c := ǫ 4 = 0.0025, (18)
where m c is the critical value of the ball mass. We also note that the balancing planes are appropriately spaced, and that the race and supports have low damping values, which match those in [START_REF] Chung | Dynamic analysis of an automatic dynamic balancer for rotating mechanisms[END_REF][START_REF] Green | Bifurcation analysis of an automatic dynamic balancer for eccentric rotors[END_REF]. Figure 5(a) shows the results of a two-parameter bifurcation analysis upon variation of the dimensionless parameters Ω and m, whilst we keep the other parameters fixed. The shaded regions correspond to the existence of a stable balanced state B, although as we shall discuss in section 4, other less desirable stable states may coexist in these regions. For masses of the balls which satisfy m ≥ m c the balanced state B always exists, but it is not necessarily stable. The curves marked H cyl and H con are Hopf bifurcations of B, at which oscillations are induced by the cylindrical and conical modes respectively. In a planar analysis only the H cyl curves are present, and as such, they are always accompanied by a change in stability. Here with the inclusion of the inclinational degrees of freedom, the stability of B is predominately determined by the H con curves, at which the conical type oscillations, are either created or destroyed.

For large values of Ω, the balanced state B stabilises in a Hopf bifurcation H con and there is also a secondary region of stability, for sufficiently large values of m, at rotation speeds just above the natural frequency Ω 2 . We also draw attention to the similarities between the shape of the two Hopf curves H cyl and H con , suggesting that there is in some sense, a duality between the lateral and inclinational dynamics of the automatic balancer.

The one-parameter bifurcation diagrams of Figs. 5(b) and (c) show the result of increasing m through m c for Ω = 1.5, which is a value of the rotation speed that lies in between the natural frequencies of the rotor. The angular displacement α 1 is used as the solution measure in Fig. 5(b), and we see that an unstable balanced state B is born at m = m c in a pitchfork bifurcation involving the unstable coincident state C1. Note that the two copies of the balanced state B come from interchanging the ball positions, for example α 1 ↔ α 2 . In Fig. 5(c) we plot the same bifurcation diagram as shown in panel (b) but instead we use the Euclidean norm L 2 . Here only one copy of each physically unique state is generated, and we can see that the balanced state B is in fact produced at a degenerate pitchfork bifurcation involving the states C1 and C3.

Next we consider the situation for the supercritical regime where Ω > Ω 2 . Figure 5(d) is an enlargement of our region of interest and panel (e) shows a one parameter scan as we increase m through m c , for a fixed rotation speed of Ω = 3. Again we find that the balanced state B is born at a degenerate pitchfork bifurcation involving states C1 and C3, but in contrast with the situation for subcritical rotation speeds, the coincident state C1 and the bifurcating balanced state B are now both stable. Note also, that the balanced state changes stability through a series of three further Hopf bifurcations, which is in agreement with the vertical slice taken through panel (d) at Ω = 3. Now let us consider some stability diagrams similar to Fig. 5(a), in which we vary other dimensionless parameters. We take m = 0.025 so that condition (18) needed for balance is satisfied, and the values for the other parameters when fixed, are again given by ( 15) and ( 17).

In Fig. 6 we display some two-parameter stability charts upon variation of the eccentricity ǫ, the support damping c, and the race damping c b , as the rotation speed Ω is increased. Here, we are again able to recover the H cyl curves which are present in the planar analyses, but as before the H con bifurcations have more influence on the stability of the balanced state. We also highlight some of the same features which were present in Fig. 5, namely: the stabilisation of the balanced state through a H con Hopf bifurcation for high enough rotation speeds, the smaller secondary stable region at rotation speeds just above Ω 2 , and most notably in panel (c) the similarity between the shapes of the H cyl and H con curves.

Figure 7 shows the eccentricity ǫ, plotted against Ω, whilst we also vary the ball mass so that m = 4m c = ǫ. Thus, for any given value of the eccentricity, the balls have a mass which is four times the minimum amount that is required to balance the system. A logarithmic scale is used for the vertical axis so that a wide range of eccentricities can be considered. Here, the main area of interest occurs for small eccentricities and supercritical rotation speeds where there is a large connected stable region. We also note that the H con curve, which bounds this region, asymptotes towards Ω = Ω 2 as ǫ → 0 so that there is no stable region in the subcritical regime. Finally, before moving on to a general dynamic imbalance, we shall briefly discuss the case of a couple imbalance, for which there is a shaft misalignment χ = 0 but no shaft eccentricity ǫ = 0. The steady states for this case are illustrated schematically in Fig. 8. The only difference between the static imbalance steady states of Fig. 4 is that here the bottom race has been rotated through an angle of π radians. Again, the symmetry between the two races is preserved by this type of imbalance, and the balanced state B is always born in a degenerate pitchfork bifurcation involving the coincident states C1 and C3. Moreover if the misalignment parameter χ, is appropriately chosen so that we again have m c = 0.0025, the results shown in Fig. 5(a) for a static imbalance can be recovered. Thus, there is an equivalence between the properties of the automatic balancer for a pure eccentricity and for a pure misalignment, and we shall not provide any explicit bifurcation diagrams for the couple imbalance case.

(c) ǫ Ω Ω Ω Ω 1 Ω 2 c c b H con H con H con H con H con H cyl H cyl H cyl H cyl H cyl

Dynamic Imbalance

A dynamic imbalance arises when there is both a shaft eccentricity ǫ = 0, and a shaft misalignment χ = 0. A phase angle β must also be included to describe the difference between the directions of the two imbalances. The parameter set which was given in (15) and ( 17) is used again, except that we change the imbalance The critical mass for the balls m c , is now given by

m c := max 1 4 ǫ cos β ± χ (J t -J p ) z 2 + (ǫ sin β) 2 ≃ 0.00261, (19)
which is a similar value to that in the static imbalance case of the previous section. In Fig. 9(a) we provide a two-parameter bifurcation diagram upon the variation of the ball mass m and the rotation speed Ω. Qualitatively the graph is similar to the static case of Fig. 5(a), in that the balanced state also stabilises in a Hopf bifurcation H con , for high enough rotation speeds in the supercritical regime. However the secondary stable region has now almost disappeared, and we see a sharp upward spike in the lower H cyl curve which is indicative of mixing between the cylindrical and conical modes. In panel (b) we scan vertically through Fig. 9(a) holding the rotation speed fixed at (20)

m m Ω Ω 1 Ω 2 PF B,Ct3 10 -3 (a) (b)
Physically this measure corresponds to the average vibration level at the points on the shaft that are a unit length from the rotor centre. We see that the degenerate pitchfork bifurcation which was present for the static imbalance case has been split into three separate non-degenerate pitchforks. For a generic dynamic imbalance this unfolding occurs because the symmetry between the two races has been broken. As a consequence, the unbalanced steady state C3 separates into two physically different cases. The configuration with coincident balls in to top race is called Ct3, whereas the corresponding state with the bottom race balls coincident is Cb3. We find that the symmetry and resulting degenerate bifurcation is only preserved when there is a pure static or couple imbalance or if the phase is such that cos β = 0. In all these cases there is no component of the shaft eccentricity which lies in the direction of the misalignment. Thus from (13) we can see that α12 ≡ α34 and so the balanced state B is born from a state where the balls in both races are coincident. Finally in Fig. 10 we provide a two-parameter bifurcation analysis as we vary the support damping c versus Ω (a), and the race damping c b versus Ω (b). Again we choose m = 0.025 so that condition (19) for the existence of the balanced state is satisfied. By comparison with Fig. 6 we note that in both cases the secondary stable regions have been lost. Also the stable region for high rotation speeds and high support damping c is dramatically reduced.

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

We build on the bifurcation analysis of the previous section by presenting some results which illustrate the important dynamics of the automatic ball balancer (ABB). Here we use the Matlab routine ode45 to perform a direct numerical integration of the equations of motion (9-10), which were derived in section 2. As a check, we note that our simulations agree with those of [START_REF] Chung | Dynamic analysis of an automatic dynamic balancer for rotating mechanisms[END_REF][START_REF] Green | Bifurcation analysis of an automatic dynamic balancer for eccentric rotors[END_REF] for appropriate choices of the parameters and initial conditions. Figure 11 shows the effect of varying the system parameters for the static imbalance case that was discussed in section 3.2.1. The positions in parameter space of these runs are highlighted in Figs. 5 and7. We plot the vibration measure Ā against the dimensionless time t, here the black curve represents the motion of the rotor with the ABB and the grey curve the motion of the rotor without the ABB. The initial conditions for all four plots are the same and are given by Therefore the rotor starts at rest in the undeflected position and the balls are initially stationary (with respect to the rotating frame) and are placed on opposite sides to each other so that they do not add

α 1 = -α 2 = α 3 = -α 4 = π/
(iv) (i) (ii) (iii) Ā Ā Ā Ā t t t t

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

to the imbalance of the rotor. Panel (i) shows the simulation for (Ω, ǫ, m) = (2, 0.01, 0.025) which, with reference to Fig. 5, lies in the secondary region of stable balanced operation, at a rotation speed just above the natural frequency Ω 2 . We note that although the autobalancing state is eventually achieved, the transients produce slightly higher levels of initial vibrations. Also, if we perturb the initial conditions so that the shaft is slightly tilted, we find that the rotor can easily destabilise into a periodic type motion with significantly worse levels of vibration (not shown but qualitatively similar to that of (ii)). This effect may be understood by considering the results of Fig. 9, where we commented on the disappearance of the secondary stable region upon the introduction of a shaft misalignment. The initial tilt of the rotor leads to a motion of the balls that produces a small couple imbalance, which in turn destabilises the operation of the ABB. Figure 11 (ii) shows the time profile for (Ω, ǫ, m) = (2.5, 0.01, 0.005). This point lies outside the region of stability for the balanced state B, however the transient response dies down to a very low vibration level before eventually destabilising into periodic motion. This type of behaviour is typical of the dynamics near a saddle type fixed point when there are symmetrical initial conditions. In panel (iii) our parameters are (Ω, ǫ, m) = (4, 0.01, 0.005), which lies in the larger primary region of stable balanced operation. For this case, balance is achieved after t ≈ 400 and the transient period for both curves are approximately the same. Moreover the dynamics here are more robust than in (i) to perturbations in the initial tilt of the rotor. Again, with particular reference to Figs. 5 (a) and 9 (a), we observe that the primary region of stability is not greatly influenced by the inclusion of a couple imbalance. Thus initial motions of the balls, which cause small inclinational motions, tend not to destabilise the device. Finally, in panel (iv) we illustrate the behaviour for (Ω, ǫ, m) = (1.5, 0.025, 0.025) which, as shown in Fig. 7, lies well outside any region of balanced operation. This is a very undesirable state as the rotor now appears to undergo chaotic motion and reaches levels of vibration which are far higher than that for the rotor without an ABB mechanism. Now we turn to the effect that the initial conditions have on the operation of the ABB. It was found that there was little sensitivity to typical variations in the initial values of the vibrational coordinates (x, y, φ x , φ y ) or their velocities ( ẋ, ẏ, φx , φy ), however the initial value of the ball positions α i can have a dramatic impact on the outcome of the system. We consider the dynamic imbalance case of section 3.2.2 where (ǫ, χ, β) = (0.005, 0.005, 1.0). In particular we take (Ω, m) = (3.5, 0.005), which is indicated in Fig. 9 and lies in the stable balanced region. We display in Fig. 12 the vibration level Ā for the rotor with and without the ABB mechanism and the angular displacement of the balancing balls α 1 , α 2 (black curves) and α 3 , α 4 (grey curves) are shown underneath. Panel (a) illustrates the results for the same initial conditions that were used for Fig. 11; namely those given by (21). Here the balls get attracted to the balanced state B with angular positions α i given by (13).

In panel (b) we take the initial conditions α 1 = -α 2 = π/3, α 3 = -α 4 = 2π/3, this has the effect of increasing the initial couple imbalance of the rotor. Here the device settles down to a periodic motion where the vibration level is about twice that of the rotor without the ABB. We find that one of the balls in the top race lags behind the rotor with α = -Ω. Thus this ball only performs one revolution for every two of the rotor. Note that a model which included partitions of the race [13] or allowed for collisions between balls would prevent such dynamics. Finally in Fig. 12 (c) we consider starting conditions where the initial couple balance is maximised so that α 1 = α 2 = 0, and α 3 = α 4 = π. In this case the rotor appears to undergo a chaotic motion and the vibration level of the rotor is an order of magnitude worse.

Conclusion

We have provided the first nonlinear bifurcation analysis of a two-plane automatic balancer for rigid rotors. It has been demonstrated that the use of rotating coordinates enables an autonomous formulation for the equations of motion. Two-parameter stability charts obtained by numerical continuation show that the considered device can effectively eliminate imbalances arising from both shaft eccentricity and shaft misalignment. However, the balancing process only works for sufficiently high rotation speeds above the second critical frequency [START_REF] Sperling | Simulation of two-plane automatic balancing of a rigid rotor[END_REF].

In our investigations, we have found all the bifurcations which appear in the planar analysis of the ABB [START_REF] Chung | Dynamic analysis of an automatic dynamic balancer for rotating mechanisms[END_REF][START_REF] Green | Bifurcation analysis of an automatic dynamic balancer for eccentric rotors[END_REF], together with additional bifurcations that occur as a result of the out-of-plane motions. Furthermore, we have highlighted how the symmetry properties of the imbalance affects the degeneracy of the bifurcation at which the balanced state is born.

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

In addition, we have presented numerical simulations which show the coexistence of the balanced state with other less desirable dynamics. Here as expected, the probability of successful operation is maximised by releasing the balls from directly opposite points on the race so that the initial imbalance is not increased.

The model introduced in this paper has a number of assumptions that may not be valid for an autobalancer operating on a real machine. We plan to continue the present work by investigating the performance of the automatic balancer on rotors which have asymmetric support characteristics or where the balancing planes are not equidistant from the midspan. Also more research needs to be carried out in order to determine the influence of ball-race interactions and rotor flexibility on the autobalancing process. However it is clear that the ABB is a highly nonlinear device, and so the regions of its stability depend in a non trivial way on the parameters of the system. As a result, we expect that the methods provided here will influence the design requirements for the ABB, in order that it may successfully emulate the two-plane balancing procedure.

Fig. 2 :

 2 Fig. 2: Definition of the co-ordinate system.

3 CFig. 3 :

 33 Fig. 3: Definition of the imbalance.

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

 4 the average and differential angular displacements within each race, we find that cos α12 cos ᾱ12 = s x + c x 4 , cos α12 sin ᾱ12 = s y + c y 4 , cos α34 cos ᾱ34 = s x -c x cos α34 sin ᾱ34 = s y -c y 4 ,

Fig. 4 :

 4 Fig. 4: Schematic diagram of the steady states for the case of a static imbalance. White surrounds denote that both balls in the race are coincident.

ǫ

  = 0.01, χ = 0, (β = 0), z = 2, c b = 0.01, and c = 0.01, so that c 11 = 0.02, c 12 = 0, c 22 = 0.18.

Fig. 5 :

 5 Fig. 5: Bifurcation diagrams in the (Ω, m)-plane showing curves of the (degenerate) pitchfork PF and Hopf bifurcations H which involve the balanced state B, (a) and (d). We take ǫ = c = c b = 0.01, with the other parameters given by (15-17). The shaded regions correspond to a stable balanced state, and the subscripts for the Hopf bifurcations denote which mode is involved. Numerical simulations at the points indicated by a • are presented in Fig. 11. One-parameter diagrams in m for Ω = 1.5 are shown in (b) and (c), with norms given by α 1 and L 2 respectively, and labels for the steady states described as in Fig. 4. Panel (e) shows the one-parameter bifurcation diagram for Ω = 3.

Fig. 6 :

 6 Fig. 6: Two-parameter bifurcation diagrams, upon variation of the parameters: ǫ against Ω (a), c against Ω (b), and c b against Ω (c). The shaded regions correspond to a stable balanced state B, and Hopf bifurcations H are shown with subscripts indicating the mode of oscillation. Parameters when fixed are m = 0.025, ǫ = c = c b = 0.01, with the other values specified in the text.

Fig. 7 :

 7 Fig. 7: Diagram showing stable regions of the balanced state B upon variation of the parameters ǫ against Ω, whilst m is also varied so that m = ǫ. The vertical scale is logarithmic and c = c b = 0.01.

Fig. 8 :

 8 Fig. 8: Schematic diagram of the steady states for the case of a couple imbalance. White surrounds again indicate that both balls in the race are coincident.

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

  = 0.005, χ = 0.005, β = 1.

Fig. 9 :

 9 Fig. 9: Panel (a) is a bifurcation diagram in the (Ω, m)-plane showing curves of the pitchfork PF and Hopf bifurcations H which involve the balanced state B. Parameters when fixed are, ǫ = χ = 0.005, β = 1, and c = c b = 0.01 with the other values specified in the text. The shaded region corresponds to a stable balanced state, and the subscripts for the Hopf bifurcations denote which mode is involved. Numerical simulations at the point highlighted by a • are presented in Fig. 12. A one-parameter diagram in m for Ω = 3 is shown in (b), with norm given by Ā.

A c c e p t e d m a n u s c r i p t Ω = 3 ,

 3 while increasing the ball mass m through m c . In this diagram we use the norm Ā defined byĀ = 1 2 (x + φ y ) 2 + (yφ x ) 2 + (xφ y ) 2 + (y + φ x ) 2 .

Fig. 10 :

 10 Fig. 10: Two-parameter bifurcation diagrams, upon variation of the parameters, c versus Ω (a), and c b versus Ω (b). The shaded regions correspond to a stable balanced state B. Parameters when fixed are m = 0.025, ǫ = χ = 0.005, β = 1, and c = c b = 0.01 with the other values specified in the text.

  2, and (x, y, φ x , φ y ) = ( ẋ, ẏ, φx , φy ) = ( α1 , α2 , α3 , α4 ) = 0.

Fig. 11 :

 11 Fig. 11: Numerical simulations showing vibration levels for the static imbalance case with low damping c = c b = 0.01. Initial conditions are fixed so that the rotor begins at rest and the balls in each race are on opposite sides with (α 1 , -α 2 ) = (α 3 , -α 4 ) = π/2. Values of Ω, ǫ and m are fixed at (Ω, ǫ, m) = (2.0, 0.01, 0.025) (a), (2.5,0.01,0.005) (b), (4.0,0.01,0.005) (c), and (1.5,0.025,0.025) (d). These points are specified in parameter space by a • in Figs. 5 and 7.

Fig. 12 :

 12 Fig.12: Numerical simulations for the dynamic imbalance (ǫ, χ, β) = (0.005, 0.005, 1.0) showing vibration levels Ā above the ball positions α i . The parameter set used is the same as that in section 3.2.2 with Ω and m fixed at (Ω, m) = (3.5, 0.005). This position is also highlighted by a • in Fig.9. Initial conditions are (α 1 , -α 2 ) = (α 3 , -α 4 ) = π/2 (a), (α 1 , -α 2 ) = π/3, (α 3 , -α 4 ) = 2π/3 (b) and (α 1 , α 2 ) = 0, (α 3 , α 4 ) = π (c).

Acknowledgements

DJR gratefully acknowledges the support from a CASE award provided by the EPSRC and Rolls-Royce plc.


Alexander Fidlin

Jon Juel Thomsen

Nontrivial effects of high-frequency excitation for strongly damped mechanical systems

Keywords: high-frequency excitation, strong damping, effective properties, control, averaging

Some nontrivial effects are investigated, which can occur if strongly damped mechanical systems are subjected to strong high-frequency (HF) excitation. The main result is a theoretical prediction, supported by numerical simulation, that for such systems the (quasi-) equilibrium states can change substantially with the level of damping. For example, a strongly damped pendulum, with a hinge vibrated at high frequency along an elliptical path with horizontal or vertical axis, will line up along a line offset from the vertical; the offset vanishes for very light or very strong damping, attaining a maximum that can be substantial (depending on the strength of the HF excitation) for finite values of the damping. The analysis is focused on the differences between the classic results for weakly damped systems, and new effects for which the strong damping terms are responsible. The analysis is based on a slightly modified averaging technique, and includes an elementary example of an elliptically excited pendulum for illustration, alongside with a generalization to a broader class of strongly damped dynamical systems with HF excitation. As an application example, the nontrivial behavior of a classical optimally controlled nonlinear system is investigated, illustrating how HF excitation may cause the controller to leave the system in an unexpected equilibrium state, quite different from the setpoint. The effects can be interesting for specialists in control of mechanical systems and structures. However the obtained results are more general. Similar effects could be expected first of all for microsystems where damping forces are typically dominating over inertia forces.

Introduction

Nontrivial effects of high-frequency (HF) excitation on mechanical systems have been investigated intensively in the last decade. The monograph by Blekhman [START_REF] Blekhman | Vibrational Mechanics -Nonlinear Dynamic Effects, General Approach, Applications[END_REF] is the main reference in this area. The main effects of HF excitation in mechanical systems can be classified into three groups [START_REF] Thomsen | Some general effects of strong high-frequency excitation: stiffening, biasing, and smoothening[END_REF][START_REF] Thomsen | Slow high-frequency effects in mechanics: problems, solutions, potentials[END_REF]

:
• Stiffening is the ability of HF excitation to increase (or decrease) the effective stiffness of the mechanical system, with respect to slow motions; one of the consequences of this is the ability to stabilize or destabilize equilibria through changing the sign of the effective stiffness. • Biasing is the ability of HF excitation to attract the mechanical system to a new steady state, e.g., a quasi-equilibrium or a steady state velocity (like in all cases of vibration induced transportation). • Smoothening is the ability of HF excitation to smoothen nonlinearities, on the average; especially important is this effect for systems with dry friction, where HF dither is often used to reduce the effective hysteretic behavior with respect to slow motions.

Different mathematical approaches can be used to obtain the equations governing the slow motions of a mechanical system under HF excitation. Standard averaging [START_REF] Bogoliubov | Asymptotic methods in the theory of nonlinear oscillations, Gordon and Breach[END_REF][START_REF] Fidlin | Nonlinear Oscillations in Mechanical Engineering[END_REF] or multiple scales analysis [START_REF] Nayfeh | Perturbation Methods[END_REF] are the general ones, whereas the method of the direct separation of motions [START_REF] Blekhman | Vibrational Mechanics -Nonlinear Dynamic Effects, General Approach, Applications[END_REF] is highly specialized and effective in this particular area [START_REF] Blekhman | Selected Topics in Vibrational Mechanics, Series A on Stability, Vibration and Control of Systems[END_REF][START_REF] Thomsen | Vibrations and Stability: Advanced Theory, Analysis and Tools[END_REF]. In [START_REF] Thomsen | Slow high-frequency effects in mechanics: problems, solutions, potentials[END_REF][START_REF] Fidlin | Nonlinear Oscillations in Mechanical Engineering[END_REF] However, the main attention in these publications is paid to weakly damped systems which are usual in mechanical applications. In that case the strong HF excitation terms are balanced with large inertia forces (i.e. acceleration terms which are proportional to squared frequency). There are at least two application areas where the situation is significantly different: First, it is quite usual to apply strong damping to control systems. The reason for that is the objective to avoid unnecessary oscillations, and to lead the controlled system to a given position as soon as possible. HF vibrations are sometimes used in control systems, e.g. to smoothen the effective dry friction [START_REF] Kovaleva | Optimal Control of Mechanical Oscillations[END_REF][START_REF] Bullo | Averaging and Vibrational Control of Mechanical Systems[END_REF]. However, in the presence of strong damping terms, a strong HF excitation may prevent correct operation of the control system, e.g. the controller may leave the system at an equilibrium that is shifted away from the desired setpoint [START_REF] Fidlin | Non trivial effect of strong high-frequency excitation on a nonlinear controlled System[END_REF]. Another and perhaps even more important area is microsystems. Inertia terms are generally proportional to the mass of the system, i.e. to the cube of its linear dimension, while damping forces are usually applied to the surface of the body and thus proportional to the square of the linear dimension. Decreasing linear dimension therefore increases the importance of the damping terms.

This paper investigates nontrivial effects that may occur if a strongly damped mechanical system is subjected to strong HF excitation. The analysis is focused on the differences between the classic results well known for weakly damped systems, and new effects for which the strong damping terms are responsible. A classic example of a weakly damped pendulum subjected to elliptic vibrations of its suspension point is considered in Section 2, and compared to the same system with strong damping. These results are generalized in Section 3. Finally, as a more complicated example showing application of the general results, a classical nonlinear system under optimal control (an up-pointing pendulum on a cart) and its (mis)-behavior under HF excitation is analyzed in section 4.

Illustrative example: Pendulum with elliptic HF excitation

Consider the pendulum in Fig. 1 having length l and mass m in gravity g, whose suspension point is excited harmonically at frequency g l

ω ω = %
, and horizontal and vertical displacement amplitudes h a l and v a l , respectively, with a relative phase γ . The equation of motion is, in nondimensional form:

( )

+ + = + + && & (1)
( ) where dots denote differentiation with respect to nondimensional time t t g l = % and linear damping with coefficient β has been added. The chosen excitation form means the suspen- sion point moves along an ellipse, whose orientation is determined by a h /a v and by γ.

sin h a l t ω γ + % % g m, l θ sin v a l t ω % % ( ) sin h a l t ω γ + % % g m, l
We suppose here, following the standard analysis [START_REF] Blekhman | Vibrational Mechanics -Nonlinear Dynamic Effects, General Approach, Applications[END_REF][START_REF] Thomsen | Some general effects of strong high-frequency excitation: stiffening, biasing, and smoothening[END_REF][START_REF] Thomsen | Slow high-frequency effects in mechanics: problems, solutions, potentials[END_REF][START_REF] Bogoliubov | Asymptotic methods in the theory of nonlinear oscillations, Gordon and Breach[END_REF][START_REF] Fidlin | Nonlinear Oscillations in Mechanical Engineering[END_REF], that the excitation frequency is much higher than the pendulum natural frequency, and that the excitation displacement amplitudes are small, so that the velocity amplitudes are of order unity, i.e.:

( ) ( ) ( ) ( ) 1 1 1 1; ; 1 ; 1 . v h v v h h a O a O a A O a A O ω ω ω ω ω - - - = = = = = = (2)

Case of weak damping

First we assume that the damping is not large, i.e.

( )

1 O = β
. Introducing a fast time t τ ω = as the new independent variable, and ρ θ = & as a new dependent variable, we convert (1) to the following system of two first-order differential equations:

( ) ( )

2 sin sin sin sin cos v h A A θ ερ ρ εβρ ε θ τ θ τ γ θ ′ = ′ = - - + + + (3)
cos sin cos cos v h A A = - - + ρ σ τ θ τ γ θ. (4)
The resulting equations for θ(τ) and σ(τ) are then in the standard form for averaging [START_REF] Bogoliubov | Asymptotic methods in the theory of nonlinear oscillations, Gordon and Breach[END_REF][START_REF] Fidlin | Nonlinear Oscillations in Mechanical Engineering[END_REF].

Considering the first-order approximation, and returning to the slow time t, we obtain the following equation governing the slow motions of the elliptically excited pendulum (index 1 here and further denotes a first-order averaging approximation):

( )

2 2 1 1 1 1 1 1 1 4 2 2 sin sin2 cos cos2 v h v h A A A A θ βθ θ θ γ θ + + =- - - && & . ( 5
)
This equation explains the well known effects of HF excitation for pendulums [START_REF] Blekhman | Vibrational Mechanics -Nonlinear Dynamic Effects, General Approach, Applications[END_REF][START_REF] Thomsen | Some general effects of strong high-frequency excitation: stiffening, biasing, and smoothening[END_REF][START_REF] Thomsen | Slow high-frequency effects in mechanics: problems, solutions, potentials[END_REF][START_REF] Fidlin | Nonlinear Oscillations in Mechanical Engineering[END_REF]:

• Stabilization of the up-pointing equilibrium by pure vertical excitation (A h =0).

• Destabilization of the down-pointing equilibrium by pure horizontal excitation (A v =0).

• Biasing of the pendulum towards the main excitation axis in case of elliptic excitation

(A h A v ≠0).
The last effect can be visualized if we introduce the corresponding effective (i.e. average) potential for (5):

( )( ) 2 2 1 1 1 1 1 1 cos 1 cos 2 cos sin 2 8 4 v h v h A A A A θ θ γ θ Π = - + - - + . ( 6
)
which is depicted in Fig. 2 for different values of the phase difference γ . The minima of the potential function correspond to stable quasi-equilibria of the elliptically excited pendulum. Biasing occurs if, in the presence of HF excitation, these equilibria are shifted away from the unexcited equilibria , j j θ π = ∈Ζ. The bias is maximal for 0 γ = or γ π = . For estimating the bias offset from the down-pointing equilibrium, we can linearise [START_REF] Fidlin | Nonlinear Oscillations in Mechanical Engineering[END_REF] in the vicinity of 0 θ = , and obtain an approximate expression for the new stationary equilibrium in presence of HF excitation:

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

1 2 2 cos 2 v h st v h A A A A γ θ = -+ - . ( 7
)
Note firstly that there is no bias if 2 γ π = ± , i.e. if the excitation ellipse has a horizontal or vertical axis, and secondly that the stationary bias does not depend on the damping coefficient β.

Case of strong damping

Again we investigate the system (1), but assuming instead the damping to be strong, i.e.

( )

2 1 , 1 k k O β ω = >> = , (8)
so that instead of (3) one has:

( )

sin sin sin sin cos . v h k A A θ ερ ρ ρ ε θ τ θ τ γ θ ′ = ′ = -- + + + (9)
Now we eliminate the non-small oscillating terms, though the strong damping term will be still present. This can be achieved by the following transformation, from ρ to a new variable q, based on the particular solution of the unperturbed system (i.e. [START_REF] Kovaleva | Optimal Control of Mechanical Oscillations[END_REF] with 0

ε = ): ( ) ( ) ( ) ( ) ( ) ( ) 2 2
sin cos , where sin sin cos cos

1 sin sin cos cos . 1 s c v h s v h c q F F kA A k F k A A k F k ρ θ τ θ τ θ γ γ θ θ θ γ γ θ θ = + + + + = + - - = - + (10)
Applying this transformation to [START_REF] Kovaleva | Optimal Control of Mechanical Oscillations[END_REF]

s c s c s c v h s v h c q F F q kq F F q F F kA A k F k A A k F k γ θ θ ′ = + + ′ ′ ′ = -- - + + + - + ′ = + + - ′ = - + (11)
where primes on F s and F c indicate derivatives with respect to θ. The structure of these equations is typical for systems with strong damping. The procedures for the corresponding asymptotic analysis (including mathematical proof and error estimation) can be found in [START_REF] Fidlin | Nonlinear Oscillations in Mechanical Engineering[END_REF][START_REF] Volosov | Method of Averaging in the Theory of Nonlinear Oscillating systems[END_REF]. The first-order approximation (direct averaging of the right hand sides of the two first equations in [START_REF] Fidlin | Non trivial effect of strong high-frequency excitation on a nonlinear controlled System[END_REF]) is:

( ) ( ) ( ) ( ) ( ) 1 1 1 0 1 sin , 2 s s c c q kq F F F F θ ε θ ε θ θ θ θ ′ = ′ ′ ′ = -- - + (12)
which, however, gives insufficient information about the system's dynamics: It fails in predicting any evolution of θ, and demonstrates only the dominance of the damping term. Thus the second-order approximation has to be considered. The appropriate general results for this can be formulated in the following theorem [START_REF] Fidlin | Nonlinear Oscillations in Mechanical Engineering[END_REF]: Consider a system with a slow master variable x and a strongly damped slave variable y:

( ) ( ) ( ) ( ) 0 0 , , , , 0 ; 0 x
X x y t y ky Y x y t x x y y ε ε ′ = ′ = -+ = = (13)
where k > 0 and 0 < ε 1, alongside with the following second-order approximation:

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 1 2 0 0 1 2 1 0 0 0 , 0, 0 ; 0 ; , 0, 1 , 0, ; , 0, , t y y k Y t x y X t X X u Y t u X d x k y ξ ε ξ ε ξ η η ε ξ ξ η ξ ξ ξ ξ ξ τ ξ τ = = ′ = Ξ + Ξ ′ = -+ = = Ξ = ∂ ∂ Ξ = + = -Ξ ∂ ∂ ∫ (14)
where denotes averaging with respect to time t. Under the usual requirements [START_REF] Bogoliubov | Asymptotic methods in the theory of nonlinear oscillations, Gordon and Breach[END_REF][START_REF] Fidlin | Nonlinear Oscillations in Mechanical Engineering[END_REF] on continuity of the right hand sides in [START_REF] Kwakernaak | Linear Optimal Control Systems[END_REF], the solutions of the systems ( 13) and ( 14) are close, with a residual error O(ε 2 ) on an asymptotically long time interval after a short transient process at the very beginning, i.e.:

( ) ( )

2 2 1 1 , , . ln x O y O O t O ξ ε η ε ε ε ⎛ ⎞ ⎛ ⎞ -= -= < < ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ ⎝ ⎠ (15)
Applying this theorem to the transformed pendulum system [START_REF] Fidlin | Non trivial effect of strong high-frequency excitation on a nonlinear controlled System[END_REF], we calculate the necessary averages:

A c c e p t e d m

a n u s c r i p t

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 1 0 0 2 2 2 2 2 2
sin cos 0; cos sin ,0, sin sin cos sin cos sin cos ; 1

sin cos 1 1 1 1 sin sin 2 cos 2 , 2 1 2 2 1 1 s c s c s c s c s c q q v h v h v h X q F F u F F Y F F q F F X X F F q A A A A A A k k k k k θ τ θ τ θ θ θ = = = + + Ξ = = - + ′ ′ = - - + + + ∂ ∂ ′ ′ = + = ∂ ∂ ⎛ ⎞ - Ξ = - - - - ⎜ ⎟ + + + ⎝ ⎠ (16)
and obtain the second-order approximation:

( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 sin cos 1 1 1 sin sin 2 cos 2 2 1 4 1 2 1 sin sin cos sin cos . v h v h v h s c s c A A A A A A k k k k k k q kq F F F F γ γ θ ε θ θ θ ε θ ε θ τ θ τ θ τ θ τ ⎛ ⎞ - ′ ⎜ ⎟ = - + + - ⎜ ⎟ + + + ⎝ ⎠ ′ ′ ′ = - - - + + (17)
The first equation is the important one, describing the slow dynamics of the strongly damped pendulum by a first-order differential equation. To estimate the bias of the downpointing equilibrium, we linearise (17) near 0 θ = , and obtain an approximate expression for the new equilibrium in the presence of HF excitation:

( ) ( ) 2 2 2 2 sin cos 2 1 v h st v h A A k k A A γ γ θ - = + + - . (18)
which approaches the solution [START_REF] Blekhman | Selected Topics in Vibrational Mechanics, Series A on Stability, Vibration and Control of Systems[END_REF] for the weakly damped pendulum as the damping coefficient k→0.

The structure of the solutions ( 7) and (18) are similar: both vertical and horizontal excitation must be present to bias the pendulum. However, there are at least two significant differences. Firstly, the variation of the stationary bias θ 2nd with the phase γ is more complex now, with the phase difference generating the strongest bias (largest |θ 2nd |) being dependent on

A v = 1, A h = 1, ω =50, γ = -π/2.)

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

damping. Secondly, the bias itself depends strongly on the damping coefficient: For the weakly damped pendulum, there is no bias in the case of elliptical excitation with horizontal or vertical major axes (γ = ±π/2). The strongly damped pendulum, by contrast, is biased in that case, with a maximum bias for damping 1 k = . This is illustrated in Fig. 3, displaying the approximate prediction (18) for parameters as given in the legend; as appears the results agree with results of numerical simulation of the full system (1). Fig. 4 shows how the bias depends on the phase difference for the same parameter values and the maximizing value 1 k = ; as appears phase control can be used to shift the biased equilibrium continuously through a range of negative and positive angles.

Note that (18) is based on linearization, assuming small θ 2s . Thus, the seeming possibility of infinite bias as the denominator approaches zero is spurious; the solution of the full system (17) remains limited.

Returning to dimensional form we can rewrite (18) as follows:

( )

2 2 2 2 2 2 2 2 sin cos 4 2 1 v h st v h g a a l l g a a l g l β ω γ γ ω θ β ω ω ⎛ ⎞ - ⎜ ⎟ ⎝ ⎠ = ⎛ ⎞ - + + ⎜ ⎟ ⎝ ⎠ % % % % (19)
This form of the solution demonstrates explicitly that if the vibration acceleration is dominant compared to gravity, i.e. 2 g l ω >> % , then the stationary tilt of the damped pendulum depends on the phase difference between the vertical and the horizontal excitations in the same way as it does in case of the weak damping (cf. ( 7)).

Generalization

Consider the following system:

( ) ( ) ( ) ( ) ( ) [ ] ( ) 1 2 , , , 1 diag ( ) ( ) ( ) , 0, 1, , . m i k k k k i m ε ε τ ε τ ε ′ = + ′ = - + = << = > = x h y X x y Ky G x h 0 0 K x x x x x L K (20)
where the vector functions h, X, G and the diagonal damping matrix K are bounded and Lipschitz-continuous in x and y, and X and G are 2π-periodic in the fast time τ.

Applying the same approach as above we obtain the second order approximation:

( ) ( ) ( ) ( ) 2 1 2 , ε ε ε τ ′ = + ′ = - + K G ξ Ξ ξ Ξ ξ η ξ η η (21)
Note that the steady state solution for y in (20) is εK -1 G(x,τ), and thus the first-order contri- bution is:

( )

1 1 ( ) ( , ) ( , ) ε τ τ
- = + h K G X Ξ ξ ξ ξ . ( 22
)
Here the first term is small, and if the average ( , ) 0 τ = X ξ (as in Section 2.2), the firstorder contribution does not give us any additional information. The second-order contribution can be calculated to (cf. ( 14)):

( ) ( ) ( ) 0 1 2 0 , , where , . d τ
τ τ θ θ - = = ∂ ∂ = + = ∂ ∂ ∫ x y X h u K G u X x y ξ Ξ ξ ξ ξ (23)
As ε→0, then ξ→x+O(ε 2 ) and η→y+O(ε 2 ) on a time scale O(1/lnε)<τ<Ο(1/ε) (cf. ( 15)).

Equations ( 20)-( 23) generalize the approach of calculating approximate equations of slow motion, for systems with several strongly damped state variables.

Application example: Failure of standard optimal control under HF excitation

The problem

Here we employ the general results (20)-(23) to illustrate how a classical control system, which is inherently strongly damped, in the presence of HF excitation may behave quite differently than expected.

Fig. 5 shows the system, which is often used to test or illustrate control algorithms [START_REF] Kwakernaak | Linear Optimal Control Systems[END_REF]. It contains a cart used to balance a pendulum in the up-pointing position, against gravity. The system state is defined by the pendulum angle θ and the cart position S, as observed from the rigid platform. The cart has mass M and linear damping coefficient d, and the pendulum has mass m and rotational inertia J about its center of gravity at distance L from the loss-free hinge. The rigid platform can be vibrated relative to a fixed inertial frame, with horizontal and vertical amplitude h a and , where non-dimensional variables and parameters are given by (cf. Fig. 5):

( )

( ) 0 0 2 2 0 0 0 ; ; 1 ; ; ; ; ; 1 ; . 2 h v h v g J U t t l L u l m L M m l a a S m L d s a a l l l M m l M ω ω ω ω ω α β ω Ω ⎛ ⎞ = = = + = = ⎜ ⎟ + ⎝ ⎠ = = = = < = + % (25)
The control u in ( 24) is intended to keep the system at zero state,

{ } T s s θ θ = = z 0 & &
, and return it to zero in case of any perturbations due to initial conditions or external forces. It is designed as an optimal linear regulator, basing its actions on a linearized model of the unforced system to be controlled, and on measurements -assumed here to be error-and noisefree -of instantaneous values of the state z i.e., u = u(z). The linearized model of the unforced system can be obtained from (24) by letting 0 v h a a = = , and linearising for small values of θ and θ & , to obtain a system of linear first-order equations in state space form:

, where

0 1 0 0 0 2 1 0 0 1 1 1 ; . 0 0 0 1 0 2 1 1 0 0 1 1 1 u β α α α α β α α α = + ⎡ ⎤ ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ - - ⎢ ⎥ ⎢ ⎥ - - - ⎢ ⎥ ⎢ ⎥ = = ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ - ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ - - - ⎣ ⎦ ⎣ ⎦ z Az B A B & (26)
The control law ( ) u z for the optimal linear state-feedback controller is:

ˆT u = -k z . ( 27
)
where the optimal gain vector k is chosen so as to minimize a weighted sum of integrated control error and control effort: The minimizing solution for k can be found in many texts on control theory (e.g. [START_REF] Kwakernaak | Linear Optimal Control Systems[END_REF]):

( ) 2 0 ( ) T Ru dt ∞ + ∫ z Qz z , ( 28
) U S θ sin A t Ω % ( ) sin t γ Ω + % U S θ sin A t Ω % ( ) sin t γ Ω + % ( ) sin h a t γ Ω + % ( ) sin v a t Ω % g M m L U S θ sin A t Ω % ( ) sin t γ Ω + % U S θ sin A t Ω % ( ) sin t γ Ω + % ( ) sin h a t γ Ω + % ( ) sin v a t Ω % U S θ sin A t Ω % ( ) sin t γ Ω + % U S θ sin A t Ω % ( ) sin t γ Ω + % ( )
1 ˆT R - = k B S, (29)
where S is the solution of the algebraic Riccati equation:

1 0 T T R - + - + = A S SA SBB S Q . ( 30
)
For given values of (A,B,Q,R), the Riccati matrix S is readily computed using, e.g., the builtfunction LQR in MATLAB [START_REF] Biran | MATLAB 6 for Engineers[END_REF] .

Thus, for a given mass ratio α and damping β of the physical system, A and B is calculated using (26), and with chosen Q and R, the Riccati matrix is given by the solution to (30). Then k is calculated using (29), and the control force u can be calculated using (27) for any state z, i.e. in explicit form:

1 2 3 4 ˆˆˆû k s k s k k θ θ = -- - -& & . ( 31
)
where the control coefficients 1 4 kare typically large, compared to the parameters of the uncontrolled part of the system (26).

Near the setpoint z=0, the controlled system is governed by ( 26) with (27) inserted, i.e.

= z Az

&

where ˆT = -A A Bk . The optimal control (29) with (30) then guaranties that all eigenvalues of A have negative real parts, so that the setpoint is asymptotically stable to small perturbations.

System response with strong HF excitation: Numerical simulation

The question is how the controlled system behave in the presence of small but rapid oscillations of the supporting platform, i.e. if 0 1 , 0 1 ,

h v a a < <
and ω 1 in (24). It was already mentioned in Section 2 how elliptic HF excitation may cause equilibrium bias, shifting the equilibrium(s) of the pendulum away from being straight down (or up). The control system, however, might be able to compensate for this just by using more power, so that the performance of the controller could be expected to be only slightly affected by the HFexcitation. The setpoint z = 0 would still remain stable, though with a small overlay of HF oscillations corresponding to the excitation. Or another stable equilibrium point would appear, deviating less from z = 0 for larger allowable control efforts. Numerical simulations contradict these expectations. As appears from Fig. 6, the controller works perfectly if there is either no (Fig. 6a), or only horizontal (Fig. 6b), or only vertical (Fig. 6c) HF-excitation, in all cases stabilizing the pendulum in the up-pointing position, and moving the cart to the zero point. The pendulum cart motions are overlaid by small HFoscillations, but on the average the control system fulfills its objective. However, as appears from Fig. 7, it fails when vertical and horizontal oscillations are applied simultaneously. The control system still stabilizes the cart and the pendulum, but the final average cart position is not even close to zero, and the pendulum is stabilized with some permanent tilt, i.e. both state variables are biased (Fig. 7a). An attempt to reduce the permanent error through increased control effort does not help, the stronger controller just amplifies the misbehavior (Fig. 7b). Next we investigate this effect theoretically.

A c c e p t e d m a n u s c r i p t 11 4.3. Asymptotic analysis

Here we derive and examine an approximate analytic prediction for the biased equilibrium state. As for the classical analysis in [START_REF] Kwakernaak | Linear Optimal Control Systems[END_REF] (for the system without HF excitation) we let α = 0, which corresponds to neglecting the horizontal reaction force between the cart and the pendulum. This restriction, assuming a very light or short pendulum, greatly simplifies the analysis, without sacrificing essential features of the system. We also neglect mechanical damping β, as compared to the dominating damping coefficient 2 k of the control system. Making allowance for domination of the control terms we assume the corresponding coefficients to be of the same magnitude order as the high excitation frequency, i.e.

( )

ˆ, 1 O ω = = k k k
. Taking (24) along with the control force (31) into account we obtain the equations of motion:

( ) ( ) ( ) 2 1 2 3 4 2 1 2 3 4 sin cos 1 sin sin . h v s ks k s k k a t k s k s k k a t ω ω ω θ ω θ ω ω γ θ ω θ θ θ ω ω θ = - - - - + + = + + + + -
In [START_REF] Fidlin | Nonlinear Oscillations in Mechanical Engineering[END_REF] the linear approximation of this system was investigated, while here we apply asymptotic analysis to the full nonlinear system.

To transform (32) to a form suitable for averaging, we first introduce a fast time τ = ωt and a small parameter ε=ω -1 1, and add new dependent variables p(t) and q(t) for converting to a system of first-order equations:

k p k q A k s k q k p k q A k s k ε θ ε τ γ ε θ θ θ τ ε θ θ θ ′ = ′ = ′ = - - + + - + ′ = + - + + + (33) where ( ) , (1), (1)
,

v v h h d d A a O A a O τ ω ω ′ = = = = =
and the variables s and θ are seen to be slowly varying in τ. Next we eliminate the large oscillating terms in the equations for p and q, by introducing a third slow variable, next to s and θ. This can be achieved by calculating first the particular solution of the unperturbed system ((33) with ε = 0), which gives that θ and s are constants, and a third equation (obtained by computing cos p q θ ′ ′ + for ε=0, which can be integrated with respect to τ to give:

( )

cos cos cos sin cos h v x p q A A θ θ τ γ θ τ = + + + - . ( 34
)
where x is a constant of integration. Then consider instead (34) as the definition of a new dependent variable x(τ). Eliminating q, the corresponding transformed form of (33) becomes:

(

cos cos cos sin cos sin sin sin cos cos cos cos cos cos sin cos sin .

h v h v h v h s p x p A A x p A A p k k p k x A A θ τ τ γ ε θ ′ = ′ = - - + + ′ ′ = - + + + ′ = - - - - + + + + - + (35)
The first three equations now have the standard form, i.e. their right-hand sides are small, but the last one contains non-small oscillating terms (those with A h , A v ), as well as non-small linear terms (those with k 2 and k 4 ) expressing strong damping due to optimal control. The oscillating terms can be eliminated if we apply the following transformation based on the particular periodic solution to the unperturbed system:

)

1
h v h v p r k x k k A k k k k A k k A k k k k k k A k k k k λ θ λ θ τ λ θ τ τ θ θ ε ′ = + + + ′ = - + + - + + + ′ ′ = - + + ′ - + + + ′ = - - + (37)
Since (37) has the general form (20), the general solution (21)-(23) applies, provided k 2 -k 4 is positive. The latter will normally hold true, since k 2 -k 4 >0 is the necessary condition for stability of the control algorithm in the linear case. The first-order approximation, obtained by retaining only terms to the lowest order of ε in each equation in (21) becomes:

(

)

1 1 1 1 1 1 1 1 1 1 1 1 1 3 2 2 2 2 1 3 2 3 2 1 2 4 1 1 cos 1 sin cos cos 2 cos 2 1 1 1 sin 2 cos sin 2 2 2 cos . h v v h v h h s x x x AA A A A A A r k k r ελ θ ε λ θ ε θ ελ λ θ ε θ γ λ ε θ λ λ λ γ λ γ θ ′ = ′ = - ′ = - - + + ⎛ ⎞ + - + + + - ⎜ ⎟ ⎝ ⎠ ′ = - - (38)
This approximation determines the stationary tilt of the pendulum, however it is insufficient for finding the position of the cart. As a matter of fact, the corresponding stationary solutions can be readily obtained by equating to zero the left-hand sides, and solving the algebraic equations for (s 1 ,θ 1 ,x 1 ,r 1 ). Substituting the expression for λ 1 from (36), one finds that the first, second, and fourth equation just yield trivial stationary solutions x 1 =0, r 1 =0, while s 1 ∈R (to first order of approximation). Thus the important equation is the third one. Substituting into this equation λ 2 and λ 3 from (36), linearising near θ 1 = 0, and letting 1 0 x′ = , we obtain an approximate prediction for the stationary tilt θ 1 = θ 1st of the pendulum:

h v st v h v a a k k k L J g mL a g G ak ak k k J L mL J G g g k k L mL γ γ θ ⎛ ⎞ ⎜ ⎟ ⎜ ⎟ Ω - - Ω ⎜ ⎟ ⎛ ⎞ ⎜ ⎟ + ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ ⎝ ⎠ = Ω - +Ω - - ⎛ ⎞ + ⎜ ⎟ ⎝ ⎠ ⎛ ⎞ = + - + ⎜ ⎟ ⎝ ⎠ (41)
Here G holds the terms that are independent of vibration amplitudes. In the important case of small gravity we obtain the following approximation:

( ) ( ) 2 4 2 1 2 2 2 2 4 ˆˆcos , ˆˆĥ v st h v a a k k g L a k a k k γ θ - = - << Ω - - (42)
This form of the solution demonstrates that if the vibration acceleration is dominant compared to gravity, i.e. 2 g l ω >> % , then the stationary behavior of the damped pendulum is similar to that discussed in section 2.2.

Comments to the theoretically predicted equilibrium state

The prediction (39)-(40) of the stationary state is identical to what was obtained in [START_REF] Thomsen | Slow high-frequency effects in mechanics: problems, solutions, potentials[END_REF] by the method of direct separation of motions. However, it differs from the result in [START_REF] Fidlin | Nonlinear Oscillations in Mechanical Engineering[END_REF][START_REF] Fidlin | Non trivial effect of strong high-frequency excitation on a nonlinear controlled System[END_REF], which was obtained for a linearization of a system similar to (32), and lacks the middle term in the denominator of (39). This deviation is yet a demonstration, beyond what was observed already in e.g. [START_REF] Thomsen | Some general effects of strong high-frequency excitation: stiffening, biasing, and smoothening[END_REF][START_REF] Thomsen | Vibrations and Stability: Advanced Theory, Analysis and Tools[END_REF], of the significance of the choice of linearization stage when calculating effects of HF-excitation for nonlinear systems: Results obtained for a linearized system may differ from linearized results of a nonlinear analysis.

The structure of the approximate prediction (39) is similar to that of a simple strongly damped pendulum (18). As appears from (39)-(40), a nonzero value for 1st θ and 2st s can occur only when 0 h v A A ≠ , i.e. when the HF excitation is not purely horizontal or vertical.

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

Further, (40) reveals that if 3 1 0 k k > then 2 0 st st s θ < , that is: with the controller adding positive (or negative) stiffness to both the cart and the pendulum part of the system, the cart and the pendulum will be biased in opposite directions. This is a consequence of the function (28) to be minimized by the optimal controller, including the choice of weighting constants Q and R. It is unrelated to the HF excitation, but merely reflects that at equilibrium ( 0)

s θ = = & &
, the function (28) will be minimal for u(z) 2 =0, which implies (k 1 s+k 3 θ) 2 = 0 and thus (40).

It also appears from (39)-(40) that the phase difference γ between vertical and horizontal excitation can be used to control the equilibrium bias. At a certain phase difference, , the up-pointing, unbiased pendulum equilibrium θ = 0 even becomes stable, as if the system were not HF-excited at all. Finally we note from (39) that the main source of bias is the part of the controller which damps oscillations of the cart: If 2 0 k = , i.e. the controller doesn't damp cart motions (cf. (31)), then the right-hand side of (39) vanishes, and 1 0 st θ = becomes a stationary, unbiased equilibrium.

Comparison with numerical simulation

Figs. 8 and 9 compare results of the approximate analytical prediction (39), (40) to results of numerical simulation of the original equations (32). MATLAB's built-in function ODE45 [START_REF] Biran | MATLAB 6 for Engineers[END_REF] was used for the numerical simulation, though other integrators were tested with similar results. The simulations were run until stationary conditions were attained, whereafter average values and oscillation amplitudes were calculated numerically from the stationary part of the time series.

Fig. 8 shows a prediction given by ( 39), (40) of how the stationary average value of the π , where the HF excitation translates the platform along a circular pattern (for a h =a v ).

The bias diminishes when the phase difference is near integer multiples of π where the platform translates along straight lines. Fig. 9 shows how the bias increases with increasing vertical HF excitation amplitude. With increased bias, it appears, the agreement between analytical prediction and numerical simulation becomes poorer, reflecting that (39) is valid only for small values of the pendulum average angle. Actually (39) predicts that the bias approaches infinity, reverses sign, and then drops off for even higher values of v A . For the "real" system, however, numerical simulation shows that the bias just increases monotonically.

Comments on the chosen asymptotic approach

The chosen asymptotic approach, based on the modified averaging procedure for systems with several strongly damped state variables, is definitely not the only one suitable for the considered problem. The same results can be obtained, for example, by the method of direct separation of motions [START_REF] Blekhman | Vibrational Mechanics -Nonlinear Dynamic Effects, General Approach, Applications[END_REF][START_REF] Thomsen | Slow high-frequency effects in mechanics: problems, solutions, potentials[END_REF][START_REF] Fidlin | Nonlinear Oscillations in Mechanical Engineering[END_REF]. However, it is supposed that for strongly damped systems it is important to distinguish between master and slave variables (cf. ( 13)), which was underlined by the chosen asymptotic procedure.

Conclusions

According to the theoretical predictions of this work, strong high-frequency excitation may change the equilibria of strongly damped systems in a manner that differs substantially from that which is well known for slightly damped systems. The analysis of this has been illustrated above for a simple example, using a slightly modified averaging approach. This has been generalized to a broader class of systems, and exemplified for a classical control problem. Numerical simulations were shown to support the theoretical predictions, and thus a natural next step is to test the predictions with a simple laboratory setup.

The effects of HF excitation for strongly damped systems, investigated in this paper, may be important at least in two areas: Control science, and microsystems.

In control science HF excitation has been known and used for a long time in the form of the so-called dither signals, for example, for smearing out discontinuities such as may be associated with dry friction and sticking behavior [START_REF] Meerkov | Principle of vibrational control: theory and applications[END_REF]. However, with strong damping terms, quite usual in linear feedback control system in general, the HF excitation can cause unexpected bias, like in the example investigated in Section 4. The problem can be also solved in the framework of the control science for example by the appropriate filtering of the velocity signal (cf. [START_REF] Fidlin | Nonlinear Oscillations in Mechanical Engineering[END_REF]).

For microsystems the effects of HF excitation is relevant, and maybe important, because of the increasing role of damping with downsizing. The authors expect that HF effects may be a key for explaining several experimental observations of effects of microwave radiation on large molecules (like proteins) and nano-structures in solvents [START_REF] Copty | Evidence for a specific microwave radiation effect on the green fluorescent protein[END_REF]. The corresponding model equations are much more complex than for the present study, and typically contain integral terms describing long-term interaction of the micro-object with the fluid. The investigation of such effects, and the development of the appropriate mathematical approaches, could be an important direction for future work.

  these

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

 2 are briefly compared, and shown to give identical or very similar results.

Fig. 1 .

 1 Fig. 1. Model of the excited pendulum.

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

 3

  small terms in (3) can be eliminated if we apply the following transformation (based on the solution to the unperturbed problem, i.e. ε =0 in (3)):

Fig. 3 .

 3 Fig. 3. Variation of stationary pendulum equilibrium bias θ 2nd with damping k, for a maximizing value of phase γ. (Parameters:A v = 1, A h = 1, ω =50, γ = -π/2.)

Fig. 4 .

 4 Fig. 4. Variation of pendulum bias θ st with phase difference γ. (Parameters as for Fig. 3, and k = 1.)

va

  respectively, and frequency Ω. The system reaction to perturbations is governed by a feedback control force, ( , , , ),

Fig. 5 .

 5 Fig. 5. Inverted pendulum, balanced by a moving cart.

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

 10 Q and R are positive definite weighting factors. These factors are chosen to reflect a suitable trade-off between the conflicting requirements of quickly reducing the setpoint error to a small value (large Q/R), without expending excessively large control power (small Q/R).

Fig. 6 .

 6 Fig.6. Results of numerical simulation of the cart and the pendulum system (24) under optimal control (31), when there is (a) no HF-excitation (a h =a v =0); (b) only horizontal HF-excitation (a h =0.02); (c) only vertical HF-excitation (a v =0.02). Solid line: pendulum tilt θ(t); Dotted line: cart position s(t). Parameters: α=0, β=0, ω=50, γ=π /2, R=1, Q = [ [10 0 0 0] T 0 [0 0 10 0] T 0] ], ˆT k ={3.1623, 9.1310, 14.1829, 13.8258}, z(0) T ={0, 0, 0, 0.5}.

Fig. 7 .

 7 Fig. 7. Results of numerical simulation of the cart and the pendulum system (24) under optimal control (31), when there is simultaneous horizontal and vertical HF-excitation (a h a v ≠0), and (a) controller gain ˆT k ={3.1623, 9.1310, 14.1829, 13.8258}; (b) stronger (twice the one for (a)) controller gain ˆT k ={6.3246, 18.2620, 28.3658, 27.6516}. Solid line: pendulum tilt θ(t); Dotted line: cart position s(t). Parameters (unless otherwise stated): as for Fig. 6.

θFig. 8 .

 8 Fig.8. Stationary average pendulum angle θ 1st and cart position s 2nd of the optimally controlled inverted pendulum, as a function of the phase difference γ between horizontal and vertical HF excitation. Parameters: α=β=0, R=1, Q=[ [10 0 0 0] T 0 [0 0 10 0] T 0] ], k = {3.1623; 9.1310; 14.1829; 13.8258} T , ω=50, a h = a v =0.02.

A c c e p t e d m a n u s c r i p t 2 π or 3 2

 22 pendulum angle θ (solid line) and cart position s (dashed) varies with the phase difference γ of the HF excitation. The agreement between the analytical predictions and numerical simulation (in marker symbols) is seen to be very good. As appears a substantial bias -an offset from the control setpoint at zero -occurs when the phase differences is near

Fig. 9 .

 9 Fig. 9. Stationary average pendulum angle θ 1st and cart position s 2nd of the optimally controlled inverted pendulum, as a function of the HF excitation amplitude ratio A v /A h , with A h kept constant. Parameters: γ=π/2, and otherwise as for Fig. 8.

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

	( ) sin	sin	(	( ) ( ) sin cos	( )	cos	) ( ) (	sin	( )	) cos ,
	( ) ( )	cos cos	( 1 sin ( sin 2 2 1	) cos sin ) cos sin	,
					we obtain:
										γ=0 γ=0
										γ=π/4 γ=π/4
										γ=π/2 γ=π/2
	Π Π
	2π -2π -					0 0				2π 2π
						θ θ
	Fig. 2. Potential function for the slow motions of the elliptically excited pendulum.

& &imposed horizontally to the cart. Motions of this system are governed by:

& && & &(32)

Substituting this into (35), the original system (32) has been transformed into its final form:

As mentioned, the first-order approximation (38) does not give information about the stationary value of cart position s 1 . The second-order approximation is needed for that, as for the uncontrolled pendulum in Section 2.2. This is obtained by including all terms in (21), and calculating the corresponding stationary value of s 1 . Omitting for brevity these elaborate calculations, we state here the final and quite simple result: