G W Blankenship

A Kahraman

STEADY STATE FORCED RESPONSE OF A MECHANICAL OSCILLATOR WITH COMBINED PARAMETRIC EXCITATION AND CLEARANCE TYPE NON-LINEARITY

A mechanical system exhibiting combined parametric excitation and clearance type non-linearity is examined analytically and experimentally in an effort to explain complex behavior that is commonly observed in the steady state forced response of rotating machines. The specific case of a preloaded mechanical oscillator having a periodically time-varying stiffness function and subject to a symmetric backlash condition is considered. A generalized solution methodology is proposed based on the harmonic balance method. The resulting non-linear algebraic equations are solved by using a direct Newton-Raphson technique, in which a closed form Jacobian matrix is computed using frequency domain methods. Analytical solutions are validated by comparison with numerical integration results and experimental measurements obtained from a gear dynamics test rig.

INTRODUCTION

In the study of rotating machinery dynamics, of primary interest is the prediction and measurement of the steady state forced response as excited by internal system components such as gears, rolling element bearings, turbines, pumps, sprags, clutches and the like. Such systems are inherently non-linear with temporally and spatially varying system parameters which arise due to periodically changing contact regimes and separations due to clearances or backlash. For instance, in a spur gear pair, the gear mesh stiffness varies periodically with the angular positions of the gears as the number of tooth pairs in contact alternates between n and n+1, where n is the integer part of the profile contact ratio. Similarly, the radial stiffness of a ball or roller bearing is time-varying due to fluctuations in the number of rolling elements in contact as the shafts rotate. In both examples, a clearance is introduced into the system in the form of gear backlash or a radial gap to ease the assembly and ensure proper operation.

In the course of the authors' experimental investigations with mechanical power transmission systems, several phenomena have been observed in the forced response of such systems which are yet to be explained analytically; these include but are not limited to multiple coexisting limit cycles, modulation sidebands in narrow band spectra, sub-and superharmonic motions and aperiodic or apparent chaotic behavior. In an attempt to develop a fundamental understanding of the physical mechanisms which give rise to such complex behavior, a governing equation is introduced based on the simple mechanical oscillator of unit mass shown in Figure 1:

x¨+2zx˙+c[u(t)]g[x(t)]=f(t), (1a)
c(u)=w 1 + s K k=1
[w 2k cos (ku)+w 2k+1 sin (ku)], (1b)

g[x(t)]= 8

x(t)-1, 0, x(t)+1,

x(t)q1, =x(t)=E1, x(t)Q-1.

(1c)

Here, ( •) denotes differentiation with respect to time t, z is an equivalent system damping ratio, f(t) is an external forcing function, c[u(t)]=c[u(t)+2p] is a spatially periodic stiffness function and g[x(t)] is a non-linear restoring function. In general, c(u) is assumed to be positive definite with w 1 =1 and can therefore be expressed in terms of the Fourier series given by equation (1b). The restoring function g[x(t)] may describe softening or hardening characteristics as well as clearance type behavior. In this study, g[x(t)] is assumed to be a symmetric, piecewise-linear clearance or backlash function as defined by equation (1c). The external forcing function f(t) may take on many different forms depending upon the system configuration, but is usually assumed to be periodic in u and can, therefore, be represented as the Fourier series

f(u)=f 1 + s L l=1 [f 2l cos (lu)+f 2l+1 sin (lu)]. (1d)
The system rotation angle u(t) is given by u(t)=Lt+bx(t), where L is a time-invariant mean angular velocity. The term bx(t) represents a linear variation in u(t) with x(t).

Hence, c(u) and f(u) are both exponentially modulated by x(t) and b is the so-called angle modulation depth [1].

1.1. literature review Previous studies in the literature have focused almost entirely on two limiting cases of equation (1). By far the most studied form of equation ( 1) is the linear system obtained by neglecting the clearance non-linearity such that g[x(t)]=x(t). This leads to a form of the Hill or Mathieu equation [START_REF] Richards | Analysis of Periodically Time-varying Systems[END_REF]. In most studies, sinusoidal or piecewise linear variations of the system parameters are assumed and, typically, only the initial value problem is solved with an emphasis on dynamic stability and parametric behavior of the system [START_REF] Richards | Analysis of Periodically Time-varying Systems[END_REF]. A general solution methodology for the steady state forced response of systems subject to combined parametric and forcing excitations was presented by Hsu and Cheng [START_REF] Hsu | [END_REF] in terms of the fundamental matrix of the homogeneous system. The harmonic balance approach has also been employed by Blankenship and Singh [1] to explain the presence of modulation sidebands in machinery spectra.

The other limiting case of equation ( 1) considers a piecewise linear oscillator with external excitation and time-invariant stiffness such that c[u(t)]=1. Such systems have been shown to exhibit an entire spectrum of nonlinear phenomena such as the dependence of steady state solutions to initial conditions, multiple coexisting solutions as well as subharmonic and chaotic motions. Due to the strong non-linearity introduced by g[x(t)], standard perturbation solution techniques are not suitable [4]. The viable solution techniques may be divided into four categories. Piecewise linear solution techniques, such as those employed by Shaw and Holmes [5] and Natsiavas [START_REF] Natsiavas | [END_REF], seek closed form linear solutions corresponding to each segment of g[x(t)]. The resulting solutions must be precisely matched at the instants when impact occurs. This procedure may be quite difficult to implement, especially with impacts that are not equally spaced. Single and multiple term harmonic balance (HBM) approaches have been employed effectively to obtain the steady state forced response by several investigators including Tomlinson and Lam [START_REF] Tomlinson | Frequency response characteristics of structures with single and multiple clearance-type nonlinearity[END_REF], Kahraman and Singh [4], Pfeiffer and Fritzer [START_REF] Pfeiffer | [END_REF], Choi and Noah [9], Kim and Noah [10] and Lau and Zhang [11]. One limitation of HBM is that the effects of initial conditions on the solution cannot be considered. A third method introduced by Pfeiffer and Kunert [START_REF] Pfeiffer | Rattling models from deterministic to stochastic processes[END_REF] uses a stochastic model applied to the well-known Fokker-Planck equation. Finally, direct numerical integration has been widely employed to study systems with clearance type non-linearities; examples include the investigations by Mahfouz and Badrakhan [START_REF] Mahfouz | Chaotic behavior of some piecewise-linear systems, part II: systems with clearances[END_REF] and Kahraman [START_REF] Kahraman | [END_REF]. Although very valuable as a tool to validate other analytical methods, numerical integration is not practical due to the computational expense required to obtain steady state solutions, especially for lightly damped systems.

The main focus of this study is the steady state forced response of preloaded systems exhibiting combined parametric excitation and clearance type nonlinearities. Although one limiting version of equation ( 1) having an external alternating forcing term was solved numerically by Kahraman and Singh [START_REF] Kahraman | Interactions between time-varying mesh stiffness and clearance non-linearities in a gear system[END_REF] and Sato et al. [START_REF] Sato | Bifurcation sets and chaotic states of a gear system subjected to harmonic excitation[END_REF], a generalized analytical treatment of equation ( 1) is yet to be developed. Since steady state solutions are of interest in this study, a multiple term HBM approach is employed. The resulting nonlinear algebraic equations are solved by using a direct Newton-Raphson technique where a closed form Jacobian matrix is computed using a frequency domain method similar to that of Kim and Noah [10]. Stability issues are examined using standard perturbation techniques along with well known results from Floquet theory [START_REF] Richards | Analysis of Periodically Time-varying Systems[END_REF].

Experimental studies having as their aim the demonstration of subharmonic and chaotic motions in mechanical systems with clearances have concentrated mostly on small laboratory bench-top set-ups which employ relatively simple non-linear elements. Examples include a beam subjected to a rigid support with clearance [START_REF] Masri | [END_REF][18][19] or an impacting pendulum [20,21]. Typically, the apparatus and system parameters are designed such that simple time domain data acquisition and analysis are possible. In this study, a very practical and considerably more complicated engineering application, namely, a precision spur gear pair with backlash and time-varying mesh stiffness characteristics, is used to demonstrate that the proposed governing equation ( 1) is indeed pertinent to the study of rotating machinery dynamics. Experimental measurements are used clearly to validate the proposed analytical model with an emphasis on describing jump discontinuities in steady state forced response curves.

scope

A generalized solution methodology is proposed to obtain the steady state forced response of equation (1) for the special case b=0 such that u(t)=Lt. Analytical solutions are compared with numerical integration results and experimental measurements obtained from a gear dynamics test rig. The main emphasis is to obtain analytical solutions which adequately describe measured forced response phenomena in the fundamental and superharmonic resonance regimes. An approximate harmonic solution is obtained by using describing functions in order to arrive at mathematical conditions for the occurrence of non-linear behavior.

GENERAL SOLUTION BY HBM

The steady state solution x(t)=x(t+T) of equation ( 1) is assumed to be periodic in t with fundamental period T=2p/L. Accordingly, g[x(t)]=g(t)=g(t+T) must also be periodic in t with the same period. Therefore, both functions can be expressed in terms of Fourier series [10]:

x(t)=u 1 + s R r=1 [u 2r cos (ru)+u 2r+1 sin (ru)], g(t)=n 1 + s R r=1 [n 2r cos (ru)+n 2r+1 sin (ru)]. (2a, b)
Equation ( 2) is substituted into the governing equation ( 1) and harmonic balance is enforced by collecting even and odd terms of like frequency to form the vector equation S=0, where the (2R+1) elements of S are given by

S 1 =n 1 + 1 2 s R i=1 [n 2i w 2i +n 2i+1 w 2i+1 ]-f 1 , (3a)
S 2r =-(rL) 2 u 2r +2rLzu 2r+1 +n 1 w 2r +n 2r + 1 2 s R i=1 n 2i [w 2(i-r) +w 2(i+r) +w 2(r-i) ] + 1 2 s R i=1 n 2i+1 [w 2(i-r)+1 +w 2(i+r)+1 -w 2(r-i)+1 ]-f 2r , r $ [1, R], (3b)
S 2r+1 =-(rL) 2 u 2r+1 -2rLzu 2r +n 1 w 2r+1 +n 2r+1 + 1 2 s R i=1 n 2i [-w 2(i-r)+1 +w 2(i+r)+1 +w 2(r-i)+1 ]+ 1 2 s R i=1 n 2i+1 [w 2(i-r) -w 2(i+r) +w 2(r-i) ]-f 2r+1 , r $ [1, R]. (3c)
In order to solve S=0, it is first necessary to express the n i in terms of the u i . This is readily accomplished in the frequency domain by making use of the discrete Fourier transform (DFT). The values of x(t) and g(t) at the discrete time t=nh are

x n =u 1 + s R r=1 [u 2r cos (2prn/N)+u 2r+1 sin (2prn/N)], n $ [0, N-1], (4a)
g n = 8

x n -1, 0, x n +1,

