-_ .-... “lax“. ~
A“ .. {:fi‘r&

"~44
1:321??? !§\1"1‘,‘ :
,1,',,.1’,f:a1;1,1€,
l‘

{Half r
3
“gr", 1 j
1' ”1"

\n 7'
:‘fryfi

"Amy‘rr'n

"Q? A.
1??!»

K .1
4';
31'“,
'.:‘.’.i’
"‘1
1

 

 

 

 

 

12%?) 59‘“

MICHIGAN STATE U VERSITY LIBRARIES

1" 3 1293 00600 0024
mat-am“ *1

Michigan State4
University

NI

 

 

 

 

 

 

 

 

 

 

This is to certify that the

thesis entitled
An Automated Brocedure for Contact
Problems With and Without Frictional Effects
Utilizing the Boundary Integral Method

presented by

Glen Curtis Bennett II

has been accepted towards fulfillment
of the requirements for

Master's degree in Mechanics

 

 

4110301002- A (Lite no

T

Major pro essor

 

Date M (

0-7639 MS U is an Affirmative Action/Equal Opportunity Institution

—~

PLACE. IN RETURN BOX to remove this checkout from your record.
TO AVOID FINES return on or before date due.

 

 

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

MSU Is An Affirmative Action/Equal Opportunity institution

I "H “-1-! a!”

AN AUTOMATED PROCEDURE FOR CONTACT
PROBLEMS WITH AND WITHOUT FRICTIONAL EFFECTS
UTILIZING THE BOUNDARY INTEGRAL METHOD

By

Glen Curtis Bennett II

A THESIS

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

MASTER OF SCIENCE

Department of Metallurgy. Mechanics and Materials Science

1988

5454404

ABSTRACT

AN AUTOMATED PROCEDURE FOR CONTACT
PROBLEMS WITH AND WITHOUT FRICTIONAL EFFECTS
UTILIZING THE BOUNDARY INTEGRAL METHOD

By

Glen Curtis Bennett II

An automated incremental formulation for the solution of contact problems with
and without frictional effects is presented. The formulation in this paper is
restricted to isotropic, linear elastic materials. The formulation uses a static
procedure to solve problems in the realm of small displacement theory but, with
modification, the incremental algorithms can be expanded to handle problems
within the area of large displacements. The numerical procedure employs a
direct Boundary Integral formulation utilizing linear displacement interpolation
functions and constant traction interpolation functions. This "mixed" formulation

allows the modeling of traction discontinuities without involving supplemental
equafions.

 

ACKNOWLEDGMENTS

The author wishes to express his sincere appreciation for the help and support of
his advisor, Dr. Nicholas J. Altiero, whose persistence during the course of alive

year period allowed me to finish this thesis.

The support of this work by Garrett Turbine Engine Company, Phoenix, Arizona

is most gratefully acknowledged.
The author also wishes to express his appreciation to his mother and father for
all their support throughout his college career and to the rest of his family for their

support and motivation that allowed him to complete his graduate course work.

Finally the author wishes to say thank you to Jim and Jean Eddy for the use of

their computer equipment and knowledge which allowed him to type his thesis.

 

TABLE OF CONTENTS

PAGE

LIST OF TABLES ............................................................................................... vi

LIST OF FIGURES ............................................................................................ vii

CHAPTER I INTRODUCTION ..................................................................... 1

CHAPTER II BOUNDARY ELEMENT THEORY ......................................... 6
”.1 DERIVATION OF THE BOUNDARY INTEGRAL

EQUATIONS ................................................................. 6

”.2 BOUNDARY ELEMENT FORMULATION ................. 12

“.3 NUMERICAL PROCEDURE ....................................... 14

”.4 COMPUTER IMPLEMENTATION .............................. 19

CHAPTER III APPLICATION OF BOUNDARY ELEMENT METHOD TO

CONTACT PROBLEMS ........................................................ 24
”L1 DESCRIPTION OF THE CONTACT PROBLEM ........ 24
II|.2 CONTACT WITH AND WITHOUT FRICTION ............ 37
”L3 THE BOUNDARY ELEMENT EQUATIONS APPLIED
TO PROBLEMS INVOLVING CONTACT ................... 47
“L4 COMPUTER IMPLEMENTATION .............................. 51

PA

CHAPTER IV RESULTS .............................................................................. 60

IV.1 EXAMPLES ................................................................. 60

|V.2 DISCUSSION OF RESULTS ...................................... 81
APPENDIX A INFLUENCE FUNCTIONS .................................................... 83
APPENDIX B CALCULATION OF THE FREE TERM .................................. 84
APPENDIX C CALCULATION OF THE MEAN NORMAL ............................ 86
APPENDIX D ORDER AND CONTINUITY OF TRACTIONS AND

DISPLACEMENTS AT A BOUNDARY NODE ....................... 89
APPENDIX E SINGULAR INTEGRALS ....................................................... 91
BIBLIOGRAPHY ............................................................................................... 92

E

 

IV.2

|V.3

IV.4

LIST OF TABLES

PA
Model description for example 1 ....................................................... 67
Model description for example 2 ....................................................... 68

Contact area versus applied load for discrete nodal
points of example 1 .......................................................................... 69

Contact area versus applied load for discrete nodal
points of example 2 .......................................................................... 70

vi

Fl RE

”.1

”.2

”.3

“.4

“.5

"L1

|II.2

III.3

III.4

|II.5

|II.6

III.7

III.8

I|I.9

LIST OF FIGURES

PA
The integration path around the point of singularity ......................... 1O
Boundary-Value problem in plane elasticity ..................................... 12

Complex geometry broken into segments (boundary
elements) .......................................................................................... 15

The actual corner (left) is modeled using two corner nodes

(right) with typical BEM formulations ................................................ 16
Partial diagram of the boundary showing the

boundary elements with singular considerations .............................. 18
Simple one degree of freedom system [21] ...................................... 25

Load-displacement diagram of the problem shown in

Figure NH [21] ................................................................................. 26
Enlarged view of the contact zone .................................................... 27
Contact problem to be solved ........................................................... 30
Contact problem before load step 1 is applied ................................. 31
Contact problem before load step 2 is applied ................................. 32
Contact problem before load step k is applied ................................. 33
Two elastic bodies in contact ............................................................ 37
Enlarged diagram of the contact area .............................................. 38

vii

“L10

”L11

IV.1

IV.2

IV.3

IV.4

IV.5

IV.6

IV.7

|V.8

IV.9

IV.10

IV.11

IV.12

IV.13

IV.6

C.1

C.2

E1

The vector V,ab ................................................... 41

 

The boundary divided into discrete elements .................................. 46
Two parallel cylinders, with radii of 1.0, in contact ........................... 60

A cylinder of radius 1.0, in contact with an elastic
foundation ......................................................................................... 61

Pressure distribution along the lenght of 2 cylinders
in contact .......................................................................................... 65

Pressure distribution of 2 cylinders in contact shown in 2

dimensions ....................................................................................... 66
Model 1 example 1 (not to scale) ..................................................... 71
Model 2 example 1 (not to scale) ..................................................... 72
Model 3 example 1 (not to scale) ..................................................... 73
Model 4 example 1 (not to scale) ..................................................... 74
Model 1 example 2 (not to scale) ..................................................... 75
Model 2 example 2 (not to scale) ..................................................... 76
Model 3 example 2 (not to scale) ..................................................... 77
Model 4 example 2 (not to scale) ..................................................... 78
Contact area b versus externally applied load P for example 1 ....... 79
Contact area b versus externally applied load P for example 2 ....... 80
Normals for an arbitrary node ........................................................... 86
Definition of the mean normal between nodal pairs a and b ............ 87

Partial diagram of the boundary showing the boundary elements with
singular considerations ..................................................................... 91

viii

INTRODUCTION

The analysis of structural systems in contact is of major importance to engineers.
The contact of objects has two major effects: the transfer of load from one body
to another and in most cases, high stress gradients in the areas of contact.
Systems of this type include bearings, gearteeth, tires, turbine blade roots, and
punches contacting sheet metal during metal forming processes. The analysis of
these contact systems becomes more important as engineers try to simulate

complex structures and how they interact with their surroundings.

The first to examine and solve contact problems was Hertz [1], in 1895. Hertz’s
studies assumed that tangential forces could be neglected and that the only
forces transmitted from one body to another are normal forces. These
assumptions imply that the contact is frictionless. A typical application which can
be simulated by a frictionless model are ball and roller bearing problems.
Based on the work of Hertz, handbooks were developed that gave the contact
area and stress, for a given applied load. One such book was published in 1963
by Lipson and Juvinall [2]. This type of procedure, using simplified analytical
methods to generate tables that could be adapted to similar problems, was used
widely until numerical techniques were developed on digital computers. Other
work in the area of contact problems has occurred since, but existing analytical
solutions are restricted to relatively simple geometries and loading conditions and
have rather limited engineering application.

1

2
To solve problems of practical interest, numerical techniques have been adopted
to solve the problem of contact between elastic bodies. These numerical

techniques allow the analyst to calculate values of interest, 9g.

1. pressure distribution in the contact region,
size of the contact area,

SII'GSSGS near the contact area,

P910!“

final deformed shapes of the bodies in contact.

The most widely used numerical technique for simulation of bodies in contact is
the finite element method (FEM). When using the finite element approach, the
total domain of the body is discretized and, in the contact region, special
interface elements (gap elements) are used. The interface elements are made
up of nodes on each body that are assumed to come into contact. The distance
(gap) between a node on one body and the corresponding node on the second
body is then measured by the interface element. When the two surfaces come
into contact, the contact conditions are imposed [3]. In the region of contact, a
fine mesh is required to predict the values of interest accurately. This can be
done by using a fine mesh in the contact area and then propagating this fine
mesh throughout the interior. This can cause the total system of equations to
become large and computationally expensive to solve. The user may also refine
the mesh in the area of contact and use multipoint constraints, which introduce
master and slave degrees of freedom, to help transition the mesh from fine to
coarse. Other transitioning techniques can also be employed, but in general this
means increased degrees of freedom or the introduction of constraint equations.
These refinement techniques allow for the accurate simulation of the

displacement in the contact region. The only other special consideration that

3
needs to be given to the analysis of contact problems is the nonlinearities that
arise due to the contact between multiple bodies. These nonlinearities are
simulated using incremental techniques that approximate nonlinear behavior in a
piecewise linear manner. This type of procedure is common when finite element

analysis is used to simulate nonlinear behavior.

Fredriksson, in 1976 [4], proposed a procedure using a finite element approach
based on an iterative incremental procedure. His emphasis was on the
development of a slip criterion and an associated slip rule for contact with friction.
Fredriksson also suggested the use of a superelement approach which would
produce a more computationally efficient iterative technique because only the
superelement which contained the contact region would be updated each
iteration. Mahmoud, Salamon and Marks in 1982 [5], developed a numerical
method to provide user convenience, simplify user input, and make the solution
of bodies in contact more computationally efficient by eliminating the iterative
portion of the contact procedure. Using the theory of linear elasticity, they
automated the finite element method to obtain the extent of contact versus the
load increment. Their program dealt with advancing contact problems without
friction. Maxurkiewicz and Ostachowicz, in 1983 [6], used spring-like interface
elements to predict rigid contact, sliding, and rigid body motion within the contact
region. Torstenfelt wrote papers in 1983 and 1984 [7,8] based on a finite
element formulation which utilized only stiffness equations, i.e. no gap elements.
His procedure treated the tangential forces as known quantities, and calculated
them from the previous load step. His second paper used an algorithm that

scaled the applied load to introduce only one change of contact status in each

4
load step. This allows the numerical procedure to be linear between discrete
contact points. In 1984, Rahman, Cook, and Wilkinson [9] developed a
technique without special interface elements that could be easily extended to
multiple contact surface problems. Their program used multiple iteration
requirements in the contact region to control the convergence criterion. Many
other papers have been written using numerical procedures based on the finite
element method, using special interface elements, flexibility formulations,
stiffness equations only, or combinations of all three. Most apply a load
increment and iterate based upon a convergence criterion. This criterion may be
based upon contact conditions remaining constant for more than one iteration or
equilibrium considerations that sum forces (internal and external) until a
tolerance level is reached. However, with both the boundary and the interior
discretized, as done with a finite element formulation, this iterative process can

become computationally uneconomical.

in the solution of a contact problem, the boundary region in contact is of primary
importance. This would indicate that a numerical solution which involves
discretization of only the boundary might be a preferable approach. A numerical
technique of this type would be the boundary element method (BEM) [10,11].
Work has been done on the application of the BEM to contact problems, notably
papers by Andersson [12,13], which explored the use of constant, linear and
parabolic boundary elements. Andersson’s first paper dealt with two-
dimensional contact problems with frictional forces. His approach used constant
interpolating functions for the displacements and tractions. The analysis
assumes the conditions for adhesion, slip, and the direction of the frictional
forces for each load step. These conditions are then checked after the load step

has been applied. if the assumptions are correct, then the next load increment

5
can be applied. It the assumptions are wrong, new conditions are chosen and
the analysis is repeated. This iterative procedure is followed until the correct
contact conditions are found and the analysis is completed. Andersson’s second
paper explored the use of linear and quadratic boundary elements to solve two-
dimensional problems of bodies in contact with friction. Abdul-Mihsein, Bakr and
Parker [14] stated that the Boundary Element Method is well-suited for analysis
of bodies in contact because the tractions are calculated with the same accuracy
as the displacements. Their analysis of frictionless contact problems used an
automated incremental procedure to solve axisymmetric problems. This
procedure used quadratic elements with an automated incremental technique,
exploiting linear elasticity and small deformation theory, to solve problems

involving bodies in contact.

In this thesis, the contributions of these authors are brought together. The
intention is to assemble the best techniques, simplify, automate, and more
accurately simulate bodies in contact. The numerical form of the boundary
integral formulation presented will be different than those discussed earlier, in
that the order of the interpolating functions will be linear for the displacements
and constant for the fractions. This will allow the simulation of corners and
discontinuous loading on the boundary in a straight—forward manner. This differs
from other formulations which use the same order of interpolation functions for
both the displacements and tractions. The present formulation will be applied to
frictionless and frictional contact problems using an automated incremental
technique. The emphasis will be to handle contact problems within the realm of

small displacement theory and linear elasticity.

CHAPTER II

BOUNDARY ELEMENT THEORY

II.1 DERIVATION OF THE BOUNDARY INTEGRAL EQUATIONS

This chapter contains the formulation of the boundary element method as it
applies to plane boundary-value problems in elastostatics. The first section
contains the derivation of the boundary integral equations from the reciprocal
work theorem. In the second section, the derived integral equations are applied
to plane boundary-value problems of elastostatics. The third section contains
the numerical procedure employed to implement the boundary integral equations
via a general computer program. The last section contains the logic used in the

current computer program to solve plane boundary- value problems.

The reciprocal work theorem as stated by Betti and Rayleigh is as follows:

if an elastic body is subjected to two different systems of external
forces, then the work that would be done by the first system of
external forces acting through the displacements associated with
the second system of forces is equal to the work that would be
done by the second system of forces acting through the

displacements associated with the first system of forces [20].

7
If the two equilibrium states (X,, t,, u,) and (X,*, t,*, u,*) exist in the region R
bounded by B then by the reciprocal work theorem the following equation can be

written:

Bt.*(x)ui(x)d8(r) + IRX,*(_)ui(x)da(x)= IBt,(x)u,*(x)ds(x_) + IRX,(x_)u,*(x)da(x),
(11.1)

where
Xi = components of the body force vector,
1. = components of the boundary traction vector,

ui = components ofthe displacement vector,

ds = a differential length on the boundary B,

da = a differential area in the domain R,

and summation over repeated indices is inferred. Next, let us choose the set of
displacements, tractions, and body forces (X,*, t,*, u,*) to be those associated with

a unit excitation at 23:; in an infinite plane, i.e.

X: = O(X_ 'g)ei(c.o)s
U'* = (UR)(X_,£)9(Q), (”2)

I I j j

I: = (IR),,-(K,§)ej(§) !

where
(tR),(x,Q) = the traction, t,(x), due to a unit force, Rj(§), applied in the
infinite elastic plane,
(uR),j(x_,§) = the displacement, u,(x_),due to a unit force, Rj(§),

applied in the infinite elastic plane,

8

6(x_— g) = dirac delta function representing a disturbance at
x=L
e,(§) = unit direction vector.

The functions (tR),,(x,§) and (uR),,(x,§) are called influence functions and are

given in Appendix A. When Equations (“2) are substituted into Equation (ll.1)

the resulting equation becomes
e,(§)I BURN-(LE )u.(x_)ds(.>s) + e.(§)I R8(2<_-§)U,(x)da(r) =
9,15.) [I B(U 81,180 t,(x)d8(x) + IR(uR).,(x.§)X,(2<_)da(r)]- (”-3)
The second term on the left hand side can be simplified as follows:
e,(£).I n5(x-§)U,(2<.)da(r) = e,(§)U,(§)
provided that g is in the domain R bounded by B. Equating coefficients of e,(§)
from equation (”3), the following equations are obtained:
U,(£) = I B(t8).,(2<_,§)U,(z<.)dS(.r) + .IB(UR),,-(r.£) t,(x)d8(x) +
I R(uRMLQX,(>_<.)da(x)- (”-4)

it can be seen from the list of influence functions in Appendix A that if the

symbols 2; and Q are interchanged the following are true:

9
(untita) = (UM-1&0 = (uRliij

(”5)
(IR), [(Lil) = '(lR), i(l,£) 1

provided that ni = n,(§), where n is a unit normal vector to a plane at L. This

yields
UM) =IB(tR),.(.>s£)U,-(Qd8(£) + .IB(UR).,-(x.£)t,(§)d81§) +

IR(UR),,(L§)X,(Qda(Q. (“-6)

where i = 1,2 and x is in the domain R. Equations (“6) can be used to solve for
the displacements anywhere in the domain R given a complete set of boundary

information and are known as the field equations.

The boundary equation can be developed from Equation (”6) by taking the point
xto the boundary B. If Equation (”6) is applied to a point x on B then special
consideration is needed as C —> x since the integrands become singular. The
influence function (tR),,(x,§) varies as 1/p and the influence function (uR),,(_x,§)
varies as ln(p), where p is the distance between the source point Q and the field

point x, i.e.
p = [(X1 ' C02 + (X2 - €2)2I1/2'

See appendix A for a listing of the influence functions (tR)ji and (UR),,. Although
the integrands containing (uR),,(x,§) are singular, the integrals are definite and

pose no problem. The integrand containing (tR),,(x,§) is also singular and the

 

10
integral is undefined as x —> §. To avoid this singularity, the integration path is
taken on a semicircle of radius 8 around the point of singularity, and the limit of

the integral as e ——> O is calculated (see Figure l|.1).

F

    
 

lnteg ration
Path

Figure ”.1

The integration path around the point of singularity.

Thus,

I B(tR)ji(L£)U,-(£)d8(£) = Hm.-. F(tR),-.(2<_.§)U,(§.ld8(§) +

—.I B(iR)j,(L§)U,-(§)d8(§)i (”-7)

where the second integral on the right is the integral contribution for the entire
boundary except at the point of singularity and is known as a Cauchy principal-

value integral. Equation (”6) can now be written as:

 

11

IR(UR),,(£,§)XJ(§)da(§). (“-8)

where

I3“ = [Imago IF(IR)j,(X,§)Uj(€)dS(C)

The left-side of equation (”8) can be written as:

where 01,,(x) = 511' [3,, and Equation (”8) becomes
081(4)U,(2<_) - IBURl-ILQUJQdSK) = IB(UR).,(L§.) t,(§)ds(§) +

I R(UR),,(2<.,£)X,(§)da(Q, (He)

where x is on B. The derivation of the free term 01,, is given in Appendix B.
Before proceeding, it should be noted that the following influence functions are
simular:

(thirst) = «wrap.

where

(uc),,(x,§) = the displacement, ui (x), due to a unit displacement

discontinuity, c,(§), applied in the infinite elastic plane.

 

12
”.2 BOUNDARY ELEMENT FORMULATION

Consider a plane linear elastic region R, with boundary B, supported and loaded

in some manner, as shown in Figure ”.2.

Figure ".2

Boundary-Value problem in plane elasticity.

Now Equation (”9) can be rewritten as:

“II-Xiujlli + IB(UC),,(L§.)U,-(Qd8(§) = IB(UR),,(&§.)1,(§)dS(Q (II-10)

where the body force term has been set to zero. This is called Somigliana’s
identity, which relates the displacements and fractions at any point P=x on B to
all tractions and displacements on the boundary. Equation (11.10) is equivalent to
Equation (”9) but has a more straightforward physical interpretation for the term

(uc)ij [22].

W

 

13

Equations (ll.10) are valid for all x on B, where

B = the boundary enclosing the region R,

x = field point,

Q = source point,

5 = distance measured along the boundary,

ti = components of the traction vector,

Ui = components of the displacement vector,

(uc) = influence function for a unit displacement discontinuity applied in

an infinite elastic plane,

(UR)ij = influence function for a unit force applied in an infinite elastic
plane,
01'- = the contribution of the integral around the singularity point,

called the free term.

A list of the influence functions (uc)ij and (U R)ij are in Appendix A. For every point
on the boundary, either the traction or displacement is known in a given direction.
Equation (ll.10) is used to solve for all unknown boundary information, thus
giving a complete set of known boundary tractions and displacements. Then at

any point x in R the displacements can be calculated from the equation:

um = I sunrise tltids-I Bluc>.,(x_.t>u,<t>ds, (Itii)

which was developed in the previous section for any point x in the domain R.
From equations (ll.11), equations can be formed that will allow the stresses
anywhere in the domain R to be calculated, provided that a complete set of

tractions and displacements are known on the boundary. This is done by

14
applying the plane stress version of Hooke’s law to Equation (ll.11). Hooke’s law

states:

511: (2G/(1 — v))(u,v1 + vu2_2),

0'22 = (2G/(1 - V))(U2'2 + yum),

012: G(U1'2 + U21), (ll.12)
where G is the shear modulus and v is Poisson's ratio. It Equations (ll.12) are

applied to Equations (ll.11) then the results are:

6,..(1) = IBIGR),k,(LQ EIQdS ' IB(GC),,,(L,£)U,(QdS, (II-13)

where

<5ik = components of the stress tensor,

(GR)..,(L

Q) = stress component, oik(x), due to a unit force, R,(§),
applied in an infinite plane,

(oc),k,(_x_,§,) = stress component, o,k(x), due to a unit displacement
discontinuity, c,(§), applied in an infinite plane.

Equations (ll.11) and (ll.13) relate the boundary tractions and displacements to
the internal displacements and stresses anywhere in the domain R. A list of the

influence functions (so) and (15R)ikj are given in Appendix A.

lkj
||.3 NUMERICAL PROCEDURE

Since Somigliana’s identity can only be solved for simple geometries with simple
traction and displacement boundary conditions, the boundary B is broken into N

segments (see Figure ”3).

15

Figure “.3

Complex geometry broken into segments (boundary elements).

Over each segment, approximate functions to the actual variations of the
displacements and tractions are chosen. The integral over the entire boundary B
then becomes a summation of the N integrals over the individual segments.

Thus Equations (ll.10) become

N
“3118983) + ZIm(UC),,(x“.§)u,(st =

N
ZIkuRtixmt) tlods (“.14)

where x", represents a segment endpoints. The variation for the displacements
and tractions can be constant, linear, parabolic or of higher order forms. The
numerical procedure developed for this thesis uses a linear variation for the

displacement function and a constant variation for the traction function over each

16
segment. The formulation of constant fractions and linear displacements is
emphasized here because it allows the model to have discontinuous tractions
and also allows the modeling of corners without any special considerations. This

differs from the formulations which appear in the literature, where:

although displacements at a corner node can be specified
unambiguously, the specified tractions can only be
represented by considering two corner nodes indefinitely close
to each other (typically 0.005 times the length of the local
boundary element apart) representing the limiting endpoints of

the two surfaces [16].

 

 

 

 

Figure ”.4
The actual corner (left) is modeled using two corner nodes (right) with typical

BEM formulations.

A further discussion on the use of linear displacement and constant traction

functions is given in Appendix D.

17
To simplify the integration over each segment further, shape functions are
introduced to approximate the variation of the displacements over each segment
(boundary element), i.e.
u, = U,"‘"N,(§) + U,'“N2(§).
ij = if"
Nit) = .5111).
N2(§) = -5(1+§).
ds = .5(sm- sm_,)d§ = .5Asmdé,

where

U,"“1 = the displacement at the beginning point of element m,
U,m = the displacement at the end point of element m,

a = local coordinate that varies from -1 to 1,

Ni = shape functions,

Asm = length of element m.

Thus Equations (ll.14) become:
N

2a'.,lr">u,.n + mZAs,[u,m-‘I m(UE),,(L”.§)N,(§)d€ +

u,mI miuc>.,<rn,§>N,<t>da] = gasmtrI m(UR),,-(r”.§)d§,

(ll.15)

where t,m is the value of the traction component on element m. This can be
written in a simpler form as

N N
201'..(x”)u." + Z,[A..mnu.m-1 + annum] = 1/Gz,c..mnl=.m, (ll.16)
Ii 1 = If J 'J l m= 'l l

 

18
where

Anni = AsmI ,(uc>.,<x".t>N.l&>dt,
B”... = As,I ,luc>.,<x".§>N,la>de
0.1"" = GI m(uR>,,-(X“.§)d§,

F,m = Asmtj'".

Note that the integrands of A“.mn and Bi,rnn are singular when m=n or m=n+1. This

is the case forthe element that begins or ends with the node n (see Figure l|.5).

n+1

element m=n+1

element m=n

n—1

Figure “.5
Partial diagram of the boundary showing the boundary elements with singular

considerations.

The singular integrals are evaluated analytically while all other nonsingular
integrals are evaluated by numerical integration using a Guass quadrature

formulation. A list of the solutions for the singular integrals are in Appendix E.

19
”.4 COMPUTER IMPLEMENTATION

The computer program is based on the numerical procedure developed in the
previous section. First the coefficient matrices [U0] and [UR] are formed from the
boundary coordinates and material information. Equation (ll.16) expressed in

matrix form becomes
[UCIG{U} = [uRi{F}-

This system of equations relates nodal displacements to resultant segment
forces. To create a well-posed problem relating nodal displacements to nodal

forces, a transformation is required. If we let the nodal force be given by

FT": .5 [F74 Ff“) ,

then a transformation matrix can be written such that :

{F} = [T] {F}. am)

where
l l O O . O
O l I O . O
[T] = .5 O O l l - O
O O O l - O

O
O
O
O
O
O

 

20

and l is the 2x2 identity matrix,

If both sides of Equation (ll.17) are multiplied by [1]1 then the following equation

isformed:

iTl"lFl = {F}

and

[uc] G{U} = [UR] [T]'1{F}, nus)

where
l -l l -l . I
l | -l l . -l
[114: -l l | —l . I

provided that the number of nodes is odd. The next step in the computer
program is to post multiply the matrix [UR] by the inverse of the transformation

matrix. Now our system of equations is as follows:

[U0] G {U} = [M] {F} . (II-19>

21

where

[M] =[uRliTl'1-

To insure equilibrium on the system of equations, a set of supplemental

equations are augmented to the system. The three supplemental equations are:
1. The summation of forces in the X direction, i.e.
N
2 F,I = 0,
i=1

2. The summation of forces in the X2 direction, ie.

.M Z
n"!

ll

.0

II
—-A

3. The summation of moments about node 1, Le.

(X1F2‘- X2F1‘) = o.

'M2

II
—L

in matrix form, these can be expressed as,
[0] {F} = {0}.

where

[01:0 1 o1~oo 1

1_ 2 2_ 1... 1_ N N_1
0 0 x2 x2 x1 x1 x x2 x1 x1

 

 

22
and the system of equations becomes
.. - _ - _ - F -
[UC] 10]T G{U} [M]
= {F} (”20)
[Oi [O] {K} [O]

L _L_ L-L-

 

 

 

 

 

 

 

 

where {it} is a small perturbation (equilibrium residuals) which should approach
zero as N —> oo. Equation (11.20) contains both unknown displacements and
unknown tractions, so to set up a system of equations that can be solved forthe
unknowns, Equation (”20) must be rearranged. All unknown forces and the
coefficients that multiply them are moved from the right hand side of Equation
(11.20) to the left hand side of Equation (11.20). The corresponding known
displacements and the coefficients that multiply them are moved to the right hand
side of Equation (11.20). Once all unknown quantities are in the left hand side
vector and all known quantities are in the right hand side vector, the two matrices
on the right hand side are multiplied togetherto produce a column matrix of real
numbers. The system of equations now looks like,

[A] {X} = {B} . (”21)

 

where

{X}

unknown tractions, displacements and three 7t, terms,

that are related to the equilibrium equations,

[A] the coefficients of [U0] and modified [UR] that multiply the

unknown matrix {x},

{B}

the resultant of the matrices that were multiplied together.

23

Equation (”21) can now be solved for the matrix {x} of unknown quantities.
When the results for {x} are obtained, a complete set of boundary information is
known. From the complete set of boundary information, displacements and
stresses anywhere in the domain R can be calculated from the numerical form of
equations (ll.11) and (ll.13).

CHAPTER III

APPLICATION OF BOUNDARY ELEMENT METHOD TO CONTACT
PROBLEMS

|lI.1 DESCRIPTION OF THE CONTACT PROBLEM

Two bodies are in contact when forces are transferred from one body to another
through an area on the boundary of each body. Here, the contact area of body A
in the deformed state is denoted CA and the contact area of body B in the
deformed state is denoted CB. Each of these areas is a subset of an area

denoted as the potential contact area for each body, CA and CB

potential potential ’

respectively. These potential contact areas are defined by the analyst before the
simulation is started and should include areas larger than the expected contact
areas. After all of the external load has been applied, the actual contact areas CA
and C8 will be subsets less than or equal to the potential contact areas of the two
bodies. The contact region can also be classified as advancing or receding. An
advancing contact region grows larger as the next load step is applied. A
receding contact region shrinks as the next load step is applied. Thus in an
advancing contact problem the initial contact area is less than the final contact
region, and in a receding contact problem the initial contact area is greater than

the final contact area.

24

25
When two bodies are in contact, nonlinearities are introduced into the system.
The response of two bodies in contact will be nonlinear for various reasons. The
first of these reasons is that the stiffness of the system is a function of the relative
displacement between the two bodies. This concept of nonlinearity can be

explained with the use of a simple example (see Figure ”1.1).

 

 

 

 

 

 

t—t/vy

U
gap

 

 

 

Figure ”1.1

Simple one degree of freedom system [21].

It can be seen that when the displacement becomes greater than the initial gap,

the stiffness of the system will change abruptly (see Figure ”1.2).

 

 

L“

26

LOAD

total

slope = (kil- k2)

gap

 

ugap utotal

DISPLACEMENT

Figure “1.2

Load-displacement diagram of the problem shown in Figure NH [21].

The nonlinearity of two bodies in contact can also be expressed in terms of the

boundary conditions. During the loading of the system, the contact area can

27
expand or contract based on the relative displacement of the bodies in contact.
This will change the force and displacement conditions within the contact region
as the relative displacement (gap) of the two bodies changes. Therefore the
force and displacement conditions in the contact region are a function of the

displacement (see Figure l||.3).

 

 

Figure ”1.3

Enlarged view of the contact zone.
This can be represented by the following equations which hold for a given

increment of load, provided there is no slip:

Ati‘i = 0 gap > 0
At, = 0
-At‘*‘n = At"n gap = O
-Atat = Atbt
(lll.1)
Au‘i'n = AU”n gap = 0
Aua, = AU”t

where the direction of the local coordinate system is defined in Figure ”1.9, and

 

Ala.

Atb.

Atan

Atbn

A13

Atb

a
Aun

b
Aun

AU3

28

components of the incremental traction vector for body A
at node a in the global coordinate system,

components of the incremental traction vector for body B
at node b in the global coordinate system,

normal component of the incremental traction vector for
body A at node a in the local coordinate system normal to
the boundary,

normal component of the incremental traction vector for
body B at node b in a local coordinate system normal to
the boundary,

tangential component of the incremental traction vector
for body A at node a in a local coordinate system tangent
to the boundary,

tangential component of the incremental traction vector
for body B at node b in a local coordinate system tangent
to the boundary,

normal component of the incremental displacement
vector for body A at node a in a local coordinate system
normal to the boundary,

normal component of the incremental displacement
vector for body B at node bin a local coordinate system
normal to the boundary,

tangential component of the incremental displacement
vector for body A at node a in a local coordinate system

tangent to the boundary,

29
Au"t = tangential component of the incremental displacement
vector for body B at node bin a local coordinate system

tangent to the boundary.

Another nonlinear consideration is when friction is present in the contact region.
When frictional forces are introduced into the system, the system becomes
nonconservative and the final equilibrium state is path dependent. For these
reasons the simulation of bodies in contact is a very complex problem and

requires special numerical treatment.

As can be seen from Figure “1.2, a problem in the numerical simulation can occur
as the relative displacement closes the initial gap and a new contact status is
achieved. To simulate this nonlinear environment the process is divided into m
steps. Each of these m steps, if properly chosen, can be approximated by linear
equations. This gives a stepwise linear approximation to a nonlinear analysis. If
the contact, adhesion, and slip areas are constant throughout a load step, then
the problem is numerically linear during that load step between discrete contact
points. This type of incremental technique, when incorporated into numerical
methods, allows the user to solve complex contact problems with and without
friction. Incremental techniques will allow the simulation of nonlinear properties
of the stiffness and boundary conditions discussed earlier, and incremental
formulations will also allow the simulation of any irreversible frictional effects
encountered in the contact region. An incremental technique for the solution of a

contact problem would consist of six steps:

30
1. assemble the required matrices,

apply the external loads for the current load step,

 

solve the system of equations for unknown quantities,

update total displacement and traction quantities,

919.00!"

check to see if the total external loads have been applied
A. YES go to step 6
B. NO repeat steps 1-5,

6. calculate the stress quantities requested.

 

The incremental procedure applied to the contact problem of Figure l|l.4, with the
basic assumptions of small displacements and linear elastic material behavior,

would be as follows.

 

Figure “1.4

Contact problem to be solved.

31

 

The externally applied load P is divided into m increments such that,

P1+ P2 + + P = P . (”I-3)

m

where

Pk = the applied external load during load step k.

LOAD INCREMENT 1

   

Figure ”1.5

Contact problem before load step 1 is applied.

The matrices for body A0 and body B0 are assembled and related
using the initial contact conditions. The load increment P1 is applied
to the system (P1 is determined by the incremental rules developed),

where P13 P. The system of equations is solved for the unknown

 

32
quantities giving the equilibrium state for load step 1. The
incremental displacements (AU), and incremental tractions (At), are
now known everywhere on both bodies A and B. We assume for this
example that the total externally applied loads have not been applied

(P,<P) and another load step will be applied.

LOAD INCREMENT 2

 

Figure “1.6

Contact problem before load step 2 is applied.

The matrices for body A, and body B, are assembled and related
through the contact conditions. The current states A, and B, are the
deformed states at the end of the previous load step. The

coordinates of the deformed states are

(X2), = (X2)o + (A02), (III.4)

33
The next load increment Pals applied to the deformed state at the
end of load step 1. Now the system of equations are solved for all
unknown quantities, giving the equilibrium state for load step 2. The

total displacements and tractions are Updated as follows:

(u,)2 = (U,), + (AU,)2. (|||.5)

KTH LOAD INCREMENT

 

Figure III.7

Contact problem before load step k is applied.

The matrices for body A,_, and body Bk, are assembled and related

through the contact conditions. The current states Ak_1 and B,1 are

34
the equilibrium states at the end of the previous load step k—1. The

coordinates of the deformed states are

(X1)k—1 = (x1)k-2 + (AU1Ik-1

(x2)k_, = (x2),,_2 + (Au2)k, (III.6)

The system of equations are solved forthe kth time to compute all

unknown quantities, and the total quantities are Updated once more,

(ti)k = (1014+ (Ati)k
(u,)k = (U,)k_, + (AU,)k_ (III.7)

This procedure is followed until the mth load step is applied at which
time all external loads have been applied to the system. The total
displacements and tractions are equal to the sum of all the

incremental displacements and tractions,

(t,) = (At,),+ (At,)2+ ---+ (At,)m
and

(u,) = (Au,),+ (Au,)2 + + (Au,) (lll.5)

rn.

Stresses can now be calculated from the total displacements and

tractions.

The incremental algorithm is used to follow the load history as closely as

possible. If the increments are sufficiently small, the total quantities are good

 

35

approximations to the exact values, but if the analyst chooses too many small
increments, the computational requirements are large. An incremental solution
can converge to a wrong answer by selecting a load increment that is to large,
and this is the case when too many contact conditions are allowed to change in
the same load step. Thus, the most important part of the incremental technique
becomes the control of how large of a load step to use. The load step should be
chosen such that it is as large as possible without sacrificing the accuracy of the
analysis. There are two types of approaches to the control of the load step

during the incremental process:

1. The load steps can be chosen prior to the analysis using
previous knowledge. Then after each load step a convergence
check is done. This convergence criterion may be based upon
contact conditions remaining constant or equilibrium
considerations that sum forces and iterate until the force
imbalance is within a tolerance limit. If the convergence test
passes, the next load step is applied and then that state is
checked for convergence. If the convergence criterion is not
satisfied, then the equations are adjusted and solved again. This
iterative technique is performed until the convergence criteria is
satisfied. Then the next load step is applied and the process

repeated until the total load is applied.

2. The load steps can be chosen during the analysis such that
within each load step the analysis is numerically linear between

discrete points. The total load is applied to the structure, and

 

 

 

36
based on the changing contact conditions, a scale factor is
determined. This scale factor is used to scale back the load to a
point where linear equations will simulate the deformation
numerically. This would then constitute a load step. The
remaining external loads are applied and the process repeated

until the total load is applied.

It can be seen that the simple example in Figure ”1.1 can be simulated using two
incremental steps, where within each step the analysis is numerically linear
between discrete points. The difference between the two incremental
approaches can be demonstrated using the simple example shown in Figure
“1.1. Using the first approach the initial load step should be as close to Pgap as
possible. This will keep the number of iterations to a minimum and assure a
good solution. If the load step was chosen such that the applied load P was
greater than P gap there is a possibility of multiple iterations to achieve the proper
convergence conditions. This may also lead to convergence to the wrong
answer. If the load is chosen properly (0 s (P - Pgap) s e, where e is small) then
the approximation will be good. Load steps are then applied until all external
loads have been applied and all convergence considerations satisfied. The
analysis is very dependent on a preconceived idea of how the system will react
under loading. If the analysis is too complex for a good estimation of the load-
displacement behavior, then multiple runs will be required to obtain satisfactory
results.

The second approach would be to apply the total load P and, based on a

total

change of contact status at U=Ugap, calculate a scaling factor. The scaling factor

would be equal to PW/Ptotal and would then be used to scale back the externally

37
applied loads and all displacement and traction quantities. Load step 1 would
then be the applied load P gap that caused the first change of contact status in the
P

contact zone. When the additional load (P was applied, there would be

total - gap)

no change of contact status, the scale factor would be equal to 1.0 and the
analysis would be completed. Thus the analysis would be completed in two
steps without any iterations. This approach is not based on any preconceived
idea of the load displacement behavior by the analyst. With this incremental
approach the largest possible load step is chosen, while still allowing the analysis
to remain numerically linear. Although this is a simple example it shows the
capabilities and benefits provided by an automated incremental load step
technique. The automated incremental load step procedure is the logic behind

the implementation of the analysis program developed in this paper.

III.2 CONTACT WITH AND WITHOUT FRICTION

To simulate bodies in contact, two sets of equations are developed for each of

the bodies in contact (see Figure ”1.8).

‘

Figure ”1.8

Two elastic bodies in contact.

38
The equations are developed for body A and body B and the two sets of
equations are then related through the contact region according to the type of
contact being simulated. All points in contact are transformed to a local
coordinate system, normal and tangent to the boundary, and the contact

conditions are imposed. (see Figure Ill.9)

Body A

 

Figure “1.9

Enlarged diagram of the contact area.

If the contact problem is considered frictionless, then only normal compressive
forces are transmitted from body A to body B. The equations that couple the two

sets of integral equations would be:

-At":1 = Atnb,

At,a _ 0,

At,b = o, (”1.6)
AU 3 = AU b,

for all nodes in contact. If Coulomb type friction is introduced into the contact
region, then normal compressive forces and tangential forces are transferred

from body A to body B. For Coulomb type friction the tangential force At, is

39

related to the normal force Atn at each node,

Att _<_ uAtn (III.7)

The direction of the tangential force Att is in the opposite direction of impending
motion. The contact conditions that couple the two sets of equations, for contact

with Coulomb friction, would be

-Atna = Atnb,
-At,3 = A1,”,
Aun":l = AUnb,
AU,a = AU,b,
(lll.8)
1,8 s utna,
t,b s utnb,

for nodal pairs that are in a nonslip state. Although the contact between two
bodies introduces nonlinear effects into the system discussed earlier, the
problem can be solved using an incremental procedure. The incremental
process is a stepwise numerically linear approximation to the nonlinearities that
arise. If the assumptions of small displacements and linearly elastic material
behavior are made, then an automated incremental procedure can be
implemented. The incremental process is used to calculate the contact area for
a given load step. Within each load step the problem can be considered

numerically linear if the increase in the contact area are small and the contact

40
conditions remain constant throughout the load step. This is achieved by only
allowing the load, during any step, to be large enough to cause the next change
of contact status. With the introduction of a load that causes only one change
of contact status per load step, the problem is numerically linear within that load

step. The three states of contact are as follows:

1. new contact,

2. gapping,

3. change from adhesion to slip.

The first of these conditions to occur controls how large the load increment will
be. During the first load step all nodes in contact are assumed to be in an
adhesion state and the total external load P is applied. Then for any load step k
the total external applied load is

k-1

Pk ___ up _ 2.2,], (”1.9)
l=1

where

= applied load after k increments,

P = total external load,
Pi = previous load increment
Yk = scale factorforthe kth increment.

For each load step k, y, is calculated to find the next change of contact status. 7,

41
is calculated for each change of contact status possible and at each node within
the potential contact zone at which that change of contact status could occur.
The smallest y, calculated will be the scaling factor for the kth load increment.
To calculate y, for the case of new contact, the relative displacement along the

unit direction vector V,“ must be computed (see Figure ”1.10).

 

 

Body A
a
/I b _\
Body B
FIGURE lll.10

The vector V,”-

The displacements along the unit vector V,” for nodes a and b are found as

follows:

Node a

(Viablk (Uia) = (D8)

k!

Node b
(Viab)k (Uib)k = (DbIK’

42

where
(V,‘=‘b)k = unit vector from node a to node b for load step k,
(ufl)k = displacement vector for node a for load step k,
(Uf’)k = displacement vector for node b for load step k,
(Da‘)k = the magnitude of the displacement vector
(uf‘)k along the vector (Vf‘b)k for load step k,
(Db)k = the magnitude of the displacement vector

(U,")k along the vector (V,ab)kfor load step k.

Since the displacements of node a and node b along the vector V,” are known,

 

the relative displacement (Dab)k can be calculated,

(Dab)k = (D"")k - (Db)k. (III.10)

The scaling factor y, for the change of contact status to new contact, during the

kth load increment, can be calculated as follows:

7": (Gab),/(Dab)k, (|||.11)

where

(Dab)k = the relative displacement between
node a and node b for load step k,

the nodal gap between nodes a and b

A
C)
m
0'
v
x
II

at the end of load step k-1.

43

The contact conditions for this set of nodes now becomes

At: = At,b
AU,“l = Aunb,
for frictionless contact and
-At,,al = Atnb,
-At,3 = At,b,
Auntl = AU,b
Au,al = Au,b,

for contact with friction and in a non-slip state.

To calculate y, for the case of a nodal pair already in contact pulling apart and
gapping, the load history of the normal force between the contact nodes must be
followed. To be in contact there must be a compressive normal force between
the nodes of the nodal pair. If the force between bodies goes from compressive
to tensile, then the two nodes will tend to gap and there is a point in this process
where the normal force is zero. The point during the application of the load
where the normal force is zero is the point where the nodal pair will begin to pull
apart or gap. The scale factor y, is calculated such that it will scale the applied
load to the point where the normal force between the contacting nodal pairs is
zero. This is accomplished by first calculating the next tangential force (At, )k. that

occurs when the remaining load
k1

pjg,i=>

44
is applied to the equilibrium state of the k-1 load step. The value (tna),,_,
represents the total compressive force between nodal pairs in contact. The value
(Atna), if negative, represents a tensile force between nodal pairs which indicates
the nodal pairs will gap. The scale factor 7, should be calculated such that the
force between nodal pairs is zero. At this point the condition of nodal contact is
removed. The point where the nodal force between the nodes is zero can be

calculated as follows:

«3),, + (Atna)k = o (”1.11)

(Atna)k = - (tna)k-1

YKIAinalt- = '(tnalk-1

y, = -(ina)k_,/(Atna)k. (I||.12)

which will be a postive number. The contact conditions for this set of nodal pairs
becomes
A$=Au=o
(lll.13)
Atna = At,b = 0.

The change of contact status from adhesion to slip will occur when the total
tangential force 1, becomes greater than the maximum allowable shear force t,’“a",

where

t,max = uAtn. (lll.14)

45
After load step k-1 has been applied and the contact is in a non-slip state the

relationship between the tangential force t,3 and the normal force tna is

(A151),1 s u(At,a),_,. ("1.15)

If during the kth load step a nodal pair is to slip then the following must be true,

(A151)k > u(Atn‘=‘)k
or (”1.16)
1.141."),- > mums), -
The scale factor y, can be calculated to find the percentage of the applied load at
which point
(tf‘)k = u(tna)k. (lll.17)
Equation (|||.17) can be rewritten in the following form
(t,"=‘),_1 + yk(At,a)k. = u[(tna),., + y,(At,,a)k. ]. (III.18)
Equation (“1.18) may now be solved for the scaling factor ykgiving
1,: fits)... - (1.3),.11/[(Atf).-- MAtnalk-l- (”09>

The new contact conditions for the nodal pair are

-At,,“‘ = Atnb,

a_ b
Aun - Aun,

46
At,a= uAtna
At,b= uAtnb,

for contact with friction.

When a nodal pair changes state from adhesion to slip, the two nodes begin to
move apart relative to each other in a tangential direction, the basic assumption
of small displacement is violated. The violation is because the arc length
between discrete nodal points must change in order to accommodate the slip
condition experienced. The boundary of each body is divided into a finite number
of discrete points, and in the potential contact region of both bodies the pair of
contact elements that are opposite and expected to come into contact should

have the same arc length (see Figure Ill.11).

BODY A

(element k)

\q
/ BODY B

Figure "1.11

I'

:> Equal arc lengths

 

 

The boundary divided into discrete elements.

47
It is then assumed that due to small displacement theory, the arc lengths of
opposite elements modeled with equal arc lengths before deformation, will have
approximately the same arc lengths after deformation due to loading. This will
allow successive nodal pairs in later load steps to come into contact. If node q
on one side of element k, is allowed to slide, while node r is in the adhesion state
with its nodal pair, then the arc length of element k must expand or contract thus
violating the assumption of small displacements that allow successive nodal pairs
to come into contact. If all nodal pairs begin to slide then a rigid body motion

occurs and the problem is no longer static.

III.3 THE BOUNDARY ELEMENT EQUATIONS APPLIED TO PROBLEMS
INVOLVING CONTACT

The implementation of the Boundary Integral Equations to solve problems
involving bodies in contact is a straightforward extension of the equations
developed in Chapter II. A set of Boundary Integral Equations are developed for

each of the bodies in contact for every load increment

BodyA
summits) - IA(UC),,(A,§)AU,(§)dS = I ,(uR).,(x.t>At,-lt)ds.

BodyB (lll.20)
oc',,(r)AU,(A) - IB(UC),,(A.§.)AU,(£)dS = IB(uR).,-(A.§)t,(§)ds.

and then related through the contact region.

 

48
Now the Boundary Integral Equations must be written in a form to be used in an
incremental formulation. The incremental form of the Boundary Integral Equation

is as follows:

oc',,(x)Au,(x) - IB(UC),,(X,§)AU,(§)ds = IB(UR),,(X,§)A1,(§)ds (Ill.21)

Then the unknown tractions At,and displacements AUi are determined for each load
step. The total tractions and displacements are equal to the sum of all the

incremental tractions and displacements,

m

t= 2, (At).

and (Ill.22)

m

ui = 2 (AU,),,.

The incremental procedure developed earlier in this chapter is based on
introducing only one change of status in each increment. Then for each load
step the simulation remains numerically linear between discrete points and the
incremental Boundary Integral Equations (Ill.21) can be used. The algorithm for
calculating the scaling factor y, for the kth load step is used to find the exact load
level at which the the next change in contact status will occur. The next load
step is applied and again the incremental Boundary Integral Equations (Ill.21) are
used to solve for the unknown displacements and tractions. This procedure is
followed until the mth load step is applied and ym = 1.0. This means that all

externally applied loads have been applied and the analysis is completed.

Since Somigliana’s identity can only be solved for some geometries with simple

49
traction and displacement functions, as stated in Chapter II, the boundary is broken
into N segments. Then over each segment an approximate function to the actual
variation of the displacements and fractions is chosen. The integral over the entire
boundary then becomes a summation of N integrals over each individual segment.
The numerical form of the boundary equation for an incremental procedure then

becomes,

01' lj(x"),AU( )N;I() (U)c,,.( (,)x"§AU,d(t;)s
N
ZxI:(UR),,(x",§)At,(C)ds. (”1.23)

where Equation (”1.23) is developed for each body. The numerical procedure
developed here uses a linear variation for the displacement function and a
constant variation for the traction function over each segment. To simplify the
integration over each segment further, shape functions are used to approximate

the variation of the displacements over each segment (boundary element),

1,. - t,"‘,
NJE) =.5(1-&).
N2(§)=.5(1+§), (Ill.24)
ds 5(s s )dé = 5Asmd§,
where
U,""1 = the displacement at the beginning of element n,
U,“ = the displacement at the end of elementn,
g = local coordinate that varies from-110 1,
Ni = shape functions,

As = length of element m.

50

Substituting Equations (”124) into Equation (“123) gives,

N
201',,(X")AU," + zAsmleufl'lIm(Uc),,(x”,§)N,(§)d§ +

N
Au,mIm(uc),,(xn,§)N2(&,)d§] = gAsmAtijm(uR),,(x",§)d§ (Ill.25)

Equation (Ill.25) is developed for each body and the supplemental contact

equations are used to relate the two systems of equations. This is the same

 

equation developed in Chapter II, but now it is applied to each of the two bodies

the and two sets of equations are coupled through the contact conditions.

For a plane problem using a Boundary Integral formulation, the number of
degrees of freedom per node is four. There are two tractions and two
displacement terms for each node. For every node except the nodes in contact,
two degrees of freedom are known, thus giving two unknowns and two equations
per node. In the contact region the displacements and tractions are expressed in
terms of a local coordinate system, which is normal and tangent to the boundary.
For contact without friction, the number of unknowns per contact nodal pair is six.
For contact with friction, the number of unknowns per contact nodal pair is eight,
while Somigliana‘s Identity provides only four equations, two for each body. This
means that for the case of contact with friction there are four more unknown
terms per nodal pair in the contact region. Therefore the equations generated
using Somigliana’s identity must be supplemented in orderto have the same
number of equations as unknowns. The number of unknown fractions and

displacements for the two bodies in contact, body A and body 8, would be,

51

unknowns = 2(non-contact nodes on body A and B) + 4(contact nodal pairs).

The number of equations that can be generated are,

equations = 2o[(nodes on A) + (nodes on B)],

for contact with friction. To supplement these equations and relate the equations

of body A to the equations of body B we use the contact conditions,

-At,,11 = Atnb,

-At,a = At,b,

Au,a = Au,b,
and (HL26)

AU,a = Aunb,
for the case of contact with friction. These contact conditions supplement four
equations for every nodal pair in contact. Thus there are the same number of
equations as unknowns, and the equations of body A are related to the equations
of body B through the supplemental contact equations. This set of equations can

now be solved for all unknown displacements and tractions.

III.4 COMPUTER IMPLEMENTATION

The flow of the incremental contact program is similar to the program described

in Chapter II. First the coefficient matrices [U0] and [U R] are formed for each

 

52
body A and B, from the boundary coordinates and material information. Each set

of equations can be expressed as follows:

Body A

[UCIAGAIAUL = [U RMAF},

BodyB (Ill.27)
[UCIAGBIAUh = [URIB{AF}B-

This set of equations can be written in the following matrix form:

 

 

 

 

 

 

 

 

[uc]A 0 GAAUA [UR]A 0 AFA
0 [uc]B CbAUB 0 [UR]B AFB _

Each system of equations relates nodal displacements to resultant segment
forces. To create a well-posed problem relating nodal displacements to nodal

forces, a transformation is required. If we let

in

m +
AFi = .5 (AF, + AF:n I, (“1.28)
then a transformation matrix can be written such that

{AF} = [T]{AF}, (Ill.29)

53

where
l | O O o O
O l l O o O
[f] = 5 O O l I 0 O
O O O l . O
I O O O - l ,

and l is a 2x2 identity matrix

[33]

Equation (”1.29) can now be rearranged by multiplying both sides by [1']1 giving

the following:

-1

[T] {AF} = {AF}
This transformation is performed on both systems of equations.
Body A

iUC].G.{AU}. = iuRl. [T13 {AF}.
Body B

[UC].G.{AU}. = luRl. lime.

54

where

provided that the number of nodes is odd. The next step in the computer
program is to post multiply the matrix [UR] by the inverse of the transformation

matrix. Now the system of equations is as follows:

 

 

 

 

 

 

 

 

" '- - - - - - 1
[uc]A 0 GAAUA [M]A 0 AFA
_ 0 [UC]B_ _GBAUB _ _0 [M]B_ _AFB_ ,
(”1.30)

where
[M] = iuRiiTif‘

This equation relates the displacements and tractions at any point P with all other
tractions and displacements everywhere on the boundary. This is true for both
sets of equations involving body A and body B. To enforce equilibrium on both
systems of equations, a supplemental set of equations are augmented that will

insure equilibrium. The three supplemental equations for each body are,

 

 

55

 

1. The summation of the forces in the X direction;

N
; AF,i =0,

2. The summation of the forces in the X2 direction;

N
Z. AFZI = o,

3. The summation of the moment about node 1;

N
2 (x1 AF2‘ - X2AF,i) = o.

 

In matrix form this can be expressed as,

[0] {AF} = {0}.

 

 

where
1 O 1 O o o o 1 O
[Q] = 0 1 O 1 o . . 0 1
Im0 0 X21“X22 X12‘X 1 o o o X21‘X2N X1N_X1‘l -

56

The system of equations now becomes,

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

"o o ‘ ".A ‘ " [01A ‘ " "
101T luc] o 1} [M] o F
A A UA A {}A
T
IUCIB -[Q]B {UIB 0 IMIB {F}B
o
o 9. [o
_ ° _ _ B _ L 13 _ _ _

 

(lll.31)

The equations for body A and body B are unrelated at this point in the program.
Therefore this system of equations needs to be coupled through the contact
region as described earlier in this chapter. In order to use the contact conditions
to couple the system of equations, the nodal pairs coming into contact must be
transformed into a local coordinate system. The local coordinate system is
defined by the mean normal between the nodes in contact on body A and body
B. The calculation of this mean normal is described in Appendix C. A
transformation matrix [L] is used to transform the coefficients of the nodes in
contact for both bodies to a local coordinate system, normal and tangent to the
boundary. The transformation matrix [L] is of the following form:

F cos 0,, sin 0,, I.

[L] =

 

 

- Sin 6m cos 0,,

After the transformation is made the system of Equations (lll.31) become:

 

 

 

 

 

 

 

 

 

 

(111.32)
—0 o _ "1,, [01,, " '
101T iuci {u}' [M]' o F.
A A A A {}A
iucl [Q]; iuié, o [Ml,'3 {F},3
0
o 0 AB [0],,

 

 

 

 

 

 

 

 

 

 

Equation (”1.32) has both unknown displacements and unknown fractions, so to
set up a system of equations that can be solved forthe unknowns, equation
(”132) must be rearranged. Also, Equation (”132) still has two sets of
uncoupled equations, which need to be modified to incorporate the contact
conditions. The contact conditions that relate the two sets of equations are of the
form

U,81 = Unb.
This condition can be inforced by moving the coefficients that multiply U,b in the
left hand side matrix so that they multiply Una. The column that contained the
coefficients is then removed. A similar process is preformed to enforce the other

contact conditions. All unknown forces and the coefficients that multiply them are

moved from the right hand side of equation (“1.32) to the left hand side of

 

58
equation (lll.32). The corresponding known displacements and the coefficients
that multiply them are moved to the right hand side of equation (lll.32). Once all
unknown quantities are in the left hand side vector and all known quantities are in
the right hand side vector, the two matrices on the right hand side are multiplied
together to produce a column matrix of real numbers. The system of equations

now looks like
[A] {X} = {B}, (lll.33)

where
{x} = contains unknown tractions, displacements and six terms
related to the equilibrium equations,
[A]: the coefficients of [U0] and modified [UR] that multiply the
unknown matrix {x},
{B} = the resultant matrix of matrices that were multiplied

together.

At this point the nodal gap between perspective nodal pairs is calculated. The
gap is calculated as the distance between potential nodal pairs along a vector
between nodal pairs in the direction from body A to body B. The set of nodal
gaps are stored for later use in determining the scaling factor for the current load

step.

Equation (”1.33) can now be solved for the column matrix {x} of unknown
quantities. When the results for {x} are obtained, a complete set of incremental
boundary information is known. The scaling factor 7,, can now be calculated

using the procedure described earlier in this chapter. The purpose of the scaling

59
factor is to find the percentage of the applied load step that will cause the next
change in the contact status. Thus providing that the equations will remain
numerically linear throughout the load step. The total displacements and
tractions are Updated by adding the incremental quantities to the total quantities
from the previous load step. If y, = 1.0 then the total external load has been

applied and the analysis is complete.

(tn)k = (tn)k-1 + A(tn)k
(tt)k = (it).-. + A(tt)k
(lll.34)
(un)k = (un)k-1 + A(un)k

(Ut)k = (ut)k-1 + A<Ut)k

If y,C < 1.0 then the coordinates must be updated to reflect the deformation from
the last load step and the complete process is repeated. Once the total
displacements and fractions have been Updated for the last time, to give the final
quantities on the boundary, the displacements and stresses anywhere in the
domain R of each body can be calculated from the numerical form of Equations

(ll.11) and (ll.13) for each body A and B.

CHAPTER IV

RESULTS

IV.1 EXAMPLES

The procedures developed in the previous chapters will be applied to example
problems with analytical solutions, so that a comparison can be established. The
first example is the solution of two parallel cylinders in contact. Both cylinders
have a radius of 1.0 and are infinite in length (see figure IV.1).

P

Figure IV.1

Two parallel cylinders, with radii of 1.0, in contact.

The examples problems solved for were made of the same material and had the

following material properties

60

61
E = 2.5

= .25
The applied load versus contact area is compared to the analytical solution found
in the Hangbggk of Stress and Strength [2]. Multiple runs were made with
different degrees of discretization to test the effect of mesh density on the a
numerical results. The boundary of the two cylinders is all that needs to be
modeled, when using a boundary integral formulation, as opposed to a finite
element solution which would discretize the total domain. It can be seen from
this simple example that when the boundary integral method is used, the data

preparation is decreased dramatically.

The second example is for the solution of a cylinder on an elastic foundation.
The radius of the cylinder is 1.0, and the width of the foundation is 4.0 (see figure

IV.2).

 

 

 

 

 

 

—1I_\
L

Figure IV.2

A cylinder of radius 1.0, in contact with an elastic foundation.

 

62
Both bodies where made of the same material having the following properties
E :25
v = .25

 

The load versus the contact area is compared to the analytical solution found in

the Handbogk of Stress and Strength [2]. Again multiple runs are made to see

the accuracy at various discretization levels.

As seen in Fiqure IV.3 the pressure distribution is parabolic and constant along
the length of the cylinders. Although Figure IV.3 is fortwo cylinders in contact,

the same type of pressure distribution occurs for contact involving a cylinder on

 

an elastic foundation. Both examples can be simulated using a plane model and
assuming a unit thickness. The applied load P verses the contact area bis used
as the comparison between the analytical results and the numerical results from
the development program. The analytical results as presented in the Handbook

of stress and strength [2] for two cylinders in contact is

b = 1.13 [PA(a(1/R, + 1/R2))]"2

and for a cylinder on an elastic foundation the equation is

b = 1.13[PRA/a]"2.

For both examples the model is for a unit thickness

63

and the radii are

R=R,=R2=1.0

The variable A is evaluated using the following formulation

A = (1-v,2)/E, + (1-v22)/E

2

where
v, = v2: .25
and
E, = E2: 2 5
This gives a value of
A = .75

for both contact problems. The final equations for the contact area b versus the

applied load are

b = .69198085P"2

for contact between two cylinders and

64
b = .9786087P"2

for the contact of a cylinder on an elastic foundation. For a definition of a, b and

R see Figures IV.3 and IV.4.

Table IV.1 and Table IV.2 contain a summary of the discription for both
examples. Table IV.3 and Table IV.4 contain the analytical and numerical results
forthe contact area b verses the applied load P for example 1 and example 2.
Figures IV.5 through IV.8 show all 4 models for example 1. Figures IV.9 through
IV.12 show all 4 models for example 2. The numerical results are also plotted
against the analytical curve for example 1 in Figure N13 and for example 2 in

Figure IV.14.

 

 

.8980 E 285:6 N 6 50:2 9: 92w 8:32:26 9:89.“.
«.2 9:9”.

65

 

 

 

 

66

.wco_w:oE_U N E 56% SmEoo E maoccio 93 B cozsnEwB 9835
v.
N

>

 

_ 939“.
r

 

 

 

__ fl..—v

68

 

: m m H
mm mm me m;
mm 5 _.N or

9:8
89:8
.6 .8822

N .68
35:56
.6 52:52

esooa
wEoEo_o
6 28:52

 

v _oUo_>_ m _ooo_>_ N _oUOE F _ono_>_

.m 29:96 6.. 83:88 Eco—2

N.>_ 058.

 

69

 

romeo. cameo. eomeo. mambo. smeeo. amass.
mmmmo. mmmmo. ommmo. somma.
moomo. mmomo. mmmmo. secs.
mmaoo. mmsoo. eaooo. meooo. mmeao. mammo.
mesoo. sauce. emmoo. Nwmoo.
momoo. Nomoo. ocmoo. mmaeo.
44000. 44000. cmooo. emooo. _mooo. emomo.
a a a a a
a _oeozv A... .282 a _oeozv : 68% 98383 a
two; 83 two; nmou boot. moi

8:92 8__&< Ba? .0932 8:84., sage

.F anmxo lo 958 .260: 2986 22 932 83% mamao> moam 63:00
m.>_ 228.

70

 

 

cameo. cmmeo. moose. «also. memmo. smear.
OFoFO. momeo. meoeo. somme.
momoo. mmmoo. eeeao. oseoe.
mFooo. Noooo. emeoo. oaeoo. mmeoo. mammo.
Nmmoo. ommoo. Naeoo. Names.
omeoo. omaoo. mmaoo. melee.
Nmooo. mmooo. mmooo. mmooo. meooo. emomo.
a o a o a
$685: 86.8% aoeozvcaoozv findings 0.
poo; poo; poo; poo; vmou moi
eo__&< 8:93 Ea? uo__&< e232 coccoo

.m oEmeo lo Egon _mno: 9986 6.. 32 8:an mamao> moan 83:00

v.>_ oEmP

71

Figure IV.5
Model 1 example 1 (not to scale).

 

72

Figure IV.6
Model 2 example 1 (not to scale).

73

7 equally spaced
nodes on both ""

/ ‘:“ \
sides of center and a!

on each body

 

Figure IV.7
Model 3 example 1 (not to scale).

 

7 equally spaced

nodes on both , ”’— “3?:
____ 3.:

sides of center and \ :.: /

on each body

Figure lV.8
Model 4 example 1 (not to scale).

75

 

 

0
D
‘f
D
I
I
I
I

 

Figure IV.9
Model 1 example 2 (not to scale).

 

76

1?

 

 

Figure IV.10
Model 2 example 2 (not to scale).

 

  
  
  

7 equally spaced
nodes on both
sides of center and
on each body

 

1p——o———0—-—o

77

 

 

Figure IV.11
Model 3 example 2 (not to scale).

  
  
 

7 equally spaced
nodes on both
sides of center and

on each body
1—7’”

 

78

 

 

 

Figure IV.12
Model 4 example 2 (not to scale).

79

CONTACT BETWEEN TWO CYLINDERS

 

 

.071
.06- O
05‘ I] o
. j .0
LOAD .04. O .
P . o o analytical
.03. 0; El model 1
. o A model 2
.02- 09 0 model 3
. 00 O + model 4
.01 ~ 0
. as 0 CO at
O ....... G ..... Q O
._01 . , . . ffi . . - , , ,

 

o .02 .04 .06 .08 .1 .1'2 '14 .136 '18
CONTACT AREA 3

Figure IV.13
Contact area b versus externally applied load P for example 1.

80

CONTACT BETWEEN A CYLINDER AND
_035. AN ELASTIC FOUNDATION

 

 

.03- O
‘ o
LOAD 025:
P O
02 0%
1 $0
.015- o<> O analytical
« OO I'_'I model 1
.01- o 03 Amodel2
. Ce 0 model 3
.005 C8 0 n + model 4
. O
0 ....... Q ..... Q 0 8
-.005 . . - . - . , .

 

0 .02 .04 .06 .08 ° .1 ' .1'2 '.14 3 .1'6 ‘18
CONTACT AREA B

Figure IV.14
Contact area b bersus externally applied load P for example 2.

81
V2 DISCUSSION OF RESULTS

From the two examples problems run, several things may be observed. The
accuracy of results at any node, especially those near the contact region, are
dependent on the element lengths in the vicinity of that particular node. This
effect can be seen by looking at the results for the last contact point (b=.14658)
for example 1. In model 1 the element lengths neighboring the last node are
.02094 and .369618,respectively. The load P required to extend the contact area
to include this node is .05398 which is greater than the analytically predicted
result of .04487. In model 2 a node is added to reduce the larger element length.
The element lengths neighboring the last node are now .02094 and .11512. The
load P required to bring the last node into contact is .04504, which is a much
better answer. In fact if another node was added, the load P required to bring the
last node into contact would be reduced again. What can be seen from this is
that the accuracy of results is highly sensitive to the mesh lenths surrounding that

node.

The results for the first node in most cases considered deviated from the
analytical values more than the values calculated for the other nodes in the
contact region. This may be the result of having two large of a mesh size around
the initial contact point. A better formulation may be to use a higher order
interpolation function for both the tractions and displacements. This type of a
formulation would allow the analyst to use a coarser mesh and still obtain

accurate results

As the mesh was refined the results tended toward the analytical values. This

type of trend would indicate that numerical values are converging as the mesh is

82
refinded. This property of convergence is very important as a test for the validity

of numerical programs.

The boundary integral formulation produces ill-conditioned matrices, which cause
numerical problems. The current solution technique used is a Guassian
elimination procedure. Other techniques that are developed for ill-conditioned
matrices should be sought to improve the accuracy of the analysis. Techniques
to monitorthe decomposition phase and report numerical accuracy loss, should
also be implemented so that the analyst will know when numerical inaccuracies
have occurred. This type of numerical problem will account for some of the

discrepancies in the results presented.

In general, the procedure works well, but numerical aspects of the method need

to be modified to provide more accurate results.

 

APPENDIX A

INFLUENCE FUNCTIONS

83
The influence functions for plane stress are as follows [18]:

(UR)ik = [-(3-v)6,,log p + (1+v)q,q,,]/(8rcG)
(UC),, = [2(1+v)(n,q,3-n2q23) + (1-v)n,q, + (3+v)n2q2]/(41rp)

(uc) = [2(1+v) (-n2q,3-n,q23) + (1+3v)n2q, + (3+V)n,q2]/(4112p)

12

(UC) = [2(1+v) (-n2q,3-n,q23) + (3+v)n2q,+ (1+3v)n,q2]/(4it:p)

010).. = [2(1+V) (-n.q.3-n.q.3) + (1-V)n2q2+ 3 + Vln.<l.l/(4vip)
(Git).11 = l-2(1+V)q.3-(1-V)Q.l/(41tp)

(0R) = [2(1+VIq23-(3+V)q.l/(4np)

(on) = [2(1+v)q,3-(1+3v)q,]/(4itp)

(OR) = [2(1+V)q23-(I+3V)q2]/(47tp)

’5
59
ll

.2. [2(1+v)q.3-(3+v)q.l/(4irp)
222 ['20 +V)q23-(1 ‘V)q2]/(4T[p)

(60) = G(1+v) [(1+4q,2-8q,4)n,+2q,q2(1-4q,2)n2]/(21tp2)

a
3’
ll

(cc)... = G(1+v) [(1~8q.2q.2)n.+2q.q.(1-4q.2)n.l/(21rpzl
(60)..1 = G(1+v) [(1 -8q.2q.2)n.+2<rl.q2(1-4q.2)n.l/(2npzl
(cc) = (so)

121

(CC) = (oc)22,

(cc) = G(1+v) [(1+4q22-8q24)n2+2q,q2(1-4q22)n,]/(21tp2)

p = [(X1-Q1)2+(X2-C2)2]°5

q1= (X1-C1I/p q2 = (X1-C2) /p

 

 

APPENDIX B

CALCULATION OF THE FREE TERM

84
The free term 06,, does not have to be calculated by solving the integral around
the point of singularity. A relationship forthe free term can be obtained by
considering what will happen during a rigid body motion. Then from information

already calculated during the analysis the free term is determined.

Recall that
[uc]G{U} = [UR]{F}.

If we now look at this equation during a free body translation of an amount UO in

the x and y directions, the following can be written,

{F} = {0}

and

{U} = {U0}-

The matrix form of the boundary integral equations can then be written as,

[UCl{U.} =10}.

This implies that
N
(UC)2,,_, 2n-1 = -2(UC)2,,-, 2m-1
m=1
n

mi

This means that the free term never has to be calculated directly, and the
coefficient containing the free term is found by adding all the terms that multiple a

displacement in the x or y direction, and taking the negative.

APPENDIX C

CALCULATION OF THE MEAN NORMAL

86

To define the local coordinate system of any set of nodal pairs a “mean normal”

is defined. Each node is bounded by two elements on each side. Each element

may have a different normal to the surface, see figure below.

Now the following is defined:

and

n, = cos 8, n2 = sin 8,

n,,= cos 0,, n2, = sm 8,

n,s = cos 6,, r12s = sm 0,,

element n+1 \,

element n

Figure 0.1
Normals for an arbitrary node.

87
where
= normal to element n,

[1,
n,s = normal to element n+1,

and
|.5( n, + D.“ = .5[( cost)f + c036,)2 + (sine, + sin0,)2]-5
= -5[2 + 2(cose,cos0s + sinef sin0,)]-5

cost) = (cost)f + coses)/ [2 + 2(cos6,cost)s + sinef sin08)]-5

sine = (sinef + sin65)/ [2 + 2(cost),coses + sinef sin6,)]-5

Next consider a set of nodal pairs that are in contact, as seen below.

 

Figure C.2

Definition of the mean normal between nodal pairs a and b.

 

88

The “mean normal” from node a to node b, can now be defined as,

Em: -5( fl... ' flb)/I'5( fl; flsll.
where
|.5( 3. - n,)| = .5[( 0050,, - 0066,)2 + (sine, - sin0,)2]-5
= -5[2 - 2(coseacoseb + sine, sin0,)]°5
and

c030m = ( cost),1 - coseb)/ [2 - 2(coseacoseb + sine, sin0,,)]-5

sinem = (sinea - sineb) / [2 - 2(cos0,,cos()b + sine, sin8,,)]-5

This defines the local coordinate normal and tangent to the surface of each body.

 

APPENDIX D

ORDER AND CONTINUITY OF TRACTIONS AND DISPLACEMENTS AT A
BOUNDARY NODE

 

 

89
When developing the numerical formulation of the Boundary Element Method a
choice of order and continuity of the traction and displacement functions must be
made. What is meant here by choice of order is not whether the tractions and/or
displacements will be simulated by constant, linear, or higher order functions, but
will the tractions and displacement functions have the same order. The two

choices are

1. both the tractions and displacements have the same order interpolation

functions (constant, linear, ...),

2. the traction and displacement functions have a different order (constant

tractions and linear displacements).

First consider the case where the interpolation functions are of the same order
(both constant, linear, quadratic, ...). To have a square system of equations
(square meaning the same number of unknowns as equations) the continuity
conditions must be the same at each node. That is, both the traction and dis-
placement components are either continuous or discontinuous. The problem
with a discontinuous formulation is that a discontinuity of displacement at a node
means that a node would have to split because the displacement on each side of
that node is multi-valued. This does not simulate a continuous medium well.

The more common formulation found in most industry standard programs uses a
continuous traction and displacement formulation. In industry this is known as an
isoparametric formulation. The problem with this formulation is that at the cor-
ners of a model a node is mUlti-valued. To simulate the corner these programs
place two nodes very close to each other. This takes care of the multi-valued

traction problem but is not a true representation of the actual problem being

90
simulated. A second problem with this type of formulation is that most load
model loading is discontinuous. That is a load is concentrated over a prescribed
area and everywhere else is zero. This poses another problem for this type of

formulation.

Instead of choosing the displacement and traction functions to be of the same
order, one may chose the displacement functions to be an order lower than the
order of the traction functions. In order to have a square system of equations for
this formulation the displacement must be discontinuous and the tractions must
be continuous. This type of a procedure has two problems. As stated above a
discontinuous displacements formulation causes problems simulating a continu-
ous medium. Also, continuous traction formulations have to have special proce-
dures for corners and discontinuous loading. The other choice for order of the
traction and displacement functions is to allowthe displacement functions to be
one order higher than the traction functions. A square system of equations may
be obtained for this formulation with continuous displacements and discontinuous

tractions.

 

APPENDIX E

SINGULAR INTEGRALS

 

 

91

The singular integrals are evaluated analytically while all other nonsingular inte-
grals are evaluated by numerical integration. The integrals become singular
when the element begins or ends with the node n (see Figure E.1).

n+1

element m=n+1

element m=n

n-1

Figure E.1
Partial diagram of the boundary showing the boundary elements with singular
considerations.

The contributions from the singular integrals are as follows:

AHnn: A11n+1n = Bunn: B11n+1 n = 0
A22n ” = A22“1 ” = B22n ” = 822"+1 " = 0
A,2” n= -A2,n ” = (1—v)/2it

8,2n+1 n= -BZ,n+1 ” = -(1-v)/21l:

A,2n+1n + 8,2“ n = A2,“ n + B2,“ n = ((1 -v)/2n)log(As,/As,,,)

 

 

 

BIBLIOGRAPHY

92

BIBLIOGRAPHY

Hertz, H., miscellaneous papers -On The Contact Of Elastic
Solids (translated by D. E. Jones) Macmillan, London, (1896).

C. Lipson and R. C. Juvinall, Handbook 91 Strese and Strength,
Macmillan, New York (1963).

ABAQUS, Example Problems Manual, Hibbitt, Karlsson and
Sorensen, lnc., Version 4.5, (1982).

Fredriksson, 8., Finite Element Solution of Surface Nonlinearities
in Structural Mechanics With Special Emphasis to Contact and

Fracture Mechanics Problems, Semputere S Structures, S, 281-
290(1976)

Mahmoud, F. F., Salamon, N. and Marks, W. R., A Direct
Automated Procedure for Frictionless Contact Problems,

lnternetional Jegrnel Fer Ngmerieel Metheds In Engineering, _1_8,
245-257 (1980).

Mazurkieewicz, M., Ostachowicz, W., Theory of Finite Element
Method For Elastic Contact Problems Of Solid Bodies,

Computers & Structures, 1_, 51-59 (1983).

Torstenfelt, 8., Contact Problems With Friction in General
Purpose Finite Element Computer Programs, W
Smarties. 16. 487-493 (1983)-

Torstenfelt, B., An Automatic Increment Technique for Contact

Problems with Friction, WM 19, 393-
400(1984).

Rahman, M. U., Rowlands, R. E., Cook, R. D., Wilkinson, T. L.,
An Iterative Procedure for Finite-Element Stress Analysis of

Frictional Contact Problems, Sempgters 8 Strgctgree, 1_8, 947-
954(1984)

 

10.

11.

12.

13.

14.

15.

16.

17.

18.

19.

20.

 

 

93

Rizzo, F. J., An Integral Equation Approach to Boundary Value
Problems Of Classical Elastostatics, 83-95 (1965).

Rizzo, F. J., The Boundary-Integral Equation Method: A Modern
Computational Procedure in Applied Mechanics, 1-5 (1975).

Andersson, T., The Boundary Element Method Applied To Two-
Dimensional Contact Problems With Friction, Proceedings of the
Third International Seminar, Irvine, California, July 1981, Ed. C.
A. Brebbia, Springer-Verlag, Berlin.

Andersson, T., The Second Generation Boundary Element
Contact Program, Proceedings of the Fourth International
Seminar, Southhampton, England, September 1982, Ed. C. A.
Brebbia, Springer-Verlag, Berlin.

Abdul-Mihsein, M. J., Bakr, A. A. and Parker, A. P., A Boundary
Integral Equation Method For Axisymmetric Elastic Contact

Problems, Computers S Structuree, 2_3, 787-793 (1986).

Kalker, J. J., On Elastic Line Contact, Journal of Applied
Mechanics. 3,1125-1132 (I972).

Banerjee, P. K. and Butterfield, R., Bpundary Element Methods
In Engineering Seienee, McGraw Hill Book Co.,London, (1981 ).

Press, W. H., Flannery, B. P., Teukolsky, S. A. and Vetterling, W.

T., Numerieel Recipes The Art pf Seientific Computing,
Cambridge University Press, Cambridge, (1987).

Garcia, J. J., Coupling of the Finite Element and Boundary
Element Methods forthe Solution of Plane Problems of Linear
Elasticity, MS. Thesis, Department of Metallurgy, Mechanics and
Materials Science, Michigan State University (1984).

Timoshenko, S. P. and Goodier, J. N., Theory of Elasticity,
McGraw Hill Book Co., New York (1953).

Wang, C. Applieg Eleetieity, McGraw-Hill, New York, (I953).

21.

22.

94
Cook, R. D., Concepte and Applieatipne of Finite Element
Analyeie, John Wiley and Sons, New York (1974).

Altiero, N., Gauazza, On a Unified Boundary Integral Equation

Method, Joprnal of Elesticity,fl, 1-9 (1980).

 

 

 

NIV I,_I IBRQRI IES

‘ lllllllllllllllllllllllllllll llllll llf l l

 

312930060