Comput Mech (2007) 41:159–174 DOI 10.1007/s00466-007-0175-9

ORIGINAL PAPER

Numerical solution of elastic bodies in contact by FEM utilising equilibrium displacement fields ˇ Lubor Fraštia

Received: 27 July 2006 / Accepted: 13 March 2007 / Published online: 30 May 2007 © Springer-Verlag 2007

Abstract In this paper a new effective formulation of the computation of the statics of elastic bodies in contact is described and demonstrated. Line contact, plane strain, and negligible friction are assumed. The formulation is an extension of the standard finite element method (FEM). With the aim to utilise the Hertz theory directly, we use the exact solution of the elastic 2-dimensional halfspace loaded by Hertzian pressure distribution and enforce the contact condition by the method of Lagrange multipliers. In numerical examples we have focussed on the demonstration and evaluation of the accuracy of the new formulation for selected applications compared with the state of the art node-to-segment contact algorithm implemented in the software system ADINA. Proposed formulation is more accurate for problems where Hertz contact dominates the strain state, especially for small number of elements, whereas we obtained a fairly good agreement with ADINA for a more general bending problem. Keywords Contact · Elasticity · Finite element · Static analysis · Variational method

1 Introduction The contact of two elastic bodies was theoretically explained by Hertz already in 1882 [1]. Results of his theory hold assuming that, see also [2, Chap. 4]:

ˇ Fraštia (B) L. Department of Applied Mechanics, SjF, ŽU, Univerzitná 1, 01026 Zilina, Slovakia e-mail: [email protected]

• the contact surfaces are continuous and nonconforming, i.e. for line contact it means a
1/K , where a is the halfwidth of the contact surface and K is the relative curvature K = K 1 + K 2 , where K i is the curvature of ith body at the centre of the contact surface. • strains are small, linear elasticity theory is sufficient; • each body can be considered as a halfspace; • friction is neglected. Hertz theory was used in the analysis of stresses and displacements in the vicinity of the contact surface [3,4]. Results of that research were used in practical applications such as approximate determination of stresses and displacements of the elastic cylinder between two bodies [5]. In general, with this approach, one tries to simplify geometry of the real problem and then to find the solution in closed form. These days the most frequently used numerical method for stress analyses is the finite element method (FEM). Engineers use FEM in cases where complicated geometry, material properties etc. causes that obtaining analytic solution of stresses is not effective and often impossible. The problem of the stress analysis of bodies in contact needs special care and one has to design proper FEM formulation for it. It is important to design a robust algorithm which incorporates contact conditions. During last 35 years several such algorithms were proposed and used. For a review, the reader may consult [6,7]. Usually the classification of numerical approaches for contact problems is on the basis of the mathematical formulation of the problem as follows: 1. the formulation is in the form of variational inequalities or inequality principle of virtual work, see e.g. [8], and usually optimization numerical methods are employed.

123

160

2. the formulation is in the form of variational equalities, or using principle of virtual work, e.g. [9]. Presented formulation is based on this approach. There are several possibilities how to enforce the contact condition, for example: penalty function method, method of Lagrange multipliers and its generalizations such as augmented Lagrangian method or perturbed Lagrangian method, barrier methods, constraint function method, mortar method. Many solution procedures, employing these techniques, were described and studied in the literature, see e.g. [9–13]. Up to now, several general purpose algorithms have been successfully implemented in commercial FEA packages. Most of the above cited works are general purpose approaches for contact with friction. A general purpose numerical algorithm has to deal with the difficulties that are inherent in the contact as a strongly nonlinear problem: there are often problems with convergence, it is important to use a very fine mesh. The aim of this paper is to propose an algorithm for special cases where the Hertz theory is a very good approximation. We want to exploit this theory directly in the formulation and to obtain better numerical scheme. Similar motivation was behind other earlier contributions in this domain, for example [14], where we used hybrid FEM formulation. In the Hybrid FEM [15], the equilibrium field, so called Trefftz field, is defined on the finite element. Trefftz fields are functions satisfying all governing equations in the strong sense, i.e. equilibrium is satisfied at each point. These Trefftz fields have in practice finite number of independent parameters and so they cannot satisfy arbitrary boundary condition in the strong sense. As a result, boundary conditions are satisfied in weak sense. For this purpose, a new independent displacement field is defined on the boundary of the element and compatibility between elements is enforced by the generalised Hu–Washizu variational principle [16]. Polynomials are often used as Trefftz functions [14]. Advantages of hybrid methods are: because the equilibrium is satisfied within the element, only boundary integrals are needed to compute numerically; high precision in the modelling of the local effects such as distribution of the stresses in the vicinity of the contact surface. We learned by studying hybrid methods also their disadvantages: significant errors were cumulated near the boundaries of the elements mainly due to the discontinuity in the primary variables, displacements, between the elements; a need to compute inversion for the singular matrix for each element. In the formulation proposed in this paper, the primary variable is continuous in the whole domain Ω despite the fact that we also have the domain Ω H , where there is defined equilibrium field uH , which is not a polynomial along the boundary of Ω H . Continuity is fulfilled by the special interpolating function uHh , which is defined for elements that

123

Comput Mech (2007) 41:159–174

share some nodes with the elements of the domain Ω H . The field uHh + uh defines approximation on those elements in such a way that continuity is preserved.1 Next, we use the standard variational principle [16,17], where the contact condition is enforced by the method of Lagrange multipliers. Later, the Lagrange multiplier is eliminated analytically. Outline of this paper is as follows. In the Sect. 2 we derive variational formulation FEM with enforced contact conditions and utilise the field uH that satisfy the Hertzian contact traction. In the Sect. 3 we deal with discretisation and linearˆ and derive iterative isation of the virtual work vector F(u) procedure for numerical solution of the contact by Newton method. In Sect. 4 we describe performed numerical experiments of selected problems with the aim to evaluate the accuracy of the proposed formulation. We discuss obtained results. Brief final conclusion is in Sect. 5 and the direction of the future research is drawn. Selected results of the Hertz theory for the case of line contact and uH field are summarised in the Appendix A. In this paper, we often use superscript for miscellaneous purposes not only as an exponent. To minimise ambiguities in formulas, when the superscript stands for exponent the base will be enclosed in parentheses, e.g. uu = (u)2 .

2 Variational formulation of the contact problem Energy functional of a deformable body Ω loaded on its boundary ∂Ωt can be written as

AdS, (1) L[u] = W dV − Ω

∂Ωt

where W is the strain energy density and A = t¯i u i is the work of the prescribed traction t¯ per unit length of the boundary ∂Ωt , assuming body in 2-dimensional space. The function u, the displacement field, is the independent variable of the functional. For linear elastostatics the stationarity condition of the L[u] with respect to u is

(2) t¯i δu i dS = 0, δL[u] = Ci jkl εi j δεkl dV − Ω

∂Ωt