x n q1, =x n =E1, x n Q-1, n $ [0, N-1], (4b)
respectively, where h=2p/(NL) and Ne2R to avoid aliasing errors [START_REF] Oppenheim | Digital Signal Processing[END_REF]. The Fourier coefficients of g(t) are determined by taking the inverse DFT equation ( 4b)

n 1 = 1 N s N-1 n=0 g n , n 2r = 2 N s N-1 n=0
g n cos (2prn/N),

n 2r+1 = 2 N s N-1 n=0 g n sin (2prn/N), r $ [1, R]. (5a-c)
The above n i are then substituted into equation (3) and the solution vector u=

[u 1 u 2 u 3 • • • u 2R u 2R+1 ]
T is determined by employing the Newton-Raphson procedure [START_REF] Thisted | Elements of Statistical Computing[END_REF]:

u (m) =u (m-1) -[J -1 ] (m-1) S (m-1) . (6)
Here, u (m) is the mth iterative solution based on the previous (m-1)th guess and J -1 is the inverse of the Jacobian matrix, the elements of which are defined in Appendix A. The iteration procedure described by equation ( 6) is repeated until the vector norm of S (m) is below a specified tolerance. The initial guess u (0) dictates the number of iterations required for convergence as well as the particular steady state solution found in multiple solution regimes.

stability of steady state solution

The stability of a given periodic solution x(t) is determined by examining the stability of the perturbed solution x(t)+Dx(t), using Floquet theory. The variational equation for the perturbation Dx(t) is given by

Dx¨+2zDx˙+c(t)f(t)Dx(t)=0, (7a)
where f(t) is a discontinuous separation function

f(t)= 6 1, 0, =x(t)=q1, =x(t)=E1, (7b)
that is consistent with the discrete formulation of equation (A2j) given in Appendix A. Equation (9a) can be written in the state-space form z˙(t)=G(t)z(t), where z(t)= [Dx(t) Dx˙(t)] T is the state vector and G(t)=G(t+T) is the periodic transition matrix given by

G(t)= $ 0 -c(t)f(t) 1 -2z % . (8)
The stability of the perturbed solution and consequently the stability of the corresponding system solution x(t) are determined by examining the eigenvalues of the monodromy matrix M=Z(T), which is obtained by solving the homogeneous matrix equation Z (t)=G(t)Z(t) given initial condition Z(0)=I 2 . Here, I 2 is the 2×2 identity matrix. One method to compute M is by direct numerical integration of the above equation by constructing c(t), x(t), and hence f(t), from the previously known Fourier coefficients u i and n i . A more efficient method to compute M has been developed by Hsu and Cheng [START_REF] Hsu | [END_REF]. The state matrix G(t) is approximated as a series of step functions G n at N discrete time intervals t=nh, according to

G n = 1 h g nh (n-1)h G(t) dt, n $ [0, N-1]. (9)
An approximate value of M is obtained from the expression

M3 t N-1 n=0 $ I 2 + s M m=1 0 hG m n m ! 1% , ( 10
)
where M is the number of terms in the approximation for the matrix exponential. This method is easily implemented in the present approach, as M can be readily constructed from the previously available DFT results. The complex eigenvalues l 1 and l 2 of M correspond to the so-called Floquet multipliers which must satisfy the relation l 1 l 2 =exp(-2zT). The solution x(t) is clearly unstable if the modulus of either l 1 or l 2 is greater than unity. Bifurcations may also be classified by observing the behavior of the Floquet multipliers as they leave the unit circle [START_REF] Guckenheimer | Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields[END_REF]. For instance, a saddle node bifurcation, which corresponds to the merger of two periodic solutions, occurs when one multiplier leaves the unit circle at +1 while the other remains inside. A flip bifurcation, which corresponds to the transition from an unstable periodic solution to a stable periodic solution having twice the period, occurs when one multiplier leaves the unit circle at -1 while the other remains inside. Finally, a stable solution will become unstable and doubly periodic whenever a complex conjugate pair leaves the unit circle.

harmonic solution

In an attempt to obtain various closed form expressions which can be used approximately to characterize the system, a single term harmonic balance approach is employed by using describing functions to account for any separations allowed by g[x(t)]. The stiffness and forcing functions are both assumed to be harmonic at frequency kL, which yields the following reduced form of the governing equation

x¨(t)+2zx˙(t)+{1+w a sin [ku(t)]}g[x(t)]=f m +f a sin [ku(t)+G]. (11)
Here, subscripts m and a refer, respectively, to the mean and alternating components of the stiffness and forcing functions with w m =w 1 =1, w a =zw 2 2k +w 2 2k+1 Q1, f m =f 1 , f a =zf 2 2k +f 2 2k+1 and relative phase G=tan -1 (f 2k+1 /f 2k )-tan -1 (w 2k+1 /w 2k ). The steady state solution is assumed to be purely harmonic: where amplitudes u m and u a and phase angle a are the unknowns. Letting 8(t)=kLt+a and substituting equation [START_REF] Pfeiffer | Rattling models from deterministic to stochastic processes[END_REF] into equation (1c) yields

x(t)=u m +u a sin (ku+a), (12)
g[x(8)]= 8 u m -1+u a sin 8, 0, u m +1+u a sin 8, 8 $ [0, 8 - 1 ) or 8 $ (8 + 1 , 2p], 8 $ [8 - 1 , 8 - 2 ] or 8 $ [8 + 2 , 8 + 1 ], 8 $ (8 - 2 , 8 + 2 ), (13)
where 8 2 1 =p+sin -1 g 2 and 8 2 2 =2p-sin -1 g 2 are transition angles that correspond to the zeros of the equation x(8)=21, respectively, with g 2 =(u m 21)/u a . Equation ( 13) can be approximated as a harmonic expression in terms of the so-called describing functions n m and n a :

g[x(8)]3n m +n a sin (ku+a). ( 14a
)
The describing functions n m and n a are given below for the three possible contact regimes shown in Figure 2: no impact (NI) where u m -u a q1; single sided impact (SSI) where =u m -u a =E1; and double sided impact (DSI) where u m -u a Q-1:

n m = 1 2p g 2p 0 g[x(8)] d8= 8 u m -1, 1 2 (u m -1)+A -/p, u m +(A --A + )/p, NI, SSI, DSI, (14b)
n a = 1 p g 2p 0 g[x(8)] sin 8 d8= 8 u a , 1 2 u a +B -/p, u a +(B --B + )/p, NI, SSI, DSI, (14c)
A 2 =(u m 21) sin -1 g 2 +u a z1-g 2 2 , (14d) B 2 =u a sin -1 g 2 +(u m 21)z1-g 2 2 . ( 14e
)
Substituting equations ( 12) and ( 14) into equation (11) and enforcing harmonic balance yields the following set of non-linear algebraic equations:

n a -k 2 L 2 u a +n m w a cos a-f a cos (G-a)=0, (15a) 2zkLu a -n m w a sin a+f a sin (G-a)=0, ( 15b
)
n m -f m + 1 2 n a w a cos a=0. (15c)
The above equations may be solved for u m , u a and a separately for each impact regime by substituting the appropriate describing functions from equations (14b, c). An expression for the transition frequencies separating linear NI and non-linear SSI solution regimes is readily obtained for the special case f a =0 by solving for the values of L and a which satisfy the condition u m -u a =1. The transition frequency at both sides of the kth superharmonic resonance is given by

L 1,2 = 1 k z 1-2z 2 3 zw 2 a -4z 2 (1-z 2 ) . ( 16
)
The above equation is accurate in regions of the kth superharmonic resonance since the steady state solution is very nearly harmonic with frequency kL in these regimes regardless of the harmonic content of c; this is demonstrated later in Section 3.5. For a special case of a lightly damped system in which zW1, equation ( 16) reduces to L 1,2 2z13w a /k. Furthermore, since the discriminant w 2 a -4z 2 (1-z 2 )q0 in order for the transition frequencies to be real, all solutions must be in the NI regime so long as w a Q2z. Hence, the existence of separations is dictated entirely by w a and z for lightly damped systems (zQ0•1) and is completely independent of preload f m .

general hbm solution in the NI regime

A general solution for the linear NI case is obtained by substituting n 1 =1+u 1 and n i =u i , i $ [1, 2R+1] into equation [START_REF] Hsu | [END_REF]. The conditions S=0 is satisfied by solving the linear equation Au=b, where the elements of square matrix A and column vector b, both of dimension (2R+1), are given respectively by

a 11 =2, ( 17a
)
a 1,2i =a 2i,1 =w 2i , i $ [1, R], (17b)
a 1,2i+1 =a 2i+1,1 =w 2i+1 , i $ [1, R], (17c)
a 2i,2i =1-(iL) 2 + 1 2 w 4i , i $ [1, R], (17d)
a 2i,2i+1 =2z(iL)+ 1 2 w 4i+1 , i $ [1, R], (17e)
a 2i+1,2i =-2z(iL)+ 1 2 w 4i+1 , i $ [1, R], (17f)
a 2i+1,2i+1 =1-(iL) 2 -1 2 w 4i , i $ [1, R], (17g)
a 2i,2j =a 2j,2i = 1 2 [w 2(j-i) +w 2(j+i) ], i $ [1, R], j $ [1, (i-1)], (17h)
a 2i,2j+1 =a 2j+1,2i = 1 2 [-w 2(j-i)+1 +w 2(j+i)+1 ], i $ [1, R], j $ [1, (i-1)], (17i) a 2i+1,2j =a 2j,2i+1 = 1 2 [w 2(j-i)+1 +w 2(j+i)+1 ], i $ [1, R], j $ [1, (i-1)], (17j)
a 2i+1,2j+1 =a 2j+1,2i+1 = 1 2 [w 2(j-i) -w 2(j+i) ], i $ [1, R], j $ [1, (i-1)], (17k) b={2f 1 +2 w 2 w 3 w 4 w 5 ... w 2K w 2K+1 0 0 0 ...} T . ( 18
)
Solving the above system of equations directly requires considerably less computational effort than the Newton-Raphson procedure defined by equation [START_REF] Natsiavas | [END_REF]. This leads to an efficient general algorithm to solve the governing equation (1). In NI regimes bounded approximately by the transition frequencies given by equation ( 16), the above set of equations is solved directly for u by using a Gaussian elimination procedure [START_REF] Thisted | Elements of Statistical Computing[END_REF], whereas equation ( 6) is employed to determine u in SSI and DSI regimes.

EXPERIMENTAL VALIDATION

In order to validate equation ( 1) and demonstrate several non-linear phenomena exhibited by the mechanical oscillator of Figure 1, a unique gear dynamics test rig was designed and constructed. Specialized instrumentation and data acquisition and signal processing software were developed to measure in real time the vibratory displacement of a gear pair operating under conditions typical of an industrial rotating machine. Parameters required by equation (1) were determined experimentally and analytical solutions obtained by HBM and numerical integration are compared with experimental steady state measurements.

