ABSTRACT A NUMERICAL METHOD FOR THE SOLUTION OF PROBLEMS IN THREE DIMENSIONAL ELASTICITY by Hotten A. Elleby A numerical method was developed based on expanded finite differences of the displacements of the nodal points for use in the solution of problems involving solid bodies. The method presented was derived using the equilibrium of the region around a nodal point of the finite difference grid as the basis of solution rather than the classical approach which uses the equilibrium of the nodal point as the basis for solution. To investigate this method, the solutions using this method were found for a concentrated line load acting on the surface of a cube, a concentrated load acting on the cube, a distributed load act- ing on the cube, an eccentric load acting on a beam column, and a beam column under constant applied moment. These solutions were also studied with respect to variations in Poisson's ratio. The solution of these problems was determined without the use of external fictitious nodal points to satisfy the boundary stress con- ditions, and also without the use of superposition of solutions when the problems contained singularities. ABSTRACT Z Hotten A. Elleby The solutions of these problems converged very rapidly with respect to decreasing grid spacing, except in the proximity of singularitie s . The equations became very unstable as Poisson's ratio approached one -half. This is the case of zero dilatancy of the material which means the average stress at a point is independent of the strains at that point. However, within the usual range of Poisson's ratio, the solutions appear to have an accuracy of at least three percent when not in close proximity to points of singularity and with a reasonable choice of grid size. A NUMERICAL METHOD FOR THE SOLUTION OF PROBLEMS IN THREE DIMENSIONAL ELASTICITY BY \' Hotten A‘. " Elleby A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Civil and Sanitary Engineering 1964 ACKNOWLEDGEMENTS The author wishes to express his sincere appreciation and thanks to his major professor, Dr. Charles E. Cutts, for his generous help and encouragement which made this thesis possible. The author would also like to thank Dr. William A. Bradley for his invaluable help and also Dr. Robert K. Wen and Dr. Edward A. Nordhaus for their guidance. Acknowledgement is made to the Computer Laboratory for their services in the use of the Control Data Corporation 3600 computer . ii TAB LE OF CONTENTS Page I. INTRODUCTION ............................. 1 II. THEORETICAL DEVELOPMENT OF THE METHOD OF EXPANDED DIFFERENCES ..... 4 III. METHOD OF SOLUTION .................... 31 IV. EXAMPLE PROBLEMS ...................... 37 Example Problem One 45 Example Problem Two and Three 63 Example Problem Four 76 Example Problem Five 88 V. SUMMARY AND CONCLUSIONS .............. 91 VI. APPENDIX .................................. 96 VII. BIBLIOGRAPHY ............................ 108 iii Table LIST OF TA B LES Page Vector components used for the determination of an octant . .............................. . ..... 15 Values of r used to determine the location of the centroid of a stress on an octant ........... . . . . 27 Values of oz for the two solutions with various grid spacings .......... . ...... . ............... 54 iv Figure LIST OF FIGURES Page Isometric view of basic model octants surrounding the reference nodal point, and also the rectangular parallelepipeds of which they are apart.. ..... .......... 10 Position of stresses as used in Equations 9 for the classical finite difference solution of the three dimensional equilibrium equations ........ . . . . . . 11 Position of the centroidal stresses on the first octant, l, of the model, and also the nodal points which determine the internal displacement func- tions withinthis octant 12 Isometric view of cube with line load for Example ProblemOne ........ 40 Isometric view of a cube showing concentrated load on the Z-Axis for Example Problem Two . . . . 41 Isometric view of a cube showing concentric distributed load over small portion of the top surface for Example Problem Three ..... . . . . . . . 42 Isometric view of a beam column showing the location of the eccentric load for Example Prob- lemFour...... ....... ........... 43 Figure 10 ll 12 13 14 15 16 17 Page Isometric view of a beam column showing the surface stress distribution for the applied moment for Example Problem Five . . . . . . ......... . . . . . 44 Outline of a thin plate with the concentrated loads and the associated boundary phi functions for the biharmonic solution for Example Problem One . . . . 50 Values of oz for the two solutions plotted with respect to the inverse of the mesh size ...... . . . . 55 Values of 0'; on the centroidal Z-Axis for various values of z plotted with respect to Poisson's ratio. 56 Values of (TX on the centroidal Z-Axis for various values of z plotted with respect to Poisson's ratio. 57 Values of oz on the centroidal Z-Axis. for z = 15 B/16 plotted with respect to Poisson's ratio. 58 O'x and cry on the centroidal Z-Axis for 'v = 0.3 . . . . 59 Lines of constant ax on the Y-Z plane with x = 0 forPoisson'sratioof0.3 . 60 Lines of constant cry on the Y-Z plane with x = 0 for Poisson's ratio of 0.3 ......... . ..... . ..... 61 Lines of constant oz on the Y-Z plane with x = O for Poisson's ratio of 0. 3 ......... . ...... . . . . . 62 vi Figure 18 19 20 21 22 23 24 25 26 27 Distributed load for Problem Two with a grid spacing of B/8 (a) and for Problem Three with a grid spacing of B/16 (b) ..................... Values of crz on the centroidal Z-Axis for various values of z plotted with respect to Poisson's ratio . Values of Oz on the centroidal Z-Axis for z = 15B/16 plotted with respect to Poisson's ratio. Values of 0x and cry on the centroidal Z-Axis for various values of crz plotted with respect to Poisson's ratio ...... . . . . . . ................... Values of crx and cry on the centroidal Z-Axis for z = B plotted with respect to Poisson's ratio . . . . Lines of constant 0x on the Y-Z plane with x = 0 for Poisson's ratio of 0. Problem Two . . . . . . . . . Lines of constant oz on the Y-Z plane with x z: 0 for Poisson's ratio of 0. Problem Two . . . . . . . . . Lines of constant O'x on the Y-Z plane with x = O for Poisson's ratio of 0. Problem Three ..... . . . Lines of constant oz on the Y-Z plane with x = O for Poisson's ratio of 0. Problem Three ........ Load distribution areas for Problem Four . . . . . . . vii Page 65 68 69 7O 71 72 73 74 75 77 Figure 28 29 30 31 32 33 34 35 Page oz stress distribution on the centroidal X-Axis with z = O for various values of Poisson's ratio. PartOne. ....... .............. 81 crz stress distribution on the centroidal X-Axis with z = O for various values of Poisson's ratio. PartTwo ...... .. ...................... . 82 oz stress distribution parallel to the X-Axis with y = O and for various values of z and with Poisson's ratio of 0.2. Part One ....... . ....... 83 oz stress distribution parallel to the X-Axis with y = 0 and for various values of z and with Poisson's ratio of 0.2. Part Two ...... 84 Lines of constant oz for the X-Z. centroidal plane with Poisson's ratio of 0. 2. Part One . . . . . . . . . . 85 Lines of constant Uz for the X-Z centroidal plane with Poisson's ratio of O. 2. Part Two . . . . . . . . . 86 Lines of constant (TX for the X—Z centroidal plane with Poisson's ratio of O. 2. Part Two ....... . . . 87 Assumed stress distribution acting on beam (a), nodal forces as determined by contributing areas, (b), for grid spacing of B/4, and nodal forces as determined by contributing areas, (c) , for grid spacingofB/8 ................. 89 Figure Page 36 Extreme fiber stress at x = B/2, y = O, and z = O for grid spacings of B/4 and B/8 for various values of Poisson's ratio . . . . . . ................ 90 ix CHAPTER I INTRODUCTION Advances in computer technology in the past few years have made possible the development of high speed, large storage digital computers. Problems which formerly required many hours of com- putation, or were limited by the available storage capacity, can now be solved with much less difficulty. Therefore, it is now possible to work problems, involving a great amount of calculation and a large amount of storage capacity, that would previously have been avoided because of the impracticability of solution. The study of three dimensional solids is an area that for this very reason has only recently been given attention outside of those problems which can be solved in closed form, or those problems which can be reduced to a less complicated level, i.e., one or two degrees of freedom instead of three, such as, problems in plane stress, plane strain, torsion, etc. This dissertation will investigate the stress distribution in three dimensional elastic solids using a special numerical method developed by the author. This dissertation will be limited to solid rectangular parallelepipeds which have at least two planes of symmetry at their centroid. This dissertation will also investigate the effect of Poisson's ratio on the stress distribution within the solid. The solution of problems in three dimensional solids has been considered by many authors . One proposed method of solution has been the application of a frame -work or a lattice type analogy as suggested by Hrennikoff (2), D'appolonia and Newmark ( l) as well as McHenry (7) . Another method has been to use a direct stiffness matrix for various basic solid shapes as suggested by Melosh (8). The method which will be presented in this dissertation will also be an analogous method based on a three dimensional grid-work system, except that this method will be concerned with the equilibrium of a region around a nodal point of the grid, as opposed to the classical approach which satisfies equilibrium based on the derivatives of the stresses at the nodal point and which is also based on the assumption that the stress condition at the point is the average representation of the region around the point. In essence, the method to be presented uses a solid model to represent the solid region around the point, and will be concerned with the equilibrium of this solid model. The solution of the three dimensional equilibrium equations is found by applying Lagrangian linear interpolation formulas . The resulting equations will be in the finitedifference form, except they will be representative of the solid instead of the nodal point and thus be appropriately called "Expanded Differences . " Chapter II will discuss the development of the theory used in the proposed method, and will compare this method to the classical approach when such a comparison will be helpful to the general understanding of the method. The technique of applying this method for the general solution of a problem will be developed in Chapter III. Chapter IV is devoted to the application of the method to a few selected example problems. These include, concentrated loads (concentric and eccentric), a distributed concentric load, and also the solution using this method for a beam—column with an applied moment and zero end shear. Chapter V summarizes the use of the proposed method. Appendix A indicates the programming in Fortran language for one of the example problems. CHAPTER II THEORETICAL DEVELOPMENT OF THE METHOD OF EXPANDED DIFFERENCES The three dimensional relationships between stress and strain have been known for many years. These relationships have been presented in detail by many authors (6), (9), (10) , (13), (14) and ( 15) . In the derivations which follow, it will be assumed that the material is isotropic, homogeneous and can be considered to be a continuum. Also, that Hooke's Law is an adequate representation of the relationship between stress and strain. It will also be assumed that the functions which will be used to represent the stresses and deformations within the body will possess piecewise continuity across the grid lines in the body. The three principal axes will be called X, Y and Z and the dis- placements of a point within the body will be called 11, v and w in the X, Y and Z directions respectively. The equations which relate stress to strain for three dimensions in elastic solids are: Bu Bu 6v crx-Xe+ZG-é-; Txy—G(5;+5;) 8v Bu 8w 2 + — : —-— — cry Xe 2G 8y sz G (az + 8x) (1) 8w 8v 8w 2 X + 2 -- : —-— —— 0.2 e G 82 Tyz G (82 + By) In the previous equations E is Young's modulus, v is Poisson's ratio, and G is the shearing modulus where, E G: 2(l+v) (2) x is Lamé's constant, vE =(1+V)(l-2v) (3' X- and e is the dilatation per unit volume of a point, 8“ Emir ‘ax By an (4) Equations 1 can be written in a somewhat more usable form for computation purposes when the following terms are defined, _ E " 2(1+v)(l-2v) E. (5) and, a=l—2v (6) Using these constants, equations 1 become, 6x : E'[(1+a)%3+(1-a)('g%+ '3'?” 0y = E'[(l+a)%§-+(1-¢)(%E’+ 3'31] (,2 = E'[(1+a)-g-E+(1-¢)(%E+ '33)] TxyzE’a(g_§+%}) (7) szzE'a (33+ ‘33) TyzzE'a(g-:-+g%) These equations relate the stress at a point to the displacement derivatives at the point. The exact displacement functions u(x, y, z) , v(x, y, z) and w(x, y, z) are usually not available as con- tinuous functions throughout the whole body except for a very few special cases, and thus, various numerical methods must be developed to attain approximate solutions. There are also series solutions available using the Neuber-Papkovitch equations, (9) and (11) . The displacement functions are derived (or approximated) by simultaneously satisfying internal equilibrium and compatibility of the body along with the external boundary conditions either in terms of stress or displacement. The equations of equilibrium are: EF=O x 23F =0 (8) Y 2F =0 z If we assume that u, v and w are continuous functions we can derive the following equilibrium equations of an infinitely small element (dx, dy, dz) as follows for zero body forces. 80x 87 x 872x 8x ‘+ By + 8z 2 0 Bay BTXX 8Tz ay + 3x + 8z=0 '9' 86 87 37 Substituting Equations 7 into these Equations 9, the equilibrium equations can be written in terms of the displacements which means that compatibility of displacement will be automatically satisfied. The se equations are: E' 33+ E'a Vzu ll 0 E'%:—+E'aV2v=O _ (10) E' 33+ E'a Vzwz o where, _..._a_i._§i, 82 — 2 2 2 8x 8y Bz V2 (11) These equations are generally known as the Navier-Stokes equations of equilibrium. These equations satisfy equilibrium at a point within the body and everywhere within the body, providing u, v and w are continuous functions in x, y and z. If u, v and w are not continuous functions throughout the entire body, such as in the case when the displacements in the body are specified only at points, as for example in a grid work of nodal points, then equilibrium is only satisfied at the nodal points where the derivatives of these discon- tinuous functions are obtained. These derivatives will only be accurate to the order of the equation assumed passing through the nodal points, and the equilibrium equations will be only accurate. to the degree in which these derivatives represent the true average conditions around the point. The accuracy of the final solution will then be dependent upon the accuracy of the displacement functions and the size of grid employed. The classical method of solving the equilibrium equations in terms of the grid work, is to substitute into the equilibrium Equations 10 the normal second order finite difference approximations for the needed derivatives (4) and then write the finite difference equations for every nodal point of the grid work in the solid. This will lead to a set of linear simultaneous algebraic equations which can then be solved by various methods (12) . The existence and uniqueness to solutions of the partial differential Equations 10 has been shown by Korn (3) and Lichtenstein (5), for first and second boundary value problems under general conditions. The method as outlined above is the classical method of solving three-dimensional problems based on the Navier-Stokes equilibrium equations. As mentioned previously, these equations only satisfy, in general, the equilibrium of a nodal point in the grid work of the body. It would seem that a solution that converges faster (with respect to the size of the grid) could be obtained by satisfying equilibrium within a region around the nodal point instead of simply at the nodal point alone . The basic philosophy of the method to be described stems from the fundamental viewPoint that the prime interest is in satisfying equilibrium throughout the entire body. It would seem then that since a grid work is being used to represent the solid body, the equilibrium equations can be written in terms of the forces in the body as repre- sented by the grid system, rather than using the derivatives of the stresses at each nodal point in the grid system, as is presented in Equations 9. The stresses at any point within the body are given in terms of the first derivatives of the displacement functions. Assuming straight line variation of displacement between nodal points, the first deriva- tive of one of these functions along a grid line and in the direction of the grid line between two nodal points will be a constant. Therefore, the point where this derivative can be best represented is at the cen- ter of the particular grid line segment. This idea can then be ex- tended to all three coordinate axes which leads to the concept that the equilibrium equations can be formed by using the summation of forces on the planes passing through the center of the grid lines and per- pendicular to them as shown in Figure 1. When there are no body forces, the sum of the forces acting on the planes surrounding a nodal point can then be set equal to zero for equilibrium. In the classical approach it can be shown that the equilibrium equations can be formed by representing the average stress on these planes as the stresses which exist at the intersection of the coordinate axes with these planes (Figure 2). In the method to be developed, the average 'stress will be computed from the expanded derivative in each 10 ' I l I I I, I, .I I I / ‘\ \ ‘\ \ ‘ \ \ I . \ \ \ I -— " x )x . ' ~ - . \ \. . ‘ \ M ‘ i \ \ \ ‘ \ ‘ ,_ \ ‘ \\\ \\' \ \ \\\ x‘ \ \\ I ‘ \ ‘ \ x \ \ \ ‘1 c I _ - Y X .. I . \.\4 \\ l"\ \ \ i \ \L \L ‘ \ \ \ \ \ M“ > b c - -’ —~k‘ . 2 E , } J . . > I . ‘ I I‘ A ‘\ b“) - - , , - ,. 1 ' ,7 , I l c r r ». l 1 A / I \. f ’ / / . i / J I ,r / / 4 7f 1 1' r4 7 a _- / i 1' A , ." I I A! ‘ - - I 0 T—- ’l I I II Figure l. Isometric view of basic model octants surrounding the reference nodal point, and also the rectangular parallelepipeds of which they are a part. ll (C’if/,"l,6) (O'hfi’h.’)) V (171,0,"1‘5) (h1,h2'h3 A r r 03/2 Figure 2. Position of stresses as used in Equations 9 for the classical finite difference solution of the three dimensional equilibrium equations. Z ( 02.9.2311 ‘ l V". -C ' , [(0,11ng Figure 3. Position of the centroidal stresses on the first octant, 1, of the model, and also the nodal points which determine the internal displacement functions within this octant. 13 quadrant of the plane, and the equilibrium equations will then be determined from these stresses (Figure 3). With linearity assumed in both cases, the total sum of the forces acting on any plane through the body will be the same, but the forces when resolved to the nodal point will be different, depending on the curvature of the stress pat- tern. The efficiency of the method will become apparent when com- paring the relative convergence of the two solutions. This method can be described by comparing the volume en- closed by the planes which pass through the center of the grid line segments around a nodal point to a solid model. The displacement functions for an octant of the model will be determined by the parallel- epiped of nodal points of which that octant is a part (Figure l) . Linearity will be assumed between the nodal points of this rectangular parallelepiped and will be used in determining the displacement functions. The displacement functions will be continuous within the parallelepiped, but they will be only piecewise continuous across the grid line boundaries . The model to be used in deriving the equilibrium equations will be, in general, composed of the eight smaller octants which surround the nodal point and each of which is one octant of the parallelepiped of nodal points which is used to determine the displace- ment functions. The dimensions of one of these smaller octants will be one —half of the dimensions of the parallelepiped used in determin- ing the displacement functions. The equations of stress within an l4 octant of the model will be determined from the first derivatives of the displacement functions within the octant. The resulting stress equations will be, in general discontinuous across the boundaries between the octants of the model because the displacement functions are, in general, only piecewise continuous across these boundaries. The eight octants which surround a nodal point will have forces acting on their faces (Figure 3) as determined by the displacements of the nodal points of which that octant is a part. These forces can then be used to determine the equations of equilibrium. The displacement functions within a particular octant will be derived by using an extension of the Lagrange interpolation formula within the boundaries of the eight nodal points needed for linearly describing the stresses and displacements between nodal points . Before the equations of displacements are determined for each octant, it will be convenient for future developments to define the method to be used in identifying the relative position of the nodal points surrounding the reference nodal point from which equilibrium equations will be formed, and also in identifying the eight separate octants surrounding the reference nodal point. The position of a functional value with reference to the centroidal nodal point of the model can be indicated by showing the function to be dependent on three arguments I, J and K which will be the vector components of the point in the x, y and z directions respectively from the reference 15 nodal point. These three arguments will uniquely describe this point. For example, u(I, J, K) would be a typical use of these arguments. This quantity would signify the u displacement of the point atT+ 3+ K with respect to the reference nodal point. A particular octant of the solid model can be identified by a set of vectors T, Tand K which are set equal to the vector components of the grid lines coinciding with the internal edges of that particular octant. There are then eight sets of identification vectors, one for each possible octant. If we number the octants (n) one through eight, and express the grid ,h,andh spacmgs hl 2 3 in the x, y and 2 directions respectively, a table of the eight sets of identification vectors for the octants can be formed as follows: TABLE I CC TANT IDENTIFICATION VECTORS n i" 3" K 1 h1 h2 h3 2 h1 h2 —h3 3 h1 —h2 h3 4 h1 -h2 --h3 5 —h1 h2 h3 6 -h1 h2 —h3 7 -h1 -h2 h3 8 -h1 -h2 -h3 16 These sets of identification vectors can be made into sets of unit vectors by dividing them by their respective grid spacings and thus, B1n : In /h1 u—i : —-I h B2n Jn/ 2 (12) B3n : Kn/h3 The x, y and z coordinates of a point within an octant can then be expressed in the following form: n — rln B1n hl _. < < n r2nB2nh2 0 .. r -1 (13) zn : r3nB3nh3 The x, y and z coordinates of nodal points governing an octant can be expressed in somewhat the same form by setting the value of r1, r2 and r3 to be equal to zero or one. The constants m1, m2 and H13 can beused to represent these values of r where, m1 = 0, 1 H12 = 0, l (14) m3 = 0, l The x, y and z coordinates of a nodal point for a particular octant can then be expressed as follows and omitting the vector notation for convenience, 17 xnml : rnl B1n h1 ynmz : m2 an ha (15) Znm3 : m3 B3n h3 The interpolation equations for the displacement functions between the nodal points will be found by using an extension of the Lagrange interpolation formula for equal intervals as found in Kunz (4). The development of the equation for a single dependent varia- ble and only one argument x is as follows. Suppose, f(x0), f(x1), . . . , . ., f(xn) are the functional values for arguments x0, x ...... , xn where the interval between the x's is 1! a constant h and, xm=xo+mh (OSmSn) (16) x = x0 .+ rh . . th . . . In finding an n -degree polynomial P(x) pas Sing through the n+1 p01nts (xo,f(x0)), (x1,f(xl)), ...... , (xn,f(xn)), let am(r) be an th . . n -degree polynomial that 18 zero (possesses roots) for all the tabulated arguments except xm, and for this argument it is equal to one, i.e., 0 forkrfm am(rk) ‘{l for k = m (17) From these equations it follows, am(r) has n zeros 0, 1, . . . . , m-l, m+l, . . . ., n and that am(rm) = 1. It is found from satisfying the se conditions that, l8 _ r(r-l)(r-2).. .(r-m+l)(r-m-l)..(r—n) — m(m-l)(m-2)...(+l)(-l)..(m-n) : (_1)n-m r[n+1] (18) m! (n-m)! (r-m) The equation for the interpolating polynomial P(x) that passes through all points becomes, 11 n-m [n+1] P(x) = >2 (,‘1' r ccx > (19) m. (n-m)! (r-m) m m=0 or, n P(x) = E a (r)f(x ) (20) m m m=0 also, P(xk) = akcrkmxk) = fuck) (21) The same development can be used to find the interpolation equation for a dependent variable which is a function of two arguments x and y. Suppose that f(x0, yo), . . .,f(x0, ynz), . . . ., f(xn1, yo), . . ., f(xn , ynz) are functional values for arguments x0, y0;. . .; xnl, ynl where the interval between the x's is a constant h and between the y's 1 is a constant h2 and, = < < xml x0 + mlh1 (O - m1 — n1) _. < < ymz‘yo‘LrIihz (0 "m2 ”‘2' (22) x = x0 + rlh1 y 2 Yo + rzhz 19 To find polynomial P(x, y), a function of two arguments, x and y which passes through all of the points one may proceed as follows. Let a (r , r mlmZ l 2 for all the tabulated arguments except x , y , and for this argu- ml ”‘2 ) be the polynomial that is zero (possesses roots) ment it is equal to one, i.e., F0 forqufm l amm (r1 , r2 )=< 0 forkzafmZ (23) 1 2 k1 kZ K 1 for kl =ml,k2‘=mZ In order to determine the two variable polynomial which will satisfy the above conditions, consider, 0 forklim 1 am (rl ) = _ (24) : V 1 R1 1 for k1 m1 0 for k # m a (r ) = 4 2 2 (25) ”‘2 2k 2 l for k2 = m2 It follows that if there is a polynomial that satisfies Equation 24 and another polynomial that satisfies Equation 25, then, am (rl)°am (r2) aamm (rl,r2) (26) 1 2 l 2 Equation 26 must be an identity because Equation 24 satisfies Equation 23 for the zeros in the x direction and Equation 25 satisfies Equation 23 for the zeros in the y direction and therefore the two multiplied together must satisfy the combined conditions of Equation 23. The 20 form of the polynomials for Equations 24 and 25 has already been determined in Equation 18. The polynomial a (r , r ) becomes, mlmz l 2 ( _1)(n1-m1)(nZ-m2)r1[nl+l]rz [n2+ 1] a (r ,r ) = ' ' ' ' (27) mlm.Z 1 2 m1. m2.(n1—m1).(n2 m2). (r1 m1)(r2 m2) The equation of the interpolating polynomial P(x, y) that passes through all of the functional values becomes, P(x, y) n1 n2 (-1)(n1+nZ-m1-mz)ri[n1+l] r2[n2+2] m1=0 m2=0 ml‘mz' (n1'm1" 'nz‘mz" (r1'm1' 'rz'mz' m1 m2 ”1 n2 _ E E l 2 1 2 also, P(xk . Yk') l 2 = a (r ’1‘ )f(x .Y )=f(x . ) (30) k I R1 2 1k1 2k2 kl k2 kl 33<2 If the dependent variable is a function of three arguments x, y and z, the polynomial that passes through all points can be determined in the same manner as in the last two cases. The subscript 3 will be used to refer to those variables in the z direction. The polynomial which possesses zeros at all of the tabulated arguments except for the unique set of arguments x , y , x becomes, m m m 1 2 3 a'rn m m (r1, r2, r3) : am (r1) .am (r2)° a'rn (r3) (31) 123 l 2 3 21 ) aLrn m m (r1, 1'2’ r3 (n +n +n —m -m —m) [11 +1] [n +1] [n +1] (-1) l 2 3 l 2 3 r1 1 r2 2 r3 3 (32) m1! m2! m3! (ml-m1)! (nZ-mz)! (n3—m3)! (rl-m1)(r2~mz)(r3-rn3) The equation of the interpolating polynomial P(x, y, z) which passes through all of the functional values becomes, P(X.y.Z) n1 n2 n3 - E Z 2‘ a (r r r )f(x ‘ .. _ _ . . :Y .z ) (33) ml-O mZ—O m3—0 mlmzm3 l 2 3 m1 m2 m3 In order to apply this polynomial to the particular problem at hand, there must be included the position of the octant, and the dis- placement function that is being determined. This can be done in the following manner. Because of linearity in x, y and z, n = n = n = l the octant identification number :3 II H. H the displacement function number then, un(X. y. Z) = f1n(X. y. Z) ll vn(X.y.Z) f2n(X.y.Z) (34) wn(X. y. Z) f3n(X. y. Z) 22 The general polynomial Pin(x' y, 2) that passes through all of the tabulated displacements for one of the above equations becomes, >1: 2': g (..1)'3 m1 “‘2" m 3)r 11(r -l)r2 (r2 -1) P. (X.y.Z)= __ _ _ 1n ml—O mZ—O m3-0 (rl-m1)(r2-m:) r (r -l) .(3___§___) fi(X m I Ynm I Z ) (35) r3 1an3. n 1 2 nm3 The derivative of this function with respect to x, y or 2 can be performed in the usual manner. For example, the derivative of this function with respect to x would be the following: l l 1 (3-m -m -m ) 2 apln - 23 2 E (-l) l 2 3 (r1 -2m1r1+ m1) 8x '— m =0 m =0 m =0 2 1 2 3 (rl-ml) (rz-m2)(r3-m3) r2r3(r2- l)(r3-l) n h fl(xnm’ym’znm) (36) 1 1 n 2 3 but, m _ 2 1‘ m1 therefore, 1 l l (3-m -m -3rm) BPln .. E Z Z (-l) l 2 rr2 r32(r -l)(r3-1) 8x m1=0 m2=0 m3: (rZ-m2Hr3-m 3) Bln . ' -—f h 1(xnm'ynm'znm) (37) 23 The derivatives of Equation 35, as can be seen from the exam- ple derivative in Equation 37, are functions of only two of the three arguments, and being that they are linear equations, they form warped plane surfaces with respect to these two arguments. If it is desired, the previous equations can also be shown in matrix form. As an example, Equation 35 is shown in matrix form on page 25. When Equation 35 is expanded, for a particular variable displacement function, as for example u, it will appear in the follow- 1, r2 and r3 and also I, J and K are the variable ing form, where r quantities for a particular octant from Table I, and with the subscript n omitted from these variables for convenience, u(r ,r ,r ) =u(0,0,0)+r 1__ 2 3 1[u(I, 0. 0)-u(o, o, 0)] + r2[u( 0, J, 0) -u( 0, 0, 0)] + r3[u(0, 0, K) -u(0, 0, 0)] + r u(I,J,0)-u(0,J,0)-u(I,0,0)+u(0,0,0)] (38) 1r2' +r u(0,J,K)-u(0,J,0)-u(0,0,K)+u(0,0,0)] 213' + r u(I,0,K) -u(I,0,0)-u(0,0,K)+u(0,0,0)] 1r3[ + rlrzr3 [u(I, J, K) —u(I, J, 0) -u(0, J,K) -u(I, 0,K) + u(I, 0, 0)+u( O, J,.D)+u( 0, 0, K) -u(0, 0, 0)] When Equation 37 is expanded in the same manner as in the above equation, it will appear as follows: 24 8u(r1, r2, r3) - B1n 8x h [u(I, O, 0) -u( 0, 0, 0)] 1 B1 r2 h“ [u(I.J. 0)-u rANQJAV G be 26 of equilibrium will be satisfied for the region of the nodal point. The conclusion can also be reached that since the body is being represented as a mesh work of nodal points with linearity assumed between the nodal points, the accuracy of the solution will be dependent on the size of the mesh chosen, since the linear lines will then be able to represent the true configuration of the displacement functions with less error. As was stated previously, the determination of the equations of equilibrium is to be developed on the basis of the integration of the forces acting on the surfaces of the octants which are formed by passing planes through the mid-point of the mesh lines (Figure 3). The average stress on any octant's surface will be the stress that occurs at the centroid of the octant's surface. This is due to the fact that the stress surface is a warped plane, and the average elevation of a warped plane surface occurs at the centroid of the area covered. The appropriate values of r for the centroid of the stresses on the surfaces of an octant are shown in Table II. To find the stress at one of these points, substitute the appropriate derivatives, as illustrated in Equation 39, into the stress function, Equations 7, and introduce the required values of the variables for that particular octant and then substitute in the values of r taken from Table II. TABLE H STRESS r1n an r3n 0 1/2 1/4 1/4 xn 0 1/4 1/2 1/4 yn 0 1/4 1/4 1/2 zn T 1/2 1/4 1/4 xyn T 1/2 1/4 1/4 xzn T 1/4 1/2 1/4 yxn T 1/4 1/2 1/4 yzn T 1/4 1/4 1/2 zxn T 1/4 1/4 1/2 zyn The positive direction of the total force on an octant should be made consistent with the major coordinate axes. The stresses are considered positive according to the usual conventions. Therefore, the vector direction of the stress must be determined, by multiplying the stresses by the appropriate B coefficient for the octant on which the stress appears. Thus, ——-b —-L 0' : 0' xn 1n xn .O'_.:E-.O' yn Zn yn 0' :E‘O’ zn 3n zn T—‘zfid'l’ xyn ln xyn “=? szn lnszn (42) -—_\ —l T :BT yzn 2n yzn rhz-d B. T yxn 2n yxn *l l B T zxn 3n zxn “ *l B T zyn 3n zyn 28 The total vector force acting on an octant can then be deter— mined by summing the individual average vector stresses by the area over which they act. This produces three total vector forces for any nun—A —L —-N . . . ' octant (Fx' Fy and F2) . If the grid d1v181ons are all the same, a, then, h :h :h :a (43) and the forces acting on an octant can be developed as follows: “ F __g ___‘ _‘ fl = 0' + T + T (44) 2. xn yxn zxn a F: 1 -——)—(Z— : E' R: [( l+a)['9U(O, 0: 0) '31.].(0, 0: K)-3U(O, J! 0) a -u( 0, J, K)+9u(I, O, O)+3u(l, O, K)+3u(I, J, 0)+u(I, J, K)] +(l-a)'2B B [-3v(0,0,0)-3v(I,0,0)-V(0,0,K) ln 2n -v(I, 0, K)+3v(0, J, 0)+3v(l, J, O)+v(0, J, K)+v(I, J, K)] +(I-a)'2B B [-3w(0,0,0)—3w(l,0,0)—W(O,J,O) ln 3n -w(I, J, 0)+3w( 0, 0, K)+3w(I, 0, K)+w( O, J, K)+w(I, J, K) ]] +aE' 1%; [-9u( O, O, 0) -3u( O, 0, K) -3u.(l, 0, 0) -u(I, O, K)+9u(0, J, 0)+3u(0, J, K)+3u(I, J, 0)+u(I, J, K) +2B1BZ [-3v( 0, O, O) —3v(0, J, O) -v(0, O, K) -v(0, J, K) +3V(I, O, O)+3v(l, J, 0)+V(I, 0, K)+v(l, J, K)] - 9u( 0, 0, O) -3u(l, O, 0) ~3u( O, J, 0) -u(I, J, 0)+9u( O, 0, K)+3u(I, 0, K) +3u( O, J, K)+u(I, J, K)+ ZBlnB3n[-3w( O, 0, 0) -3w( 0, O, K) -w( 0, J, 0) -w( 0, J, K)+3w(I, O, O)+3w(l, 0, K)+w(I, J, O) +w(I, J,K)fl 29 or, F" E, :31: TZZ' [—(9+Z7a)u(0, o, 0)+( -3+3a)[u(0, O,K)+u(0, J, 0)] +(-1+5a)u(0,J,K)+(9+3a)u(l,0,0)+(3+‘5a)[u(I,J,O) +u(I, 0,K)]+(1+3a)u(I,J,K)+2- BlnB2n['3V(o’ 0,0) +(3-6a)v(O,J,0)-v(0,0,K)+(l-Za)v(0,J,K) (45) +(-3+6a)v(I,0,0)+3v(I,J,O)+(-1+Za)v(l,0,K)+v(I,J,K)] +2-B B [-3w(0,0,O)+(3-6a)w(O,O,K)-W(O,J,O) ln 3n +(l-Za)w(0,J,K)+(-3+6a)w(l,0,0)+3w(I,O,K) +(-l+2a)w(I,J,O)+w(I,J,K)]-J Similarly, F E, .123:T6;[-(9+27a)v(0,O,O)+(-3+3a)[v(0,O,K)+v(l,0,0)] a +( -1+5a)v(I, O,K)+(9+3a)v(0, J, O) +(3+5a)[v(I, J, O)+v(0, J,K) ]+(l+3a)v(I, J,K) +2.B B [-31.1(0, 0: 0)+(3-6C)11(I, 0: 0)-U(O, 02K) ln 2n +(1-2a)u(I, 0, K)+( —3+6a)u(0, J, 0)+3u(I, J, O) (46) +( -l+2a)u(0, J,K)+u(l, J, K)]+ Z ' BZnB3n[-3w(0,0,0) +(3-6a )w(0, O, K) -w(l, O, 0)+(1-2a)w(l, O, K) +,( -3+6a)w(0, J, 0)+3w(0, J, K)+( -l+2a)w(I, J, 0) +w( I. J. K) 1] 3O and, f" Zn El . 2 =T6;[—(9+27a)w(0,0,0)+(-3+3a)[W(I,O,O)+W(O:J:O)] a +( -1+5a)w(l, J, O)+(9+3a)w(0, O,K) +(3+5a)[w(0, J,K)+w(1, O,K)]+ (1+3a)w(I, J,K) +2' BZnB +(1-2a)v(l, J, 0)+( -3+6a)v(0, O,K)+3v(0, J, K) (47) 3n[-3v( O, O, 0)+( 3-6a )v( 0, J, O) -v(I, O, 0) +(-l+2a)v(l,O,K)+V(I,J,K)]+2’B B [-3u(0,0,0) n 3n 1 +(3-6a)u(I, 0, 0) -u(0, J, 0)+( l-Za)u(I, J, O) +( -3+6a)u(0, 0,K)+ 3u(I, O,K)+ ( -l+ Za)u(0, J,K) +(u(I,J,K)]] These three equations represent the three principal vector forces acting on each octant. The general total force acting on an octant would then be the sum of these three vectors . The total force acting on the planes surrounding a nodal point would then be the sum of the forces on the individual octants. If the octants used are to be actually part of the physical solid body, then the number of octants surrounding a nodal point can vary from one octant for a corner, two octants for an edge, four octants for a surface to eight octants for an interior nodal point. The method for combining these forces into the proper equilib- rium equations will be discussed in the following chapter. CHAPTER III METHOD OF SOLUTION The general solution of solid problems using the octant force equations developed in the last chapter is obtained by simultaneously solving a set of equilibrium equations formed from the nodal points of the whole body. Each nodal point in the body, in general, has three degrees of freedom of displacement. There is available an equation of equilib- rium in the direction of each unknown displacement component. Therefore, there is one equation for each unknown displacement throughout the whole body, forming a complete and solvable set of equations with no ambiguity. In the interior of the body, the equilibrium equations are formed by summing the vector forces on all eight octants to zero in each of the principal directions. The equilibrium equations formed on the boundary, however, require fairly close attention, especially in those regions where the curvature of the displacements along the boundary is comparatively large. In the interior of the body the equilibrium equations are formed by essentially taking the difference of two first order equations, thus, forming a second order difference. On the boundary, there are two alternatives to the solution. The first method would be to eliminate the octants that are not actually 31 32 part of the body and use the boundary forces as part of the equilib- rium equations for the affected nodal points. This will then lead to a first order solution in the neighborhood of the boundary, which will lead to small error, for small curvatures of the boundary displace- ments. The second method would be to create external fictitious nodal points one grid spacing removed from the boundary. These fictitious nodal points can be determined by substituting the boundary stress conditions for their finite difference equivalents in terms of these fictitious points. The stresses on the boundary will be correct and the equilibrium equations will be of the second order throughout the body. The first method will be applied to the solution of the example problems in the following chapter because of their relative simplicity. Using the first method, the general equation for the equilibrium of a nodal point can be shown as follows: 8 8 F:O=ZTF+ZY (48) n: : 8 8 F=O=ZTF +Z‘Z n= : If there are no body forces and the nth octant exists, T =1 and n X = Y = 2n = 0. If the nth octant does not exist and there are no 33 body forces, then Tn = 0 and Xn, Yn and Zn are the appropriate boundary forces in terms of the boundary normal force and shearing forces. The identification subscripts which occur in the equations for the forces on the octants, Equations 45 to 47, are at this point only identified with respect to the reference nodal point where the equilibrium equation is being formed. For a complete set of equations, the reference nodal point for the equilibrium equation must be referenced to the main centroidal axis system of the whole body. This can be accomplished if we define the distance to the nodal point in the following manner, XL = La L = 0, I, Z, 3, ..... yM = Ma M = O, 1, 2, 3, ..... zN = Na N = O, l, 2, 3, ..... Therefore, the identification of the position of a nodal point for a particular displacement function can be indicated as follows: u(I, J, K) L, M, N = u(La+I, Ma+J, Na+K) v(I, J, K) L, M, N = v( La+I, Ma+J, Na+K) ( 50) w(I, J, K) L, M, N = w( La+I, Ma+J, Na+K) The general equations of equilibrium of a nodal point with respect to the main centroidal axis system then become: 34 8 8 z T F + z xn :0 “=1 nL,M,N an,M,N n‘l L,M,N 8 8 z: T F + >3 Y =0 (51) “=1 nL,M,N ynL,M,N n=1 n1., M,N 8 8 2 Tn Fz + 2 z :0 “=1 L, M, N nL, M, N “=1 nL, M, N In general, the system of equations produced for a given body with even a fairly small number of grid spacings, would yield a tremendous number of equations to be solved simultaneously. This would then eliminate as a possible method of solution the methods using the inverse of the coefficient matrix. The remaining choices are those which use relaxation and iteration. These methods are comparatively simple for the types of equations being used for solu- tion, because there are only 57 basic types of equations. The only variables of this system of equations will then be the location of the reference nodal point (L, M and N), the type of equation to be used (one of the 57 basic types) and the boundary forces, if present at that particular nodal point. The equations can, therefore, be written with these variables built into them. In order to converge to a solu- tion by iterating on the equations thus formed, a relaxation or even an over relaxation procedure can be adopted. If this is done, the relaxation of the equations can proceed until the change in the values of the displacements between-iterations is a small acceptable value. 35 This value represents the convergence error of the iteration solution. When this method of iteration is applied to the previous general equations of equilibrium, the equations can be represented in the following form, 8 8 R(L,M,N)=ET F +23x x ml n1.,M,N xnL,M,N n=1 nL,M,N 8 €- 2 T (9+27a) (52) n=1 nL,M,N u(L,M,N) =u(L,M,N) + c R (L,M,N) (53) new old r x 8 8 R(L,M,N)=ET F +2Y y n=1 nL,M,N ynL,M,N n=1 nL,M,N 8 E- E T (9+27a) (54) n=1 L,M,N v(L,M,N) =v(L,M,N) + c R (L,M,N) (55) new old r y 8 8 Rz(L,M,N) = z Tn an + 2: zn n=1 L,M,N L,M,N n=1 L,M,N 8 .1. >3 T (9+27a) (56) ml nI_.,M,N w(L,M,N) =w(L,M,N) + C R (L,M,N) (57) new old r z Where, the acceptable error 36 =e 2Rx(L,M,N) 2Ry(L,M,N) 2Rz(L,M,N) (58) and, Cr is the over relaxation constant that can be used when desired. When the error in the displacements has converged to the small error, e, the stresses at the nodal points can be determined using Equations 7. These stresses can then be used for the usual purposes and applications . CHAPTER IV EXAMPLE PROB LEMS The example problems that will be shown in this chapter are for the purpose of demonstrating the application of the method and the behavior of the equations involved under different types of loading conditions. The first is an example problem which will be used to compare the stress distribution in a solid cube supported on a plane surface under the loading of a line load acting perpendicular to the top sur- face along the centroidal X-Axis (Figure 4) when computed with the proposed method and Poisson's ratio equal to zero, to the stress distribution computed using the plane stress and plane strain classical solution of a thin plate, with a concentrated load acting at the center of one edge. The classical solution will be solved using the Airy's stress function (15) and finite difference methods. The stress distribution for various other values of Poisson's ratio will also be shown using the proposed method. The stress distribution, for Poisson's ratio equal to zero, is unique because the normal stress is dependent only on the normal derivative. Thus, when the boundary conditions are constant in the direction of one of the principal axes, the stress distribution on any plane perpendicular to that axis is a constant. When Poisson's ratio is not zero, the normal stresses are 37 38 dependent on all three normal derivatives and the stress distribution will vary along this axis . The second problem will consider the stress distribution in a solid cube due to the application of a concentrated load acting in the direction of the Z-Axis at the intersection of the centroidal Z—Axis and the surface of the cube (Figure 5) . The effect of Poisson's ratio on the stress distribution will also be investigated. The third problem will find the stress distribution in a solid cube due to a distributed load over a small concentric area on the surface of the cube (Figure 6). This stress distribution will be com- pared to the solution determined in problem two. This problem will also demonstrate whether some of the effects of the changing of Poisson's ratio are due to the nature of the concentrated load, or due to the behavior of the equations when the ratio approaches a value of one -half. When Poisson's ratio is one —half, the dilatation of the material is zero, and in such a case, the average stress for a point in the body is indeterminate, and only the differences between the average stress and normal stresses are determinate. Therefore, a solution of problems using the displacements of the body is impossible. The fourth problem will show the stress distribution in a beam- column due to the action of an eccentric load (Figure 7). This will correspond to the case of a post-tensioned member, or a column with an eccentric load. 39 The fifth problem will find the stress distribution in a beamw column subjected to the theoretical stress distribution of a beam under constant moment (Figure 8). Since the solution of this problem is already known, the problem will demonstrate the relative error associated with the use of this method's first order equilibrium equations at the boundary surfaces in conjunction with a variable Poisson' 3 ratio. 4O 00)“: >< Figure 4. Isometric View of cube with line load for Example Problem One. 41 >< Figure 5. Isometric view of cube showing concentrated load on the Z—Axis for Example Problem Two. 42 Til; I r)‘\f’*——-—————*N \\ "\ 14)“ NF Mm X Figure 6. Isometric view of cube showing concentric distributed load over small portion of the top surface for Example Problem Three. 43 t I a 9-,“), ,(" e. , 1 l a '9 l I 1 . u | I ' I i 1 l 1 i s E 4 m : ‘V : Z i / / allN 1%....»‘tv ! /, . / . _) -1. / / . —- ——~ - — y/ / ..=. / Figure 7. Isometric view of a beam-column showing the location of the eccentric load for Example Problem Four. 44 .— efi._ . 28 ,5 i l I i L Figure 8. Isometric View of a beam-column showing the surface stress distribution for the applied moment for Example Problem Five. ‘ 45 EXAMPLE PROBLEM NUMBER ONE In this first problem, we will find the stress distribution in a cube due to a line load acting along the X-Axis on the surface of the cube as is shown in Figure 4. The solution will be found for three different grid spacings. The grid spacing in any particular solution will be the same in all three principal directions. Each solution using the proposed method, with Poisson's ratio equal to zero, will also have a parallel solution computed using the classical biharmonic approach using finite differences. For the first grid, we will use a grid spacing of B/4, for the second grid, we will use a grid spacing of B/8 and for the third grid, we will use a grid spacing of B/l6. For purposes of comparison between the various solutions, we will specify that the average stress on any X-Y plane be unity. Therefore, the computed stresses will be stress concentration factors for the assumed unit average stress. If we consider the line load to be distributed over a one -half grid spacing on either side of the loaded line, then we can develop the following relationship between the speci- fied unit stress on the cross section, the stress on the distributed line load, and the size of the grid spacing. a = grid spacing B = dimension of cube N = number of segments in B length The refore, 46 B = Na Q = stress on loaded area then, QaB = 1B or, o = 113- = 1 - N a By inspection, the w displacements will be symmetrical with respect to the centroidal Y-Z and X-Z planes on any X-Y plane, and the u and v displacements will be anti-symmetrical with respect to these same planes. The numerical solution is, therefore, required in only one quadrant of any X-Y plane. All boundary forces can be set equal to zero except for those nodal points which lie along the loaded area. In the equations of equilibrium, for the nodal points, the boundary forces can be expressed in the equation in terms of the boundary stresses because the grid spacing is a constant. Therefore, the equations can be divided by the area over which the stresses act. It will be convenient for programing purposes to write the equations such that the boundary forces can be included in the equa— tions. This can be done by entering the equation with the proper nodal displacement incremented by some value and when leaving the equation subtracting this incremental value, leaving the nodal displacement at the correct value. This technique and the method 47 of determining stresses can be demonstrated as follows: As an example, part of the equation of equilibrium for a quadrant on the Z surface (Equations 47 and 48) can be written in the following form, E! boundary _ 16a [-(9+Z7¢)W(0.0.0)+ etc.'°"]+ Kcrz -— Owhere K is an arbitrary constant for the boundary stresses. Defining, __ E'(9+ 27a) ’ 16a _ E'( 9+ 27a) - 16a __ E'(9-l-Z7a)u - 16a u! the equation can be written in the following form, E-(94'27a )w'(0, o, 0)+ . .. etc. ]/(9+ 27a] +bi°undaryz 0 Setting w'(0, 0,0) = w'(.0, 0,0) = w'(0.0,0) - C then C : Kcrzboundary as a correct substitution. When the equation has been relaxed (Equations 56 and 57), C can be eliminated as part of the w displacement, w(O, O, 0) = (w(O, 0, 0) -C) + C This is a very simple means of expressing the boundary forces, since C can be chosen as any convenient constant consistent with the 48 boundary stress distribution. When the stresses are determined from the primed values of the displacements, the equations must be corrected for the value of K and for the elastic constants. As an example of the procedure, the az stress at a nodal point L, M, N 'would be determined in the following manner: <7 (L,M,N) = E' [(1+a)9-‘3+ (1-a)a_“+ .31. z 32 8x By E' 168. aw' 311' 8V. : -—- —— 1- —— ......_ K E'(9+27a)[(1+a)az+( a)8x +3y The normal finite difference equations for these derivatives, Kunz (4), would be substituted giving the correct stresses for the given loading condition. The solution to the problem as a plane strain problem for Poisson's ratio equal to zero, was found by applying the familiar biharmonic solution as explained by S. Timoshenko and J. N. Goodier (15) , 2 2. 2 Y 32 By Y Byaz which satisfies equilibrium for zero body forces . Then, v4¢ = O for compatibility. For the boundaries, or, . 2 2 4) :y2?-+za—- “Q'fy .33- ‘ fir?) 23 "J A ) ____._-._..-._ . 7 l l l l r i if.) ' '3‘). ( l M! g .. ll I "\ \\ x.-. . ) (’3 N) 3 , : , l \ [I l t ,7 <3)”- : w - ) .., i l . l3i ”3.3T \ a .. y . i H R I_,’ ,4, r\ ‘ . -.4_— m“. -fll—b'---—w-"—- -«<—.- ”a- - Figure 9. Outline of a thin plate with a concentrated load and the associated boundary phi functions for the biharmonic solution for Example Problem One. 51 Poisson's ratio and for the grid spacings of B/8 and B/16. The values of oz at z = 15 B/l6 for the same parameters are shown in Figure 13. From these figures, it can be seen that as Poisson's ratio approaches one -half, the divergence between the solutions of the two grids greatly increases even when comparatively far from the loading surface. This would seem to indicate that the convergence of the equations to the true solution with increasing Poisson's ratio is not dependable for large grid spacings. It will also be noted, that within the usual range of Poisson's ratio for structural materials, 0.1 to 1/3, that the discrepancies between the solutions for the two grid sizes are small when comparatively far from the loaded surface where the singularity exists. For solutions of the stress distribution near the singularity, increasingly small grid spacings would be necessary. This is clearly evident in Figure 13 which shows the crz stress at z = 15B/l6 for various Poisson's ratios and the two smaller grid sizes. The values for the larger of the two grid spacings were computed on the basis of straight line interpolation between z = 7B/8 and the assumed distributed stress on the top surface. The lower line represents an assumed stress over one grid spacing of B/8 or a stress of 8 at the top surface and the upper line for the smaller grid length of B/l6 or a stress of 16 on the top surface. The Flamant (15) solution for a concentrated load for z = 15B/16 and z = 7 B/8 also is indicated in Figures 11 and 13 for Poisson's ratio of zero. The 52 agreement between the Flamant solution and the numerical method solution seem to indicate that as far as the vertical stress distribu- tion is concerned, the solution in the neighborhood of the concentrated load is the Flamant solution. However, this will not be true for the 0x and (TY stress distributions. The 6x stress does not exist as part of the Flamant solution because the problem is assumed to be that of the solution of a flat plate in plane strain or plane stress. In the Flamant solution, the 6y stress is always of the same sign. In the real problem under investigation, we are not dealing with the semi- infinite plane and it is required in this problem that the sum of the normal stresses on any vertical section be zero, and therefore, the cry stress will be of varying sign so as to satisfy this condition. The Flamant solution for 0y on the Z-Axis is zero and the 0.x stress on the Z-Axis would be the oz stress multiplied by Poisson's ratio (assuming the plane strain solution). Figure 14 shows the 0.x and or stress on the Z-Axis for a Poisson's ratio of 0. 3 and a grid spacing of B/16. The Flamant solution can be used in the neighborhood of the concentrated load for the vertical stresses because the boundary conditions on the surface where the load is applied agree in both cases. The Flamant solution can not be used for the horizontal stresses because the boundary conditions on the vertical planes do not agree. The relative importance of this consideration would be 53 unimportant for design purposes in materials which are equally as strong in tension as in compression.“ However, in materials which are not equally as strong in tension as in compres- sion, such as concrete, the presence of tension in the material could be of concern. The presence of these tension stresses will be even more marked in the example problem to follow which is the case of the concentric concentrated load. As an example of the general distribution of the stresses throughout the block, lines of constant stress for crx, 0'y and O'z on the Y-Z plane with x = 0 are shown in Figures 15, 16 and 17 respectively for the grid spacing of B/16 and with Poisson's ratio of O. 3. .v .N eofimfiommhuxo 80.3 poudeoU as“. 350m dungeon “83.309" 3“: «ES ”Emwwuwm gosh a. 54 woo.: mom.: >531. owmd >.mm** combs. o.mm** 380$”... Qnm m MIN owmg: $4.2 meat. 2&5 . 93.... Mame... New... 8:... 94> m: mmm.m Swim mi. 3&5 8.2. 284. 0.3 2%.? 9m m m .33.... :35 9.1. moat... w.$ m3... mic... mama... as.» a mood coed 9:: god oi: 2:8 52. come 9m N . m $8M $8M 5:: 238. 0.2: $5M Too SYN a...» m 284 mi; 90.3 «3.4 who“: :8; ~53 cum; 9m Em $84 New; 82: Re; 8:2: is; of: €34 94> $0; $0; 0.2: moo; ~12: $0; A.12: 234 GA o moo; So; «.2: 3o; m.:: to; 52: M24 as». as .uuxm axe .Huxm axe .uuxm aim are; Nine 35 nine Em nine 3m @2302 N Gofldaom N poumfiommtxm b mswm nmoz N mZOHHDJOm 03H. ”NEH. MOh mMNHw mmmuz mDOHm<> 39:5 b .m0 wHDA<> HHH HA m——-/~’ ,9 V - ’ ‘ C’J“ '—--~ID-~-—— — -4 ~"""¥( Cir“ ‘—' '1" .. -'- r --~' I I -*—"“ — IV 0 0.1 0.2 0.3 0.4 O. 5 Figure 20. Values of 62 on the centroidal Z-Axis for z = 15B/16 plotted with respect to Poisson's ratio. 70 I? 60 ~ ______ grid spacing B/8 Prob. Two / grid spacing B/l6 Prob. Two / ---—--- - grid spacing B/l6 Prob. Three x" q _ 30 — .-4/ v/u , ,/ . z = 7 B/8 ./ ‘ A / I , , , , ’ , ' . ”I, F /'/I/ ‘-_I‘ x. ”’23) I O 5L1: :w- r“ "if; 4") - - d: T ’— * "’ 1 U '0 O l O 2 0.3 0.4 0 5 20 *i i 10 1‘ . l _ l I q // I / o :7.» .. - - ! ._.. —-—-— 7 ~ "w: ore-Arv—OOSII/ O 0.1 0.2 0.3 0.4 If; 2 g 1 4 l , i 0 J ’ " J" r "" ' ‘ *' ' — ' ‘ -- __‘; _:‘.'_- -- 111244;— U . -.. - . --__..__....___.... . ._ — ‘- ._ _ ‘2. _.-. _..:..._.... 4.]. —-”’" 0 0.1 0.2 0.3 0.4 0.5 Figure 21. Values of “x and cry on the centroidal Z-Axis for various values of z plotted with respect to Poisson's ratio. I 71 I. -------- grid spacing B/8 Problem Two i grid spacing B/16 Problem Two r}, I ——- - -—-——- grid spacing B/l6 Problem Three / i — -— - ——- Dist. Theory B/8 ,' l / // ——- - - - --— Dist. Theory B/l6 I, , / / I, l / I 200 / /J //‘ I / z“ I/ I // / . 100 «’ /’ , z = B / . /" I ITO/VI / V. // ’ riaflj "f” /” , - - ““’ {/3/ 1’ ‘. - "y "' f?” ’ f/ 1.," ’d_ .. ‘ ”r- ’/.:(:" / M.._ ._._—-——- - -~“\fI ,, ’ .. ’ ‘ Q" ‘ .. _. .. o-—+-- --.-.----....- - ». -, . ‘Z/ 0 0.1 0.2 0.3 0.4 0. 5 Figure 22. Values of 0x and cry on the centroidal Z-Axis for z = B plotted with respect to Poisson’s ratio. m _.—.—.__..- +..—___..._.-_.--7—. _ ._ -.._ __._ ...._ ..__ -‘_..__.-. __-‘._-.---_.__. -_- ...-..._ .-_._.. 1 Figure 23. Lines of constant “x on the Y-Z plane with x = 0 for Poisson‘s ratio of O (+ = comp. ) . Problem Two. Figure 24. Lines of constant crz on the Y-Z plane with X: 0 for Poisson‘s ratio of 0 (+ = comp. ) . Problem Two. *--—_—----_ -- .-- --_r.-.. . _.---- A I T Figure 25. Lines of constant “x on the Y-Z plane with x = 0 for Poisson's ratio of 0 (+ = comp.) . Problem Three. Figure 26. Lines of constant 0'z on the Y-Z plane with for Poisson's ratio of 0 (+ = comp.) . Problem Three. II 76 EXAMPLE PROBLEM FOUR In this problem (Figure 7), the stress distribution will be investigated in a beam column which is acted upon by an eccentric concentrated load. The problem will be divided into two parts. In the first part, the load will be considered to be concentrated on a square area centered around the nodal point, whose dimensions will be equal to the grid spacing which for this first part will be B/4 (Figure 27a). In the second part, the load will be considered to be concentrated on a rectangular area centered around two consecutive nodal points along the upper surface X-Axis, whose dimensions will be one grid spacing in the Y-Axis direction, and two grid spacings in the X-Axis direction. The grid spacing for this second part will be B/8 (Figure 27b). The assumed eccentricity in the first part will be B/4 and in the second part 3 B/16. Symmetry will be assumed to exist across the centroidal Z-X plane, therefore, eliminating one -half of the structure for solution. As in the previous problems, the w displacements on the plane z = 0 will be set equal to zero. Displacements of nodal points lying on the plane of symmetry will also be assumed zero in the normal direction. The equations of equilibrium and method of solution will be handled in the same way as in problem one. The solution will in each case be computed for Poisson's ratio of 0.0, 0.2, 0.3 and 0.4. The solutions for values of Poisson's ratio 77 s- - —-~ )— Y . _ LVII /// I. . WM/hef II- -. ”—— .— -7— ... . . m . . . p . _ _ _ I II...IIL-.II-I . I Us _ _ ,- m _ + M - n 4 4. n M .I _ _ _ h n . _ m _ m III_-I IIIIIPII h ”1. I_- II- n _ . I . I. I I . _ c A . . I __ . m . ,_ vvvvvv ilil llfl (‘ %“ l Ll loll. “I IIIIII4I'I- -- .TIVIVA Load distribution areas for Problem Four. Figure 27 . 78 nearer 1/2 will not be computed because of the instability of the equations as evidenced in the previous problems. At a distance from the load, the stress distribution in the beam column will be governed by the equations: ”13.1322 Uz‘A I or, P 2 UzzA (1 leZX) B If we substitute into the equations above the value of the eccentrici- ties for part one and two, the stress distribution across the beam would be: for part one: P 3x (I. 2: (1+?) for parttwo: P 9x crz"A (1+Z—B) In both parts, the average stress on the cross section of the beam column was set equal to 500. The stress distribution associated with this average stress would be the computed stress distribution from the above equations if the beam column was infinitely long, and if there was no error introduced because of the approximate nature of the equations . 79 Figure 28 shows the axial stress along the X—Axis in part one with 2 equal to zero and for Poisson's ratios between 0.0 and 0.4. The error in the extreme fiber stress with respect to a true straight line variation varies from -4%to +2070. Figure 29 shows the same stress distribution for part two. The error in the extreme fiber stress for this case varies from -2%to +2% which is a smaller range than in part one. The theoretical stress distributions for part one and two are indicated in Figures 28 and 29 by the broken lines. These stress distributions are based on the stress distribution in the infinite beam column. Figures 30 and 31 show the vertical stress distribution on the centroidal X-Axis for part one and two respectively with Poisson's ratio of 0. 2. Only the distribution in proximity to the loaded surface is shown since the lines from the base to approxi- mately B/2 are parallel. Figures 32 and 33 show the lines of constant vertical stress on the centroidal X-Z plane for part one and two respectively with Pois son's ratio of 0. 2. In both cases, it can be seen that the trajectories are virtually vertical between the base of the column and the center of the column. Figure 34 shows the lines of constant stress 6x on the centroidal X-Z plane with Poisson's ratio of 0.2 for part two. The maximum tension stress is approximately one -tenth of the maximum theoretical compression fiber stress or two -tenths of the average stress which would be of considerable interest in design practice if the material happened to be concrete. 80 It can be seen from these figures that the stress distribution agrees with the theoretical stress distribution when comparatively far from the loaded surface. The error varies with respect to Poisson's ratio and the range of the error decreases with decreasing grid spacing. Tension stresses that are small in comparison to the compression stress exist and may in some circumstances cause Concern . Error 81 Q / q “5% @\\\f— Theoretical 1200 ‘ U = C13 - . v = o 800 600 \-\ - \\ 400 200 m-_-..--l..nm- n1-.. _--- , I I _ F. -___._,.._---_.\.\§ei_-----.-_._.-, X - B/2 - B/4 I B/4 \\-\«\.;\ B /2 i \\\ I \:\‘~ . j \ \ \‘\ ‘\\\-‘ - zoo _ \o I {gt/l Figure 28. oz stress distribution on the centroidal X-Axis with z = 0 for various values of Poisson's ratio. Part One. Error +-:°I- F\ I/—— Theoretical \w -570“. \\‘ —-—~—U:C . \‘C. “L/: C l} \\ \ \S; r" -‘~~y- -—-.~-— ~'--r - - B/Z - B/4 82 i I 1200 I I IoooI I I I I sooI I I I I - 200“ Figure 29. 0' stress distribution on the centroidal X-Axis with Z 2 = 0 for various values of Poisson's ratio. Part Two. 83 5000 1 I" I I I I 1 I ‘ 4000 I I" I I ' I III I I I 1 3000 « ,' z = 14 B /. 8 I /— I . I I 3’ -z=6B/4.I 2000 4 .' [:I"O\\I :'\‘:\ z = 10 B / 8 I I ‘ . .- I I I §‘\ i‘Xeq/fi ‘1 O “I .’ ‘.'\‘ II ‘ l 5 “\‘Ir "r4 -_\¢I I I"'” W ‘“ s I fTTL‘“? ‘\...:’;" - B /2 " B /4 E ’1 v I I - 1000‘ Figure 30. oz stress distribution parallel to the X-Axis with y = 0 and for various values of z and with Poisson's ratio of 0. 2. Part One. 84 CY; I12,000 +10,000 Xx, / 1 4 8,000 l J ‘ I 6,000 4.913. 16 I ! [3"‘\.I I 4,000 \I \a 1 ,14B/s\ I\ // \\ 2,000 - \. fi__ -_ .1. , IV\ ?‘\\\: , T \4‘ if; ‘x I - B/Z - B/4 B/4 13/2 Figure 31. 02 stress distribution parallel to the X-Axis with y = 0 and for various values of z with Poisson's ratio of 0. 2. Part Two. 85 Lines of constant oz for the X-Z centroidal Figure 32. Part One . plane with Poisson's ratio of 0.2. C 5 5 ,. J J #00 303 201‘) II if! {I » ' 1 ‘ I I i 2 f i ' - ' ‘ I ‘ . 1 ‘ .__._..A_ _._ V . ..-._l _ . , ,.L__n ._,._.A-. -_,j.._-,it___4 Figure 33. Lines of constant oz for the X-Z centroidal plane a with Poisson's ratio of 0. 2. Part Two. _- _‘--- I C‘) c l ! l i 5—--- --- .-., ,..._ - . .. -. _..- .. - .. . .. . . A Figure 34. Lines of constant (TX for the X-Z centroidal plane with Poisson's ratio of 0. 2. Part Two. 88 EXAMPLE PROBLEM FIVE In this problem (Figure 8), the accuracy of this proposed method will be studied by comparing the known stress distribution in a beam under constant moment and zero shear to the solution which results from the numerical method using the first order boundary equations. Figure 353. shows the assumed stress distribution on a line parallel to the X-Axis on the upper surface of the beam. Fig- ure 35b and Figure 35c show the assumed boundary nodal stresses for grid spacings of B/4 and B/8 respectively. Symmetry was assumed to exist only with respect to the X-Z plane. The solution will be handled as in the previous problem. Figure 36 shows the outside fiber stress for solutions with various values of Poisson's ratio for the two grid spacings. In Figure 36, the error associated with a base stress of 500 is also shown. Within the usual range of Poisson's ratio (0.2 to 0.3), the error ranges from -Z.761% to -O. 782% for a grid spacing of B/4 and from -0.38% to +0.870% for a grid spacing of B/8. As in the previous problems, the accu- racy of the solution is dependent on the value of Poisson's ratio, and instability of the equations becomes progressively worse as Poisson’s ratio approaches 1/2. 89 - /‘-{/1’_/ +/CC (a) -7 + 7 ._ ~ I , 57:- i “4 +4 +? . 1“?“ I i A . l I ' >~~\I. - - ‘ - .: a ..._ \ 7 I ' f ' 5' I C..- ...__ ' I ‘ __ Y- - - l -..l . l. - - _ i. J (b) -/5 HS _/( . . A F—‘Lo . "HG I . l V,» I ' r‘x'xq I -£f 0' +4 I -’ It”. ..: l I ; I "\.\\ l ./ I i ‘ ‘““ : A 177%: l was ' 3 ' 3 ’ T‘TSJ - . -- :"A' i ' I ’6 i ~ ‘ l I ?\\- I k/ ‘ I l ‘I I -y _ I. -1 -_I.---Z..\...L.2-:'r.-l 1--.--.“ .- ..-_..l (C) Figure 35. Assumed stress distribution acting on beam (a), nodal forces as determined by contributing areas, (b), for grid spacing of B/4, and nodal forces as determined by con- tributing areas (c) for grid spacing of B/8. 90 1 Us 1300 I x = -B/2 I ‘f y = 0 I]. 1 f I z = O i l I‘ I/ ‘ I 1200 4 I/ I I!" I 5/ I f," i [If/I i Error I/ 1100 «I + 10% I I l i I" + 5% I . Theoretical value 1000 '4 -- _,..,..« “— /A'-”- 900 1 n - V —-'I2 0.4 0 5 0.2 0.3 :B/Z, y=0, andz=0for Figure 36. Extreme fiber stress at x grid spacings of B/4 and 8/8 for various values of Poisson’s ratio CHAPTER V SUMMARY AND CONCLUSIONS A numerical method for the solution of the equilibrium equations for a three-dimensional body in terms of its displacements has been presented. The method is designed to satisfy the equilibrium equa- tions for a region around a nodal point of the grid system employed to represent the body. The applicability of the equations was studied by applying them to various representative problems. The example problems for this study were chosen to show the disadvantages as well as the advantages of this method. The problems were also solved without any use of refinements that would tend to give better results. This objective was accomplished by first, not making use of fictitious nodal points to satisfy the boundary equations, which leads to first order equations with respect to the boundary, and second, by solving problems with concentrated loads on the boundaries which leads to singularities. Since, solutions of problems in three—dimensional elasticity which make use of the displacement functions are impossible when Poisson's ratio is one-half, because the dilatancy of the body is zero, the solutions were also studied with respect to this ratio approaching one-half. In problem one, when U = 0 it was found that the solution of the numerical method presented was better than the 91 92 classical solution with respect to grid size. It was also found that as Poisson's ratio approached one -half, the solutions became unstable. In problems two and three, it was found that the instability of the equations was a function of the grid size. This is clearly evident because the solutions for the fine grid in problem one and problem three agree very well when not in close proximity to. the load, but the coarse grid diverges sharply from this solution when Poisson's ratio approaches one -half. In problem four, the solution of. the concentrated eccentric load, when compared at a distance from the concentrated load, agreed with small error to the theoretical solution of the infinite beam column, which agrees with St. Venant's theorem. The maximum tension stress was approximately 20% of the average stress. In problem five, it was found that in using the first order boundary conditions, the error within the usual ranges of Poisson's ratio was comparatively small. In reviewing the results for these problems, a number of conclusions can be drawn: 1. A generalmethod has been presented which can be used to solve problems in three dimensional elasticity. Up to this time very little literature has been available concern- ing the numerical solution of stress distribution problems 93 in three dimensions. The method is applicable to all problems that can be described by a suitable mesh, and is limited only by the capability of the computer. This method includes as a possible variable the use of Poisson's ratio. The primary influence of Poisson's ratio is in its effect upon the stress distribution in the region of a concentrated load. When the equilibrium equations on the boundary use the boundary forces rather than the boundary displacement derivatives, there is a secondary effect of Poisson's ratio that is introduced into the solution. If Poisson's ratio is limited to values less than 0. 4, the accuracy of the solution is good. 3.’ The solutions appear to converge monotonically with respect to decreasing grid spacing. It has also been shown that when l/ = O the method presented in this thesis compares more favorably to the exact solution than the classical numerical solution with respect to grid spacing. There are a number of refinements that can be included in the solution of a problem in three dimensional bodies that were not used in the presentation of the method. First, fictitious points once re— moved from the boundaries could be used to satisfy the boundary stress conditions. This would then lead to equilibrium equations on the boundary of the second order rather than that of the first order. 94 There are, however, a number of drawbacks in this approach. At the intersection of surfaces, there are more unknown nodal displace- ment values involved in the fictitious points than there are stress equations available. This leads to approximations as to the curvatures of the boundary surfaces at these points, and thus some error may be introduced. Second, the well known subtractive process can be used to eliminate singularities from the numerical solution. This also would present some difficulties. The boundary conditions to be used in such a process would involve defining the displacements on some surfaces, and the stresses on other surfaces which would lead to some compli- cations in writing the computer program necessary for the solution of the problem. It would, however, lead to a very accurate solution in the region of the singularity. There are many avenues for the future development of this method. The application of the refinements to some problems should be attempted to evaluate their potential as useful inclusions in a general solution. The application of the use of these expanded differ- ences should be attempted with plate problems to determine their applicability in this situation. In conclusion, a numerical method for solution of three dimen— sional elasticity problems has been presented, and by all implications 95 seems to be a method which is easy to use, and which gives good results for comparatively small grid spacings for the usual structural materials . VI APPENDIX FORTRAN COMPUTER PROGRAM FOR EXAMPLE PROBLEM FIVE 96 97 pnoann BEND FORMAT STATEMENTS 500 F09MAT15X.16HPO1SSON-S RAT10-.F7.41 100 FORMAT15x.19HNuusER OF CYCLES - .15.5x.7HRMAx . 0614.8) 101 FORflAT(315) . 102 F0RMAT12X.312.2X.91E11.5.1X)1 103 FORMAT (BX-5H! J K.6x.2Hu3.10x.2Hv3.10x.2Hw3.10x.3stx.9x.3Hsvv. 19x.3Hszz.9x.3HTxv.9x.3Hsz.9x.3HTv2.14./1H01 104 FORMATtlHloSXml9HPROB FOR H A ELLEBY./1H01 DIMENSION u113771.v¢13771.w113771.nx113771.RYt13771.RZ113771o13181 1.Rx1¢13771.Rv1¢13771.Rz1113771 COHHON UQVQUORXORYQRZQSXXQSYY0$ZZITXYoTXZoTYZORl0R20R30940R50960 1R70R8 GOVERNa.0000001 NUMBER-O nz1¢11=0.07 R21121-0.04 R21131c0.0 921(41--0.04 ~RZI(5)--Oo07 3x1111-0.0 RX1(2).°02 9x1131-0.3 ,RX1141-0.4 RX1£51.0045 9x116130.49 RXH7)=3.0 Rx1(8,)3300 RXlt9)85oO 9x11101s5. RX1(11)83.' 9x11121-9. Rx1<131s9. Rx11141~5. Rx1115139. N431 Mass 600 COflTINUE 00 11 L81.1377 u+U(L+891+U(L-89)+U(L+91)+U(L-91))+12.*(V(L+101 4+V1L-IO)‘V(L+B)‘V(L‘8)+U(L+82)+U(L*82)‘W(L+BO)“U(L*BOI)+20*(V(L+91 5)+V(L-91)+V(L+7I)+V(L*71)‘V(L+73)“V(L-73)_V(L+89)'V(L’891+UIL+91I 6+"(L-9I)+U(L+73)+U1L‘73)‘UIL+71)‘WCL‘71I‘U(L+89)‘U(L‘89III/RIE UILI=UILI+1025*RXIL) CONTINUE SURFACE PLUS X AXIS I‘N L‘LAKE(I¢J0K) CALLZS(L071989o-9I9‘73000000001 CALLYSIL0710899‘919‘73900090001 CALLXS(L9710390‘9IC‘TSQOQOQOQOI CONTINUE J3N2 SYN Z AXIS ON Y SURFACE I31 L=LAKEIIOJOKI CALLZSIL9730‘890000000000001 CALLXSIL9730*89009000000000) CALLYSIL0730-89000000000000) CONTINUE SURFACE PLUS Y AXIS DO 154 I329NI LzLAKEIIOJOK) CALLZSIL07I0739-890‘91000000001 CALLYS(L9710730-899~919000©OOOI CALLxstL.71.73.-s9.-91.0.0.0.01 CONTINUE , OUTSIDE EDGE z AXIS I3N LfiLAKEIIodoK) CALLZSCLQTIo‘QI009090900000) CALLYSIL0719’9I000090900000) CALLXSIL07I0-91000000000000) CONTINUE K31 J31 INSIDE X AXIS DO 334 1:29N1 L3LAKEII0J0K) CALLXSILogl089.090.090.000) CONTINUE 23 32 324 234 24 244 28 123 400 106 121 107 102 X AXIS AT THE SURFACE IBN L'LAKEII0J0K) CALLXSCL08900000000000000) DO 234 J=20N12 INSIDE Y AXIS 131 L=LAKEIIOJOKI CALLYSIL091073000000000000) CALLXS(L0910730000000000001 INSIDE 2 SURFACE DO 324 1320N1 L=LAKE(IOJOK) CALLYS(L07107308909100000000) CALLXSIL07I07308909100000000) CONTINUE SYN Y AXIS ON X SURFACE I‘N L=LAKECIoJ0K1 CALLYS(L07I0890000000000001 CALLXSIL071089000000000000’ CONTINUE J=N2 Y AXIS AT THE SURFACE I=1 L3LAKEII0J0K) CALLYSIL07300000000009000) CALLXSIL073000000000000001 SYM X AXIS ON Y SURFACE DO 244 I320N1 L'LAKEIIOJOK) CALLYSIL071073000000000000) CALLXSIL0710730000000000001 CONTINUE OUTSIDE Z EDGE AT SYM INTERSECTION I=N L8LAKE110J0K) CALLYSIL07100000000000000) CALLXSIL07100000000000000) CONTINUE CONTINUE RMAX=00 DO 110 I=10N DO 110 J=10N2 DO 110 K=ION3 L'LAKEII0J0K) IFCABSF(RMAX)“ABSF(RX(L))) 10601210121 RMAXzRXIL) 1F(ABSF(RMAX)“ABSF(RY(L))1 10701220122 RHAX‘RY(L) 103 I22 IFIABSFIRMAXI’ABSFIRZ(LI11 10801100110 108 RMAX‘RZILI 110 CONTINUE IF‘ABSF(RMAX)“GOVERN)12501250401 125 N4=N4+l NUMBER=0 NU‘Z IF(N4-61 70007000701 700 GO TO 116 701 N5=N5+1 N4=1 RZIIII'0015 R21(2I=0012 R2113I=0008 RZIC4I=0004 R2 I 15130 0 O RZII6)=‘0004 RZI(7)=‘0008 RZII813‘0012 RZI(9)=“0015 IFIN5-3170207020703 703 NU=1 702 GO TO 116 401 GO TO 115 116 PRINT 104 REUIND 5 WRITE TAPE 50U0V0W0NUMBER0N40N50RZI(I)0RZII210R211310R211410 IRZI(5)0RZI(6)0RZII710R211810R2119) REUIND 5 PRINT 1000NUMBER0RMAX PRINT 1030 N PRINT 1010N9N20N3 PRINT 5000P1 DO 120 K=I0N3 DO 120 I=IoN DO 120 J=I0N2 L‘LAKEIIOJOK) CALL STRESS (I0J0K0N0N20N31 120 PRINT 102 I0J0K0UCL10VIL)0U(LI0SXX0SYY05220TXY0TX20TYZ PAUSE 100 GO TO (504013OI0NU 504 CONTINUE STOP END SUBROUTINE STRESSIIOJOKONON20N31 DIMENSION U(1377)0V(I37710U(1377)0RX(137710RYII377)0R211377101313) COMMON U0V0U0RX0RY0RZOSXX0$YY0SZZ0TXY0TXZ0TYZ0RI0R20R30R40R50R60 1R70R8 AL=RB+R6—20 ALIRAL+10 45 11 12 13 104 ALMI‘IoVAL L=LAKEIIOJOK1 NI‘N-I ANI=N1 NIZRNZ‘I ANIZ’NIZ AN=N AN11=AN_105 C1=4000000/(AN1*R11 CA=CI*AL1 CBRCI*AL CC=C1*ALM1 L1=L+1 L2=L-I L3=L+9 L43L-9 L53L+81 L63L-81 LX=L1 LXX:L2 LY=L3 LYY=L4 L23L5 LZZ‘LO UC=10 VC=10 UC310 IFII-1120203 LXX‘L UC=2o IFIN’II40405 LX=L UC820 IFIJ‘1160607 LYY=L VC=2o IFIN2“J)80809 LY‘L VC320 IFIK-1110010011 LZZ8L UC=2o IFCNB-K) 12012013 LZ=L UC‘Z. CONTINUE UXX‘UC*(U(LX)’UILXX11 UYYBVC*(U(LY1“U4LYY11 UZZ‘UC*IU(LZI-UILZZII VXXIUC*(VILXI‘VILXX11 60 105 VYYBVC*(V(LY1‘V(LYY)1 V22=WC*(V(LZ)*V(LZZ)1 UXX=UC*(U(LX1-U(LXX)) UYYSVC*(U(LY)‘W(LYY)1 UZZ‘WC*(U(LZI‘W(LZZI) SXXSCA*UXX+CC*(VYY+UZZ) SYY=CA*VYY+CC*(UXX+UZZ) SZZ=CA*WZZ+CC*(UXX+VYY) TXY=CB*(UYY+VXX1 TXZ=CB*(UZZ+UXXI TYZSCB*(VZZ+WYY) RETURN END FUNCTION LAKE’IIOJQK) LAKE 3 I + J*9 + K*81 - 90 END SUBROUTINE XS (L0IBI01320183018401550186018701881 DIMENSION U(1377)0V(1377)0U1137710RXI1377)0RY1137710R211377101318) COMHDN U0V0U0RX0RY0RZ0SXX0SYY0SZZ0TXY0TXZ0TYZ0R10R20R30R40R50R60 1R70R8 IBIII=IBI 13(2)=IB2 IB(3)=I83 18(413184 13(513155 13(613136 18(71‘187 15(815188 XSUM=00 DO 9 M3103 AM=M-I IA=1 JA=9 KA=81 A1810 A2=10 A3310 IFIIB(M))20703 KA=~81 A3=-10 IFIIB(M)-KA)40705 JA=-9 A23-10 IFIIB(M)—KA-JA)60708 IA=-1 A1=~10 CONTINUE XSUNBI’RI’UIL)+R2*(U(L+KA)+U(L+JA11+R3*UIL+JA+KA1+R4*U(L+IA1+R5*I 1U1L+IA+JA)+U(L+KA+IA1)+R6*U(L+1A+JA+KA1+20*A1*A2*(’30*V(L)+R7*V1L 2+JA)-V