Contact force algorithm in explicit transient analysis using finite-element method

Finite Elements in Analysis and Design 43 (2007) 580 – 587 www.elsevier.com/locate/finel

Contact force algorithm in explicit transient analysis using finite-element method Fu-Jun Wang a,∗ , Li-Ping Wang a , Jian-Gang Cheng b , Zhen-Han Yao b a College of Water Conservancy and Civil Engineering, China Agricultural University, Beijing 100083, China b Department of Engineering Mechanics, Tsinghua University, Beijing 100084, China

Received 6 May 2006; received in revised form 21 December 2006; accepted 30 December 2006 Available online 20 February 2007

Abstract The contact force calculation has significant effect on the accuracy and efficiency of finite-element analysis for contact–impact problems. In this paper, an algorithm for contact force calculation in explicit finite-element analysis is presented. The classical penalty method and Lagrange multiplier method are combined by interpreting the impenetrability between two bodies as a control parameter in this algorithm. The constraints of no-penetration condition and sticking condition are imposed accurately without the requirement to choose a proper parameter as in the penalty method. With this algorithm, no coupled system of equations needs to be solved. The interaction of contact pairs applied on same segment is treated with iteration procedure. This algorithm is totally compatible with the explicit time integration scheme. Numerical examples demonstrate that this algorithm performs well for contact–impact problems. In contrast to the penalty method, this algorithm increases the stable time step because the effect of penalty parameter to the stability does not exist. The new algorithm provides an accurate and efficient way to study contact–impact problems with explicit finite-element method. 䉷 2007 Elsevier B.V. All rights reserved. Keywords: Contact force algorithm; Finite-element method; Explicit transient analysis; Contact–impact problem

1. Introduction Dynamic contact–impact problems have received wide attention due to their universality and importance in many engineering fields, e.g. the automobile industry, power industry, etc. The computations of the dynamic impact and contact of interacting bodies are notoriously difficult, due to the strong nonlinearity of the associated system. The usual numerical methods for dealing with this kind of problems are based on the finiteelement method. In transient contact–impact analysis, the time integration of the equilibrium equations is usually performed using some form of an explicit scheme since high frequencies are involved [1,2]. In the aspects of dealing with contact constraints, the existing methods can be mainly classified into two major branches. One is the penalty methods [3,4], which allow ∗ Corresponding author. Tel./fax: +86 10 62736972.

E-mail addresses: [email protected], [email protected] (F.-J. Wang). 0168-874X/$ - see front matter 䉷 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.finel.2006.12.010

penetration to occur but penalize it by applying surface contact force models. Another is the Lagrange multiplier methods [3,5,6] which exactly preserve the non-interpenetration constraint. The penalty approach satisfies contact conditions approximately and the accuracy of approximate solutions depends strongly on the penalty parameter. This method is widely used in complex three-dimensional contact–impact problems since it is simple to use in explicit finite-element system. However, there are no clear rules to choose the penalty parameter; it depends on particular problems considered [7]. On the other hand, the penalty method will affect the stability of the explicit analysis, which is conditionally stable, when the penalty parameter reaches a certain value [8]. Unlike the penalty method, the Lagrange multiplier method does not use the algorithmic parameter and it enforces the zeropenetration condition exactly. Thus, this method can give out very accurate displacement fields in the analysis of static contact problems [8]. However, for dynamic contact problems, this

F.-J. Wang et al. / Finite Elements in Analysis and Design 43 (2007) 580 – 587