application example: a gear pair

As a practical application of the proposed governing equation (1), consider a purely torsional two-degree-of-freedom (DOF) model of a gear pair supported by rigid mounts as shown in Figure 3. The rotational angle of each gear u i (t)=V i t+u i a (t) is written as the sum of a nominal rotation angle V i t and a vibratory component u i a (t) which arises due to the gear meshing action. Here t represents real time, and V i is the mean rotational velocity of gear i. The V i and V j are related kinematically by r i V i =r j V j , where r i is the base radius of gear i. Since the system is semi-definite, the dynamic model may be reduced to a single DOF in terms of the so-called dynamic transmission error co-ordinate j(t)=r i u i a (t)+r j u j a (t), which represents relative displacement across the gear mesh interface. Here, subscripts i and j denote gears i and j respectively. The corresponding equation of motion is given by m eq j0+c eq j'+k[u*(t)]gˆ{j(t)-e[u*(t)]}=F, (19a) Here, (') denotes differentiation with respect to t, m eq =I i I j /[I i (r i ) 2 +I i (r i ) 2 ] is the system equivalent mass determined from the polar mass moment of inertia I i of each gear, c eq is an equivalent viscous damping coefficient which describes frictional losses of the system, F is the mean load carried by the gear teeth and 2b is the total backlash. The instantaneous mesh stiffness of the gear pair k[u*(t)] varies periodically according to the number of teeth in contact for a given system rotation angle u*(t)=V*t+bj(t). The term e[u*(t)] describes manufacturing errors associated with mating gear teeth and must also be periodic in u*.

gˆ{j(t)-e[u*(t)]}= 8 j(t)-e[u*(t)]-b, 0, j(t)-e[u*(t)]+b, j(t)-e[u*(t)]qb, =j(t)-e[u*(t)]=Eb, z(t)-e[u*(t)]Q-b. (19b)
Thus, e gives rise to external forcing which may be neglected for precision gear pairs with no intentional tooth flank modifications. The system rotational angle is written in terms of the system velocity V*=V i /n i and modulation depth b=n i /r i , where n i =N i /GCF(N i , N j ) is an integer and GCF(N i , N j ) is the greatest common factor of gear tooth numbers N i and N j . Hence the fundamental frequency V* of k and e is n i N i times lower than the so-called gear meshing or tooth engagement frequency V=n i N i V i . However, if tooth-to-tooth variations of each gear are kept small by precision manufacturing, then both k and e may be assumed to be periodic at the meshing frequency V.

Let t=v n t be a dimensionless time parameter, where v n =zk m /m eq is the undamped natural frequency of the corresponding linear time-variant system and k m is the mean value of k[u(t)]. Accordingly, the equation of motion ( 21) may be reduced to the governing equation (1) for the case e[u*(t)]=0 by defining x(t)=j(t)/b, z=c eq /2zm eq k m and c[u(t)]=k[u(t)]/k m . The dimensionless forcing frequency is given by L=V/v n and f [u(t)]=f 1 =F/k m b in accordance with the mean and alternating convention defined earlier.

The above single-DOF model is valid for precision gear pairs provided that the equivalent translational stiffness of shafting, bearings and supporting structures are sufficiently large relative to the equivalent mesh stiffness [4]. In cases in which the amplitude of j(t) is relatively small (less than 10 mm), angle modulation effects are negligible and the assumption b=0 is valid [1].

experimental gear dynamics test rig

The four-square gear test rig shown in Figure 4 was designed precisely to measure j(t) for a gear pair operating under a range of speed and load conditions. The lubrication bonnet and several guards have been removed for clarity in this photograph. Several unique features are required in order to ensure that the behavior of the test rig is indeed described by equation (1). These features are labelled in the schematic diagram shown in Figure 5. A precision test gear pair (i) is straddle mounted between oversized ultra-precision spherical roller bearings (ii) which are themselves supported by rigid bearing pedestals (iii). The span of the test gear shafts (iv) is kept to an absolute minimum and the shafting diameter is sufficiently large to ensure that any translational or shaft bending modes are well above the primary torsional resonance of the test gear pair so that the SDOF model is valid. The test gear pair is dynamically isolated from the slave gear box (v) by compliant elastomer couplings (vi) and flywheels (vii). Precision helical slave gears (viii) are machined with a different number of teeth from the test gears so that the meshing frequencies are not commensurate. Furthermore, the slave gear tooth flanks have specially ground topological modifications to ensure minimal excitation over a range of operating loads. A power recirculation arrangement is employed to ensure precise torque regulation and efficient operation. System torque is manually adjusted by loosening a split hub coupling (ix) and applying a static torque using a locking pin (x) and a removable torque arm with dead weights that are not shown. Once system torque is applied, the split hub is tightened and the torque arm removed. Power is supplied by a DC motor (xi) equipped with electronic speed control. In order to keep damping to a minimum, the system is lubricated with low viscosity automatic transmission fluid under controlled elevated temperatures. Lubrication is applied to the gear mesh exit and all four bearings at controlled flow rates to ensure consistent film thickness and hence constant damping ratio. Experimental investigations indicated that the system damping is dictated primarily by the bearings under the above lubrication conditions; thus employing a linear viscous damping term in equation (19a) is valid even when tooth separations take place.

The test gears are of a standard 50 tooth, 3•0 mm module, 20°pressure angle unity ratio design such that r=70•477 mm and were ultra-precision ground with no intentional flank modifications such that e=0 and f a =0. The spectral content of k a is dictated primarily by the so-called involute contact ratio which may be varied over a suitable range by selecting mating gears of the same nominal design but having different major diameters. The mean stiffness k m depends mostly on the gear face width and involute contact ratio. Accordingly, Figure 5. A schematic of the test rig: i, test gear pair (AGMA Class 14); ii, ultra-precision spherical roller bearings; iii, rigid bearing pedestals to reduce structural dynamics; iv, test gear shafts that are very stiff in bending and torsion; v, slave gear box necessary for power recirculation; vi, compliant elastomer couplings and vii, flywheels for torsional isolation between test and slave gears; viii, slave gears with topological modifications for reduced excitation; ix, split hub; x, locking pin for torque adjustment; xi, 20 hp DC motor with electronic speed control; xii, location of torsional accelerometers; xiii, slip rings to connect accelerometer signals to conditioning amplifiers; xiv, optical encoders (18 000 lines) to control data sampling rates; xv, drive belt; xvi, polymer granite base. k m may be adjusted for a given contact ratio by altering the gear face width. In this way, parametric studies can be conducted with f m and c as the independent variables with z constant.

instrumentation

The measurement of z(t) to the necessary accuracy of 0•01 mm over a typical range of operating speeds from 200 to 4000 shaft rpm proved to be a very challenging task. This is achieved by attaching four tangentially mounted piezoelectric accelerometers (xii) to the gear wheel as shown in Figure 6a. Slip rings (xiii) are used to connect the accelerometers to the necessary signal conditioning amplifiers. Signals from all four accelerometers z i (t), i=1 to 4, are added vectorally to cancel gravitational effects as well as any transverse vibrations which may occur. Although two diametrically opposed accelerometers are sufficient, four are used in order to achieve higher measurement accuracy by averaging random errors. The resulting angular acceleration of each gear is then added in the time domain to obtain j0(t) and subsequently integrated twice to yield j'(t) and j(t), respectively, by using the analog circuitry shown in the block diagram of Figure 6(b). This method allows direct measurement of j(t) and its time derivatives to an extremely high accuracy of less than 0•001 mm at 1000 Hz. The measurement accuracy is highly dependent upon the excitation frequency, since displacement is directly proportional to the measured acceleration amplitude by the square of the excitation frequency.

data acquisition and signal processing

Analog signals representing j(t) and j'(t) are captured and analyzed in real time by using a high speed multiple channel signal analysis workstation. Synchronous time averaging is employed by controlling the analyzer sampling rate from a clock pulse obtained from a high precision 18 000 line optical encoder (xiv) mounted on the test gear shaft. This not only results in a high signal to noise ratio but, more importantly, allows all analyses to be easily conducted with respect to dimensionless time t once v n is known from an experimental modal analysis. Averaging also diminishes the influence of any slight differences among individual teeth that exist due to unavoidable manufacturing errors; the effects of which were otherwise neglected by the periodicity assumptions in Section 3.1. All Fourier transforms are performed in the so-called orders domain rather than the frequency domain. The measured j(t) and j'(t) are normalized to obtain x(t) and x'(t) respectively. The rth harmonic amplitude A r of x(t) is determined by

A r = X s N i +BW/2 s i =N i -BW/2 X (rs i )X *(rs i ), (22a)
where X (s i ) is the complex valued discrete Fourier coefficient of x(t) corresponding to shaft order index s i , N i is the number of teeth on gear i, BW is the analysis bandwidth in shaft orders and (*) denotes complex conjugate. A typical value of BW=4 was chosen in order to account for sidebands and leakage and bin errors inherent to experimental spectrum analysis. An equivalent root-mean-square (r.m.s.) amplitude A rms of the alternating component of x(t) is given by

A rms = X s R r=1 A 2 r , (22b)
where the maximum number of harmonic terms R available to approximate A rms depends on the operating mesh frequency range 200EVE3400 Hz and the upper frequency limitation of the angular accelerometer VE5000 Hz. Although the alternating component of x(t) can be measured very accurately, the mean value of x(t) cannot be obtained from accelerometer measurements.

solution validation

Experimental values of A rms are compared with equivalent r.m.s. amplitudes obtained from a three-term harmonic balance solution in Figure 7 over the frequency range 0•2ELE1•1, which includes the fundamental and first two superharmonic resonances. Stiffness values w i were determined experimentally for a normalized preload of f 1 =1•0, given v n =19 500 rad/s and b=0•015 mm. The value of z was estimated to be 0•01 from experimental data. Also shown in Figure 7 are corresponding solutions obtained by direct ; and (e) phase portraits j' 1 (t) versus j1(t) and j' 2 (t) versus j2(t). The system parameters are the same as in Figure 7.

time domain numerical integration of equation ( 1) by using a fifth-sixth order variable step size Runge-Kutta algorithm that are in excellent agreement with the HBM solution. Overall, the analytical solution matches quite well with the experimental data which exhibited only NI and SSI responses for the particular physical configuration tested. The response of the physical system was extremely well behaved and repeatable. In regions of multiple solutions, the lower branch was obtained by slowly increasing the system speed until jump-up occurred. Similarly, the upper branch was obtained by slowly decreasing system speed until jump-down occurred. The frequencies at which these transitions occurred were very distinct and could be detected audibly by a 10-40 dB change in sound level. Differences between analytical and experimental results in the upper solution branches are attributed to inaccuracies in the simplified mesh interface model as well as slight manufacturing variances that caused e and f a to be other than zero under actual operating conditions.