which is the well-known variational principle used in FEM. εi j = (u i, j + u j,i )/2 are Cartesian components of the strain tensor, Ci jkl is the tensor of elastic constants of the material and t¯i is the prescribed traction on the boundary ∂Ωt . We use the Einstein summation convention. Subscript indices have values 1, 2 if another range is not specified. The indices in subscript after a comma denote partial derivatives with respect to the spatial coordinates. 1 We use uh to denote the interpolating functions on the finite element of the characteristic size h, which is common in the literature.

Comput Mech (2007) 41:159–174

161

2a

Let us consider several bodies in contact. Energy functional of this configuration will be the sum of the energy functionals of the bodies and new terms that enforce the contact conditions by method of Lagrange multipliers. We assume contact without friction so contact conditions will enforce only impenetrability of the outer boundaries of the bodies. Relative distance between the surfaces of two bodies Ω I and Ω J can be expressed with the help of the gap function g(x, t), for which holds: g = 0 in case that the point x (selected point on one surface so called contactor) is in contact with the other surface so called target; g > 0 in case that x is located outside the target body; g < 0 if x penetrated inside the target body. t is pseudotime. The gap function can be defined, following [9], as     g(x, t) = x − y (x, t) · en y (x, t) ,

e1n (0) h

Hh

H

H

Hh

h

Hh

H

H

Hh

h

Hh

Hh

Hh

Hh

h

h

h

h

h

h

h

h

h

h

∂Ωc

∂ Ω HE

∂ Ω HO

where en is the unit normal of the target profile at point y on the target profile which satisfies the condition

Fig. 1 Subdomains on the FEM discretised domain Ω

||x − y (x, t)|| = min{||x − y|| : y ∈ ∂Ω J }.

into the formulation. This field is described in Appendix A. Figure 1 shows one body in contact. Let us denote by Ω HI the domain, in which we assume the solution in the form

It can be shown that if there is unique y (x, t), then the gradient of the gap function is   (3) ∇g(x, t) = en y (x, t) . For the sake of clarity, in the following we will deal with the contact of two bodies only: Ω 1 , Ω 2 . Generalisation to the multibody contact is straightforward. Modified functional for two bodies is ⎡ ⎤


2  ⎢ ⎥ c I 2 L [{u } I =1 , λ] = AdS ⎦ − λgdS. ⎣ W dV − I =1

ΩI

∂Ω I

∂Ω c

The bodies in contact will be in equilibrium if the stationarity condition of the L c is met ⎡ ⎤

2  ⎢ ⎥ δL c = t¯i δu i dS ⎦ ⎣ σi j δεi j dV − I =1




− ∂Ω c

ΩI




δλgdS −

∂ΩtI

λδgdS = 0.

(4)

∂Ω c

We assume, that the boundary ∂Ω c , i.e. the set of points on the actual contact surface, is known, dependent on the primary variables of the problem, in fact ∂Ω c (λ). 2.1 Equilibrium functions in the approximation of displacements Now, let us include the equilibrium field  uHI = uHI x; a, K , K I , I = 1, 2

(5)

u I = uh I + uHI = Nh uˆ h I + uHI ,

(6)

where Nh is the matrix of standard FEM interpolating functions on the finite element, uˆ h I is the vector of nodal values of displacements. Note, that in this case the nodal value is not total displacement at its nodal point. The domain Ω HI is some close neighbourhood of the contact surface Ω c . We choose Ω HI in such a way that its boundary goes along the element edges. One part of this boundary, ∂Ω HEI inter(e)lement, separates it from the rest of the body and other part, ∂Ω HOI (o)uter, is part of the outer boundary of the body. Then we have ∂Ω HI = ∂Ω HEI ∪ ∂Ω HOI Certain part of the boundary ∂Ω HOI is actual contact strip ∂Ω c , with nonzero contact traction. It is important to choose Ω HI in such a way that ∂Ω HOI contains whole ∂Ω c , so ∂Ω HOI ∩ ∂Ω c = ∂Ω c . Two profiles are shown on Fig. 2 in two positions: as they are approaching and during contact. We assume that centres of the contact, s1 at ∂Ω 1 and s2 at ∂Ω 2 , are a priori known and independent of t, denoting selected material points. Similarly, we assume that unit normal vectors at the centres of the contact en1 and en2 are independent of t. The halfwidth of the contact strip is a, K 1 and K 2 are curvatures2 of the profiles at the centre of the contact in the referential configuration, t = 0, and K = K 1 + K 2 is the relative curvature. It seems that for our formulation the following definition of the gap function g is useful.     g y H , t = x2 y H , t − x1 y H , t · en1 2

We define that convex surface will have positive curvature.

123

162

Comput Mech (2007) 41:159–174

1/K2

1/K 1

For the domains Ω HI of both bodies, we can write the variation of the L cH as follows δL

e1t (0) x

1

H1

s (0)

1

2

s (0)

e2n (0)

e1n (0)

 u

hI

I =1

∂Ω HEI




yH1

−

Fig. 2 Contact profiles




gδλdS −

∂Ω c

   = x1 y H , t − x2 y H , t · en2 

 = u2 y H , t + x2 y H , 0     − u1 y H , t + x1 y H , 0 · en1 ,     where points x1 y H , 0 , x2 y H , 0 are the points at profiles in the referential configuration. The unit normal vectors, in case of the smooth profiles and K > 0 at the closest points of the approaching profiles and in case of contact at the centres of the contact, are en1 = −en2 . It means that this definition of g treats both profiles in the same way. There is no distinction between contactor and target. The whole situation is shown at Fig. 2. Thevectors {−en1 , −et1 } define local coordinates xH = x H , y H . We can write



For t = 0 last term is zero and  we have the gap in the referential configuration g y H , 0 . Last term in (7) is the difference between normal components of the displacements.3 For sufficiently smooth profiles, e.g. arcs, which are loaded a  insuch 3 H for way that a
1/K , we can neglect the term O y





I σiHI j δεi j dV

=

I σiHI j δu i, j dV Ω HI




tiHI δu iI dS ∂Ω HI

We neglect the difference between the tangential displacements.

123

−

I σiHI j, j δu i dV, (10)

Ω HI

where the last integral is zero because the field σ HI satisfies homogeneous equilibrium equations. Substituting from (10) into volume integral in (9), extending the domain of t¯ t¯ I |∂Ω HOI \∂Ω HOI ∪∂Ω cI  ≡ 0 t

and using the assumption that ∂ΩtHOI ∩ ∂Ω cI = ∅ we obtain cH

=

2 
 I =1

 dV σihj I δ εihjI + εiHI j

Ω HI




− 3




=

δL (8)

(9)

where the integral over ∂Ω HEI is the virtual work of traction tEI due to cohesive interaction with the rest of the I th body. On the boundary ∂Ω HEI we assume that condition of the continuity of the displacements is fulfilled. Straightforward numerical computation of the volume integrals in (9), e.g. using Gauss quadrature, demands to use high order formula compared with the order usually employed for standard finite elements with polynomial interpolation functions. The following procedure enables us to obtain boundary integrals from them. Using the symmetry of the tensor Ci jkl and Gauss theorem we can write the following

1 g y H , t = s2 (0) − s1 (0) · en1 + K y H 2      3 + u2 y H , t − u1 y H , t · en1 . +O y H (7)