method requires the solution of implicit augmented systems of equations, which can become computationally very expensive for large problems [1]. The first application for dynamic analysis was presented by Hughes et al. [5]. Chaudhary and Bathe [9] developed a general Lagrange multiplier method for dealing with three-dimensional static and dynamic contact problem. In their method, the combination of the Lagrange multiplier method and the implicit Newmark integration algorithm has been used. To make it useful in explicit analysis, some improved approaches [1,6,10,11] have been proposed. Although the first two methods [6,11] are consistent with explicit scheme, they are yet time consuming [12]. In the last two methods [1,10], the term Lagrange multiplier method is usually used to imply that the impenetrability condition is satisfied, but further assumptions are involved to update the global position of the contact surface and to determine the tractions, velocities and accelerations over the contact surface [7]. In this paper, the classical penalty method and Lagrange multiplier method are combined to make a new contact algorithm for the contact force computation in explicit analysis. In this method, both no-penetration condition and interaction of neighbor contact points are considered. This algorithm is totally compatible with explicit time integration scheme. The contact constraint condition and the equilibrium equation can be satisfied sufficiently in each time step. It provides an accurate and efficient way to study contact–impact problems with explicit finite-element method as compared to penalty method. The outline of this paper is as follows. In Section 2, we summarize the features of the contact problem and finite-element procedure with explicit time integration. In Section 3 we describe the new algorithm in detail. The normal and tangential contact forces and their computing procedure will be given. In Section 4, we examine the algorithm by two examples to demonstrate the performance of the algorithm. 2. Basic description of the contact–impact problem 2.1. General definitions of the contact–impact problem Consider two continuous bodies (with volumes V 1 and V 2 and surfaces S 1 and S 2 , respectively) which come into contact with each other at some instant of time t. The algorithm to be developed is applicable to any number of bodies but we restrict formal development to two bodies for simplicity of notation. The contact problem is defined as an initial-boundary problem where each body obeys the condition of contact equilibrium and constitutive equations of the material. According to the principles of continuum mechanics, the contact conditions can be described as follows: (1) Geometrical condition: At time t, the surfaces on two contact bodies should be the impenetrability condition: V1 ∩V2 =0

and

S 1 ∩ S 2  = 0.

(1)

For a point on the contact surface (i.e. common boundary), the position vectors (x1 and x2 ) must be the same, and the

581

outside surface unit normal vectors (n1 and n2 ) in respective to the two body surfaces must be in the opposite direction: x1 = x2

and

n1 = −n2 .

(2)

(2) Kinematics condition: This condition restrains the velocities (u˙ 1 and u˙ 2 ) of the points on the common surface of the contact bodies to satisfy the impenetrability: u˙ 1 · n1 + u˙ 2 · n2 = 0.

(3)

(3) Kinetic condition: The condition for the balance of momentum on the boundary requires: f1 + f2 = 0,

(4)

where f is the Cauchy traction acting on the contact boundary. In the same time, no tensile normal force on the contact boundary should be obeyed: f1 · n1 0

and

f2 · n2 0.

(5)

2.2. Finite-element formulations For dynamic contact problem, the governing equations of the system motion are obtained using the principle of virtual work as [13,14]:    t T Tt [t e]  dV = [ u] b dV + [t u]T T dS tV

0V

0S



−

tV

 +

tS c

[t u]T t u¨ dV ([t uc ]Tt fc ) dS,

(6)

where t e is the infinitesimal strain,  is the Cauchy stress, b is the body force per unit volume, T is the applied external force over the surface of the body,  is the mass density, t u¨ is the inertial force, t u is the virtual displacement, fc is the contact force resulted on the contact surface Sc , uc is the gap between the two bodies on the contact surface, the superscripts t and 0 represent the time t and time 0, respectively, and superscript T represents transpose of a vector. The fourth item on the right side of the equation denotes the function associated with the contact boundary condition. By using the finite-element discretization procedure, the discretized equation of motion to be solved in explicit finiteelement analysis may be put into the following form [13]: Mt u¨ = t Q − t F + t Rc ,

(7)

where M is the mass matrix, u¨ is the acceleration vector, Q is the prescribed external load vector, F is the internal force vector and Rc is the force vector associated with contact forces. The derivation of contact force vector Rc will be discussed later in Section 3. Integrating Eq. (7) over time by the central difference method, the displacements at time t + t are obtained according to their values at time t − t and t as follows: t+t