The measured time domain response of the system within double solution regimes in the vicinity of the fundamental (k=1) and first superharmonic (k=2) resonance is shown in Figures 8 and9, given L=0•9 and L=0•45, respectively. The measurements are presented in dimensional form in these two figures in order to demonstrate how small the amplitude of j is in practice. In both cases, the measured response is nearly harmonic with frequency kL. The harmonic content of the analytical solution presented in Figure 7 is examined by plotting the forced response curves A rms , A 1 , A 2 and A 3 versus L for NI and SSI solutions only, as shown in Figure 10. Clearly, the steady state forced response in the vicinity of the kth resonance is nearly harmonic so that A rms 3A k . The amplitude A k is dictated almost The system parameters are the same as in Figure 7.

entirely by the kth harmonic stiffness amplitudes w 2k and w 2k+1 and is influenced only slightly by stiffness harmonics other than those of frequency kL. Hence, for practical purposes, the system response in regions near the fundamental and superharmonic resonances may be described by a single-term harmonic balance solution at frequency kL. However, additional harmonic terms must be included in the corresponding NI solution in off-resonance regimes. Also shown in Figure 10 are the corresponding r.m.s. and harmonic amplitudes obtained from an LTV analysis by neglecting backlash completely. Note that the LTV solution becomes unstable in the vicinity of the fundamental (k=1) resonance for the lightly damped case considered here. This occurs because of an unstable parametric resonance attributed to the second harmonic stiffness terms w 4 and w 5 that coincides with the otherwise stable k=1 resonance [1]. Hence, the clearance non-linearity has a stabilizing effect on the forced response which is in exact agreement with the observed behavior of the physical system.

PARAMETRIC STUDIES: HARMONIC SOLUTION

The system response has been shown to be nearly harmonic in resonance regimes where the presence of the clearance type non-linearity dictates the system behavior. Accordingly, the harmonic solution of Section 2.2 is employed to study the effects of w a , z and f m on the steady state forced response of a lightly damped system, zE0•1, with an emphasis on the fundamental and superharmonic resonance regimes.

effect of stiffness variation

The effect of w a on the mean u m and alternating u a components of x(t) is shown in Figure 11 for a given preload f m =0•5 and damping ratio z=0•05. Both stable and unstable SSI and DSI solutions are shown along with the linear NI portion of the steady state forced response over the frequency range 0•5EkLE1•5. For w a Q2z, the system behaves in a linear fashion and no jumps or unstable solutions occur. As w a is increased beyond 2z, two stable solutions coexist. The shape of the forced response curve for w a Q0•25 is similar to that of a forced Duffing's equation of the softening type. Both jump-up and jump-down frequencies and the amount of frequency overlap between the upper and lower stable solutions are highly dependent upon w a . The lower jump-up frequency is nearly equal to the transition frequency L 1 given by equation ( 16) and this corresponds to a fold bifurcation. While the upper transition frequency L 2 predicted by equation ( 16) is evident from a corresponding inflection in the forced response curve, the actual jump-down frequency occurs at a much lower frequency. This behavior is consistent with the experimental forced response curves shown in Figure 12, which were obtained by overlaying several measured response curves A k versus kL, with each harmonic having a distinct value of w a for a given gear involute contact ratio. As w a is increased further such that u m -u a Q-1, a DSI solution is revealed which folds back to the right as shown in Figure 11

effect of damping ratio

The effect of z on the steady state forced response is illustrated in Figure 13, given f m =0•5 and w a =0•25 over the same frequency range. For 2zqw a , the system behaves in a linear fashion. The effect of decreasing z is similar to increasing w a with the exception that the jump-up frequency is only slightly affected by z, which is in agreement with equation ( 16). As z is decreased further the solution bifurcates to a DSI solution.

effect of preload

The effect of f m on the steady state forced response is illustrated in Figure 14, given z=0•05 and w a =0•25 over the same frequency range. Since w a q2z, the system will exhibit jump phenomenon regardless of f m . Such behavior is quite different from that of a clearance-type oscillator with time-invariant stiffness, where increasing the preload was shown to prevent any separations from occurring [4]. The jump-up frequency associated with the lowest solution branch is independent of preload and is governed entirely by w a and z. However, the jump-down frequency is highly dependent on preload becoming lower as f m is decreased. This behavior is consistent with the experimental forced response curves shown in Figure 15, although some variation in the jump-up frequency does occur as f m is varied. This may be attributed to the fact that e and f a are not exactly zero under actual operating conditions. The frequency overlap between the upper and lower SSI solutions tends to zero as f m :a. DSI solutions exist only for higher values of f m and were not observed in the physical test system because of loading limitations.

CONCLUSIONS

A new differential equation is introduced having periodic stiffness variation and clearance type non-linearity which has been shown experimentally to describe accurately non-linear behavior observed in the steady state forced response of a class of rotating machinery in the fundamental and superharmonic resonance regimes. An efficient HBM solution methodology is presented that employs a direct Newton-Raphson procedure to solve the resulting system of equations in the non-linear solution regimes and a Gaussian elimination procedure in linear solution regimes. This procedure may be extended to multiple-DOF systems. Analytical solutions are validated by comparison with numerical integration results and experimental measurements obtained from a gear dynamics test rig. An approximate harmonic solution is obtained which yields mathematical conditions for the occurrence of non-linear behavior in the system response.

Well behaved coexisting stable response solutions have been clearly demonstrated using a relatively common mechanical system that has significant relevance to automotive, aerospace, marine and industrial power transmission and gearing applications. The gear pair system used as the engineering application is inherently non-linear and exhibits a range of complex behavior beyond that reported in this publication. For instance, subharmonic resonances, as well as aperiodic motions and chaotic behavior, have been observed and reported in another recent paper [START_REF] Kahraman | Interactions between commensurate parametric and external excitation in systems with clerance type non-linearity[END_REF].

Figure 1 .

 1 Figure 1. The physical system: a mechanical oscillator of unit mass.

Figure 2 .

 2 Figure 2. Non-linear restoring function g(8) for (a) no-impact (NI), (b) single sided impact (SSI) and (c) double sided impact regimes (DSI).

Figure 3 .

 3 Figure 3. An application example: a single-DOF torsional model of a gear pair.

Figure 4 .

 4 Figure 4. The experimental test rig.

Figure 6 .

 6 Figure 6. (a) The arrangement of linear accelerometers on a gear wheel and (b) the block diagram of the analog circuitry used to compute j(t) from measured accelerometer signals zi (t).

Figure 7 .

 7 Figure 7. A comparison of HBM solutions (--) with numerical integration (q) and experimental (w) results. R=3, K=3, with w2=0•15, w4=0•075, w6=0•06 and all other wi=0, f1=1•0 with all other fi=0, and z=0•01.

Figure 8 .

 8 Figure 8. Measured coexisting steady state solutions j1(t) and j2(t) at L=0•9: (a) j1(t) versus Vt/2p; (b) FFT of j1(t); (c) j2(t) versus Vt/2p; (d) the FFT of j2(t); and (e) phase portraits j' 1 (t) versus j1(t) and j' 2 (t) versus j2(t). The system parameters are the same as in Figure7.

Figure 9 .

 9 Figure 9. Measured coexisting steady state solutions j1(t) and j2(t) at L=0•45: (a) j1(t) versus Vt/2p; (b) the FFT of j1(t); (c) j2(t) versus Vt/2p; (d) the FFT of j2(t); (e) phase portraits j' 1 (t) versus j1(t) and j' 2 (t) versus j2(t).The system parameters are the same as in Figure7.

Figure 10 .

 10 Figure 10. A comparison of the harmonic content of multiple term HBM solutions (W) with the corresponding LTV system results (-----); (a) Arms versus L; (b) A1 versus L; (c) A2 versus L; (d) A3 versus L. The system parameters are the same as in Figure 7.

  for w a q0•25. Transitions from NI to SSI are readily detected by the deviation of u m from its otherwise constant value in NI regimes, as shown in Figure 11(b).

Figure 11 .

 11 Figure 11. The effect of stiffness variation wa on the (a) alternating ua and (b) mean um components of the steady state forced response for wa $ [0•1, 0•5], fm=0•5, and z=0•05. --, Stable solutions; • • • • •, unstable solutions.

Figure 12 .

 12 Figure 12. The effect of stiffness variation wa on the measured alternating component ua of the steady state forced response for wa20•20, 0•15, 0•10, 0•05 and 0•02, z20•01 and fm=0•5.

Figure 13 .

 13 Figure 13. The effect of the damping ratio z on the (a) alternating ua and (b) mean um components of the steady state forced response for z $ [0, 0•1] with an increment of 0•02, fm=0•5 and wa=0•25. --, Stable solutions; • • • • •, unstable solutions.

Figure 14 .

 14 Figure 14. The effect of the preload fm on the (a) alternating ua and (b) mean um components of the steady state forced response for fm $ [0•25, 1•25] with an increment of 0•25, z=0•05 and wa=0•25. --, Stable solutions; • • • • •, unstable solutions.

Figure 15 .

 15 Figure 15. The effect of the preload fm on the measured alternating component ua of the steady state forced response for fm21•0, 0•5 and 0•25, wa20•20 and z20•01.

ACKNOWLEDGMENTS

The authors wish to thank the Powertrain Division of General Motors Corporation for supporting this research effort. Mr S. Moorman and Mr T. Vogel are acknowledged for their assistance in developing signal processing software and performing experimental analyses, respectively.

APPENDIX A: ELEMENTS OF THE JACOBIAN MATRIX

The elements of the Jacobian matrix J ij =1S i /1u j are given by

where

is the Krohnecker delta. Meanwhile, the partial derivatives 1n i /1u j in equations (A1) are determined directly from equations ( 4) and ( 5) as [10] 1n

where f n is a discrete separation function given by


Nicholas J Carpenter

Robert L Taylor

Michael G Katona

LAGRANGE CONSTRAINTS FOR TRANSIENT FINITE ELEMENT SURFACE CONTACT

A new approach to enforce surface contact conditions in transient non-linear finite element problems is developed in this paper. The method is based on the Lagrange multiplier concept and is compatible with explicit time integration operators. Compatibility with explicit operators is established by referencing Lagrange multipliers one time increment ahead of associated surface contact displacement constraints. However, the method is not purely explicit because a coupled system of equations must be solved to obtain the Lagrange multipliers. An important development herein is the formulation of a highly efficient method to solve the Lagrange multiplier equations. The equation solving strategy is a modified Gauss-Seidel method in which non-linear surface contact force conditions are enforced during iteration. The new surface contact method presented has two significant advantages over the widely accepted penalty function method: surface contact conditions are satisfied more precisely, and the method does not adversely affect the numerical stability of explicit integration. Transient finite element analysis results are presented for problems involving impact and sliding with friction. A brief review of the classical Lagrange multiplier method with implicit integration is also included.

INTRODUCTION