λδgdS,

∂Ω c

Ω HI

2

y H ∈ −a, a. Then the variation of g is     δg y H , t ≈ δ u2 y H , t − u1 y H , t · en1   = δ uh2 + uH2 − uh1 + uH1 · en1

,λ

  tiEI δ u ih I +u iHI dS






I =1

 t¯iI δ u ih I +u iHI dS

−

−





2

Ω HI




e2t (0)

∂ΩtHOI



+u

HI

2 
   hI HI σihj I + σiHI δ ε dV + ε = j ij ij

2

s (t) = s (t)

2a

cH

∂Ω HOI \∂Ω cI

 t¯iI − tiHI δu iI dS

Comput Mech (2007) 41:159–174






−

163

tiEI − tiHI δu iI dS

solution in the form

∂Ω HEI


 I I I HI I −λn i − ti δu i dS u i n i δλdS−




+ ∂Ω cI








(11)

∂Ω c

  where we substituted for g y H , t from (7), for δg from equaT  tion (8) and used notation enI = n 1I , n 2I for vector of the unit normal. We see that the Lagrange multiplier λ stands for the contact traction. Next, using the symmetry of Ci jkl and then Gauss theorem we get this identity

h I HI σi j δεi j dV = εihjI δσiHI j dV Ω HI

Ω HI




= Ω HI




u ih I δtiHI dS.

(12)

∂Ω HI

After substituting (12) into (11) and rearranging terms we finally obtain δL

cH

=

2 
 I =1

σihj I δεihjI dV

Ω HI




  t¯iI − tiHI δu iI − u ih I δtiHI dS

− ∂Ω HOI \∂Ω cI




∂Ω HEI

+




u iI n iI δλ + u ih I δtiHI dS

∂Ω cI

−




∂Ω cI




−

−λn iI

− tiHI







 δu iI dS

g y , 0 δλdS. H

(15)

• uHh (x) = uH (x) along such element edge, which is part of the ∂Ω HE ; • uHh (x) = Nh (ζ )uH (ˆx) along such element edge, which has only corner node xˆ ∈ ∂Ω HE ; • uHh (x) = 0 along such element edge, which does not share nodes with ∂Ω HE . Such a field can be constructed on the given 2-dimensional Lagrange element as a tensor product of the standard 1-dimensional interpolating functions with the field

and possibly also uH (ˆx)-products of the interpolating functions of a certain corner nodes xˆ of the element. This procedure is in analogy with the construction of the classical interpolating functions for standard 2-dimensional Lagrange element used in the domain Ω h . We will not go to details of their derivation for the sake of brevity. In order to be compatible with a neighbouring element from Ω H , δuHh depends on δuH so existence of the field uHh does not enlarge the number of free parameters of the variational problem. We can also write   uHh x; a, K , K 1 = uHh Nh (ζ )ˆx; a, K , K 1 , where ζ = (ξ, η) are the natural coordinates of the given isoparametric element. Finally, we observe that because of fields uH and uHh elements of the domains Ω H and Ω Hh are not isoparametric. For the above mentioned rest of the body, Ω Hh ∪ Ω h , the variation of the energy functional can be written as



σihj δεihj dV − t¯i δu ih dS δL cHh uh , uH =

  tiEI − tiHI δu iI − u ih I δtiHI dS

−

on the domain Ω ,

u=u

h

uH |∂Ω HE

u i,h Ij δσiHI j dV

=

(14)

h

where the field uHh satisfies the compatibility conditions

∂Ω cI

g y H , 0 δλdS,

−

u = uh + uHh on the domain Ω Hh ,

Ωh

(13)

∂Ω c

∂Ωth


  σihj + σiHh δ εihj + εiHh dV + j j Ω Hh




2.2 Compatibility in displacements and formulation on the rest of the body

−

In this section we will deal exclusively with the I th body so we will drop the redundant index I . We can partition the rest of the body as follows, see Fig. 1. Elements that are outside of Ω H but share some nodes with those in Ω H belong to the domain Ω Hh The remaining  Ω belong to the  elements of domain Ω h , so Ω h = Ω \ Ω H ∪ Ω Hh . We assume the

−

∂ΩtHh




 t¯i δ u ih + u iHh dS  tie δ u ih + u iH dS,

(16)

∂Ω HE

where te , defined on the ∂Ω HE , is a traction due to the cohesion with the Ω H domain. Obviously te = −tE . It is assumed that displacements are continuous along ∂Ω HE .

123

164

Comput Mech (2007) 41:159–174

2.3 Formulation of the contact problem

2.4 Elimination of the Lagrange multiplier

Final form of the stationarity condition (4) is

Let us consider domains Ω QI that consist of elements, whose edges contain some part of the contact strip ∂Ω c . Symbol Ω QI stands for the e(q)uilibrium domain as will be seen below. For the domains Ω QI the condition of equilibrium is ⎡ ⎤

2 ⎢ ⎥ σiIj δεiIj dV − A I dS ⎦ ⎣

δL c =

2 
 I =1

+


σihj I δεihjI dV −

ΩI

2 
 I =1

Ω Hh I

∂ΩtI

 I Hh I I hI σiHh + σiHh j δεi j j δεi j

I dV + +σihj I δεiHh j




+

∂Ω cI




−

∂Ω cI




−

 t¯iI δu iI dS

I =1


 tiHI δu iI + u ih I δtiHI dS

u iI n iI δλ + u ih I δtiHI dS  −λn iI − tiHI δu iI dS (17)

which was obtained be adding the contribution (13) and two contributions (16); and further rearranging terms with the aim to show how contact condition modifies the standard FEM. We note that: • Terms in the first sum are almost identical as in case of the standard FEM without contact, compare with (2). The only difference is that in (17) the field u I is defined differently on different parts of the boundary ∂ΩtI : equations (6), (14) and (15). • Next term is the contribution of the non-equilibrium and non-isoparametric approximation on Ω Hh to the volume integral and so it is needed to compute it only for the elements of Ω Hh . • Remaining terms are due to contact conditions and due to rewriting of some volume integrals in terms of the boundary ones, see Sect. 2.1. Direct numerical computation of the integrals on ∂Ω c is ineffective because the contact strip, on which the contact traction is applied, is changing with the load ∂Ω c = ∂Ω c (λ). We know that λ has the physical meaning of the contact traction. According to Hertz theory this traction is elliptic distribution of pressure (30). Variation of this pressure with respect to change in amplitude of the load is discontinuous at the endpoints of the contact strip. On the other hand the equilibrium field uH satisfies boundary condition of the Hertzian pressure along the line boundary. In the following we will use this fact and will focus on the elimination of λ and on analytic solution of the integrals on the boundary ∂Ω cI .




λδgdS −

∂Ω c

∂Ω c

123

−

∂Ω HI \∂Ω cI

 g y H , 0 δλdS = 0,

Ω QI




∂Ω QEI

δλgdS = 0,

(18)

∂Ω c