u = M−1 [t 2 (t Q − t F + t Rc ) + 2Mt u − Mt−t u],

(8)

582

F.-J. Wang et al. / Finite Elements in Analysis and Design 43 (2007) 580 – 587

where t is the time increment which should be less than that given by Courant criterion in order to ensure the stability of calculation. With a diagonal mass matrix and using a method to explicitly determine the contact force, the equations associated with contact bodies are then uncoupled with each other, Eq. (8) can be written for each degree of freedom.

where  and  are the isoparametric coordinates defined for master segment, j and j take on their nodal values at j

(±1, ±1), and xi is the nodal coordinate of the jth node in the ith direction, see Fig. 1. Note that r is at least once continuously differential and that the normal is non-zero: jr jr ×  = 0, j j

(14)

3. Contact force algorithm (c , c ) are the contact point coordinates on master segment M. They can be determined by solving Eqs. (9), (10), (11). The details can be found in [16].

3.1. Determination of the contact point To study the interaction mechanism and contact force computation, let us take a slave node and a master element as an example of contact pair, as shown in Fig. 1. Assume that slave node is not identified as lying on the intersection of two master segments. Consequently, the contact point C is defined on the master segment M as the point closest to the slave node S. The contact point is calculated by directly solving the minimization problem. The function J to be minimized is defined as [15]: Min: J = {[t − r(c, c )][t − r(c , c )]}1/2 with

(9)

jJ jr = [t − r(c, c )] (c , c ) = 0, j j

(10)

jJ jr = [t − r(c, c )] (c , c ) = 0, j j

(11)

where t is the position vector drawn to slave node S, r represents the master segment M. Each master segment can be described by a bilinear parametric representation [3,12]: r = f1 (, )i1 + f2 (, )i2 + f3 (, )i3 , fi (, ) =

4 

j

j x i ,

j =1

(12)

1 j (, ) = (1 + j )(1 + j ), 4 (13) 

N4

master segment S M

n

N1

r

N3

 ∂r ( , ) c c ∂

C

t x3

∂r ( , ) c c ∂

N2

x2 x1 Fig. 1. Slave node S contacts with master segment M.

3.2. Normal nodal force The normal direction at the contact point can be calculated either at the contact point on the master surface or at the slave node. In general, the two choices will yield different directions. Here, the normal at the contact point on the master segment is taken as the normal direction for the contact force determination. This direction is calculated by taking the cross product of the derivative with respect to the isoparametric coordinates: n=

jr/j × jr/j . |jr/j × jr/j|

(15)

Analogize to the conventional penalty method. Each slave node is checked for penetration through its master segment at this method. If the slave node does not penetrate, nothing is done, but if it does, an interface force as described below is applied between the slave node and its contact point. Consider the motions of slave node S and contact point C in the normal direction. The equations of motion of the two nodes are described by 

M t u¨ n = t Fˆ n + t Rcn

( = I, II),

(16)

where subscript n denotes a component in the normal direction, left superscript t denotes the time t, superscript  could be I or II, I denotes that the associated variable is defined at nodal S, and II at contact point C, M is the nodal mass, u¨ is the acceleration of the node, Fˆ is the nodal force without contact force, Rc is the nodal contact force. Rc can be decomposed into two parts: t

Rcn = R¯ cn + fn ,

(17)

where f is the nodal contact force increment contributed by this contact pair, R¯ c is the nodal contact force contributed by other contact pairs. In Eqs. (16) and (17), fnII (or fnI ) is to be found according to contact conditions, and M II is to be constructed according the contact point position and its master segment M.

F.-J. Wang et al. / Finite Elements in Analysis and Design 43 (2007) 580 – 587

To evaluate M II , here we use the defence node algorithm developed by Zhong [13]. We suggest the following: u¨ II n =

n 

Nk u¨ nk ,