Surface contact kinematic conditions can be enforced by prescribing displacement constraints to prevent structural or continuum domains from overlapping and to control surface contact sliding. Surface contact also involves contact force conditions, typically consisting of a tension limit condition for normal forces and a friction limit condition for tangential forces. Because of the non-linearity associated with surface contact force conditions, an iterative strategy is generally required to obtain a precise solution.

Lagrange multiplier methods and penalty function methods are the two most common used approaches to enforce finite element surface contact displacement constraints. Lagrange multiplier methods are alternatively referred to as mixed or hybrid variational methods by some authors, and penalty methods are commonly referred to as 'contact', 'gap', or 'joint' element methods. For transient analyses by explicit integration, penalty methods have received the most attention in the literature and in commercial finite element programs. Some of the recent work involving the use of Lagrange methods and penalty methods for finite element surface contact is found in References 1-3, 6, 8, 10, 13 and 14. And an extensive reference to the literature on this subject is given in Reference 14.

The primary focus herein is on Lagrange multiplier methods. The presentation begins in Section 2 with a preliminary discussion of the finite element equation of motion and two dimensional surface •contact. A brief review of the classical Lagrange multiplier method is presented in Section 3, where it is shown that the classical method is not compatible with explicit integration operators. In Section 4 an alternative formulation that is compatible with explicit operators is presented, which is referred to as the 'forward increment Lagrange multiplier' method. Compatibility with explicit operators is established by referencing the Lagrange multipliers one time increment ahead of the associated surface contact displacement constraints. A one dimensional impact example involving a single contact constraint is presented in Section 5. Finite element analysis results are compared for alternative methods of enforcing the contact constraint. Enforcing the constraint by the forward increment Lagrange multiplier method leads to a highly accurate and well behaved solution. By contrast, the performance of the classical Lagrange method with an implicit integration operator is shown to be poor. Results obtained using a penalty function method to enforce the constraint are also presented.

A two dimensional finite element surface contact formulation based •on the forward increment Lagrange multiplier method is developed in Sections 6, 7 and 8. Kinematic conditions and displacement constraints are considered in Section 6. In Section 7 the Gauss-Seidel method is introduced to solve the coupled forward increment Lagrange multiplier equations. The Gauss-Seidel method is then modified in Section 8 to allow for the enforcement of contact force conditions during iteration. In Section 9 the forward increment Lagrange multiplier method and the modified Gauss-Seidel method are exercised to solve a two dimensional surface impact example. A two dimensional finite element sliding problem is then presented in Section 10, followed by closing remarks in Section 11.

The forward increment Lagrange multiplier method is an extension of the ideas presented in Reference 14. The most significant contribution in the present paper is the formulation of an efficient method to solve the coupled forward increment Lagrange multiplier equations that arise in two dimensional surface contact.

CONSTRAINED EQUATION OF MOTION

The finite element semi-discretized equation of motion is expressed in general form as Mii + F(u, it) = R [START_REF] Bathe | A solution method for planar and axisymmetric contact problems[END_REF] in which M is the mass matrix, u is the vector of displacement degrees of freedom, li is velocity, ii is acceleration, F is the internal force vector and R is the external force vector. In addition to the usual prescribed boundary conditions, it is assumed that the solution of equation ( 1) is also subject to surface contact displacement constraints.

An illustration of two dimensional surface contact between bodies that are spatially discretized using low order continuum finite elements is shown in Figure 1. Contactor nodes are denoted by a 'C' and target nodes by a 'T'. Displacement constraints are prescribed to prevent the contactor nodes from penetrating the target domain and to control tangential sliding of contactor nodes along target surfaces. These constraints may be expressed as G{u+X}=O [START_REF] Belytschko | A three-dimensional impact-penetration algorithm with erosion[END_REF] where X is the material co-ordinate vector, the sum of u and X is the spatial co-ordinate vector, and G is a surface contact displacement constraint matrix. The components of G are typically unknown a priori and generally change as displacement and deformation occur. Starting from a configuration in which the surfaces are separated, the motion of contactor and target nodes must be tracked so that displacement constraint components can be introduced in G as contact occurs. During contact the components of G may change with time as required to ensure that the associated contact force reactions satisfy contact force conditions. For example, if a force component normal to a target surface approaches a tension limit force condition, then the associated displacement constraint must be eliminated to allow surface contact separation. Similarly, if a tangential force component approaches a friction limit force condition, then the associated tangential displacement constraint must be relaxed to allow sliding. The components of G change with time when sliding occurs.

It is convenient to first consider the less complicated problem of treating equations (2) as known linear equality constraints that do not change during an integration time increment. This assumption simplifies the initial discussion of Lagrange multiplier methods as presented in the following two sections. A method for relaxing displacement constraints to enforce two dimensional surface contact force conditions is formulated in Section 8.

LAGRANGE MULTIPLIER METHOD

Lagrange multipliers may be introduced into the equation of motion to give

Mii + F(u, u) + GTA = R (3)
where the components of the Lagrange multiplier vector A are the surface contact forces. The Lagrange multiplier method proceeds by treating A as unknown and solving equations ( 2) and (3) simultaneously.

For an elementary small displacement problem in which internal forces are strain-rate independent and proportional to displacement, the constrained equation of motion referenced to time

tn + 1 is Gn+l {un+l +X}= 0 (4a) (4b)
Equations ( 4) may be solved by direct time integration, see for example Hughes et a/. 7 and Bathe and Chaudhary. 1 Herein, the following second order direct time integration operator is considered,

Un + 1 = qo + boAiin + 1 Dn +l = ql + btAiin+l On+ 1 = q2 + b2Aiin+ 1 (Sa) (5b) (5c) where h. l. h 2 "" qo = Un + Un + 2 Un (5d) ql = On + hiin (5e) q2 = Dn (5f) and bo = !h 2 Po (5g) bt = hPt (5h) b2 = 1 (5i) h = tn+ 1 -tn (5j)
This temporal discretization is equivalent to the well-known Newmark method. 12 However, the operator form presented here is formally known as the Beta-2 method, which is a subset of the generalized Beta-m method developed by Katona and Zienkiewicz. 9 Two well-known Beta-2 methods are:

(i) the constant-average-acceleration method, also known as the trapezoidal rule; which corresponds with Po = Pt = !;

(ii) a single step version of the central difference method, which corresponds to Po= 0 and

P -.1 1 -2.
Moreover, a Beta-2 operator is referred to as implicit if Po =I= 0, and explicit if Po = 0. Substituting equations (5) into equations (4) leads to the following incremental equation of motion:

[ [b2M + boK] G~+t] {~iin+1} = {Rn+1-{Mq2 + Kqo}} boGn+ 1 0 An+ 1 -Gn+ 1 {qo + X} (6)
For the surface contact constraints, the rows of Gn+ 1 are linearly independent. Therefore the above system of equations is non-singular if

[b 2 M + b 0 K] is non-singular and b 0 =I= 0.
Conversely, the system is singular for b 0 = 0, thereby excluding admissibility of explicit integration operators.

If the non-linearity associated with surface contact force conditions is considered, then an iterative form similar to equations (6) arises. A general two dimensional non-linear surface contact formulation of this type has been developed by Bathe and Chaudhary, 1 wherein a system of equations similar to equations ( 6) is re-solved for each iteration. The behaviour and accuracy of the method are good for static and slow transient problems. However, the method is ill behaved when inertial forces are relatively large. This is demonstrated by a one dimensional finite element example in Section 5.

FORWARD INCREMENT LAGRANGE MULTIPLIER METHOD

It is easy to show that the contact forces An+ 1 directly influence the forward differences ~iin + 1 , but act one time step too late to influence the history terms, Un , On or iin . Furthermore, the forward difference terms ~iin+ 1 are absent in the second order incremental equation

h • .l h2- b A- Un+l = Un + Un + 2 Un + oL1Un+1 (7)
for explicit integration because b 0 = 0. Therefore, the contact forces A.,.+ 1 have no influence on displacements u,.+ 1 , which explains the explicit operator singularity of equations ( 6).

An alternative formulation that is compatible with explicit integration is developed by relating displacement constraints at time t,.+ 1 with Lagrange multipliers at time t,.. This approach is referred to here as the forward increment Lagrange multiplier method, for which the incremental equation of motion is Mu,. + F(u,., ii,.) + G~+ 1 A.,.= R,.

G,.+l {un+l +X}= 0

(Sa) (8b)
For explicit integration of equations ( 8), consider the well-known multi-step central difference method, which is given by

u:+ 1 = h 2 M-1 {R,.-F(u,., ii,.)} + 2u,.-u,.-1 A.,.= [h 2 Gn+tM-1 G~+1]-1 G,.+ 1 {u:+t +X} c h2M-1GT '\ Un + 1 = - II+ tAn (9a) (9b) (lOa) (lOb) (lOc)
(10d) (10e)
The incremental solution of equations [START_REF] Landers | An augmented Lagrangian formulation for the finite element solution of contact problems[END_REF] proceeds by first calculating. the usual explicit displacement update for all degrees of freedom, including surface contact nodes. These displacements u:+ 1 are denoted with a right superscript asterisk to indicate that the effect of surface contact forces A.,. is neglected. Following the explicit update, the spatial co-ordinates (x: + 1 = u: + 1 + X) are processed by a search algorithm to identify all surface contactor nodes that have penetrated target surfaces. Based on the penetration search information a constraint matrix G,.+ 1 can be assembled_and equations (10d) can be solved for the contact forces A.,.. The incremental displacements u~ + 1 associated with these contact forces are then calculated from equation (lOe), and the total displacements u,.+ 1 are calculated from equation (lOa).

An efficient Gauss-Seidel iterative approach to combining and solving equations (lOd) and (lOe) for two dimensional surface contact is developed in Section 7. A modification is then introduced in Section 8 which relaxes the displacement constraints as required to ensure that A.,. satisfies surface contact force conditions.

An incremental algorithm based on equations [START_REF] Landers | An augmented Lagrangian formulation for the finite element solution of contact problems[END_REF] is given in the Appendix as Algorithm 1. The incremental algorithm for Beta-2 explicit methods is similar to the multi-step central difference algorithm and is discussed briefly in Section 10. Two advantages of the Beta-2 explicit approach are that numerical damping can be prescribed by simply selecting an appropriate value of p 1 greater than t, and only a slight modification is needed for the more general case of non-uniform time increment integration.

ONE DIMENSIONAL IMPACT EXAMPLE

A simple one dimensional surface contact problem is solved by the classical Lagrange multiplier method of equations ( 6), and by the forward increment Lagrange multiplier method of equations [START_REF] Landers | An augmented Lagrangian formulation for the finite element solution of contact problems[END_REF]. A penalty method solution is also presented for comparison.

The problem consists of two linear elastic prismatic rods moving with equal speed in opposite directions. The configuration and properties are given in Figure 2. The rods are initially undeformed and the problem is symmetric about the point at which the two rod faces impact.

The exact solution for displacement, velocity and stress at the contact surface of the left rod is given by