where u iI = u iHI + u ih I and A I is virtual work per unit length due to cohesive forces on the boundary ∂Ω QEI , which separates Ω QI from the rest of the I th body. We assume that displacements are continuous on ∂Ω QEI . Integrating volume integrals by parts we get ⎤ ⎡

2  ⎢ ⎥ A I − tiI δu iI dS ⎦ σiIj, j δu iI dV − ⎣− I =1

Ω QI

∂Ω QEI



− −λn iI − tiI δu iI dS − δλgdS = 0, ∂Ω c

(19)

∂Ω c

T  where tiI = σiIj n Ij and enI = n 1I , n 2I are the components of the unit normal vector of the outer boundary. We can see that if the displacement field uiI on Ω QI is equilibrium then the volume integral will be identically zero. If we define the virtual work A I as follows A I = tiI δu iI , on ∂Ω QEI ,

(20)

also second integral in the sum in (19) will vanish. Then satisfying condition of the impenetrability of the surfaces of the bodies  g y H , t = 0, on ∂Ω c , (21) the condition (19) will be reduced to the form
 −λn iI − tiI δu iI dS = 0, ∂Ω c

which can be satisfied if we put λn iI = −tiI at each point of ∂Ω cI . In other words, following this procedure we satisfied the equilibrium in the strong sense on the Ω QI domain. According to Hertz, the distribution of the contact pressure is given by (30). Equilibrium field uHI satisfies this boundary condition on the line shaped boundary ∂Ω c . Furthermore tHI is zero outside the contact strip along the whole boundary of the halfspace so that we can extend the domain ∂Ω cI in the integrals (17) and (19) also outside the contact strip so that it will be equal to the whole (o)uter boundary, ∂Ω QOI , of the

Comput Mech (2007) 41:159–174

165

Ω QI . However, the shape of the ∂Ω QOI has to be linear, and because of tiI = tiHI +tih I to comply with the Hertz boundary condition we get the following conditions tih I = 0 λn iI

−tiHI

=

on ∂Ω QOI ,

(22)

on ∂Ω

(23)

QOI

.

Details of the derivation of the equilibrium polynomial field uh I , satisfying boundary condition (22) are not shown in this paper. The procedure is based on [14]. Finally, the equilibrium condition (17) can be rewritten into the form

 δL c uh1 , uH1 , uh2 , uH2 =

2 
 I =1

+

∂ΩtI

 I Hh I I hI h I Hh I σiHh dV + σiHh j δεi j j δεi j + σi j δεi j

Ω Hh I




+

 tiHI δu iI +u ih I δtiHI dS

∂Ω HI \∂Ω QOI






−

t¯iI δu iI dS

−

ΩI

2 
 I =1




σihj I δεihjI dV

u iHI δtiHI dS


+

 g y H , 0 n i1 δtiH1 dS = 0,

∂Ω QO

∂Ω QOI

(24) assuming that the conditions (20), (21), (22), (23), and that uh I is equilibrium on the Ω QI are satisfied. Integrals on the boundary ∂Ω QOI , with the use of (7), can be calculated analytically.

3 Numerical formulation of the contact problem by Newton method In this section we will outline a numerical procedure that can be employed to obtain the numerical solution of the contact problem. The procedure is based on the standard isoparametric formulation of the FEM [16]. For the sake of simplicity we will focus on the case where the approximate solution in the Ω H is

written in the matrix notation, is  T ˆ = 0, uˆ = uˆ h , b , δL c = δ uˆ T F(u) ˆ is the virtual work vector. The scalar product of where F(u) two symmetric tensors which occurs in (24) can be expressed in matrix notation as σi j δεi j = δ T τ , where τ = (σ11 , ˆ can be written σ22 , σ12 )T ,  = (ε11 , ε22 , 2ε12 )T . Then F(u) by substituting the following formulas into the (24) uh uH uHh h H  Hh τ t

= Nh uˆ h , = uH (x; b, K , K 1 ), = uHh (ζ ; b, K , K 1 ), = Bh uˆ h , =  H (x; b, K , K 1 ), =  Hh (ζ ; b, K , K 1 ), = C, = ZC,

H

h h

   T  δL c = δ uˆ h , b F uˆ h , b = δ uˆ hT F1 uˆ h +δbF2 uˆ h , b +δ uˆ hT F3 (b) + δ uˆ hT F4 + δbF5 (b). ˆ is i.e. vector F(u)  F uˆ h  + F (b) + F  1  4 h 3 ˆ = F uˆ , b = , F(u) F2 uˆ h , b + F5 (b)

(25)

where  2 
  h hT h h F1 uˆ = B CB dV uˆ ≡ Khh uˆ h , I =1

 F2 uˆ h , b =

ΩI

2 
 I =1

BHhT CBh dV uˆ h

Ω Hh I






+

H

where the equilibrium field uH will depend only on one parameter: the square of the halfwidth of the contact strip b = (a)2 , see (5) and Appendix A. The use of b instead of a is more convenient because the contact force is proportional to b, see (29), similarly as load is proportional to displacement in linear elasticity. Condition of the equilibrium (24),

= Nh δ uˆ h , = ∂b uH δb = NH δb, = ∂b uHh δb = NHh δb, = Bh δ uˆ h , = ∂b  H δb = BH δb, = ∂b  Hh δb = BHh δb, = Cδ, = ZCδ,

where in the last two rows we assume  =  h ,  H ,  Hh . C is the matrix of the elastic coefficients and Z is the so called matrix of the unit normal. With respect to the dependence on the primary variables T  uˆ = uˆ h , b , the variation of the functional L c can be expressed in this form

u = u + u = N uˆ + u (b), h

δuh δuH δuHh δ h δ H δ Hh δτ δt

HT

B

CZ N dS uˆ T

∂Ω HI \∂Ω QOI

F3 (b) =

2 
 I =1

F4 = −




BhT C Hh dV +

Ω Hh I

2


h

h

= ∂b F3T uˆ h ,

 NhT ZC H dS ,

∂Ω HI \∂Ω QOI

NhT t¯ dS,

I =1 I ∂Ωt

123

166

Comput Mech (2007) 41:159–174

F5 (b) =

2  




I =1




NHhT t¯ dS −

−

∂ΩtHOI

∂ΩtHh I




+

Iterative procedure will be   K b(i−1) ∆uˆ (i) = −F uˆ (i−1) , i = 1, 2, 3 . . . .

NHT t¯ dS

We note the following facts.

NHT ZC H dS + α(b),

∂Ω HI \∂Ω QOI

• The matrix K is symmetric as can be expected. • Furthermore K is dependent only on the primary variable b, K = K(b), where dependence on b is only in the last row and column of the K. This can be exploited in the effective solution procedure for the resulting equations (28). • The matrix Khh is the standard stiffness matrix for the isoparametric FEM without the contact conditions and equilibrium fields, which can be obtained as a sum of the stiffness matrices for all the elements regardless from which domain they are Ω H , Ω Hh or Ω h . Only exception is that in Ω Q the field uh has to be equilibrium.

where
2 

α(b) = −

uHI · ∂b tHI dS

I =1 QOI ∂Ω




 g y H , 0 en1 · ∂b tH1 dS

