THESIS i . . 4-7 ‘ .v f} ". “h - ."-.'. ‘ ... - é . 5 ‘3'. ,‘.. ‘5‘“. ,IR Y i -t§"1 S" : .icricirn333 if; k ' q I lea :‘?r‘- v E at} “J‘- V' Lami- ’ 3' ‘WMW'W +l This is to certify that the thesis entitled IMPROVEMENT OF THE BOUNDARY INTEGRAL METHOD FOR PLANE ELASTOSTATICS USING COMBINATIONS OF SING- ULARITIES AND ANALYTIC INTEGRATION presented by ENAYAT . H . MAHAJERIN has been accepted towards fulfillment of the requirements for ~Pi+.—B.——— degree in _Mechanics_ Date F? g/ 0-7 639 MSU LIBRARIES “ RETURNING MATERIALS: Place in book drop to remove this checkout from your record. FINES will be charged if book is returned after the date stamped below. IMPROVEMENT OF THE BOUNDARY INTEGRAL METHOD FOR PLANE ELASTOSTATICS USING COMBINATIONS OF SINGULARITIES AND ANALYTIC INTEGRATION BY Enayat H. Mahajerin A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Metallurgy, Mechanics and Materials Science 1981 hi ABSTRACT IMPROVEMENT OF THE BOUNDARY INTEGRAL METHOD FOR PLANE ELASTOSTATICS USING COMBINATIONS OF SINGULARITIES AND ANALYTIC INTEGRATION BY Enayat H. Mahajerin Starting with a general formulation for the Boundary- Integral Equation (BIE) method in plane elastostatics, a wide range of options are investigated. To avoid individu- al programming, a general influence function is replaced by a suitable linear combination of fundamental components of singular solutions which are classified in appr0priate forms. Accordingly, a useful analytic scheme is applied to both straight and curved boundaries, wherein the crude discretization is modified using a new technique which leads to highly accurate results. It is also shown that the components of the so called "indirect" approach always lead to a closed-form solution which makes the technique both efficient and independent of numerical integration. The method requires far less operations compared with the most efficient numerical techniques such as the Romberg-Richardson method. The idea is also extended to the case of a variable-singularity-density. Enayat H. Mahajerin Based on the general scheme derived, a compact com- puter program which can be applied to any arbitrary in- fluence function, is developed which requires minimum Storage and computation, The program is self sufficient in that it uses no computer library subroutines. ACKNOWLEDGEMENTS I wish to express my deepest appreciation to my major advisor, Professor Nicholas J. Altiero, under whose in- spiration, constant supervision and unfailing interest this study was made possible. Grateful appreciation is exPressed to the members of my doctoral committee, Professors George E. Mase, Larry J. Segerlind and David Yen who have contributed to the aca- demic background necessary for the preparation of this dissertation. Thanks are extended to Professor David L. Sikarskie, Chairman of the Department of Metallurgy, Mechanics and Materials Science and Gary Burgess my fellow graduate student and friend with whom many profitable discussions were held. I also acknowledge my indebtedness to Ms. Michelle Ward for the care and patience she displayed in typing this manuscript. Finally, special thanks is due to my family, especially to my wife Tahera, for their patience and moral support. ii LIST OF TABLES . . . LIST OF FIGURES . . LIST OF APPENDICES . LIST OF SYMBOLS . . CHAPTER I CHAPTER II CHAPTER III TABLE OF CONTENTS INTRODUCTION . . . . . . . . . . . BACKGROUND AND PRELIMINARIES . . . II.1 II.2 II.3 II.4 THE BASIC EQUATIONS OF ELASTOSTATICS . . . . . . . FUNDAMENTAL ELASTOSTATICS PROBLEMS O O O O O O O O O O PLANE ELASTOSTATICS PROBLEMS BOUNDARY INTEGRAL EQUATION METHOD 0 O O O O O O C O O O II.4.l THE DIRECT APPROACH II.4.2 AN INDIRECT APPROACH II.4.3 GENERALIZATION OF THE INDIRECT APPROACH . II.4.4 REMARKS . . . . . . NUMERICAL SOLUTION OF THE INTEGRAL EQUATIONS . . . . . . . . . . . . III.l III.2 III.3 INTRODUCTION . . . . . . . THE STRUCTURE OF THE GENERAL INFLUENCE FUNCTION . . . . SIMPLIFICATION OF THE INTEGRAL EQUATION FOR NUMERICAL TREATMENT . . . . . . . . . iii Page vi viii ix 10 10 12 14 16 18 18 19 22 CHAPTER IV CHAPTER V III.4 APPROXIMATE INTEGRATION OVER A DISCRETIZED BOUNDARY . . . . III.5 EXACT INTEGRATION OVER THE III.6 DISCRETIZED BOUNDARY . . . . . III.5.1 INFLUENCE FUNCTION OF THE FIRST GROUP... III.5.2 INFLUENCE FUNCTION OF THE SECOND GROUP.. III.5.3 INFLUENCE FUNCTION OF THE THIRD GROUP. . NUMERICAL RESULTS . . . . . . BOUNDARIES OF PIECEWISE-CONSTANT CURVATURE . . . . . . . . . . . . . . IV.l IV.2 SOURC FORMULATION OF THE METHOD . . . (l) 0,0 ' ' IV.l.2 INFLUENCE FUNCTION OF THE FIRST GROUP . . . . IV.l.l EVALUATION OF P IV.1.3 INFLUENCE FUNCTION OF THE SECOND GROUP . . . IV.l.4 INFLUENCE FUNCTION OF THE THIRD GROUP . . . . NUMERICAL RESULTS . . . . . . . IV.2.1 EXAMPLES . . . . . . . E-DENSITY VARIATION OVER A BOUNDARY SUBDIVISION . . . . . . . . V.l LINEAR VARIATION OF DENSITY OVER A STRAIGHT BOUNDARY . . . . . . V.l.l INFLUENCE FUNCTION OF THE FIRST GROUP . . . . . V.l.2 INFLUENCE FUNCTION OF THE SECOND GROUP . . . . VARIABLE DENSITY OVER A CURVED BOUNDARY . . . . . . . . . . . . iv Page 24 27 32 33 33 35 44 45 56 57 64 66 68 70 80 80 81 83 84 Page CHAPTER VI CLOSURE . . . . . . . . . . . . . . . . 88 APPENDICES o o o o o o o o o o o o o o o o o o o o o 90 REFERENCES . . . . . . . . . . . . . . . . . . . . . 107 Table 3.1 3.3 3.4 3.5 4.1 LI ST OF TABLES Results for Example Problem 3.1 Case (1) Loading, u=l.,v=.25 Plane Strain Singularity Applied: (uR) (HR) a.y' a8.y Results for Example Problem 3.1 Case (1) Loading, u=l.,v=.25 Plane Strain Singularity Applied: (uc)a B . . . . . . . . Results for Example Problem 3.1 Case (1) Loading, u=l., v=.25 Plane Strain Singularity Applied: (ug)a BE (uc)a B + 2p :3 e V5 YB (ur)0‘.8 . . . . . . . . . . . . Results for Example Problem 3.1 Case (2) Loading with n=l. Singularity Applied: (HF)B E (HR)B with v=.5 . . . . . . . . .°Y . . . .°Y . Results for Example Problem Case (3) Loading with u=l., v=.25 Plane Strain Singularity Applied: (uR)a B and (HR)a B . . Results for Example Problem 3.2 Loading: Traction Corresponding to the Stress Functiin ¢(x1,x2) = x1 x2 - x12x23 Singularity Applied: (HF) 5 (HR) with v=.5, u=l. . . . . ' . Comparing the Cost and Accuracy of the Present Method with Romberg-Richardson Method (NM)B Y: (HR)B y, (HF)B Y, n=16, u=l . Results for Example Problemm4.l (NM) : (HR) , (HF) Case (1) with n=60§'5=1 . .Bzy. . . 8'? . . . . . . . . Results for Example Problem 4.2 (NM)B Y: (HR)B Y’ (HF)B Y' Case (1) with n=60,'u=1 . . I . . . .°. . . . . . . . . . . vi Page 38 39 4O 41 42 43 69 72 74 Table Page 4.4 Results for Example Problem 4.2 (NM)3.Y: (HR) .Y, (HF)B,Y, Case (2) with n=60, u=§ . . . . . . . . . . . . . . . 75 4.5 Results for Example Problem 4.3 (NM)B.y‘ (HR)B.Y, (HF)B.Y' Case (1) with n=60, u—l . . . . . . . . . . . . . . . 77 4.6 Results for Example Problem 4.3 égfisaléo,(Eiie'l'.(?F§B:*I Iaie.(?). . . 78 4.7 Results for Example Problem 4.3 (NM)B Y: (HR)B.Y, (HF)B.Y' Case (3) 79 with fi=6o, u=l . . . . . . . . . . . . . . . vii LIST OF FIGURES Figure 2.1 Elastic Boundary Value Problem . . . . . . 2.2 Region D Embedded in the Infinite Plane. . 3.1 Geometry and Coordinate System for an Arbitrary Region D . . . . . . . . . . . . 3.2 Region D of Figure 3.1 After Discretization 3.3 Geometry and Coordinate System for Exact Integration Over a Straight Segment . . . 4.1 An Irregular Region with Curved Boundary Embedded in the Infinite Plane . . . . . . 4.2 Boundary Field Point 5 and Source Point x’ 4.3 Geometry and Coordinate System Related to Formulation . . . . . . . . . . . . . . . 4.4 Concave Boundary Segment KC=+1 . . . . . . 4.5 Convex Boundary Segment KC=-1 . . . . . 4.6 The Case w # 6’ 4.7 The Case w = e’ 4.8 The Singularity 4.9 Example Problem 4.10 Example Problem Case (9 4.1 . 4.2 . 4.11 Example Problem 4.3 . 5.1 Linear Variation of Straight Segment . . . O O C O O O Source Density Over a 5.2 Linear Variation of Source Density Over a Curved Segment viii Page 13 20 26 29 46 47 49 52 52 60 61 63 71 73 76 82 85 Appendix A LIST OF APPENDICES Page Derivation of Some Important Influence Functions . . . . . . . 90 The Romberg-Richardson Method Applied to Boundary Integrals . . ... . . . . . . 93 Computer Programs ix LIST OF SYMBOLS Constants Linear combination matrix Quadrative weights Constants Boundary of region D A boundary segment Constants Edge dislocation Vector generated from dipole of CY Elastic constants Constant Region of interest Submatrix of the coefficient matrix of linear system equations Young's Modulus of Elasticity Massonet's singularity Vector generated by dipole of FY Functions Rieder's singularity Vector generated by dipole of GY Collective denotation for integrals Mesh length for straight segments Hm m,n I(k) (NM)B.Y mus. mm B .Y.(um)8 v,vlthrough v6 W through W 1 5 ’ Ix Ix I X K: Collective denotation for integrals Collective denotation for integrals Convexity factor Collective denotation for any arbitrary singularity (Ry,Cy,.....etc.) Collective denotation for vector gen- erated by dipoles of MY (rY,Cy,...etc.) Total number of subdivisions on the boundary Components of the normal at the boundarv field point and sourcegxfixuzrespectively Collective denotation for elastic fields Influence function corresponding to N and MY 8 Nonintegral term corresponding to (NM)B Collective denotation for integrals Constants Collective denotation for integrals Distance between field point and source point Body force singularity Vector generated by dipole of RY Components of displacement Displacements DB generated by MY and mY respectively Functions of Poisson's ratio Fundamental boundary integrals Field point and sourcegxflanrespectively Relative coordinates xi Y (X0 1Y0) z (m) E-m' —k vluyE |m 8' 08 (HM)B.y’(HM)aB,Y Center of curvature Integration points Radius of curvature Poisson's ratio the shear modulus and Young' 3 mo ulus respectively Internal field point Airy stress function Tractions and stresses respectively Influence function corresponding to H8’ H08 due to MY respectively Kronecker delta (611:522=1,512=621=0) Permutation tensor (e =0,e =1) -8 11:822 12= 21 Differentiation with respect to field point and source point respectively Cauchy Principal Value (CPV) xii CHAPTER I INTRODUCTION Many problems of practical interest have been solved by boundary integral equation (B.I.E.) methods. The simplicity of these methods, their numerical efficiency, and their extreme accuracy make them very appealing for application to a wide range of engineering problems. The groundwork for integral-equation methods was es- tablished in the classical work of Fredholm [l], Betti [2], Somigliana [3], Lauricella [4], Muskhelishvili [5], Mikhlin [6], Miche [7], Weinel [8], Kupardze [9] and others. With the advent of high-speed digital computers, these methods have taken on great practical Significance. In general, B.I.E. methods can be classified as "direct" as formulated by Rizzo [12] and Cruse, [13], or "indirect", as established by Massonet [11], Altiero and Sikarskie [14], Banerjee [21], Hiese [l9], Altiero and Gavazza [15] and Others. The idea of the "direct" approach is based on the application of Green's theorem, converting an integral over the entire domain into one over the boundary of that domain. One must then solve a system of boundary equations, the un- knowns of which are real boundary values, i.e. tractions and displacements. For traction boundary-value problems, these equations are singular. For mixed boundary value problems, some of the equations are non-singular. Once the entire set of boundary values is known, the field solution is obtained by integration. On the other hand, the "indirect" approach is based on the principle of superposition, in which the region of interest is embedded in another region of the same material for which the various influence functions are known. Typically the region into which the embedding is done is taken to be infinite, but this need not be so. A layer of singularities e.g. body forces, dislocations, dipoles, etc., is then applied to the embedded boundary and the elastic fields (stresses and displacements) can be eXpressed as the integrated influence of the singularity layers. The unknown layer is obtained by enforcing the boundary condi- tions at n boundary points. Since there are two boundary conditions associated with each boundary point, components of the layers of singularities are solutions of Zn linear equations. The characteristics of the integral equations generated by the "indirect" method depend upon the type(s) of singularities used. For example, the use of a body force layer leads to singular equations for traction boun- dary-value problems. For mixed problems, however, some of the equations are non-singular as with the "direct" method. The use of dislocation dipoles, say, leads to singular equations for displacement problems. Other possibilities exist which will be discussed in this dissertation. Although both the "direct" and the "indirect" methods have been compared to one another in terms of numerical efficiency and accuracy, which of them is the "best" for- mulation has yet to be established. It has been suggested by Hiese [24] and Altiero and Gavazza [15] that singular integral equations may be preferable to non-singular ones. In this dissertation, a general formulation is estab- lished which includes a wide range of integral-equation methods. In Chapter II, the aforementioned techniques are reviewed and generalized. Some restrictions are imposed on various formulations and the groundwork is laid for de- velopment of a suitable scheme for both numerical and analytical treatment. In Chapter III, using some funda— mental relationships, a general scheme for straight boundary segments is derived which avoids individual programming. In Chapter IV, closed-form solutions of the boundary integrals of the indirect type are derived for curved boundary segments, thus avoiding crude discretization techniques. A variety of problems with curved boundaries are examined and highly accurate results are obtained using less Operations than are required by the method of Romberg-Richardson. In Chapter V, analytical eXpressions are obtained for curved boundary segments subjected to variable singularity den- sities. In Chapter VI the closure and conclusions are pre- sented. Finally, in the appendices, a general computer program which includes any arbitrary choice of influence functions of elastostatics for both straight and curved boundaries is presented. CHAPTER II BACKGROUND AND PRELIMINARIES In this chapter, the fundamental equations of linear elastostatics are presented so that the notation to be used in this dissertation is made clear. Also, the various boun- dary integral equation methods are reviewed and presented in such a way to make clear their similarity in form. II.1 THE BASIC EQUATIONS OF ELASTOSTATICS The equations of static equilibrium can be written as: H + b = 0 (0,8 = 1,2,3) (2.1) where 1% and bu are the components of stress and body force 8 respectively. The comma notation denotes partial differen- tiation with respect to the coordinate variable whose sub- script follows the comma, and summation is implied over 1 to 3 when an index is repeated. The stresses Na and the 8 strains ea are related by: 8 H = , a8 CaBY5 6Y5 (2 2) where the elastic moduli Casyd have the symmetry properties: CGBYd = CYOGB = CBayd = Cdde . (2.3) K]; For an isotropic medium: Cm8Y5 = A 608 6Y5 + u (Gav 585 + Gas 68y) (2.4) where A, u are the Lame'Constants and 608 is the Kronecker delta. The strains are given by: _ 1 £08 — 2 (ua,8 + uB,a) (2.5) where ua are the components of the elastic displacement. The components of the boundary traction H are given by: B = = : (2'6) Ha HaB nB CaByd 8Y6 nS CoByd uy,6 nB where nB are the components of a unit normal directed out- ward from the surface of the boundary. II.2 FUNDAMENTAL ELASTOSTATICS PROBLEMS In the absence of body force, the equations of equili- brium are: C = 0 . (2.7) u aBY5 7.58 To solve (2.7) in a domain D, one considers the following types of boundary conditions: (i) ua is prescribed everywhere on the boundary of D, ua = 5a (the geometric problem); (ii) Ha is prescribed everywhere on the boundary of D, Ha = Na (the statical problem): (iii) ua is prescribed on Bu' Ha is prescribed on Bt' where B = Bu + B is the boundary of the region t D (the mixed problem). Figure 2.1 Elastic Boundary value problem. :w-fi .uvv 9'... ,. Osobu v... 1 ‘00-! up . u»... -A- ...C. . ‘LA LA.: I .'V ‘-a§ 1 Pa ‘5‘ a" “P 5- ...~ u.“ I... v ‘- )4 d. For the geometric problem, (2.5) may be substituted into (2.2) and the result, in turn, substituted into (2.1). This leads to the following equations: '* (X‘+ u)‘u + b = 0 (2.8) u u 8.8a a 0.88 which are well known as the "Navier-Cauchy" equations. The solution of this problem is given in the form of a displace- ment vector ua satisfying (2.8) throughout D and fulfilling the boundary condition ua = Ea on B. If e is given ex- a8 plicitly an; a function of the coordinates, then (2.5) im— plies six independent equations for the determination of three displacements. Thus, the system is over-determined and does not have a solution for an arbitrary choice of EaB' Therefore, if the displacement components are to be single- valued and continuous, the compatibility equations: + _ - = 2. Ea8.76 €75.a8 an,BG €66.av 0 ( 9) must be satisfied. There are eighty one equations in (2.9), six of which are distinct. For the statical problem, (2.9) may be combined with Hooke's Law and the equilibrium equations (2.7) to produce the governing equations: H + —l— H + b + b aBrYY l+v yy.a8 a,8 8.0 (2.10) v — + l-v 6a8 prY — 0 which are well known as the "Beltrami-Mitchell"equations of compatibility. The solution of this problem is obtained by specifying the stress tensor which satisfies (2.10) through— out D and fulfilling the boundary conditions Ha = fa on B. For the mixed problem, the system of equations (2.2), (2.5) and (2.1) must be solved. The solution gives stress and displacement components throughout D. The stress com- ponents must satisfy 1% = Na over B and the displacements t! must satisfy the boundary conditions 1% = Ed on Bu- II.3 PLANE ELASTOSTATICS PROBLEMS Many problems in elasticity may be treated with the two dimensional theories known as plane stress and plane strain. In plane stress, the geometry of the body is essentially that of plate with one dimension much smaller than the other two and 128 = 1&8 (x1, x2) a,B = 1,2 and (2.11) II31 = H32 = H33 = 0- Accordingly: _ 1+v _ 2 £08 — E I28 E Gas INY' (2.12) For plane strain problems, the geometry of the body is that of aprismatic cylinder with one dimension much larger than the other two and ua = ua (x1, x2), u3 = 0, a = 1,2. (2.13) Accordingly: = 16 H v (II11 + H22). (2.14) H e + 2 a8 a8 YY U a (18' 33 = 10 The compatibility condition for plane problems reduces to: an e85 EGBIY5 = O, a,B,Y,6 = 1,2 (2.15) where e11 = e22 = 0, 1. e12 = ‘ezi = In the absence of body force, the solution of a plane problem is obtained from the biharmonic equation: v4 o = o, (2.16) and consequently: H = e e o, (2.17) dB aY By Y5 where o is the Airy's stress function. II.4 BOUNDARY-INTEGRAL EQUATION METHOD Boundary-Integral Equation (BIE) methods are the ex- tension of the works of Betti [2 ],Somigliana [3 ],Fredholm [l] and Lauricella [4 ]. During the past two decades, there has been rapid development in the solution of engineering problems using integral-equation techniques. There are a number of such techniques and, in general, these can be classified as either "direct" or "indirect". The best known of these are presented here as background. II.4.1 THE DIRECT APPROACH Rizzo [12] employed Somigliana's identity: 0: u (g) = {Humms (935’) HB(§’)— (IIR)OL.B(§_,§’)UIB (>_<_’)ldS: 2.18 58D ( ) 11 where (uR)a 8 (3,5’) denotes the a component of the dis- placement at 5 due to a unit force in the 8 direction at x’ in an infinite medium, and (HR)a (§,x’) are the 8 tractions which correspond to the displacements (uR)a (F. ’5’) . In three dimensions ds denotes an element of area of B, and in two dimensions, an element of arc length of B. If g approaches a point x on B from inside D, then (2.18) becomes: é—ua (x) +Zz:(HR)a.B (1935’) u8(§’) ds = f(UR)a.B (gay) ”3 (5') as B (2.19) where the integral on the left-hand side is to be inter- preted in the Cauchy Principal value (CPV) sense. For the geometric problem where the displacement is prescribed everywhere, the left side of (2.19) is known and the re- sulting integral equations for Ha are of the first-kind (not singular). On the other hand if the boundary conditions are of type (ii), i.e. traction is prescribed everywhere,1flua right-hand side of (2.19) is known, and one solves a system of singular-integral equations of the 2nd kind for Uh. This arises because (uR)“ is 0 (Zn le) in two dimensions and B is O (Igl-l) in three dimensions, and these singularities are integrable. However (HR)a is O (Ifil-l) in tWO dimen- B sions and is O (le-Z) in three dimensions and consequently the resulting integral must be considered as a CauchY F"__TF1 12 Principal Value integral. For mixed boundary-value problems, only some of the integrals are singular. In any case, after all boundary values are obtained, the displacement at any point 5 2D can be obtained from (2.18). II.4.2 AN INDIRECT APPROACH Benjumea and Sikarskie [20] , Altiero and Sikarskie [14] and Banerjee [21] extended the work of Massonet [11]anuideve10ped a different form of the boundary integral method. Massonet solved the traction problem by embedding the real body in a series of fictitious half planes, sequentially tangent totflma real body. The integral equation solution of the traction problem then has the following form _ L . , HB(§) — 2 68.7FY (x) +‘.IgL(HF)B_Y (§.§) I“Y (g) ds (2.20) where (HF)B Y(3935’) denotes the 8 component of traction at x due to a point force F in the y direction at x’ on a half plane. However, (HF)B Y(§,x’) for an anisotropic half plane is not independent of the coordinate system and, thus, the method becomes extremely cumbersome for anisotrOpic problems. Benjumea and Sikarskie[2o]remedied this by replacing the half plane by the infinite plane. Altiero and Sikarskie [14] and Banerjee [21] obtained independently a formulation for mixed boundary value problems. The idea is to embed the region D in a infinite plane (in two dimensions) or infinite space (in three dimensionsh see Fig. (2-2). A layer of body force RY is then distributed over the surface B in the 13 cue-Plane Fi g . p 14 infinite medium. The stresses and displacements throughout the medium are determined according to the principal of superposition by: umB (g) = Juana“ (5g) Ry(§’) ds (2.21) uB (_€_) =1]; (UR)B-Y (épfi) RY (15’) ds (2.22) where (HR)a8 Y are stresses associated with nun B y. Note that; (HR)B.Y (§.§ ) = (HR)aB.Y (3.5’) na (5). (2.23) As a + x a B, (2.21) and (2.22) lead to: 1 n (x) = — a R (x) + fUIR) (x.x’) R (x') as B — 2 dB 8 - B 8-Y — — Y - (2_24) uB (3E) =‘1[(uR)B-Y (x,x’) RY (x’) ds. (2.25) For the geometric problem the left hand side of (2.25) is known and one solves an integral equation of the lSt kind for the unknowns Ry. For the statical problem the left- hand side of (2.24) is known, and one solves the singular . nd integral equations of the 2 kind for the unknowns Ry. As with the direct approach, not all of the integral equations are Cauchy singular for mixed type boundary conditions. In any case, the solution at field point 5 2D is obtained from (2.21) and (2.22). II.4.3 GENERALIZATION OF THE INDIRECT APPROACH Equations (2.21) and (2.22) generate the elastic field in the region D due to the application of single force RY 15 on the embedded boundary. In a similar way, Lauat [22] utilized a layer of dislocations, Cy, to generate the elastic fields. In [23], the groundwork has been laid for another indirect formulation in which dislocation dipoles, cy, are distributed over the embedded boundary. Other possibilities also exist, such as linear combinations of RY and Cy, or dipoles of RY and CY. Two well known singularities, Massonet's, Fy, and Rieder'sIliH, GY, can be shown to be linear combinations of RY and Cy. If MY is used to denote any of these singularities, then one can write: uB (g) =.£.(uM)B.Y(§,x’) MY (x’) ds, (2.26) HaB (§)=f(IIM)aB. (533’) MY (95’) ds (2.27) where (uM)B Y §,§’) denotes the B component of displace- ment of field point 35D due to application of a singularity MY in the Y direction at source point xfeB,anui(HM&B.y(§,x’) are the stresses which correspond to the displacement (uM)B Y (§,§’). In the case of the application of the dipole N2 of My, one can write (2.26) and (2.27) as: uB (g) =.l[.(um)8.Y (5,5’) mY (33’) ds, (2.28) n (a) = B (111m)“ as _ (5.1211!)Y (95") ds. (2.29) BoY If in (2.26) and (2.27) MY is Ry, Cy, FY, or Gy, then in Y Y Y ar di oles of R C , F and G res ectivel . g'Y e p Y' Y Y Y p y (2.28) and (2.29) mY is ry, cy, f , or gy, where rY, c , f 16 II.4.4 REMARKS Since the dipole of any scalar singularity is a vector, and the dipole of a tensor of rank K is a tensor of rank K+1, one increases the order of singularity by employing (2.28) and (2.29)*. This is due to the fact that: (um) (§,§’) = a B-Y (uM)8 x, (gli’)l (2.30) Y (“New (Ly) ax} (m) .6 8.5 (3.5’) (2.31) which means that the gradient with respect to source point co-ordinates of an influence function of a singularity M6 represents the corresponding influence functionoftfluadipole m» of this singularity. But for plane problems (UM)8 Y is _1 . )1 30 (um) BOY W111 be 0 (lx-x’l-l) and (Him)B Y will be 0 (Ix-x’l-z). Conse- 0 (2n léfé’l) and (HM)B Y is O (l§-§’l quently the boundary integral equation for (2.29) does not exist. This can be verified easily using Hooke's law in the following form: _ 2v HaB — u(Say 686 + 606 58y + 1-2v 6a8 yd (2.32) Now the fundamental solutions (111MB.Y and (um)B.Y must satisfy (2.32), but (111108.Y is O (lg-g’I-l) so that (MUS.Y will be 0 (]x—x’l-z) and its Cauchy Principal-Value (CPV) does not exist. Stresses are economically obtained by numerical * Note that the m used above is technically the vector of the dipole of M1. 17 differentiation of the displacements (see'Table 3.2). Thus the boundary integrals associated with (2.26), (2.27) and (2.28) can be written in the following forms: “a (35) =f(uM)Cl 8 (Eng) MB (g’) ds. (2.33) B . 1 , , um (35) = - 5 GaBmB (35) + {(um)a.8(§.§) InB (g) ds, (2.34) _ l . . Ha (x) — + -2- 6018MB ()5) + immafiqé) MB (33) ds. (2.35) Using various types of singularities, one can form integral equations of the first kind which behave logarithmically, i.e. equation (2.33), or equations of the 2nd kind which are singular and of O (Ix-x’l_l), i.e. (2.34) or (2.35). In any case, the resulting integral equations are solved numerically. However singular integral equations of the Cauchy type may be desirable because they lead to a diagonally dominant system of linear algebraic equations which can be solved by iteration techniques. Such techniques require fewer operations and provide computational savings. Based on equations (2.34) and (2.35), one needs to choose the following singularities to achieve singular integral-equa- tions of the second kind: (i) For geometric problems: (uc) (ug)a B or any linear combinations of these. c.8' (ii) For statical problems: (HR) (HF)a or any linear combination of these. 0.8' B In [24] it has been shown that the various singular- ities do not lead to the same numerical efficiency. CHAPTER III NUMERICAL SOLUTIONS OF THE INTEGRAL EQUATIONS III.1 INTRODUCTION Closed-form solutions of the integral equations of types (2.33-2.35) can be found only for special geometries and boundary conditions, e.g. a circular cavity subjected to hydostatic pressure. Thus it is important to develop approximation techniques for the solution of these equations. The principle of most approximate methods is to discretize the boundary of the region; i.e. replace the real boundary by an n-sided polygon. Integration over the polygon can be carried out either numerically or exactly. The idea of re- placing integral equations by finite sums was first employed by Nystron [25] using the fundamental ideas of Fredholm [l ] and Hilbert [26]. Rectangular formulas were employed and other quadrature formulas were proposed. In this chapter, a general scheme for any arbitrary influence function will be presented. The technique will be useful for developing a general computer program for any choice of influence functions. 18 19 111.2 THE STRUCTURE OF THE GENERAL INFLUENCE FUNCTION The influence function represents the connecting link between the elastic field and the singularities which generate this field. The known boundary conditions are the known parts of the integral equations, and the density of the singularities are the unknowns. Influence functions (which are the integrands of the integral equations) de- pend on the relative coordinates of the field point §=(€l.€2)and the source point x=(xi,xé) i.e. E = (x,y) where X = 51 - xi and Y = $2 - x5, see Figure (3.1). Thus one can write: (NM) (5,5’) = (NM) (3 - x’) = (NM) (X,Y) (3.1) where (NM) (§,x’) denotes the elastic field N at g due to a singularity M applied at source point x’. From Figure (3.1) it can be seen that r = I; - x’l and from equation (3.1) one can write: 35a (NM) (3. 5’) = - axé (NM) (5,5’) (3.2) where 350 and axé denotes differentiation with respect to field point Ea and source point x’a respectively. To develop the numerical scheme, it is necessary to look at the structure of two important fundamental solutions, namely (uR)B Y and (uC)B Y which are given in [19] as: V1 3‘ 6 £nr - V r r } (3.3) u (”mew = 5 BY 3 8 Y {V r -v {V2 e £nr+v 6 arctg (172-) (UC)B-Y 1 BY 6 BY 1 (3.4) + V3 rB eYA rx} 20 ~eb-X >X] Figure 3.1 Geometry and coordinate system for an arbitrary region D. (I In 21 where B.Y.6 and A are indices referring to the coordinate directions, ra=(xa-x;)/r, a=l,2 and V1 through V are con- 6 stants depending on Poisson's ratio and the shear modulus. Other influence functions are either linear combina- tions of (uR)8 Y and (uC)B Y or can be derived simply by differentiation or integration. For example, via Hooke's law, equations (3.3) and (3.4) lead to: - 2.1:; (HR)B.Y ‘ V1 {an r [ V2 58y + 2V3 rBrY] * (3.5) 4- V2 (e6K nK ré/r) eBY} ia‘ _ BE 1 (He’s. ‘ V4 {an r [en r8 + 8A8 rY”) (3.6) + (e6K nK ra/r) day} where ar/an = rlnl + an2 It can be seen that (UR)B.Y and (uC)B.Yare linear com- binations of the terms log r, XY/rz, xz/rz, Y2/r2 and arctg (£4 . Thus any linear combination of (UR)B-Y and (uC)B.Y will have these same terms. Similarly it can be seen that (HR)B.Y and (HC)B-Y are linear combinations of X3/r4, XZY/r4, XY2/r4 and Y3/r4. The functions (UC)B-Y and.(uR)B.Y are also linear combinations of this second set of functions. Differ- entiating (HR)B.Y and (HOB.Y creates a third group of in- fluence functions which are linear combinations of X4/r6, X3Y/r6, XYB/r6 and Y4/r6. The functions (HC)B.Y and (IIr)B.Y belong to this group. As can be seen, members of the third group are strongly singular. They occur in the influence functions of stresses due to dipoles. It can be seen that 22 (uc)B Y is O (%) in two dimensions. Thus the corresponding tractions. (HC)B.y are of 0(1/r2) and their Cauchy Principal‘values (CPV) do not exist. A general method to derive any arbitrary influence function of elastostatics is shown in Appendix A. 111.3 SIMPLIFICATION OF THE INTEGRAL EQUATIONS FOR NUMERICAL TREATMENT The general integral equations (2.33-2.35) which arise in the indirect formulation can be written as follows: NB(§) = a 68y MY (x) + j£(NMg.Y(x,x’)MY (x’)ds (3.7) For N E u the integral is weakly singular and a = O, B B (if MY is replaced by mY the integral is singular and a = -%). For N8 2 Us the integral is singular and a = t, (if MY is replaced by mY the integral is strongly singular and a may not exist). In general one can expand (NM)B Y as follows: HNM)1.1 '311 a12 a13 a14 ais' r Fl ‘(NM)1.2 b a21 a22 a23 a24 a25 F2 (NM)2.1 == a31 a32 a33 a34 a35 ‘ F3 ’ (3'8) (NM)2.2J _a41 a42 a43 a44 a45J F4 . F5 ] 23 whereA= (ak£)' k = 1,...4, (i = 1,...5, represents the coefficients of the linear combinations, and F = F1, (195’) is selected according to the following cases: 1 T g = {log r, XY/rz, XZ/rz, Yz/rz, tg‘ (§)} = u u . . . when (NM)B.Y ( R’B-Y' ( C)B-Y or a linear combination of these; 3 = {x3/r4, XzY/r4, XYz/r4, y3/r4, 1}T when (MUS-Y = (HR)B.Y, (HC) , (ur)8 y, (uc)B.Y or a B.Y linear combination of these; and E = {X4/r6, X3Y/r6, XY3/r6, Y4/r6, 1}T w = . . . f hen (NM)B.Y (HC)B.y' (Hr)8.Y or a linear combination 0 these. The arbitrary constant, unity, is added to the arguments of the two later cases for computational purposes. This addition makes the dimension of A equal for all choices of influence functions. For (uR)B Y: ’ v5 0 -v3 0 o ‘ A = 2% W3 0 0 0 (3.9) 0 —v3 0 o 0 LV5 0 0 -v3 0 - . For (UC)B-Y: ' 0 v3 0 0 v6' .5 = -vl -:2 Z -:3 : Z (3.10) 2 3 . 0 -v3 0 o V64 . ~ 24 For (HR)B : (V5+2)nl (V5+2)n2 V2n1 V2n2 0 v + - .A = V1 2n2 (V3 2)n1 (V5+2)n2 ‘V2nl 0 (3 ll) -V2n2 (V5+2)n1 (V3+2)n2 Vznl 0 _ Vzn1 Vzn2 (V5+2)nl (V5+2)n2 OJ . F‘ : or (RC)B.Y F n1 n2 -n1 n2 0" n 0 0 -n 0 liz: v4 1 2 (3.12) n2 0 0 -n2 0 . n2 -n1 n2 n1 OJ In equations (3.9 - 3.12), Vl = -l/4HV) V2 = V - 1, V3 v +‘i, v4 = -2uvlv3, v5 = 3v - 1 and v6 = 2v where v elv for plane stress, V = (l-vL/v for plane strain and v = A/Z (A+u) is Poisson's ratio. . . I ff' ' c The linear combination coe iCient matrix‘A,for other types of influence functions can be obtained using the method outlined in appendix A. 111.4 APPROXIMATE INTEGRATION OVER A DISCRETIZED BOUNDARY The boundary B is divided into n intervals and the general integral equation (3.7) is discretized as foll ows: i _ NB (BE ) -- aOBY M8 (Xi ) + i f ' (NM)81'X’ ‘ ; j=l B. Y — )R%((§ ) dSI i=l,,.n (3.13) 3 25 where Bj are sides of the polygon obtained by discretiza- tion of B: see Figure 3.2. Equation (3.13) leads to a 2n system of .linear equations with its off diagonal elements obtained by integrating the influence function over each interval and.its diagonal elements determined by the argument givem.in the discussion between equations (3.7) and (3.8). Since in the singular integral equations constants can be taken outside of the integral sign,[ 5 ], one can use (3.8) and reduce (3.13) to a linear combination of the following integrals: “IF (x,x’) {Ml (x’) + M2 (g’)} ds, a = 1....5. (3.14) a _ If the resultant boundary data (i.e. NB (xi) are defined at the mid—point of each interval, the weight of the singula— rities (i.e. MY (5’)) can be considered constant over each interval. Consequently (3.13) can be written as a linear combination of simpler integrals, b 4f Fa (5. 5’) as = Wu ' (3.15) To evaluate the integrals Wu numerically, a general quad— rature formula is considered. One can write; m We = K221 AK Pa (:5, ix) + e, (3.16) where AK' EK are weights and integration nodes respectively and e is the quadrature error. Depending on the type of quadrature used, one can employ the following AK and Z . _K° 26 A Integration nodes lk Figure 3.2 Region D of Fig. 3.1 after discretization, 27 (i) Rectangular Formula: Z = ___ _ b-a -1 Er £2 a + —-. ---Z = §_+ (m-1)=fi=, and (ii) Tangential Formula: _ 2-9. b—a Ei‘e+‘2—firéz=e+3=zn=vw _ 12-2 lb-al ER) — 3 + (2m-l) —fi_' and AK = ‘m- . (iii) Chebyshev's Formula: 2'2 (m) -Z-K T+T§K “‘3sz 3| (m) - where ER are Chebyshev's pOints, the roots of Chebyshev's polynomial of degree m, (note that m i 10 [29]) (iv) Gauss' Formula: I -e| ___, (m) _——— _ =§_ + —j— EK and AK - 2 A )6 (m) K (In) where EK are the roots of the Legendre polynomial of (m) . 9 0 degree m, and AK are the GaUSSian coeffiCients for the interval (-1,l)- III.5 EXACT INTEGRATION OVER THE DISCRETIZED BOUNDARY The integrals of equations (3.15) can be analyt' 1cally evaluated over the discretized boundary ° Except for two terms, log r and arctg (Y/X)) other terms have s' .1 lml ar structures. Thus it is possible to generate a s cheme to A 28 evaluate these integrals. It can be seen that all of these integrals are contained in the following form: n 100 = [gig—ds, K=l,3,m,n=l,...4 and m,n B r2K (3.17) m + n _<_ 4. Relating X, Y and r to the arc length 5 will facilitate the method. Thus the following relationships are considered: X=a +aS o 1 Y = b0 + blS (3.18) 2 __ 2 r _C0+CIS+CZS where a0, al,... , C1 and C2 are constants which are not lyre- sently known. Substituting (3.18) into (3.17) leads to the following equations: m n k) _ f(a0 + alS) (b0 + b1S) ( d5 ° (3.19) Im'n B (C0 + cls + c252)K Considering an arbitrary interval Bc, see Figure (3.3) and using the triangle ABC obtained by connecting a field point, A, to the extreme points B and C of the interval, one can define the angles a, (p and o as follows: 2 2 h +h _ 2 X_XB h a=arccos (ch )) V=arccos(—l—2h\2___) h 2 2 2 12 hl +h _hz and (b = arc cos (—————2—h-IH—_____) where O _<_ ‘1), (1’ f. " and XB’ Xc' h1' h2 and h are defined in Figure (3.3) . Consequently, one can write; 29 O >x' Figure 3.3 Geometry and coordinate system for exact j4ntegration over a straight segment. 30 X = x1 - x 1 = x1 - xB - Scosa Y = x2 - x’2 = x2 - YB - SSina (3.20) r2 = h 2 - 2h S cos¢ + 82 1 1 Comparing (3.20) to (3.18) one can find: a0 = x1 — xB a1 = - cosa bo = X2 ' YB bl = - sina (3.21) _ 2 Co ‘ hi C1 = - 2hl cos ¢ C2 = l . Substituting the constants obtained in (3.20) into (3.19) one obtains: m . n I(k) = Jf(xl - xB - cosa S) (x2 - yB - Sina S) ds. m'n (s2 - 2h cos ¢ 5 + h 2)K (3.22) IBxpanding the numerator in l 1 (3.22), one can write the in- izegrals in (3.22) as a linear combination of the following auxiliary integrals: Sm ds. (3.23) m,k (S2 - ZhIScos uplle integrals in (3.23) are ()Ile can also use a suitable all integrals in (3.23) are jLI'ltegral: 2 ¢ + h1 )K easily evaluated using [28]. change of variables and show that obtainable from the following I— 31 Q =./” ds 0'1 (s2 - 2hlScos ¢ + hlz) [us-h d5 2 2 . 2 1 cos 9) + hl Sin 9} (3.24) 1 S - hl cos o h = ————*——— arctg ( . ) hl Sln ¢ hl Sln ¢ 0 hl Sln ¢ Consequently: _ h Ql,l — h1 cos o Q0,1 + in ( 2/hl) (3.25) h - h1 cos ¢ cos 9 2 2 (3.26) Q0,2 = i h 2 + TQ1,1} / 2h]. 511'?! d) 2 _ l l __ 1 01,2 - h]. cos it QO'Z + f (_2 2) (3.27) h h l 2 h2 2 (3.28) Q2,l = h + 2 hl cos ¢ Kn ( /hl) + hl cos 2 ¢ Q0,l Q = Q - h 2 Q + 2 h cos ¢ Q (3 29) 2,2 0,1 1 0,2 l 1.2 ' _ _ 2 Q3,2 - Q1,2 h1 Q1,2 + 2h1 cos ¢ Q2,2 . (3.30) Cine should note that if sin ¢ = 0, then: 1 l l l 1 Q = |—— - ——| and o = — |——— - ———|. (3.31) 0,1 h2 hl 0,2 3 h23 hl3 IQ(D‘w the w integrals are easily evaluated by decomposing (13.22). 32 III.5.1 INFLUENCE FUNCTIONS OF THE FIRST GROUP When (NM)6 Y = (uR)B Y or (uC)B Y, one obtains: ) 1 ‘= (X1 ’ XB) (X2 ' YB) Q0,1 (3.32) - {(xl - xB) sina + (x2 - yB) coso} 01,1 + sina COSa Q2 1. I _ (1) _ _ 2 _ _ W3 “ I2,0 " (X1 XB) Q0,1 2 (X1 XB) COS“ Q1,1 (3.33) + COS d QZ,1° _ (1) _ _ . W4 _ 10,2 _ h W3 (3.34) IJsing integration by parts and a suitable change of variables one obtains h cos ¢ Kn (hl/h2) + h (£n h W = - 1) + h w sin ¢ 1 1 1 1 (3.35) W5 = h GB + {(x1 - xB) sina - (x2 - yB) coso} Ql,1- (3.36) Chne should note the following special cases of equation (3.35) (:1) When ¢ = 0 or TT: wl = (hl zn h2 (n hz) + h2 - hl- (3.37) (:11) If the field point locates on the interval: W1 = h (in h - 1). (3.38) 33 III.5.2 INFLUENCE FUNCTIONS OF THE SECOND GROUP Similarly for (NM)B Y = (HR)B.Y' (HC)B-Y' (UC)B.Y' (ur)B Y one can obtain: _(2)_ _3 _ -2 Wl — 13,0 — (x1 x3) Q0,2 3 (x1 XB) cosa Q1,2 2 3 (3.39) +~3 (xl - xB) cos a 02,2 - cos a Q3,2 2 _ (2) _ - _ W ‘ I ' (x1 XB) (X2 YB) Q0,2 2 2,1 - (XI—x3) {(xl-XB) sina + 2(x2-YB) 0059} 01,2 2 . + {(x2 - yB) cos a + (x - xB) SinZd} 02,2 1 . 2 - Sino cos Q3,2 (3.43) W = 1(2) = (x - x ) Q - cosa Q - W a: 3 1,2 l B 0,1 1,]. l (3.21) W = 1(2) = (x - y) 0 - sino Q - w (3.42) 4 0,3 2 B '0,1 1,1 2 13131.5.3 INFLUENCE FUNCTIONS OF THE THIRD GROUP For (Hm) where m are r , f , c ,...etc. one obtains 8.Y Y Y Y Y 3 2 2 _ 3 w _ I + 4a0 al Q1,3 + 6a0 al Q2’3 ( ) _ 4 1 4,0 ‘ ao Q0,3 3 4 3.43) + 4a0al Q3,3 + a1 Q4,3 ( _ (3) _ 3 2 w2 ' I3,1 ‘ ao bo Q0,3 + a0 (3alb0 + aObl) Q1,3 2 3 3 1,3 = aobo Q0,2 o) Q1,2 + albl Q2,2 - W2 (3.45) W4 = 16?; = Q0,1 + W1 ao2 Q0,2 lQl,2+a1202,2) (3.46) where 00,3 = (c1 + 28)(l/2r4+ 3/(4c0 - c12)2r2}/(4c0 - C12) + 6 Qo’l/(4c0 - cl ) (4.37) Qm,3 = ‘{Sm-l/r4 + (3 ' m) c1 Qm-1,3 ' ‘m‘l)Con-2,3}/ One can see first group (90,1' (s-m), m=l,2,3,4 that all integrals except W (4.38) and W of the 5 are derivable from the fundamental integral Subroutine "INT" (see Appendix C) evaluates the boun- Ciary integrals both numerically and analytically, Some eexamples were tested, and better results were obtained Ilsing the analytic method with less time and fewer opera- tions . .- .1. fffi-d-v—u u'l 35 III.6 NUMERICAL RESULTS In this section, some numerical results are presented for a few test problems. Regions with straight boundaries are examined. Example 3. l A rectangle generated by the points (.5, .75), (-.5, .75), (-.5, -.75) and (.5, -.75) is examined. The follow- ing boundary conditions are considered: Case (1) = 1 .25 E u on x1 = i .5 _ (l-v) - .5 u x2 = -.5 3 x1 u on x = + .75 2 _ = + .375 ‘1'“) - u 0 everywhere and H = i 1 CH1 x2 = i .75 U ==() on x1 = i 5 H2 = O on x1 = i .5 -52Xl u on x = + .75 2 _ + .375 (l'V) - u uplle following formulations are employed: 36 For boundary conditions of type (1): Method (i) - Equation (2.33) with (uM)a B a (uR)a B is employed for evaluation of R Internal displacements are 8. obtained using the same equation. Stresses are obtained from (2.27) with (HM)aB Y = (HR)aB Y. See table (3.1). Method (ii) - Equation (2.34) with (um)a B E (uc)a B is used for evaluation of CB' Internal displacements are obtained using (2.28). Stresses are obtained using numeri- cal differentiation of displacements which avoids the highly . See expensive process of (2.29) with (Hm)a E (HC)aB Y B-Y table (3.2). Iflethod (iii) - Equation (2.34) with (um) E (ug)a is 8 Iised which corresponds to applying a linear combination of the a.B (iislocation dipole c and the force dipole r Internal dis- 8 8' Inlacements are obtained using (2.28). Stresses are obtained \ria numerical differentiation of the displacements. See table (3.3). IPor boundary conditions of type (2): Equation (2.35) with (HM)a s (HF)a is used for the B .8 Eivaluation of F8' Then stresses are obtained via (2.27) With (HM)(YB Y s (HF)a . See table (3.4). B-Y I-"'~:>r boundary conditions of type (3) Equations (2.33) with (uM)a (uR)a on Bu and .B .B (12.35) with (HM)a B s (HR)a B on Bt are used for evaluation 37 of RB' Displacements and stresses at field points are obtained via (2.26) and (2.27) respectively. see table (3.5). Example 3.2 A square plate is generated by the points (1., 1.), (-l., -1.), (1.,-l.) and (-l., 1.). The variable boundary conditions corresponing to the stress function: 4 - x 4x - x 2x 3 ' e ‘ 1 2 1 2 1" H — 6x1 x2nl + (6xlx2 4xl )n2 _ 2 _ 3 2 _ 3 H — (6xlx2 4xl )nl + (12xl x2 2x2 )n2 are applied to the boundary. Equation (2.35) with (HM)a 8: (HP) is used for the evaluation of F Then stresses are 0.8 8' obtained Via (2.27) Wlth (HM)aB Y= (HF)aB.Y see table (3.6). The boundary of each region is divided into 40 equal segments. All integrals are evaluated analytically. Com- putation is done on a CDC 6500 computer and the times in parenthesis indicate the required CP execution time. 38 5.00m 00.0 m0 000.H 00.0 00.0 00.0 mNN.0 mNN.0 w0N.0 mhm0.0l ohm0.0I 0.0 m.0 v 000.H 00.0 00.0 00.0 mNN.0 mNN.0 00.0 000.0 000.0 0.0 0.0 m 000.H 00.0 00.0 00.0 000.0 000.0 00.0 mhm0.0l mhm0.0l 0.0 m.0 N 000.H 00.0 00.0 00.0 000.0 000.0 00.0 000.0 000.0 0.0 0.0 H HOHHMm uooxm .Ebz monumw uomxm .Edz fix .Pa .02 mm NH Ha usaom = = = m: H5 mwumcwpuoou macaw >. r. maAm=0 . eamsv upmflammm wpfluwasmcfim cwmuum mamam mm.u> ..Hn: suw3 unflpmoa AHV mmmo H.m stances mememxm new mueemom H.m memee 39 A.omm N0.00 m0 5mm.o oo.o oo.o om.o m-.o vmm.c oa.~ mnmo.o: hwmo.0I 0.0 m.o v mmm.o 00.0 00.0 om.o mmm.o vm~.o 00.0 000.0 000.0 0.0 0.0 m mmm.o 00.01 00.0 00.0 000.0 000.0 om.a mnmo.o| memo.0| o.o m.o m mmm.o 00.0 00.0 00.0 ooo.o 000.0 00.0 ooo.o 000.0 0.0 0.0 A HOME? uomxm .Ecz Houumw uooxm .Ecz Nx HR .02 mm: NH: Ha: H ucwo ms 5 mmumcflpuooo came n.5Aosv undefined muflumaswcflm cflwuuw ocean mm.u> eat .Hue sues essence Ids mmeo H.m aoflnonm weasexm hoe mueemom ~.m mqmae 4O it.‘¢§ 0.. .r H. 1.66m em.ov mo mmm.o 00.0 00.0 om.o mmm.o vmm.o om.o memo.ot memo.0| 0.0 m.o v mmm.o 00.0 00.0 om.o mmm.o v-.o 00.0 000.0 000.0 0.0 0.0 m mmm.o 00.0 00.0 00.0 000.0 000.0 oo.o memo.o| memo.ot 0.0 m.o m mmm.o 00.0 00.0 oo.o 000.0 000.0 00.0 000.0 000.0 0.0 0.0 H HOHHMw uooxm .EDZ HOHNMw pooxm 352 mx HR .02 mm: NH: Ha: ms H: mmumcflpuoou WMWWM n.5AHsvmrm WW :m+m.eaosvmm.aimsv ”pmwammfl muflnmasmCHm cemuum mcoam mm.u> 0cm .Hu: :ufl3 ocflpooa Adv ommu H.m eoenoue madamxm hoe muesmmm m.m memee 41 TABLE 3.4 Results for Example Problem 3.1 Case 2 loading with u=l. Singularity used: (HF)B YE(HR)B Y with v=.5 Field Coordinates P§:Pt x X H11 H12 H22 1 2 1 0.0 0.0 0.00 0.00 0.996 2 0.0 0.2 0.00 0.00 0.997 3 0.0 0.4 0.00 0.00. 1.000 4 0.0 0.6 0.00 0.00 1.000 5 0.2 0.2 0.00 0.00 0.994 6 0.2 0.4 0.00 0.00’ 0.998 7 0.2 0.6 0.00 0.00 1.000 8 0.3 0.2 0.00 0.00 0.990 9 0.3 0.4 0.00 0.00 0.994 10 0.3 0.6 0.00 0.00 1.000 Exact 0.00 0.00 1.000 CP (.69 sec.) \ “'19". 42 A.oom mm.ov mo mmm.o 00.0 00.0 om.o mmm.o vmm.o 00.0 mnmo.ot memo.ol 0.0 m.o v wam.o 00.0 00.0 om.o mmm.o vmm.o 00.0 000.0 000.0 0.0 0.0 m 500.0 00.0 00.0 00.0 000.0 000.0 00.0 memo.o: memo.ot 0.0 m.o m ham.o 00.0 00.0 00.0 000.0 000.0 00.0 000.0 000.0 0.0 0.0 H HOHHMw 0.0me .852 Houummw uomxm .55 «x ax .02 mm NH HH ucflom a = = m: as mmumcwcuoou pamflm m.5 m.5 Am=v cam Amcv "wwwammd mumeHsmcfim aflouum mcoam mm.u> 0cm .Hun sues mcflpooH Amy wmou H.m smenonm meeeexm hoe muesmmm m.m mgm.H mm.m w>.H uonum .xmz 0mH.m 00H.m om0.H 000.0 me0.m| HOH.m| 0.0 0.0 0H mHm.v 00~.e omw.0 mvm.0 Hmm.ml vmm.m| me.0 me.0 0 «m0.HI VHO.H| 000.0 000.0 000.0 000.0 0.0 0.0 0 000.H mmm.H www.ci mmm.0i 000.0) He0.0| v.0 0.0 e 0mv.m mmv.m 000.0 050.0 0m0.~| 0e0.m| 5.0 v.0 0 00H.~ emH.N mmv.0 Hmv.0 00m.HI NHm.HI 0.0 0.0 m 0m~.H 0m~.H 0m~.0 mv~.0 0me.0i m0e.0i m.0 m.0 0 000.0 000.0 «00.0) v00.0) 000.0 000.0 0.0 0.0 .m 000.0 000.0 0mm.0n mmm.0| 000.0 000.0 0.0 «.0 m 000.0 000.0 000.0 000.0 000.0 000.0 0.0 0.0 H pomxm .Ecz uomxm .Edz uomxm .Esz mx Hx umwmm mm: NH: HH: moumchuoou pHWHm W H": .m.n> nuH3 0.5Amzv m 0.5Amcv "poms muHHMHs0ch mmx Hx I «x Hx u ANx.Hx0e :oHuocsm mmouum may on mCHpcommwuuoo mcoHuomue "mCHpmoq m.m EmHnonm mHmmem mom MUHsmmm 0. m mHmdfi CHAPTER IV BOUNDARIES OF PIECEWISE-CONSTANT CURVATURE Replacing the real boundary by a polygon simplifies the calculation of the integrals and enables one to gene- rate the scheme developed in Chapter III. Obviously the method is suited only for boundaries composed of a finite number of straight segments along which the unit outward normal is considered to be constant. For irregular boun- daries (curved) the method becomes crude and increasing the number of subdivisions can not guarantee improvement in re- sults because of round-off error and the numerical insta- bility of a large system of equations. The most suitable numerical technique for curved boun- daries is the Romberg-Richardson extrapolation, which is based on error analysis [29]. The method converges to the exact solution of the non-singular integrals but it is re- latively costly. This method was used by Heise [19] to evaluate the boundary integrals for elasto-static problems. In this chapter the general boundary integrals of elastostatics are evaluated analytically over irregular boundaries using the assumptionmof piecewise constant curvature. 44 45 IV.1 FORMULATION OF THE METHOD Consider an irregular region D (which is in general multiply connected), the boundary of which is denoted by B (Figure 4.1). Any internal field point is denoted by In (51,52) and any boundary field point by x = x (s) = (Xlfliz) where s is a measure along the boundary. When it is necessary to distinguish between a boundary field point and a source point, the latter is denoted by x’ = g’ (s) = (x’1.x’ The unit outward normal and unit tangent at 2). a boundary field point are denoted by nu = na(6) and ta = ta (0) respectively (Figure 4.2). To calculate the elastic field N(§) (which could be the displacements 1%3(§) or the stresses HaB(§)).OneCKNP- structs the following integral equation: N(€) =.]. (NM) (5.x’(s’)) M (8’) ds’ 36D (4.1) B 07 — Y As g + x(s) 63, the above equation becomes one of the fol- lowing*: uB(s) =J£ (uM)B Y (s,s’) MY (3’) ds’ ; or (4.2) u8(s) = [um]B.YmY(s) + fl; (um)8.Y (s,s’) mY (SOS-:33) H8(s) = [HM18-Y MY(s) + (E(HM)B-Y (s,s’) MY (3’) 32:4) where the influence functions depend on the relative coor- dinates §_= xfis) - x (s) and the orientation of the normal * Recall equations (2.33) through (2.35). 46 CDC—Plane Figure 4.1 An irregular region with curved boundary embedded in the infinite plane. 47 +X] C) Figure 4.2 Boundary field point x and source point xi 48 at angle 9 (S)- If the parametric equation of the boundary B is not known, one must replace the boundary by a finite number of arcs. Then the general influence function (NM)(s,s’) is a function of the center of curvature, radius of curva- ture and angle 6, i.e.: (NM) (315’) = (NM) (EEO, Dr 63126 I 036’) (4-5) where x0 = (xo,yo) and p are the center of curvature and the radius of curvature of the arc respectively. These are related to the field point 5(8) and the source point §’(s’). The non-integral factors which arise in (4.3) and (4.4) de- pend on the normal to the boundary [10], [16]. A suitable number of circular arcs has many advantages over a variable radius of curvature because, in the first case, the radius of the curvature is along the normal to the boundary and one does not introduce any new coordinates (Figure 4.3). Since all influence functions can be derived from essential components using a suitable linear combinations, it is more advantageous to consider these essential com- ponents and justify them according to the notation used in Figure (4.3). To simplify the calculations, one can move the origin of coordinates to the field point and obtain: X = x0 + p cos 6 (4.6) Y yO + 0 sin 8 (4.7) 49 x Y (3 H M2 9, ‘0 (l'!7 Ml ' 5' l I 6// / g. / 5. 0 i . /‘// / S // ‘V /’ "I ,5 >x g D N2 8 0 ”I Figure 4.3 Geometry and coordinate system related to the formulation. 50 2 _ 2 2 2 . r — x0 + yo + p + 2pxo cose + 2p yO Sine. (4.8) Letting f(0.0) = (X02 + Y 2 2 .. + in0 0 + p + prO cose 2pyo s ) one can obtain: 3 = {glogfflo.6). (x0 + pcoseHyO + osin6)/f(0.6). (x0 + pcose)2/f(p.e). (yo + psine)2/f(p.e). (4.9) “l [(yQ + psin8)]}T (XO + pcose) t9 when (NM) = (uR) , (uC) or a linear combination 8 Y B-Y B.Y of these; _ 3 2 2 E — {(x0 + pcosa) /f (0.6). (x0 + 0cose) (yo + pSin8)/f2(p,8), (X0 + pCOSG) (4.10) (yO + 05in9)2/f2(0.9). (yO P psine)3/f2(o.e).1}T when (NM)B-Y = (HR)B.Y' (HC)8.Y' (uc)8{Y,(ur)B.Yor a linear combination of these; and 4 3 3 E = {(XO + pcose) A? (0.9). (XO + pccse) (yo + psin6)/f3(p.6). (x0 + 0cose) (4.11) (yo + psine)3/f3(p.e), (y0 + psine)4/f3(p.e).1)T when (NM)B Y = (Hc)B Y, (Hr)B Y or a linear combination of these. In order to simplify the calculations and to be able to present the results in a suitable and clear form, one can consider the integrals of log h3.8) and 51 (x + psine)]'separate1y and write the following 0 t9 form for the other components: G(K) = (Xo + pcose)m (Yo + psine)n ds m,n B (x02 + yo2 + p2 + chose + 2psin8) where, m,n = l,2,3,4; m + n i 4; K = 1,2,3 and ds=K pde. (4.13) c The angle 0 and the convexity factor KC must be defined according to Figure (4.4). To evaluate the integrals given by (4.12), the inte- grand is decomposed. Neglecting Kc momentarily one obtains: G(1) 1) +X 2 (1) _ ( 2 (1) 3 (1) 1,1 ' pxoyo H0,0 0p H1,0 + yop H0,1 + p H1,1 62:, = as); + 421+ 03 4%; as}; = as}; + H132. + .3 4%) Gif0 = px03 Héf0 + 3°2X02 3631+ 3p3x0 H0?i + 4 H0?i Gi?i = °X§Y03020 + 02"02 Hi?0 + 202x0Y0 H0?i + 203x0 Hifi + 3y0 “6 i + ”4 Hifi 43; = 422, + 432, + 2 .2 as + 1421+ 43:, + 4 4?; G62; = py03 H030 + 3p2y02 Hi?0 + 303),0 H230 + p4 H(2) (4.12) (4.14) 52 ucmemmm >Hmpccon xm>cou m.v musmHm . o H+ u x ucwEmmm mumpcson w>mocou «.0 oustm ,«xm 1< / .\ m 0 0 0 .a \ 4‘ / \\ l. e\./\0. 53 where H(K) _Jf +sinm 6 cos n6 (a0 cose + co sine)K d9 (4'15) and a = x + y 2 + op2 b = 20x c = 20y Using 0 02 0 ' 0 0' 0 0' trigonometry and partial fractions, one can show that: H(1) _ ;L_ _ H(l) _ (1) H2 0 — c0 {cose aOH1 O , b0 Hl,l} (4.16) H(1) __ (l) _ (1) HO 2 — H0,0 H2,0 (4.17) (2) _ c0 cose — bQ sine (l) H0,0 - {(a0 + bO cose + cO sine) a0 H0,0} (4.18) 2 2 2 / (aO - bO - cO ) (2) _, (2) _ (2) H0,2 — HO,0 H2,O (4.19) (2) _ H(2) _ (2) 2’1 - HO 1 H0’3 (4.20) H(2) __ H(2) _ (2) H1, 2 — Hl’o H3’0 . (4.21) To be able to evaluate the other necessary integrals in (4.15) one must define an angle a such that: ct): 8 +0 (4.22) and a0 + b0 cose + c0 Sine = a0 + dO Sin 0 (4.23) where 2 _ 2 2 _ -1 b d0 — b0 + c0 and a - tg ( O/co). (4.24) Now the required integrals of (4.15) become: (l) 54 H0,0 = P0,0 Hi}; = c050 Pifg - sina Pg}: Hg}; = Sing P138 + COSa PST: Hi}; = - % sin 20 Péfé + cos 20 Pifi + sin a P H{?g = cosa P{?g - sina Pé?i Héfi = sina P{?é + cosa Pé?i H{?1 = - % sin 20 P5?5 + cos 20 P{?i + sin 20 P; Hé?é = sinza Pé?g - sin 20 Pifi + cos 20 P H§?é = 35in20 coso P335 - 3sin0 c0520 Pé?i + cos 3o Pé?g + sin 30 Pg?; Héfg = 35in0 coszo P333 + 3Sin20 coso Péfi - sin 30 P§?$ + cos 30 Pg?; where P(K) = sinmo cosne do m, (a0 + dO Sln¢)K The following integrals are easily evaluated: pél’l = 2n 'f(¢)/d0 06?; = -1/dO f(6) The remaining integrals have the following form: 2) ,0 (4.25) (4.26) (4.27) (4.28) 55 pifg = (4 - aOP P‘l’vd0 Pi}; = (sin¢ - aOP P(1))/dO Péfé =-(cos¢ + aOP(l))/d0 Péfg =-(d0 cos¢/f(¢) + aoP P(l))/El P{?% = (a0 cos0/f(0) + dO P(1))/El P{?1 = (Péfi - sin¢/f(¢))/d0 éié= 4 - P5?) - P§€é =-(cos0 + 2a0d0 P230 + a02 :(2))/d02 Pé?; =-(co52¢/f(¢) + 2 p{}1)/do Where El =_(a02 ' doz) f(¢)=a0 + dosin0. (4.30) It can be seen that P(1) (2) nd P(2) (1) from P0,1 (l) 2,0 1,1’ P1,1 a 0,3 are evaluated (see 4.27). One can also show that =_(cos¢ + aoo/do - a 2P(1)/d0 )/d0 0 56 (2) _ 2 2 3 2 (1) 2 (a0 " 2‘3‘0‘30 ) P0,0 /Eld0 Pm - -{E ¢ + 25(3 «p - a 2d ooscb/f(¢>)) + a 3cos(¢)/f(¢>)}/E d 2 3,0‘ 1 1 o o o 1 o . 3 2 . (l) 2 -{2a(a0 -2aod0 ) + aoa} P0,0 /Eld0 (4.31) Equations (4.31) shows that all necessary integrals are (l) obtainable from P0 0. (l) IV.l.l EVALUATION OF P0 0 I Consider the triangle obtained by connecting a field point g to the center of curvature g and the source point 0 g’. From Figures (4.6) and (4.7) one can find: 2 2 _ 2 2 2 2 2 2 a - d — (D + X0 + YO ) ' 49 (X0 + YO ) (4.32) _ 2 2 2 which indicates that the following cases exist for the boundary integral equation method: Case (1) a0 = dO In this case r, p and the unit normal n’ are collinear, (1) Figure (4.7). P0,0 is easily integrated as follows: (1)__l_ dd) ____1_. P0,0 _ a0 1 + sin¢ — a0 tan ( ) . (4.33) hfle '1 4 57 Case (11) a0 > d0 In this case r, p and the unit normal n' are not collinear, Figure (4.6). P(l) is again integrable [28] 0,0 and the result is: P0,0 _ arCSln {a0 + dO sin¢}/(ao ‘ do )- (4-34) Once P613 is evaluated, the P I (k) ' (4.31). Substitution of the Pm n into (4.25) leads to the P(k) (k). Finally the am n can be obtained from I I (k) are obtained using values of the H (4.14) o The fundamental integrals defined by (3.15) can now be evaluated for arcs of piecewise - constant curvature: IV.1.2 INFLUENCE FUNCTIONS OF THE FIRST GROUP _ (l) _ 3 . _ 2 _ . _ W2 — Gl,l —_Kco COSZGSln¢/d0 ch (yocosa XOSlna . _ 2 . (l) apc052a) 1n f(¢)/dO * {(xoyO p Sinacosa) P0,0 (l) I 2 ' (l) + p(x0cosa+y051na) P1,0 + p 51n2a P2,0 } ch _ (l) _ 3 . . _ 2 _ W3 — 62,0 — ch SinZaSin¢/d0 ch (2x0cosa (l) ‘. 222.2 apsanQ) 2n f(¢)/d0 KCD{(XO +0 +0 Sin a) P0,0 . (l) __ 2 (l) + 2x0p31na Pl,0 o cosZa P2,0} 58 _ (l) __ W4 — 60,2 — che - W3, (4.35) where a’= aO/do. (4.36) W1 and W5 are non-elementary integrals. This is one of the disadvantages of the direct formulation. W1 is evaluated using suitable recursion formulas: K Wl =-7§ 0.] Kn (a0 + do sin¢)d¢ K =-7§ p (w(i) + W(i)) where (4.37) (l) _ W l - ¢£n a0 and (4.38) n n- sinnfi cos¢ In = Jf31n ¢ d¢ = —E— In-l - n . (4.40) To evaluate W5, first integrate by parts: _ -1 W5 —-ch .ltan Yq_+ psin¢ {X0 + pCOS¢} d6 _ -1 yfl + psin¢ _ 2 (l) _ (1) _ 0X T(l)} (4.41) where: 59 T(1) = f ede 0,0 (a0 + b0 cose + co Sine) 'r =f a sine de 1,0 (a0 + b cose + c sine) 0 0 (1) _Jr 6 cose d6 T0,1 — (a0 + b0 0056 + c0 sinG)’ (4'42) The integrals Tél% and Tfla are non-elementary [28]. One I I can show that: (l) _ _ (l) (l) To,o " “ P0,0 + s 0,0 (l) _ _ (l) _ . (1) (l) Tl,0 - a H1,0 Sina 81,0 + cosa 80,1 (1) _ _ (l) (l) . (1) T0,l — a HO,1 + cosa 81,0 + Sina 80,1 (4.43) where S(1) =j' ¢>d¢ . 0,0 a0 + d0 Sin¢ 5(1) =f ¢sin ¢ d4: 1,0 a0 + d0 sin¢ (4.44) (1) _f ¢cos ¢ deb - - ._ 50,1 — a0 + do sin¢ — [¢ in (a0 + d0 Sin¢, 2W1]/d0. It can also be shown that: 2 8(1) = [L _ 8(1) . 1’0 (.2 a0 0’0)/d0. (4 45) Hence, one must evaluate 831g for the following cases: I Case (1), a0 = d0 In this case 5010 is easily evaluated as follows: I s<1>=_1_f___¢.d¢._____ _ 1-1 0,0 a0 1 + sin¢ { ¢ tan 4 2) n ¢ (4.46) + 2 Zn cos (4 - 5)}/a0 *‘I 0 Figure 4.6 The case 4i¢ 9’. 61 x Y ()2 J) A B n, / .— / 3' 6 / ‘3’ / 9 / ’ A // ./ 6% z /<%\ ¢ 1o 3’ ‘0 \3) +x l +4 0 I Figure 4.7 The case (‘1 : 9', 61 x Y 112 J) J) B n, / '- / u' 0 / ‘3’ / 9 / /’ A / / '/ 9s ’ ’(9‘ is 1. > (o ‘3! 1 IX 0 >111 Figure 4.7 The case #1 _-_ 9' 62 Case (ii), ao > do (1) In this case S0 0 is not evaluated in terms of I elementary functions [28]. When a suitable recursion formula is used, one can write: (1) _42 °° _ n (1) So,o"‘2" +331 (1) In (4.47) where I‘i) = ./4 sinn ¢ d¢ = (sinn ¢ - n4 cos¢ sinn-lw/n2 5.1;, n > 2 (4.48) 1(0) = ¢2/2, 1(i) = sin¢ - ¢cos¢ and I(%) = (sin2¢ - ¢sin2¢ + ¢2)/4 . The case 6 = 6’, in Fig.(4.8) must be considered separately. One can obtain: 6 __ __ B . 2 W1 —J{log r ds — KC 0 j" Kn (2 0 Sin 2) de 6A _ (3) —-2Kc p {60 Zn 2p + w 1 } (4.49) 90 Where W(i)= Kn (sin %) d8. 0 Since 2: 9—9-8333 = - Zn (2 sin %) one can write: n=1 w‘i’=- Z §1—:42—‘E 016111 (4.50) n l 63 _|'_|_ / m‘::e: \ IPQ/ / // Figure 4.8 The singularity case (9:8,)- 64 eo _ sine cose f.l~—7 d3 --2 Kc p sinZIe/Z d6 0 (4.51) =-KC p (Km (1 + cos 60) - cose0 + l - Zn 2} 6o 2 _ x2 __1 cos a WB—IFst-ZKCD{ sin2672 d8 =-KC p {sineo + tan e0/2 - 60} (4°52) w = nyz ds =-l K 60 Sinze d6 4 :7 2 c °.f sin2 672 0 =_ ' (4.53) KC 0 (Slneo + 60) 8 _ -1 Y _ o -1 W5 —./ tan (3) ds ——2KC pdf tan (tan 6) d6 0 --K 62 4 4) '- Cpo' (.5 IV.l.3 INFLUENCE FUNCTIONS OF THE SECOND GROUP In this caseIClosed-form solutions exist for all integrals: w = G(2) _ 3 2 2 . 2 . _ 1 3,0 ‘ ch{(xo +3X 0 3X 0 Sin a)(b Slne 0 0 0 c cose)/Elf(¢) + (3x02pcosa+p3sina)/d0f(4) 0 +3xopzsin2a(£n f(¢)-dosin¢/f(¢))/d02 - + _ 3 2_ 2 . 2 (1) ch{ a0(x0 +3xoo 3x00 Sln a) P0,O 2 . 3 . 2 (2 ) (3X0 oSlna+3o Sinacos a) P1 0 - W 65 3 . (2) 2 (2) - o Sin3a P + 3x0p cosZa P2,0 3,0 03cos3a Pézg} (4.55) 2 2 . 2 . 2 2 -ch{ (xO yO—xop SinZa—yoo Sln a+yop )(bosine- c cose)/Elf(¢) - (2x0y pcosa-xozo-o3+ 0 0 3osin3acosza)/dof(¢) + (2x 02c032a+ 0 yoozsinZa)(zn f(¢)/dO-sin¢/f(¢))/dO - p3sin3a(-cosz¢/f(¢)-2(sin¢-a£n f(¢))/dO)/d0} + 2 2 . 2 2 . 2 -Kco{-a0(xO yO-xoo Sin2a+yoo -yoo Sln (1)/El P0,0 + ((x0 + (2x0 I o O O 0') w Q "U 2o+o3)cosa+2x 0Y0 2 . _ 2 ( o Sin2a yoo cosZa) P2,0 (l) psina-3p3sin2acosa) Pizg 2) -KC02(¢sina+cosa£n f(¢>)/d0 - _ o . ( ch(xo a051na) P0, 1) 0 -W chz(sina£n f(¢)-¢Cosa)/d0 - - ( Kco(y0 apcosa) P0 1 I ) 0 (4.56) l (4.57) _w . (4.58) 2 66 IV.l.4 INFLUENCE FUNCTIONS OF THE THIRD GROUP For this case: W1 = Géég 43: W3 = GiB; and W4 = G53: (4.59) Using r2 = (x0+pcose)2+(y0+psin6)2 the following relationships can be obtained: _ P(2) ._ (2) 4 1 Po 0 2G2,0' 2“. I 2 I _ G(2) w2 + w3 — G1 1' (4.60) Consequently one needs to evaluate W1, W2 and the two I (2) and 3(2) 2 0 1 1 only. Using equation I I auxiliary integrals G (4.15) one obtains: (2)__ 2 (2) 2 (2)+ H(2) G2,0 ” “Kc {OX0 P0,0 + 20 x0 H0,1+ Ho,2} (2) _ _ (2) 2 (2 ) H(2) G1,1 ‘ Kc {pxoyb P0,0 + p X0H1, o + p 2yoH o, 1 + P3 H{?1} (4.61) 67 Fran (4.12) one obtains: _ _ 4 (3) 2 3 (3) 3 2 (3) W1 ‘ Kc {9X0 It30,0 + 4‘3 xo H0,1 + 6" Xo H0,2 4 (3) 5 (3) + 40 x0 H0'3 + p HO,4 (4.62) _ _ 3 P(3) 2 (3 ) H(3) W2"Kc{‘”‘oyop,oo W3OX0Y001+30X0Y0 4 3 2 3 3 3 2 3 +9Yoflé3+p:oH1(,()+3p X011191 4 (3)+ (3) + 30 x0 H]. 2+ 5H0'3} , (4.63) P(3) = H(3) = -{d0cos(¢)/f(¢) + 2a0 Péfg - dO Pizg}/El- (4.64) Accordingly all integrals in the right-hand side of equations I (4. 62) and (4.63) are derivable from the fundamental integral, P(l) P,0 0 Subroutine "'ICRVW" (see Appendix C) is developed for this method. 67 Fn:n(4.LD one<1fi2dns: __ 4 (3) 2 3 (3) 3 2 (3) W1 ‘ K {”0 P0,0 + 49" xo H0,1 + 6“ Xo H0,2 4 (3) 5 (3) + 40 x0 H0,3 p H0,4 (4.62) __ _ 3 P(3)+ 2 (3) H(3) W2 ' KC {”0 yoP o, 0 +23" X0 Y0 H0,1 + 393 x0Y0H o, 2 + (34 + 13:. + H13. 4 (3)+ (3) + 30 x0 Hl 2+ p5 110,3} , (4.63) where (3) (3) _ _ (2) P(2) P0,O — Ho'o — {docos(¢)/f(c)) + 2aO P0,0" dO P1 ()}/E21 (4.64) Accordingly all integrals in the right-hand side of equations I (4. 62) and (4.63) are derivable from the fundamental integral, P(1) P0, 0 Subroutine "TCRVW" (see Appendix C) is developed for this medmfl. 68 IV.2 NUMERICAL RESULTS Before considering examples, the accuracy and cost of the present method is compared to the method of Romberg- Richardson. Consider the coefficient matrix of a system of linear equations which is an n block of 2x2 matrices each obtained by integrating the influence function corresponding to a field point 5(5) and source point x’(s), D8” = £3 (NM)B.Y (Sl,S’) ds’, i,j=l,...n (4.65) 2) BIY=1121 i#j° For a segment of an elliptic region (see Fig. 4.11) corres- ponding to i = 1 and j = 8 (say) the results are obtained and compared to the Romberg-Richardson method. See Table (4.1). 69 .mapfimmn mmmmmuocfl ponume compumaoflmumumneom may now confiswmu mEHu no on» .womusoom once Mom comm mov.o mmH.o «Ho.o .wmm mo mmmommammmoo.o mmoammammmoo.o mmvommammmoo.o man mhvmmmmwmaao.o moammqumaao.o msvmmmmwmaao.o Hmo msvmmmmvmaao.o moammmmwaaao.o mnvmmmmvmaao.o man mmvmmmoaoooo.o mammmaomoomo.o mmvmmmomoomo.o Han acousood mica mom >0mnsoo< mloa mom >H .9650 m0 xflnwmz acmpnmnowmlmumneom comoumsoflmlmumnsom oonumz oaumamcfi uCWEmHm 0H": s @HHC 17me 12.me "4.4221 .conumz comoumnofimumumneom may sue: ponumz ucmmmum on» mo homuooofi can umou may mcwnmmfioo H.¢ mqm<fi 70 IV.2.1 EXAMPLES In this section, some numerical results are examined for a few test problems. The following regions with curved boundaries are considered: Region 1 - A limacon generated by Z = a-2bsinw, a>2b, see Figure (4.9). Region 2 - A region generated by Z = a+coszw, see Figure (4.10). 2 Region 3 - An elliptic region generated by Z2 = I:3355§7W e SYZ/r4, (qo + qls) Y3/r4, 1} (5.5) Substituting for X and Y from (3.18) and decomposing each component leads to: _ 3 2 W1 ‘ q0 W1 + q1 {a0 Q1,2 + 3a0 a1 Q2,2 2 3 + 3a0a1 Q3,2 + a1 94,2} _ — 2 2 W2 — q0 W2 + q1 {aob0 01,2 + (aobl + 2aO albl) Q2,2 + (a2 b + 2a + a2b 1 0 Oalbl) Q3,2 1 1 94,2} 84 W3 = aoqo Q0,1 + (aoqi + alqO) Q1,1 + aiqi Q2,1 ‘ W1 W4 = bqu Q0,1 + (qul + biqo) Q1,1 + blql 92,1 ' W2 w = 1. (5.6) V.2 VARIABLE DENSITY OVER A CURVED BOUNDARY A variable density of the form M(6)=c%)sink8+-coske 1&5 considered, where qO and q1 are constants re— lated to geometry and k is an integer. Furthermore,k = 0 corresponds to a constant load and K=an denotes no elastic deformation, due to the fact that sin m =cos—— *X () ~s;x' Figure 5.2 Linear over a curved segment. variation of source density 86 It is clear that sin (e - ej-l) and sin (e - ej) are of the form qo Sin 6 + ql cos 6, in which q0 = cosej_l, ql==-Sinej_1 amui q0 = cos ej, ql = — sin ej respectively. Letting qO sin 6 + q1 cos 6 = §;(6), the general components of the augmented kernel of the second group are: E = {g(e)(x0 + pCOSB)3/f%6) , g(e)(x0 + pcose)2 2 (Y0 + psine)/f(e) , g(e) (x (5-9) 0 + pcose) l... 2 (yo + psine)2/f2(6), g(e)(y0 + psine)3/f(e) ,l}T. I 0 0 $1119. Decomposing each component leads to the following where f(9) = a0 + b cose + c integrable functions: 3 (2) 2 2 (2) 3 1 quo p H1,0 + 3q0X0 p H1,1 + ino p H0,1 2 2 (2) 3 (2) 3 (2) + 3qlx0 p H0,2 + 43q0x0p Hl,2 + quxop H0,3 w = 4 (2) H(2) + q0p H1,3 + q1p4 H0, 4 _ 2 (2) (2) W2 ‘ qu0 yo0 H1, 0 + X00 2(2q0y0 + ino) H1,1 (2) 2 2 <2) 2 (2) + q1X02Y0p H0,1 + qu 0 p H2,0 + 2q1X0Y0p H0,2 (2) H(2) + p 3(qoyo + Zino) :1, 2 + 2q0X0p3 H2, 1 + q0p4 H222 + qu00 3H623 + qip4 Hi2; w3 = qu 0p P310 + in 0p Héli + qo02 Hili + q102 H012 W1 W4 = q0Y0p Pélb + qiyop H611 + q102 “iii + H(l) q192H2,0 W2 87 .As was shown in Chapter IV; all of the H(K) inte rals in m,n g the above equations are derivable from the fundamental (l) integraLyPo'o. CHAPTER VI CLOSURE The solution method which has been described in this dissertation is the result of a continuing study on the be- havior of the components of fundamental solutions in the theory of elastostatics. As a result of this study, a general formulation which includes several versions of the so called "indirect approach" was formulated. Accordingly, these fundamental solutions were classified in appropriate forms, leading to the development of a corresponding numer- ical scheme which applies to straight segments or any boun- dary which is replaced by straight segments. Due to the fact that the entries of the coefficient matrix of the Zn linear system of equations are obtained from the integration of the components of the influence functions over individual segments, the straight segment approach may prove to be too crude for curved boundaries. Thus for curved boundaries, discretization was replaced by dividing the boundary into a finite number of arcs. Analog- ous to straignt segments, an analytic scheme for curved segments was developed. It was shown that, in any case, the components of the "indirect" approach always lead tc»exact integrations. Accordingly this makes the technique both efficient and independent of numerical integration. 88 89 For displacement boundary value problems, a suitable linear combination of the force dipole and the dislocation dipole leads to this type of formulation. For traction problems, the single force or the Massonet singularity is used. In chapter V the idea of analytic integration was exten- ded to suitable singularity variations over individual segments. Based on the formulation explained in this dissertation, a compact computer program for arbitrary influence functions of linear elastostatics was developed in which all of the aforementioned schemes were incorporated into)Sub- routines for minimum storage and computational effort. The program is self sufficient in that it uses no computer library subroutines. Based on the results obtained, the computer program requires a minimum CP execution time. On the CDC 6500 computer, it was shown that this time varies from less than 0.7 seconds for simple geometries and boundary conditions with 40 segments, to less than 3.7 seconds for a complicated geometry, variable boundary conditions, and 60 segments. It is important to note that the major portion of the CP time for these cases is consumed in determining the den- sities of singularities (i.e. in solving the Zn system of linear equations). Computation of the elastic fields takes very little of the total CP time. The program is coded in FORTRAN and is included in the Appendices. APPENDICES APPENDIX A DERIVATION OF SOME IMPORTANT INFLUENCE FUNCTIONS The fundamental influence functions, (uR)B.Y, (uC)B.Y (HR)B.Y and (HC)B-Y are given in (3.3) through (3.6). In- fluence functions related to GY' FY and their dipoles, i.e. gY and fy, can then be obtained from the following relation- ships [19]: (NG)B Y = (NC)B-Y + 20 g? e6y (NR)B.6 (A.l) (NF)B-Y = (NR)B-Y + 5% g: e5Y (NC)B,5 (A.2) (Ng)8.Y = (NC)B-Y + 2n ;; 86y (Nr)8.6 (A.3) (mos.Y = (was.Y +.g: éi edy (NC)B,5 (A.4) Influence functions related to stresses “08 are easily obtained from Hooke's law: _ 2 (HM)aB.y—u(6a>1 686 + 58A 605 + a 6A6 6018) axA (uM)A.0 (A.5) where M is any of R, F, G or C*. For dipoles: = , NM A.6 (NIII)B.Y 3x ( )Boa ( ) * Similar relationships exist for mypthe dipole of MY' 90 91 Linear combination matrices for any other influence functions can be evaluated using (A.l) through (A.5). For examplelthe corresponding linear combination matrix related to the displacement due to a dislocation dipole, i.e. hufls is obtained from (A.6) and (3.4): (V5+2)n’l (V5+2)n 2 V2 n 1 V2 n 2 0 A _ _V -v2 n’2 (V5+2)n’l (V3+2)n’2 v2 n 1 0 .1.‘ 1 , , , _ , V2 n 2 (V3+2)n l (V5+2)n 2 V2 n 2 0 L V2 n’1 V2 n’2 (V5+2)n’l (V5+2)n’2 OJ I For (ug)B the matrix is: V6nl V6n2 V2nl V2n2 O 0 V n' V n' 0 0 A: l 31 32 (A.7) "’ TTV5 , . 0 V3nl V3n2 O 0 vani V2n2 V6nl V6n2 0 " I I 2 source point x'. where ni and n are components of the unit normal at the 92 A.l SPECIAL CASE In the numerical solution of many problems of elasticity, the effect of Poisson's ratio is relatively unimportant. In plane traction boundary value problems without body forces, Poisson's ratio has no effect on the stress distribution at all. If one then chooses v=.5, for which v2=v-1=1;V - 1 = 0, the following useful relationships are obtained from (A.l): (NG)B.Y = (NC)B-Y (A.8) (NF)B.Y = (NR)B.Y (A.9) (Ng)B.Y = (NC)B Y (A.lO) (Nf)B.Y = (Nr)8.Y (A.ll) where N8 is either 118 or H8' The above equations are also valid for N8 replaced by the stresses H 8' Cl APPENDIX B THE ROMBERG-RICHARDSON METHOD APPLIED TO BOUNDARY INTEGRALS 6B It is required to evaluate I=%A!~F(e)de, where F(6) is chosen according to the cases outlined in this dissertation. - . . 9 -9 USing the composed trapeZOidal rule With 091 = Em A where m is the number of subdivisions of the interval (9A,6B): A6 m-l 11 = —12 {F(6A) + F(eB) + 2 El F(6A + 1801) + ()(A61)2 (B'l) is an approximation to I. An approximation to I can also be determined with 861 halved, i.e. with A62 = 3%; : A6 2m-l 12 = —72 {F(6A) + F(eB) + 2 iél F(eA + 1882) 2 (B.2) + C)(Aez) . The first term of the error series can now be eliminated by taking a suitable combination of I1 and 12: I ~ 412 - 11 (B.3) ~ 4 _ 1 . One can see that the leading error term is now'O(Ael)4. Assume for the sake of simplicity of notation, that m is K a power of 2, i.e. m = 2 and let the approximation given 93 94 by (3.1) and (B.2) be deSignated by T0,k and T0,k+l respectively, i.e.: 4T -T 1 0,k+1 0,k I ~ 4 _ 1 . (8.4) Designate the right-hand side of (B.4) by Tl,k' Successive halving of the interval will give a sequence of values T0,k and each successive pair can be combined to give the values Tl,k' The sequence of values Tl,k can then be com- bined in a similar manner to remove the leading error term C(89)4 according to: = 3' - j_ Tj’k (4 Tj_1,k+l Tj_l’k)/(4 l), (3.5) j = 1,2,... k = 0,1,2,... This process can be continued to form a sequence of columns with error terms of increasing order. Using the fact that AB = (BB-6A)/2k, one can see that the error term for the approximation Tj,k values in successive columns converging rapidly to the true is of order [(eB-GA)/2X] (2j+2) with.the values of the integrals. To generate the computer program, the integrals are written as follows: INTEGRALS OF e:‘HE FIRST GROUP: = 3 W1 ZfeAKn f(e) d9 W2 = pfegflxo + pcose) (y0 + psin6)/f(e)} de 6A 95 6B 2 W = 0 {(x0 + pcose) /f(e)} de A 9B W4 = p f {(y0 + psin6)/f(6)} de A 9 Y + ' = B An pSlne W5 0,] e arctg {x0 + pCOSB} d INTEGRALS OF THE SECOND GROUP: 8 (3.6) 9 w = p .I'B (x0 + pcose)3/f2(e) d8 l 1 9A 9 2 W2 = p ./.B (x0 + pcose) (y0 + psin6)/f%e) de 6A ‘ 6 . 2 2 W3 = p .[.B (x0 + pcose)(y0 + pSine) /f(e) de 6A 9 . 3 2 W4 = p ./.B (y0 + pSine) /f(6) d6 (3.7) 6A INTEGRALS OF THE THIRD GROUP: 6 W1 = p Jf 3 (x0 + pCOSG)4/f3(9) de 6A 9 3 . 3 W2 = p f B (x0 + pcose) (v0 + pSine)/f (6) de BA BB . .3 3 W3 = p Jf (x0 + pcose)(y0 + pSine) /f (e) de 6A 6 W4 = p Jf B (yo + psine)4/f3(e) de (3.8) 6A Where f(8) = a0 + b0 cose + c0 sine and a0, b0 and c0 are given in (4.15). Subroutine RMBGW is developed for this method, (see Appendix C). A function statement F(i,X0,yO,Rk,Z) is used for the integrals to be evaluated over a curved segment. APPENDIX C COMPUTER PROGRAMS Two general computer programs are developed: (i) Program BIEM This program is based on the formulation of Chapter III, i.e. when the boundary is replaced by a finite number of straight segments. (ii) Program TCURVE This program is based on the formulation of Chapter IV when the boundary is replaced by a finite number of of arcs. A. INPUT DATA The following information must be provided as input N Total number of subdivisions on the boundary. PR Poisson's ratio for the material. G Shear modulus of the material. IP IP = 0 for plane strain, IP = l for plane stress. IB IB 0 for displacements prescribed all over the boundary, IB = l for traction prescribed all over the boundary. IBMX Indicates the number of nodes at which displace- ments are prescribed. 96 INUM NSB IW NF X(2N+2) IBC (2N) B . AUXILIARY NOTAT ION XN(2N) C(10,6) Tl’TZ Xo'yo RK,K II‘A'TB DA(N,6) 97 Integration factor. If INUM = 0, the integrals are evaluated analytically. If INUM = l, the integrals are evaluated numerically using a nine point Chebyshev quadrature formula. Number of subdivisions for each segment, (for the case INUM = 1 only). If IW = l, the values of the fundamental integrals a will be printed. XL“ Number of internal field points. Coordinates of the boundary segment end-points. ~~~ Subsequently stores the solutions. Type of boundary conditions; i.e. IBC = 0 implies that a displacement is prescribed and IBC = 1 implies that a traction is prescribed. Stores coordinates of boundary segment mid-points. Stores linear combination matrix elements in C(10,5) and stores integration nodes (roots of Chebyshev’s polynomial) in the last column. Components of unit outward normal to the boundary. Center of curvature of each curved segment. Radius of curvature and convexity factor for each curved segment. Angles 6 defined in Figures (4.4) and (4.5). A'eB Stores (x K,Rk,TA,TB), with each of these olYol variables occupying one row. D(2N,2N+l) T(20,20) 98 Stores coefficients of the linear system of equations in the first 2N columns. Stores boundary conditions in the last column. Stores extrapolated values from the Romberg- Richardson technique. SUBROUTINES INT TCRVW RMBGW DTA TR UTR MTX Calculates integrals over straight segments using the analytic method if INUM = 0 and using the numerical method if INUM = l. Calculates integrals over curved segments analytically. Calculates integrals over curved segments using the Romberg-Richardson technique. Replaces an arbitrary (curved) boundary by a finite number of arcs and calculates x0,y0, Kc,p,0A, and 6 for each segment and stores B them in DA(N,6). Stores the linear combination matrix for (HR)B and (HR) For v=.5, the results are also a8.y' valid for (HF)8.Y and (HF)aB.Y' Stores the linear combination matrices for (uR)B.Y ’(HR)B.Y and (HR)GB.Y for a general mixed problem. For v=.5, it also contains the corresponding (HF)B.Y' (HF)B.Y' and (HF)o8. values. Solves a system of 2N linear equations using Gaussian elimination with partial pivoting. 99 It should be noted that: (i) If only stresses are required, the smaller sub- routine SINT is used in place of INT; (ii) If tractions are prescribed all over the boundary the smaller program TRAC is used in place of BIEM; (iii) For boundaries composed of both straight and curved segments, program TCURVE employs both subroutines SINT (for straight segments) and TCRVW (for curved segments); (iv) NSB, IW and INUM are used in the subroutines INT and SINT only: (v) To avoid confusion, integrals of the second group are designated by W6, W7, W8 and W9. (vi) Since for each field point the solution ul, u2, H11, H12 and n22 is stored in X(2N+2), make 2N+2:5NF. D. OUTPUT DATA N, PR, G, IP, IBMX, INUM, (see INPUT DATA). ll’ H12, n22 at internal field points. 00 00 *E C 200 400 450 500 550 600 650 100 ***§******fié***§§*******§******§********************************** PROGRAM BIEM(INPUT.OUTPUT.TAPE5=INPUT.TAPE6=OUTPUT) *#**§*s****+****************************************************** DIMENSION X(122).XN(120).IBC(120).F(24).C(10.6).D(1201121) NAHELIST/PROP/N1PR.G.IP.IB.IBMX.INUM READ(5)*) N.PR.G.IP:IB:IBMX:INUM:NSB)IN 5 NRITE(6:PROP) READ(5)*) (X(2*I*1):X(2*I):I=1:N) 5 X(2*N‘1)=X(1) 5 X(2*N+2)=X(2) READ(5.*) (IBC(2*I-1)aIBC(2*I):I=1:N) READ(5:*) NF.(F(2*I 1):F(2*I))I=1:NF) 5 N2=2*NF 5 N5=5*NF IF(IP.EG.1) PR=PR/(1.+PR) 5 PI=3.14159 26535897 5 NN=2*N 5 M=NN+1 DO 10 I=11N 5 XN(2*I-1)=(X(2*I-1)+X(2* I+1))/2. N(2*I)=(X(2*I)+X(2*I+2))/2. CALCULATION OF ELEMENTS OF D(2N)2N+1) MATRIX DO 20 I=1pNN 5 IBC(I)=1 5 IF(IB.EG.O) IBC(I)=O COIJTINUE 5 IF(IBNX.EQ.O) GOTO 30 DO 600 I=11N 5 12=2*I 5 Il=12-1 5 IO=IZ+1 5 I4=I3+1 S= SGRT((X(13)-X(I1))**2+(X(I4)-X(12))**2) XF=XN(II) 5 YF=XN(I2) 5 T1=(X(I4)-X(I2))/S 5 T2=(X(Il)- X(IS))/S DO 600 d=11N 5 J1=2*d-1 5 J2=2*d 5 J3=J2+1 5 J4=J3+ 1 XA=X(J1) 5 YA=X(J2) 5 X8=X(J3) 5 YB=X(J4) 5 IF(I. E0. J) GOTO 200 IF30 5 CALL INT 5 YA=X(J2) 5 XB=X(J1+2) 5 YB=X=.1679061842148 5 C(5.b)=0. 5 C(4.b)=-C(é,é) 5 C(3,é>= +*C(7.6) 5 Cl2.6)=—C(B,é) 5 C(1.b)=-C(9.6) DX=(XB-XA)/NS: 5 DY=(YBwYA)/NSB 5 DO 35 I=1.NSB XA=XA+(I—1)*DX 5 XB=XA+DX 5 YA=YA+(I—1)*DY 5 YB=YA+DY AX=