Boundary element analysis of 3-D elasto-plastic contact problems with friction

Computers and Structures 82 (2004) 555–566 www.elsevier.com/locate/compstruc

Boundary element analysis of 3-D elasto-plastic contact problems with friction H. Gun

*

Usak Engineering Faculty, Afyon Kocatepe University, Ordu bulvary, Usak 032 00, Turkey Received 18 July 2003; accepted 1 February 2004

Abstract In this paper, a boundary element method for three-dimensional elasto-plastic stress analysis of frictional contact problems in the framework of small strains and small displacements is presented. The contact conditions are used to unite the different system of equations for each body in contact. Several types of contact interface conditions are covered including infinite friction, frictionless and Coulomb friction. In order to include plastic deformation, the initial strain formulation and von Mises yield criterion are employed. The proposed algorithm is applied to the elasto-plastic analysis of the frictionless rigid punch, a large plate with a rigid, circular, cylindrical inclusion under the increasing load and the ceramic femoral head frictional contact problem. Ó 2004 Elsevier Ltd. All rights reserved. Keywords: Boundary element; Plasticity; 3D-frictional contact

1. Introduction The boundary element (BE) method is well established as an accurate numerical tool, particularly well suited for linear elastic problems. Due to its high resolution of stresses on the surface, the BE approach has been shown to be accurate in problems involving stress concentration, fracture mechanics and contact analysis. However, its extension to non-linear problems including material and geometric non-linearity is not widespread and is under-developed when compared to the finite element (FE). In many non-linear BE formulations the interior of the solution domain has to be discretised, thus losing the main BE method advantage of surfaceonly modeling and resulting in a significant increase in computational effort. Despite this fact, a significant aspect is that the number of simultaneous algebraic equations which have to be solved is only dependent on

*

Fax: +90-276-26341-96. E-mail address: [email protected] (H. Gun).

the boundary discretisation and not on the interior discretisation. Because of the unknown contact area and––in the presence of frictional contact––the unknown contact conditions, elastic contact problems are inherently nonlinear. By including plastic deformation, the numerical analysis gets more complicated. In most engineering applications, severe deformations may take place in the contact area, thus producing rapidly changing stress increments in solution domain. It is evident that the application of the BE approach to contact problems is efficient and suitable. For elastic and thermo-elastic contact problems, the application of the BE method has been widely discussed (see, for example, Refs. [1–7]). To cover the plastic contact analysis, the BE formulation would require a carefully designed robust numerical algorithm in order to incorporate nested iterations and load increments that are capable of monitoring contact developments as well as marching the solution along the elasto-plastic material path. Introducing frictional stickslip behaviour would further complicate the numerical algorithms. In such formulations, two different incremental-iterative procedures are incorporated in order to

0045-7949/$ - see front matter Ó 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruc.2004.02.002

556

H. Gun / Computers and Structures 82 (2004) 555–566

find the proper contact conditions for each load increment. One is to calculate the plastic deformations. The other is to find the correct contact area. Some algorithms using fixed load increments are presented [8,9]. Alcantud et al. [10] presented the extension of the load scaling approach to plasticity, in which linear interpolation is employed to find a permissible load increment. Kong et al. [11] presented a similar mathematical approach in order to treat friction. Huesmann and Kuhn [12] proposed an automatic incremental technique, in which contact conditions change at only one node at the end of the increment, for two-dimensional elasto-plastic contact problems including friction. Their algorithm takes into account the elasto-plastic material behaviour over a fast iterative scheme. For plastic analysis, in general, it is more practical to reduce the number of the load steps to a minimum and to attempt to find average values for the stress and strain increments which are reasonably representative of the particular load step. Although experience is needed in order to specify the optimum size for the load steps, it is possible to roughly determine the size of a particular load step by considering the maximum allowable deviation from proportional loading for the present load step. Such elasto-plastic BE algorithm requires another iterative scheme which does not resort to load incrementation to cover the contact analysis. In this paper, a new extension of the BE formulation to a combined three-dimensional plasticity-contact problems is presented and the method is shown to be an accurate alternative technique to the FE method in handling the plastic deformation under contact loads.

2. BE analytical formulations By considering the plastic strain increments as initial strain increments then modifying Betti’s reciprocal theorem to include plasticity, the pseudo-boundary integral equation for initial strain approach can be written as follows (by neglecting body forces): Z Cki ðP Þu_ i ðP Þ ¼  Tki ðP ; QÞu_ i ðQÞ dSðQÞ S Z þ Uki ðP ; QÞt_i ðQÞ dSðQÞ S Z ð1Þ þ Rkij ðp; qÞ_epi dV ðqÞ V

with the following notations: S is the surface of the body and V is the volume of the solution domain. Uki and Tki are fundamental displacement and traction tensors in the k direction at the field point Q or q due to an orthogonal unit load at the variable point P or p in the i direction. Capital letters are used to indicate that the point concerned lies on the surface S. Also Rkij denotes

