THE UNIVERSITY OF MI CHI GAN COLLEGE OF ENGINEERING Department of Aerospace Engineering Technical Report ON THE SOLUTION OF PLANE, ORTHOTROPIC ELASTICITY PROBLEMS BY AN INTEGRAL METHOD Rafael Benjumea David. L. Sikarskie ORA Project 03616 supported by: NATIONAL SCIENCE FOUNDATION GIANT NO. GK-17101 WASHINGTON, D.C. administered through: OFFICE OF RESEARCH ADMINISTRATION ANN ARBOR August 1970

V V ''j7 Y~~ / bi

ABSTRACT The present paper is concerned with the application of integral equation techniques to problems in plane orthotropic elasticity. Two approaches for solving such problems are outlined, both of which are characterized by embedding the real body in a "fictitious" body for which the appropriate influence functions are known. "Fictitious" tractions are then introduced such that the boundary conditions on the real body are satisfied. This results in a coupled set of integral equations in the fictitious traction components. Once these are found the unknowns, i.e., stresses, etc., are found in a straightforward manner. The difficulty is in introducing the fictitious traction field such that the resulting integral equations are useful computationally, i.e., are Fredholm equations of the second rather than the first kind. A sufficient condition for this is that the fictitious traction field is applied to the boundary of the real body. The two approaches mentioned above differ in the choice of influence function used, in one case the influence function being singular in the field and the other singular on the boundary. A solution method already exists in the isotropic case using the boundary influence function [3]. An alternate formulation is presented using an internal influence function which is shown to have computational advantages in the anisotropic (orthotropic) case. To illustrate the methods, the stress field is found in a "truncated" orthotropic quarter space, under the condition of a given traction on the truncated surface, traction free elsewhere. This problem is of interest in certain ii

Rock Mechanics calculations, e.g., to a first approximation the stress field is that due to a rigid wedge penetrating a brittle, orthotropic elastic solid (prior to chip formation). iii

NOMENCLATURE a. - elastic compliances ij B - bounding surface of the real region B* - bounding surface of the fictitious region C. - constants, i = 1,...,8 D - interior of the real region D*- interior of the fictitious region E,E,E - elastic moduli in x,y,z directions xx yy zz F - body force field (vector) f =f i + f j - real traction vector (i,j - unit vectors in x,y directions) x y f*- fictitious traction vector GxG,G, = shear moduli xy' xz' yz G.. - influence function (boundary); stress field due to a unit load in the mj;q q direction G. ij - influence function (boundary) for a moving fictitious region ij;q H. - influence function (internal) lj;q k - ratio of constants (Aj -) -2 L - length of truncated boundary in example problem m - form parameter n =n i + n j - unit normal vector x y p - scalar pressure magnitude p;p * - real, fictitious dimensionless tractions. P;PP - field oundary points 1'P P f,ield, boundary points iv

P,P - real, line load components in the ith interval xi yi P* P*. - fictitious, line load components in the ith interval xi' yi r,r = (x-x'),(y-y') - x,y lengths between boundary, field point x y x-x'y y -y' r,r - dimensionless x,y lengths between boundary, field point ( L ' L) x y L L s - distance along the boundary t = t i + t j - real traction vector x y t* - fictitious traction vector T = T i + T j- line load x y u,v - displacements in x,y directions i = r + Air (i =1,2), complex variable i x i y a- convergence parameter a1~,2 - real or complex roots - angle of deviation from a qarterspace ( = O - quarterspace) pij - elastic constants for orthotropic, plane strain 7 - angle of orthotropic symmetry from y axis - constants, see pi below (i = 1,2) 6. - constants, see p. below (i = 1,2) E.. - strain components 13 T = v2 y - stretched y coordinate 9 - angle of truncation 1 - angle of field point from the vertical = tan 4, coefficient of friction 4 1~2., j2 S- roots of a fourth order equation (complex) v,v,...,v - Poisson's ratios xy' yx' zy v

- coordinate along truncated boundary - angle of friction, also, stress functions Qi(zi) - analytic functions of the complex variable z. (i = 1,2) f - reference angle, fictitious boundary and y axis vi

1 I. INTRODUCTION In the present paper a study is made of the application of integral equations to plane problems in nonisotropic elasticity. An excellent general review of the plane problem has been done by Teodorescu [1]. The use of integral equation techniques in elasticity has been extensive, both for solving a variety of specific boundary value problems and in the development of general approaches valid for larger classes of problems. In the latter case the work of Muskhelishvilli [2], and Massonnet [3,4,5] is of special pertinence. In a numerical sense, casting the problem in a form requiring the solution of a single or a set of integral equations can be viewed as an alternate approach to such methods as finite element, finite difference, etc. It has advantages, however, which justify its use in certain problems over other methods. A major advantage is that the final system of algebraic equations is of lower order, i.e., the equations hold on the bounding surface of the region rather than in the field. This is offset somewhat, however, by the fact that the resulting matrix of coefficients is dense and techniques available for the band (and sparse) matrices which generally appear in the finite element, etc., approaches are usually of little help. A major problem is formulating the integral equations themselves, in particular, deriving equations in a form suitable for numerical computation, i.e., Fredholm equations of the second kind. The utility of the integral equation approach is, in fact, highly dependent on the details of this formulation. This is illustrated in the present paper in which two ways of formulating integral