t ~ timpact timpact < t < trelease t ~ trelease where and Ut = g Ut = 0 A= -voJEP U1 = 2vo(gfvo + L JPfE-t/2) timpact = g / Vo (lla) (llb) (llc) (lld) trelease = gfvo + 2LJPjE (lle)
which is plotted in Figures 3(a) to 3(c). The three finite element solutions also plotted in these figures are described in the following.

For the finite element analyses each rod is discretized by 20 equal length uniform strain elements and corresponding diagonal mass matrix. All three finite element solutions are obtained by direct time integration using a uniform time step increment, h = 0•2226E-5 sec.

The 'Lagrange-Implicit' finite element results presented in Figures 3 correspond with the classical Lagrange multiplier method given by equations ( 6). The integration is performed using ,, ' ' ' ' ' ' ' ' ' \ .::::::- § . the implicit operator recommended by Chaudhary and Bathe, 5 for which {3 0 = 1 and P 1 = ! . The solution commences without a contact constraint and continues until the following condition is satisfied:

g -100 w > .... V.CRANGE-INPUCIT -• v.GRANGE -EXPUC!f -200 --PENAllY -EXPUC!f -El1ICT -300 .0 2.5 5.0 ' ' ', ' ' , ' ' , ' ' ' ,
Ut + X 1 > Uz + X 2 n+ l n + l
which indicates penetration. When this occurs the constraint

Ut + X 1 = Uz + X 2 n+l n + l (12a) (12b)
is introduced and On+ 1 is re-solved while concurrently solving for the contact force An+ 1 • The constraint remains active until the contact force becomes tensile, at which time the constraint is inactivated and On + 1 is re-solved to allow for contact surface separation. Fortunately these two conditions for activating and inactivating the displacement constraint are never mutually exclusive. It is observed in Figures 3 that the Lagrange-Implicit results are generally poor with regard to velocity and contact force. And the results worsen if a smaller time step increment is used. Impact and release conditions can be used to improve the Lagrange-Implicit method, as proposed by Hughes et al. 7 The more recent theoretical developments by Moreau 11 might also serve as a basis to formulate impact and release conditions; however, such efforts are not easily extended to the more general problem of two and three dimensional finite element surface contact. As a final comment regarding the Lagrange-Implicit method it is noted that use of the trapezoidal rule method leads to much greater inaccuracy than that observed with the implicit operator employed here.

The 'Penalty-Explicit' finite element results are obtained using a linear penalty function constraint that is active whenever inequality (12a) is satisfied. The contact force is calculated by multiplying the depth of penetration by a penalty parameter. The parameter used is 7•5E + 6lb/in. The transient solution is obtained using the central difference method given by equations ( 9). The ratio of the time integration increment to the critical time step size is 90 per cent for the unconstrained mesh. When the penalty constraint is active the numerical stability ratio increases to 99 per cent. A larger penalty parameter would reduce the contact surface penetration observed in Figure 3(a) but would also further stiffen the system and therefore require use of a smaller time step increment to maintain numerical stability. Overall the Penalty-Explicit results are a good approximation to the exact surface contact solution.

The 'Lagrange-Explicit' finite element results are obtained by the forward increment Lagrange multiplier method with central difference integration as given in equations [START_REF] Landers | An augmented Lagrangian formulation for the finite element solution of contact problems[END_REF]. The results obtained with this method are the most accurate in predicting the contact force and virtually exact for displacement and velocity at the contact surface during contact. It is noted that following surface contact release (i.e. t > tretease) the Lagrange-Explicit solution exhibits a noticeable error in centre of mass velocity, as observed in Figure 3(d). A closer study of the forward increment Lagrange multiplier method reveals that, during the time step in which impact occurs, the collision of surface nodal masses corresponds with the mechanics of a plastic collision and thus an associated loss of kinetic energy. However, in the mesh refinement limit this kinetic energy loss vanishes because surface nodal mass vanishes. Overall the Lagrange-Explicit results are an excellent approximation to the exact surface contact solution, and the accuracy of the method is found to improve with mesh refinement. Furthermore, the forward increment Lagrange multiplier method has no adverse affect on the numerical stability of explicit integration; the stability ratio is actually less than 90 per cent during contact.

The forward increment Lagrange multiplier method of equations [START_REF] Landers | An augmented Lagrangian formulation for the finite element solution of contact problems[END_REF] is outlined in the following for the one dimensional impact problem. Surface contact penetration is in which

Gn+ 1 = [ -1, 1]
If p:+ 1 is negative, then penetration exists and the contact force is calculated in which

M= [~' :,]
The displacements due to this contact force are then 1 m2

Given the symmetries, m 2 = m 1 , X 2 = X 1 + 2g, and ut = -ut , leads to

n+ l n+l P:+ 1 = 2(-ut +g) n+l and if p: + 1 is negative, then * A = mrPn+l n 2h2 u~ -1. p* n+l - 2 n+ 1 (13a) (13b) (13c) (13d) (13e) (14a) (14b) (14c) (14d)
Note that the displacement in equation (14d) corresponds with the exact displacement solution of equation (11 b).

TWO DIMENSIONAL SURFACE CONTACT KINEMATIC CONDITIONS

The explicit displacement update given by equation (10c) allows contactor nodes to penetrate target surfaces. This is illustrated in Figure 4

G, j, =0 (16a)
• +1 • + 1 where -[<1 -ex, ) 0 a I 0 -1 -~] G •+I • +1 (16b) I = • +1 0 (1 -a1 ) 0 a, 0 n+l n+l and - { . }T Xr = XM At Xsr Br Xcr CJ •+1 •+1' Y. + t' •+t' Y.+t' •+t' Y.+t (16c)
The relative position of node CI along the target surface is associated with the value of a, . In

• + 1
the following a, is prescribed to maintain the relative position of node C I, thereby imposing

• +! a tangential constraint to control sliding. If the node C 1 is in contact with the target surface at time t,., then a1 is simply

•+ 1 where v A.BI = XBI -XAI • • V A.Cl = Xcr -XAI • •
The above geometric relations are illustrated in Figure 5. The distance between the node C 1 and the target surface at time tn , measured in the direction normal to the target surface, is

A _ VAClx VA~ly-VACly VABix ilNJ - IVABII ( 18
)
If L\Nr is greater than zero, then node C 1 is not in contact with the target surface. at time tn . In this case the relative position of node C 1 and the associated value for ar should correspond with the

• + 1
point at which contact occurs during the time increment. However, this is quite difficult to determine precisely. Therefore an approximation for the point of contact is based on the assumption that nodal displacements are linear with respect to time during the increment tn to t: + 1 . This leads to

where

L\~1 (19b) a= 6.~1 -ANI a*= 1 -a (19c)
and

V* •V* A* - ACI ABI (19d) n-IV~BII A* -V~clx v~Bly -v~Cly v~Blx Nl- IV~BII (19e) V~BI = x.fr •+1 -xfr • + 1 (19f) V~cr = X~r •+1 -xf1 • + 1 (19g)
The surface contact displacement constraint form given by equations ( 16) is typical and can be found for example in References 1, 4 and 8. The usual implementations are with low order finite elements that employ linear spatial and displacement interpolation functions along target and contactor bodies. With regard to the evaluation of a1 , it might be argued that equations (19) n+l provide a sufficiently accurate approximation for the point of contact even when the contactor node is not in contact with the target surface at time tn. Nevertheless, equations (19) provide a more accurate approximation of the point of contact without significant computational cost or complexity. At the other extreme, attempting to calculate the exact point of contact occurring between time tn and tn+ 1 is generally quite complicated, and may provide no real improvement in accuracy considering the limitations of the finite element and time integration discretizations.

GAUSS-SEIDEL ITERATION

An equation solving strategy based on the Gauss-Seidel method is developed to solve equations (lOd) and (lOe). In the following it is assumed that the 2 x 6 constraint element matrices, given by equation (16b), are assembled to form the global matrix Gn+ 1 • The two rows of Gn+ 1 corresponding with Gr are given by Gr . Equation (10d) is then expressed as

n+ 1 ~t + 1 for I = 1, 2, ... , N N " h 2 Gr M-1 GJ A1 = Gr {u:+l +X} L. ,.-+1 n+ l " n + l (20) J = 1
where M is the global mass matrix, N is the number of contactor nodes, and l1 is the 2 x 1 vector

• consisting of surface contact force components A. .. 1 andA.,1 associated with the J th contactor node.

• •

The i to i + 1 Gauss-Seidel iteration of equation ( 20) is,

for I = 1, 2, ... , N N i,i+lPI = G!+t {u:+l +X}-L h2G!+.M-1G!+tiJ..! J=l I -1 -" h 2 Gr M-1 GT i+ 1 J..1 L.. n+l o+ l n J=l AI.. = [h2G, M-lG]' ]-1 i, i + lp I n + 1 n+ 1 J (21a) (21b) (21c)
Equations ( 2lb) and (21c) may be efficiently evaluated by element level calculations. However, equation (21a) is complicated because coupling generally occurs between the Ith and Jth displacement constraints. Fortunately a more efficient expression for calculating i, i+ 1 PI can be obtained by combining equations (lOd) and (lOe). From equation (lOe), the value of u;+ 1 at the beginning of iteration i to i + 1 is and at the end of the iteration

i c - Un+l- J = l i+1 0 c _ n+l- J = 1
Defining the difference between i + 1 u~+ 1 and iu~+ 1 as

Auc _ i + 1 0 c _ i 0 c L.l - n+1 n+l
and substituting from equations (21c) and ( 22) leads to

N Auc = L Au) J=l (22a) (22b) (23) (24a)
where

~Uc = -h 2 M-1 GJ ~A. J n+i J (24b)
During iteration ito i + 1, just prior to solving for ~11 , the current estimate of u~+ 1 , denoted as

i. i+ 1 0 c is n+ 1' I-1 i, i+ t 0 c _ i 0 c + " A 0 c n+ 1 - n+ 1 L. Ll J (25)
J=l Substituting from equations (22a) and (24) into equation ( 25), and also substituting the difference, e + 1 A.1 -iJ..J ), in place of ~A.J, leads to Based on equations ( 27) and (28), the i to i + 1 Gauss-Seidel iteration can be expressed as,

for I = 1, 2, ... , N i,i+tp =G 1 {i,i+t 0 c +u* +X} I n+i n+l n+l i + 1 AI = iAI + ~~..1 n n ~uc = -h 2 M-1 GT ~A. I n+ 1 I i,i+t 0 c ~ {i,i+t 0 c + Auc} n+l n+l Ll I (29a) (29b) (29c) (29d) (29e)
where the arrow in equation (29e) means: replace i,i+tu~+t by e,i+tu~+t + ~un. In equations (29) all references are to the Jth constraint element only. Therefore the Gauss-Seidel iteration can proceed by simple element level calculations, without assembling any global matrices. The coupling between the Jth and Jth constraint elements that appears in equations (21a) is accounted for very efficiently by the global vector update in equation (29e). Also, the matrix in equation ( 29b) is a simple 2 x 2 diagonal matrix given by

