- Email: [email protected]
Constrained variational approach for dynamic analysis of elastic contact problems
Finite Elements in Analysis and Design 10 (1991) 125-136
125
Elsevier
Constrained variational approach for dynamic analysis of elastic contact problems G y o u n g Jae H u h and Byung M a n K w a k Department of Mechanical Engineering Korea Advanced Institute of Science and Technology, P.O. Box 150, Cheongryang Seoul, South Korea
Received November 1990 Revised June 1991 Abstract. A solution method for dynamic analysis of elastic contact problems with rigid body motion under small deformation is presented. The contact surface is assumed unbonded and frictionless. A variational statement constrained by the geometric compatibility conditions on the contact surface is formulated as the basic principle of the dynamics and its equivalence to the governing equation of equilibrium is shown. Introducing increments in rigid body motion, the variational statement is described in an incremental form and the geometric compatibility conditions are linearized. The .finite element method is adopted as a numerical approximation technique. For the time integration of dynamic response, the displacements are approximated by admissible functions and discontinuities of the velocities due to contact are considered. The resulting discretized system is described as a form of linear complementarity problem, suitable for numerical solution. The formulation is illustrated by means of two numerical examples.
Introduction T h e dynamic contact problems often encountered in m a n y machine parts can cause stress amplification due to impact and lead to a failure of the bodies. It is, however, very difficult to efficiently obtain solutions of d y n a m i c contact problems due to the inherent characteristics that the b o u n d a r y conditions of the bodies under consideration are not k n o w n a priori. The finite element a p p r o a c h has been used for the problems and various descriptions of contact suggested [1-5]. Hughes et al. [1] have p r o p o s e d a scheme taking into account the d y n a m i c c o n t a c t effects. The c o m p o n e n t s of the traction vector at the contact nodes are included as u n k n o w n s and a discretized f o r m of impact and release conditions are used. A s a n o [2] has formulated an approximated hybrid type of virtual work principle b y using Lagrange multipliers and subsidiary conditions for various contact and separate states. Martins and Oden [3] have presented a numerical scheme based on the variational formulation for the d y n a m i c analysis of elastic bodies subjected to friction and prescribed normal tractions. O s m o n t [4] has used contact springs to model the contact and a two-step N e w m a r k algorithm for time integration. C h a u d h a r y and Bathe [5] have employed a Lagrange multiplier a p p r o a c h to enforce the constraints of compatible surface displacements and presented a solution m e t h o d for the d y n a m i c analysis of three-dimensional contact problems with friction. The works cited above are ad hoc in nature in the dealing of contact conditions and the c o n t a c t - i m p a c t p r o b l e m still remains complex. A more systematic formulation with basic principles needs to be developed. In this paper, d y n a m i c analysis of elastic contact problems u n d e r small deformations is considered using a complete description of contact conditions. The contact surface is assumed u n b o n d e d and frictionless. The p r o b l e m is first described b y partial differential equations with 0168-874X/91/$03.50 © 1991 - Elsevier Science Pubfishers B.V. All fights reserved
G.J. Huh, B.M. Kwak /Dynamic analysis of elastic contact problems
126
inequalities. A variational statement constrained by the geometric compatibility conditions on the contact surface is then proposed for numerical analysis of the dynamic contact problems. The equivalence of these two descriptions is shown by considering the necessary conditions of the variational statement. This variational statement is described in an incremental form and the geometric compatibility conditions are linearized. For numerical analysis, the deforming body is discretized by the finite element method, while the displacements are approximated by admissible basis functions over time increments. A possible discontinuity in the velocity due to contact is allowed in the representation. The Lagrange multiplier technique is employed to enforce the geometric compatibility condition of contact. For numerical implementation, several contact check points are selected conveniently at nodal points on the contact surface. Two examples are considered to show the implementations. The first example is the longitudinal impact of two elastic bars for which an analytical solution for the contact force is available. The second example is the impact of an elastic sphere with a rigid plane. The results are compared with quasi-static solutions by the Hertz model.
Statement of the problem For notational simplicity, consider the problem of two elastic bodies 121 and 122 which can impact and contact. Denote by X the position vector of a generic material particle with respect to an inertial coordinate frame (0, x) at time t and X 0 the position vector of the same particle in the non-deformed state at initial time t o. Then,
w = X - Xo,
(a)
where w is the displacement vector of the same particle. The two bodies are governed by the momentum equation, the kinematic relations, and the constitutive equation, which are, respectively
o,~,j + b, = p/~,
in 121 U 122 ..= $2,
, ~ = ~ (1i . W j+wj.,)
in
W12z,
(3)
oij = Cijk,,k*I
in 121U122,
(4)
~'~1
(2)
/
where Cijkl denotes the elastic constitutive coefficient. The velocity, the density and the body force are denoted by fi,(x, t), p and b, respectively. The stress state is described by the Cauchy stress oi~ and the strain tensor c~ is used to measure the deformation. Index notation with the usual summation convention is used. The notation ( ).j denotes differentiation with respect to the coordinate x j and (") denotes twice differentiation with respect to time. The initial conditions are w(x, to) = % ( x ) ,
(5)
if(x, to) = ri0(x),
(6)
where % ( x ) and ri0(x) are the prescribed initial displacement vector and velocity vector, respectively. Following the usual approaches used in contact analysis [6,7], the boundary surface F of each body is composed of three disjoint parts: r w on which displacements are prescribed, F v on which tractions are prescribed, and Fc, the potential contact surface. Displacement boundary conditions for each body are w, =
on rw,
(7)
G.J. Huh, B.M. Kwak / Dynamic analysis of elastic contactproblems
127
y/
~ u n d e f ° r m e d
"--'"~
~-'------
state
g'tX~) = 0
/y.,,\
...
X
g'(x~)= o
L
Fig. 1. Schematic representation of contact condition.
0 xt Fig. 2. Coordinate systems and notation.
where ~ is the prescribed displacement. Traction boundary conditions are
o,;,,j
=
on rF,
(8)
where n = { nj } is the outward unit normal vector on the boundary surface and if, is the prescribed traction. The basic conditions of contact along the contact surface are that no material overlap should occur, and contact forces developed cannot be in tension when the contact surface is assumed unbonded. N o friction is included such that tangential components of the contact force are zero. Then the condition on the normal traction N is,
N=oij ninj ~0, ~tPI
m
ill
(9)
where the superscript m is used to designate quantities associated with the ruth body. The geometric compatibility condition has been generally described by Kwak [6,8]: when the geometry of the potential contact surface of the non-deformed body shown in Fig. 1 is described by g'(X~') ffi 0, (10) then the impenetration compatibility condition states that, at time t,
gt(X~+w2-,,')~O
onr¢.
(11)
Now, since the normal traction must vanish ( N -- 0) whenever the gap opens (gl < 0), or vice versa, the following condition should be satisfied: N-g1=0
o n F c.
(12)
Another geometric compatibility condition using g2 for the potential contact surface of body 2 can be expressed as follows:
g2(X~ +w 1 - w 2)~0
o n F c.
(13)
This condition is the same as condition (11) and so the simpler of the two descriptions can be used.
Variational formulation of the problem In this paper the variational statement [9-11], based on Hamilton's law of varying action, constrained by the geometric compatibility on the contact surface is proposed for the numerical
G.J. Huh, B.M. Kwak /Dynamic analysis of elastic contact problems
128
analysis of the dynamic contact problem. That is, 8 f " J , ( w, ri,)dt +
"to
f0w-8wdI2]/o' = 0
(14)
subject to
g'[X2o+WZ-wl]<~O
o n e c,
(15)
where
Jl(w,
-
b-w) o a - f g . wdr .
(16)
Introducing the Lagrange multiplier p corresponding to inequality constraint eqn. (15), the necessary conditions [12] for eqns. (14) and (15) are obtained as follows: %7.J + bi = p/0i
in ~2,
(17)
o,Tn j = if,.
on F F,
(18)
o,;'n) - p ( V g ' ) i = 0
on Fc,
(19)
0,,'2n,2 + ix(Vg,), = 0
on Fc,
(20)
/~. gl = 0,
(21)
/~ >/0,
(22)
where (Vgl)i denotes the xi-directional component of the gradient of gl(x) at x = X 1. On the contact surface, the outward unit normal vectors of two bodies are opposite to each other at contact. Thus, from eqns. (19), (20) and (9), t~lXTg~ I = - N ,
(23)
where I Vg~] is the magnitude of XTg~ and is positive. The necessary conditions for this variational statement are equal to the governing equations described in the previous section. This variational statement just proposed is very suitable for numerical implementation of the dynamic contact problem, as described below. Again for simplicity of description, a two-body contact problem is considered. Any of the two bodies is allowed to have a rigid body motion. We consider that these bodies are two-dimensional and elastic with rigid body motions. Let P be a generic material particle on the body as shown in Fig. 2. The position vector X of the point P with respect to an inertial coordinate frame (0, x) is given by X = X~ + r,
(24)
where X c is the position vector of the origin G of the moving coordinate system (G, ~) that moves with the body and r is the position vector of the same particle P with respect to the moving coordinate system. If r0 is the vector that defines the position of P in the undeformed state, the position vector r is expressed as follows:
r=ro+v,
(25)
f pi~ d$2" = 0
f o r m = 1, 2,
(26)
f0[(,0 +v)×e] dl2m = 0
for m = 1, 2,
(27)
where v and 6 are the relative displacement and velocity vector, respectively, due to the deformation. Conditions (26) and (27) are used to determine the location G and rotation of the
G.J. Huh, B.M. Kwak / Dynamic analysis of elastic contact problems
129
coordinate system attached to the moving body, such that the rigid body motion and deformation can be decoupled. The problem considered herein deals with small elastic deformations, so the relative displacement(v) due to the deformation can be assumed small, and eqn. (27) can be linearized. From eqns. (3) and (4), the stress tensor (ou) and strain tensor (¢u) are expressed as follows; ',,=
l(avi + av'l
(28)
o, i = C,jk/Ckt,
(29)
fo[(ro+v) xe] d~m---fo[r0xeld~
"
form=l,2.
(30)
The rigid body motion, however, cannot be assumed small. It is assumed that the solution up to a discrete time t~ is known and that the solution for the discrete time t~ + At is required, where At is a suitably chosen time increment. We can define the increments (ua) for the rigid body motion from time t~ to time t~ + At as u a = [ua,, ua 2, ua3] T, (31)
uc,(* ) = X~j(ti + ' r ) - Xcj(ti) = 0(t, +
(32)
for j = 1, 2,
(33)
o(ti),
where • = t - t i, Xaj denotes the x j-directional component of the vector X a and 0 is the rotation of the moving coordinate system with respect to the inertial coordinate system. A time increment is chosen such that the increments in the rigid body motion are assumed small. Let the geometry of the potential contact surface of the non-deformed body at time ti be denoted by smooth functions, gl(al) =0
and
g2(a2) = 0 ,
(34)
where
(35)
a ' = X~'( ti) + r~'( ti).
Then the geometric compatibility conditions can be finearized [7,8], and the variational statement is rearranged as follows; r t + At
8J'
(36)
Jz(uG, k G, v, /~) dt + 8 j 3 ',"+a'= 0
tt
subject to
ga(a2) +
(Vg'),[a
2 u2
,
,
on Fc ,
(37)
f pi~ dI2" = 0
for m = 1, 2,
(38)
fo[r ox~]dI2 m=0
form=l,2,
(39)
where
J2= f[½o, 2
G/ k=l
0 G,
(40)
130
G.J. Huh, B.M. Kwak / Dynamicanalysis of elastic contactproblems
2 8J3 = E (/~k ) TmoSU k ak + f p r . 8v dO, k=l
m~=
m 0
0
I1]
m
12
I1
12
IG
,
= fo(,o
(42)
bodyk
m=fpd ~ k, Ij = fo(i
(41)
× ,o)j dok
(43) for j = 1, 2,
,0) dl2k,
(44)
(45)
and i 3 is the unit normal vector to the x~-x 2 plane. The coefficient matrices a u, flu and ~'u represent the rigid body displacements of points of Fc, F v and ~2 in the /-direction due to a unit displacement in the j t h rigid body degree-of-freedom, respectively, and ti a is the velocity vector of the rigid body motion.
Finite dement discretization and computational method For the numerical analysis of the variational statement, a usual finite element discretization [13] is used. The relative displacement v is approximated by using a nodal displacement vector V: v(~, t) = H(~,) V(t), (46) where H denotes the assembled shape functions. In discretizing condition (37), a finite number of points from the mating surfaces are selected. In general contact check points can be chosen arbitrarily. For the purpose of description, these points are chosen as the node points on the contact surface, and condition (37) is then described as follows: G('r) = Do( ti)U0 + D( t i ) V + Go( ti) <~0,
(47)
DT(ti)Ua = [ V g l ( a } ) ] n • [a~kU2k--Ct~nkulk],
(48)
O f ( t i ) V = V g l ( a 2 ) • [ 02 - v~],
(49)
o0j(t,) =
(50)
where
and Uc denotes the vector that represents the increments in the rigid body motion. In these equations, subscript j is used to designate quantities corresponding to the j t h contact check point, and Daj, Dj are the j t h row vectors of matrices D a, D. In writing eqn. (47), we assume that the shape function on the contact surface is of low order (piecewise constant or linear), and so eqn. (47) is implied by eqn. (37). The discretized variational statement corresponding to eqns. (36)-(39) is obtained as f0 At
(51)
subject to
M,¢'= o,
(52) (53)
G.J. Huh, B.M. Kwak /Dynamic analysis of elastic contact problems
131
where
A = ½VTKV-- ½t)rMoU~- ½~'TM2p-- VrFc- VrF,
(54)
8J5 = (]gMoSV G + I/TM28V ,
(55)
M:= ( ML
(56)
K = f B T C B dO,
(57)
and
Mo=
m~'
M2 = f o n T n dO,
(59)
vgr~ = f b,v,ju~, dO + f F=
f b . H dO +
,a, juo, drF,
f~.n dFF,
(61)
M. = f p H dO, = fp(rXH)
(60)
(62) dO.
(63)
In the above equations, /)o and l;" are the velocity vector of the rigid body motion and the relative velocity vector, respectively, B is the strain-displacement transformation matrix, and C is the stress-strain material property matrix. Suppose that the j t h contact check point comes into contact at time t i ( t ~ < ti <~ti + At) and stays in contact from time t f to time ti + At. The inequality condition (52) corresponding to the j t h contact check point is rearranged in the time domain as follows: gj(¢) < 0
0 < ~"< ,j,
(64)
gj(~j) =0, ~j(¢) = 0
(65) 1-j < ~-~< At,
(66)
where 9 = t f - t~, and gj is the j t h component of the condition (52). Thus, the g j ( , ) may be discontinuous at contact time ~ and a jump in ~j(¢) is non-positive as the point ¢ = ~ is crossed in the positive direction. The discontinuity of the velocity vectors [13] at time ~ due to the jump in ~j(r) is obtained from eqns. (51) and (66) as follows: ~f]G( ~ + O) = If]G ( 7j -- O) -- s j M 01Dgj,
(67)
~'(~ + O) = ~'(~ - O) - s j M f l D f ,
(68)
sj = [D%Mo1D~ + D j M f 1 D f ] - l
j(1"j-0),
(69)
where sj is the jump that shows the discontinuity of the velocity vectors due to the jump in ~j, and is non-negative. D% and Dj are the j t h row vectors of matrices D e and D, respectively. For the numerical analysis of the discretized variational statement in time domain, the nodal displacement vector II(T) is approximated by the admissible basis functions, V(~) = N,(~) V',
(70)
G.J. Huh, B.M. Kwak / D y n a m i c analysis of elastic contact problems
132
where Ni(~- ) denotes a basis function and V ~ is a coefficient vector. In this paper, three basis functions [14] are taken as follows: N, = - ½ 7 ( 1
-,)),
(71)
U 2 = (1 - 71)(1 + ~7),
(72)
U 3 = ½~1(1 + ~1),
(73)
where = ( 2 ~ / a t ) - 1.
(74)
Thus, V l, V 2 and V 3 in eqn. (70) denote the values of V(~) at time ~"= 0, ½At and At, respectively. The variation of the nodal displacement vector satisfying the two given initial conditions (displacement and velocity) can be expressed in the form,
~V= W(n)~V 3,
(75)
where W())) is an arbitrary function and is differentiable with respect to time in order that the variation of the velocity vector can be obtained. The increments in the rigid body motion are approximated by the same basis functions. In the time domain, the time cheek condition (52) is chosen at time ~"= At and the j u m p in gj(~-) is considered at this time in practice. If this j u m p is positive at time ~-= At, this j u m p is considered occurred at time T = 0. These j u m p conditions are expressed as follows:
S° . Gfl = 0,
(76)
G3>~0,
(77)
S° >_-0,
(78)
G3 = ~d(At -- 0)
(79)
where
and S° denotes the j u m p sj at time z = 0. Introducing the Lagrange multipliers P * and ~* corresponding to eqns. (52) and (53), respectively, the necessary conditions for the discretized variational statement eqns. (51)-(53) and the j u m p conditions (76)-(78) at time r = 0 can be summarized as a linear complementarity problem as follows;
[7:]:rQ,,
[:,.,].,
(80)
P . G 3 = 0,
(81)
s o" (~3 = o,
(82)
e >/o,
(83)
S O>~ 0,
(84)
a 3 ~ o,
(85)
~3>/0
fora 3=0,
(86)
S° = 0
forG 3<0,
(87)
G.J. Huh, B.M. Kwak / Dynamic analysis of elastic contact problems
133
where Qll = - ½( At )2 DGMo 1DT -- D M *D T,
(88)
Q13 = - ( A t ) D ~ M o 1DT - ( A t ) D M * [ K * - ~(23' + 1)K] M ; 1 D v,
(89)
Q2~ = (2/at) Q . ,
(90)
Q22 = D c M o 1DT + D M ; ' D r + ( 2 / A t ) Q , 2 ,
(91)
R 1 =DG[U ~
+ (At)U~ + ½(At)2Mo1F~] + G o
+ DM* I F * + ( K * - K ) V I +
( A t ) ( K * - -~(2T + 1)K)IY'],
(92)
R 2 = ( 2 / A t ) [ D G U ~ + D V 1 + Go] - [ 06/.)"~ + DI)"] + ( 2 / A t ) R , ,
(93)
M * = ( K * ) -1 - ( K * ) - 1 M T [ M1( K * ) - 1 M T] - ~ M I ( K * ) -1,
(94)
K* =
#= 3'
½[(2/at)2M2 + (B +
½f) wn(1
= fw ( .
½~ + ~ ) K ] ,
+n)d./~lWd~, l
+ ½) drl/f_~]W drl,
F* = fI_WF dT//f_l' W d~/, FG* =
fI_ WFG d ~ / ~
x= f_'lWX * d,/fLw
(95) (96) (97)
(98)
W d,,
(99)
aT,
(lOO)
P= £, we* d,/ f
(101)
In the above equations, V 1, ¢1, Ud and/)~ denote the values of II, l?, UG and UG at time ¢ = 0, respectively. G 3 denotes the value of G at time ~"= At. A modified simplex method [15] is used to solve this linear complementarity problem. From the solution of P and S °, the nodal displacements at time ~-= At, V 3 and U3, can be calculated as follows: K * V 3 + MTX = F *
+ (K* -K)V'-DTp
+ ( A / ) [ K * - ¼(2), + 1)K] [ l;'1 - M-~D+S°I 2 J, M , V ~ = 0,
MoU2 = Mo[ VJ
(102) (lO3)
+ ( At )UJ] - (At)DTS ° + ½( At)2[ Fc, _ DT p] .
(104)
Numerical solution
Longitudinal impact of two elastic bars
An elastic bar moving with a uniform velocity V0 is considered to impact another elastic bar. Figure 3 shows the finite element mesh using four-node 2-D continuum elements and a table of
G.J. Huh, B.M. Kwak /Dynamic analysis of elastic contact problems
134
()
()
()
()
() () ()
()
©
C)
C)
C)
C:f.,~,
()
()
()
()
(]//
BAR 1
BAR 2
Velocity ( m / s ) Length (m) A r e a ( m 2) Young's modulus (Pa) Poisson's ratio Density ( k g / m 3)
Bar 1
Bar 2
0.01 0.08 0 . 2 x 1 0 -3 0.206 × 1 0 1 2 0.0 0.785 × 104
0.0 0.16 0 . 2 × 1 0 -3 0.206 × 1012 0.0 0.785 × 104
Fig. 3. Finite element mesh and data for longitudinal impact of two bars.
constants used. Numerical data [14] fl and y for the variation of the nodal displacement in the time domain and the time increment At are taken as fl = 0.25, y = 0.5, At = 0.002 ms. The obtained contact force is shown in Fig. 4 and is compared with the corresponding analytical solution [16]. In spite of the coarse mesh, the result shows good agreement with the analytical solution.
I m p a c t o f an elastic sphere against a rigid plane
An elastic sphere moving with a uniform velocity g0 is considered to impact a rigid plane. Figure 5 shows the finite element mesh and the data used. Because of symmetry, the finite element model consists of 4-node 2-D axisymmetric elements. It consists of 72 elements and 95
80.0 o
F.E.M. 1-Dim.
Exoct
Solution v hJ o n~
40.0
¢
z o o o 0.0 0.0
i
6
6
, 0.1
TIME
(rnsec)
Fig. 4. Longitudinal impact of two bars. Contact force vs. time.
\
¢ ¢ / / /
,o
i..-
/ / / / /
p=O.O1 u=0.3 ~=0.25 y=0.5 E=IO00 R=5 At = 0 . 5 X l O -3
/ / / / /
I
V0
/
//////////// Fig. 5. Finite element mesh and data for impact Of an elastic sphere.
G.J. Huh, B.M. Kwak / Dynamic analysis of elastic contact problems
135
Table 1 Maximum contact force and duration of impact for an elastic sphere Velocity (V0)
Maximum contact force
0.1 0.3 3.0
Duration of impact
Hertz
Present
Hertz
Present
4.96 18.52 293.8
5.05 19.50 331.2
0.388 0.311 0.197
0.380 0.303 0.183
nodes. Table 1 shows the maximum contact forces and durations of impact. Figure 6 shows the change of contact force with respect to time. It is observed that the numerically obtained maximum contact forces are larger than that of the quasi-static Hertz's rule solution [17] and that the calculated duration of impact is shorter than the one predicted by the Hertz's static rule. This comparison is presented only to show the general tendency because no other choice is available.
Conclusion
A computational method for dynamic analysis of elastic contact problems with 2 rigid body motion under small deformation has been presented and applications of the method to two-dimensional contact problems have been described. The contact surface is assumed unbonded and frictionless. A variational statement, based on Hamilton's law of varying action, constrained by the geometric compatibility condition on the contact surface has been proposed as a convenient means to numerical analysis of the dynamic contact problems and its equivalence to the governing equation of equilibrium has been shown. This variational statement is described in an incremental form to consider rigid body motion and discretized by usual finite elements in space domain. To enforce the geometric compatibility condition, contact check points are chosen as the node points on the contact surface. In time domain, the displacements are approximated by admissible basis functions and the discontinuities of the velocities due to contact are considered. Although emphasis in the presentation has been on the analysis of two-body contact, the method is directly applicable to multibody contact problems with arbitrarily selected check points. To apply the method to more practical problems, however, frictional conditions must be included with additional complexity.
25.0
hi (..) (Y
400.0
(a)
..,,..
20.0
o •
Vo: Vo=
(b)
0.1 0.3
300.0
o
tY
0 15.0
o
%o
o 0
200.0
(0
~
o
ob -
o
°o o
~10.0 0 (J 5.0
Z
o
o (j 100.0
0.0 0.0
0.1
0.2 0.3 TIME
0.4
0.5
0.0 0.0
o
o ,o
o o
i
0.1
% %1
0.2 TIME
Fig. 6. Impact of an elastic sphere: contact force versus time. (a) Vo = 0.3, 0.1; (b) Vo = 3.0.
0.3
136
G.J. Huh, B.M. Kwak / D y n a m i c analysis of elastic contact problems
References [1] T.J.R. HUGHES, R.L. TAYLOR, J.L. SACKMAN,A. CURNIER and W. KANOKNUKULCHAI,"A finite element method for a class of contact-impact problems", Comput. Methods Appl. Mech. Eng. 8, pp. 249-276, 1976. [2] N. ASANO, "An approximate hybrid type of virtual work principle for two elastoimpact contact bodies", Bull. J S M E 26 (221), pp. 1849-1856, 1983. [3] I.A.C. MARTINS and J.T. ODEN, "A numerical analysis of a class of problems in elastodynamics with friction", Comput. Methods Appl. Mech. Eng. 40, pp. 327-360, 1983. [4] D. OSMONT, "Computation of the dynamic response of structures with unilateral constraints (contact)--Comparison with experimental results", Comput. Methods Appl. Mech. Eng. 34, pp. 847-859, 1982. [5] A.B. CHAUDHARY and K.J. BATHE, "'A solution method for static and dynamic analysis of three-dimensional contact problems with friction", Comput. Struct. 24 (6), pp. 855-873, 1986. [6] B.M. KWAK,"Complementarity problem formulation of three-dimensional frictional contact", J. Appl. Mech., pp. 134-140, 1991. [7] G.B. LEE and B.M. KWAK, "Formulation and implementation of beam contact problems under large displacement by a mathematical programming", Comput Struct. 31 (3), pp. 365-376, 1989. [8] B.M. KWAK, "Numerical implementation of three-dimensional frictional contact by a linear complementary problem", KSME J. 4 (1), pp. 23-31, 1990. [9] O.P. AGRAWAL and S. SAIGAL,"A novel, computationally efficient, approach for Hamilton's law of varying action", Int. J. Mech. Sci. 29 (4), pp. 285-292, 1987. [10] C.D. BAILEY, "Hamilton, Ritz, and elastodynamics", Trans. ASME, J. Appl. Mech., pp. 684-688, 1976. [11] C.D. BAILEY and R.D. WITCHEY, "Harmonic motion of nonconservative, forced, damped systems subjected to nonpotential follower forces", Comput. Methods AppL Mech. Eng. 42, pp. 71-88, 1984. [12] C. LANCZOS, The Variational Principals of Mechanics, Univ. Toronto Press, 4th edn., 1970. [13] K.J. BATHE, Finite Element Procedures in Engineering Analysis, Prentice-Hall, Englewood Cliffs, N J, 1982. [14] O.C. ZIENKmWICZ, The Finite Element Method, McGraw-Hill, New York, 3rd edn., 1977. [15] B.C. LEE and B.M. KWAK, "A computational method for elastoplastic contact problems", Comput. Struct. 18 (5), pp. 757-765, 1984. [16] J.A. ZUKAS,T. NICHOLAS,H.F. SWIFT, L.B. GRESZCZUKand D.R. CURRAN, Impact Dynamics, Wiley, New York, 1982. [17] S.P. TIMOSHENKOand J.N. GOODIER, Theory of Elasticity, McGraw Hill, New York, 3rd edn., 1970.
125
Elsevier
Constrained variational approach for dynamic analysis of elastic contact problems G y o u n g Jae H u h and Byung M a n K w a k Department of Mechanical Engineering Korea Advanced Institute of Science and Technology, P.O. Box 150, Cheongryang Seoul, South Korea
Received November 1990 Revised June 1991 Abstract. A solution method for dynamic analysis of elastic contact problems with rigid body motion under small deformation is presented. The contact surface is assumed unbonded and frictionless. A variational statement constrained by the geometric compatibility conditions on the contact surface is formulated as the basic principle of the dynamics and its equivalence to the governing equation of equilibrium is shown. Introducing increments in rigid body motion, the variational statement is described in an incremental form and the geometric compatibility conditions are linearized. The .finite element method is adopted as a numerical approximation technique. For the time integration of dynamic response, the displacements are approximated by admissible functions and discontinuities of the velocities due to contact are considered. The resulting discretized system is described as a form of linear complementarity problem, suitable for numerical solution. The formulation is illustrated by means of two numerical examples.
Introduction T h e dynamic contact problems often encountered in m a n y machine parts can cause stress amplification due to impact and lead to a failure of the bodies. It is, however, very difficult to efficiently obtain solutions of d y n a m i c contact problems due to the inherent characteristics that the b o u n d a r y conditions of the bodies under consideration are not k n o w n a priori. The finite element a p p r o a c h has been used for the problems and various descriptions of contact suggested [1-5]. Hughes et al. [1] have p r o p o s e d a scheme taking into account the d y n a m i c c o n t a c t effects. The c o m p o n e n t s of the traction vector at the contact nodes are included as u n k n o w n s and a discretized f o r m of impact and release conditions are used. A s a n o [2] has formulated an approximated hybrid type of virtual work principle b y using Lagrange multipliers and subsidiary conditions for various contact and separate states. Martins and Oden [3] have presented a numerical scheme based on the variational formulation for the d y n a m i c analysis of elastic bodies subjected to friction and prescribed normal tractions. O s m o n t [4] has used contact springs to model the contact and a two-step N e w m a r k algorithm for time integration. C h a u d h a r y and Bathe [5] have employed a Lagrange multiplier a p p r o a c h to enforce the constraints of compatible surface displacements and presented a solution m e t h o d for the d y n a m i c analysis of three-dimensional contact problems with friction. The works cited above are ad hoc in nature in the dealing of contact conditions and the c o n t a c t - i m p a c t p r o b l e m still remains complex. A more systematic formulation with basic principles needs to be developed. In this paper, d y n a m i c analysis of elastic contact problems u n d e r small deformations is considered using a complete description of contact conditions. The contact surface is assumed u n b o n d e d and frictionless. The p r o b l e m is first described b y partial differential equations with 0168-874X/91/$03.50 © 1991 - Elsevier Science Pubfishers B.V. All fights reserved
G.J. Huh, B.M. Kwak /Dynamic analysis of elastic contact problems
126
inequalities. A variational statement constrained by the geometric compatibility conditions on the contact surface is then proposed for numerical analysis of the dynamic contact problems. The equivalence of these two descriptions is shown by considering the necessary conditions of the variational statement. This variational statement is described in an incremental form and the geometric compatibility conditions are linearized. For numerical analysis, the deforming body is discretized by the finite element method, while the displacements are approximated by admissible basis functions over time increments. A possible discontinuity in the velocity due to contact is allowed in the representation. The Lagrange multiplier technique is employed to enforce the geometric compatibility condition of contact. For numerical implementation, several contact check points are selected conveniently at nodal points on the contact surface. Two examples are considered to show the implementations. The first example is the longitudinal impact of two elastic bars for which an analytical solution for the contact force is available. The second example is the impact of an elastic sphere with a rigid plane. The results are compared with quasi-static solutions by the Hertz model.
Statement of the problem For notational simplicity, consider the problem of two elastic bodies 121 and 122 which can impact and contact. Denote by X the position vector of a generic material particle with respect to an inertial coordinate frame (0, x) at time t and X 0 the position vector of the same particle in the non-deformed state at initial time t o. Then,
w = X - Xo,
(a)
where w is the displacement vector of the same particle. The two bodies are governed by the momentum equation, the kinematic relations, and the constitutive equation, which are, respectively
o,~,j + b, = p/~,
in 121 U 122 ..= $2,
, ~ = ~ (1i . W j+wj.,)
in
W12z,
(3)
oij = Cijk,,k*I
in 121U122,
(4)
~'~1
(2)
/
where Cijkl denotes the elastic constitutive coefficient. The velocity, the density and the body force are denoted by fi,(x, t), p and b, respectively. The stress state is described by the Cauchy stress oi~ and the strain tensor c~ is used to measure the deformation. Index notation with the usual summation convention is used. The notation ( ).j denotes differentiation with respect to the coordinate x j and (") denotes twice differentiation with respect to time. The initial conditions are w(x, to) = % ( x ) ,
(5)
if(x, to) = ri0(x),
(6)
where % ( x ) and ri0(x) are the prescribed initial displacement vector and velocity vector, respectively. Following the usual approaches used in contact analysis [6,7], the boundary surface F of each body is composed of three disjoint parts: r w on which displacements are prescribed, F v on which tractions are prescribed, and Fc, the potential contact surface. Displacement boundary conditions for each body are w, =
on rw,
(7)
G.J. Huh, B.M. Kwak / Dynamic analysis of elastic contactproblems
127
y/
~ u n d e f ° r m e d
"--'"~
~-'------
state
g'tX~) = 0
/y.,,\
...
X
g'(x~)= o
L
Fig. 1. Schematic representation of contact condition.
0 xt Fig. 2. Coordinate systems and notation.
where ~ is the prescribed displacement. Traction boundary conditions are
o,;,,j
=
on rF,
(8)
where n = { nj } is the outward unit normal vector on the boundary surface and if, is the prescribed traction. The basic conditions of contact along the contact surface are that no material overlap should occur, and contact forces developed cannot be in tension when the contact surface is assumed unbonded. N o friction is included such that tangential components of the contact force are zero. Then the condition on the normal traction N is,
N=oij ninj ~0, ~tPI
m
ill
(9)
where the superscript m is used to designate quantities associated with the ruth body. The geometric compatibility condition has been generally described by Kwak [6,8]: when the geometry of the potential contact surface of the non-deformed body shown in Fig. 1 is described by g'(X~') ffi 0, (10) then the impenetration compatibility condition states that, at time t,
gt(X~+w2-,,')~O
onr¢.
(11)
Now, since the normal traction must vanish ( N -- 0) whenever the gap opens (gl < 0), or vice versa, the following condition should be satisfied: N-g1=0
o n F c.
(12)
Another geometric compatibility condition using g2 for the potential contact surface of body 2 can be expressed as follows:
g2(X~ +w 1 - w 2)~0
o n F c.
(13)
This condition is the same as condition (11) and so the simpler of the two descriptions can be used.
Variational formulation of the problem In this paper the variational statement [9-11], based on Hamilton's law of varying action, constrained by the geometric compatibility on the contact surface is proposed for the numerical
G.J. Huh, B.M. Kwak /Dynamic analysis of elastic contact problems
128
analysis of the dynamic contact problem. That is, 8 f " J , ( w, ri,)dt +
"to
f0w-8wdI2]/o' = 0
(14)
subject to
g'[X2o+WZ-wl]<~O
o n e c,
(15)
where
Jl(w,
-
b-w) o a - f g . wdr .
(16)
Introducing the Lagrange multiplier p corresponding to inequality constraint eqn. (15), the necessary conditions [12] for eqns. (14) and (15) are obtained as follows: %7.J + bi = p/0i
in ~2,
(17)
o,Tn j = if,.
on F F,
(18)
o,;'n) - p ( V g ' ) i = 0
on Fc,
(19)
0,,'2n,2 + ix(Vg,), = 0
on Fc,
(20)
/~. gl = 0,
(21)
/~ >/0,
(22)
where (Vgl)i denotes the xi-directional component of the gradient of gl(x) at x = X 1. On the contact surface, the outward unit normal vectors of two bodies are opposite to each other at contact. Thus, from eqns. (19), (20) and (9), t~lXTg~ I = - N ,
(23)
where I Vg~] is the magnitude of XTg~ and is positive. The necessary conditions for this variational statement are equal to the governing equations described in the previous section. This variational statement just proposed is very suitable for numerical implementation of the dynamic contact problem, as described below. Again for simplicity of description, a two-body contact problem is considered. Any of the two bodies is allowed to have a rigid body motion. We consider that these bodies are two-dimensional and elastic with rigid body motions. Let P be a generic material particle on the body as shown in Fig. 2. The position vector X of the point P with respect to an inertial coordinate frame (0, x) is given by X = X~ + r,
(24)
where X c is the position vector of the origin G of the moving coordinate system (G, ~) that moves with the body and r is the position vector of the same particle P with respect to the moving coordinate system. If r0 is the vector that defines the position of P in the undeformed state, the position vector r is expressed as follows:
r=ro+v,
(25)
f pi~ d$2" = 0
f o r m = 1, 2,
(26)
f0[(,0 +v)×e] dl2m = 0
for m = 1, 2,
(27)
where v and 6 are the relative displacement and velocity vector, respectively, due to the deformation. Conditions (26) and (27) are used to determine the location G and rotation of the
G.J. Huh, B.M. Kwak / Dynamic analysis of elastic contact problems
129
coordinate system attached to the moving body, such that the rigid body motion and deformation can be decoupled. The problem considered herein deals with small elastic deformations, so the relative displacement(v) due to the deformation can be assumed small, and eqn. (27) can be linearized. From eqns. (3) and (4), the stress tensor (ou) and strain tensor (¢u) are expressed as follows; ',,=
l(avi + av'l
(28)
o, i = C,jk/Ckt,
(29)
fo[(ro+v) xe] d~m---fo[r0xeld~
"
form=l,2.
(30)
The rigid body motion, however, cannot be assumed small. It is assumed that the solution up to a discrete time t~ is known and that the solution for the discrete time t~ + At is required, where At is a suitably chosen time increment. We can define the increments (ua) for the rigid body motion from time t~ to time t~ + At as u a = [ua,, ua 2, ua3] T, (31)
uc,(* ) = X~j(ti + ' r ) - Xcj(ti) = 0(t, +
(32)
for j = 1, 2,
(33)
o(ti),
where • = t - t i, Xaj denotes the x j-directional component of the vector X a and 0 is the rotation of the moving coordinate system with respect to the inertial coordinate system. A time increment is chosen such that the increments in the rigid body motion are assumed small. Let the geometry of the potential contact surface of the non-deformed body at time ti be denoted by smooth functions, gl(al) =0
and
g2(a2) = 0 ,
(34)
where
(35)
a ' = X~'( ti) + r~'( ti).
Then the geometric compatibility conditions can be finearized [7,8], and the variational statement is rearranged as follows; r t + At
8J'
(36)
Jz(uG, k G, v, /~) dt + 8 j 3 ',"+a'= 0
tt
subject to
ga(a2) +
(Vg'),[a
2 u2
,
,
on Fc ,
(37)
f pi~ dI2" = 0
for m = 1, 2,
(38)
fo[r ox~]dI2 m=0
form=l,2,
(39)
where
J2= f[½o, 2
G/ k=l
0 G,
(40)
130
G.J. Huh, B.M. Kwak / Dynamicanalysis of elastic contactproblems
2 8J3 = E (/~k ) TmoSU k ak + f p r . 8v dO, k=l
m~=
m 0
0
I1]
m
12
I1
12
IG
,
= fo(,o
(42)
bodyk
m=fpd ~ k, Ij = fo(i
(41)
× ,o)j dok
(43) for j = 1, 2,
,0) dl2k,
(44)
(45)
and i 3 is the unit normal vector to the x~-x 2 plane. The coefficient matrices a u, flu and ~'u represent the rigid body displacements of points of Fc, F v and ~2 in the /-direction due to a unit displacement in the j t h rigid body degree-of-freedom, respectively, and ti a is the velocity vector of the rigid body motion.
Finite dement discretization and computational method For the numerical analysis of the variational statement, a usual finite element discretization [13] is used. The relative displacement v is approximated by using a nodal displacement vector V: v(~, t) = H(~,) V(t), (46) where H denotes the assembled shape functions. In discretizing condition (37), a finite number of points from the mating surfaces are selected. In general contact check points can be chosen arbitrarily. For the purpose of description, these points are chosen as the node points on the contact surface, and condition (37) is then described as follows: G('r) = Do( ti)U0 + D( t i ) V + Go( ti) <~0,
(47)
DT(ti)Ua = [ V g l ( a } ) ] n • [a~kU2k--Ct~nkulk],
(48)
O f ( t i ) V = V g l ( a 2 ) • [ 02 - v~],
(49)
o0j(t,) =
(50)
where
and Uc denotes the vector that represents the increments in the rigid body motion. In these equations, subscript j is used to designate quantities corresponding to the j t h contact check point, and Daj, Dj are the j t h row vectors of matrices D a, D. In writing eqn. (47), we assume that the shape function on the contact surface is of low order (piecewise constant or linear), and so eqn. (47) is implied by eqn. (37). The discretized variational statement corresponding to eqns. (36)-(39) is obtained as f0 At
(51)
subject to
M,¢'= o,
(52) (53)
G.J. Huh, B.M. Kwak /Dynamic analysis of elastic contact problems
131
where
A = ½VTKV-- ½t)rMoU~- ½~'TM2p-- VrFc- VrF,
(54)
8J5 = (]gMoSV G + I/TM28V ,
(55)
M:= ( ML
(56)
K = f B T C B dO,
(57)
and
Mo=
m~'
M2 = f o n T n dO,
(59)
vgr~ = f b,v,ju~, dO + f F=
f b . H dO +
,a, juo, drF,
f~.n dFF,
(61)
M. = f p H dO, = fp(rXH)
(60)
(62) dO.
(63)
In the above equations, /)o and l;" are the velocity vector of the rigid body motion and the relative velocity vector, respectively, B is the strain-displacement transformation matrix, and C is the stress-strain material property matrix. Suppose that the j t h contact check point comes into contact at time t i ( t ~ < ti <~ti + At) and stays in contact from time t f to time ti + At. The inequality condition (52) corresponding to the j t h contact check point is rearranged in the time domain as follows: gj(¢) < 0
0 < ~"< ,j,
(64)
gj(~j) =0, ~j(¢) = 0
(65) 1-j < ~-~< At,
(66)
where 9 = t f - t~, and gj is the j t h component of the condition (52). Thus, the g j ( , ) may be discontinuous at contact time ~ and a jump in ~j(¢) is non-positive as the point ¢ = ~ is crossed in the positive direction. The discontinuity of the velocity vectors [13] at time ~ due to the jump in ~j(r) is obtained from eqns. (51) and (66) as follows: ~f]G( ~ + O) = If]G ( 7j -- O) -- s j M 01Dgj,
(67)
~'(~ + O) = ~'(~ - O) - s j M f l D f ,
(68)
sj = [D%Mo1D~ + D j M f 1 D f ] - l
j(1"j-0),
(69)
where sj is the jump that shows the discontinuity of the velocity vectors due to the jump in ~j, and is non-negative. D% and Dj are the j t h row vectors of matrices D e and D, respectively. For the numerical analysis of the discretized variational statement in time domain, the nodal displacement vector II(T) is approximated by the admissible basis functions, V(~) = N,(~) V',
(70)
G.J. Huh, B.M. Kwak / D y n a m i c analysis of elastic contact problems
132
where Ni(~- ) denotes a basis function and V ~ is a coefficient vector. In this paper, three basis functions [14] are taken as follows: N, = - ½ 7 ( 1
-,)),
(71)
U 2 = (1 - 71)(1 + ~7),
(72)
U 3 = ½~1(1 + ~1),
(73)
where = ( 2 ~ / a t ) - 1.
(74)
Thus, V l, V 2 and V 3 in eqn. (70) denote the values of V(~) at time ~"= 0, ½At and At, respectively. The variation of the nodal displacement vector satisfying the two given initial conditions (displacement and velocity) can be expressed in the form,
~V= W(n)~V 3,
(75)
where W())) is an arbitrary function and is differentiable with respect to time in order that the variation of the velocity vector can be obtained. The increments in the rigid body motion are approximated by the same basis functions. In the time domain, the time cheek condition (52) is chosen at time ~"= At and the j u m p in gj(~-) is considered at this time in practice. If this j u m p is positive at time ~-= At, this j u m p is considered occurred at time T = 0. These j u m p conditions are expressed as follows:
S° . Gfl = 0,
(76)
G3>~0,
(77)
S° >_-0,
(78)
G3 = ~d(At -- 0)
(79)
where
and S° denotes the j u m p sj at time z = 0. Introducing the Lagrange multipliers P * and ~* corresponding to eqns. (52) and (53), respectively, the necessary conditions for the discretized variational statement eqns. (51)-(53) and the j u m p conditions (76)-(78) at time r = 0 can be summarized as a linear complementarity problem as follows;
[7:]:rQ,,
[:,.,].,
(80)
P . G 3 = 0,
(81)
s o" (~3 = o,
(82)
e >/o,
(83)
S O>~ 0,
(84)
a 3 ~ o,
(85)
~3>/0
fora 3=0,
(86)
S° = 0
forG 3<0,
(87)
G.J. Huh, B.M. Kwak / Dynamic analysis of elastic contact problems
133
where Qll = - ½( At )2 DGMo 1DT -- D M *D T,
(88)
Q13 = - ( A t ) D ~ M o 1DT - ( A t ) D M * [ K * - ~(23' + 1)K] M ; 1 D v,
(89)
Q2~ = (2/at) Q . ,
(90)
Q22 = D c M o 1DT + D M ; ' D r + ( 2 / A t ) Q , 2 ,
(91)
R 1 =DG[U ~
+ (At)U~ + ½(At)2Mo1F~] + G o
+ DM* I F * + ( K * - K ) V I +
( A t ) ( K * - -~(2T + 1)K)IY'],
(92)
R 2 = ( 2 / A t ) [ D G U ~ + D V 1 + Go] - [ 06/.)"~ + DI)"] + ( 2 / A t ) R , ,
(93)
M * = ( K * ) -1 - ( K * ) - 1 M T [ M1( K * ) - 1 M T] - ~ M I ( K * ) -1,
(94)
K* =
#= 3'
½[(2/at)2M2 + (B +
½f) wn(1
= fw ( .
½~ + ~ ) K ] ,
+n)d./~lWd~, l
+ ½) drl/f_~]W drl,
F* = fI_WF dT//f_l' W d~/, FG* =
fI_ WFG d ~ / ~
x= f_'lWX * d,/fLw
(95) (96) (97)
(98)
W d,,
(99)
aT,
(lOO)
P= £, we* d,/ f
(101)
In the above equations, V 1, ¢1, Ud and/)~ denote the values of II, l?, UG and UG at time ¢ = 0, respectively. G 3 denotes the value of G at time ~"= At. A modified simplex method [15] is used to solve this linear complementarity problem. From the solution of P and S °, the nodal displacements at time ~-= At, V 3 and U3, can be calculated as follows: K * V 3 + MTX = F *
+ (K* -K)V'-DTp
+ ( A / ) [ K * - ¼(2), + 1)K] [ l;'1 - M-~D+S°I 2 J, M , V ~ = 0,
MoU2 = Mo[ VJ
(102) (lO3)
+ ( At )UJ] - (At)DTS ° + ½( At)2[ Fc, _ DT p] .
(104)
Numerical solution
Longitudinal impact of two elastic bars
An elastic bar moving with a uniform velocity V0 is considered to impact another elastic bar. Figure 3 shows the finite element mesh using four-node 2-D continuum elements and a table of
G.J. Huh, B.M. Kwak /Dynamic analysis of elastic contact problems
134
()
()
()
()
() () ()
()
©
C)
C)
C)
C:f.,~,
()
()
()
()
(]//
BAR 1
BAR 2
Velocity ( m / s ) Length (m) A r e a ( m 2) Young's modulus (Pa) Poisson's ratio Density ( k g / m 3)
Bar 1
Bar 2
0.01 0.08 0 . 2 x 1 0 -3 0.206 × 1 0 1 2 0.0 0.785 × 104
0.0 0.16 0 . 2 × 1 0 -3 0.206 × 1012 0.0 0.785 × 104
Fig. 3. Finite element mesh and data for longitudinal impact of two bars.
constants used. Numerical data [14] fl and y for the variation of the nodal displacement in the time domain and the time increment At are taken as fl = 0.25, y = 0.5, At = 0.002 ms. The obtained contact force is shown in Fig. 4 and is compared with the corresponding analytical solution [16]. In spite of the coarse mesh, the result shows good agreement with the analytical solution.
I m p a c t o f an elastic sphere against a rigid plane
An elastic sphere moving with a uniform velocity g0 is considered to impact a rigid plane. Figure 5 shows the finite element mesh and the data used. Because of symmetry, the finite element model consists of 4-node 2-D axisymmetric elements. It consists of 72 elements and 95
80.0 o
F.E.M. 1-Dim.
Exoct
Solution v hJ o n~
40.0
¢
z o o o 0.0 0.0
i
6
6
, 0.1
TIME
(rnsec)
Fig. 4. Longitudinal impact of two bars. Contact force vs. time.
\
¢ ¢ / / /
,o
i..-
/ / / / /
p=O.O1 u=0.3 ~=0.25 y=0.5 E=IO00 R=5 At = 0 . 5 X l O -3
/ / / / /
I
V0
/
//////////// Fig. 5. Finite element mesh and data for impact Of an elastic sphere.
G.J. Huh, B.M. Kwak / Dynamic analysis of elastic contact problems
135
Table 1 Maximum contact force and duration of impact for an elastic sphere Velocity (V0)
Maximum contact force
0.1 0.3 3.0
Duration of impact
Hertz
Present
Hertz
Present
4.96 18.52 293.8
5.05 19.50 331.2
0.388 0.311 0.197
0.380 0.303 0.183
nodes. Table 1 shows the maximum contact forces and durations of impact. Figure 6 shows the change of contact force with respect to time. It is observed that the numerically obtained maximum contact forces are larger than that of the quasi-static Hertz's rule solution [17] and that the calculated duration of impact is shorter than the one predicted by the Hertz's static rule. This comparison is presented only to show the general tendency because no other choice is available.
Conclusion
A computational method for dynamic analysis of elastic contact problems with 2 rigid body motion under small deformation has been presented and applications of the method to two-dimensional contact problems have been described. The contact surface is assumed unbonded and frictionless. A variational statement, based on Hamilton's law of varying action, constrained by the geometric compatibility condition on the contact surface has been proposed as a convenient means to numerical analysis of the dynamic contact problems and its equivalence to the governing equation of equilibrium has been shown. This variational statement is described in an incremental form to consider rigid body motion and discretized by usual finite elements in space domain. To enforce the geometric compatibility condition, contact check points are chosen as the node points on the contact surface. In time domain, the displacements are approximated by admissible basis functions and the discontinuities of the velocities due to contact are considered. Although emphasis in the presentation has been on the analysis of two-body contact, the method is directly applicable to multibody contact problems with arbitrarily selected check points. To apply the method to more practical problems, however, frictional conditions must be included with additional complexity.
25.0
hi (..) (Y
400.0
(a)
..,,..
20.0
o •
Vo: Vo=
(b)
0.1 0.3
300.0
o
tY
0 15.0
o
%o
o 0
200.0
(0
~
o
ob -
o
°o o
~10.0 0 (J 5.0
Z
o
o (j 100.0
0.0 0.0
0.1
0.2 0.3 TIME
0.4
0.5
0.0 0.0
o
o ,o
o o
i
0.1
% %1
0.2 TIME
Fig. 6. Impact of an elastic sphere: contact force versus time. (a) Vo = 0.3, 0.1; (b) Vo = 3.0.
0.3
136
G.J. Huh, B.M. Kwak / D y n a m i c analysis of elastic contact problems
References [1] T.J.R. HUGHES, R.L. TAYLOR, J.L. SACKMAN,A. CURNIER and W. KANOKNUKULCHAI,"A finite element method for a class of contact-impact problems", Comput. Methods Appl. Mech. Eng. 8, pp. 249-276, 1976. [2] N. ASANO, "An approximate hybrid type of virtual work principle for two elastoimpact contact bodies", Bull. J S M E 26 (221), pp. 1849-1856, 1983. [3] I.A.C. MARTINS and J.T. ODEN, "A numerical analysis of a class of problems in elastodynamics with friction", Comput. Methods Appl. Mech. Eng. 40, pp. 327-360, 1983. [4] D. OSMONT, "Computation of the dynamic response of structures with unilateral constraints (contact)--Comparison with experimental results", Comput. Methods Appl. Mech. Eng. 34, pp. 847-859, 1982. [5] A.B. CHAUDHARY and K.J. BATHE, "'A solution method for static and dynamic analysis of three-dimensional contact problems with friction", Comput. Struct. 24 (6), pp. 855-873, 1986. [6] B.M. KWAK,"Complementarity problem formulation of three-dimensional frictional contact", J. Appl. Mech., pp. 134-140, 1991. [7] G.B. LEE and B.M. KWAK, "Formulation and implementation of beam contact problems under large displacement by a mathematical programming", Comput Struct. 31 (3), pp. 365-376, 1989. [8] B.M. KWAK, "Numerical implementation of three-dimensional frictional contact by a linear complementary problem", KSME J. 4 (1), pp. 23-31, 1990. [9] O.P. AGRAWAL and S. SAIGAL,"A novel, computationally efficient, approach for Hamilton's law of varying action", Int. J. Mech. Sci. 29 (4), pp. 285-292, 1987. [10] C.D. BAILEY, "Hamilton, Ritz, and elastodynamics", Trans. ASME, J. Appl. Mech., pp. 684-688, 1976. [11] C.D. BAILEY and R.D. WITCHEY, "Harmonic motion of nonconservative, forced, damped systems subjected to nonpotential follower forces", Comput. Methods AppL Mech. Eng. 42, pp. 71-88, 1984. [12] C. LANCZOS, The Variational Principals of Mechanics, Univ. Toronto Press, 4th edn., 1970. [13] K.J. BATHE, Finite Element Procedures in Engineering Analysis, Prentice-Hall, Englewood Cliffs, N J, 1982. [14] O.C. ZIENKmWICZ, The Finite Element Method, McGraw-Hill, New York, 3rd edn., 1977. [15] B.C. LEE and B.M. KWAK, "A computational method for elastoplastic contact problems", Comput. Struct. 18 (5), pp. 757-765, 1984. [16] J.A. ZUKAS,T. NICHOLAS,H.F. SWIFT, L.B. GRESZCZUKand D.R. CURRAN, Impact Dynamics, Wiley, New York, 1982. [17] S.P. TIMOSHENKOand J.N. GOODIER, Theory of Elasticity, McGraw Hill, New York, 3rd edn., 1970.