(18)

k=1

M II

n 

Nk fnk /Mk = fnII ,

fnk = k fnII

(k = 1, n),

fnk = fnII ,

k=1

n 

k = 1,

Nk ( k M /Mk ) = 1.

(30)

t

u¨ =

1 t−t ( u − 2t u + t+t u), t 2

(31)

(21)

t

u˙ =

1 t+t ( u − t−t u). 2t

(32)

II

(22)

k = Mk /M II

(23)

However, k may generally be written as [13]:   n  l−1 II l M (Nk ) , k = Mk (Nk )

(24)

where l is an integer equal to or greater than 1. When l = 1, Eq. (24) becomes Eq. (23). The more reasonable choice for k is l = 2. Thus, we obtain:   n  II 2 M (25) (Nk ) . k = Mk Nk k=1

Substituting Eq. (25) into the second item of Eq. (21), we obtain:   n n   M II (26) M k Nk (Nk )2 = 1, k=1

which gives n 



M k Nk

k=1

n 

(Nk )2 .

(27)

k=1

mk ,

k=1

mk = Mk Nk



n 

 (Nk )

2

.

(28)

k=1

If the contact point is located at the segment node k with mass Mk , we have M II = Mk .

= (2t un − t−t un ) +

(t)2 t ˆ   ( Fn + R¯ cn + fn ) M

(33)

or t+t  un

(t)2 t ˆ   = t t−t U˙ n + ( Fn + R¯ cn + fn ), M

(34)

where 1 U˙ = (2t u − t−t u). t

(35)

On the other hand, the no-penetration condition gives that t+t I gn = t gn + t+t uII un = 0, n −

(36)

where gn is the gap in normal direction between slave node S and contact point C. This condition will be used to construct the new contact algorithm as implemented in Lagrange multiplier method. Substituting Eq. (34) into (35), and considering that the nodal contact forces increment at node S and contact point C satisfy Eq. (4), we obtain:  t Fˆ II + R t Fˆ I + R II ¯ cn ¯I M I M II n I fn = I − n I cn II II M +M M M  t−t U t−t U tg ˙ nII ˙ nI n , (37) + − − t t (t)2 where fnI is the normal contact force increment on slave node S, which is contribution of the current contact pair. It can be described by vector form as



This can be written as n 

t+t  un

t+t

k=1

k=1

Substituting Eqs. (31) and (32) into Eq. (16), we obtain:

t−t

One choice of k is

M =

Nk Fˆnk .

(20)

k=1

II

n 

To integrate Eq. (16), we adopt a central difference approximation [17]. The accelerations and the velocities are written as

k=1

M II =

FˆnII =

(19)

where n is the total number of the master segment nodes, Nk , Mk , are shape function and mass of node k, respectively. fnk denotes the part of the normal contact force increment contributed to node k by the contact point. k is a function to be chosen. Substituting Eq. (20) into Eq. (19), we have n 

For the calculation of FˆnII in Eq. (17), the same method adopted in Eq. (18) is used here. When the Fˆnk is known, we have

k=1

k=1

n 

583

(29)

fIn = |fnI |n.

(38)

The normal contact force increment on contact point C can be determined by fnII = −fnI . The normal contact force increment on contact point, fnII , is then distributed over the master segment nodes according to Eq. (20) and (28), i.e.: fnk =

mk II f . M II n

(39)

The complete normal contact force in which other contact pairs are considered will be discussed in Section 3.4.

584

F.-J. Wang et al. / Finite Elements in Analysis and Design 43 (2007) 580 – 587

3.3. Tangential contact force

performed. With an explicit time integration method employed, the following procedure may be used:

Assume that slave node S sticks with master segment M on contact point C, also as shown in Fig. 1. The tangential contact force, or friction, can be determined according to the similar procedure to the normal contact force computation: ftI

M I M II = I M + M II +

ftI

t−t U ˙ II t

