THE APPLICATION OF THE MIXED FINITE ELEMENT METHOD TO THE ELASTIC CONTACT PROBLEM by JORCITO TSENG (B.A. Sc., The University of British Columbia, 1977) A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE STUDIES Department of Civil Engineering We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA January 1980 (c) Jorgito. Tseng, 1980 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the Head of my Department or by his representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Civil Engineering The University of British Columbia 2075 Wesbrook Place Vancouver, Canada V6T 1W5 Jorgito Tseng ii ABSTRACT The finite element method is applied in conjunction with Reissner's mixed variational principle to the investigation of two-dimen-' sional elastic contact problems. The versatility of the mixed principle in incorporating boundary conditions pertinent to the contact problem is demonstrated. Contact conditions are modelled by appropriate manipulations of boundary variables. In cases where the contact region is independent of the applied loading, an iterative procedure is used to establish the sliding and adhering portions in the contact region. Numerical results for displacements and stresses are independently confirmed by the finite element analysis in conjunction with the minimum potential energy principle. In cases where the contact region is a function of the applied loading, or progressive contact, an incremental formulation is employed to describe the discretized contact stages. In the example of a frictionless contact between a long cylinder and a rigid base, good confirmation is obtained from Hertz's analytical solution. Criteria for taking one contact stage to another are also outlined for frictional progressive contact. iii TABLE OF CONTENTS Page ABSTRACT iii TABLE OF CONTENTS i i i LIST OF TABLES v LIST OF FIGURES vi ACKNOWLEDGEMENTS vii! NOMENCLATURE ...:-.=: >* CHAPTER 1. INTRODUCTION 1 1.1 Background 1 1.2 Purpose and Scope 4 2. THEORETICAL BACKGROUND 5 2.1 Reissner Principle 5 2.1.1 Background 5 2.1.2 The Variational Principle 6 2.2 The Finite Element Method 11 2.3 Application of the Finite Element Method to Reissner's Principle 13 3. APPLICATION TO THE CONTACT PROBLEM 16 3.1 Contact Conditions , 16 3.2 The Contact Model 19 iv CHAPTER Page 3.3 Rigid Body Modes 27 3.4 Progressive Contact 30 4. COMPUTER PROGRAM 35 5. SAMPLE PROBLEMS 41 5.1 Sliding Block 41 5.2 Stick-Slip Problem 56 5.3 Frictionless Progressive Contact Problem 73 6. CONCLUSIONS 80 REFERENCES 82 APPENDIX A — The Finite Element Model 85 APPENDIX B — Constraints 94 APPENDIX C — ....................... 99 C.1 Derivation of the Separating Distance between Two Nodes in a Contact Nodepair 99 C.2 Derivation of the Increment Load Factor k' , 103 m+1 V LIST OF TABLES Table Page 5.1 Numerical Results for Problem A: A Frictionless Contact Problem 44 5.2 Numerical Results for Problem B: A Sliding Contact Problem without Normal Stress 47 5.3 Numerical Results for Problem C: A Sliding Con-tact Problem with Normal Stress 53 5.4 Numerical Results for Problem C: Convergence of Displacements with Grid Refinement 54 5.5 Boundary Condition Results for Stick-Slip Contact: Comparison between QST and Mixed Element in the Determination of the Adhering-Sliding Changeover Point with Iteration Convergence 60 5.6 Boundary Condition Results for Stick-Slip Con-tact: Comparison between QST and Mixed Element in the Determination of the Adhering-Sliding Changeover Point with Grid Refinement Convergence .. 61 5.7 Numerical Results for Stick-Slip Contact: Conver-gence of Displacements with Grid Refinement 71 5.8 Numerical Results for Stick-Slip Contact: Conver-gence of Stresses with Grid Refinement 71 vi LIST OF FIGURES Figure Page 3.1 Adjacent adhering and sliding nodes 26 4.1 Flowchart for computer program to solve contact problems with a fixed contact region 38 4.2 Flowchart for computer program to solve pro-gressive contact problem . 5.1 A sliding contact problem 42 5.2 The Ix i Mixed finite element grid for Problem A 44 5.3 The 1x1 Mixed finite element grid for Problem B 46 5.4 The Mixed finite element grids for Problem C ^9 5.5 Free body diagram of the rigid block problem 50 5.6 Strain energy versus grid size for Problem C 52 5.7 Stress distributions on the contact face for Problem C 53 5.8 The deformed block in Problem C 55 5.9 A stick-slip contact problem 56 5.10 The QST and Mixed finite element grids for the stick-slip problem 59 5.11(a-b) Strain energy versus grid size for the finite element model for the stick-slip problem 62-63 5.12(a-d) Stress distributions on the contact face for the finite element models for stick-slip problem 65-68 5.13(a-b) The T x x distribution on the centre line of the block for the stick-slip problem 67-70 5.14 The deformed block in the stick-^jjp problem 7 2 vii Figure Page 5.15 A progressive contact problem 73 5.16 The CST and Mixed finite element grids for the progressive contact problem 75 5.17 Applied loading versus the size of contact width 77 5.18 Applied loading versus relative approach 78 5.19 Normal stress distributions for the progressive contact problem 79 Appendix A.1 A Mixed triangular element 86 A. 2 A Mixed element in a local coordinate system 92 C.I Contact nodepairs 99 viii ACKNOWLEDGEMENTS The author expresses his gratitude to his advisor, Dr. M.D. Olson for his patient guidance and invaluable advice during the research and the thesis preparation. He would also like to thank Dr. N.D. Nathan and Dr. D.L. Anderson for their valuable advice during this work. Much grateful appreciation is due to the National Research Council of Canada for the financial support of this work. NOMENCLATURE A list of symbols used repetitively in this thesis is given here. Symbols for summation indices and other conventional mathematical symbols are not included. Symbol Description A the flag indicating adhering condition in the computer program C the flag indicating sliding condition in the computer program d the distance between the nodes in a nodepair E the elastic matric relating a to e F a functional f_ the force density vector in linear elasticity f ,f the consistent loads on the variables T , T , T T xy yy xy yy 7 7 7 respectively, of the mixed finite element model , J a functional k an elastic spring constant k m the multiplicative factor for A Q_ to invoke a change in the contact condition at the ith node pair with the mth increment L the area coordinates of a triangle X Symbol Description ml the direction cosines of the outward normal on S with respect to the x-y cartesian coordinates p surface traction AQ_ a test incremental consistent load vector Ag the test incremental solution vector corresponding to AQ R the cumulative consistent load vector r the cumulative consistent solution vector corresponding to R S a boundary of a domain T the differential operator matrix in the equilibrium equations of linear elasticity u the displacement vector in linear elasticity u,v the components of u in the x,y directions, respectively, of the x-y cartesian coordinates e the strain vector in linear elasticity eX ' £y'^Xy * r , e components of e with respect to the x-y cartesian coordinates y Coulomb's coefficient of friction Coulomb's coefficient of kinetic friction u Coulomb's coefficient of static friction p s v Poisson's ratio $ a global coordinate function cj) a local coordinate function xi c Symbol Description o_ the stress vector in linear elasticity T , T , x the components of a with respect to the x-y cartesian xx yy xy coordinates fi a bounded continuum fi* an approximate domain of fi fi a subdomain of fi* or an element domain General Subscripts A,B,etc. points on a finite element model C that part of the boundary lying within the contact region F that part of the boundary lying outside the contact region n,t denotes the normal and tangential directions, respectively T that part of the boundary where surface tractions are prescribed U that part of the boundary where displacements are prescribed x in the x-direction of the x-y cartesian coordinates y in the y-direction of the x-y cartesian coordinates ,x differentiation with respect to the variable x ,y differentiation with respect to the variable y denotes a vector quantity xii Symbol Description General Superscripts denotes the prescribed value of the variable e in the element 1 1. INTRODUCTION 1.1 Background In the design of many engineering structures it is very common to have components such as connecting rod and shaft, gear teeth in mesh, coming into physical contact with each other. It is important for the designer of these assemblies to determine the deformations and stresses at and near the contacting surfaces. Mathematical analysis of the contact problem for ideal geometric configurations has long been in existence. Hertz [1] first solved the displacement and stress distributions in the smooth contact region between two ellipsoidal bodies pressed together with the assumption that the contact region is small compared with the principal radii of curvature of the undeformed ellipsoids. Goodman and Kerr [2], later in their investigation of an elastic sphere indenting an elastic cavity, removed this restric-tion while retaining the assumptions that displacements and stresses are small so that the analysis was still within the framework of the linear theory of elasticity. Tu and Gazis [3] investigated the frictionless contact problem of a plate pressed between two spheres using the Hertz assumptions except that the radius of the contact region may not be small in comparison with the thickness of the plate. Such analytical efforts involved long mathematical developments which solved, in the realm of linear elasticity, unbonded contact problems between bodies of ideal geometric shapes. Towards the understanding of unbonded 2 contact between surfaces of more general configurations, Blackketter and Christensen [4] employed an approximate method in their study of contact between two two-dimensional nearly rectangular smooth elastic bodies. Although the magnitude of error in the iteration process was difficult to establish, the numerical values were in close agreement with experimental results. Goodman [5] applied the incremental formulation along with the Hertz assumptions to analyze bonded contact between normally loaded rough spheres. Relative slip of the surfaces in contact is prevented by the shear stresses de-veloped on the interface in contrast to the Hertz problem in which relative slip occurs without resistance. Many of these analytical solutions, especially those with incremental formulation, lend themselves to numerical methods. For practical situations where geometries are more complicated and the contact phenomenon is most likely of the stick-slip type, that is neither totally bonded nor unbonded, the need for more general numerical approaches is obvious. The recent advent of high-speed digital computers rendered the application of the finite element approach [6,7,8] possible in the study of the contact problem. The finite element method is very versatile in that it has the capability to handle various boundary conditions, to accommodate dissimilar materials and geometries, and to model different types of loads such as body forces and thermal stresses. Parsons and Wilson [9], with a potential energy formulated constant strain element, solved the frictionless contact problem between elastic bodies. Ohte [10] extended this method to contact problems 3 with friction where the irreversibility of the stick-slip phenomenon was modelled. At each increment of load, unknown contact forces were determined by an iterative procedure. Chan and Tuba [11,12] used a similar technique to study the effects of clearance, load and friction on turbine blade fastenings. Later, with a potential energy formulated linear strain element, Caertner [13] was able to model the frictional stresses by forcing the strain variables to the appropriate values according to the contact conditions on the interface. The results of this method agreed well with those obtained experimentally in the case of a connecting rod with a loose-fitting bushing. More recently, Okanoto and Nakazawa [14] introduced a technique in which contact elements are used to determine the contact stresses and de-formations. Numerical results of this approach showed a reasonable agreement with experimental data. In these studies the finite element method is applied in conjunction with the minimum potential energy principle. Investi-gations in this area with the complementary energy formulation or the mixed formulation have been lacking. In the complementary energy formulation, stresses, instead of displacements are the in-dependent variables and the mixed formulation has both the dis-placements and stresses as independent variables. In studies by Dunam and Pister [15], Mirza [16], Mirza and Olson [ 17,18] it was observed that in plane elasticity problems and particularly in stress sing-ular situations, finite element analyses applied with the mixed variational 4 principle yielded results more accurate than those from corresponding potential energy formulated models. Since the contact problem is a stress concentration problem which, exhibits a stress singularity at the edge of the contact region as Goodman noted, further studies on this problem employing the mixed method seems warranted. 1. 2 Purpose and Scope In this work, the mixed finite element method is applied to two-dimensional elastic contact problems. Contact constraint con-ditions pertinent to the stick-slip behaviour of Coulomb friction are developed. In contact problems where the contact region is independent of the loading, an iterative scheme is devised for the determination of adhering and sliding portions of the region. In progressive contact problems where the contact region is a nonlinear function of the applied loading, an incremental approach is employed to model the geometric non-linearity. Fortran computer programming is used to obtain numerical results of displacements and stresses in each problem investigated. 5 2. THEORETICAL BACKGROUND 2.1 Reissner Principle 2.1.1 Background In the linear theory of elasticity, the three quantities stresses, strains and displacements denoted by a , e and u are related by three field equations:-Equilibrium equations: - T T a = f . . . . .(2.1) Strain-displacement relations: T u - e . . . . . ( 2 . 2 ) Stress-strain relations: a = E_ e . . . . (2.3) In these equations, T is a differential operator matrix, f is a vector of force density components and E is the elastic matrix governed by material properties. It is common to combine Equations 2.2 and 2.3 to form the stress-displacement or constitutive relations: a = E T u . . . . (2.4) 6 The three widely used variational principles in elasticity problems are 7 the principle of minimum potential energy, the principle of minimum complementary energy and the Reissner [19] mixed variational principle. In the potential energy theorem, the constitutive relations are sub-stituted into the equilibrium equations, by which the variational equations become equivalent equilibrium equations in terms of displace-ments. In the complementary energy theorem, however, the equilibrium equations restrict the class of admissable stresses and the variational equations become equivalent to the constitutive relations. The Reissner principle is a more general variational principle in that the resulting variational equations are both the equilibrium and constitutive equations. 2.1.2 The Variational Principle In plane elasticity, the linear theory furnishes explicitly the equilibrium equations: - T XX, x T = f, x - T - T = f .(2.5) xy,x yy,y y and the constitutive (stressTdisplacement) equations: 7 u, - -F (T - VT ) = 0 x E xx yy V , - = ( V T + T ) y E xx yy u, + v, y x |(1 + V) T E xy (2.6) For a region A in the plane continuum bounded by a curve S, Equations 2.5 and 2.6 are solved subject to these boundary conditions on S: on S-p or u = u v = v on S U (2.7) where p = Six + mx *x xx xy p = £x + mx y xy yy and l,m are the direction cosines of the outward normal on S. A functional F is defined as follows: F = If [T U , + X v, + x (u, +v, ) - T ^ - ( X 2 + x2 i 1 xx 'x yy 'y x y 1 'y 'x' 2E xx yy 2vx x + 2(1 +v)x 2 )] +uf +vf ]dA xx yy xy x y [up x + v p y ] dS [(u-u)p x +(v-v)Py]ds 'U . . . .(2.8) 8 Reissner's theorem states: Among all states of stress and displacement which satisfy the boundary conditions, the actually occurring state of stress and displace-ment within the small displacement theory of elasticity renders 6 F to zero. To illustrate the theorem, the functions u, v, T , T xx yy and T are given arbitrary variations 6u, 6v, 6T , <5 T and ST X y => i ' ' xx yy xy yielding the first variation of F: <5F = // [6T ( U , -1 (T - VT )) ^ xx x E xx yy + 6 T ( V , - i ( T - VT ) ) yy y E yy xx + 6T ( U , + V , - 2 T ) xy y x E xy Suf - 6vf x y + Su, T + Su, T + Sv, T . + S 'V, T ]dA x xx y xy y yy x xy / (p 6u+p 6v)dS - / (p 6u+p 6v + (u-u) 6 p v +(v-v) Sp) ds . . . . .(2.9) x y Subjecting the terms containing variations of displacement derivatives, that is 6u, , 6u, , 6v, , and 6v, , to integration by x y x y parts yields: 9 6F = ff [6T (u, - -i- (T - VT ) A xx 1 'x E xx yy' + 6 T (v, — F - ( T - V T ) yy 'y E yy xx ' + 6T (U, + v, - 2 ^ 1 ^ T ) xy y x E xy + 6u (-T - T - f ) xx,x xy,y x + 6 V ( - T - x - f ) dA YY>Y xy,x y' + f~ i C [6U ( £T + mx )+ 6V ( £ T +ITIT )]ds S T +S U xx xy xy yy + / - (p 6u+p 6v)dS - /_ (p 6u+p 6v + ( u - u ) 6p + (v - v) 6p ) ds . . . . . (2.10) x y According to Lagrange's lemma, the vanishing of 6F for arbitrary variations 6u, <5v, 6T , <5 T , and 6T furnishes the followinq xx yy xy 3 differential equations in the domain and conditions on the boundary: u, - i (T - VT ) = 0 x E xx yy v, - i (.T - VT ) = 0 y E yy xx u ' y + v ' y " I ( 1 + v ) Txy = ° ' ' ' ' , ( 2 ' 1 1 ) 10 or •T - T - f xx,x xy,y x T xy,x - Tyy,Y - f y = 0 • ' • • - --C2.12) IT + mx - p = 0 xx yy "x on Sj IT + mx - p = 0 . . . . . (2.13) xy yy y u - u = o on Sy v - v; = 0 . . . . .(2.14) Therefore the vanishing of the functional 6 F is equivalent to the plane, elasticity problem posed by Equations 2.5, 2.6, and 2.7. Note that this implies only a stationary value for F, in general neither a minimum nor a maximum. In the context of calculus of variations, the Euler equations of the area integral in functional F are the governing differential equations in plane elasticity and the line integrals in functional F are added for the satisfaction of prescribed boundary conditions. 11 2. 2 The Finite Element Method The finite element method is a technique of numerical analysis which provides approximate solutions to continuum problems in which appropriate differential equations and boundary conditions are imposed on unknown functions, Oden [6], Zienkiewicz [7], and Oliveira [8]. In the formulation of the method, the bounded continuum fi is replaced by a domain fi* such that fi* may be subdivided into a number, say E, of non-overlapping subdomains, fi , called elements. The construction of elements is such that adjacent elements are to share a common boundary, the sum of the subdomains fi equals domain fi* and that fi* tends to the bounded continuum fi as the number of fi becomes large. Each element contains a family of fields which will constitute the approximate solution in that corresponding part of the continuum. For each function u sought in the problem, an approxi-mate function is assumed in each element which can be expressed as: n u G = I uf ?(x.) = 1 if k = j K ' ; k,j = 1 , 2 , - , n . . . .(2.15) = 0 if k ± j 12 This ensures the linear independence of the in fie. It is common, however, to first choose any u containing n linearly independent t terms, for instance a polynomial, and by means of a linear transforma-g tion, cast u into the form expressed in Equation 2.11. In the formulation of the finite element method, it is only necessary to assume coordinate functions defined over individual elements. However, for the implementation of the method, it is convenient to introduce functions defined over fi* such that e e e $ k (x.) = k (x.) for x.efi = 0 for x . £ f i e . . . . .(2.16) Then for the entire domain, the approximate function is now E m u = • E £ uf $f . . . .(2.17) e=1 k=1 The adjacent elements are made compatible by specifying that the values of the approximate functions be the same at coincident nodes (and thereby providing continuity along all common element boundaries). The nodal values, or the degrees of freedom, of each element is related to those of the domain fi by the relationship: M u k = . E , F fk u i • • • - ( 2 - 1 8 ) 13 where if node u. is coincident with u. i 0 otherwise and M number of degrees of freedom in domain fi u. nodal value at nodes in domain fi . This is directly analogous to the well known displacement method in structural analysis relating local to global degrees of freedom, as described in detail by Bathe and Wilson [20]. Finally, in terms of the global degrees of freedom, the approximate function is expressed as: E m M Refined approximations of the sought functions may be achieved by subdividing the domain fi* into smaller elements, by increasing the number of coordinate functions per element or both. 2.3 Application of the Finite Element Method to Reissner's Principle methods showed encouraging results although their mathematical properties were not yet well understood. Recently, generalizations, convergence and completeness criteria of the mixed methods have been studied by Oden [21], Reddy and Oden [22]. Mirza [16] extended u .(2.19) Early applications of the finite element technique to mixed the investigation and established the convergence and completeness criteria for the mixed finite element method. On these grounds, the basic procedures in applying the finite element technique to the Reissner principle will be presented here. The approximate solution which the finite element method furnishes for the differential Equations 2.5 and 2.6 and boundary conditions 2.7 is determined by rendering the functional F in expression 2.8 stationary. e e Let T and u take on approximate functions within an element as denoted by the superscript e : m T = £ tr. . ; e = 1,2, . . . / ( 2 . 2 0 ) k=1 k k Substitutions of these into functional F yields a function F in terms e e of the element degrees of freedom x ^ and u^ . For stationarity the derivatives i f - and K 9 u k are set to zero; resulting in m+n equations. The matrix of coefficients, or the stiffness matrix, is symmetric and indefinite. This process is demonstrated in Appendix A for specific coordinate functions. 15 To form the global structure, the element stiffness matrices are assembled in a manner identical to the displacement method in structural analysis as mentioned in the previous section. In the potential energy approach, the addition of columns during the assem-blage of element matrices corresponds to equating the nodal displace-ments of adjacent elements and the adding of rows corresponds to the summing of equilibrium equations. In the mixed approach, the same holds true for the rows and columns corresponding to the displacements variables. For the stress variables, the addition of columns corresponds to equating the nodal stresses of adjacent elements and the adding of rows corresponds to the summing of constitutive conditions. From another angle, the process simply transforms the element degrees of freedom to the structure degrees of freedom. Boundary conditions are applied as suitable constraint equations or consistent right-hand-side vector entries. This is explained further in the appendices. The indefiniteness of the matrix of coefficients poses no real difficulty in solving the equations by the method of Gaussian elimination with partial pivoting. 16 3. APPLICATION TO THE CONTACT PROBLEM 3.1 Contact Conditions Most of the everyday physical systems have forces transmitted to them through a mechanism commonly known as "physical contact." For instance, consider a ladder leaning against a wall. The ladder is supported by the sum total of countless electromagnetic interactions among the atoms in the adjacent surfaces of the ladder, the ground and the wall. From an engineering point of view, a submicroscopic analysis of such interactions is unnecessarily complicated and it is sufficiently accurate to lump them into a single force, the compressive contact force. Confining the discussion to contact between solids, the contact forces are very short ranged, in fact, they are negligible when the objects are more than a few atom diameters apart. In macro-scopic investigations, the assumption that adjacent points on surfaces in contact occupy the same space is therefore reasonable. In studies of the contact phenomenon, the problem may be classified into one of three categories: unbonded contact, bonded contact and stick-slip contact. These are just empirical descriptions found to hold for some surfaces and are really manifestations of the same event. Contact between well lubricated metallic surfaces may be modelled as unbonded. In the contact regions, the shearing com-ponent of the contact force vanishes, leaving only the compressive 1 7 normal pressure-.': Itfollows then relative slip can occur without impediment. Static equilibrium and kinematic compatibility are maintained in the normal direction within the contact region. For perfectly rough surfaces coming into contact, bonded conditions may govern. When resistance to relative slip occurs, shearing stresses will develop on the interface. If these shearing stresses are sufficient to prevent any relative slip between adjacent points on the surfaces after they have come into contact, the conditions developed are those of a bonded contact. Within the contact region, equilibrium and compatibility are maintained: in all directions. In other words, continuum is established between the surfaces in contact and all governing differential equations for a continuum must likewise hold true across the contact interface. The stick-slip type of contact is the one mostly encountered in daily experiences. Consider a book resting on a table. If a small horizontal force is applied to the book, the book may remain stationary. As the force increases, the book remains at rest until a critical value is reached at which the book slips. From static equilibrium considerations, it can be said that the friction force rises with the applied force up to the critical value called the maximum static friction force. Once in motion, the friction force typically drops to a lower level and is maintained independent of the speed of the slip. For a given pair of surfaces, the magnitudes of these frictional stresses are proportional to the normal stresses pressing the two surfaces together and may be expressed in these relations: 18 T . ^ U - T nt H s nn s n t k ^ k nn | y k l ^ |us | £ 1 , typica l ly . . . .(3.1) where T ^ is the static frictional shear ing stress pr ior to slip s x x is the kinetic frictional shearinq stress du r i nq slip nt. a => r-k y , y. are coulomb coefficients of static and kinetic fr ict ion S K x is the normal pressure between the surfaces in nn contact. In the normal direct ion within the contact reg ion, static equil ibrium and kinematic compatibil ity are maintained. In the tangen-tial d i rect ion, the absolute value of the shearing stress is limited to a fraction of the compressive normal stress with its d irect ion opposite to the direct ion of relative motion. In static systems, it is plausible that some region of the contact are adher ing while others are s l id ing. Th i s s t i ck - s l ip type of contact is character ized by contact fr ict ion stresses va ry ing from zero to the allowable maximum. 19 3.2 The Contact Model The problem of contact with Coulomb friction is in essence a boundary variable constraint problem. A manipulation of the Reissner principle on the contact boundary to suit this purpose is presented here. In the plane stress problem with zero body forces, the differential equations are: 0 .0 9 9x 0 9 ~9y u i •) 0 0 0 0 3 9y 9 9x V 0 9 3x 0 1 ,E V E 0 T X X 0 0 9 9y V E 1 E 0 T yy 0 '' 9 9y 9 9x 0 0 2(1+v) E T xy 0 >. . . . .(3.2) The functional principle for Equation 3.2 with homogeneous boundary conditions can be written as: J = II [x U , +T V , +T ( U , + V , )~-iU(T2 + T 2 -.. ^ xx x yy y xy y x 2E xx yy " 2 v T x x T y y + 2 ( 1 + V ) T x y ) ] d A . . . .(3.3) 20 Taking the first variation of the functional in Equation 3.3 and integrating by parts the variations of derivatives yields: 6 J = / / [6T (U, - J - (T - VT )) . xx x E xx yy + 6T ( V , - -JL- ( T - V T )) yy y E yy xx x i 2(1+v) , + 6 T (U, + V , \ - T ) xy y x E xy ' * 6 U ( - T - T ) xx, x xy,y + 6 v (-T -* T. .::)•] dA yy.y x y * x + / s [ 6 u ( A T x x + m T x y ) + 6 v ( £ T x y + m T y y ) ] ds . . . .(3.4) where S is the total boundary of the domain containing both Sy and V On Sy where the displacements are prescribed, the variations 6u and 6 v must be zero and therefore, the integral over Sy vanishes Hence the stresses on Sy are unknowns and emerge as part of the solution. On S T where the variations 6u and 6v are not zero, the vanishing of 6J requires the stresses to satisfy the assumed homo-geneous conditions as natural boundary conditions. In the contact problem, where the contact boundary shear stress at a point is bounded by a fraction of the compressive normal 21 stress at the same point, it is desirable to have stresses as forced boundary condit ions. Hence, over the contact boundary the conditions are neither exclus ively displacement nor stress but mixed. The d i s -placements on the contact boundary will become natural boundary condit ions. An inspection of the variation of the functional 3.4 shows that an addition of a line integral - /g 5 T n n u * D + 6 T n t u * 1 + T n n 5 u * - + T nt^ u " - d s " ' ' ^ 3 , 5 y c where S c is the part of boundary in contact will achieve this end . T h u s a functional F def ined as: J - fs p u .+ p v ds . . . .(3.6) c y will require displacements be natural boundary conditions over S c and forced elsewhere and stresses be forced boundary condit ions .over S and natural elsewhere, c The Reissner pr inc ip le with non-homogeneous boundary c o n -dit ions; may now be written as: F = J " fS P x U + P v V d s c y fc up +vp d S - / _ [ (u -u)p + (v -v )p Ids* 3 F T y FU y + / S u ( P x ~ P j + v ( p - p ) d S + /«. up +vp d s • • .(3.7) 3 C T y y ^ C U y 22 where Sp is that part of boundary not in contact. S C T U S c u " S C S F T U S F U S F S c [I S p =-. S . . . .(3.8) subscript T denotes that part of the boundary where surface traction is prescribed, and subscript U denotes that part of the boundary where surface displacements are prescribed. In forming the finite element matrix equations, the element stiffness matrix K is the matrix of coefficients formed from differentiat-ing the domain integral with respect to the variables u, v, T x x » T yy ' and T x v ' ' n the case of an element lying on the contact boundary, there is an additional contribution from the boundary .integral S .^; which is just one side of the element. Along Sp, the surface tractions are consistent loads on appropriate displacement variables as indicated by the boundary integral over Spy. Let f . be the con-i sistent load on the displacement variable u. then f is qiven by: i u. i f u. = TOT, [ / S p T " P x + v p y d s l where u = u.(J).(s) and r i • • • • (3 • 9) v = y.<|>. (?) . 23 On S p y , the displacements are constrained to suit the p r e -scr ibed condit ions. S imilarly, along S^,, the displacements are c o n -sistent loads on appropriate stress variables over S ^ y . Let f be T x y j the consistent load on the stress variable x then f is q iven xyj -r 3 by : x y j f = — — [/_ u(£x> + mx ) + v ( £ x +mx )ds] H x Y i s c u x x x y x y y y where x = -T, 4>.(s) ; x = x x x , x y y , x x y . . . . .(3.1 On S^-j., the stresses are constrained to suit the fr ict ion condit ion. If the element is oriented in a local normal-tangential co-ordinate system, then the constraint condition is simply sett ing the tangential shear stress to a fraction of the compressive normal s t ress . Fur ther demonstrations are g iven in the Appendices A and B. In the bonded, or adher ing type of contact where the tangential shear stress is less than the allowable fraction of the com-press ive normal s t ress , relative slip does not occur between the surfaces of contact. Continuum is establ ished across the interface. In the case where the body is against a r ig id sur face, the normal and tangential displacements within the contact region are both zero. In the finite element model, this results in zero consistent loads on the stress variables on the contact boundary . In the s l id ing type of contact where relative slip does occur , the tangential shear stress is constrained to the compressive 24 normal stress multiplied by Coulomb's coefficient of friction y. When y is zero, the contact is unbonded. It is convenient to work in local normal-tangential coordinates so that both displacements and stresses are defined with respect to the boundary normal. Continuum is maintained across the region of contact in the normal direction. The tangential displacements of the nodes slipping relative to each other enter the matrix equation as consistent loads on the tangential stress variables. The need for an iterative procedure is obvious as the tangential displacements in the equation formulation should agree with those in solution. An iterative procedure where the displacements in the solution of the previous iteration are used to calculate the con-sistent loads in the present iteration is employed. It will be observed that convergence with such iterations is rapid. After the friction con-straint is applied on the stress variables of the nodes in contact as described in Appendix B, the system of equations is solved iteratively. Schematically, the iterative procedure may be put into this form: [K] [ 6 ^ = [f,] DK] [ 5 n ] = [ f n ] . . . .(3.11) where [f ] = function of [6 ,] . n n-1 The stopping criterion is the convergence of [<$n] or [6 n ] " " " [ f ] within a tolerable limit. 25 In the stick-slip type of contact, an iterative scheme is employed to determine the unknown adhering and sliding regions. The designated contact nodes in the finite element representation are assumed to be adhering nodes in the first solution of the matrix equation. The tangential stresses at these nodes are then compared with the normal stresses to check against the assumption of adhesion. If at any of these nodes the tangential shear stress exceeds u times the compressive normal stress, the friction condition is applied at that node constrain-ing the stress variables. These are nodes with the 'sliding' condition. The re-formulated equations are solved iteratively for the displacements in the sliding contact region. When convergence is within tolerance, the conditions at the contact nodes are examined for comparison with the assumptions made in the previous solution step. At the nodes assumed to be adhering, the stresses are compared with the Coulomb friction condition as described above. At the nodes assumed to be sliding, the tangential relative slip should be in the opposite direction to the assumed tangential shear stress. If not, the node is released from the friction constraint and re-defined to be an adhering node. A revision in the contact condition of any contact node prompts the next step in the iterative scheme. A solution to the problem is obtained when the contact conditions in the solution coin-cide with those in the previous assumption. The accuracy in determining the locations of points at which adhering conditions change to sliding conditions is governed by the size of the finite element grid along the boundary of contact. 26 Suppose that on the contact boundary, there are two consecutive nodes at one of which the contact condition is adhering and at the other slid-ing. Since the stick-slip behaviour is a discrete change in the boundary condition, and the enforcement of contact conditions is performed at the nodes, it can only be said that in the finite element model such a stick-slip changeover point lies between these two nodes. An examin-ation of the stresses and the displacements between the two nodes shows contradictory conditions. Consider body R resting against a rigid boundary C. Also consider points A and B of R which lie on the contact boundary and that the variations in the stresses and dis-placements are" linear between them as shown in Figure 3.1. Figure 3.1 Adjacent adhering and sliding nodes 27 If the contact condition at A is adhering and at B is sliding then IT IT I A < ii,-1' IT IT L = u 1 xy yy 'A K ' 1 xy yy 'B H and u A = 0 , | U b | > 0 . . . . .(3.12) It follows then | T X y / Tyy I < ^ a n c ' ! u I = u a * a n y point between A and B. Thus while the stresses indicate adhering behaviour, the dis-placements do not. The region over which this anomaly occurs may be reduced only by refining the finite element grid along the contact boundary. 3. 3 Rigid Body Modes The following three independent rigid body modes will reduce the functional J in 3.3 denoting strain energy to zero: (iii) u = -cy, v - ex; T X X = o. T yy = o, T xy = 0 : T X X = o. T yy = o. T xy = 0 T X X = o. T yy = o. T xy = 0 (3.14) The finite element approximations should exhibit the rigid body modes of 3.14. By specifying kinematic boundary conditions, the rigid body modes may be removed. A detailed eigenvalue analysis of stiffness matrices of mixed elements may be found in Mirza [16]. 28 For plane triangular mixed elements, the eigenvalue distribution finds three positive, three zero and three negative. The eigenvectors corresponding to the three zero eigenvalues are linear combinations of the rigid body modes in conditions 3.14. For a contact element, one side of which has the added line integral L; /g p x u + p vdS , the compliance with the zero strain energy rigid body modes is not affected. An eigenvalue test on the stiffness matrix of such an element, however, will yield no zero eigenvalues. This follows directly from the condition that prescribed displacements on are consistent loads on appropriate stress variables: Consider the eigenvalue equation (K - XI) 6 = 0 for zero eigenvalues, if any _K 6 = 0 . . . . .(3.15) Since any non-zero displacement pattern 6 will cause consistent loads (entries in the right hand side vector), there cannot exist any zero eigenvalues. The reduced stiffness matrix of an element having constrained stress variables at the nodes ending line integral L in the manner described in Appendix B also includes rigid body modes with zero strain energy. It is interesting to note, however, that it possesses a zero eigenvalue, that is for particular rigid body modes, the required consistent load vectors are zero. This is a direct result of symmetriz-ing the constraint equations which leads also to a reduction in the right hand side vector. For example, consider an element with line integral L 29 between nodes 2 and 3 and that the tangential direction is defined by side 2-3. Suppose that at both nodes 2 and 3, the tangential shear stress variables are constrained to the normal stress variable multiplied by y. For the particular rigid body translation when the normal dis-placement is - y times the tangential displacements, the consistent load on the normal stress variables at nodes 2 and 3 is - y times that on the tangential stress variables. When the constraint is applied to the system of equations, the right hand side vector reduces to zero which is the zero eigenvalue equation. For the element shown in Figure A. 2 the eigenvectors corresponding to the zero eigenvalues when y is 0.2 and 0.3 at both nodes 2 and 3 are respectively: y = 0. 2 y = 0.3 " 0.566139 " " -0. 553001 v 1 -0. 113228 -0.165900 T xx 1 0 0 T yy. 0 0 ^ 1 0 0 u 2 0.566139 -0.553001 V 2 -0. 113228 -0. 165900 T = 0 1 0 T y y 2 0 0 ET ) x y 2 u 3 0.566139 -0.553001 v 3 -0.1 1 3228 -0.165900 T x x 3 0 0 T 0 0 (T ) x y 3 . . . .(3.16) 30 To eliminate these rigid body modes, kinematic boundary con-ditions are applied through consistent loads on the stress variables. While the kinematic stress variables present no real difficulty in incorporating prescribed displacement condition, it does warrant an iterative procedure in solving the matrix equation when the natural dis-placement variables are sought in the solution. As mentioned pre-viously, this follows from the need for reconciliation of the consistent loads on the stress variables and the displacements in the solution. 3. H Progressive Contact In the previous sections, the types of contact are described with the assumption that the region of contact is known. The pre-determined contact region could be a stage of contact of a more general contact phenomenon, progressive contact. Consider a rubber ball squeezed between two boards. If the two boards are to approach each other, the contact regions formed between the ball and the boards will become larger. Formally, pro-gressive contact is the contact phenomenon where the contact region is a function of the relative approach of two bodies in consideration, or a function of the forces pressing them against each other. In the finite element method, the portions of the boundaries expected to come into contact are designated as the contact boundary on which lie designated contact nodes. In the case of two bodies coming into contact, these nodes on opposite sides of the contact inter-face constitute contact nodepairs. Continuous progressive contact 31 becomes discretized successive closing or opening the gaps between the nodes of contact nodepairs and changing the type of contact at nodes already in contact. The problem reduces to the determination of the load increments required to bring the two bodies from one con-tact stage to the next. Suppose that at the end of the previous load vector incre-ment A R , the contact condition at the ith contact nodepair is C 1 m ^ m denoting open or closed and A* denoting adheringv or sliding conditions at the mth stage and that the cumulative load vector is R and the 3 m cumulative solution vector is r . The problem, then, is to establish an incremental load A R , that would cause a chanqe in the contact m+1 3 condition of one and only one nodepair. Assuming a test load increment A Q, the system with contact conditions C 1 , A 1 is solved to obtain a test increment solution vector m m A q. Assuming linear elastic behaviour, the test load increment AQ must differ from the required AR , by the same factor k , goveminq M m+1 7 m+1 a a the Aq to Ar , ratio. Suppose k1 , is the multiplicative factor for M m+1 ^ K m+1 ^ AQ to invoke a change in the contact condition at the ith nodepair, then k , , the multiplicative factor for AQ to invoke a change in the m+1 ^ 3 contact condition of the system, must be the smallest of the positive k1 , of all nodepairs. m+1 r The method of evaluation of the k1 , depends on the contact m+1 r conditions C 1 , A 1 . For C 1 open, k1 , is the factor required to m m m ^ m+1 ^ close the distance dj^ separating the ith nodepair. The calculation of .i d m is given in Appendix C. Then k m + 1 is given by 32 d" (r ) :m m m+1 d 1 (r ) - d ' ( r +Aq) m m m M ' (3.17) where d (r m +Aq) is the distance between the nodes of the ith nodepair evaluated after the test increment solution vector Aq. For C* closed, k1 is the factor required to reduce the m m+i ^ normal stress at the ith node to zero. Then k' , is qiven by m+1 a 7 m+1 T (r ) nn m m AT' (r + Aq) nn m 1 (3.18) where the subscript n denotes the direction normal to the contact boundary. For C 1 closed and A 1 adhering, k ! , is the factor required m m a m+1 M to satisfy Coulomb's friction condition, T given by nt nn Then k , is m+1 m+1 U T (r ) + T . (r ) nn m nt m m m tiATnn(Aq) +AT n t (Aq) for and T . (r ) A T „ ( A q ) > T (r ) A T .(Aq) nt m nn 1 nn m nt m m m+1 T . (r ) - y T (r ) nt m nn m m m [ A T n t ( A q ) - y A T n n ( A q ) j nn for 33 Tnn ( r m ) A T n t ( A ^ > T n t ( r m ' A W A ^ ( 3 ' 1 9 ) m For C m closed and A m sliding, k m + 1 is the factor required to reduce the tangential slip to zero. Then k m + 1 is given by u! (r ) . 1 tm m m+1 " -. . . . . .(3.20) A u ' ( r m + Aq) The derivations of the various k1 's are qiven in Appendix C. m+1 a ^ Having determined the multiplicative factor k , , hence the a ^ m+V increment load required to bring the progressive contact to the (m+1)th stage, the contact conditions over the contact region are updated. At the ith nodepair where the open or closed condition is to change, the appropriate boundary conditions are applied. If the nodes are to leave the contact region, that is the contact condition of the nodepair has become open, then all displacements and stresses are set to be in-dependent variables. If the nodepair has just come into contact, the nodes are first assumed to be adhering. That is, continuum is established between the nodes. In the case of contact against a rigid face, both the normal and tangential displacements relative to the rigid face are constrained to zero. Nodepairs coming into adhering conditions from sliding conditions and vice versa will also have their boundary conditions accordingly constrained as discussed earlier. To account for the geometrically nonlinear nature of the problem, the coordinates of the nodes are updated at each contact 34 stage in addition to the rectification of boundary conditions. The nodes in a nodepair that are in sliding contact are forced to occupy the same point in space by averaging their coordinates at the end of each incre-ment. The inaccuracy in determining the cumulative solution vector involved in this process depends firstly on the relative magnitudes between the incremental displacements and the size of the boundary elements. If the ratio of the incremental displacement to the length of a contact element is small, then taking the Ar to the new nodal coordinates is a reasonable approximation. Secondly, the coordinates of the nodepair not yet in contact should be such that the effect of coordinate averaging will be small when the nodes come into contact. A judicial choice of the initial nodepair coordinates should serve to minimize the amount of mis-match when the nodes come into contact. A solution to the progressive contact problem, is obtained when the sum of the load increments reaches the given load level. 35 4. COMPUTER PROGRAM A computer program is written in Fortran to solve two-dimensional contact problems where the contact region is known. The mechanics of the program is illustrated in the flowchart of Figure 4.1 followed by brief descriptions of the program components. The program is essentially a standard linear elastic finite element program with the added iteration loops LOAD - DGBAND and KNOCK - MAKE. In the LOAD - DGBAND loop, the matrix decomposition is performed once only. The decomposed matrix is retained for subsequent iterations for each of which a different consistent load vector is generated by the subroutine LOAD. The computing time saved is sizable. The criterion for convergence is that the consistent loads on appropriate stress variables ought to agree within a tolerable limit with the displacements in the solution. A good indicator of such convergence is the energy stored in the body at each iteration step. In the KNOCK - MAKE loop, the boundary conditions in the contact region are revised according to the agreement between the conditions on the contact boundary in the solution and the assumed ones. This may lead to new constraint equations on or removal of old ones from the previous matrix equation. The global stiffness matrix assembled originally is stored without any Coulomb friction : 36 friction constraints along the contact boundary. Any new set of stress constraints required is applied on this stiffness matrix and the awkward programming of removal of constraints is avoided. In the problems solved, the coefficient of friction u is taken as a constant. For many combinations of engineering materials, the difference between static and kinetic coefficients of friction is small. Furthermore, in static problems where displacements are small and kinetic effects are negligible, it is reasonable to use the static co-efficient of friction in the constraint model. Another computer program is coded in Fortran to solve two-dimensional progressive contact problems. The flowchart in Figure 4.2 highlighting the program is followed by brief descriptions of some program components. Note that several subroutines appear also in Figure 4.1. 37 Descriptions of subroutines in Figure 4.1 LAYOUT BANDWH CALER LININT BUILD MAKE DGBAND LOAD EXPAND TOUCH COFRIC KNOCK Geometric, material and loading data input Calculates bandwidth of global stiffness matrix Calculates the stiffness matrix of each element Adds to the stiffness matrix of each contact element the contribution from the line integral L Assembles the global stiffness matrix in a two-dimensional array and discards the rows and columns corresponding to homogeneous displace-ment boundary conditions Re-assembles the global stiffness matrix from a two-dimensional array form to a one-dimensional triple bandwidth form Solves symmetric noh:-positive definite matrices by Gaussian elimination using partial pivoting Calculates consistent loads on the stress variables using displacements in the solution of the previous iteration step Inserts into the solution the previously removed homogeneous boundary conditions Examines the contact conditions against those assumed and revises them for the next loop if necessary Multiplies the rows and columns corresponding to the T variables on sliding nodes by y and adds them to^the rows and columns corresponding to the x variables at the same nodes v v Discards the rows and columns corresponding to the x variables on slidinq nodes, xy 3 Figure 4.1 Flowchart for computer program to solve contact problems with a fixed contact region Desc r i p t i on of Sub rou t i ne s in F i gu re 4.2 SET PDIST PDIST2 KN CUMU Sets cumulat ive loads, cumulat i ve s t res ses to zero Ca lcu lates the d i s tance between the nodes of contact nodepa i r s not yet in contact Ca lcu lates the new d i s tance between the nodes of contact nodepa i r s not yet in contact a f te r the app l i cat ion of the load increment AQ Ca lcu lates the mul t ip l i ca t i ve factor k n for the test load increment to take the system from one state of contact to the nex t Ca lcu lates the cumulat ive load v e c t o r , cumulat ive s t res ses and updates the coord inates of the nodes of the model to account for the geometric n o n -l i nea r i t y 40 START STOP SET LAYOUT B A N D W H G A L E R L IN INT S E T U P EXPAND C U M U K N PDIST2 DGBAND PD IST Figure 4.2 Flowchart for computer program to solve progressive contact problem 41 5. SAMPLE PROBLEMS Results from investigations of several contact problems are presented here. The method employed is the mixed finite method described previously in Section 3.2. The mixed element used herein has linear interpolating functions for both its displacement and stress variables. The structure in each problem is modelled by finite element grids. Elements not having one of their sides forming part of the contact region follow the original formulation of Resissner's principle. Those lying on the contact region have added contribution to their stiffness matrices arising from the boundary integral along the contact boundary. Some results obtained from investigations using the potential energy formulation are also presented for comparison. 5.1 Sliding Block Consider the problem setup in Figure 5.1. H2 y / -/ / k -AAA-- W -Q / / / ? <— I -"7 7" h E = 21000.ksi v = 0.3 k = 200. ksi/ in. Figure 5.1 A sliding contact problem The unit thickness block has a Young's modulus of elasticity E and a Poisson's expansion factor v. The coefficient of friction between the block and the rigid base is y . The elastic foundation attached to the left side has a continuum spring constant k. Suppose y is sufficiently small such that when the uniform loads Q and P are applied, the frictional stress developed on the bottom face is not enough to prevent displacement. The problem is then a contact problem of the sliding type. Three different examples of this problem are considered in the following. 43 Problem A The sliding contact problem is first investigated with Coulomb's coefficient of friction set to zero, that is, a frictionless compressive contact problem. This problem is very simple but was included essentially to debug the program. The data chosen are: y = 0 . 0 Q = 40.0 ksi P = 0.0 ksi h = 1.0 in I = 1.0 in. . . . . ( 5 . 1 ) Because the exact solution is expressable by the linear inter-polating functions used within each element, the simplest 1 x 1 finite element grid shown in Figure 5.2 is used. The uniform load Q is lumped at nodes A and C as consistent loads. Similarly, the elastic spring foundation becomes lumped springs at nodes A and B. For this problem, the exact solution is obtained with no iteration. Consider the equation of constraint, T = y x , on the M ' xy K y y ' stress variables lying in the contact region. When y is zero, the x x y variables are zero and are eliminated altogether in the process of static condensation along with the consistent load formed from the u y \ V\A 10Qk/ i n ^ x y = ^ y y Figure 5.2- The 1x1 mixed finite element grid for Problem A displacements at the same nodes. The results which are exact are given in Table 5.1. A , B , C , and D refer to the corners on the finite element model in Figure 5.2. Table 5.1 Numerical Results for Problem A: A Frictionless Contact Problem A B C D u (in.) 0.0 0.0 1.904761 9 x 1 0 _ 3 1. 9047619x10~ 3 v (in.) -5. 71 42857x1O" 4 0. 0 - 5 . 7 1 4 2 8 5 7 X 1 0 " 4 0.0 0. 0 0.0 0.0 0. 0 T y y ( k s i ) -40.0 -40. 0 -40.0 -40.0 x x y (ksi) 0. 0 0. 0 0. 0 0. 0 45 Problem B The following problem has a block displacing relative to the rigid base. It is not sliding against it in the sense that no compressive stresses exist to keep the block in compressive contact with the base. The purpose of this exercise is to exemplify the convergence in the iterative process to determine the displacements along the contact boundary. The data chosen are: y = 0 . 2 Q = 0.0 ksi P = 40.0 ksi h = 1.0 in 1 = 1.0 in . . . . .(5.2) Again, because the exact solution can be expressed by the assumed functions in the finite element, the 1x1 grid is used in the investigation as shown in Figure 5.3. The results to the sliding problem B are tabulated in Table 5.2. The subscripts A, B, and D refer to the corners on the finite element model in Fiqure 5.3. The columns of values for the f ' s are the COn-lS T sistent loads on the particular stress variable calculated from the appropriate displacements from the previous iteration. For instance. y 46 \ W A 100.k/in A w /I A A ^ - 2 0 k 10O.k/in I x y z P l y y 2 - 2 0 . k Figure 5.3 The 1 x 1 mixed finite element grid for Problem B f B in the 3rd iteration is evaluated from /.-.^u ds where u is the x BD xy linear interpolation of u^ and u^ in the second iteration. Note also that f on side BD has a zero contribution from the v displacement Tyy because it is a prescribed zero displacement but has a y times f Txy contribution resulting from the static condensation process. It is observed that it takes only a few iterations to obtain a stationary approximation of the solution vector. Even with the un-realistic initialization of the f^'s to zero, the rate of convergence is rapid. TABLE 5.2 (a) Numerical Results for Problem B: A Sliding Contact Problem Without Normal Stress Iter. No. Strain Energy ( k - i n u.(in.) v B ( i n . ) T (ksi.) x x c x (ksi.) V V B x (ksi.) x y D 0 1 2 3 4 4.038091703588185 4.038095233362261 4.038095238095234 4.038095238095234 4.038095238095237 -0.20018556 -0.20000012 -0.20000000 -0.20000000 -0.20000000 ^0.03996418 {0.00003709 0.00000002 0.00000000 0.00000000 -40.039067 -40.000026 -40.000000 -40.000000 -40.000000 0.11133697 0.00007457 0.00000005 0.00000000 0.00000000 0.02226739 .0.00001291 0.00000001 0.00000000 0.00000000 Exact 4.038095238095238 -0.20000000 0.00000000 -40.000000 0.00000000 0.00000000 TABLE 5.2 (b) Iter. No. u B (in.) u Q (in.) f (kip.) x y B f (kip.) T f (kip.) f (kip.) 0 -0.19981444 -0.20171792 0.0 0.0 'o .o 0.0 1 -0.19999988 -0.20190464 -0.100224465 -0.100541715 -0.020044893 -0.020108343 2 -0.20000000 -0.20190476 -0.100317460 -0.100634860 -0.020063480 -0.020126972 3 -0.20000000 -0.20190476 -0.100634920 -0.100634920 -0.020063492 -0.020126984 4 -0.20000000 -0.20190476 -0.100317460 -0.100634920 -0.020063492 -0.020126984 Exact -0.20000000 -0.20190476 -0.100317460 -0 . 100634920 -0.020063492 -0.02012984 48 Problem C Combining Problems A and B, a general sliding contact problem is generated. The data chosen for this problem are: u = 0 . 2 Q = 40.0 ksi. P = 40.0 ksi. h = 1.0 in. I = 1.0 in. . . . .(5.3) This problem has no simple solution and therefore it is solved with a number of finite element grid of progressive refinement as shown in Figure 5. 4. For each finite element grid, the solution is obtained iteratively. The convergence with iterations is illustrated in Table 5.3 where the numerical results obtained for the strain energy in the block from each finite element grid are tabulated. Again, a rapid rate of convergence with iteration is observed. To show the convergence of strain energy with finite element grid refinement the values of strain energy from the fifth iteration are plotted against the number of elements per unit height in Figure 5.6. The somewhat oscillatory converging manner in the grid refinement process is characteristic of results obtained with the mixed method as }—w 20.k } W 1 0 0 . k / i n X x y ^ I y y 10.k 20.k 1Qk 1 W-5 0 . k/in r X y = p r y y 20. k V - 1 0 . k « * - 2 0 . k -« -10 .k 49 1 X 1 2 X 2 A X 4 l \ IM\ l \ l \ l\KK 8 X 8 Figure 5.4 The mixed finite element grids for Problem C 50 opposed to the monotonic convergence associated with the boundedness of potential energy formulations. The normal and shear stresses on the bottom face of the block are shown in the Figure 5. 7 . At the trailing end of the block where there is a stress reduction and near the leading end where there is a stress concentration, the finer grids reveal the stress gradients more dramatically. Although this problem has no exact solution, a rough compar-ison can be obtained by considering the limiting case of a rigid block. Consider the free body diagram of the rigid block in Figure 5.5. -± ± ±e-spr ing react ion R Q t x y z JJlyy tyy = A-I y y = B Figure 5.5 Free body diagram of the rigid block problem 51 Since there are only three equilibrium equations, there can only be a linear approximation for the unknown normal stress distri-bution on the bottom face. For a unit block, the unknowns are given by R = Q - y P A = P( 1 + 3 * F ) B = P( 1 - 3y ) . . . . . ( 5 . 4 ) Included in Figure 5.7 is a plot of this stress distribution (shown dashed) for the purpose of comparison. In the finite element model this block is relatively rigid in comparison to the spring, hence the 1x1 grid solutions is in good agreement with the rigid block solution. Although the exact continuum solution to the problem is unknown, the boundary conditions that must be satisfied by such a solution are definite. For instance, along the bottom edge BD, the vertical dis-placement v should be zero. However, in the present method, this condition is a natural boundary condition and therefore is only achieved in the limit of zero element size. The actual numerical convergence of the displacements at B and D are shown in Table 5.4. Also shown for comparison are the vertical displacements at A and C. It is seen that V g and V p are about 8 percent of v^ and v^,. The deformed shape of the block as obtained from the 8x8 grid is shown in Figure 5.8. 52 y N u m b e r of e l e m e n t s H i gh Figure 5.6 Strain energy versus grid size for Problem C 53 x ( i n ) Figure 5.7 Stress d istr ibut ions on the contact face for Problem C TABLE 5.3 Numerical Results for Problem C: A Sliding Contact Problem with Normal Stress Iter.No. Strain Energy (k-in.) ' 1x1 2x2 4x4 8x8 0 1 2 3 4 0.52467207 0.052482976 0.052482998 0.052483002 0.052483003 0.012431459 0. 016296002 0. 016318660 0. 016318819 0.016318826 0.012253379 0.016706516 0.016749069 0.016749315 0.016749316 0.011973734 0.016672901 0.016715101 0.016715100 0.016715100 TABLE 5.4 Numerical Results for Problem C: Convergence of Displacements with Grid Refinement GRID v B ( 1 0 4 in.) v D ( 1 0 V l . ) v A ( 1 0 4 in.) v c ( 1 0 4 in.) 1x1 2x2 4x4 8x8 - 4 . 0388415 -2.4680971 -2 .2220836 -1.9649921 -5 .1975942 -1 .1409362 -0 .9702678 -0 .5083620 -19 .686294 -17.795228 -17 .754885 -17.963476 -8 .3809920 -9 .2390954 -9.1502081 - 9 .1980224 EXACT + 0.0 0.0 - -4= Figure 5.8 The deformed block in Problem C 56 5. 2 Stick-Slip Problem Consider the problem set up in Figure 5.9. y Q * y *V >y y y y E , ^ I X y z L l t y y Tin. ~7 7—7—7—7—7—7—7—/ J> ' / < 2.in. » Q - 20. ks i E = 21000. k s i V = 0-3 Figure 5.9 A stick-slip contact problem The unit thickness block under the uniform load Q will expand in the x-direction as a result of Poisson's effect. Depending on Coulomb's coefficient of friction, frictional stresses will develop on the bottom of the block. The stresses may be sufficient to prevent slippage in some parts of the contact region (adhering portion), while at other parts slippage does occur (sliding portion). This problem then as different from the sliding problem, has unknown boundary conditions in the contact region. 57 For comparison, the problem is investigated by both the mixed finite element model and a potential energy formulated finite element, namely the quadratic strain triangular element (QST). The method of iteration for the mixed element described in Section 3.2 is used to determine the points at which adhering conditions change to sliding conditions and the displacement and stress patterns thereof. The method for the QST is the one used in previous works for potential energy formulated finite element models. In this pro-cedure, frictional stress is applied as consistent loads on the tangential displacement variables. The stresses are a fraction, as the Coulomb friction condition requires, of the normal stresses extracted from the solution of the previous iteration. Schematically, the iterative pro-cedure may be put into the same form as Equations 3.11: [K] [ 6 , ] = [f,] [K] [ 5 n ] = [ f n ] • • . .(5.5) where [f ] = function of [ 6 , ] . 1 n n-1 When [ 6 n ] T [ f n ] has converged within a tolerable limit, the iteration stops. Since the displacement variables are associated with the forced boundary conditions, an adhering condition at a node is modelled by an elimination of the normal and tangential displacement variables and sliding condition is modelled by the elimination of the normal displacement variable while retaining the tangential displacement variable as a degree of freedom. 58 Noting the line of symmetry in the problem, only half the block needs to be modelled. The grids 1x1, 2x2, 3x3, 4x4, and 6x6 for half the block are used for both the mixed element and the QST investigations. The 1x1 grids are shown in Figure 5.10. For ease of comparison, the grids 1x1, 2x2, 3x3, 4x4, and 6x6, are used for both the QST and the mixed element investigations (for one half of the problem). The iterative approximations of contact conditions on the bottom face of the block with the 6x6 finite element grid for u =0.10 are tabulated in Table 5.5. The points A, B, C, D, E, F, and C refer to the consecutive grid points lying on the contact boundary with point A at the centre line of the block and point G at the outer edge of the block. From the third iteration on, as in the case of the mixed element and fourth in the case of the QST, adjacent adhering and sliding nodes are found to be nodes B and C. In general it is observed that it takes fewer iteration for the mixed element than the QST to reach the same results using the same grid.,. Also tabulated-is the \i = 0.15 case. The adhering-sliding changeover point is found to be a function of Coulomb's coefficient of friction \i. The adhering and slid-ing regions approximated with different finite element grids are tabulated in Table 5.6. The results show good agreement between the mixed element and the QST. In Figures 5.11, the strain energy in the block is plotted against the grid size at different values of the two methods. In the QST plot, the adhering curve showed a monotonic convergence from >x ,u I MIXED MODEL 1 a d h e r i n g : u = 0 , s l i d i n g : X x y = P " C y y . t v= 0 QST M O D E L I a d h e r i n g : u = 0 , v= 0 s l i d i n g : f u = jJTyy The QST and mixed finite element grids for the stick-slip problem TABLE 5.5 Boundary Condition Results for Stick-Slip Contact: Comparison between QST and Mixed Element in the Determination of the Adhering-Sliding Changeover Point with Iteration Convergence QST MIXED Iter. y = 0.10 No. a b c d e f g a b c d; e f g 0 A A A A A A A A A A A A A A 1 A A A A S S s A A A S S S S 2 A A A S S S s A A S S S S S 3 A A S S S S s A A S S S S S 4 A A S S S S s Iter. y = 0.15 No. a b c d e f g a b c d e f g 0 A A A A A A A A A A A A A A 1 A A A A A S s A A A A A S A 2 A A A A S S s A A A A S S s 3 A A A A S S Sv A A A A S S s TABLE 5.6 Boundary Condition Results for Stick-Slip Contact:CConfparisoh between QST and Mixed Element in the Determination of the Adhering-Sliding Changeover Point with Iteration Convergence QST MIXED \ x u N. 0 .1 .2 . 3 .4. :. 5 .6 .7. .8 . 9 1.0 • 0 . 1 .2 . 3 .4 .5 .6 .7 • 8 .9 1.0 . 05 A s s A S s /10 A s s A S s 2 X .15 A s s A A s 2 * 20 A A s .05 A S S s A S S s .10 A S S s A S S s 3 X .15 A A S s A A S s 3 .20 A A A s .05 A S s S : s A S S S s .10 A S s S s A A s s s 4 X .15 A A s s s A A A s :S 4 .20 A A A A s A A A A s .05 A S '• S S S s s A s S S s S s .10 A A A A S s s A A A A S ' S s 6 x .15 A A A A S s s A A A A S S s 6 .20 A ' A A A A A s A A A A A S A 62 MIXED 3.82 3.80 i o o m >-c n t_ o> c UJ c o CO 3.78 3.76 3.74 3.72 3.70 LI=0.0 • O J J = 0 . 0 5 - op =0.10 ^ JJ - 0.15 ig JJ=0.20 a d h e r i n g Figure 5.11 1 2 3 4 5 6 Number of E l e m e n t s H igh Strain energy versus grid size for the finite element model for the stick-slip problem: (a) mixed element 63 Figure 5.11 Strain energy versus grid size for the finite element model for the stick-slip problem: (b) QST 64 below as predicted by the potential energy theorem. When frictional sliding is allowed to occur, such a behaviour is not observed. In the mixed element plot, apparent convergence from above is observed for all levels of y . A reasonable agreement is found between the two sets of results. The stress distributions on the contact face at different values of y obtained from the 6x6 grid analyses with both the OST and the mixed method are plotted in Figures 5.12.. Also shown in Figure 5.13 is a plot of a x x x distribution on the center line of the block at different values of y . Again, they indicate good agreement between the two methods. Again, as it is with the sliding problem, the exact continuum solution to the problem is unknown, but the boundary conditions that must be satisfied by such a solution are definite. The numerical con-vergence of the displacements associated with natural boundary conditions at B and D are shown in Table 5.7 for y = 0.10. It is seen that v D and are, respectively, about 3 and 4 percent of v ^ and v^,. The convergence of the stresses x and x to the stress free condition xx xy at C and that of x x y to zero at A and B by argument of symmetry are shown in Table 5.8 also for y = 0.10. They are about 2 percent of x which are also tabulated for the purpose of comparison. The x x B deformed shape of the block as obtained from the 6x6 grid is shown in Figure 5. 14. 65 V 40. k s i 12. 10. 1.0^g>f o.o t > _± ± ±_ / / / / / / 0.0 1.0 -* x 8.: o o MIXED • • QST a : J J - 0 . 0 5 b : J J =0.10 6. 4, • — — • 0.0 Figure 5.12 0.2 0.4 0.6 0.8 1.0 x ( i n ) Stress d istr ibut ions on the contact face for the finite element models for the st ick-s l ip problem: (a) shear stress 66 x ( i n ) Figure 5.12 Stress distributions on the contact face for the finite element models for the stick-slip problem: (b) shear stress x ( i n ) Figure 5.12 Stress distributions on the contact face for the finite element models for the stick-slip problem: (c) normal stress 68 Figure 5.12 Stress distributions on the contact face for the finite element models for the stick-slip problem: (d) normal stress 69 Figure 5.13 M I X E D y 40. ksi 1.0£>- >t- i ^ f o-o p ^xy= j j t y y "7~7—7—7—7 0.0 1.0 The T x x distribution on the centre line of the block for the stick-slip problem: (a) mixed element Figure 5.13 The T x x distribution on the centre line of the block for the stick-slip problem: (b) QST TABLE 5.7 Numerical Results for Stick-Slip Contact: Convergence of Displacements with Grid Refinement GRID v B ( 1 0 4 in.) v D ( 1 0 4 in.) v A ( 1 0 4 in.) v c ( 1 0 4 in.) 2x2 - 0 . 026837 -1 .010900 -18.219505 -19 .626840 3x3 - 0 . 0 3 9 8 2 7 - 0 . 837304 - 1 8 . 291 595 -19.640766 4x4 0.114818 -0 .076521 . -18 .255615 -19.641686 6x6 0.049041 -0 .672582 -18 .279924 -19 .605969 EXACT 0.0 0.0 - -TABLE 5.8 Numerical Results for Stick-Slip Contact:. Convergence of Stresses with Grid Refinement GRID x (ksi.) i x x c x (ksi.) x y c T (ksi.) X>A T (ksi.) x y B T (ksi.) x x B 2x2 -1 .171252 - 0 . 1 8 7 8 5 4 0.799319 -1 .147168 - 7 . 951464 3x3 -0 .002095 - 0 . 6 5 2 5 3 4 -0 .449692 -0 .818124 -9.773741 4x4 -0 .593755 -0.228151 0.134665 -1 .002770 -10 .501745 6x6 -0 .152562 - 0 . 2 2 3 8 2 4 -0 .233440 -0 .082713 -10 .65878 EXACT 0. 0 0.0 0.0 0.0 -73 The system of equations of the finite element model is linear with respect to the loading and nonlinear with respect to Coulomb's coefficient of friction y . Hence the position of the adhering-sliding changeover point on the contact face is independent of the applied loading but dependent on y . 5. 3 Frictionless Progressive Contact Problem The problem of a long, solid cylinder under a uniform line load resting on a rigid base is considered here. The contact formed between the deformed cylinder and the base is assumed frictionless. Shown in Figure 5.15 is a schematic axial view of the problem setup. I P / / / / / / / / I x y = P^yy '• JJ = 0.0 Figure 5.15 A progressive contact problem 74 This problem belongs to the class of geometrically non-linear problems, Martin [23], Hartz and Nathan [241, Kawai [25]. For large and small deformations alike, the conditions on the contact boundary will vary as a function of the applied loading P. It is the concern of this exercise to determine the size of the contact region, the stress distributions thereof, and the relative approach between the cylinder and the base in relation to the applied loading. The method described in Chapter 3 is used in the investiga-tion. Although the incremental formulation can accommodate large de-formations, the analysis is restricted to small deformations for which comparison solutions are more readily available. The linear displacement, linear stress mixed element is applied with the incremental formulation to the plane strain problem. Since the problem in question is frictionless, the contact boundary is strictly a Sy boundary involving the normal displacement and a S-p boundary in-volving the tangential stress. Hence the original form of Reissner's principle may be employed in the element stiffness formulation. For numerical comparison results, the mixed element is replaced by the corresponding displacement element, the constant strain triangle, CST and the problem is re-solved. For analytical comparison results. Hertz's solution to the frictionless contact between two identical long, solid cylinders is presented. By symmetry, the two cylinder problem may be reduced to the one shown in Figure 5.15. The assumption in Hertz's solution is that the contact region is small in comparison to the radius of the sylinder. The grids for the finite element analyses are shown in Figure 5.16. 75 Figure 5.16 The CST and mixed finite element grids for the pro-gressive contact problem 76 The finite element solutions are shown graphically against _ 3 the Hertz solution. The deformations are of the order 10 times the diameter of the cylinder. The contact region created under such deformations has its width up to 9 percent of the cylinder diameter. In Figure 5. 17 , the contact width, b, obtained with grid B, is plotted against the applied loading P. The agreement is reasonable since the number of elements used is few. The relative position of the CST curve with respect to the Hertz's agrees with the results obtained by Yamada [26]. Figure 5.18 shows the plot of the relative approach between the cylinder and the rigid base as a function of the applied loading. Again, the mixed element seems to yield closer numerical results to the Hertz solution than the CST solution. The normal stress distributions in the contact region at two values of P plotted against a contour of the Hertz solution are shown in Figure 5.19 . Comparison curves are difficult to show in this instance due to the large increments of load taken as a consequence of the coarse nature of the finite element grids. Also, the validity of the Hertz solution becomes doubtful when the contact width is beyond 10 percent of the cylinder diameter. Therefore, on the one hand, the small deformations will keep the number of contact nodes few and on the other hand, the larger deformations will reduce the accuracy of the Hertz solution. Bearing these in mind, the results obtained in-dicate that the mixed element applied in conjunction with the incre-mental formulation yields reasonable numerical results to the geometric non-linear progressive contact problem. 77 160. 140. 120. .8. 100. Q_ cn c -a a o 80 "O CD CL < 60 MIXED CST HERTZ A : GRID A B : GRID B .04 .08 .12 .16 Contact Width b ( in . ) .20 Figure 5.17 Applied loading versus the size of contact width 78 Figure 5.18 Applied loading versus relative approach 16. 79 14 12 to 8. i H .00 .02 -»x / / / / y ' MIXED HERTZ a : P = 110. k i p s .04 .06 x (in.) .10 Figure 5.19 Normal stress distributions for the progressive contact problem 80 6. CONCLUSIONS The problem of contact between elastic bodies requires differ-ent boundary conditions on the contact interface depending on the type of contact under investigations. In frictionless and adhering types of contact, the boundary conditions contain either prescribed displacements or prescribed stresses, respectively. In the case of contact against a rigid base, the boundary conditions are homogeneous either in dis-placements or stresses. The more general problem, contact with friction, is best modelled by a constraint on the surface shear stress in relation to the normal stress. This calls for, in the variational principle, the association of stress variables with forced boundary con-ditions. Also implied in the mixed formulation is the association of displacement variables with natural boundary conditions. Reissner's Principle is able to accommodate such require-ments by a suitable manipulation of the boundary integrals in the energy functional. When an approximation solution is sought using the mixed finite element method, the friction constraint leads to linear constraint equations on the tangential shear stress and normal stress variables on designated contact nodes on the contact boundary. The non-homogeneous displacements will enter the right-hand side vector, rendering the necessity for an iterative scheme. For the problems investigated, only several iterations were needed to achieve good agreement between the displacement input and output. 81 In stick-slip problems where the conditions on the contact boundary are sought, again, iterations are required. It is observed that only several iterations are needed to attain the final boundary conditions. It is also observed that it takes more iterations for the displacement or potential energy model in which the friction condition cannot be applied as a constraint condition. Nevertheless, the two formulations yielded comparable results. In contact problems where the contact region is known and independent of the applied loading, the positions of points at which adhering conditions change to sliding conditions are also independent of the applied loading but dependent on Coulomb's coefficient of friction u . In the geometrically non-linear progressive contact problem, an incremental formulation shows that the size of the contact region, the displacement and stress patterns are all non-linear functions of the applied loading. The iterative method for two-dimensional progressive contact may be extended to general contact between bodies of dissimilar materials. Conceptually, the technique could well be extended to three-dimensional general contact problems. 82 REFERENCES 1. Hertz, H., Miscellaneous Papers (translated by Jones, D.E., and Schott, G.A.) , MacMillan and Co. Ltd., London, 1896, pp. 246-183. (contains English translations of the papers from Journal Reine und angewandte Mathematik (Crelle), vol. 92, 1881, pp. 150-171, and Verhand Vereins Beforder Gewerberfleisses,, Nov. 1882). 2. Goodman, L.E., Keer, L.M.s, "The Contact Stress Problem for an Elastic Sphere Indenting an Elastic Cavity," International Journal of Solids Structures, vol. 1, 1965, pp. 407-415. 3. T u , Y., Gazis, D . C , "The Contact Problem of a Plate Pressed between Two Spheres," Journal of Applied Mechanics, Dec. 1964, pp. 659-666. 4. Blackketter, D.O., Christensen, H.D., "Contact Stress between Two-Dimensional Finite Elastic Bodies," Journal of Applied Mechanics, Sept. 1969, pp. 397-402. 5. Goodman, L.E. , "Contact Stress Analysis of Normally Loaded Rough Spheres," Journal of Applied Mechanics, Sept. 1962, pp. 515-522. 6. Zienkiewicz, O . C , The Finite Element Method, McGraw Hill Book Co., (UK) Ltd. , 1977. 7. Oden, J . T . , "Finite Element Applications in Nonlinear Structural Analysis," Proceedings of the Symposium on Applications of the Finite Element Method in Civil Engineering, edited by W.H. Rowan J r . , R.M. Hackett, 1969, pp. 419-442. 8. De Arantes E. Oliveira, E.R., "Theoretical Foundation of the Finite Element Method," International Journal of Solids Structures, vol. 4, 1968, pp. 929-952. 9. Wilson, E.A., Parsons, B., "Finite Element Analysis of Elastic Contact Problems Using Differential Displacement," International Journal of Numerical Methods in Engineering, vol. 2, 1970, pp. 387-395. 10. Ohte, S., "Finite Element Analysis of Elastic Contact Problems," Bulletin of the Japanese Society of Mechanical Engineering, vol. 16, no. 95, 1973, pp. 797-804. 83 11. Chan, S.K., Tuba, I.S., "A Finite Element Method for Contact Problems of Solid Bodies - Theory and Validation," International Journal of Mechanical Science, vol. 13, 1971, pp. , 615-625. 12. Chan, S.K., Tuba, I.S., "A Finite Element Method for Contact Problems of Solid Bodies - Application to Turbine Blade Fastenings," International Journal of Mechanical Science, vol. 13,':?. 1971, pp. 627-639. 13. Gaertner, R., "Investigation of Plane Elastic Contact Allowing for Friction," Computers and Structures, vol. 7, 1977, pp. 59-63. 14. Okamoto, N., Nakazawa, M., "Finite Element Incremental Contact Analysis with Various Frictional Conditions," International Journal for Numerical Methods in Engineering, vol. 14, 1979, pp. 337-357. 15. Dunham, R.S., Pister, K.S., "A Finite Element Application of the Hellinger-Reissner Variational Theorem," Proceedings of Con-ference on Matrix Methods in Structural Mechanics, AFFDL-TR-68-150, 1968, pp. 471-487. 16. Mirza, F.A., Convergence of Mixed Methods in Continuum Mechanics and Finite Element Analysis, Ph.D. thesis. University of British Columbia, Vancouver, 1 977. 17. Mirza, F.A., Olson, M.D., "Energy Convergence and Evaluation of Stress Intensity Factor K-| for Stress Singular Problems by Mixed Finite Element Method," International Journal of Fracture, vol. 14, no. 6, 1978, pp. 555-573. 18. Mirza, F.A., Olson, M.D., "The Mixed Finite Element Method in Plane Elasticity," International Journal for Numerical Methods in Engineering, 1 980 (in press). 19. Reissner, E., "On a Variational Theorem in Elasticity," Journal of Mathematics and Physics, vol. 29, 1950, pp. 90-95. 20. Bathe, K .J . , Wilson, E.L., Numerical Methods in Finite Element Analysis, Prentice Hall Inc., New Jersey, 1976. 21. Oden, J . T . , "Generalized Conjugate Functions for Mixed Finite Element Approximations of Boundary-Value Problems," The Mathematical Foundations of the Finite Element Method - With Applications to Partial Differential Equations, A.K. Aziz, edition. Academic Press, New York, 1 972, pp. 629-669. 84 22. Reddy, J . N . , Oden, J . T . , "Convergence of Mixed Finite Element Approximations of a Class of Linear Boundary-Value Problems," Journal of Structural Mechanics, vol. 2, 1973,. pp. 83-108. 23. Martin, H . C , "A Survey of Finite Element Formulation of Geo-metrically Nonlinear Problems," Recent Advances in Matrix Methods of Structural Analysis and Design, edited by R.H. Gallegher, J . T . Oden, Y. Yamada, University of Alabama Press, 1971. 24. Hartz, B .J . , Nathan, N.D., "Finite Element Formulation of Geo-metrically Nonlinear Problems," Recent Advances in Matrix Methods of Structural Analysis and Design, edited by R.H. Gallagher, J . T . Oden, Y. Yamada, University of Alabama Press, 1971. 25. Kawai, T . , "Finite Element Analysis of the Geometrically Non-Linear Problems," Recent Advances in Matrix Methods of Struc- tural Analysis and Design, edited by R.H. Gallagher, J . T . Oden, Y. Yamada, University of Alabama Press, 1971. 26. Yamada, Y. , "Recent Japanese Developments in Matrix Displace-ments Method for Elastic-Plastic Problems in Japan," Recent Advances in Matrix Methods of Structural Analysis and Design, edited by R.H. Gallagher, J . T . Oden, Y. Yamada, University of Alabama Press, 1971, pp. 283-316. 85 APPENDIX A T H E F IN ITE ELEMENT MODEL T h e derivat ion of the matrix equations for a two-dimensional isotropic tr iangular element with linear d istr ibut ions in both d isp lace-ments and stresses is given here. Cons ider the functional in the mixed variational pr inc ip le in any cartesian x - y coordinate system when body forces are zero: F (U ,T ) = / / [ 2 T T T U - T T C x ] d A - 2.4 p T u d s A J 2 0 _ (u - u) p d s (A. where xx T T ) y y xy u = (u v) ( P x V p = IT + mx " x xx xy IT +mx xy y y ' Si, m are the direct ion cosines of the outward normal. 3x " 0 f 9y 3_ 3_ 9y 9x I E 1 -v 0 -V 0 1 0 0 2 (1+v ) for plane stress and 1 - v 1 0 1 - v 1 - v 0 0 1 - v for plane strain Figure A .1 A mixed triangular element 87 For the triangular element shown in Figure A. 1 the function expansions of u and x are: u = (L T L 2 L 3 ) u. (A.2) v = (L 1 L 2 L 3 ) v. (A.3) xx ( L 1 L 2 L 3 } xx. X X , X X . (A.4) yy 0-., L 2 L 3 ) yy-yy. yy-(A.5) xy L 3 ) x y 2 l T x y 3 (A.6) where L 2 L 3 are area coordinates. Substituting A.2 to A.6 into A.1, functional F becomes 88 F ( U , T ) = 2x 'a c -vc • 0 b ' u + c -vc c b a 0 0 - 2 (A.7) Taking the first variation of Equation A. 7 and setting it at zero yields the stationary value condition as represented the following matrix equation: 0 0 T a 0 b 0 0 0 b a u d ' a 0 c - v c 0 • 0 b -vc. c 0 . 0 . b a 0 0 2(1+v)c (A.8) where the submatrices are: u ( U 1 U 2 U 3 V1 V 2 V 3 } ( T T T T T T T T T ) xx 1 x x 2 x x 3 yy} yy2 yy3 xy., x y 2 x y 3 [a] y 2 y 2 y 2 r.2 y 2 y 2 = / / U L dA A 1 J ' . . ( A . 9 ) 3 2 1 3 2 1 II L.L., dA I 1 i y [c] 12 E 1 1 II - ^ - L.L. dA A I for plane stress. {d} d. = / p L. ds 1 S x 1 {e} = / p L. ds s y' Rearranging the equations to attain nodal order yields 0 0 a. 1 0 b11 0 0 a21 0 0 0 b11 311 0 0 0 b21 C1 1 f11 0 312 0 C12 f12 c11 0 0 b 12 f12 C12 9l1 b n a l 2 0 0 0 0 0 3 22 0 C22 Q b 22 f22 symmetric C22 where [f] = [g] < 2(1+V) [c] or [K] {u} = {f-y "21 a21 0 0 g 1 2 b 2 2 322 0 0 g 2 2 13 13 '23 23 0 0 0 b, 3 1 0 0 0 "13 "23 331 0 b31 U1 d i 0.... b31 331 V1 e i C13 f13 0 T X X . 0 f13 C13 0 X YYi 0 0 0 9l3 T x Y l 0 a"32 0. b 32 u 2 d 2 0 b32 a32 V 2 6 2 C23 f23 0 T x x 2 0 f23 C23 0 T y y 2 0 0 0 923 T x y 2 0 a33 0 b 33 U 3 d 3 0 b 33 a33 V 3 6 3 C33 f33 0 T x x 3 0 C33 0 C33 T y y 3 T x y 3 0 0 (A 91 The right-hand side vector arises only from the line integral over the part of boundary S-p where stresses are prescribed. The line integral over Sy where displacements are prescribed has a value of zero but represents additional constraint equations u =u. This will be discussed further in Appendix B. On an element where constraints on stresses are sought in the solution on side 2-3 of the element, the functional F in A.1 takes the following form: F (U ,T ) .// [2x' T T U ' - T , T C T' ] dA - 2 3 S T p , T u1 DS - 2 p '"^ u 1 ds 1 S T 3 S U [u' - u ' ] T p'ds - 2 [u 1 - u']^p' ds - 2 2 S T + S U p , T u'ds + 2 2 S T [p 1 - p']"^u' ds + 2 S 2^U u p 1 ds (A.15) where the variables are primed to indicate the local x ' - y ' coordinate system as shown in Figure A. 2: 92 ^ x x 3 ^ y y 3 ^ x y 3 U2.V2 ^ x x 2 ^ y y 2 ^x\/2 Figure A. 2 A Mixed element in a local coordinate system Substituting in the linear expansions for stresses and displacements and taking the stationary value of the first variation result in a similar governing matrix equation: [K1] {u'} = {f} . (A.16) Matrix [K1] differs from matrix [K] in A. 14 in that it contains additional terms arisen from the line integral p'^ u' ds in this form of functional F. Only submatrix [b] is affected and the new submatrix is denoted by [b 1]: 93 x 3 X 2 X1 x 3 X 2 X1 x 3 X 2 X1 X 3 X 2 X1 x 3 X 2 X1 x 3 X 2 X1 0 0 X 3 X 2 x 3 X 2 X 3 X 2 X 3 X 2 / / L , L ) , y d A • L.L. ds ' J (A.17) The right-hand side vector {f1} still contains entries d and e from line integrals /L.pds taken from node 3 to 1 and from node 1 to 2. If side 2-3 lies on Sy, then 2 5 U u '"^ p' ds .(•A .18) yields consistent loads on the stress variables x , x at nodes 2 and yy x y 3: (0 0 0 0 0 0 0 0 h 2 k 2 0 0 0 h 3 k 3 ) T where L. u' ds i . . . .(A.19) 2 3 L. v' ds i If side 2-3 lies in Sy then 2 / [p1 - p ' ] T u' ds (A.20) represents constraint equations: p 1 = p1 (A.21) 9H APPENDIX B CONSTRAINTS The finite element method numerically approximates the boundary value problem by discretizing both the continuum and the boundary. The application of discretized boundary conditions is facilitated by applying appropriate constraint equations, most commonly linear, to the governing algebraic equations in the finite element model. A method of such an application is described here. The general form of i linear constraint equations on n degrees of freedom of an unconstrained finite element model is [c] {u} - {d} = {0} . . . . (B. l ) where [c] is an nxi matrix and {u} and {d} are vectors of length n representing degrees of freedom and constraint values,; respectively. Prescribed displacements and skew boundaries are particular applica-tions of Equation B.1. In the method of static condensation. Equation B.l is rearranged to the partitioned form: [c. |c.] \ J- - {d} = {0} . . . .(B.2) ' u. where j U j j represents the variables to be preserved and | U ' | those to be eliminated and [c] is an ixi invertible matrix. The choice of u. and i i u. is not unique, however. I H 95 Then {ujO> is expressed in terms of {u.}: . .(B.3) where and [T.] = - [c ] 1 [c ] j i L j [ T d ] = [ C j ] -1 Finally, {u} may be expressed in terms of | uj j {u} = u. J u. T. J u. {d } (B.4) where I is the j x j identity matrix. In cases where {d} is a zero vector, the process of static con-densation may be done by a variable transformation followed by a variable elimination, analogous to the skew rigid boundary problem where in the displacement approach the variables are transformed into normal and tangential directions with the normal displacement variable set to zero. Thus, with {d} = {0} , Equation B.4 becomes {u} u. J u. u. J = [T] u. J (B.5) 96 In the finite element model, the governing matrix equation is [K] {u}. .= {f} (B.6) Substituting B.5 into B.6 gives [K] [T] l U j j = {f} (B.7) Symmetrizing the Equation B.7 yields [ T ] T [K] [T] i U j \ = [ T ] T {f} (B.8) In the solution of B.8, the eliminated variables \ u i } are recoverable from the equation U= 1 = [T.] ^1 (B.9) The Coulomb type friction constraint is applied in this manner on the nodes which lie in the region of contact. The unconstrained vector of variables at such a node is {u} f ^ t u n Ttt Tnn CB.10) 97 where the subscripts n a n d t denote the normal and tangential directions at the node. The constraint condition is nt y T nn . . (B.1T) where y is the coefficient of Coulomb friction, the sign of which depends on the direction of frictional slip. Let be the variable to be suppressed and [T.] is readily identified: J 7 or | T n t J - = [0 0 0 y] nt [ T ] u. Ltt nn (B.12) U t u n \t T nn . (B.13) This furnishes the transformation matrix [T ] : u t ' 1 0 0 0 ' r u n 0 1 0 0 T t t • = 0 0 1 0 Tnn 0 0 0 1 Tnt 0 0 0 u Ltt nn (B.14) or 98 u n tt nn nt - [T] u tt nn . (B. 15) Substituting into the governing Equations B.8 and symmetrizing yields: [ T ] ' [K] [T] u tt nn = [ T ] ' {f} . . .(B.16) After the solution vector for the node is obtained, T . is recovered from ' nt Equation B.12. In this particular case where the constraint equations involves so few variables, the process of static condensation may be accomplished more expediently by firstly multiplying the rows and columns of matrix [K] corresponding to variable T by y and adding them to the rows and columns corresponding to variable T and secondly discardinq the x x ^ » nn 7 a nt rows and columns. 99 APPENDIX C C.I Derivation of the Separating Distance between Two Nodes in a Contact Nodepair Consider nodes i, j , k and m, n, o lying on the boundaries of bodies A and B, respectively. Nodes j and n form a contact nodepair. Figure C.1 Contact Nodepairs : •':>;..-* The normal unit vector at node j , nt has the same direction as the vector joining the center of the circle containing nodes i, j , k and node j . Let k (q.b,) (x - a.) 2 + (y - b.) 2 r. 2 . (C. 100 be the equation of such a c i r c le . Subst i tut ing (x., y . ) , (x., y.) and (x, , y . ) into C.I gives I I J J K K where ( x i (ys -b.) = r. J J (*. b.) = r. J J . ,2 2 b.) = r. J J . 2 yields a.. b. : J C 2 C 6 ' C 5 C 3 a . J C 1 C 5 " C 4 C 2 b. J C 6 C 1 " " C 3 C 4 C 2 C 4 " " C 1 C 5 c 1 = 2(x. -- x . ) C 2 = 2(y. • C 3 2 - Y J " 2 y i + 2 2 x. - X . J ' = 2 ( x k - x . ) C 5 = 2 ( y k - V C 6 2 2 y k 4 2 2 x. - X . J k ( C 2 ) (C.3) 101 Hence, nt is given by: for for for * = ' x i - V i ' b i > . / (x. -a. ) 2 + (y.-b.) 2 jj x kj > 0 i.e., c .c 5 - c 2 c 4 < 0 ; (a. - x., b. - y.) nt = —1 ) J Vf i /(x.-a.) 2 + (y.-b.) 2 J J 1 J J i x k j < 0 ; <"c4' c i } n. = i / 9 2 c. .-> .-»• ji x kj = 0 . . . . (C.4) Similarly, n^ may be determined-Let nt = (xT, yi) and n = (x , y ) . . . . .(C.5) Then the averaqe unit normal n. for nodes i and n is 3 jn ' or is Let then where 102 -y n. -> •> n. - n J n • -> -> n. - n 1 j ni ( x J - x n , y J - y n ) n = - • . . . .(C.6) / ( x ' - x ) - ( y ' t y ) By definition, the separating distance d between nodes j and n n. • n. J J " n j n = ( x ' y } d = (x . -x n ) ( x j n ) + (yj-y n ) ( y j n ) . . . . ( C 7 ) ( x J - X ) in / J ru ' 2 " j n 2 /(x J-x ) - (y'-y ) , j n, (yJ - y ) /(x J-x ) - (y'-y ) and x' and y' are evaluated from C U depending on ji x kj . 103 C.2 Derivation of the Increment Load Factor k , m+1 In elastic problems, the load to solution ratio is constant. Hence in the incremental approach to the progressive contact problem, A* R the ratio AQ/Aq is equal to m+1 where AQ, Aq are the test m+1 increment load and solution and AR , , Ar . are the increment load m+1 m+1 and solution sought to bring the contact to the (m+1)th contact state. Cumulatively, hence R , = R + AR . , = R +k , A Q m+1 m m+1 m m+1 r , = r + Ar ' = r + k , Aq . . . . .(C.8) m+1 m m+1 m m+1 M (a) From open to close condition: d' = d 1 + k' , Ad 1 = 0 m+1 m m+1 therefore i d m .k - _ ~~ — — . . . . (C • 9) Ad where Ad 1 = d'(r +A ) - d' m q m (b) From close to open condition: i i , i A i x = T + k . A T nn , nn m+1 nn +1 m 104 therefore i " I nn , i _ m_ km+1 " Ax nn (c) From adhering to sliding condition: T . + y T = K * + K , A T , +y T +k . A T ) : 1 nt , 1 H nn , 1 nt m+1 nt 1 H nn m+1 nn m+1 m+1 m m for T . A T nt nn m m+1 > T A T , nn .nt m nn + T m nt m p A T n n + A T . nn nt for "nn A T m nt nt A T m nn y T m+1 hn nt m m A T . - I IAT nt K nn (d) From sliding to adhering condition: u ' = u + k' , Au. = 0 t „ t m+1 t m+l m therefore u t . i m k m+1 A u t