-1 T -2 2 [1 OJ Gr M Gr = ((1 -llr ) f mAI + llr fmsi + 1/mci) 1 n + l n+l n + l n+i O ( 30
)
where mAI, m 8 I and mc 1 are the nodal masses of nodes A 1 , B 1 and CI for the Jth constraint element. The Gauss-Seidel algorithm for a two dimensional surface contact based on equations (29) is given in the Appendix as Algorithm 2. This algorithm corresponds with step 2 of Algorithm 1.

In addition to the steps outlined above, an iterative convergence check is necessary. The following is recommended and appears in step 2.3 of Algorithm 2 (31

)
where IAAI is the mean-square norm of contact force changes, and li+ 1 Anl is the mean-square norm of the i + 1 estimate of An.

TWO DIMENSIONAL SURFACE CONTACT FORCE CONDITIONS

The procedure developed in Section 7 for solving equations (10d) and (10e) leads to surface contact force reactions An that satisfy the prescribed displacement constraints. However, satisfying the displacement constraints may result in contact force reactions that fail to satisfy prescribed surface contact force conditions. When this occurs the associated displacement constraints are relaxed to obtain a modified solution that satisfies the contact force conditions.

In the global co-ordinate system the components of At are A. d and Ayr. These components may be n n n

transformed to obtain the components Art and ANI, which correspond with the directions tangent

TWO DIMENSIONAL IMPACT EXAMPLE

A two dimensional surface contact example is presented to demonstrate the application of Algorithms 2 and 3. The example consists of a contactor surface impacting a stationary target surface during an increment of motion from time tn to tn + 1 • The target and contactor surfaces are illustrated in Figure 6. The target surface includes three segments spanning from nodes 1 to 4, which are motionless at tn. Nodes 5 to 16 are contactor surface nodes, which are in motion at tn . The contactor nodes penetrate the target segments during the increment from tn tot:+ 1 , as shown in Figure 6(a).

The spatial co-ordinate positions of target and contactor nodes at t,. and t:+ 1 are given in Table I. The information in Table I Two different friction coefficient cases are considered. In the first case the friction coefficient Jl is large enough to prevent all contactor nodes from sliding on the target segments during the time III andIV. The time step increment for the calculations is h = 1•0.

Node-A 1 1 1 1 1 2 2 2 2 3 3 3 3 Node-B 1 12
An overview of the convergence behaviour of Algorithms 2 and 3 is given in Figure 7. Mean-square norms are used to measure relative contact force error. The actual error is

I i + 1 l.n -ln J actual error = lln I (33)
where i+ 1 ln is the contact force estimate at completion of the current iteration, and ln is the converged solution. The estimated error is where ;ln is the estimate from the previous iteration. Equation ( 34) is similar to the expression used for the convergence check in step 2.3 of Algorithm 2.

It is observed in Figure 7 that the order of convergence is approximately linear. Nevertheless, an accurate solution is obtained quickly. For example, the actual contact force error is reduced from 100 to 1 per cent in only 3 iterations. Furthermore, the relaxation of tangential displacement constraints in step 2.2.3c of Algorithm 3, which occurs in the frictionless sliding case, has only a slight effect on the convergence behaviour.

Note that contactor nodes 6 and 12 are not in contact with the target at time tn + 1 . Enforcing the displacement constraints initially assumed for these contactor nodes results in tensile normal force reactions. This feature of the example serves to exercise the relaxation of normal displacement constraints in step 2.2.3c of Algorithm 3.

In this example the incremental displacements are large and the pattern of contactor node motions is unusually chaotic. For smaller incremental displacements and a more orderly pattern of motion, Algorithms 2 and 3 convergence more quickly. Also, increasing the number of target segments and contactor nodes has been found to have little effect on convergence behaviour.

A significant factor affecting convergence is the relative magnitude of contactor and target nodal masses. If target node masses are significantly greater than contactor node masses, the Gauss--Seidel method will converge more rapidly. Conversely, if the target nodes are relatively less massive, the convergence rate is slower. This occurs because equation (lOd) becomes more diagonally dominant as the magnitude of target masses increases. This effect should be considered when assigning contactor and target surface designations.

TWO DIMENSIONAL EXAMPLE

As a further demonstration of the forward increment Lagrange multiplier method formulated in Algorithms 2 and 3, two finite element example problems are presented. The finite element Table III. Forward increment Lagrange multiplier solution for surface contact forces 1.. 11 • No sliding for J-t > 0•88, no resistance to sliding for problems are compared with an elementary single degree of freedom (SDOF) sliding contact problem, which is illustrated in Figure 8(a). The SDOF problem consists of a rigid object that slides horizontally on a rigid base, while connected to a fixed point by a linear elastic spring. The mass of the object is specified as 30, the spring stiffness is 10 000 and the undeformed spring length is 10. At time t 0 the object has an initial horizontal velocity of 75, and is located directly above the fixed point such that the spring is undeformed. The first of the two finite element problems is illustrated in Figure 8(b), which features a two dimensional linear elastic object. The object is spatially discretized by 12 bilinear quadrilateral plane stress elements, which are initially square at time t 0 . The modulus of elasticity of the object is 10 6 , Poisson's ratio is 0•25 and the mass density is 2. The spring is modelled by a two node uniform strain .truss element having a modulus of elasticity of lOS, a cross section area of 1 and density of 1•2. Large rotation of the truss is accounted for by incrementally updating spatial co-ordinates during the finite element internal force calculations.

J-t = 0 Contact }-t > 0•88 J-t=O element Contact I node ATI AN I ATI AN I • n n • 1 5 -0•5773 -0•9155 0•0000 -0•9155 2 6 0•0000 0•0000 0•0000 0•0000 3 7 -0•2486 -1•9289 0•0000 -1•9298 4 8 -0•7111 -0•9423 0•0000 -0•9441 5 9 0•6817 -0•9349 0•0000 -0•9316 6 10 -1•0628 -1•8470 0•0000 -1•8454 7 11 -0•4890 -1•8032 0•0000 -1•8025 8 12 0•0000 0•0000 0•0000 0•0000 9 13 0•7551 -1•7067 0•0000 -1•7060 10 14 -0•2748 -1•5813 0•0000 -1•5809 11 15 -0•7292 -1•4872 0•0000 -1•4871 12 16 -1•2441 -1•4244 0•0000 -1•4246
The second finite element problem is illustrated in Figure 8(c), where an elastic base is introduced and discretized by 15 bilinear quadrilateral plane stress elements. The base elements are initially square at time t 0 • The material properties and mass density ofthe base are the same as given above for the elastic object.

For the SDOF system, the equation of motion at time tn is miin -F. + Ax = 0

• • (35a)
where m is the mass of the sliding object, u is the horizontal displacement of the object and Ax is the contact friction force. The horizontal and vertical components of the spring force are where and F; =F.(~:) F: =F.(:.)