t

M I M II = I M + M II +

t−t U ˙ tII

t

t

II FˆtII + R¯ ct 

M II − t

t−t U ˙I t

t

−

II FˆtII + R¯ ct 

M II −

t−t U ˙ tI

t

−

−

I + R¯ ct 

t Fˆ I t

MI

tg

 t 2

,

(t) −

I + R¯ ct 

t Fˆ I ty

tg

(40)

MI  t 2

(t)

,

(41)

where subscripts  and  denote the tangential directions with isoparametric coordinates defined on a master segment. ftI and ftI are the tangential components of contact force increment of slave node S. Let us define a slip criterion: = |ft | − |fn | =



ft + ft − |fn |,

(42)

where is the dynamic frictional coefficient. For discussion and determination of , see [4]. Then, the vector of actual tangential contact force increment on slave node S is determined by ⎧ jr jr ⎨f I ( ,  ) + ftI (c , c ) if 0, t  j c c j ft = ⎩

|fn |nt if > 0,

(43)

where nt is the tangential direction which is determined by nt =

 jr jr ftI (c , c ) + ftI (c , c ) j j    I jr  I jr f  ,  ) + f ,  ) ( ( c c . t  t j c c j

(44)

3.4. Procedure of contact force computation Since one master segment may contact with more than one slave nodes, adjacent contact pairs may exist. The nodal contact force caused by a contact pair will affect the contact force computation of other contact pairs when the above contact algorithm is used. This is why the contact force caused by one contact pair is called the contact force increment. Therefore, an iterative procedure has to be performed after contact force increment is determined by the above method in order to consider the interaction of adjacent contact pairs. Assume that all relevant quantities are known at time t − t and t, and contact force t Rc at time t + t in Eq. (8) is to be

1. Find all contact pairs in the entire computational domain by calling global contact searching algorithms, e.g. HITA in [13] and algorithm in [18]. 2. Find all slave nodes, i.e. all contact pairs, for each master segment by calling local contact searching algorithms, e.g. FFS algorithm in [16]. 3. Set j = 1. 4. For each contact pair, calculate contact force increment fn and ft according to (37)–(44). 5. For each master segment in the domain, loop over its all contact pairs, and set (R¯ cn )(j ) = (R¯ cn )(j ) + fn , (R¯ ct )(j ) = ¯ ct )(j ) + ft . Thus, the contact forces have been distributed (R ¯ cn has onto corresponding nodes. It should be noted that R to be set to 0 if its value is larger than 0 because no tensile normal force has to be obeyed. ¯ Ic )(j ) − 6. For each slave node in the domain, check |(R (j −1) (j ) I I ¯ ¯ (Rc ) |/|(Rc ) | < . If not, j = j + 1, and then goto step 4. Otherwise terminate. Here,  is a control factor that is preassigned with a value of 0.05. It should be mentioned that the maximum iteration counter jmax should be set up to avoid deadlock of iteration. Usually, jmax = 5 can make the results accurate enough. 4. Numerical examples 4.1. Impact of a square-section beam against a rigid wall The square-section beam is a typical example of crashworthiness analysis. As the beam is impacted on the rigid wall, buckling and severe folding take place and extensive self-contacting areas are formed. The large mass at the beam end contributes to the continuous buckling of the beam. The system is shown in Fig. 2. The beam has an initial velocity of 15.6 m/s and is free from any external loading. A lumped mass of 1440 kg is attached to the right end of the beam to simulate crash loading. Due to the symmetry, only one quarter of the beam is modeled. This problem has been investigated with different contact algorithms by several researchers, including Benson and Hallquist [15], Zhong [13], Belytschko [1], Bathe [19] and the authors [16]. This example was performed to demonstrate the capability of the new algorithm to simulate the contact of shell elements. The contact analysis was performed by using the new algorithm and the classical penalty method, respectively. The global

Rigid Wall

Attached Mass, M Beam v

Fig. 2. Impact of a square-section beam against a rigid wall.

F.-J. Wang et al. / Finite Elements in Analysis and Design 43 (2007) 580 – 587

contact search and local contact search are performed by using the author’s previous algorithms [16,18,20]. Fig. 3 shows the deformed configurations of the beam at time 2.16 ms. Fig. 3(a) is the result with a fixed time step of t =0.4×10−5 s, which is obtained with the new algorithm. Fig. 3(b) presents the result of the penalty method while the same other settings are adopted as the new algorithm. It is clear that a penetration was found for the penalty method. The penalty method can only give acceptable results when small penalty parameter and small time step size are used. This demonstrates that high accuracy can be obtained by using the new algorithm.

585

55 nodes of the middle right of the second transversal beam to simulate the engine loading. The impact–contact problem was analyzed by using the present contact algorithm. The deformed configurations of frontal part of the vehicle frame at time 12, 30 and 52 ms are shown in Fig. 6. Fig. 7 gives the longitudinal permanent deformation (i.e. the longitudinal displacement) of the frame end point A (indicated in Fig. 5) as compared to experimental result [12]. The problem was also analyzed by using penalty method. The numerical results show that the present algorithm increases the efficiency of total computation by 20% as compared to the penalty method when a stable integration and the similar accuracy are obtained.

4.2. Impact of vehicle frame against a rigid wall 5. Conclusions and discussions This example concerns the impact between a vehicle and a rigid wall. The crashworthiness experiment system is shown in Fig. 4. Here, we focus our study on the vehicle frame because it is the main element for deformation and energy absorption of the vehicle. The finite-element model used for the vehicle frame is as shown in Fig. 5. The vehicle frame has an initial velocity of 48 km/h. A lumped mass of 550 kg is distributed on

a

b

Analyzed with the new algorithm.

Analyzed with the penalty method.

Fig. 3. The deformed configurations at time 2.16 ms.

A new contact algorithm for contact force calculation in explicit finite-element analysis has been discussed. It can be viewed as a combination of the classical Lagrange multiplier method and the penalty method. The constraints of no-penetration condition and sticking condition are imposed accurately without the requirement to choose a proper parameter as in the penalty method. The constraints are satisfied by interpreting the impenetrability between two bodies as a control parameter in a special formula resulted from the Newton’s second law. In some aspects, this algorithm is similar to DENA algorithm [10], but it is more accurate than DENA because the interaction of adjacent contact pairs is considered in the new algorithm. Although the local iterations are needed in this algorithm to calculate contact force, it is quite simple to implement this algorithm in an explicit integration procedure. With this algorithm, no simultaneous equation system need to be solved, thus it is totally consistent with explicit scheme. The efficiency of the explicit integration method does not decrease when the new

Fig. 4. Crashworthiness experiment for a vehicle.

586

F.-J. Wang et al. / Finite Elements in Analysis and Design 43 (2007) 580 – 587

20

Displacement (cm)

A

15 computational

10

experimental 5 0

Fig. 5. The finite-element model of the vehicle frame.

a

0

10

20

30 Time (ms)

40

50

60

Fig. 7. The longitudinal displacement of the frame end.

algorithm is adopted. In contrast, the overall efficiencies increase for the two numerical examples examined in this study as compared to the penalty method. The reason is that the stable time step size increases when the effect of penalty parameter does not exist. The sufficient accuracy can be obtained by this algorithm. This method can be adopted as a new approach to study contact–impact problems with explicit scheme. Currently, research is underway for different kinds of contacts. t=12ms

b

Acknowledgments This work is supported by the National Nature Science Foundation of China (Grant no. 10372114), and by the Program for New Century Excellent Talents in University of China (Grant no. NCET-04-0133). References

t=30ms

c

t=52ms Fig. 6. The deformed configurations of frontal part of the vehicle frame.

[1] T. Belytschko, M.O. Neal, Contact–impact by the pinball algorithm with penalty and Lagrangian methods, Int. J. Numer. Methods Eng. 131 (1991) 547–572. [2] F. Cirak1, M. West, Decomposition contact response (DCR) for explicit finite-element dynamics, Int. J. Numer. Methods Eng. 64 (2005) 1078–1110. [3] J.O. Hallquist, G.L. Goudreau, D.J. Benson, Sliding interfaces with contact–impact in large-scale Lagrangian computations, Int. J. Numer. Methods Eng. 51 (1985) 107–137. [4] P. Wriggers, V.T. Vu, E. Stein, Finite-element formulation of large deformation impact–contact problems with friction, Comput. Struct. 37 (1990) 319–331. [5] T.J.R. Hughes, R.L. Taylor, J.L. Sackman, A. Curnier, W. Kanoknukulchai, A finite element method for a class of contact–impact problems, Comput. Methods Appl. Mech. Eng. 8 (1976) 249–276. [6] N.J. Carpenter, R.L. Taylor, M.G. Katona, Lagrange constraints for transient finite-element surface-contact, Int. J. Numer. Methods Eng. 32 (1991) 103–128. [7] J.G. Malone, N.L. Johnson, A parallel finite-element contact/impact algorithm for nonlinear explicit transient analysis, part 1: the search algorithm and contact mechanics, Int. J. Numer. Methods Eng. 37 (4) (1994) 559–590. [8] N. Hu, A solution method for dynamic contact problems, Comput. Struct. 63 (6) (1997) 1053–1063.

F.-J. Wang et al. / Finite Elements in Analysis and Design 43 (2007) 580 – 587 [9] A.B. Chaudhary, K.J. Bathe, A solution method for static and dynamic analysis of three-dimensional contact problems with friction, Comput. Struct. 24 (1986) 855–873. [10] Z.H. Zhong, L. Nilsson, Lagrange multiplier approach for evaluation of friction in explicit finite-element analysis, Commun. Numer. Methods Eng. 10 (3) (1994) 249–255. [11] D.S. Sha, K.K. Tamma, M.C. Li, Explicit computational developments for elasto-plastic large deformation frictional contact–impact problems, ASME, Comput. Eng. Div. (1994) 95–109. [12] F.J. Wang, Parallel computation of contact–impact problems with FEM and its engineering application, Ph.D. Thesis, Tsinghua University, Beijing, China, 2000. [13] Z.H. Zhong, Finite Element Procedures for Contact–impact Problems, Oxford University Press, Oxford, 1993. [14] K.J. Bathe, Finite Element Procedures in Engineering Analysis, PrenticeHall, Englewood Cliffs, NJ, 1996.

587

[15] D.J. Benson, J.O. Hallquist, A single surface contact algorithm for the post-buckling analysis of shell structures, Comput. Methods Appl. Mech. 78 (1990) 141–150. [16] F.J. Wang, J.G. Cheng, Z.H. Yao, FFS contact searching algorithm for dynamic finite-element analysis, Int. J. Numer. Methods Eng. 52 (7) (2001) 655–672. [17] D.R.J. Owen, E. Hinton, Finite Elements in Plasticity: Theory and Practice, Pineridge Press, Swansea, UK, 1980, pp. 388–389. [18] F.J. Wang, J.G. Cheng, Z.H. Yao, A contact searching algorithm for finite element analysis of contact–impact problems, Acta Mech. Sin. 16 (4) (2000) 374–382. [19] K.J. Bathe, J.W. Walczak, O. Guillermin, P.A. Bouzinov, H.Y. Chen, Advances in crush analysis, Comput. Struct. 72 (1–3) (1999) 31–47. [20] F.J. Wang, Y.T. Feng, D.R.J. Owen, Parallelisation for finite-discrete element analysis on distributed-memory environment, Int. J. Comput. Eng. Sci. 5 (1) (2004) 1–23.