the stresses of the corresponding fundamental solution. u_ i , t_i and e_ pi are displacement, traction and plastic strain rates respectively. Cki is the free-term tensor, whose components depend on the geometry. The integral expression for the total strain rates at an interior point p can be obtained by differentiating the corresponding identity for the displacement rate. In the initial strain approach, the convected differentiation of the related domain integrals is employed and at internal points the correct expression for the total plastic strains can then be written as follows: Z e e_ kh ðpÞ ¼  Sikh ðp; QÞu_ i ðQÞ dSðQÞ S Z þ Deikh ðp; QÞt_i ðp; QÞ dSðQÞ S Z Reijkh ðp; qÞ_epij ðp; qÞ dV ðqÞ  V

8  10m p e_ ðpÞ þ 15ð1  mÞ kh

ð2Þ

e in which Dekij , Skij and Reijkh are the derivatives of the aforementioned fundamental solutions [13,14] and Reijkh is given by:

Reijkh ¼

1 1 3ð1  2mÞdij r;k r;h þ 3dkh r;i r;j 8pð1  mÞ r3  15r;i r;j r;k r;h þ 3tðdjh r;i r;k þ djk r;i r;h þ dih r;j r;k þ dik r;j r;h þ ð1  2mÞðdjk dih þ dik djh   dij dkh Þ

ð3Þ

For a material obeying the von Mises yield criterion the plastic strain increments are given by the following incremental elasto-plastic flow rule: ! S_ kl e_ kl S_ ij 3 _epij ¼ ð4Þ 0 2 1 þ H =3 l ðr_ eq Þ2 in which S_ ij and r_ eq denote the current deviatoric stress tensor and the equivalent stress respectively. e_ pij is the  is the shear modulus, and H 0 total strain increments, l represents the plastic modules. The formulations can be also extended to cover kinematic, isotropic or mixed hardening material behaviour [8,13].

3. Numerical implementation The detailed formulation of the BE method for the plastic analysis using isoparametric quadratic elements is well covered in the literature (see, for example, Gao and Davies [14,15]) and will be summarized here. Since interior modeling is required, both boundary surface elements and volume cells are necessary to perform the integrals. For the boundary surface element, the geo-

H. Gun / Computers and Structures 82 (2004) 555–566

metry can be described in terms of quadratic shape functions in a local co-ordinate axes system as follows: xi ðn1 ; n2 Þ ¼

8 X

Nc ðn1 ; n2 Þðxi Þc

ð5Þ

c¼1

where Nc is the quadratic shape function and n is the local co-ordinate. Similarly, the displacement and traction vectors can be expressed in terms of quadratic shape functions. For the interior cells, the geometry can be defined in terms of quadratic shape functions, which are described in local co-ordinates, n1 , n2 and n3 as follows: xi ðn1 ; n2 ; n3 Þ ¼

20 X

Nc ðn1 ; n2 ; n3 Þðxi Þc

ð6Þ

c¼1

Similarly, the displacement increments can be expressed in the domain cells. Therefore, in discretised form, for each body in contact, the elasto-plastic boundary integral equation in the initial strain approach (by neglecting body forces) can be written as follows: M X 8 X

Cki ðP Þu_ i ðP Þ þ 

Z

þ1

Z

1

¼

m¼1 c¼1 þ1 1



1

e_ pij ðqÞ

m¼1 c¼1 þ1 X

Z

Z

þ1

1

Z

Z

M X 8 X

þ1

1 M X

1

e Sikh ðp; QÞNc ðn1 ; n2 ÞJ ðn1 ; n2 Þ dn1 dn2

8 X

t_i ðQÞ m¼1 c¼1 Z þ1 Z þ1

þ

1

ðP ; qÞNc ðn1 ; n2 ; n3 ÞJ ðn1 ; n2 ; n3 Þ dn1 dn2 dn3



ð7Þ

þ

kij

where P denotes the node where the integration is performed. Q indicates the cth node of the mth boundary surface element and q indicates the cth node of the mth volume cell. It should be noted that the integration process is performed separately for each domain in contact. A set of linear algebraic equations is obtained for each contacting body after performing the necessary numerical integration over the boundary elements and domain cells, which can be formed as follows: h i _ ¼ ½B ½t_ þ ½R e_ p ½A ½u ð8Þ in which the matrices ½A , ½B and ½R contain the integrals of the displacement, traction and domain kernel integrals, respectively. For three-dimensional problems, if the total number of boundary nodes is N and the total number of the interior cell points is H , then the solution matrices ½A and ½B will be square matrices of size

u_ i ðQÞ

m¼1 c¼1 þ1 Z þ1



Uki ðP ; QÞNc ðn1 ; n2 ÞJ ðn1 ; n2 Þ dn1 dn2

1

D X 20 X

1

¼

t_i ðQÞ

m¼1 c¼1 Z þ1 Z þ1





Tki ðP ; QÞNc ðn1 ; n2 ÞJ ðn1 ; n2 Þ dn1 dn2

M X 8 X

þ

