LIBRARY MIchIgen State UnIversIty PLAC E. IN RETURN BOX to remove this checkout from your record. TO AVOID FINES return on or before due due. DATE DUE DATE DUE DATE DUE J T—TI JJ MSU le An Afflrmetlve AcIIon/Equel Opportunlty Inetituflon ammo-n1 APPLICATION OF THE BOUNDARY ELEMENT NIETHOD TO FINIT E ELASTICIT Y By Husain Jubran Algahtani A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Materials Science and Mechanics 1992 ABSTRACT APPLICATION OF THE BOUNDARY ELENIENT METHOD TO FINITE ELASTICITY BY HUSAIN JUBRAN ALGAHTANI A boundary element method solution for the problem of finite plane deformation of elastic compressible and incompressible solids is presented. Two types of constitutive relation are considered: the Blatz-Ko model for compressible materials and the Mooney-Rivlin model for incompressible materials. The finite elasticity solution is obtained by solving the nonlinear boundary element equations using an iterative procedure. The iterative procedure is implemented in two computer codes which can be used to solve the problem of plane finite elasticity of Blatz-Ko and Mooney-Rivlin materials, respectively. The computer codes were tested using several numerical examples. ii ACKNOWLEDGEMENTS All praise to be to Allah, Lord of the Worlds, for his providence and divine direction throughout my life and especially during my stay in the United States. I am indebted to King Fahad University of Petroleum and Minerals (KFUPM) in Dhahran, Saudi Arabia which gave me the opportunity to pursue my graduate education and supported me during my study in the U.S.A. Many people deserve special thanks and appreciation for their effort and their positive influence toward the completion of this research. I am very grateful to my advisor Professor Nicholas J. Altiero, Associate Dean of Engineering, for his professional direction, guidance and support throughout this work. I also would like to thank, the members of my guidance committee, Professor Thomas J. Pence, Professor Iwona Jasiuk and Professor David H. S. Yen, for their contributions and helpful comments. Many thanks and deep appreciation go to my family members in Saudi Arabia, my father Jubran who had to sacrifice his iii education to take care of his family, my mother Noorah for her continuous prayers, my brother Said for his moral support and encouragement and the rest of my brothers and sisters. Last but not least, I would like to express my appreciation and gratitude to a wonderful companion, my wife Noorah, the person who stood by me and was a trusted friend and a compassionate partner in all my ups and downs during my study in the U.S.A. I am also grateful to my son, Mohammad, and daughter, Sarah, for their countless kisses they shower me when I leave or come back from school. iv TABLE OF CONTENTS LIST OF FIGURES CHAPTER 1 CHAPTER 2 CHAPTER 3 INTRODUCTION 1.1 Statement of the Problem ............ 1.2 Background ............... ........... FORMULATION OF BOUNDARY INTEGRAL EQUATIONS 2.1 Linear Elasticity ... ...... .......... 2.1.1 Governing Equations . ..... ..... 2.1.2 Integral Formulation ........... 2.2 Finite Elasticity ................... 2.2.1 Governing Equations ..... ...... 2.2.2 Constitutive Law .............. 2.2.3 Integral Formulation .......... BOUNDARY ELEMENT FORMULATION 3.1 Boundary and Domain Discretization ... 3.2 Singular Integrals .... ..... ......... PAGE vii l 3 11 11 13 19 19 21 24 28 34 CHAPTER 4 CHAPTER 5 CHAPTER 6 BIBLIOGRAPHY COMPRESSIBLE MATERIALS 4.1 Plane Stress .. ....... ............... 4.2 Plane Strain ......... . ..... ......... 4.3 Iterative Procedure ... ...... ........ 4.4 Computer Program BLATZKO ............ INCOMPRESSIBLE MATERIALS 5.1 Plane Stress ................. . ...... NUMERICAL EXAMPLES 39 4O 41 44 6.1 Incompressible Plane Stress Deformation 56 6.2 Incompressible Plane Strain Deformation 71 6.3 Compressible Plane Stress Deformation 6.4 Compressible Plane Strain Deformation vi 81 85 90 LIST OF FIGURES Discretization of the boundary and the domain Flow chart for the program BLATZKO Flow chart for the program RIVLIN Uniaxial stretching of an elastic sheet Various models used for the elastic sheet problem Total edge force versus the number of boundary elements Total edge force versus the extension ratio Comparison of the BEM and the FEM solutions for the total edge force Uniaxial stretching of a sheet with a central circular hole Model for a quarter of the sheet Total edge force versus the extension ratio Comparison between the BEM and the FEM results for the total edge force Deformed shape of a quarter of the hole boundary Infinite cylinder under uniform internal pressure vii PAGE 38 46 55 59 60 61 62 63 66 67 68 69 70 74 Axisymmetric model for the infinite cylinder Applied internal pressure versus the interior node displacement Displacement profile for various internal pressure values Average hydrostatic pressure profile Radial stresses in the cylinder Circumferential stresses in the cylinder Model for a quarter of the sheet Total edge force versus the extension ratio Comparison between BEM and FEM solutions for the total edge force Simple extension of an elastic sheet Extension along X2(Kn) versus extension along X1(>\11) all versus A“ viii 75 76 77 78 79 80 82 83 84 87 88 89 CHAPTER 1 INTRODUCTION 1.1 Statement of the Problem The possible nonlinearites which may be present in boundary value problems of continuum mechanics are 1. a nonlinear strain—displacement relationship (geometric nonlinearity), 2. a nonlinear stress-strain relationship (material nonlinearity), and/or 3. nonlinear boundary conditions. In finite elasticity, the finite strain is a nonlinear function of the displacement. The constitutive law (stress- strain relationship) for a hyperelastic material is derived from a strain energy function, and the result is a nonlinear relation. Furthermore, due to the large deformation there is a distinction between the undeformed and deformed coordinates. This means that in finite elasticity, all of the above types of nonlinearity are present. 2 In addition to the above, the analysis for incompressible materials involves a constraint on the solution. That constraint is characterized, for finite deformation, by the fact that I,==1, where I,is the third principal invariant of the Cauchy-Green strain tensor C. While such a constraint is favored in analytical solutions since it leads to certain simplifications in some problems, this is not the case for numerical techniques. Moreover,the stress tensor for an incompressible material is determined within a scaler function called the hydrostatic pressure. Particularly, in the displacement-based finite element method, numerical difficulties arise when implementing the incompressible constraint since Poisson’s ratio equals 0.5 and therefore the bulk modulus is infinite and the relationship between volumetric stress and strain is indeterminate [1-3]. Many engineering materials belong to the class of nonlinear elastic materials. These materials are commonly subjected to large strains and deformations. Examples of these materials are synthetic rubbers, polymers, solid. propellant. rocket grains and biological materials such as human skin, brain tissue, and papillary muscles. Even in normal physiological functioning, these biological. materials are subjected 'to finite strain. As examples, the skin over the middle joint of a finger elongates about 20% when the joint is bent [4]. Likewise, the brain tissue may be subjected to strains of the 3 order of 30% in the situation of head injury [5]. Although many problems of nonlinear solid mechanics such as plasticity, elastoplasticity, and finite deflection of plates have been treated by the boundary element method, there has been only a very limited application of this method to finite elasticity. The main objective of this dissertation is to present a general boundary element formulation of the equations governing the finite deformation of compressible and incompressible plane elastic bodies for given constitutive laws. The formulation is then applied to obtain solutions for special cases. 1.2 Background In solving boundary value problems in finite elasticity, one must deal with highly nonlinear partial differential equations. Analytical solutions exist for only a few problems involving simple geometries and constitutive laws. This difficulty has been overcome by the use of numerical techniques. The finite element method (FEM) is the most widely used 4 numerical technique for solving problems in finite elasticity. Since most of the highly elastic materials such as rubbers and polymers are considered to be incompressible under finite deformation, most studies utilizing FEM deal with incompressible bodies. In some of these studies, however, the formulation can be reduced to the compressible case. The first applications of FEM to finite elasticity were due to Herrmann [6], Becker [7], Oden [8], Peterson et al [9], Oden and Sato [10], Oden and Kubitza [11], Oden and Key[12], and others. In all of the references cited above, the "mixed formulation" was used to incorporate the incompressibility constraint. The general procedure involved deriving the stiffness equations based on the principle of stationary potential energy. The resulting nonlinear algebraic equations were then solved by a suitable algorithm such as the Newton—Raphson method. Oden [13] also used.the method of incremental loading where, during each increment, the body was assumed to behave linearly. At the end of each increment, the nonlinear equations of equilibrium were satisfied using iterations. The above procedure was applied to solve some problems of plane stress and plane strain for the so-called Mooney-Rivlin material. Murakawa and Atluri [14,15] developed a hybrid finite element 5 formulation based on a complementary energy principle. Their variational.principlezhad ,as independent variables, the first Piola-Kirchoff stress tensor, a point-wise rotation tensor, the hydrostatic pressure and an interelement boundary displacement field. The incompressibility constraint was introduced through a Lagrangian multiplier applied to the strain energy function of the material. One of the approaches which has been employed to incorporate the incompressibility constraint into the FEM formulation is the reduced/selective integration technique [16-18]. In this technique, the volumetric terms in the stiffness matrix are integrated using a lower-order Gauss integration than is used to integrate other terms. Malkus and Hughes [16] showed that this selective integration can be equivalent to a mixed formulation, with the number of integration points used to evaluate the volumetric terms equal to the number of pressure degrees of freedom used in the mixed formulation. Argyris et al. [19] presented a method for incompressible analysis in which the strain energy function was divided into two parts: a deviatoric part and a dilatational part. The usual finite element discretization was applied to the deviatoric part of the strain energy such that the nodal displacement became the single field. variable. The dilatational strain energy was, however, assumed to be 6 associated with another continuum behavior which was defined by the nodal displacements as well as another unknown vector representing the lack of compatibility between the two continua. More recently, in a paper by Chang et a1. [20] a finite element formulation for the large strain analysis of nearly incompressible materials was presented. The finite element equations were derived on the basis of a perturbed Lagrangian variational principle from which both the displacement and the hydrostatic pressure fields were independently approximated by appropriate shape functions. There are many other papers on the application of FEM to finite elasticity, but.a complete bibliography on this subject is beyond the scope of this dissertation. The first application of the boundary element method (BEM) to geometrically nonlinear elastostatics was due to Novati and Brebbia [21]. The general statement of the integral equations was presented but no numerical examples were given. Kamia and Sawaki [22] gave a BEM formulation for the finite deflection of elastic plates with immovable edge conditions. Such boundary conditions allowed the use of Berger’s equation which is a simplified form of the equation governing the 7 nonlinear deflection of plates [23]. Tosaka and. Miyake [24] presented a BEM formulation for geometrically nonlinear problems of elastic bodies with extremely nonlinear behavior such as the snap-through phenomenon” Two numerical exampleS‘were[consideredz bending of a one dimensional shallow arch, and bending of a shallow spherical shell. In both examples, the nonlinear system of equations obtained from the nonlinear integral equations was solved using the Riks-Wempner method. The earliest BEM application to elastoplasticity was due to Swedlow and Cruse [25]. This was based on an "initial strain" approach. Kumar and Mukherjee [26] used a similar initial strain formulation and presented solutions for simple problems of expanding cylinders and spheres. The first formulation of BEM for 3—D elastoplasticity was due to Banerjee et al. [27]. Their formulation was based on the "initial stress" distribution within the yielded zone. Both boundary and domain elements were chosen to be constant (i.e., constant shape functions) so as to simplify the analysis. The deformation gradients inside the domain were calculated from the nodal displacements as in the FEM in order to avoid singularities. This analysis was improved in another study [28] by the use of isoparametric quadratic boundary elements 8 and a superior numerical integration scheme over the boundary elements and the elastoplastic volume cells. In all the elastoplastic applications discussed above, the numerical solutions required an iterative procedure which was not efficient in some cases. For that reason, Raveendra [29] and Banerjee and Raveendra [30] presented a non-iterative two dimensional elastoplastic analysis which is similar to the variable stiffness method used in FEM. This procedure was extended to the three dimensional case by Henry and Banerjee [31]. There has been very limited work on the applications of BEM to finite elasticity. In a paper by Phan-Thien [32], a constitutive law based on a micro-structural model was derived. The resulting 2-D integral equations involved domain integrals due to the nonlinear part of the constitutive law. The numerical solution of the integral equations was obtained by discretizing the boundary into constant elements (for displacement and traction) and the domain into triangular elements. The solution of the resulting algebraic equations was then obtained by iteration. A computer program based on the above procedure was tested by two examples. In the first example, the problems of homogeneous simple shear and uniaxial extension were considered. The results were satisfying for deformations up to 300% and 185% in the shear and uniaxial 9 extension cases respectively. In the second example, the problem of finite deformation of a circular elastic slice perfectly bonded to two parallel rigid end plates was considered. Phan-Thien extended the above analysis to include the three dimensional case [33]. In both studies, the computational problem of singular domain integrals was avoided by calculating the deformation gradients indirectly from the domain nodal displacements as in FEM. The hydrostatiijressure inside the domain was not determined and therefore, the stress obtained was only a part of the total stress. In.the present study, a boundary element formulation for large elastic deformation of compressible and incompressible plane bodies is presented. The constitutive law is divided into linear (Hooke’s law) and nonlinear parts. A linear solution is first obtained by neglecting the nonlinear part of the constitutive law. The linear solution is then corrected by considering the nonlinear terms through an iterative procedure. Unlike reference [32], the domain deformation gradients are calculated directly from the domain integral equations. Also unlike :reference [32], the hydrostatic pressure is; obtained inside the domain and therefore, the total stress can be calculated. The above formulation is tested using numerical examples for the following cases: 1. plane 2. plane 3. plane 4. plane stress strain stress strain deformation deformation deformation deformation 10 of a compressible body; of a compressible body; of an incompressible body; and of an incompressible body. CHAPTER 2 FORMULATION OF THE BOUNDARY INTEGRAL EQUATIONS 2.1 Linear Elasticity The formulation for the case of linear elasticity is given here in order to introduce the boundary integral method and because the linear solution will be utilized as a part of the nonlinear solution later. 2.1.1 Governing Equations Consider a plane linear elastic, isotropic, homogeneous body that occupies a region Q and is bounded by a boundary P. The equations of equilibrium can be written as: ”11.: +19, =0, er, (2.1) 11 12 where caissthe stress tensor, biis the body force vector, and i=1,2, j=1,2. Summation on repeated indices is assumed and the comma indicates differentiation. The stress tensor is related to the strain tensor through Hooke’s law which can be expressed as a” = 15116“: + 2pc”, (2.2) where 6“ is the Kronecker delta whose properties are _ 0 1*1 511 ‘ {1 1:1: (2-3) A and u are Lame’s constants which can be expressed in terms of the more familiar shear modulus, G, the modulus of elasticity, E; and Poisson’s ratio, v, tar the following formulae = = _______ = vEI " G A (1+v)(1-2v)’ (2‘4) 13 and efi is the strain tensor which is defined by 1 e” = —2- (u1'j+uj'1) , (2.5) where u,is the displacement vector. The boundary conditions considered here are such that at each point on the boundary either a displacement or a traction is specified in each coordinate direction. i.e: u, = u, or gun, = t1 = (:1, (2.6) where niis the unit vector normal to the boundary. 2.1.2 Integral Formulation The weighted residual statement correspondimg to equations (2.1) can be written as ’ = , 2.7 four,“ +191) u, do 0 ( ) 14 where u; is a weight function which will be chosen later. Integrating equation (2.7) by parts, one obtains Equation (2.8) is usually' called the "weak" variational statement and it is the basic equation needed to derive the stiffness equations of the finite element method. Integrating the first domain integral in equation (2.8) by parts again results in foobd u, do - Lo}, 11, ui d[‘ (2.9) +fp011 111 u; dI‘ +[0u1'b1 d0 = 0, where a; is the stress due to uf. Equation (2.9) is sometimes called the "inverse" variational statement and it is the basic equation needed to obtain the boundary element equations. Let us choose the weight function, uf, such that 011,, + A (x.£) e, = o, is infinite plane, (2.10) where eiis a unit load applied in the i direction at f and A 15 is the Dirac delta function which has the following properties Anni) =0. xsfi. (2.11) fof(x) A(x,€) do = HE) . tea. The displacement and traction solutions associated with equations (2.10) can be written in the form C H. H U,, (x. E) e]. (2.12) = T11 (xlfi) ejl where U0. represents the displacement at point x in the 1 direction due to a unit load applied at E in the j direction in infinite space, and I} is the corresponding traction. The expressions for U0. and Ti; can be written as _ 1 u 1 . U11 "’ aflflh[(3—v ) O1jln; + (1+V )Pipj] r 1 u u T1;: = ‘zn—ruu’v )511 + 2(1+v)pipj]pknk (2.13) -(1—v') (pjn, — 91111)}. where r is the distance from S, the point of application of 16 the unit load, to any point x under consideration, i.e. r qua-51)? + (29—5,)2. (244) piis given by P: = _, (2.15) n,is the normal to the boundary, u’is given by v for plane stress, (2 . 16) —; for plane strain, where v is Poisson’s ratio, u is the shear modulus, and h is the body thickness. This solution is usually called the "fundamental" solution. Substituting equations (2.10) and (2.12) into equations (2.9), we obtain u,<£) = — fP Tkj(X:€) ukm arm +frUkj(x,£) 1:,r dI‘(X) (2.17) + [a U,,(X.£) bkm dam. £60. 17 where we have employed equations (2.11) and (2.13). Equations (2.17) express the displacement at any point inside I2 in terms of the boundary displacements, u*, the boundary tractions, tk, the known domain body forces, bk, and the known fundamental solution. Note that only one half of the information on the boundary values of displacements and tractions is available, i.e. specified. The other half of the information can be obtained if we apply equations (2.17) at the boundary, i.e. take £ to P. In doing so, the boundary integrals in equations (2.17) become singular. This singularity and other singularities presented in the finite elasticity formulation are discussed in Chapter 3. The resulting boundary integral equations can be written as a,,(£)u,,(£)+[r T,,(x.£) ukm arm =[P U,,(x.£) txm d[‘(X) +[0U,,,(x.£)b,m dam. teI‘. (2.13) 18 where the integral on the L.H.S is interpreted in the Cauchy principal value sense and the components of aKi depend on the boundary geometry at 5. For a smooth boundary aKi = k6”, where 5H is the Kronecker delta. Once equations (2.18) are solved for the non-specified boundary displacements and tractions, the solution at any point inside the domain can be obtained by means of equations (2.17). Unlike the finite element procedure, the displacement gradients inside the domain can be obtained directly by taking the proper derivatives of equations (2.17). The result is auj(£) __ 623,, aa, — fr ——(x,£) ukun arm 35, an” +fp a—E...’(x’£’ txm arm .fo aU*1(x.§) mm mm. gen. 35. (2.19) 19 2.2 Finite Elasticity 2.2.1 Governing Equations The formulation presented here will be based.on the undeformed configuration (total Lagrangian formulation). The position vector of a particle in the undeformed configuration is denoted by x. After deformation, the particle takes a new position denoted by x. The deformation gradient tensor is defined by - 6x1 F11 — 5171' (2.20) The left and right Cauchy-Green deformation tensors can be found from F0. according to 311 = Fur-'11:! _ (2.21) The physical meaning of these kinematic quantities can be found in any standard text book, (e.g. [34]). 20 If the material is incompressible, then the following constraint must be satisfied detF= 1. (2.22) The equations of equilibrium can be written as 011“, + bi = 0' (2.23) where a“ is the stress tensor per unit undeformed area (known as the Piola-Kirchoff stress) which is, in general, a nonlinear function of the deformation gradients (and a hydrostatic pressure, in the case of an incompressible material). The boundary conditions are the same as described in Section 2.1. Let us assume that a“ can be split into linear and nonlinear parts, i.e 1 011 = 011 + 031, (2.24) where 0% is the stress encountered in linear elasticity which 21 is related to the strain via Hooke’s law and onij is the remaining part. Then, the equilibrium equations become 011.: + ”’1 + 031.1) = 0' (2'25) where the nonlinear part has been added to the body force, resulting in a "fictitious" body force. Note that, in general, 0“ cannot be directly split into "Hookean" and "non-Hookean" parts as inferred by equations (2.24). However, this can be accomplished simply by replacing 0"" by (“ii - a'ij) . 2.2.2 Constitutive Law The constitutive law for a hyperelastic material is derived from a strain energy function, W. The general form of the constitutive law can be written as __ 3W 01:, - “61:71, (2.26) where “fl is the stress per unit undeformed area and Ffi is the 22 deformation gradient tensor. If the material is incompressible, another term involving a hydrostatic pressure due to the incompressibility constraint is added to the R.H.S. of equation (2.26) as follows: BW' -T = — F I . 011 31..” P 11 (2 27) where p is the hydrostatic pressure and F57 is the transpose of the inverse of I}. For an isotropic material, W is usually given as a function of the jprincipal invariants of the Cauchy-Green. deformation tensor , CU I i e e e W=W(Il,Iz,Ia), (2.28) where 11,12 and I3 are the invariants of C“ and are given by 11 = bijcaj' ”H n 1 E (aijaklcijckl _ bijbklcikcfl) ’ (2.29) I3 = det CU . 23 Using equation (2.28) in equation (2.26) and applying the chain rule for derivatives yields the following 6W 3W OIJ=Z[—Fij+-a_i-: BI; 6W '“53*(1}F31"1HF31CE1"IQmCLH%u)]- a ( 11F11"53k‘%u ) (2.30) Equations (2.30) suggest that one should know the form of the strain energy function, w, in order to obtain a stress- displacement relation. Various forms of W have been proposed for specific materials. Although the methodology for the solutions given here is valid for the various forms of w, only the following two forms will be considered for compressible and incompressible materials, respectively. 1. The Blatz-Ko model [35]: 1: W:.E( 2 .I +2¢I§—5), (2.31) U where u is the shear modulus. The constitutive law for this material then becomes 24 I I 3 3 3 2. The Mooney-Rivlin model [36]: W=C1(I1-3) +Cz(Iz-3), (2.33) where Cl and C2 are elastic constants called the Mooney-Rivlin constants. Note that the third strain invariant, I“ does not appear in this model because for an incompressible material, I3== 1. This constraint will be used later to determine the unknown hydrostatic pressure which was introduced in the constitutive law. The constitutive law for a Mooney-Rivlin material becomes 011 = 2 [ (C1 + 0211) F1: ‘ 02311511] ‘ PFu—T' (2'34) 2.2.3 Integral Formulation If the body force in equations (2.17) corresponding to the linear case is replaced by the fictitious body force, then we 25 obtain the following domain integral equations u,(£> = —fr rumt) aux) arm + [P mama) aim dI‘(X) + f‘1 mama) ohm) dam + [a ovum) him dour), 550' (2.35) where Tkj and UKi are the same influence functions derived for the linear elasticity case and tk is the linear part of the traction, which is given by t; = a,“ n1. (2.36) Note ‘that the second. boundary integral on ‘the lR.H.S of equations (2.35) involves only the linear part of the boundary traction, while in a finite elasticity problem the total traction is prescribed. The nonlinear part, however, can be retrieved by applying the divergence theorem to the domain integral in the expression. Equations (2.35) then become 26 u,<£) =—fr T,,(x.£) ukm arm +f U,,(x.£) t,(X) arm ’faa 6155”,“ 07mm dD(X) +1; Huang) bkur) dQIXI. 560- (2.37) Equations (2.37) express the displacement at an internal point, f, in terms of the boundary displacements, the boundary tractions, the body force and a domain integral whose integrand consists of nonlinear functions of the deformation gradients (and a hydrostatic pressure in the case of an incompressible material). By applying equations (2.37) at the boundary, one can obtain the following boundary integral equations a,,(£)u,(£)+fr T,,(x.£) ukm arm = [P U,,(x.£) gm arm -06]. 60” (X15) 072m dO(X) . [a Uk,(x.£) bk(x) doom. 6613 (2.38) Note that the domain integral in equations (2.38) contains 27 additional unknowns compared to the linear elasticity case, equations (2. 18) . These unknowns are the deformation gradients (and a hydrostatic pressure in the case of an incompressible material) which are contained in the expressions for oflm. The displacement gradients can be obtained by differentiating equations (2.37) to obtain au,(£) ____ aTk 35.. P —a—€" (LE) u),(X) dI‘Ur) U' +fa a £10125) tax) are!) an; n Mfr—Kali (X,E) 0),..(1') arm 6U” 6E? (2.39) E) b,,(X) (10 (X) 560. CHAPTER 3 BOUNDARY ELEMENT FORMULATION 3.1 Boundary and Domain Discretization In order to derive a system of algebraic equations, let us divide the boundary into N elements and the domain into M cells as shown in Figure 3.1. In the case of no real body forces (i.e b=0), the boundary integral equations (2.38) become N a,,(£)u,,(£)+ {f r,,(x.£)u,,mdrm} 1;}. I“ = U x, t dl‘ (3.1 gun ,,( z) ,m m} > in: an” (x no” (X) dO(X) EeI‘ 1,1 max, ' k” ' ' Let.us next approximate the displacements and tractions on the boundary by assuming linear variations over each boundary element, i.e. 28 U: = 9. (n) uj‘H’ + 0. (n) u)“, _ (3.2) t1 = ¢1(n) t;21 1) + Q: (n) t;31), over element. i, ‘where 12‘ is the displacement in the j direction at node i, tfim and t?”” are the tractions in the j direction "before" and "after" node 1 respectively, and‘P.and ‘P.are linear interpolating functions which can be written in terms of a local coordinate n as @101) =(—1‘2‘-9-)—. 92m) =—(—1—;fll. 451151. (3.3) where n ranges linearly from a value of -1 at node i-1 to a value of 1 at node i over element 1. Furthermore, let us assume that the fictitious body force is constant over each domain cell and assumes the value at the centroid of the cell. The global coordinate system X is transformed linearly to the local coordinate n using the same linear interpolating functions x = o1(n)x<1-n + «b, (mx‘i’. <3-4> over element 1. Similarly, 29 30 d8(X) = (5(1) -28(1-1))dn .__ 3381“, (3.5) where Asiis the length of element i. Using equations (3.2) through (3.5) in equations (3.1), we obtain ” As ak,(5)u,,(€)+ { 4‘ afifnuml T,,,(X(n).E)dn} =1 +" ”“1 u‘[ (1— ) 1' m ) 0dr: 2: 4 I: I‘ '1 Id 7| I 1+1 4 czi'lfnu-manum)man} (3.6) s { ‘ti‘f (1+n) U,,(x(n).£)dn} 1a. 4 " 60,: -Z ( of...) iafa (x step size I Yes . . No Total boundary conditions apphed? l Yes (Stop) Figure 4.1 Flow chart for the program BLA'FZKO wL CHAPTER 5 INCOMPRESSIBLE MATERIALS In this Chapter, the boundary element solution is obtained for a material obeying the Mooney-Rivlin constitutive law. The solution is obtained using an iterative procedure for solving the nonlinear equations. The implementation of the iterative procedure leads to a Fortran computer program called RIVLIN. 5.1 Plane Stress Using’ the incompressibility constraint. given by (equation (2.22), the deformation gradient tensor for the plane stress case becomes A11. A12 0-! F = 121 A22 0 , (5.1) o o l A. where A is given by 47 43 A = 111122 — 112121. (5.2) Recall that the constitutive law for a Mooney-Rivlin material is given by “11 = 3 [ (C1 + 9211) F1: ‘ czaungj] - pFij‘T. (5.3) The hydrostatic pressure can be expressed in terms of the plane deformation gradients by using the plane stress assumption 033 = 0. (5.4) Substituting equations (5.1) and (5.3) into equation (5.4), we obtain the following equation for the hydrostatic pressure: 2 P=?[C1+Cz(lii+lgz +412+121HI (5-5) where C, and C, are the Mooney-Rivlin constants, and A is given by equation (5.2). Note that we have eliminated.the hydrostatic pressure from the 49 equations. The hydrostatic pressure ,however, can be calculated from the solution of the deformation gradients by the use of equations (5.5). 5.2 Plane Strain The deformation gradient in this case is given by A11 A12 0 F = 121 12, o. (5-6) 0 O 1 The incompressibility constraint becomes 111122 - 112121'-1 = O. (5.7) Note that for small deformations (linear elasticity), the above constraint reduces to the following linearized version of the incompressibility constraint in + 122 — 2 = 0. (5.8) 50 5.3 Iterative Procedure Recall from Chapter 3, that the BEM equations for an incompressible material are given by A {x} = {a} + {flVuml}. (5.9) (det F)1 - 1 = o, 1=1,M. (5.11) Note that equations (5.11) are not necessary for the plane stress case since the hydrostatiijressu e was eliminated from the equations and the incompressibility constraint was satisfied.by the(deformation.gradient.given.in equation (5.1). The solution for the above nonlinear equations can be obtained using the following procedure. 1. Apply an increment of the boundary traction (or displacement) and use equations (5.9) and (5.10) to calculate the initial values for boundary unknowns, (x}fi and the plane displacement. gradients, {uU}H based. on linear' elasticity 51 assumptions (i.e: {f}={g}=0). Thus H 1 1 1 {Va} = B u + D {t}. (5.13) 2. Obtain initial values for the hydrostatic pressure as follows a. For the plane stress case, the initial hydrostatic pressures are calculated by substituting the results for the initial plane displacement gradients, obtained in equations (5.13), into equations (5.5). b. For the plane strain case, the initial values for the hydrostatic pressures can be obtained by satisfying the linearized. incompressibility' constraint, i.e. substituting equations (5.13) into equations (5.8) . The resulting algebraic equations have the following form 9 = qr (5°14) 52 where [Q] is an M by M matrix, and {q} is an M by 1 column matrix. The components of [Q] and {q} are functions of the initial values of the boundary displacements and tractions, and the displacement gradients obtained in step 1. {p} is an M by 1 column matrix containing the values of the unknown initial values of the hydrostatic pressure in the domain cells. 3. Use the results of the displacement gradients from step 1 and the results of the hydrostatic pressure from step 2 to calculate the components of [ff and {g}K This completes the first iteration. 4. Use the results for [ff from step 3 in equations (5.9) to get new values for the boundary unknowns. II p 4. b A H) ...: 5. Use the results for {g}‘from the first iteration and the updated values of the boundary unknowns obtained in step 4 to 53 calculate new values for the displacement gradients, i.e. 2 2 2 1 6. Obtain new values for the hydrostatic pressure by satisfying the nonlinear incompressibility constraint given by equations (5.7) over every domain cell. This can be accomplished by using the updated values for the boundary unknowns and the displacement gradients in equations (5.7) and (5.10). The resulting nonlinear algebraic equations have the following form f,(u , t , Vu , p) = 0, 1=1,M. (5.15) The above nonlinear equations require initial values for the hydrostatic pressure as required by any iterative solution. The results of step 2 can be used for this purpose. 7. Continue iterations until convergence is reached. 8. Add another increment of boundary traction (or displacement) and repeat the above iterative procedure. 9. Add.more increments until the total boundary conditions are applied. 54 5.4 Computer Program RIVLIN The iterative procedure described above is implemented in a computer program, RIVLIN, written in Fortran to solve finite plane deformation problems involving Mooney-Rivlin materials. The approximations utilized for the boundary and domain elements are the same as those assumed in the program BLATZKO. The flow chart for the program is shown in Figure 5.1. 55 lInput Data] Calculate [H] [G]. [B]. & [D] .4 Set lfi=lgi=0 Apply a small increment & calculate the boundary unknowns, 3x2, from eq. (5.9) 1 Calculate the plane disp. gradients _-_-J inside the domain using eq. (5.10)- Plane strain Calculate P from eq (5.15) l , Reduce step size l Plane stress or plane strain Plane stress Calculate p from eq (5.5) q Calculate if] and lg] l Calculate new values for 8c 3vu$ using eqs (5.9) &(5.10) N Converged? 1 Yes Total boundary conditions apphed? Yes > NO Figure 5.1 Flow chart for the program RIVLIN )l CHAPTER 6 NUMERICAL EXAMPLES The two computer programs described in Chapters 4 and 5 for the compressible and incompressible finite plane deformation cases, respectively, are employed in this chapter to obtain numerical solutions for several example problems. 6.1 Incompressible Plane Stress Deformation Plane Stretching of an Elastic Sheet Consider an 8" square rubber sheet, 0.05" thick, subjected to a uniform plane stretching as shown in Figure 6.1. As noted in [13], this problem corresponds to the so-called bi-axial strip test commonly used to characterize ultimate properties of materials such as rubbers and polymers. The sheet is made of a Mooney-Rivlin material with elastic constants Ch and C5 of 24.0 and 1.5 psi, respectively. The boundary conditions are: 56 57 u1(:):4.yl — 14(1-1). u3(t4..v) - 0. t1(x,1:4) — 0. 83(1):,14) — o, (6-1) where A is the extension ratio along the xq direction such that for A=1, the sheet is undeformed. In order to see the effect. of the mesh refinement, the four boundary-domain element models shown in Figure 6.2 were considered. The total edge force required to stretch the sheet to twice its length was obtained for each model as shown in Figure 6.3. The total edge force was calculated by simply summing the products of the x1 components of the tractions on the boundary x1 = 4" times the lengths of the corresponding boundary elements. Figure 6.3 shows that the total edge force converges to approximately 35.9 lb as the boundary-domain mesh was refined. FEM solutions are available for this problem [13,15]. In reference [13], the solution was based on 72 elements and yielded a value for the edge force of 36.0 lbs. The second FEM solution [15] used a 6x6 mesh of 4—nodded quadrilaterals to model quarter of the sheet and yielded a value of 161 N (36.2 lb). The finite elasticity solution is compared to the linear elasticity solution in Figure 6.4 where the total edge force 58 is obtained for different values of the extension ratio, A. Figure 6.4 shows that for A=2, the edge force solution by linear elasticity is about 1.8 times the finite elasticity solution. A comparison of the present results to the FEM results from reference [15] are given in Figure 6.5 where the total edge force is plotted as a function of the extension ratio, A. This figure shows a good agreement between the two solutions. 59 8" L r I'T 1 X2 __ T , A / 2 I é 2 6 8" fl :— / r / ' § x1 / 2 2 4 ? .1 2 / 8 A 7i- 8" Figure 6.1 Uniaxial stretching of an elastic sheet. 6O Figure 6.2 Various models used for the elastic sheet problem. lbs Total edge force, 40 39 38 37 35 35 34 61 l T l l l l C)\ \ \4‘ ~_\\ --‘ \\\\O I L l L l L 4 8 1 2 l 6 2 O 2 4 Number of boundary elements Figure 6.3 Total edge force versus the number of boundary elements 28 lbs Total edge force, 70 50 50 4o 20 62 O T Linear elasticity Finite elasticity l l Figure 6.4 Total edge force versus the extension ratio 1.4 Extension ratio, >\ 1.6 1.8 l\) Toto! edge force, lbs 63 40 r m I I I V O BEM /O V FEM [)5] / 3O — i // 20 r / 1 .63 10 r — DE; 1 l l l l 1.0 12 1.4 1.6 1.8 2.0 Extension rotio, )\ Figure 6.5 Comparision between the BEM and the FEM solutions for the total edge force 64 Uniaxial Stretching of a Sheet with 3 Circular Hole Consider a 6.5" square rubber sheet, 0.079" thick, containing a central 0.5" diameter hole and subjected to a uniform uniaxial stretching as shown in Figure 6.6. The boundary conditions are 13-25 (1‘1): u,(13.25,y) u,(13.25,)d -0. u,(x,13.25) = O, t1(x,13.25) =0, (6-2) where A is the axial extension ratio. In this problem, the sheet is assumed to be made of a Mooney-Rivlin material with elastic constants Cl and C2 of 27.02 and 1.42 psi. FEM solutions for this problem are available [13,15]. Due to the symmetry of the problem, we consider a quarter of the sheet. The boundary was modeled by 24 unequal elements and the domain by 36 cells as shown in Figure 6.7. The results for the required edge force as a function of the extension ratio, A, are given in Figure 6.8 along with a comparison with the linear elasticity solution. For an extension ratio of 3, the linear solution is almost three times the finite elasticity solution. A comparison between the present results and the FEM 65 results [15] for the edge force are given in Figure 6.9. The FEM solution in [15] was based on a 6x6 mesh of 4-nodded quadrilaterals similar to the one used in the present analysis. The deformed.profiles of the initially circular hole are shown in Figure 6.10 for various values of A along with a comparison to the FEM results. 66 6.5" X2 W/fl/fl 7/2 7/ ‘ //7/////////////////////////% x 6.5" / % X1 WWW/22253772127777; 7 7 7/7/72“ /// C/ ////////////// 6.5A Figure 6.6 Uniaxial stretching of a sheet with a central circular hole. X2 Figure 6.7 Model for a quarter of the sheet. X1 Edge for (1 e 68 I I T I T I I I I 220 L _ —— Linear elasticity 200 P a 0—0 Finite elasticity 180 ~ 2 lGO r ‘1 / I40 — _ I30 - _ 100 r — Q C; I 1 l 1 l l l I l l 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 Extension rotio ,.\ Figure 6.8 Total edge force versus the extension ratio lbs Total edge force. 100 90 80 7O 69 l l l l l l l l l '— V O BEM 0 ,/v/ P V FEM [15] W/ ’ /O 4 ' I— "7// //,.O’ ////’ / v”/ H C E," I L 4 I u I I :M, a I .0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.3 :28 3.0 Extension ratio, >\ Figure 6.9 Comparision between the BEM and the FEM results for the total edge force of (lofm‘II‘ch hole boundary Xz—(Jool'd. 70 I I I 1.50 — -— O BEM 1.25 _ V FEM [15] __ 1.00 — __ 0.75 —- -— 0.50 5'? — 0.25 undef. 1-3'0— (1-1.0) \fito \ 0.00 L o» 1 5+ 0.0 0.5 1 O 1.5 XI—coord. of deformed hole boundary Figure 6.10 Deformed shape of quarter of the hole boundary 71 6.2 Incompressible Plane Strain Deformation An Infinite Cylinder under Internal Pressure Consider an infinitely long thick—walled cylinder subjected to an internal pressure, p,, as shown in Figure 6.11. The cylinder is assumed to be made from a Mooney-Rivlin material with elastic constants C1 and C2 of 80 and 20 psi , respectively. This problem has been analyzed by several investigators, e.g. [37,38], due to the availability of its analytical solution [34,36]. The analytical solution as obtained in [34] is given by“ _ R2 tr: ‘1’ + 2C2 + 2(C1+C2) I2, I2 Too "‘1’ + 202 + 2(01-1-02) 123' p = -p1. - zc2 - 2(c1+c,)(1nir— + §(—1———i) - 1nR£ + £2), 1 r2 122 1 r2 122 +1) R R -R2 m =(01 + C2) (1n : - 2111—1 +b °= ‘ no +b Ra (123 +1») (1234-12) 13 = 212vr + urz, 'r=R+ur'. 72 In these equations Ri and R0 are the inner and outer radii, respectively, R and r are the undeformed and deformed radial distances to the point of interest, respectivelyg mg is the displacement in the radial direction, p is the hydrostatic pressure due to the incompressibility constraint, 5% is the internal pressure,1rr is the radial Cauchy stress, 199 is the circumferential Cauchy stress, and C5 and Cg are the Mooney- Rivlin elastic constants. The BEM solution is obtained based on the axisymmetric model shown in Figure 6.12. The results for the internal pressure, p,, versus the radial displacement of an interior node are given in Figure 6.13 along with the exact solution. The results for the radial displacement.profile for various values of internal pressure are given in Figure 6.14. This figure shows that the BEM and exact solutions are the same for an internal pressure of 42.5 psi. For an internal pressure of 131 psi the BEM solution is about 5% higher than the exact solution. The results for the average hydrostatic pressure profile are given in Figure 6.15. The stresses obtained by the present analysis are based on the undeformed configuration "Piola-Kirchoff stresses" while the stresses given in equations (6.3) are based on the deformed configuration "Cauchy stresses". The two types of stresses are related by the following equation 73 where J=Iy=1 for an incompressible material, I} is the deformation gradient and a". and 19- are the Piola-Kirchoff and Cauchy stress tensors, respectively. Since the deformation gradients were obtained inside the domain as part of the iterative solution, we can use equation (6.4) along with the stress transformation procedure to obtain the radial and circumferential Cauchy stresses in the cylinder. The results are given in Figures 6.16 and 6.17, respectively. X2 74 Pi Figure 6.11 Infinite cylinder under uniform internal pressure X2 75 Figure 6.12 Axisymmetric model for the infinite cylinder psi Internal pressure , 76 160 I I I I I ' I I _ Exact 140 — — O BEM o 120 — i — 100— O a 20— _ Radial displacement of an interior node, in Figure 6.13 Applied internal pressure versus the interior node displacement Radial displacement (in) .[x 77 I‘\ —— EXACT O BEM Pi=131.1 psi Pi=42.5 psi l l l L l l l l l l l T T T T 1'7 I I I T I T I 7 8 9 10 11 12 13 14 15 16 17 18 Undeformed radial distance (in) Figure 6.14 Displacement profile for various internal pressure values Average hydrostatic pressure, psi 240 230 220 200 190 180 78 I I I 1 I I 8 10 12 14 16 18 Undeformed radial distance. inches Figure 6.15 Average hydrostatic pressure profile —Radial stress, psi 79 —— Exact 120 - 100 — 80—- Undeformed radial] distance, inches Figure 6.16 Radial stresses in the cylinder Circumferential stress, psi 600 550 500 450 400 350 300 250 200 150 100 50 80 O BEM Exact Undeformed radial distance, inches Figure 6.17 Circumferential stresses in the cylinder 81 6.3 Compressible Plane Stress Deformation Plane Stretching of an Elastic Sheet Consider an 8" square sheet, 0.06" thick, subjected to a uniform plane stretching. The sheet is made of a Blatz—Ko material with a shear modulus of 40 psi. The same problem was discussed in Section 6.1 but the sheet was assumed to be incompressible. A FEM solution is available for this problem [14]. In reference [14], the solution was based on 6x6 mesh of 4-nodded quadrilaterals to model a quarter of the sheet as shown in Figure 6.18. The same model was used in the present study. The finite elasticity solution is compared to the linear elasticity solution in Figure 6.19 where the total edge force is obtained for different values of the extension ratio, A. A comparison between the present results and the FEM results from reference [14] are given in Figure 6.20. 82 X2 1? X1 Figure 6.18 Model for a quarter of the sheet lbs Total edge force, 35 30 25 83 —— Linear Elasticity Q Finite elasticity Extension ratio, A Figure 6.19 Total edge force versus the 1.8 2.0 extension ratio 84 l I l l 1 l l l l V 10 _ V/v/V/ ‘ 8 — / _ /o 2 ":3 6 - / - 3 o 3:0 Q) 4 _ _ ,2 :3 /’) o BEM / FEM [14] 0 / L I I I l I I I I 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 Extension ratio, A Figure 6.20 Comparision between BEM and FEM solutions for the total edge force. 85 6.4 Compressible Plane Strain Deformation Simple Uniaxial Tension Consider simple uniaxial tension of a Blatz-Ko sheet as shown in Figure 6.21. The exact solution can be obtained for this problem as follows. The deformation gradient is given by where A22 can be obtained as a function of A” by using the equation This gives the following equation for An ._ -1/3 Use equations (6.5) and (6.7) in equation (2.32) to get the following 011 = '1 (xii/3 .. 1;: The boundary-domain mesh used to model the sheet consists of 86 4 linear boundary elements to represent the four sides of the sheet and 1 domain cell to represent the area. The BEM results for a" and An are given in Figures 6.22 and 6.23 as functions of A" along with the exact solutions. C1) 87 :3 7“ Av-v—v \ . \.\\.L 3A.\u\ ‘.“.‘: v \ . .XL u.»n;k.. T“. l §\\\V\I~.\..\i\_l :2 \\'\Z\‘\.\I\'I\I>§I . EX Figure 6.21 *.I/ #I‘ , "/1////////////////////Z Simple extension of an elastic sheet. 88 I I I I O BEM 1.0C\ h Exact _ M 0.9 — — 0.8 — 7) I I I I 1.0 1.2 1.4 1.6 1.8 2.0 MI Figure 6.22 Extension along X302.) versus extension along X,()\,). 89 0,, (p31) 4:. I \O\ O BEM Exact 1.0 1.2 1.4 1.6 1.8 Figure 6.23 011 versus Au. BIBLIOGRAPHY 1. Key, S. N., A variational principle for incompressible and nearly incompressible anisotropic elasticity, International Journal of Solids Structures 5, 951-964 (1969). 2. Fried, 1., Finite element analysis of incompressible material by residual energy balance, International Journal of Solids and Structures 10, 993-1002 (1974). 3. Malkus, D. S., A finite element displacement model valid for any value of the incompressibility, International Journal of Solids and Structures 12, 731—738 (1976). 4. Danielson, D. A., Biomechanics Symposium, (Fung and Brighton, ed.), Vol 2, New York (1973). 5. Atluri, S. N., and Kobayashi, A., Numerical modelling of nonlinear behavior of soft. biological materials, Applied Mathematical Modelling 3, 228-231 (1979). 6. Herrmann, L. R., Elasticity equations for incompressible and nearly incompressible materials by variational theorem, 9O 91 AIAA Journal 3, 1896-1900, (1965). 7. Becker, E. B., A numerical solution of a class of problems of finite elastic deformation, Ph.D. Dissertation, University of California, Berkeley (1966). 8. Oden, J. T., Numerical formulation of nonlinear elasticity problems, Journal of Structural Engineering 93, 235-255 (1967). 9. Peterson, F. E. , Campbell, D. M. and Hermann, L. R. , Nonlinear plane stress analysis applicable to solid propellant rocket grains, Bulletin of 5th Meeting of Interagency Chemical Rocket Propellant Group-Working Group on Mechanical Behavior, CPIA Publication 119, vol 1, pp 421-455 (1966). 10. Oden, J. T., and Sato, T., Finite strains and displacements of elastic membranes by the finite element method, International Journal of Solids and Structures 1, 471-488 (1967) . 11. Oden, J. T., and Kubitza, W., Numerical analysis of Nonlinear pneumatic structures, Proceedings of International Colloquium on Pneumatic Structures, 82-107 (1967). 92 12. Oden, J. T., and Key, B., Analysis of finite deformation of elastic solids by the finite element method, Proceedings of the IUTAM Colloquium on High Speed Computing for Elastic Structures, Liege (1971). 13. Oden, J. T., Finite elements of nonlinear continua, McGraw-Hill, New York (1972). 14. Murakawa, B., and Atluri, 8., Finite elasticity solutions using hybrid elements based on a complementary energy principle, Journal of Applied Mechanics 45, 539-546 (1978). 15. Murakawa, B., and Atluri, 8., Finite elasticity solutions using hybrid elements based on a complementary energy principle, Part II: Incompressible :materials, Journal of Applied Mechanics 46, 71-77 (1978). 16. Hughes, T. J. R., Equivalence of finite elements for nearly incompressible elasticity, Journal of Applied Mechanics 44, 181-183 (1977). 17. Hughes, T. J. R., and Malkus, D. 8., 0n the equivalence of mixed finite element methods with reduced/selective integration displacement methods, Proceedings of the Symposium on Application of Computer Methods in Engineering, Los Angeles, California, Volume 1, 23-32 (1977). 93 18. Noor, A. K., and Andersen, C. M., Mixed models and reduced/selective integration displacement models for nonlinear shell analysis, International Journal for Numerical IMethods in Engineering 18, 1429-1454 (1982). 19. Argyris, J. H., Dunne, P. C., and Muller, M., Isochoric constant strain finite elements, Computer Methods in Applied ‘Mechanics and Engineering 13, 245-278 (1979). 20. Chang, T. Y. P., Saleeb, A. F., and G. Li, Large strain analysis of rubber-like materials based on a perturbed Lagrangian variational principle, Computational Mechanics 8, 221-233 (1991). 21. Novati, G., and Brebbia, C., Boundary element formulation for’ geometrically' nonlinear elastostatics, Applied Mathematical Modelling 6, 136-138 (1982). 22. Kamia, N., and Sawaki, Y., An integral equation approach to finite deflection of elastic plates, International Journal of Nonlinear mechanics 3, 187-194 (1982). 23. Berger, H., A new approach to the analysis of large deflections of jplates, Journal of Applied Mechanics 22, 465-472 (1955). 94 24. Tosaka, N., and Miyake, 8., Integral equation analysis for geometrically nonlinear problems of elastic bodies, Proceedings of lst Japan-China Symposium on Boundary Element Method, 251-260 (1989). 25. Swedlow, J., and Cruse, T., Formulation of boundary integral equations for three dimensional elastoplastic flow, International Journal of Solids and Structures 7, 1673-1683 (1971). 26. Kumar, V., and Mukherjee, S., A boundary integral equation formulation for time dependent inelastic deformation in metals, International Journal of Mechanical Sciences 12, 713-724 (1975). 27. Banerjee, P., Cathie, D., and Davies, T., Two and three-dimensional problems of elastoplasticity, Developments iJiBoundary Element Methods-1, 65-95, Elsevier Applied Science Publishers, London (1979). 28. Banerjee. P., and Davies, T., Advanced implementation of boundary element methods for three-dimensional problems of elastoplasticity and viscoplasticity, Developments in Boundary Element Methods-3, 1-26, Elsevier Applied Science Publishers, London (1984). 95 29. Raveendra, 8., Advanced development of BEM for two and three dimensional nonlinear analysis, Ph.D. Dissertation, State University of New York at Buffalo (1984). 30. Banerjee, P., and Raveendra, S., A new boundary element formulation for two-dimensional elastoplastic analysis, Journal of Engineering Mechanics 2, 252-265 (1987). 31. Henry, D., and Banerjee, P., A thermoplastic BEM analysis for substructured axisymmetric bodies, Journal of Engineering Mechanics 12, 1880-1900 (1987). 32. Phan-Thien, N. , Rubber-like elasticity by boundary element method: finite deformation of a circular slice, Rheological Acta 27, 230-240 (1988). 33. Phan-Thien, N., Boundary element method for finite elasticity, Computational Mechanics 6, 205-219 (1990). 34. Green, A. E., and Zerna, W., Theoretical elasticity, 2nd edition, Oxford University Press, London (1968). 35. Blatz, P., and Ko, W. L., Application of finite elasticity theory to the deformation of rubber, Transactions of the Society of Rheology, Vol. VI, 223-251 (1962). 96 36. Rivlin, R. S., and Saunders, D. W., Large elastic deformations of isotropic materials VII. Experiments on the deformation of rubber, Philosophical Transactions of the Royal Society, vol. A243, 251-288 (1951). 37. Scharnhorst, T., and. Pian, T., Finite element analysis of rubber-like materials by a mixed model, International Journal for Numerical Methods in Engineering 12, 665-676 (1978). 38. Sussman, T., and Bathe, K., A finite element formulation for nonlinear incompressible elastic and inelastic analysis, Computer & Structures 26, 357-409 (1987). ‘IIIIIIIIIIIIIIIIIII