+ ∂Ω QO 2 




= ∂b

uh I · tHI dS.

I =1 QOI ∂Ω

The term α(b), using (7), neglecting the term O α(b) = −

π K d0 4



(1 − (ν1 )2 ) (1 − (ν2 )2 ) + E1 E2



yH

3

is

−1 ,

(26)

where E i and νi are Young modulus and Poisson number for ith material, respectively. Initial gap is denoted by d0 . It results that α is independent of b, in fact. ˆ = 0 by We solve the set of the nonlinear equations F(u) the Newton method. For this purpose we need to compute the ˆ Using the special Hessian of the energy functional, ∂uˆ F(u). ˆ see (25), we can write structure of the vector F(u), ˆ = (∂uˆ F)(b) ∂uˆ F(u)   ∂b F3 (b) Khh ≡ K(b), (27) = ∂b F3T (b) ∂bb F3T (b)uˆ h + ∂b F5 (b) where ∂bb F3 (b) =

2 
 I =1

Ω Hh I

∂b F5 (b) =

NhT ZC∂b BH dS ,


2   − I =1

∂ΩtHh I




+

∂b NHhT t¯ dS −

Initial approximation can be uˆ = 0, including b. Nonzero initial penetration d0 < 0 is prescribed.




To demonstrate the formulation and to evaluate its accuracy the following problems were considered. ∂b NHT t¯ dS

∂ΩtHOI

  ∂b NHT ZC H + NHT ZCBH dS .

∂Ω HI \∂Ω QOI

Note, that using b instead of a in discretisation results in ∂b α(b) = 0, as can be seen from (26), whereas if we had α(a) then ∂a α(a) = 0.

123

1. For each element compute standard element stiffness matrix Kehh and r.h.s. vector −F1e − F4e , see (25). 2. If the element belongs to the domain Ω Q , eliminate some of its uˆ h DOFs in order to make it equilibrium, thus obtaining reduced equilibrium Kehh and Fe . 3. If the element corresponds to some of Ω HI , Ω Hh I , ∂Ω HI , ∂Ω Hh I then corresponding contributions F3 (b), F5 (b) and their partial derivatives are calculated and added according to (25) and (27). 4. Assembly the global system K, F and perform the iteration. Repeat steps from the 3rd one until convergence is reached.

4.1 Problem descriptions






∂Ω HI \∂Ω QOI

Flowchart of a typical implementation is as follows

4 Numerical examples

BhT C∂b BHh dV

+

(28)

Block: This problem represents the case of squeezing of the infinitely long foundation of the height D with the array of rigid cylinders arranged periodically in the horizontal direction with the period 2D. We are interested in contact half-strips that are much smaller than D so the solution of the problem should approach uH solution. Geometry, constraints and loading can be seen from Fig. 3, where it is seen that only half-period, in

Comput Mech (2007) 41:159–174

167

Z

TIME 1.000 C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C B B B B B B B B B B B B B B B B B

PRESCRIBED DISPLACEMENT TIME 1.000 0.5000 B B B B B B B B B B B B B B B B B

U U 2

BC-

Z

TIME 1.000 Y

X

3

-

Constraints are – B: ux = 0, uy is free. – C: ux = uy = 0. The rigid cylinder is fully constrained. Constant displacement uy is prescribed at the bottom of the considered domain. Note, that in the figure, produced by ADINA, y, z coordinates were used in place of our x, y.

Fig. 3 Block: geometry and boundary conditions

the horizontal direction, of the foundation was discretised. Dimensions of the modelled square block were D = 10 × 10. Cylinder: An elastic cylinder, only 1 quarter of which was discretised, was squeezed between two equal rigid cylinders. The radius of the elastic cylinder was fixed, R = 10, but the radius of the rigid cylinders was varied with the aim to explore the different relative curvatures K . Geometry, constraints and loading are given in Fig. 4. For this problem the approximate analytic solution can be constructed, according to [5], as two opposite fields uH applied at the contact strips with the relative distance of the centres of the contact 2R and the isotropic pull with the amplitude σ = P/(π R), where P is the resulting force of the contact load. In the limit, where elliptic contact traction approaches the Dirac distribution Pδ(x) this isotropic pull causes that on the circumference of the cylinder outside the contact strips is zero traction as is needed for the free surface. Beam: The beam is supported in such a way that its left and right end can freely deform in y, vertical, direction while in x direction the prescribed displacement is applied. The beam bends over fixed rigid cylinder, see Fig. 5. Due to symmetry, only half of the beam was discretised, with the horizontal dimension 10. In our computations,

C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C B B B B B B B B B B B B B B B B B B B

X

Y

PRESCRIBED DISPLACEMENT TIME 1.000 1.700

B C

U U 2 3 - -

Constraints B ,C are the same as in Fig.3. Constant displacement uy is prescribed at the right side of the beam. Note,that in the figure, produced by ADINA, y, z coordinates were used in place of our notation x, y.

Fig. 4 Cylinder: geometry and boundary conditions

the height h of the beam was a parameter that assumed selected values, see the Sect. 4.2 for the specific values. This problem demonstrates the use of our formulation in the case, where locally, in the neighbourhood of the contact, there exists also significant strain due to another reason than the contact loading: the bending loading in this case. For all described problems we used Young modulus E = 1, Poisson number µ = 0.3 and plane strain in our computations. For all problems the calculations by the proposed formulation were done for mesh configuration with only Ω H domain, the boundary ∂Ω QO was defined by the minimum number of the element edges so that the whole contact strip was contained in Ω QO also for the largest given displacement of the rigid body. For the beam problem we have done also the variant, with Ω Hh , Ω h , and corresponding boundary ∂Ω HE as will be introduced later in the Sect. 4.2.4. 4.2 Evaluation of the numerical solutions 4.2.1 Criteria used in evaluation In this section we compare the numerical solutions obtained by the proposed formulation with numerical solutions obtained by: the well-known ADINA software system, version 8.3, that implements state of the art node-to-segment contact algorithm based on the constraint function method

123

168

Comput Mech (2007) 41:159–174

TIME 1.000

Z PRESCRIBED X DISPLACEMENT TIME 1.000

Y

2.900 C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C B B B B B B

B C

U U 2 3 - -

Constraints B, C are the same as in Fig. 3. Constant displacement uy is prescribed at the right side of the beam. Note, that in the figure, produced by ADINA, y, z coordinates were used in place of our notation x, y.

Fig. 5 Beam: geometry and boundary conditions

[9]; and approximate analytic solution of the problem cylinder, described in the Sect. 4.1. For all problems we computed the displacements, and derived quantities such as P, τ1max , and σ3min for selected values of prescribed relative curvature K and prescribed displacement v at the contact centre. In this section, by solution, we mean the numerical solution obtained by the proposed formulation. We used the following criteria for evaluation of the accuracy of the obtained solutions. • Relative error in the displacements in the nodal points for two solutions: evaluated one and referential one is ∆u =

 2 N    ˆ i − ur xˆ i i=1 u x ,  2 N ˆ x u r i i=1

  where u xˆ i is the computed value of the displacement   in the nodal point, ur xˆ i is the value of the referential solution at the same given nodal point. The sum is only over N nodes in the region close to the contact strip so that the influence of the farther boundary conditions is minimal and resolution for the evaluation of the local quality of the solution is higher. As a referential solution can be used the numerical solution obtained by the ADINA or solution in a closed form. • Important criterion for the evaluation of the solution for all solved problems is the total amount of the strain energy