2 equations are presented, both valid for general elastic material behavior, twodimensional geometries, and traction boundary conditions. Both approaches involve "embedding" the actual body in a fictitious body with a known influence (Green's) function. These influence functions are then superimposed in such a way that the appropriate boundary conditions are satisfied. The approaches differ essentially in the choice of influence function or singular solution used. While both lead to well defined numerical results, there are decided advantages to one or the other depending on whether the material is isotropic [3] or orthotropic. A similar approach (a fictitious surface deflection is introduced) has been used by Wu and Chiu [6] for the solution of the contact problem of layered elastic bodies. Also, it is of interest to note that a similar method has been used for melting and solidification problems in heat conduction [7,8]. Detailed derivations for both approaches are presented in Section II for a general problem. The specific problem used to illustrate the method is the determination of stresses in a rock mass consisting of a "truncated," orthotropic, quarter plane (see Fig. 5). Of ultimate interest in this problem are stress distributions sufficient to cause fracture although only the stress distributions are discussed in the present paper. This is done in the third section. The details are summarized in the last section and possible extensions to other problems such as displacement, mixed boundary value problems, and nonhomogeneous materials are discussed.

5 II. DERIVATION OF THE INTEGRAL EQUATIONS Consider a properly supported two-dimensional region, having an interior D and a bounding surface B, shown schematically in Fig. 1. In outlining the method we begin with the simplest problem namely, the traction specified everywhere on the bounding surface of a simply connected body. Comments concerning the displacement and mixed boundary-value problems follow in a. later section. Two singular solutions are now identified which will be useful in superposing solutions to more complex problems. In Fig. 1 we consider a line load T per unit length which acts at some point P1 on the boundary. The stress at some field point P for this loading can be represented by*: a ij(P) = Gij;q(P, P1) Tq(P1) i, j x, y (1) where G.. (P, P1) is an influence (Green's) function, i.e., it represents the mJ;q 1 stress state at P due to a unit line load at P in the q direction. For G.. (P, P ) known the stress field for any traction distribution t on B (see ij;q 1 Fig. 1), is given in the form: (P)= f G (P, P1) tq (P1)ds (2) B ij;q 1 q 1 Similarly, a line load T acting in the interior of the region D as shown in Fig. 2 can be considered. In this case the stress at some field point P can be written, similar to equation (1): *Repeated indicies imply the summation convention.

4 ij(P) = Hij;q(P, P1) Tq (P ( where H.. (P, P ) is again an influence function and represents the stress 11;q state at P due to an interior unit line load at P1. Knowing the influence function H.ij;(P, P1) the stress state due to a body force field F is given by:. (P) = Hi (P, P1) F (P1) dV (4) JjD 1 q 1' We will make use of equation (4), however, in the following way. Consider the following body force field (see Fig. 2): F in D' F = F = (5) elsewhere Equation (4) becomes a. (P) = H. (P, P) F (P) dtds (6) B A 3j;q 1 q 1 B' At Defining the following limit, lim / Fq(P1) dt f (P1) At At + 0 F + o q equation (6) thus represents the stress at son "traction" distributed over the surface B', i.., oj(P) = H;(P f P1) f (P1) ds (7) ijB;q. 1 q 1 The functions Gij, H.. are general at this point in the sense that all field equations, including the generalized Hooke's law, plus the boundary

5 conditions are satisfied. Note that the two singularities introduced are fundamentally different; one appearing on the boundary, the other in the field. Except for a few regions of simple geometry the functions Gi. and H.. are ij;q ij;q not available, To solve boundary value problems for regions in which these functions are not available, we consider the following. The problem to be solved consists of a region D with traction t specified over the bounding surface B (see Figo 3). The real region D is now embedded in a "fictitious" region D* for which we assume G.. and H are known. Each of the functions leads to an independent 1o;q ij;q formulation of the problem. To formulate the problem in terms of G. we inij;q troduce a fictitious traction t* on B*, as shown in Fig. 3. The stress field in D* is now given by equation (2). The boundary value problem in D is now satisfied except for the boundary conditions on B. These are satisfied as follows n = t.; on B (8) where nj, ti are the componenets of the unit normal and traction vector respectively. Substituting equation (2) into (8) we obtain, f G (P, P) t* (P ) n (P) ds = t.(P); P on B (9) B* ij;q 1 q 1 j B* Equation (9) represents a Fredholm integral equation of the first kind in the unknown fictitious traction t*. If equation (9) can be solved, our original problem is solved with stresses given by problem is solved with stresses given by

6 ij(P) = G (P, P) t* (P ) ds; P in D (10) This is the complete solution since all conditions of the original boundary value problem are satisfied. Note that a stress field is obtained in the entire region D*, but, of course, only has physical meaning in D. Due to difficulties in dealing with Fredholm integral equations of the first kind [9], however, equation (9) is not in a useful form. Sufficient conditions will be shown by which equation (9) can be reduced to a Fredholm equation of the second kind, for which extensive results are available [9]. We now formulate the problem in terms of the function H.. Consider an ij;q unknown fictitious internal traction f*, as previously defined, applied on the bounding surface B within the fictitious region D* as shown again in Fig. 3o The stress field due to f* is given by equation (7) and again satisfies the original boundary 'value problem except for the boundary conditions on B, i.e., equation (8). Substituting into equation (8) we obtain f H. i;(P, P') fq (P') nj (P) ds = ti(P); P on B (11) B iJ3;q q 1 This equation is similar to equation (9) with one important difference, namely, the kernel H.. (P, P') is singular when P, P' coincide. Because of this equaij;q tion (11) can be cast in an alternate form, i.e,, f H f* n ds = f H. f* n.ds + f H. f* n ds (12) ij; q q AB ij;q q j Bj;q q j B AB- ~B where the element of srface B contains thae point of coincidence of P, P' It can be shown that,

7 ff* lim J H.i; f* n.ds - (15) AB )ij;q q j 2 AB -+ 0 where the factor 1/2 in the limit comes from the fact that the load is applied internally. Substituting back into equation (11), the final form is f* (P) -+ f H, (P, P') f* (P') n. (P) ds = t(P); P on B 2 - i;q q where the integral is evaluated at all points on B except where P = P'. Equation (14) is a Fredholm equation of the second kind, particularly suitable for numerical evaluation. The reason for the different types of integral equations, i.e., equations (9) and (14), comes from the fact that in the latter case the fictitious traction is applied on the real boundary. This generates a singularity, the limit of which is just the dependent function itself, see equation (13). A sufficient condition, in fact, for generating Fredholm equations of the second kind in this approach is that the fictitious traction introduced is applied to the real boundary. Integral equations of the second kind can also be generated using the first approach, i.e., using the G.. functions, in the following way. Consider, again, two regions; the real region of interest D embedded in a fictitious region D* for which the function G.. is known. We further require that the 1ij;q boundaries of these two regions B, B* be tangent at some point, i.e., such that the point P1 on B* coincides with P" on B as shown in Fig. 4. At this point we

8 apply an unknown fictitious line load t* ds (ds is an element of length of the boundary B). This creates a stress field in the entire region D*. da.ij(P) = Gij.(P P1) t (P) ds (15) where P is a field point. At this point no boundary conditions on B have been satisfied. We now consider a second point P"' on the boundary B. The region D* is now translated and rotated such that the points P on B* and P'" on B are tangent. We again apply a fictitious line load at P'" giving an additional contribution to the stress field at P similar to equation (15). It is important to note that the coordinates in the function G.i have been translated and rotated*. This translated form of G. is referred to as... If this is lj;q 1ij;q done for all points on B, then the stress.field at P is. (P) I;q(P' P ) t* (P ) ds (16) B ij;q 1 q 1 This stress field satisfies everything except the boundary conditions. Substituting into equation (8) and again extracting the singularity that appears when P, P are coincident [3], we obtain; ti (P) + j i;q(P, P1) t* (P) nj(P)ds = ti(P); P on B (17) 1 B ij;q q 1 q 1 *L Again, equation (17) represents a Fredholm integral equation of the second kind. This approach has been outlined by Massonet [3] for a specific choice of fij;q namely, the field due to a line load on the isotropic half plane. Again, once *The coordinate system x, y is fixed in the real body, hence as the fictitious region changes position its coordinates must be transformed.

9 the fictitious traction t* is known, the solution is complete with the stress field given by equation (16). Two approaches have outlined above which lead to Fredholm equations of the second kind, i.e., the equations (14) and (17). In the form presented both are valid for a generalized Hooke's law. In the isotropic case both approaches are of comparable numerical effort since the elastic constants are independent of the coordinate system.. In the nonisotropic case, however, difficulties arise in generating the function G.. since for each position of the fictitious region the elastic constants must be transformed. No such problem exists for the singular solution Hij since it remains fixed with respect to the embedded real iJ;q body. In Appendix A specific examples of the functions G.. and H.. are tabuij;q ij;q lated for both isotropic and orthotropic materials. In the following section a. specific example is solved (for an orthotropic material) making use of these functions. The problem is formulated in terms of the internal load approach, i.e., equation (14). The alternate approach, equation (17), is formulated in Appendix B and the results of the two solutions discussed in the following section.

III. EXAMPLE PROBLEM A. INTEGRAL EQUATIONS The example considered is that shown in Fig. 5. The surface AC and BD are traction free while the surface AB has a specified traction t. The material is orthotropic with one axis of symmetry perpendicular to the plane of the figure (plane strain is assumed) and the other two defined by the angle y, y + 90. This problem is of interest since the stresses represent to a first approximation those due a rigid wedge penetrating an orthotropic half plane (stresses just prior to chip formation [10,11]). The traction on the surface AB is assumed to be of the form (see Fig. 5): t = p(/L)m l-cos2it V/L); 0 < < L t = tn= t p(t/L)n (l-cos2 tn/L}; 0 < t < L (18) t = t =0; elsewhere n s where p is a scalar magnitude and m is simply a form parameter. t and t are n s the normal and shear stress distributions on AB with p representing a coefficient of friction. There is no loss of generality in assuming a specific traction distribution such as that given by equation (18). Two formulations for the problem were outlined in Section II resulting in equations (14) and (17) respectively. Equation (14) represents a new formulation while the basic idea contained in equation (17) has been given by Massonnet (although no solutions for orthotropic materials have been given). For the