(35b) (35c) (35d) (35e)
in which the undeformed spring length isLand the spring stiffness is k. Since the SDOF object is constrained from moving vertically, the vertical contact force is given by -Fy. Equations ( 35) constitute a non-linear initial value problem that is readily solved by direct time integration using the central difference operators given in equations ( 9). An incremental relation for the friction force Ax is obtained by substituting equation (9b) into (35a) and rearranging so that m A:= -h 2 (un +1-2un + Un-d-F:

(36)
This force is defined as a 'stick' force if it prevents sliding (i.e. un+ 1 = un) and is given by

m Ast~t = h 2 (un-Un+ 1) -F~ (37a)
Observing that the vertical force given by equation (35c) is non-negative and introducing a coefficient of friction, p, leads to the following expression for the friction force:

Ax = { A.~~k, I A..t~kl ~ pF~ } • pF ~ sign ( A•t~k) , I Ast~•k I > pF ~ (37b)
The central difference integration solution to equations (35) and (37) is designated as 'SDOF' in Figures 9 and10 solution is calculated using a uniform time step increment of h = 0•000850, which is two orders of magnitude less than the critical time step size.

For comparison to the SDOF solutions, finite element results are also presented in Figures 9 and10. The problem description of Figure 8(b) is referred to as 'FE object', while that of Figure 8(c) is referred to as 'FE object/base'. The finite element solutions are obtained using the Beta-2 explicit integration method with Algorithms 2 and 3. ------------------------------------------------------------------------------------ -------------~----------------~--------------~----------------~--------------~ -------------~--------~~~--------------~------~~~--r-------------~ ------------T--------------r-------------~------------~------------~ ------------T--------------r-------------~------------~------------~ The forward increment Lagrange multiplier method with Beta-2 explicit integration is similar to the multi-step central difference integration in Algorithm 1. The modification of Algorithm 2 for Beta-2 explicit integration is accomplished by simply replacing h 2 , by h 2 ({3 1 + 1/ 2), in steps 2.2.3 and 2.2.5. The Beta-2 method is preferred here because it allows for high frequency numerical damping, which improves the finite element results. The Beta-2 integration parameters used are: {3 0 = 0 and {3 1 = 1. The FE object and FE object/ base problems are solved by the Beta-2 method using a uniform time step increment of h = 0•000850, which is 85 per cent of critical for both finite element problems.

~ N .. ' ' ' SDOF--- FEOBJECT --------- o FE OBJECT/BASE • ••••••••••• 0~-
' ~ CXl ' 0 N SDOF--- FEOBJECT --------- FE OBJECT/BASE •••••••••• ...... .

7~-

~~-

The contact surface of the finite element object is designated as the contactor surface, which consists of seven contactor nodes. The vertical contact force resultant given in Figures 9(b The finite element object sliding on the rigid base is a special case in which target surface mass is assumed to be infinite. Algorithm 2 is easily modified to accommodate this assumption by setting the terms involving the reciprocal of target nodal masses to zero in step 2.2.3. In this special case the matrix in equations (lOd) is diagonal. Consequently, the exact solution for An is obtained by one iteration, at each time step.

By comparison, the problem in Figure 8(c) with the finite element base is solved by the more general methodology of Algorithms 2 and 3. The convergence tolerance in step 2.3 of Algorithm 2 is specified as tol = 0•001. Satisfying this rather strict tolerance requires not more than four iterations at any time step in the transient analyses, and the number of iterations required is nearly independent of the friction coefficient.

The horizontal centre of mass displacement results for the FE object and FE object/base frictionless contact are quite similar, hence the dotted line in Figure 9(a) is visually obscured by the dashed line. It is observed in Figures 9(a) and lO(a) that the horizontal centre of mass motion of the elastic finite element object is nearly identical to the motion of the SDOF object in all cases. This is true even when friction, J1. = 0•5, is prescribed. The high frequency oscillations in the contact force resultants are expected as a consequence of the elastic compliance of the finite element object and the rapid changes in friction forces acting on the object. This is observed in Figure lO(c), where the horizontal contact force resultant is large and suddenly changes direction.

The sliding contact problem presented here is a demanding test of the forward increment Lagrange multiplier method and the modified Gauss-Seidel solution strategy in Algorithms 2 and 3. The forward increment Lagrange multiplier method proves to be highly accurate for enforcing the coupled normal and tangential contact constraints and friction limit force condition. Furthermore, the modified Gauss-Seidel strategy converges rapidly.

CONCLUSION

It should be emphasized that the coupled system of equations involved in the forward increment Lagrange multiplier method is typically small. For two dimensional surface contact, the number of unknowns is simply equal to two times the number of active contactor nodes. Moreover, the incremental displacements of interior nodes and surface nodes away from contact regions are calculated by the usual explicit operator update without loss of efficiency. The computational cost per increment for these degrees of freedom is not affected by the forward increment Lagrange multiplier contact forces.

The modified Gauss-Seidel equation solving strategy formulated in this paper requires little computational effort per iteration, and converges quickly for the forward increment Lagrange multiplier two dimensional surface contact problems. Also, it is quite convenient that non-linear surface contact force conditions can be enforced concurrently in the modified Gauss-Seidel formulation. This combined methodology generally proves to be a very accurate and efficient alternative to the penalty function method. Since the forward increment Lagrange multiplier method does not adversely affect numerical stability, a slightly greater computational cost per increment with the method may be more than compensated in many cases by the advantage of a larger allowable time step size.

For the sake of simplicity, no distinction between static and dynamic friction coefficients is made in this paper. A full range of friction law and contact condition refinements remains a topic for further investigation.

Figure l .

 l Figure l. Two dimensional finite element surface contact illustration

. 1 u

 1 ,. = 2h {un+1-u,._t} u,. = : 2 { u,. + 1 -2u,. + u,. _ t} Direct substitution from equations (9) into equations (8) leads to where 1 ii,. = h { u,. -u,.-d

Figure 2 .

 2 Figure 2. Problem description for one dimensional impact example

F igure 3 .

 3 ;-~~-r,-~~-r~r,~-r~~-r~~~-r.-~rT--~~~ . . :. l / :: u \f ' : f : f . One dimensional impact example: (a) contact surface displacement u 1 versus time; (b) contact surface velocity u 1 versus time; (c) contact force A. versus time; (d) centre of mass velocity versus time

  Figure 3 (continued)

Figure 4 .

 4 Figure 4. Illustration of Ith constraint element: (a) configuration at time t,; (b) configuration at time t:+ 1

Figure 5 .

 5 Figure 5. Relative contactor positions: (a) configuration at time t"; (b) configuration at time t:+ 1

  corresponds with completion of step 1 in Algorithm 1. Step 2 in Algorithm 1 is replaced by Algorithms 2 and 3 to solve for An and u,. + 1 .

4 Figure 6 .

 46 Figure 6. Two dimensional impact example: (a) contactor node positions denoted by circles at t" and asterisks at t!+ 1 ;(b) node positions at t"+ 1 for J.L > 0•88; no sliding; (c) node positions at tn + 1 for Jl = 0; no resistance to sliding

  0•700 0•150 0•380 0•640 0•900 0•160 0•440 0•680 0•820 increment. A value of Jl greater than 0•88 leads to this result. In the second case the friction coefficient is zero. The converged results for these two cases are illustrated in Figures6(b) and 6(c), with numerical results for surface contact forces and spatial co-ordinates given in Tables

.

  li+ ll.n -il.n l estimated error = I ; + 11 n I (34)

-Figure 7 .Figure 8 ,

 78 Figure 7. Convergence behaviour for two dimensional impact example

  Figures 9(a) and 9(b) correspond with frictionless contact, Jl = 0. Results for a friction coefficient of 11 = 0•5 are given in Figures 10. TheSDOF solution is calculated using a uniform time step increment of h = 0•000850, which is two orders of magnitude less than the critical time step size.For comparison to the SDOF solutions, finite element results are also presented in Figures9 and 10. The problem description of Figure8(b) is referred to as 'FE object', while that of Figure8(c) is referred to as 'FE object/base'. The finite element solutions are obtained using the Beta-2 explicit integration method with Algorithms 2 and 3.

  SDOF

  FE OBJECT ---------FE OBJECT/BASE•••••••

  ----------------------------------------------------------------------------------~ UJ ()

TIMEFigure 9 .

 9 Figure 9. Two dimensional sliding example. Frictionless contact, Jl. = 0: (a) centre of mass horizontal displacement of sliding object versus time; (b) vertical contact force resultant versus time

  ----------------------------------------------------------------------~

  in,-------------------------------------------------------------------------, BASE................ . ••'~ . 1.\.

TIMEFigure 10 .Figure 10

 1010 Figure 10. Two dimensional sliding example. Friction contact, JL = (}5: (a) centre of mass horizontal displacement of sliding object versus time; (b) vertical contact force resultant versus time; (c) horizontal contact force resultant versus time

  ) and lO(b) corresponds with 7 vertical contact force resultant = L ly 1 I=l (38a) where ly 1 are calculated by Algorithms 2 and 3. Similarly the horizontal contact force resultant for the friction case in Figure lO(c) is 7 horizontal contact force resultant = L A.x 1 I=l (38b)

  for a two dimensional target surface spanning between target nodes A 1 and B 1 which is penetrated by the contactor node C 1 during the time step from t,. to t:+ 1 . The penetration is eliminated by requiring nodes A 1 , B 1 and C 1 to lie on a common line at timet,.+ 1 . This condition may be expressed as Xcr = (1 -a1 )XAI + a1 XBI

	n+l	n+l	n+l	o+ l • + 1	(15)
	where a1 is a dimensionless constraint parameter and XAI , Xs1 and xc, are the spatial n+ 1 n+ 1 >1 + 1 n+ 1
	co-ordinates of nodes AI, B 1 and C I at time t,. + 1 . The constraint element matrix corresponding
	with equation (15) is

  and normal to the I th target surface at time tn + 1 • An elementary set of surface contact force conditions applicable to Art and ANI are, Seidel strategy is easily modified to account for inequalities (32) during iteration. The necessary modifications are given in the Appendix in Algorithm 3, which is a sequence of steps to follow step 2.2.3 of Algorithm 2. The procedure in Algorithm 3 begins in steps 2.2.3a and 2.2.3b by calculating the i + l iterative estimate of the An and ANI. In step 2.2.3c i + 1 ANI is set to zero

		•	•
	•	n
	for I = 1, 2, . . . , N
			(32a)
			(32b)
	Inequality (32a) is a tension limit force condition, while inequality (32b) is a friction limit force
	condition in which p. is a friction coefficient.
	The Gauss-

ft • •

if it is tensile, and i+

1

An is reduced in magnitude if it exceeds the friction limit. Following the modification of i + 1 A!t a~d i + 1

A:r in step 2.2.3c, the corresponding modification to AAx 1 and AAy 1 is calculated in step 2.2.3d before proceeding to step 2.2.4 of Algorithm 2. In the above strategy the displacement constraint element matrices G1 are not altered during n+l iteration. However, enforcement of these constraints is relaxed by the effects of step 2.2.3c. For example, if ANI is zero because of step 2.2.3c, then the Ith contactor node may separate from the • target at time tn + 1 . Similarly, sliding occurs when An is reduced in magnitude as a result of the n friction limit condition in step 2.2.3c. Since the displacement constraints are re-evaluated by equations (17) and (19) at the beginning of each time increment, the constraint parameters a.t will n+l change at each time step when sliding occurs.

Table I .

 I Surface node spatial co-ordinates at time tn and t:+ 1

		Co-ordinates at tn	Co-ordinates at t:+ 1
						Nodal
	Node	Xn	Yn	x:+l	* Yn + l	mass
	1	0•000	0•000	0•000	0•000	2•000
	2	10•000	0•000	10•000	0•000	4•000
	3	20•000	0•000	20•000	0•000	4•000
	4	30•000	0•000	30•000	0•000	2•000
	5	0•000	-2•000	2•000	2•000	1•000
	6	2•000	-2•000	3•000	0•500	1•000
	7	4•000	-2•000	5•000	3•000	1•000
	8	6•000	-2•000	8•000	2•000	1•000
	9	12•000	-2•000	11•000	2•000	1•000
	10	13•000	-2•000	15•000	3•000	1•000
	11	16•000	-2•000	17•000	3•000	1•000
	12	19•000	-2•000	19•000	0•500	1•000
	13	22•000	-2•000	21•000	3•000	1•000
	14	24•000	-2•000	25•000	3•000	1•000
	15	26•000	-2•000	28•000.	3•000	1•000
	16	27•000	-2•000	30•000	3•000	1•000
	A constraint element is activated for each contactor node that penetrates a target segment

during the increment from tn to t: + 1 . In this example 12 constraint elements are activated, corresponding with contactor nodes 5 to 16. The constraint element connectivity and constraint parameters ar , evaluated using equations (19), are given in Table

II

. n + l

Table II .

 II Constraint element connectivity and parameter values

	Element-/	1	2	3	4	5	6	7	8	9	10	11

Table IV .

 IV Forward increment Lagrange multiplier solution for spatial co•ordinates at time t 11 + 1 . No sliding for J-t > 0, no resistance to sliding for J.i = 0

		}-t > 0•88	.u=O
	Node	Xn + l	Yn + l	Xn + l	Y n+ l
	1	0•4410	1•0915	0•0052	1•0939
	2	10•2177	1•0477	9•9880	1•0465
	3	20•1130	1•2646	19•9654	1•2608
	4	30•7083	1•5971	29•9508	1•5732
	5	1•4186	1•0871	1•9957	1•0845
	6	3•0000	0•5000	3•0000	0•5000
	7	4•7427	1•0722	4•9908	1•0702
	8	7•2847	1•0608	7•9955	1•0559
	9	11•7020	1•0803	11•0200	1•0686
	10	13•9779	1•1301	15•0396	1•1550
	11	16•5507	1•1865	17•0387	1•1979
	12	19•0000	0•5000	19•0000	0•5000
	13	21•8082	1•3178	21•0533	1•2949
	14	24•7749	1•4109	25•0494	1•4199
	15	27•3178	1•4907	28•0465	1•5136
	16	28•8012	1•5372	30•0445	1•5761

ACKNOWLEDGEMENT

This work was carried out under a TRW Independent Research and Development Program, project number 89268533, in cooperation with research investigator John E. Crawford.

APPENDIX

Algorithm I. Forward increment Lagrange multiplier method using multi-step central difference integration

Step 1 Given Un-l and Un, predict velocities and displacements 1

Step 2 Calculate contact forces and displacement corrections

Step 3 Calculate velocities and accelerations

Step 4 Next increment: n +-n + 1, return to step 1.

Algorithm 2. Gauss-Seidel iteration algorithm for two dimensional surface contact. This algorithm corresponds with step 2 of Algorithm 1

Step