123

accumulated in the elastic body. Exact solution of the problem has the property of minimal accumulated strain energy. For our cases the strain energy is equal to the work of the contact traction during deformation. Displacements of the contact strip are given by the prescribed displacement of the rigid cylinder while squeezing the elastic body. Then good measure of the work of the contact traction is the amplitude of contact force P. Because in our formulation the impenetrability is exactly satisfied the lesser the obtained P, the better the solution. • The dependence of the contact force on the displacement P(v) is called the stiffness characteristics of the structure. We obtained the stiffness characteristics for each solved problem for ten evenly selected values of K from 0.02 up to 0.2 and for ten evenly selected values of v from 0.1vmax (K ) up to vmax (K ). In this way we obtained the dependence P(K , v). The maximum displacement vmax (K ) was chosen in such a way that a of the developed contact strip should not be higher than some given value (1.5 for all problems). Due to this reason the maximum displacement for the smaller values of K is smaller than the maximum for the larger curvatures as can be seen from the figures. • Next evaluated quantity is the largest maximum shear stress in the vicinity of the contact strip τ1max = max(τ1 ) = max(σ1 − σ3 ), where σ1 , σ3 denotes first and third principal stress, respectively. τ1max is an important quantity for the contact problems in general. On the basis of the known solution uH , given in Appendix A, we know that it is located below the centre of the contact in the dis. tance x H = 0.786151a. For the problems of block and cylinder τ1max ≈ τ1Hmax especially in case of the cylinder. Because of that we used τ1Hmax in our evaluation, which was computed for the same values of K , v as force P. For the beam we did not evaluate τ1max . • Next quantity was σ3min , which is the minimum 3rd principal stress. For the field uH it is the maximum contact pressure (with the opposite sign) located at the centre of the contact. As in case of τ1max also σ3min was evaluated only for block and cylinder for the same K , v.

4.2.2 Block Results for block are presented for these cases: configuration of the mesh containing only Ω H that has 8×8 linear elements as is shown on Fig. 3; and configuration with 4×4 9-node elements, where also all of them belong to Ω H . The dependence of the relative error of displacements, ∆u(K , v), is shown on Fig. 6. The solution by ADINA for the same mesh is used as a referential solution. K is on the first axis, v is on the second axis. Maximal values of error are for the smallest v. These are due to division

Comput Mech (2007) 41:159–174

169 0.35 0.3 0.25 0.2

0.4

0.15

0.3 0

0.1

0.05

0.05

0.2 0.1

0.1

0 0

0.15

0 0

0.05

0.1

0.15

0.2

0.25

0.3

a) 4-node Lagrange elements

0.35

0.4

0.1

0.2

0.2 0.05

0

0.1

0.15

0.2

0.25

0.3

0.35

0.4

a) 4-node Lagrange elements 0.35 0.3 0.25 0.4

0.2 0.3 0

0.15

0.05

0.1

0.2 0.1 0.1 0 0

0.05

0.15 0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.05 0 0

0.2

b) 9-node Lagrange elements

0.1 0.2

For both pictures: 1st axis: K, 2nd axis: v, the vertical axis: ∆ u(K, v).

Fig. 6 Block: displacement error ∆u

by small numbers: the total referential displacements. With the increase of v error increases, it reaches a maximum, and decreases or converges to some lower value. Increase and decrease is caused by the interaction of the contactor node with target segment in ADINA. This can be seen also on Fig. 7, where stiffness characteristics are shown and it was confirmed also after reviewing of results from ADINA. On the Fig. 7a, for linear 4-node elements, the top smooth characteristics is obtained by ADINA for very fine mesh of 100×100 elements, close to it is characteristics obtained by the proposed formulation. Lower characteristics with well at about v = 0.7vmax (K ) is the solution by ADINA at the same mesh as was used for proposed formulation. On the first sight results are against our expectation: weakest characteristics should be that closest to exact. But if we take into account that node-to-segment contact algorithm in ADINA prevents to penetrate into target only at the nodal points of contactor while element edges of contactor can freely penetrate, we can understand why coarser model is more weaker. This effect

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

b) 9-node Lagrange elements For both pictures: 1st axis: K, 2nd axis: v, the vertical axis: P(K, v).

Fig. 7 Block: stiffness characteristics P(K , v)

diminishes for the finer mesh so the stiffness is higher for the finer mesh. Finally with further refinement the characteristics is expected to slightly decrease and converge to the exact solution. From this viewpoint our solution should be the most accurate one. Error of the maximum contact force of our solution for 9-node elements with respect to our solution for 4-node elements was about 0.3%. For 9-node elements were again the stiffest and the weakest characteristics that obtained by ADINA. τ1max (K , v), σ3min (K , v) are shown on Figs. 8, 9, respectively. In all these figures the smoothest surface is obtained by our formulation, slightly rough is that obtained from ADINA for finer mesh and more rough is by ADINA for the first mesh. The uH field is almost the exact solution of this problem. We see, that by adding uH field into the approximation in our

123

170

Comput Mech (2007) 41:159–174

0 −0.02

0.06

−0.04

0.05

−0.06 −0.08

0.04

−0.1

0.03

0

0.02

0.05 0.1

0.01

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

a) 4-node Lagrange elements

−0.14

0

−0.16 0.1

0.15

0

−0.12

−0.18 0

0.2

0.05

0.1

0.4

0.15

0.2

0.25

0.3

0.2 0.35

0.4

a) 4-node Lagrange elements

0.07 0

0.06

−0.02

0.05

−0.04

0.04

−0.06

0.03

−0.08

0.02

−0.1

0.01

−0.12

0

−0.14

0

−0.16

0.05

−0.18

0

0.1 0.2

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

b) 9-node Lagrange elements For both pictures: 1st axis: K, 2nd axis: v, the vertical axis: τ1max (K, v).

Fig. 8 Block: stresses τ1max (K , v)

formulation we reach high accuracy even for relatively small number of elements. The error ∆u on the Fig. 6 is in fact the error of the solution by ADINA with respect to our solution, both for the same mesh. 4.2.3 Cylinder The results presented for the cylinder are for the configuration of the mesh as is sketched on Fig. 4, i.e. linear elements. 9-node elements give qualitatively the same results. The error ∆u is shown on Fig. 10. As the referential solution we used approximate solution described in Sect. 4.1. The error is in the whole region under 5%. Largest error is for largest v, where the analytic solution itself deviates from the exact solution mainly in the region in the vicinity of the contact strip. The stiffness characteristics P(K , v) is shown on Fig. 11. Because for this problem we computed ADINA solution for larger vmax (K ) also the shown characteristics is for different