11 present orthotropic problem both approaches have been carried through to numerical results for purposes of comparison. The first approach is outlined in detail in this section while the second, due to the complexities mentioned in Section II, is only discussed.* In the first formulation we chose as the singular solution the stress field in an infinite, orthotropic, two-dimensional elastic solid, due to a concentrated load acting internally. Thus the fictitious region D* in Fig. 3 is infinite. With no lose in generality we chose the axis x,y coincident with the axis of symmetry of the material (see Fig. 6). The singular solution H.. is derived 1ij;q in detail in Appendix A and is given by**: H f.* -c *1Cx (x-x')+C fy*-y) C fx(x-x' ) C4fy(y-y') C f*(x-x')+C f*(y-y') C f*(X-X')+C8f(y-y') H f* - 2 l 2y 5x Yx;q q 2A(a2-a ) (xx 2 + (yl)2 (x- 2 + (y_)2 xy;q q 2(a -a) L2(x-x.)')2 ( a(x-x) + (yy) traction f*[. The constants Cf - C*( are: H f*1 l 3 x 6yx 1xy 2 it 1 - 22 5 - 2 L1 22 The primed coordinates in equation (19) refer to the position of the unknown traction f*. The constants C - C are: *A detailed derivation of the equations is found, however, in Appendix B. **Actually, it is more convenient to present Hj f*. Also, for convenience;qthe function is presented in component form. the function is presented in component form.

12 1/2 12 5/2 2 =2"1 1 al PT C= 6 C C = 11/2 2 - a12 C a2 C3 = al / Ua2 - 12 C7 = 13/ 7 c4 = 12 l /2 = 3/2 The constants j.. are those that appear in the orthotropic, strain relations, e.g., E = xx +, 1see [12]. al',2 are related to the ij 's as follows, 212 + 66 a 2 + a2 P22 121 a11 p121 - 12 2 ll 1 plane stra The two co (20) in, stressnstants (21) 1o2 = p1/22 The reason for embedding the real region (that shown in Fig. 5) in the seemingly awkward way, as shown in Fig. 6, is that we want the axis of material symmetry in the real region to also be coincident with the x,y axis. This results in the simplest form of the orthotropic stress strain equations. The problem now reduces to finding the fictitious traction f* along this boundary such that q the correct boundary condition is satisfied, i.e., equation (18). The following dimensionless variables are introduced, x-x' r x L r y L p* = q q P t q p = q q P S = s/L (22)

13 Substituting equations (18), (19), and (22) into equation (14), the final form of the integral equations is, P~~ C C x 1 1 x+ 1 -fp*(nr +n r 1. 2 22 2 2 2 x 2 2y 1 Qr+r cxr +r DBA 2 x y 1 x y * y,/C r n - C r n C r n - C r + - 2 yxx * 1 6x y x x y x 7 x y 2 22 2d = (23) rn - Cr n Cr n-Y x n C C d -p*(n r + n r dS p y xx yy 2r2 2 2 2 y 2x y lx.y. where nx, n are the components of the unit normal to the surface DBAC. p x y x and p are the dimensionless components of the traction in the x,y direction. Using the relation, = tan O, and equation (18), p,p can be written, xy 0; on DB = (/L)1 - cos 2r(/L)] sin (e +; on BA r; on DB P = -(/tL)m [ - co1o s 2 n(L) cos ); on BA; onos AC 0; on AC

14 As mentioned previously, once the fictitious traction is determined from equations (23) the stress field is found by an integration, i.e., equation (7). Equations similar to equations (23) have been derived for the second approach, see Appendix B, using as the singular solution the line load on the boundary of an orthotropic half plane. Numerical results have been obtained for both approaches and are discussed in the following sections. B. NUMERICAL RESULTS Equations (23) are now evaluated numerically. The specific material of interest for numerical purposes is rock (granite) which is assumed to behave as an orthotropic elastic solid with the following elastic constants* [13]. E = 6.47 x 10 psi v =.168 xx xy E = 6.15 x 10 psi v =.176 yy yx E = 4.40 x 10 psi v =.141 zz xz 6 (25). G = 2.84 x 10 psi v =.209 xy zx G = 2.50 x 10 psi v =.146 xz yz G = 2.45 x 10 psi v.203 yz zy *An orthotropic material has nine independent constants, i.e., the following relations hold. v v v V v V rxy _yx. xz zx. yz zy E E 'E E 'E E yy xx zz xx zz yy

15 The relationship between the elastic constants given above and the 's in equations (20) and (21) can be shown to be (see Appendix D), 1 - vv 1-v v _ xz zx _yz z y 11 E 22 E xx yy (26) v + V V p = xy xzs z p6 1 12 66 G E xy yy It is of interest to investigate the case of 7 = 90~ (see Fig. 6). The other parameters are chosen to be: 3 = 0, 9 = 45~, m = 1, = 0 (only the normal component of the traction on AB is considered). Equation (23) can be evaluated numerically in a straightforward manner [14]. The boundary is divided into a finite number of divisions of equal or unequal length and a concentrated line load, the resultant of the traction on each subdivision, is assumed applied at the center of the division (see Fig. 7), i.e., P. = f p*pds P = f pds xi A xi s. As. px i i (27) P*' = f p*ds P f p ds yi As. yi s. y i i The trapezoidal rule is used to approximate the integral. This yields a system of algebraic equations to solve for the concentrated fictitious loads P*., P*.. X1 yi

16 +p 1 C C xi +1 x ( 1 2 + 2(-2 ) j=l Xj xi xij yiryj 2x 2 2 2 ("2 1 jX1 yJ r2 lrxij yj 1 Cr.n - Cr. n. C r n r..n: +y 2 yiaj xi 6 xij yi 4 yij xi 8 x yi i yj 2 2 2 A 2 i i xi r1 a r..+ r r.. ar X + a2 xij yij 1 xij yij J N r*_/Cr C rn Cr Cr n - C r n + 1 yij xi 5 1xij yi yij xi 7 xi i yi 2 2i(a -a) j=1 xj 2 2 2 2 / 2g(24) j oi x j di o r + r ij r + r ji L_ \ 2 xij yij 1 xij yij / C6 C8 -P* r (n r + n r ) AS yj xi xij yij a r2 2 + 2 i yi y y2rxij + y ijx (28) i=l,...,N Note that in the summation j does not take on the value i. N is the number of subdivisions counted as shown on Fig. 7. P. and P are found by intexi yi grating equations (24) according to definitions (27). r r * are x,y xij ' yij distances between the centers of the i,j subdivisions. The solution of the set of algebraic equations (28) can be done by either matrix inversion or iteration. Iteration was chosen since machine time required is less. Some difficulties can be encountered with the convergence of the process, however, these can usually be circumvented by introducing a convergence parameter a, i.e., see [3]. Equation (28) can be written in the general form,

17 N P* = 2P. -A A. P* + B..P*. x i x j= ij xj ij Yj (29) N P* = 2P - Ao Z D.P. + E..P* yi yi i xj iJ yj j=l where A, Aij, B ij, C., and D.i are found by comparing with equations (28). Following Massonnet [5], these equations can be iterated in the following form (n+l) _ (n) N (n) (n) p*(n+l) = p + (l-ai)P*(n- Ao D A P*) + B P* xi x i x i j ij yj j=l (30) (n+l) _ (n) N (n) (n) p.(* l) = 2aP + (1-)P*(n)- A o Ea D.P*( + E.P* yi yi yi.1 1i xj iJ yj j=l where the superscripts in parentheses are the number of iterations. The convergence parameter a lies between 0 < a < 1. The value a = 0. 8 was used in all computations and convergence was within 12 cycles for the initial choice (the prescribed traction field) and the desired accuracy (1Co) which is required at every point of the boundary. Once the fictitious traction is found from equations (30), the stresses can be found anywhere in the field from equation (7). Using equations (19) and (27), the numerical approximation for equation (7) becomes,* *Only one stress component is given, the others follow in the similar fashion.

18 a, 11N C P* r CP*.r + C P*.r. xxi _ 1 z lxjxij 2yjyij C jrxij+ C p 2irc-(aa) 2 2 2 r2 1 j=C1 l 2 xij yij 1 xij yij Equation (31) gives the normal stress component in the x direction at some field point i (i has x,y coordinates). The boundary points j are just those points at which the fictitious traction has been found from equations (30). Computer programs have been written both for the solution of equations (30) and for the computation of the stresses, equations (31), see Appendix E. Numerical results for the above parameters are given in Tables 1 and 2. Table 1 gives the fictitious tractions along the boundard DBAC. In dimensionless form the length AB is unity. Fictitious tractions were computed for 3, 5, 8 units along the boundaries AC, BD to determine the effects of the infinite boundary. In Table 1, results are for boundaries CA, BD 5 units long. Note that each unit is divided into four segments with the x,y fictitious traction components computed at each segment. Stresses are computed at the field points 1 - 15 as shown in Fig. 7. Results of these computations in the form of principal stresses and directions appear in Table 2. It is of interest to note that along the 45~ line from the origin, the principal directions are all close to 45~, see points 4, 7, 10. This is due to the relatively small anisotropy and the applied traction distribution considered.* *An additional case was run for a strongly orthotropic material;E = 14x106 psi E = 8xl06ps6 xy- 1.5xl06 XX psi,Eyy = 8x10 psi Gxy = 1.5x psi, other constants given by equations (25). As expected, much larger deviations (from 45~) of the principle directions of points 4, 7, 10 occurred.

19 Several additional comments on the numerical results are pertinent. Firstly, regarding the infinite boundaries; it was found that for field points of the order of L from the boundary AB little change in the stresses resulted for the inclusion of boundary points beyond 5L. This is shown in Table 3 for field point number ten, see Fig. 7. Similarly, the effects of mesh size are important. Mesh sizes of L/4, L/8, L/16 were investigated. It was found that for field points of the order of two mesh lengths from the surface, further subdivision caused little difference in the results. This is again shown in Table 3 for the field point number ten. As mentioned earlier, numerical results for the present example were also obtained using the approach given by equation (17), taking as the singular solution the line load on an orthotropic half space, see Appendix B. While this approach has advantages in the isotropic case, in the orthotropic case the required transformation of the material constant tensor for each boundary point makes the numerical computation considerably more cumbersome. Numerical results in both approaches were essentially the same.

20 IV. SUMMARY Two methods have been outlined above for the solution of two-dimensional stress boundary value problems. The methods differ in the general type of singular solution or influence function considered. When singular on the bounding surface of the region, the influence function is referred to as Giq and when singular internally (in the field) it is referred to as H.ij In order to generate Fredholm integral equations of the second kind, necessary for numerical purposes, it was shown in Section II that the fictitious tractions introduced must be applied to the real boundary. In order to do this for the Gijq formulation the fictitious region introduced must be shifted such that it is tangent at each boundary point considered, This is not the case with the H, formulation. For orthotropic problems this represents a distinct advanij;q tage since in the former case one must transform the material constant tensor as each new boundary point is considered. Clearly, this is not a problem for isotropic materials, however. In this case there is some advantage to' the G.. formulation since the singular solution used is simpler, i.e., the stress field due to a unit line load on an isotropic half space versus that due to a line load in an infinite isotropic space, which, in effect, represents the two simplest choices of G.. and H, respectively. ij;q ij;q Numerical evaluation of the integral equations is a relatively straightforward task. One advantage of the integral methods over finite element, etc., methods is that the fictitious traction and stress field are separate calculations. Once the fictitious tractions are found the stresses can be found

21 anywhere in the field without resolving the problem, subject to the mesh size, boundary length discussion in Section III. As expected, the numerical results indicate that the stress field is a strong function of the magnitude of the anisotropy. Several extensions of these ideas appear possible although no results presently are available. Displacement and mixed boundary value problems can be formulated provided the displacement field corresponding to the singular solution is available. Formulations for nonhomogeneous materials (voids, rigid, or elastic inclusions) can be formulated by introducing two sets of fictitious tractions as well as two ficitious regions. Clearly, the problem of a multiply connected region is a special case of the nonhomogeneous problem, i.e., a void. Lastly, all of this work might be extended to three dimensions. This simply increases the number of coupled integral equations, e.g., from two to three in the stress boundary value problem.

22 APPENDIX A TABULATIONS AND DERIVATION OF SINGULAR SOLUTIONS We tabulate in this appendix some of the singular solutions G.. and H.. that are available. We first present the isotropic case and then the ij;q orthotropic case.* The cases which are considered are shown in Fig. A-1. Some of these cases are available in the literature, if so, the source will be cited, if not, those cases will be derived in this appendix. a. ISOTROPIC CASE a-1. Line Load On the G t* = xx;q q G t* = yy;q q G t = xy;q q Boundary of the Half t* r + t* r 2 x x y y-2 -2 -22 x (r + r ) x y t* r + t* r 2 x x y y -2 r t — 2 -2 2 y (r + r) x y t* r + t* r 2 xx yyTC -2 -2 2 x (r + r) x y Plane (see Fig. A-i, a-l) where r = x - x', r y - y'. The unprimed coordinates represent a field x y point and the primed the position of the acting load. *For later convenience G. t* and H.. f* are given. ij;q q ij;q q

25 a-2. Line Load At a Point of an Infinite Plane H f* xx; q q = - 2 —J ----7. |f* r (a, r -2 -22 xrx -21 x 2i(r2 + r2) x y -2 - -2 + a2 ry) + f r (a3 r -2 -2 a3 r) - f* r ( a r + y y y 2 x 1 F*- -2 yy; q q - 2)2 -22 x x2 r - 2 x(r +yr) x y -27 - a2 r i 1 y+ al r) H f* xy;q q 1,, _, (-2r(r +-2r) x y x* x -, -2 r( a r y 1 x -2 - -2 + a r ) + f* r (a r 2 y x2 x where a1 = (3 - 2v)/(l - v) a2 = (1 - 2v)/(l - v) a = (1 + 2v)/(l - v) and v is Poisson's ratio. Cases (a-3) and (b-3) are discussed later. b. ORTHOTROPIC CASE b-l. Line Load on the Boundary of an Orthotropic Half Plane-z Axis of Material Symmetry (see Fig. A-l, b-l) Axis Is an G = RL( 2 2 G t* = -Re 1i ( 1 -_ xx;q q fi(l( - 2) ) 2 - t*. + 2/ X P1 41 42 z 2-2 2) ti YJ G t* yy;q.. q = -Re 1 ( i 1 2i(l l -[~ '~~ - t* 2) x [A 4 + 21 C122 t '..... -, G t* = xy;q q l 1 - l +tR e lfi ( 1 '- l2 ) Z1 t* + 2 x 21 C11 112 z 2/ 2/

where Re means "real part of" and where, = r + i r; i = 1, 2 1 x y l'' "2' 1' 2' are the roots of the following equations (il' i2 are the complex conjugate roots of p1, |2 respectively) [12], for plane strain*, 114 3 2 66) 26 + 311 ~ 2116 4 + (2'12 + 266 - 26 [L + P22 = 0 for plane stress*, all - a16 3 + (2a12 + a66)42 _- a26 + a 2 &11 16C'2 6626 12 0 where the..'s and a.'s are related to the elastic constants (see Appendix 00 ij C). b-2. Line Load in the Interior of a Orthotropic Infinite Plane H f* xx;q q H f yy;q q H f* xy;q q ff - T C f* r + C f* r 1 1 x x 2 y y 2Tr(a -2 - )2 2 2 (2 1) a 2 rx ry 2 x y - Cf* r + C f* r 5 x x 6 yy 2r(C2- -l2 -2 2 " - 1a) a ~2 r2 + r C'2 x y C f* r - C f* r 1 x y 6 y x 2itr(a- a) -2 - 2 2 1 a2 r + r L 2 x y C f* r + C f* r 3 x x 4 y yi -2 -2 a r + r 1x y C f* r + C f* r 7 x x 8 y yi -2 -2 a r + r 1x y C f* r - C f* r 3 x x 8 y y -2 -2 a r + r 1 x y *For the special case of the x, y axis coincident with the other axis of material symmetry we have, '16 = p26 = a6 = a6 = 0 16 = 2616 26

25 where = 1/2 12 53/2 a [2 C-1LB2 C,= 1 1 /2L l ] - 2 C = 12 1 -L 12 C8 a3/2 1 - a 2 2C~~ +2 4 a2 a12 + 66 1ll and l + C2 a= 1;l,212 2... 22 22 The singular solution (a-l) is found in many elasticity texts, e.g., [15] and will not be elaborated further here. The singular solution (a-2) can be found in Sneddon [16]. The orthotropic half plane singular solution (b-l) exists for several special cases. Conway [17] has the solution for a special relationship between the elastic constants, see equation (A-29), and one plane of symmetry parallel to the bounding surface. Maki [18] extends this solution using a general set of orthotropic elastic constants (same symmetry however). Lekhnitskii [12] and Green and Zerna [19] do the general problem (arbitrary orientation of axis of symmetry) for a vertical load. In this appendix the general problem is outlined for an arbitrary load direction. The singular solution (b-2) is given by Conway [20] for the same relationship between the elastic constants and the condition of plane stress. In this appendix the general case is also worked out.

26 The singular solution (a-3) is given by Conway [20]. Due to its complexity and limited use it is not presented. Its major use is for such problems such as inclusions or cavities in a half plane. The singular solution (b-3) is given by Conway [20] for the special case of one axis of symmetry parallel to the bounding surface and the special relationship between the elastic constants i.e., equation (A-29). The major use of this solution is for problems similar to those mentioned for (a-3). Again, because of its complexity, it is not presented. c. DERIVATION OF THE SINGULAR SOLUTION (b-l) The derivation of the stress field due to a line load on an orthotropic half space (see Fig. (A-l, b-l) follows the approach outlined by Lekhnitskii [12]. Two analytic functions of a complex variable zi = x + i. y are introduced; l1(z ), $2(z2), having derivatives*, 1 22 d1 (z 1 0 12oo Ny x ) dz1 1 1 2iri(1 12 - 0 -1 1 + N = - f - - ^ -^ c \ -^ (A-l) dz 2 2 i(p 4 2 1 2 where N (e), N (S) are the x and y traction components and p.i' ii are the roots (complex conjugates) of the fourth-order equations for either plane strain or stress given in (b-l) above. The derivation of these equations as well as their relationship to the stress components, equations (A-2) below, is discussed in detail in [12]. *These &. 's assume the z axis is an axis of symmetry. 1

27 The stress components in terms of the. 's are, 2 2Re xxa = 2Re J l(Zl)+ 2 2( 2 a = 2Re f|(zl) + 2( z2 (A-2) aY = -2Re L l(zl) + 2 z2(Z2 xy 2" 2_J^ To obtain the singular solution we need the stresses due to a concentrated line load T = T i + T j. This can be done as the limit of equation x y (A-l) i.e., we define, 0 -00 < < -E 0 -o < < -e T T N (- -E< E -e< N( = - < (A-3) L_0 o C < 0 << <o00 Introducing equations (A-3) into (A-i) with T t* ds and T = t* ds and x x y y taking the limit as e + O, we obtain, (2 t* + t*....... _y X ds 1 ) - - 2ti(I - 2) ds 1l t* + t* 1 y x (A-4) (z) = + ds Xz 2 2 2ir1i(,1 -12)z2 Substituting equations (A-4) into equations (A-2) we obtain the stress field due to the concentrated load t* ds. We can now introduce our previous notation, namely, Gi t* where again G i is the stress due to a unit line load in the q direction. Therefore, i q direction. Therefore,

28 Gx; - 2Re ~ (1 2 t - -) + - _))} xx; q 2e1(1(l - 2) 1 2 Y \1 z2 ( z G ca 2Re + th (A-t i /( sine + cos 21 ) G t* = 2Re 2 — t --- ------- xy;q q 4) y cos91 + sin G t* = O ' Introducing a coordinate transf;qrmation the components in polar coordinates wcan ber ithanla shown ito be, A b 2 Introucinr r A-we t rnsfo eratei the o nentK ipl cosed +ina1ts rI(r s + cose ) + 2(1 ssine1 + cose)2 -.......1.) -.c..+t* -1 - sin- (A-6) -2 cos1 + sine / x1 cose + sine ( 2 sine + cose) ~-2 coseG + sineG1 G t* = 0 ee;q q G t* = 0 re;q q where e1 is the angle as shown in Fig. A-i, b-1. In equations (A-6) we now separate the real and imaginary parts (only the 1i's are complex). It is convenient to define, k J= y +i (A-7) After separating, the following quantities appear and are defined for simpli city,

29 1 )((sYi lsine )(-cose ) -cS 1inel)cosO)+sin -2sinecos sin+ Re1 (- cose.++sinel) +2cos 2 + 5lcos12 2 2 26 sine (y sinE) +cose )(.-y cose +sin9e )+ cose) ((y sine +cosO) 8 sin 2) 1 1 1 1 _ 1 1 1, 1 1 -. 1- — s h e261sinel( sinel+C~OSl)( lc~sSl+sinel)+'lc~Sel((W 11inel+c~S el) 2 2 in ) Im1 = --—.....2 2 2'.. (-Y csel +sin9l) +62cos 1 222 2 2 262sin1(7 2sin1+cos91 )(7 Ycose1+sine1) cose ((Y2sin91+cos91)2 62sin 90) Introducing equation (A-8) into (A-6) G tsi can now be written as rr;q g follows, rr;q q R2((Y2 - y) + (6 - ) ) s2 1 2 Re 2 2 2 2 7 Res - Im +c ) Im2) - (7 - s1)(62 Re+c 61 Re2 + y Im - ly Imn ) + stxRe - Re2)(82 1) (Im -(2 Y1iI} 3(A-9) It is important to point out that y1, y2' 61' and 62 depend on the relative orientation of the axis of material symmetry with respect to the boundary of the half plane (see Fig. A-, b-l). d. DERIVATION OF ITHE SINGULAR SOLUTION (b-2) We now derive the expressions for the components of the stress tensor in an infinite orthotropic plane due to a line load applied at an internal point.

30 The strain-stress relations when the axis of symmetry coincides with the coordinate axis are [12], xx xx 11 xx 12 yy E = a + a (A-10) yy 12 xx 22 yy (.E _ f xy 66 xy Introducing the Airy stress function 0, the equilibrium equations are automatically satisfied if, 2 xx 2 2 a = 2 (A-ll) 2x ay - - - xy 8x`y The only nonzero compatibility equation in plane strain or generalized plane stress is, 2 2 2 xx + 2 - = Xy (A-12) by2 2 bx y Substituting (A-ll) into (A-10) and then into (A-12), we arrive at the following a fourth-order partial differential equation in ~, 4 4 4 -+ (2P + P 2 2 22 4 = (A-13) 11 4 (2 12 + 66) x2 + 22 (-) cy ax:y 2s axy

or in alternate form -t-2 - + a )a '2 ~+a a (A-14) where al and a2 are the roots of the following second-order algebraic equation* 2 12 +66 p l -2 -2 66a +.l = 0 (A-15) 22 22 Note that the roots aC, a2 can be either real or complex conjugates. Introducing the following change of variable proposed by Maki [18], 1I = ^T y (A-16) equation (A-14) becomes, 2 2~ ( 2> a~ 2 + a+ k2 a t + a2+ =kO (A-17) where, 2 1 k - 52 Any stress boundary value problem in plane orthotropic elasticity reduces then to the solution of equation (A-17) in the field, subject to an appropriate set of boundary conditions. Consider the problem shown in Fig. A-l, b-2 where the line load T is acting at the origin in the positive direction of the x axes. Consider also the infinite plane divided into two semi planes along the x axis with a *This is seen by putting equation (A-14) into the same form as equation (A-13) and comparing coefficient.

32 distribution of shearing stress such that the total force on each plane is T /2. Clearly the limit of this, as the stressed area goes to zero is the solution of the original problem provided the proper symmetry conditions are met, i.e., no vertical deflection v along the y = 0 surface. The solution is derived for the horizontal force and results then stated for the load in an arbitrary direction. We introduce the following stress function given by Conway [17], o = J 1 - sa x Ae +Be sinx =0 2 L a (A-18) where A and B are two unknown constants at this point. In terms of this stress function, the stresses are, 2 1 K -ca' B - -T xa 1f IAe J+ e + e sina x da xx 2 a k2 2 0 2 La~Tk - c00 --- -= cra = 2- - Ae + Be ki sina x da (A-19) yy 02 a_ xy 2 Fll _ _ 1,00 A - B - == - -T = - - e + e k cos x d yax f O k jcosaxda CY ^? J We consider a rectangular distribution of shearing stress on the boundary of length 2c and magnitude T /2c (we want the limit as c - 0). This can be represented as follows,

33 O; -oo < < x < T — x f coscx x da 2T- 0 Equating equation (A-20) with a of equation (A-19) on y = 0 ( r = O) we obtain, xy A + + - Cosa x dxa = 0 (A-21) 0 k 2 2_J which implies that, B -T A + a0 0o (A-22) k 2 2r2 Equation (A-22) represents one relationship between A and B. We now satisfy the condition of v = 0 along r = 0. We first find the expression for v in terms of 4, i.e., v(x, y) = / E dy YY Making use of equation (A-10) and (A-19), we obtain for v(x, q), v(x, ) = 22 o Ae + + k Be k sina x da L 12 00 1kB + 12 1 Ae'n + - e k sinc x da dy (A-23) for y = 0 (: = 0), ~ = 2 fo 12 B v(x, 0) = A + B) - 12(A + da = aa(~~~~~~~~~~~~~~~__

Which again implies, 1(2 B 22(A+ k ^ B) - (A +) = o 2 Solving for A and B from equations (A-22) and (A-24) we obtain, T 2 2 T x k 22 2 12 A 2 B = _}2 2 22 - 2^ C^a " kR2)~ (A-24) (A-25) With A and B known the stress = H T xx xx;x x a = H T yy yy;x x components* are given by equations (A-19), i.e., T xa k2 T x e2 k2 822 1 \2 2 22 - 2 1 )kl x1 2 22 2 2 r 22 2 2 2 522(1 - k2) I k2 x + y k2 k2 + J x 2 2 22 12 ( 2 822 - 12)k 2 2 2 2 2 2 2iTr P2(1-k ) c x + y a2 k x + y x PP O - = H T = 2 2n/ c 22(1-k) 2 x + Y k x + y 2 22 - axy (A-26) Following a similar procedure, expressions for the components of the stress tensor can be obtained for a line load T coincident with the positive direction y of the y axes. For a line load in a general direction, T, we can form Hj; T. iJ;q q Following the results of equation (4)-(7) we define T = F dV where F is a body force field of equation (5) and the subsequent limit, we see that *From previously defined notation, it is clear that axx here is Hxx;x Tx where Hxx;x is the normal component of stress in the x direction due to a unit load in the x direction. I

35 H T = H f* ds when ds is an element of length in the infinite plane. ij;q q ij;q q Therefore, C f* r + C f* r H f* 1 1 x x 2 y y xx;q q 2i(a2-a1) L 2 y 2 1 2 x y H f=:2 yyxx;q q 2 1) a2 r +r 2 x y C f* r C f* r 3 x x 4 y y -2 -2 ac r + r l x y _ C f* r + C f* r _7 x x 8 y y -2 -2 ca r + r 1 x y _ C f* r - C f* r 3 x y 8 y xi -2 -2 a1 r + r (A-27) where r = x - x' r y - y ', valid for the load at an arbitrary point x y x, y. And, C1 2 1 P22 5 2 L P22' C 1/2 a 7 12= C3/2 12 45 = "1 "2 - i 22 8 1 2 22 ^ 04 = 1/2L| a _2C — j - c 12 equations (A-27), (A-28) are the ones used in Section III, i.e., equations (19), (20). From a computer programming point of view, it is important to distinguish the two cases, a1, a2 real, and al, a2 complex conjugates. Equations (A-27), (A-28) work in the real case. An alternate form of equations (A-27), (A-28) is derived in Appendix D when aO, a2 are complex conjugates.

36 It is interesting to point out that equations (A-27) reduce to ones given by Conway [17] if a specific relationship between the constants is introduced, E E xy E x+E (1 + 2 vyx) See Appendix C for the relationship between the..'s and the elastic constants. iJ

37 APPENDIX B DERIVATION OF THE ALTERNATE SET OF INTEGRAL EQUATIONS (Using as the Singular Solution the Line Load Acting On the Surface of an Orthotropic Half Plane) We have obtained in Appendix A expressions for G;, Gr G;, the rr;q r;q 9;q' cynlindrical components of the stress tensor due to a unit line load in the q direction acting on the bounding surface of an orthotropic half plane. Following the approach given by equation (17) we now derive the set of integral equations for the specific choice of G.. function (above) and a real geometry ij;q defined by the equation f(x, y) = 0, see Fig. B-1. This consists of three piecewise continuous segments for the geometry of the specific problem of Section III. The major difficulty here is defining the function ij. lJ;q The cartesian components of the stress for a line load on an orthotropic half plane are given by, 2 G t* x G t* xx;q q 2 2 rr;q q x +y G t xy G t (B-1) xy;q q 2+ 2 rr;q q x +y 2 G t* = G t* yy;q q 2 2 rr;q q where G is the radial stress due to a unit line load in the q direction and rr;q is given by equation (A-9). The definition of G.. requires that the region Laq B* move around B such that the two are tangent at at least one point. In equation (A-9) are several expressions which change as B* moves around B, i.e.,

38 c's, 5's, P's, Re's, Im's. Expressions for these changes will be derived in terms of their values in a reference direction (x axis), which is chosen to be one of the axis of material symmetry, fixed in B. Referring to Fig. B-l, as the fictitious boundary B* is moved over the real boundary B (such that the two are tangent), the axis of material symmetry of the fictitious region D* changeso This creates the fundamental difficulty in the definition of ij; We first derive expressions for the changes in pi' i. [ 1' 2' i1' and [2 are the roots* of either of the fourth-order equations of Appendix A, section b-l and enter the expression r. The coefficients of these equations rr;q depend on the relative orientation of the axis of symmetry. Hence, it is necessary to solve these equations for each point of coincidence of B*, B. It can be shown, however, that there exists [12] a relationship between the t's for two orientations an angle r apart (one orientation is the x axis) (see Fig. B-1). r p1 cost + sint 1l r cost * 1 sinf r g2 cost + sint 2:r --- -(B-2) cost- -2 sint r r where ~1' p2 are the roots for the reference position. With respect to the fixed set of coordinates (x, y) in B, the unit normal to B at P is given in terms of t, n = sint n cos (B-3) x y *Note that the roots 1i' Y2 l' 21' are complex,

39 Substituting equation (B-3) into equation (B-2) we obtain the following expressions for l', 2', r ~ln + n 1 y x HIl- r n - l n y 1 x r 2 n + n p ' = 2 Iy x (B-4) =2 r ny-11 nx Separating the real and imaginary parts we end up with four expressions for 1, 72, and &2 in terms of 71, 72 &r and 62 and the direction cosines n, n at each point on the boundary, where the y's and b's are defined as follows, r r r k = Yk + i; k = 1, 2 These expressions are r 2 2 r2 r2 [n - n ] + n n [1-(y1) - ( ) ] = y x xy 1 1 '1 [lyn Y- ni + (6 n) y x x... 1 [n r n ]2+ n ) y l x lx [n - n2] + n n l-(7r)2- (6J)2] 2 y x x y 2 2 (B5) 2 r r 2 r 2 I[n - y n] +(6 n) [ny 2 x 2 x 1r S2 2 r- r 2 6 = 2 [n - 7Yr nx 2+ (2 n) In the expressions for Rel, Re2, Iml Im2 the angle 01 is as shown in Fig. B-l. It represents the angle between the normal to the boundary B and the radius vector to the field point P. This changes as P1 moves and must be

40 defined in terms of J, i.e., nx, n. Defining r = x - x', r y - y', the expressions for cos 01 and sin e appearing in equations (A-8) are, n r + n r x x y y cose =iosl -2 -2 1/2 [r + r x y n r -n r sine - y x x y (B-6) 1 — 2 -2 1/2 [r + r x y With the changes in the angle 01 and the 7's, 6's defined previously the new expressions for Rel, Re2, Iml, Im are found by substituting directly into equation (A-8). In equation (A-9) t* and t* represent the normal and tangential fictitious y x traction components as shown in Fig. A-1, b-1. As B* moves over B the normal and tangential components of t* are no longer coincident with t* and t*. We y x can, however, write, t* = t* n + t n; on B n x x y y t* = t* n - t n; on B (B-7) s x y y x t*, t* now replace t*, t* in equation (A-9). n s X y The definition of $ t* is now complete if the equations defined above rr;q q are substituted into equation (A-9). The resulting integral equations therefore, are,

41 r n +r n B x y y - t* f X _ _ + __ r_ X -2 -22 (62 x ) 2 x [2 2 )2] ftn ((6 -61)(7' Re - y Re2 62 Im + 6 Im) tn 2 1 2 1 1 2 2 I2 1 i2 - (2 1)( 2 Re - R + 7 I y Im)) + t* ((Rel-Re )(62-b) (I -Im2)(7-)) ds = t r n + r n t* + x ay y -_ _ _ _ Y — 2 -2 2 r 2.. Y B (r2+ rP)2 [(y _7 )2 + (6 6)2] { n((2-1)( Rel Y1 Re2 62 I + 1 Im2) - ( )( Re1 - 6 Re + )) 2 1 2 1 1 Re2 + y I y ) + t* ((Rel - Re )( 62-6 - ds (B-8) s 2 NY - I B- s where t*, t*, 's, 6's, Re's, Im's, r, r, n, n are given above. n s x y x y Equations (B-8) represent the alternate way to the one given in Section III for solving an orthotropic two-dimensional problem. As can be seen, if equation (B-8) is compared with (23) the amount of complexity involved for the orthotropic case is greater. A computer program for the solution of the problem presented in Section III using this approach has been written and results obtained, however, for brevity it is not included. Note that one more step is involved in casting equations (B-8) in the final form for the problem of Section III, i.e., direction cosines must be defined for the specific geometry.

42 In case of APPENDIX C ELASTIC CONSTANTS FOR PLANE STRAIN orthotropic elasticity the stress-strain relations are given for the the axis of symmetry coincident with the axes of coordinates by [12]; xx 11 xx 12 yy 13 zz E = a a + a a + a yy 12 xx 22 yy 25 zz zz = a + + aa azz (C-l) zz 13 xx 23 yy 33 zz E = a ar xz 44 xz E = a ar yz 55 yz E = a a xy 66 xy In plane x xx strain we have the conditions E = = E = O, hence zz xz Yz = c + 12 a 11 xx 12 yy e = 6 a + aC yy 12 xx 22 yy E = - 66 xy xy 66 xy (C-2) We define P. P66 P,. as 1J... = - a.. l0 follows a. a 3; i, j = 1 2 33 (C-3) = a66

The relationship between the a. 's and the Young's moduli E, E, E 13 xX yy ZZ Poisson's ratios v,v, v, v, and v and the shear moduli G, xy yx xz yz zx zy xy G, G are xz yz v v 1 yx xy 1 a a = - a1= 11 - E E; 44 G xx xx yy xz v v 1 xz zx 1 22 = -' = - 55 G= (-4) yy zz xx yz v v 1 a zy yz 1 a33 2 E E E; 66 G zz yy zz xy hence substituting (C-4) into (C-3) the expressions for the..'s in terms of 13 the elastic constants will be 1 1 E - (1- v zx XX 1 22 E (=1 - v v ) (C-5),22 E yz zy yy 12 E yx zx yz E y xz zy xx yy 66 = xy These are the expressions given in equation (26). Note that if we were solving a plane stress problem, we have to substitute in place of the ij's in (A-28) the corresponding constants for the plane stress case, i.e, a Equations (C-4) are stillvalid case, i.e., a.. Equations (C-4) are still valid. 13

APPENDIX D ALTERNATIVE WAY OF CASTING EQUATIONS (A-27) FOR a 1 a2 COMPLEX CONJUGATES The two roots of equation (A-15) can be either real or complex conjugates. For both cases the final expressions given by equation (A-27) are real. For the sake of programming we include here an alternative way of presenting equations (A-27) when a1 and a2 are complex conjugates. We first define r i a = r *i ~^ (D-l) 2 2 2 r i a = a -ia 1 2 2 r i where a2 and a2 are the real and imaginary part of a2. The modulus of the complex number a2 is = r2 i21/2 1ax1 = [(ar)2 + (a)2]12 (D-2) and the argument arg[a ] = tan - (p5) 2 r hence we can write = 1a, ei[arg(a2)] (D-4) Substituting (D-4) into (A-28) and defining

C2 = cos 2 [arg (a2) ] S2 = sin [arg (cx2)] (D-5) C3 3 = COS 2 [arg (c2)] S3 = sin 3 [arg (ox2)] We obtain the following expressions for 1 = 2l~2(2 1 222 the constants C1 - C8; 1- i 2 - i S2 a 2 jS I - p 1 c2 2 la1 C2 I p12 2%^) + i S2( + oIx I 12\ a ) — I 12 22/ c3 = C1 C4 = C2 53 (D-6) 21.12 C2 - 12 C3 22 12 2 11 3 J2 + i IIJ 2 + i (S3 I + s: 12 I s2 -, 22 - -o 22 - c6 C7 = 5 8 = c6 where C1 C, C2, C C6 are complex conjugates of C1, C2, C, C6. Defining further;

46 r -2 -2 A 2 rr + r 2 x y i -2 B = 2 rx (D-7) r i C. = C + iC j = 1, 2, 5, 6 J J J r i where C. and C. are the real and imaginary parts of C. defined in equation (D-6) above. Equation (A-27) will now read; I r i r - 1 (C A- C' B) f*r + (C A- C B) f* r x x 2 2 y y H fy 2 xx;q q L2 2 y;q^ 2c2r A + B i A r B) f* ( H -'*' f* = --- iv -5 ----x x — 6 6 y Bf] Equations(D8)ands (-) are equivalent. The former is for c(C A- C1 B) f*r ' (C A ra B) f* 1 1 16 x yx xy;q q i 2 2 1 2nMr L A + B Equations (D-8) and equations (A-27) are equivalent. The former is for Ql, a complex conjugates and the latter for al, a2 real.

47 APPENDIX E COMPUTER PROGRAM IX this appendix we include the computer program for the numerical solution of equations (28). The program computes the values of the constants C1,...,C8, equations (20) and then for any two locations (I,J) on the boundary the coefficients of equations (28). These now represent a system of algebraic equations whose unknowns are P*, P* (i = 1,...,N). This system xi yi is solved by iteration using equations (30). Once the final P* and P* are x y known, the program is prepared to compute the stresses at any point within the field through the use of equation (31). As input information to the program the following items are necessary: 1. Coordinates of the center point of each of the subdivisions of the boundary referred to a set of coordinate axes coincident with the axes of symmetry of the material, i.e., X(I), Y(I). The origin can be set anywhere. In the case solved in Section III the origin was taken, for convenience, at the intersection of the boundaries BD, AC. See Fig. 7. 2. Components of the unit normal at each X(I), Y(I); ANX(I), ANY(I). 3. Total number of subdivisions on the boundary, (N). 4. Numberof points where the stresses are to be found in the field, (M). 5. Traction on the boundary, represented by the cartesian components P' P i; (TX(I), TY(I)) at all points X(I), Y(I), I = 1,...N. xi, '''

6. Degree of accuracy required in the iteration, (PC). 7. Value of the convergence parameter, (AL). 8. Material properties E, (EXX), Eyy (EYY), G, (GXY), v, (PXY), Vyv (PYX), v, (PXZ), (PZX), v yz (PYZ), v (PZY). As output the program will yield the following information. 1. Components of the fictitious traction represented by the concentrated fictitious load at the center of each one of the subdivisions, TXS(I), TYS(I), I = 1,...N. 2. Components of the stress tensor a,, a, at each of the points xx yy Nxy desired, (SGMAXX(I), SGMAY(I), SOMAXY(I), I = 1,M). 5. Principal stresses ol2 and their directions at the same prescribed locations (SGMAl(I), SGMA2(I), Zl(I), Z2(I), I = 1,M). The computer programs follow.

FORTRAN IV G COMPILER - MAIN PROGRAM DIMENSION* TX(-),TY(-),LOCT(-),S(-),X(-),Y(-),ANX(-), 1 ANY(-),F1(-,-),F2(-,-),F3(-,-),F4(-, -),T(-), 2 TYS(-),TXSN(-),TYSN(-),SUMX(-),SUMY(-),XN(-),YN(-), 3 SGMAXX(-),SGMAYY(-),SGMAXY(-),SGMA1(-), SGMA2( -), 4 ZETA1(-),ZETA2(-) 201 FORMAT ('CICLE =',I4/'0' /LOCT', 1X,'TXS',30X,'TYS'/'Ot/ 1 (I4,8X, E20.8,16X, E20.8)) 300 FORMAT ( 'LOCT', 10X, TSGMAXX',30X,'SGMAYYT,30X,' SGMAXY'/ '/ 1 (I4,8X, E20.8, 16X, E20. 8, 16X, E20. 8)) 400 FORMAT (' LOCT',8X, SGMA1',20X, SGMA2 ', 20X, 'ZETA1T, 20X, 1 tZETA2'/' O/( I4,8X, E20. 8,6X, E20. 8, 6X, E20.8,6X, E20. 8)) NAMELIST/DATA/PC,AL, IN, M, XN, YN, ID, X, Y,ANX,ANY, S, TX, TY, EXX, EYY, GXY, PXY, PYX, PZX, PXZ, PYZ, PZY NAMELIST/OUT/A1,A2 NAMELIST/CONST/C1, C2, C3, C4, C5,C6, cC7, C8 NAMELIST/OUTI/A2R,A2I NAMELIST/COSTI/ C1R,C1I,C2R,C2I,C5R,C5I, C6R,C6 1 READ (5,DATA) WRITE (6,DATA) PI=3.1416 Bll=( 1. O-PXZ*PZX)/EXX B12=- ( PXY+PXZ*PZY) /EYY B22=( 1. 0-PYZ*PZY) /EYY B66=1. O/GXY P1=(2. O*B12+B66) /B22 P2=B11/B22 800 A2R=P1/2.0 A2I=SQRT(4 O*P2-P1**2)/2. 0 WRITE (6,OUTI) TM=SQRT(A2R**2+A2I**2) AR=ATAN(A2I/A2R) A22=B12/B22 A11=B12/B11 SM=SQRT(TM) CR=COS(AR/2. 0) SR=SIN(AR/2. 0) CR3=COS(3. O*AR/2. 0) SR3=SIN(3. O*AR/2. 0) C1R= SM*CR*(TM-A22) C1I=-SM*SR*( TM+A22) C2R= SM*CR*( 1. O-All*TM) C2I=SM*SR*( 1. O+All*TM) C5R=( SM*3 )*( TM*CR-A22*CR )

50 C5I=( SM*35)*( TM*SR-A22*SR ) C6R=( SM**3 )*( CR3-All*TM*CR) C6I=( SM**3)*( SR3-All*TM* SR) WRITE(6,COSTI) R=1. O/(2. O*PI*A2I) DO 81 I=1,N DO 71 J=1,N RX=X(I)-X(J) RY=Y(I) -Y(J) A=A2R* (RX**2) +RY**2 BA=A2I*(RX**2) AL1=RX*ANY(I) AL2=RY*ANX(I) PN=RX*ANX( I) +RY*ANY ( I) PM=C1I*A-C1R* BA PQ=AL2*(C2I*A-C2R*BA)-AL1*(C6I*A-C6R*BA) PR=AL2*( C1I*A-CR*BA) -AL1*(C5I*A-C5R*BA) PS=C6I*A-C6R* BA PG=A**2+BA**2 IF (I. NE. J) GO TO 61 F1(I,J)=O.0 F2(I,J)=O. O F3(I, J)=0. 0 F4( I, J)=0. 0 GO TO 71 61 F1(I, J)= R*PN*PM/PG*S(I) F2(I, J)= R*PQ/PG*S(I) F3(I, J)= R*PR/PG*S(I) F4( I, J)=-R*PN*PS/PG*S(I) 71 CONTINUE 81 CONTINUE 125 DO 9 I=1,N TXS(I)=TX(I) TYS(I)=TY(I) 9 CONTINUE NCICLE=1 GO TO 12 10 IF (NCICLE.GT.20) GO TO 22 I=1 NCICLE=NCICLE+1 DO 11 I=1,N TXS(I)=TXSN(I) TYS(I)=TYSN(I) 11 CONTINUE 12 DO 13 I=1,N SUMX(I)=O.0 SUMY(I)=O. 0

14 15 17 22 18 900 119 120 CONTINUE DO 15 I=1, N DO 14 J=1, N SX=TXS( J)*F1( I, J) +TYS( J)*F2( I, J) SY=TXS( J)*F3 (I, J) +TYS( J)*F4 ( I, J) SUMX(I)=SUMX(I)+SX SUMY(I)=SUMY(I)+SY CONTINUE CONTINUE DO 16 I=1,N TXSN(I) =AL*2. O*TX(I)+(1. O-AL)*TXS(I)-2. O*AL*SUMX(I) TYSN(I)=AL*2. O*TY(I)+( 1. O-AL)*TYS(I) -2. O*AL*SUMY(I) CONTINUE WRITE(6,201) NCICLE,(LOCT(I),TXSN(I),TYSN(I), I=1,N) DO 17 I=1,N EPSX=PC*ABS(TXS(I)) EPSY=PC*ABS(TYS(I)) DTXS=ABS( TXSN( I) -TXS(I)) DTYS=ABS( TYSN(I)-TYS(I)) IF (DTXS.GT.EPSX) GO TO 10 IF (DTYS.GT.EPSY) GO TO 10 CONTINUE WRITE(6,201) NCICLE,(LOCT(I), TXSN(I),TYSN(I), I=1, N) DO 18 I=1,M SGMAXX(I)=O. O SGMAYY(I)=O. 0 SGMAXY(I)=O. CONTINUE DO 120 I=1,M DO 119 J=1,N RXN=XN(I)-X(J) RYN=YN(I)-Y(J) AN=A2R* ( RXN**2) +RYN**2 BN=A2I*( RXN**2) C 1N=C I*AN-C1R*BN C2N=C2I*AN-C2R*BN C5N=C5 I*AN-C5R* BN C6N=C6 I*AN-C6R*BN PGN=AN**2+BN**2 SNXX=R* (C1N*TXSN( J)*RXN+C2N*TYSN( J)*RYN)/PGN SNYY=-R*(C5N*TXSN( J)*RXN+C6N*TYSN( J)*RYN)/PGN SNXY=R* (C 1N*TXSN( J) *RYN-C6N*TYSN( J) *RXN) /PGN SGMAXX( I)=SGMAXX(I) +SNXX SGMAYY(I) =SGMAYY( I) +SNYY SGMAXY( I)=SGMAXY( I)+SNXY CONTINUE CONTINUE

52 126 WRITE (6,300) (LOCT(I),SGMAXX(I),SGMAYY(I),SGMAXY(I), I=1,M) DO 60 I=1,M SGMA1( )=( SGMAXX( I) +GMAYY( I) )/2. 0+ 1 SQRT( ((SGMAXX(I)-SGMAYY(I) )/2. 0)**2+SGMAXY(I)**2) SGMA2(I)=( SGMAXX(I)+SGMAYY(I) )/2. 0 -1 SQRT(((SGMAXX(I)-SGMAYY(I))/2. )**2+SGMAXY(I)**2) IF (ABS(SGMAXX(I)-SGMYY(I)). LT. 0. 00001) GO T 62 ZETA1(I)=ABS( O. 5*ATN( 2. O*SGMAXY()/(SGMAYY(I) -SGMAXX(I) ) ) )* 1 (180.O/PI) GO TO 63 62 ZETA1(I)=45.0 63 ZETA2(I)=90. O+ZETA1(I) 60 CONTINUE WRITE (6,400) (LOCT(I),SGMA1(I),SGMA2(I),ZETA1(I), 1 ZETA2(I), I=1,M) GO TO 1 END *The dash in the arguments of the dimension statement represent the number of desired points on either the boundary or in the field.

Table 1. Components of the Fictitious Loading, see Figure 7. Note that points 1-20 on AC. are on BD, 21-24 on AB, 25-44 Boundary P*. P* Boundary P* P*. xi yi xi YLyi Location Location 1 -0.005320 -0.003425 23 0.299086 -0.194660 2 -0.005083 -0.004384 24 0.091765 o.031661 3 -0.004926 -0.004999 25 0.043434.oo005890 4 -0.004819 -0.005521 26 0.023564 -0.002726 5 — 0.004740 -0.006020 27 0.019134 -0.002559 6 -0.004677 -. 006327 28 o. 016662 -0.001143 7 -0.004620 -0.007063 29 0.014785.00013 4 8 -0.004563 -0.007644 30 o0.013256 0.00012 9 -0.004497 -0.008288 31 0.011985 0.001809 10 -0.004412 -0.009015 32. 010917 0.002329 11 -0.004294 -0.009852 33 0.010006 0.002716 12 -0.004121 -0. 010829 34. 009220. 003012 13 -0.003862 -0.011993 0.008552 0.00245 14 -0.003462 -0.013400 36 0.007921 0.003435 15 -0.002837 -0.015134 37 o.007569 0.003594 16 -0.001849 -0.017309 38 0.006864 0.003735 17 -0.000289 -0.020104 39.oo006390 0.003866 18 0.002073 -0.023863 40. 0059o 4 0.003996 19 0.004909 -0.030013 41 0.005480 0.004134 20 0.003860 -0.052446 42 0.004999gg o.004294 21 0.051697 -0.182251 43 0.004427 0.004502 22 0.362376 -0.398540 44 0.003521 0.004797

54 Table 2. Principal Stresses, Direction (principle direction is the angle between y axis and direction of maximum principle stress). Note that AC=BD=5L and mesh size = L/4. Field Max. Prin. Min. Prin. Prin. Location Stress Stress Direction 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 0. 001820 o. 003678 0.004902 o. oo4895 0. 004580. 007178 O. 007077 0. 005026 -0. 041126 -0. 000722 0.008279 o. 004030 -o. 010163 o. 005674 0.001998 -0. 191869 -0. 208151 -0. 197443 -0. 168405 -0. 290822 -0. 309108 -0. 248241 -0. 185699 -0. 613534 -0. 479524 -o. 270639 -0. 179982 -0. 460534 -0. 222955 -0.153953 10.84~ 22 38~ 36.58~ 43. 09~ 13.96~ 30.66~ 42. 29 ~ 32.532 23.91~ 40. 170 24.49~ 18.36~ 13. 64 8.66~ 7.99~

55 Table 3. Effect of Infinite Boundary, Mesh Size on the Stresses at Location 10. Max. Prin. Stress Min. Prin. Stress Prin. Direction Boundary Length AC,BD (Mesh Length = L/4) 3L 5L 8L Mesh Length (AC=BD=3L) L/4 L/8 L/12 -o. 00004774 -0. 0007224 -0. 0007728 -0. 00004724 -o. 0006542 -0. 0012311 -0. 479372 -0.479584 -0. 479749 -0.479372 -0. 506492 -0.515493 40. 16 40.17~ 40. 19~ 40. 16~ 4o. 00~ 39.97~

56 Acknowledgment The authors would like to express their sincere appreciation to the National Science Foundation (Grant No. GK17101) for the support of this work.

57 BIBLIOGRAPHY 1. Teodorescu, P. P. "One Hundred Years of Investigations in the Plane Problem of the Theory of Elasticity," Applied Mechanics Surveys, Ed. by Abramson, H. N., et al., Spartan Books, 1966, pp. 245-262. 2. Muskhelishvili, N. I. Some Basic Problems of the Mathematical Theory of Elasticity, P. Noordhoff lt., 1963. 3. Massonnet, Ch. "Numerical Use of Integral Procedure," Chapter 10 of Stress Analysis-Recent Developments in Numerical and Experimental Methods, Ed. by Zienkiewicz, O.C., and Holister, G.S., J. Wiley and Sons, 1965. 4. Massonnet, Ch. "Resolution graphomecanique de problemes generaux de l'elasticite' plane," Bull. Centre Etudes Rech. Essais Sci. Construc. Genie Civil, Liege, IV, 1949. 5. Massonnet, Ch., et al. "Application des machines a calculater electroniques a la solution du probleme aux tensions de l'elasticite' plane," 6th Congress of the I.A.B.S.E., Stockholm, 1960, p. 95. 6. Wu, T. and Chiu, Y. P. "On the Contact Problem of Layered Elastic Bodies," Quart. Appl. Math., Vol. 25, 1967, pp. 233-242. 7. Boley, B. A. "A Method of Heat Conduction Analysis of Melting and Solidification Problems," Journal of Math. and Phys., Vol. 15, 1961, pp. 300-313.

58 8. Sikarskie, D. L. and Boley, B. A. "The Solution of a Class of TwoDimensional Melting and Solidification Problems," Int. J. Solids Structures, Vol. 1, 1965, pp. 207-234. 9. Tricomi, F. G. Integral Equations, Interscience Publishers, New York, 1957. 10. Paul, B. and Sikarskie, D. L. "Preliminary Theory of Static Penetration by a Rigid Wedge into a Brittle Material," Trans. AIME, Vol. 232, 1965, pp. 572-383. 11. Benjumea, R. and Sikarskie, D. L. "A Note on the Penetration of a Rigid Wedge into a Nonisotropic Brittle Material," Int. J. Rock Mech. Min. Sci., Vol. 6, 1969, pp. 3453-352. 12. Lekhnitskii, S. G. "Theory of Elasticity of an Anisotropic Elastic Body " Holdenday, Inc., San Francisco, 1963. 13. Duvall, W. I. "Effects of Anisotropy on the Determination of Dynamic Elastic Constants of Rock," Soc. Petr. Eng., Trans. AIME, Vol. 236, 1965, pp. 309-316. 14. Hildebrand, F. B. Introduction to Numerical Analysis, McGraw-Hill, Inc., 1956. 15. Timoshenko, S. Theory of Elasticity, McGraw-Hill, 1951. 16. Sneddon, I. N. Fourier Transforms, McGraw-Hill, 1951.

59 17. Conway, H. D. "Some Problems of the Orthotropic Plane Stress," Journal of Applied Mechanics, Vol. 75, 1953, PP. 72-76. 18. Maki, A. C. "Stress Distribution in an Orthotropic Half Plane Subjected to a Concentrated Load," U. S. Forest Products Lab. Report No. FPL54, February, 1966. 19. Green, A. E. and Zerna, W. Theoretical Elasticity, Oxford, 1965. 20. Conway, H. D. "The Stress Distribution Induced by Concentrated Loads Acting in Isotropic and Orthotropic Half Planes," Journal of Applied Mech., Vol. 75, 1953, PP. 82-86.

FIGURE CAPTIONS Figure 1. Two-Dimensional Solid Submitted to a Boundary Traction t. 2. Two-Dimensional Solid Submitted to a Body Force Field F. 3. Real Region D Embedded in a Fictitious Region D*. 4. Real, Fictitious Regions Tangent at P'', P''',.... 5. Example Problem (general). 6. Choice of Coordinate Axis. 7. Specific Geometry, Calculation Points for Numerical Example. A-1. Singular Solutions for Isotropic, Orthotropic Regions. B-l. G.. Formulation - Plane Boundary B* Moves Over All Points of B. la;q

I I LyU-,\ * Ttds ~P ~~ ~ d Fig. 1 (" (011

Tr FdV B Fig. 2

Iy ~r '~ ^ II ''! / n '~ I I I D \ D* B B* \% x 'I*I - I Fig. 3

_zt \ID t? *BTFd / - / s I I / _" f I / \ \% - q~^ *4

- e8 L ~ ~y Planeof, '. o?<ort hotropic symmetry CO Fig. 5 0\ k

-___ Plane of orthotropic symmetry Fig. 6 \. oh

y -- n } 13 14 15- x / 9 O 10 12 4 5 - - Plane of L.2!_1_1 i / orthotrox 2 3 4 symme C Fig. 7 Yic ty 0\

t-V 'STI (~-q) (Z-q) (I-q ) x —1 A.^- x- A /dx A (~-D) (Z-D) (I-o) X p -.4!rl 99

y %%Ob %**% ".. I - ^ l/. r ix,y/ =u / P/one of orthotropic symmetrv - / also reference plane ^%^~ '~s^ __ F~NO 00 Fictiious, orthotropic 0 — W.... am hofalf plane Fig. B-1 0\ O0

UNIVERSITY OF MICHIGAN 3 9015 02IIIII229 3172 3 9015 02229 3172