3N  3N, whereas the matrix ½R will be a rectangular matrix of size 3N  6H . All matrices are fully populated. Since the solution procedure consists of solving Eq. (8) in iterative manner, the first step is to assume zero plastic strain rates and to obtain the first approximate solution for the primary unknowns, which are traction and displacement rates. The second set of approximate solutions for e_ pij , acting as the secondary unknowns, is next found by applying Eq. (6) and the procedure is repeated until convergence is achieved. The total strain rates at all the nodes have to be computed before the flow rule can be employed. As far as the boundary nodes are concerned, these rates are readily available from the local traction and displacement gradient rates. For the interior nodes, the displacement rates are first found from an integral equation similar to Eq. (7). One way of obtaining the stress and total strain rates is to utilize the integral identities for the strain rates and to subsequently apply the stress–strain relationship to yield the stress rates. The discretised form of Eq. (2) is e_ kh ðpÞ þ

u_ i ðQÞ

557

1 D X

1 20 X

Deikh ðp; QÞNc ðn1 ; n2 ÞJ ðn1 ; n2 Þ dn1 dn2

e_ pij ðqÞ

Z

þ1

Z

þ1

1 1 m¼1 c¼1 þ1 Reijkh ðp; qÞNc ðn1 ; n2 ; n3 ÞJ ðn1 ; n2 ; n3 Þ dn1 dn2 dn3 1

Z

8  10m p e_ ðpÞ 15ð1  mÞ kh

ð9Þ

The evaluation of the equation coefficients pertaining to the integrals over the boundary and in the solution domain is well covered in the literature [14].

4. Frictional contact conditions When two three-dimensional bodies A and B come into contact, each body has an associated normal direction x3 and two tangential direction x1 and x3 which are orthogonal to each other and follow the right-hand side screw rule. Fig. 1 shows the local axes directions at a single point for two arbitrary bodies in contact. There are six contact variables for each contacting node, namely displacement and traction rates in three directions. To unite the system of equations obtained for each body in contact, the contact condition must be imposed

558

H. Gun / Computers and Structures 82 (2004) 555–566

sticking or slipping occurs in the local direction of node (a). The integer variables Ex1 and Ex2 take value of +1 or )1 if slip occurs in the same or opposite direction as the local x1 or x2 directions of node (a). The angle b is measured between the local x1 direction and the direction of the resultant slip, t_r .

5. Contact iterations

Fig. 1. 3D contact.

on the contacting node pairs. Let node (a) from body A and (b) from body B make a node pair; then if slip occurs under coulomb friction conditions, by taking into account the local co-ordinate orientation, see Fig. 1, constraint equations based on slip in a resultant tangential direction can be expressed as follows: ðaÞ

ðbÞ

u_ x3 ¼ u_ x3 þ d ðaÞ ðbÞ u_ x1 ¼ u_ x1 ðaÞ

ðbÞ

u_ x2 ¼ u_ x2 ðaÞ ðbÞ t_x3 ¼ t_x3 ðaÞ ðbÞ t_x1 ¼ t_x1 ðbÞ a ¼ t_x2 t_x2

ð10Þ

t_ra ¼ Klt_nðaÞ ðaÞ t_x1 ¼ Ex1 t_rðaÞ cos b ðaÞ t_x2 ¼ Ex2 t_rðaÞ sin b qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðaÞ ðaÞ t_rðaÞ ¼ t_x1 þ t_x2 "

b ¼ tan1 abs

In contact problems, each domain is treated separately to form Eq. (7), and the resulting matrices [A] and [B] are coupled together according to the relevant contact conditions, with the number of unknowns remaining equal to the number of equations. The compatibility and equilibrium equations are exactly modeled by coupling the relevant variables directly. The implementation of the contact constraints in the system of equations is relatively more efficient and accurate than the FE method because both the tractions and displacements exist as nodal variables whereas, the FE equations involve only the displacements as the nodal degrees of freedom. In most practical contact problems, the extent of the contact area and the stick-slip conditions are initially unknown and must be determined by an iterative procedure and, if friction is present, the behaviour may be dependent on load history. An efficient automatic iterative scheme [4–6] can be employed in which interpenetration of the contacting elements or overlap just outside the contact area is prevented, and elements with tensile stress are released from the contact area. The iteration procedure also allows the nodes with the shear stress exceeding the coulomb friction limit to slip in the next iteration with a limit imposed on the value of the tangential traction. The most difficult task in three-dimensional contact algorithm is to determine the direction of slip which occurs along contacting surface. During the first contact iteration the bodies are assumed to have no slip between contacting surfaces. In the second iteration, the direction of the slip is set opposite to the direction of the first iteration over relative slip acting in a single resultant tangential direction between node pairs i.e. frictional traction and displacement oppose each other [6].

6. Coupling the system equations ðaÞ t_x2 ðaÞ t_x1

!#

where d is the initial gap between the contact nodes and l is the coefficient of friction. The integer variable K is an on–off switch which is set to 0 or 1 respectively if

Multi-domain formulation requires the treatment of each domain separately and coupling sets of equations for the nodes at the contact interface, prior to the solution of all of the equations. With, say, two bodies in contact, the two set of simultaneous equation can be expressed in the matrix form as follows:

H. Gun / Computers and Structures 82 (2004) 555–566

½A ðaÞ 0

0 ½A ðbÞ



_ ðaÞ ½u _ ðbÞ ½u



¼

½B ðaÞ 0

þ

 ½t_ ðaÞ ½t_ ðbÞ 2 h iðaÞ 3  p 0 6 e_ 7 ðbÞ 4 h iðbÞ 5 ½R p e_

0 ½B ðbÞ

½R ðaÞ 0



ð11Þ The matrices ½A , ½B matrices as follows: 2 ½A 11 ½A 12 ½A 13 6 ½A 21 ½A 22 ½A 23 6 6 ½A 31 ½A 32 ½A 33 4 .. .. .. . . .

and ½R can be expressed as sub-

3 32 _1 ½u ... 6 ½u _ 27 ...7 7 76 6 ...7 ½u _ 37 7 56 4 . 5 . ... . 3 2 32 ½B 11 ½B 12 ½B 13 . . . ½t_ 1 6 ½B 21 ½B 22 ½B 23 . . . 76 ½t_ 7 7 6 76 6 _ 27 ¼ 6 ½B 7 ½B ½B . . . ½t 3 7 32 33 4 31 56 4 . 5 .. .. .. .. . . . 2h i 3 2 3 e_ p ½R 11 ½R 12 ½R 13 . . . 6 h i1 7 6 ½R 21 ½R 22 ½R 23 . . . 76 e_ p 7 7 6 76 6 h i2 7 þ 6 ½R 7 ½R ½R . . . 7 32 33 p 4 31 56 6 e_ 7 .. .. .. 4 35 . . . ... .. .

ð12Þ

Where the submatrices of ½A , ½B and ½R are defined as follows: 2 3 2 3 Axx Axy Axz Bxx Bxy Bxz 6 7 6 7 ½A ij ¼ 4 Ayx Ayy Ayz 5; ½B ij ¼ 4 Byx Byy Byz 5;

½R kij

Azx Azy Azz 2 Rxxx Rxyy Rxzz 6 ¼ 4 Ryxx Ryyy Ryzz Rzxx

Rzyy

Rzzz

Bzx

Bzy Bzz 3

Rxxy

Rxyz

Rxxz

Ryxy

Ryyz

7 Ryxz 5

Rzxy

Rzyz

Rzxz

ð13Þ

_ ½t_ and ½e_ p represent the displacement, The vectors ½u , traction and plastic strain rates in the Cartesian axes as follows: 2 3 2 3 u_ x t_x 6 7 6_ 7 _ _ i ¼ 4 u_ y 5 ; ½t i ¼ 4 ty 5 ; ½u u_ z i t_z i h i h iT ð14Þ e_ p ¼ e_ px e_ py e_ pz e_ pxy e_ pyz e_ pxz i

559

displacement and three tractions) for each node pair within the contact region. Since the contact conditions for node pair are unknown, substitution for six of the unknowns in Eq. (11), with contact variables from Eq. (10) results in a set of six equations with six unknowns. To arrive the solution matrix, coefficients of the related variables are combined within the matrix and the equations are rearranged such that the unknown variables are on the left-hand side of the equation and the known variables are on the right-hand side of the equation. The form of the linear equations for the combined bodies can be expressed as follows: h i _ þ ½R e_ p ½A ½_x ¼ ½C ð15Þ in which only unknowns are in vector ½_x . Considering the contact conditions in Eq. (10) and assuming elastic deformation, Eq. (15) becomes 2 3 2 3ðaÞ A11 A12 A13 6 7 4 A21 A22 A23 5 6 7 6 7 A A A 31 32 33 6 7 2 3ðaÞ 6 7 6 7 C11 C12 C13 6 7 4 C21 C22 C23 5 6 7 6 7 6 7 C31 C32 C33 62 7 3 6 ð1  KÞA11 ðK  1ÞA12 A13 ðbÞ 7 6 7 6 4 ð1  KÞA21 ðK  1ÞA22 A23 5 7 6 7 6 ð1  KÞA31 ðK  1ÞA32 A33 7 6 7 2 3ðbÞ 6 7 6 7 D11 D12 D13 6 7 4 D21 D22 D23 5 4 5 D31 D32 D33 2 3 2 3 ðaÞ _ u x1 6 ðaÞ 7 6 7 6 u_ 7 6 7 4 x2 5 6 7 ðaÞ 6 7 _ u x3 62 37 7 6 6 7 ðaÞ ðbÞ 6 6 ð1  KÞt_x1 þ K u_ x2 7 7 66 7 ðaÞ ðbÞ 6 4 ð1  KÞt_x2 þ K u_ x2 7 57 4 5 ðaÞ t_x3 2

3 0 6 7 0 6 7 6 7 0 ðbÞ 7 ¼6 6 d½A13 7 6 7 4 d½A23 ðbÞ 5 d½A33 ðbÞ

ð16Þ

in which

i