123

0.1

−0.2

0.15

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.2

b) 9-node Lagrange elements For both pictures: 1st axis: K, 2nd axis: v, the vertical axis: σ3min (K, v).

Fig. 9 Block: stresses σ3min (K , v)

range of v. Characteristics from our formulation and from analytic solution are almost indistinguishable in the figure. τ1max (K , v), σ3min (K , v) are shown on Figs. 12, 13, respectively. Displayed data shows analogical properties as that for block. We conclude that our solution is very close to analytic one as can be expected. But in fact it is not obvious because our solution was obtained by solving the FEM equilibrium equations on the one quarter of cylinder and analytic solution was obtained using the function of one variable a, and from the condition that relative compression of the opposite points on the cylinder has to be equal to prescribed value 2v we determined the value of a. Analytic solution converges to exact one for a/R → 0. Good agreement can be seen also with ADINA solution. 4.2.4 Beam The results presented on Figs. 14 and 15 are for the mesh configuration of 8×2 9-node elements. Solutions obtained

Comput Mech (2007) 41:159–174

171

0.05

0.12

0.045

0.1

0.04

0.08

0.035

0.06

0.03

0.04

0.025

0.02

0.02 0

0 0

0.05

0.05

0.1

0.1

0.15 0.2

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.15

0.4

0.2

1st axis: K, 2nd axis: v, the vertical axis: ∆ u(K, v).

0.2

0

0.4

0.8

0.6

1.6

1.4

1.2

1

1st axis: K, 2nd axis: v, the vertical axis: τ1max (K, v).

Fig. 10 Cylinder: ∆u, 4-node Lagrange elements

Fig. 12 Cylinder: stresses τ1max (K , v), 4-node Lagrange elements

1.4 1.2

0

1

−0.05

0.8

−0.1

0.6

−0.15

0.4

−0.2

0.2

−0.25

0 0

−0.3 0.1 0.2 0

0.5

1

1.5

2

−0.35 0

1st axis: K, 2nd axis: v, the vertical axis: P(K, v).

Fig. 11 Cylinder: P(K , v), 4-node Lagrange elements

0 0.05 0.1 0.15

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

0.2

1st axis: K, 2nd axis: v, the vertical axis: σ3min (K, v). Fig. 13 Cylinder: stresses σ3min (K , v), 4-node Lagrange elements 0.06

by presented formulation and by ADINA are shown. Results for 4-node linear elements are given in Table 1 only. ∆u for beam is shown on Fig. 14. There are 4 functions ∆u(K , v), each for one height h of the beam. They look the same only slightly vertically shifted. Largest error is for the most thick beam, then follows more thin, more thin and finally the most thin one. The more thick the beam, the larger amount of strain due to contact compared with strain due to bending, and vice versa. We guess, that small increase of error with the increase of h can be due to increase of error in contact strains in ADINA solution. On the other hand for linear elements error increases with decreasing of h.4 This can be explained in the way, that for linear elements, see discussion of Table 1 below, our formulation is not able sufficiently 4

The corresponding figure is not included in the paper.

0.05 0.04 0.03 0.02 0.01 0 0 0.05 0.1 0.15 0.2 0

0.5

1

1.5

2

2.5

3

1st axis: K, 2nd axis: v, the vertical axis: ∆ u(K, v). Fig. 14 Beam: ∆u, 9-node Lagrange elements

123

172

Comput Mech (2007) 41:159–174

Table 1 Contact forces P(K , v) = P(2/10, 2.8) for beam Parameters of

h=1

the problem

P

PA

P

PA

P

PA

P

PA

4, 16×4, Ω H

0.109314

0.007037

0.124951

0.012970

0.149106

0.028989

0.206134

0.088359

4, 16×4,

Ω Hh

h = 10/8

0.118868

4, 160×40

h = 10/6

0.133092

0.156208

0.006056

9, 8×2, Ω H

0.008751

9, 8×2, Ω Hh

0.008816

9, 32×8

0.011679

0.006121

h = 10/4

0.015419

0.011780

0.015504

0.210399 0.026989

0.031968

0.027119

0.032079

0.006050

0.011675

0.084846 0.086663

0.085154

0.086645 0.026982

0.084864

Parameters of the problem: (number of nodes of the used element; number of elements: horizontally × vertically; domains that are used in the model); h height of the beam, length of the modelled part of the beam was 10; P is the contact force computed by the proposed formulation; PA is the component of the contact force in the direction of the normal to the contact profile computed by the ADINA

0.09

32×8 for the quadratic 9-node Lagrange elements), see rows 3 and 6 in the table. Let us remind that exact value of the contact force should be the smallest one, see also Sect. 4.2.1. Discussion of the results:

0.08 0.07 0.06 0.05 0.04 0.03 0.02 0.01 0 0 0.1 0.2

0

0.5

1

1.5

2

2.5

3

1st axis: K, 2nd axis: v, the vertical axis: P(K, v). Fig. 15 Beam: P(K , v), 9-node Lagrange elements

to capture the bending deformations in the vicinity of the contact strip. Qualitative insight about stiffness characteristics for different h can be gained from Fig. 15. Computed values can be found in Table 1. On Fig. 15 the stiffness characteristics for three values of h are shown, for the clarity the case of the most thin beam is not shown. For each h two functions P(K , v) can be seen that are close to each other. Slightly more stiff is always that computed by our formulation. We observe good agreement with the ADINA solution on the same mesh especially for the most thick beam. The contact forces for the beam are summarised in the Table 1. All of them are computed for K = 2/10, v = 2.8, the largest values for relative curvature and displacement used in the computation of the stiffness characteristics of the beam. As a good approximation to the exact solution we can assume the solution obtained by the ADINA program for the finer mesh (160×40 for the linear 4-node Lagrange elements and

123

• First row of Table 1 corresponds to the configuration 16×4 linear elements, all in Ω H . This case is shown on the Fig. 5. Linear equilibrium field uh used in Ω Q (first two elements along the contact profile) satisfying (22) represents only pure pull in the tangential direction to the contact profile, σth = const. Because the field u h has to have continuous first derivatives at any point in Ω Q , pure pull σth of the same magnitude has to be defined on both elements of the Ω Q . This approximation is not sufficient to capture the real strain state in the vicinity of the contact strip, where also significant bending strains are present, so that computed value P is very different from the exact one especially for more thin beams. • Fourth row in Table 1 corresponds to configuration 8×2 of the 9-node Lagrange elements, all belonging to Ω H . Here the equilibrium domain Ω Q was only one element in the direct contact with the rigid cylinder. The uh field also satisfied the symmetry condition about the axis given by normal at the contact centre. However, in this case the equilibrium uh field compatible with the 9-node Lagrange element contains also bending term which was able better to capture the strain state in the region. For the most thin beam, h = 1, where the deformations are mainly due to the bending and contact force is small we get the relative error 43%. For the thickest beam, h = 10/4, where the contact strains are significant relative error in contact force is under 2% compared with the solution obtained by the ADINA for the same configuration of the mesh. • Second row corresponds to the configuration 16×4 linear elements. As was mentioned above the mesh is shown in Fig. 5. In this case, however, only first two columns (eight elements) under the contact profile belong to Ω H

Comput Mech (2007) 41:159–174

then follows the column of elements that make up domain Ω Hh and rest of the elements belong to Ω h . Equilibrium domain Ω Q was the same as for the case without Ω Hh . It can be seen that presence of Ω Hh worsens the solution mainly for more thin beams, where relative error with respect to solution without Ω Hh is under 9%, for the case of the most thick beam it is about 2%. • Fifth row corresponds to the configuration 8×2 9-node Lagrange elements. First column of two elements belongs to Ω H , then follows one column of Ω Hh elements and finally domain Ω h . Equilibrium domain Ω Q was the same as in the case without Ω Hh . Relative error with respect to solution without Ω Hh was for all heights h of the beam under 1%. Moreover for the most thick beam we obtained slightly better solution than in the case without Ω Hh . We conclude that in this case Ω Hh has not significant influence on the solution. 4.2.5 Common remarks Note, that results of the computation by the ADINA shows that the impenetrability condition was not always exactly satisfied. The node on the block, cylinder or beam, the contactor in ADINA, located at the centre of the contact strip sometimes penetrated into the interior of the rigid cylinder, the target. Furthermore, as was mentioned earlier, node-to-segment algorithm implemented in ADINA will not consider as a contact if edges of the contactor penetrate into the target. Only position of the nodes of contactor is checked against the edges of the target elements. These facts contributed to smaller computed contact force P and led to apparently better results. Finally, note that in all cases of solved problems, solution obtained by the proposed formulation converged almost always within 6–14 full Newton iterations. For simplicity and space constraint of our presentation, all demonstrated problems are a contact of a elastic body with a rigid body. A more general contact of two elastic bodies was also solved including limiting cases where one body approaches a rigid body E 1 E 2 , for example E 1 = 1010 E 2 . Results agree with the results presented in this paper. Number of full Newton iterations were again within the above mentioned range.

5 Conclusion In this paper we described and demonstrated the main features of a new formulation of the contact problem within the framework of the FEM assuming: small displacements, rotations and strains; line contact of the bodies, i.e. the plane strain or plane stress; and validity of Hertz theory.

173

From the discussion in Sect. 4.2 we conclude that our formulation is more accurate for problems, where Hertz contact dominates the strain state, especially for small number of elements, compared with standard node-to-segment contact algorithm implemented in ADINA. On the other hand, for more general problem of bending of the beam, results, for the case of second order approximation uh in Ω Q , are in fairly good agreement with ADINA only if the thickness of the given beam was large enough. More detailed study of this more general case is needed. Further research is needed mainly in the following directions. • To extend applicability of the formulation by dropping some of the above mentioned restrictions and assumptions; • To improve the formulation for the cases where other than contact strains are significant in the contact region; • To investigate the numerical properties of the formulation including stability, consistency and convergence.

Appendix A: Some relations of the Hertz theory and the field uH In this appendix, superscript will be used as an exponent as is convenient. Exception is the letter H which marks the quantity that has some connection to the Hertz theory. On the basis of the assumptions of the Hertz theory we can derive for the case of line contact these formulas. • Halfwidth of the contact strip is      4P 1 − µ2 1 − µ22 1  + , a= πK E1 E2

(29)

where P is the contact force on the unit of the length of the line of the contact, E i , µi are elastic coefficients for the ith profile. • Hertzian pressure at the contact strip is given by 





p y H = p0 1 −



yH a

2 ,

p0 =

2P . πa

(30)

H • Stresses in the halfspace,   x ≥ 0, corresponding to the Hertzian pressure p y H are

⎡

σxHx

σxHy

⎤ H 1 + ya e−ζ ⎦, = p0  ⎣ sinh ζ  H  x −ζ a e = − p0  , sinh ζ

⎡ H σ yy

= − p0  ⎣

⎤ + e−2ζ ⎦, sinh ζ

y H −ζ a e

123

174

Comput Mech (2007) 41:159–174

where z = −y H + ix H = a cosh ζ = a cosh(ξ + iη) is the conformal mapping, which maps the upper halfspace of a complex z-plane to the band  0 ≤ ξ, 0 ≤ η ≤ π in the complex ζ -plane. Coordinates x H , y H are introduced in the Sect. 2.1, see Fig. 2. • Equilibrium displacement field uH that corresponds to the introduced stresses is defined by the formula

−

 H π E1 u H x − iu y 2P(1 + µ1 )



 ∗ e−2ζ = (3 − 4µ1 ) ζ + 2     2y H −ζ e−2ζ , −2 e−2ζ + e + ζ+ 2 a ∗

where the symbol ∗ denotes the operation of complex conjugation. More detailed treatment of Hertz theory can be found in [2] and derivation of the field uH and stresses σ H is given in [4].

References 1. Hertz H (1882) Über die Berührung fester elastischer Körper. J reine und angewandte Mathematik 92:456–171 2. Johnson KL (1985) Contact mechanics. CUP, Cambridge

123

3. Huber MT, Fuchs S (1914) Spannungsverteilung bei der Berührung zweier elastischer Zylinder. Phys Z 15:298 4. Poritsky H, Schenectady NY (1950) Stresses and deflections of cylindrical bodies in contact with application to contact of gears and of locomotive wheels. Trans ASME Ser E J Appl Mech 17:191 5. Timoshenko S, Goodier JN (1951) Theory of elasticity, 3rd edn. McGraw-Hill, USA 6. Zhong Z-H (1993) Finite element procedures for contact-impact problems. Oxford University Press, Oxford 7. Wriggers P (2006) Computational contact mechanics. Springer, Berlin 8. Kikuchi N, Oden JT (1988) Contact problems in elasticity: a study of variational inequalities and finite element methods. SIAM Publication, Philadelphia 9. Eterovic AL, Bathe K-J (1991) On the treatment of inequality constraints arising from contact conditions in finite element analysis. Comput Struct 40(2):203–209 10. Simo JC, Wriggers P, Taylor RL (1985) A perturbed Lagrangian formulation for the finite element solution of contact problems. Comput Meth Appl Mech Eng 50:163–180 11. Zavarise G, Wriggers P (1998) A segment-to-segment contact strategy. Math Comput Modell 28:497–515 12. El-Abbasi N, Bathe K-J (2001) Stability and patch test performance of contact discretizations and a new solution algorithm. Comput Struct 79:1473–1486 13. Wohlmuth B, Krause R (2004) Monotone methods on non-matching grids for non linear contact problems. SISC 25:324–347 14. Kompiš V, Fraštia Lˇ (1997) Polynomial representation of hybrid finite elements. Comput Assist Mech Eng Sci 4:521–532 15. Jiroušek J (1978) Basis for development of large finite elements locally satisfying all field equations. Comput Meth Appl Mech Eng 14:65–92 16. Bathe K-J (1996) Finite element procedures. Prentice-Hall, Englewood Cliffs 17. Lanczos C (1970) The variational principles of mechanics, 4th edn. Dover, New York