For contact equations only, it is necessary to substitute the global contact variables of Eq. (11) with local contact variables. Hence, it is necessary to convert the submatrices ½A , ½B from global expressions to local expression using the three-dimensional direction cosines. While there are only three unknowns for each node outside the contact region, there are six unknowns (three

C11 ¼ ðK  1ÞB11;

C12 ¼ ðK  1ÞB12

C13 ¼ B13  KEx1 l cos bB11  KEx2 l sin bB12 C21 ¼ ðK  1ÞB21;

C22 ¼ ðK  1ÞB22

C23 ¼ B23  KEx1 l cos bB21  KEx2 l sin bB22 C31 ¼ ðK  1ÞB31;

C32 ¼ ðK  1ÞB32

C33 ¼ B33  KEx1 l cos bB31  KEx2 l sin bB32

ð17Þ

560

H. Gun / Computers and Structures 82 (2004) 555–566

D11 D12 D13 D21 D22

¼ ð1  KÞB11 þ KA11 ¼ ðK  1ÞB12 þ KA12 ¼ B13 þ KEx1 l cos bB11  KEx2 l sin bB12 ¼ ð1  KÞB21 þ KA21 ¼ ðK  1ÞB22 þ KA22

D23 D31 D32 D33

¼ B23 þ KEx1 l cos bB21  KEx2 l sin bB22 ¼ ð1  KÞB31 þ KA31 ¼ ðK  1ÞB32 þ KA32 ¼ B33 þ KEx1 l cos bB31  KEx2 l sin bB32

ð18Þ

For a given set of approximation to plastic strain rates, in the presence of plastic deformation, the contents of _ ½R are multiplied with e_ p and the results are added to ½C and the final form of the equations can be expressed as ½A ½_x ¼ ½C_

ð19Þ

where the vector ½C_ is an updated known vector and the matrix ½A is the final coefficient matrix. Because of the nature of the both boundary integral identities and contact conditions, the equations need to be solved in an iterative manner. Therefore, two different iterative procedures are incorporated in order to find the proper contact area and the plastic deformations [16]. Within the first load increment, the solution procedure consists of solving the system of equations in iterative manner by initially assuming zero plastic strain rates and obtaining the first approximate solution for ½_x . The second set of approximations is found next by applying Eq. (4) and the procedure repeats until convergence is achieved [16]. For the next load increment, the initial strains can be assumed to be the strain rates of the previous load step. The resulting simultaneous equations can be solved using Gaussian elimination with partial pivoting. The stresses at the boundaries are obtained from the local displacement gradient and traction rates.

of the model is analysed with zero displacements prescribed in the x direction along the axis of symmetry and the y direction along the base of the foundation. The punch is compressed by means of a uniform arbitrary displacement of 0.18 in the y directions. BE mesh design having 1322 boundary surface elements and 3908 nodes are employed, as shown in Fig. 2. The regions where plastic deformation is expected are discretised using 512 20-noded volume cells. The foundation is assumed to have the following material properties: Young’s modulus E ¼ 1000, yield stress r ¼ 1 and Poisson’s ratio m ¼ 0:3. The elastic modules of the rigid punch is assumed to be 105 times greater than that of the foundation with the same Poisson’s ratio. For comparison purposes the contact is assumed to be frictionless. This problem was the subject of a benchmark finite element study carried out on behalf of NAFEM and reported by linkens [17]. Calculations are carried out, assuming perfect plasticity. The stresses in the x and y directions at the specific point A, which is located just below the edge of the punch, at the co-ordinate ð80; 140Þ, are plotted in Fig. 3. The strains in the x and y directions at the specific point A are also plotted in Fig. 4. It can be seen that there is a very good agreement between the present BE results and FEM (reference), which is believed to have emerged from a convergence study with much finer mesh discretisations. 7.2. Circular inclusion This BE application, which finds many engineering applications in connections and in joints, represents the plane strain analysis of a circular inclusion contact problem between a rigid shaft and an infinite body with cavity in which the shaft is assumed to be a perfect fit in the hole before loading.

7. Numerical applications In order to asses the accuracy of the present BE algorithm, the frictionless rigid punch, a large plate with a rigid, circular, cylindrical inclusion under the increasing load and the ceramic femoral head contact problem are studied. 7.1. Rigid punch contact problem A punch of height 100, half width 80 and thickness 10 is pressed against a foundation of height 160, half width 200 and thickness 10. Dimensions of this problem are given in consistent units such as N and mm. Plane strain is simulated in the three-dimensional model by constraining the relevant faces of the punch and the foundation in the z direction. A symmetric half

Fig. 2. BE mesh design for the rigid punch.

H. Gun / Computers and Structures 82 (2004) 555–566

Direct stresses at point A

0 -0.2

-0.4

0.1

561

σm

0.2 σx

FEM (Reference ) BEM

Y

-0.6

α -0.8

σy

R

X

-1

Applied deflection of punch

Fig. 3. Stresses at the specific point A.

εx

εy

Fig. 4. Strains at the specific point A.

For the modeling purposes, a large plate of the half height H and of the half width W is used instead of an infinite matrix. A remote tensile stress rm is applied to the matrix causing the inclusion to separate from the shaft. The contact interface is assumed to be frictionless. The stresses are measured at the angles around the edges of the hole, as defined in Fig. 5. The numerical values of the modeling of the geometry are R ¼ 25 mm, H =R ¼ W =R ¼ 8 and thickness t ¼ 10 mm. Plane strain is simulated in the three-dimensional model by constraining the relevant faces of the plate and the shaft in the z direction. A half of the model is analysed with zero displacements prescribed in the x and y directions along the axis of symmetry. BE mesh design having 1140 surface elements with 3101 nodes are employed. For the plastic analysis, the plate is discretised in two subregions, as shown in Fig. 6. Only the first region, covering the root of the plate has 56 20-noded interior cells. It should be worth mentioning that there are two major differences between element

Fig. 5. Shaft in a hole in an infinite plate.

modeling strategies of the BE and FE method, which are retained in the present BE formulation. The first difference is that the BE mesh for 20-noded interior elements (volume cells) is designed to cover only the region where yielding is expected to occur. This makes the BE method more attractive for the problems displaying limited plastic deformation. Secondly the domain and boundary meshes in BE method are independent. This allows more flexibility in BE mesh design. To assess the present BE contact algorithm, only elastic deformations under the frictionless contact conditions are considered. For comparison purposes, the shaft and the plate are assumed to have the same material properties E ¼ 73 GPa and m ¼ 0:3. It should be noted that only the boundary discretisation of the bodies are needed for the elastic deformations. The radial contact stresses, rr , at the edge of the hole against the angle, a, defined in Fig. 5, are indicated in Fig. 7. The circumferential stresses, rh , also against the angle, a, are plotted in Fig. 8. Both the present BE results and the analytical results [18] are in a very good agreement. The elasto-plastic analysis of this problem was investigated by using FE formulation [19] and 2D BE formulation [8]. The shaft is assumed to be rigid with the material properties of the plate E ¼ 73 GPa, ry ¼ 386 MPa, m ¼ 0:3 and plastic modulus, H 0 ¼ 620 MPa. The plate starts to become plastic at the load of 2.98 MN/m. Therefore three values of 3, 3.4 and 4.4 MN/m, remotely applied uniform tensile load are employed with 1, 5 and 7 load increments respectively. From comparison purposes, two more values of uniform tensile load, 0.3504 and 1.41912 MN/m, are also applied. Fig. 9 shows the variation of the circumferential stress, rh , versus the angle, a, for the three values of the remotely applied

562

H. Gun / Computers and Structures 82 (2004) 555–566

Fig. 6. BE mesh design for the inclusion problem.

450

0.7

Circumferential Stresses (MPa)

3 MN/m

0.6

σr / σm

0.5 0.4 0.3

Analitical solution

0.2

BEM

0.1

4.4 MN/m 2D-BEM [8] 1.41912 MN/m

250

0.3504 MN/m FEM [19]

150

50 -50 0

20

40

60

80

-150

0 0

5

10

15

20

25

Circumferential Angle, α Fig. 7. Radial stress distribution over the contact region.

3.5 3 Analytical solution 2.5

Alpha, α ,in Degree

Fig. 9. Circumferential stresses on elasto-plastic plate around rigid shaft.

uniform tensile load. The radial stresses against the angle, rr , are plotted in Fig. 10. It is clear that the present BE results at the contact arc of the plate are satisfactory when compared with FE results [19] and 2D BE results [8].

BEM

2

σθ /σ m

3.4 MN/m

350

7.3. Compression of a ceramic femoral head

1.5 1 0.5 0 20

40

60

80

-0.5 -1

circumferential angle, α

Fig. 8. Circumferential stress distribution at the edge of the hole.

This example is chosen to demonstrate the applicability of three-dimensional structure with two contact regions involving three dissimilar materials. The details of the axisymmetric profile and co-ordinates of the titanium spigot, femoral head and brass ring are given in Fig. 11. The outside of radius of the femoral head is specified as 14 mm, which is a standard size for many hip joint replacements. Dimensions of this problem are given in consistent units such as N and mm. The titanium spigot is assumed to have the material properties E ¼ 1:14  105 and m ¼ 0:33 and yield stress, r ¼ 800. The femoral head is made of an alumina

H. Gun / Computers and Structures 82 (2004) 555–566 3 MN/m 3.4 MN/m 4.4 MN/m 2D-Bem[8] 1.41912 M N/m 0.3504 MN/m FEM [19]

140

Radial stresses (MPa)

120 100 80

K6

KEYPOINTS : Z, R (mm) K1 : 6.8718, -10.25 K2 : 9.5361,-10.25 K3 : 9.798 ,10 K4 : 16 ,10 K5 : 16 ,17.5 K6 : 7 ,17.5 K7 : 7 ,12.1244 K8 : 0 ,14 K9 : 6 ,9.5 K10 : 6.5 ,-3 K11 : 0 ,-3 K12 : 0 ,6.75

40 20 0 5

10

15

20

25

K4 R2 K3

K12

R1

K11 K10

RADIUS MAGNINUTE (mm) R1 R2

: 0.75 : 1.2 K1

Fig. 10. Radial stresses on contact surface.

ceramic with E ¼ 2:0  105 , m ¼ 0:3 and yield stress, r ¼ 200 and the brass ring has the material properties E ¼ 1:05  105 , m ¼ 0:35 and yield stress, r ¼ 420. A uniform displacement of 0.1 mm in y direction is applied to the brass ring, compressing the femoral head. The spigot is restrained in all directions between key points K10 and K11 and in the z direction between key points K11 and K12, which is axis of symmetry. The femoral head is restrained in the r direction between key points K8 and K9. It is clear from Fig. 11 that there are two distinct contact regions between the three parts. The contact region between the inside of the femoral head and the titanium spigot is assumed to be glued, while the contact region between the outside of the femoral head and the ring is assumed to have slip-stick contact conditions. Therefore three different values of coefficients of friction, l ¼ 0:0, l ¼ 0:1 and l ¼ 0:2 are employed in the analysis to monitor the effect of friction on the von Mises stresses. To depict how von Mises stresses vary around the inside and the outside of the femoral head, the inside angle, al and a0 are defined with respect to datum ð0; 0Þ, as shown in Fig. 11. The two BE models shown in Fig. 12 are employed in the present analysis. In the first part of the analysis, only elastic deformations are considered. The von Mises stresses around the inside of the femoral head are shown in Fig. 13, revealing that the stresses peak in the two distinct locations at the edges of contact between the spigot and the femoral head and that they remain uniform in the rest of the contact region. The elastic analysis of this problem was investigated by using 3D and axisymmetric BE formulations [6,20]. The von Mises stresses around the outside of the femoral head are shown in Fig. 14. The peak stresses appear at the edges of the contact between ring and the head. The effect of increasing the coefficient of the friction is to increase the contact stresses.

K5

K7

K8 K9

60

0

563

K2

α0 αi

Z

( 0,0 )

r

Fig. 11. Geometric parameters of the hip joint contact problem.

In the next stage of the analysis elasto-plastic deformations are considered under the elastic–perfectly plastic material behaviour. The BE mesh design has 1062 more 20-noded volume cells and the finer BE mesh design has 1014 20-noded volume cells for the ceramic femoral head since plastic deformations occur only in the femoral head. It is clear that the two BE mesh design results are in a very good agreement. The von Mises stresses around the inside and the outside of the femoral head are shown in Figs. 15 and 16 respectively. The equivalent plastic strains around the inside of the femoral head are shown in Fig. 17. The peak equivalent plastic strains appear at the edges of contact between the spigot and the femoral head and they remain uniform in the rest of the contact region. The equivalent plastic

564

H. Gun / Computers and Structures 82 (2004) 555–566

2500

µ=0.1 µ=0.2 BEM mesh-2

2

Von Misses ( N/mm )

µ=0.0 2000

1500

1000

500

0 0

20

40

60

80

100

120

140

Alpha, α i ,in degree

Fig. 13. Elastic stress distribution around the inside of the femoral head.

350

µ=0.0 µ=0.1 µ=0.2

2

Von Misses ( N / mm )

300 250

BEM mesh-2 200 150 100 50 0 0

20

40

60

80

Angle,α °, in degree

Fig. 14. Elastic stress distribution around the outside of the femoral head.

350

2

Von Mises (N/mm )

300 250 200

µ=0.0

150

µ=0.1

100

µ=0.2 BEM mesh-2

50

Fig. 12. Two BE mesh designs for femoral head assembly. 0 0

20

40

60

80

100

120

140

Alpha,α i , in degree

strains around the outside of femoral head are shown in Fig. 18. The shear stress distribution around the inside

Fig. 15. Elasto-plastic stress distribution around the inside of the femoral head.

H. Gun / Computers and Structures 82 (2004) 555–566

14

2

Shear stresses ( N/mm )

BEM mesh-2

2

Von Mises (N/mm )

16

µ= 0.0 µ= 0.1 µ= 0.2

200

150

100

50

12 10

µ=0. 2

8

µ=0. 1

6

µ=0. 0

4 2 0

0 0

20

40

60

45

80

65

85

105

Angle,α i, indegree

Angle,α o , in degree

Fig. 16. Elasto-plastic stress distribution around the outside of the femoral head.

Fig. 19. Shear stress distribution around the inside contact interface of the femoral head.

0.2

1

0.18

µ= 0.0 µ= 0.1

0.14

µ= 0.2

0.5 2

0.16

Shear stresses ( N/mm )

Equivalent Plastic Strain

565

0.12 0.1 0.08 0.06 0.04 0.02

0 30

35

40

-0.5

µ=0.2 -1

µ=0.1 µ=0.0

-1.5

BEM mesh-2

0 0

20

40

60

80

100

120

140

Alpha,αi , in degree

Fig. 17. Plastic strain distribution around the inside of the femoral head.

-2

Angle,α o, in degree

Fig. 20. Shear stress distribution around the outside contact interface of the femoral head.

and the outside contact interface of the femoral head are also shown in Figs. 19 and 20 respectively.

Equvialent Plastic Strain

0.014

µ= 0.0

8. Conclusion

µ= 0.1

0.012

µ= 0.2

0.01 0.008 0.006 0.004 0.002 0 25

30

35

40

45

50

Angle, α o ,indegree

Fig. 18. Plastic strain distribution around the outside of the femoral head.

A quadratic boundary element formulation for threedimensional elasto-plastic stress analysis of frictional contact problems in the framework of small strains and small displacements is presented. The contact conditions are expressed as algebraic relationships between the displacements and tractions of the bodies in contact and directly incorporated in the equation solver. This direct coupling of the contact variables means that it is not necessary to use any form of gap or interface elements between the contacting surfaces. The purpose of the volume cells is to evaluate the interior volume integrals and, unlike the FE method, volume cells do not have to be interior-connected with the boundary elements. The

566

H. Gun / Computers and Structures 82 (2004) 555–566

proposed algorithm is applied to the elasto-plastic analysis of the frictionless rigid punch, a large plate with a rigid, circular, cylindrical inclusion under the increasing load and the ceramic femoral head frictional contact problem. The results are compared with the published FE and BE results, where possible, and show an excellent agreement. The method is convenient for handling the plastic deformation under the contact loads with friction.

[10]

[11]

References [1] Anderson T. Boundary element in two-dimensional contact and friction. Linkoping: Linkoping University Report Dept. of Mech. Eng., 1982. [2] Karami G. Boundary element method for two-dimensional contact problems. Berlin: Springer; 1989. [3] Karami G, Fenner RT. A two-dimensional BEM for thermo-elastic body forces contact problems. In: Brebbia CA, Wendland WL, Kuhn G, editors. Boundary element IX stress analysis applications. Berlin: Springer; 1987. [4] Olukoko OA, Becker AA, Fenner RT. A new boundary element approach for contact problems with friction. Int J Numer Methods Eng 1993;33:1513–22. [5] Olukoko OA, Becker AA, Fenner RT. A review of three alternative approaches to modeling frictional contact problems using the boundary element method. Proc Roy Soc London 1994;444(A):37–51. [6] Leahy JG, Becker AA. A quadratic boundary element formulation for three-dimensional contact problems with friction. J Strain Anal 1999;34:235–51. [7] Eck C, Steinbach O, Wenland WL. A symmetric boundary element formulation for contact problems with friction. Math Comput Simul 1999;50(1–4):43–61. [8] Karami G. Boundary element analysis of two-dimensional elastoplastic problems. Int J Numer Methods Eng 1993;36:221–35. [9] Tsuta T, Yamaji S. Boundary element analysis of contact thermo-plastic problems with creep and the numerical

[12]

[13]

[14]

[15] [16]

[17]

[18]

[19]

[20]

technique. In: Boundary Elements Proceedings of 5th International Conference. Berlin: Springer; 1993. Alcantud M, Doblare M, Gracia L, Dominguez J. Analysis of the contact problems between elastoplastic 2D bodies by means of the BEM. In: Brebbia CA, Chaudonet-Miranda A, editors. Boundary element in mechanical and electrical engineering. Shouthampton: Computational Mechanics Publicationm; 1992. Kong X-A, Gakwaya A, Cardou A, Cloutier L. A unified BEM formulations for elasto-plastic contact analysis by mathematical programming. In: Brebbia CA, Domingues J, Paris F, editors. Boundary elements XIV stress analysis and computational aspects. Southampton: Computational Mechanics Publications; 1992. Huesmann A, Kuhn G. Automatic load incrementation technique for plane elasto- plastic frictional contact problems using boundary element method. Comput Struct 1995;56:733–45. Lee KH, A boundary integral equation method for two dimensional elastoplastic analysis. Phd thesis, University of London, 1983. Gao X-W, Davies TG. An effective boundary element algorithm for 2D and 3D elastoplastic problems. Int J Solids Struct 2000;37:4987–5008. Gao X-W, Davies TG. Boundary element programming in mechanics. Cambridge University Press; 2002. Gun H. An effective BE algorithm for elastoplastic frictional contact problems. Eng Anal Boundary Elem, in press. (doi: 10.1016/j.enganabound.2003.09.003). Linkens D. Selected Benchmarks for Material Non-linearity. National Association for Finite Element Method and Standards Report R0026, East Kilbride, UK, 1993. Stippes M, Wilson HB, Krull FH. Contact stress problem for a smooth disc in an infinite plate. In: Proceeding of the 14th US National Congress on Mechanics, 1962. p. 799– 806. Vojiadjis GZ, Navaee S. Finite strain contact problems of cylinder embedded in body. J Eng Mech ASCE 1984;110:1597–609. Gun H. Hip joint applications using the boundary element method. Boundary Elem Commun 2001;12(2/3):10–7.