ABSTRACT APPLICATION OF NUMERICAL MAPPING TO THE MUSKHELISHVILI METHOD IN PLANE ELASTICITY by Vincent L. Plsacane A method is demonstrated for treating plane elasticity problems for simply connected regions. it consists of use of numerical mapping methods to map the simply connected region onto the unit circle in order to apply the Nuskhelishvili complex variable method. This approach now makes the complex variable method susceptible to auto- matic solution on a digital computer. It provides an advantage over the solution of the biharmonic equation by finite-difference approxi- mations, whose accuracy is limited by the fact that automatic iterative solutions of the resulting linear algebraic equations do not converge; the system of equations must then be solved by direct elimination methods, and this severely limits the number of equations a given computer can solve. Consequently, it limits the solution to a very coarse mesh. The numerical mapping, on the other hand, can be accomplished by an iteration procedure using a large number of mesh points. To demonstrate that this method is not difficult to use in practice and to make some estimate of the accuracy. an example is considered which has a large stress gradient and for which the exact solution was obtainable. The stress distribution is obtained by both the complex variable method with numerical mapping and the finite- dlfference method. A comparison of the results shows that the complex variable method furnishes better results than the finite-difference method and agrees remarkably well with the exact solution except on the boundary (where two of the three components of stress are known) and in the neighborhood of the corners. Conducting-paper techniques for the conformal mapping were also used in this study mainly as a check on the numerical procedure. However, the very good accuracy achieved by the electroanalogic method suggests that it may be used to advantage, especially when the region under consideration has curved boundaries. APPLICATION OF NUMERICAL MAPPING TO THE MUSKHELISHVILI METHOD IN PLANE ELASTICITY by Vincent L. Plsacane A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of MetallurQY. Mechanics, and Materials Science i962 ACKNOWLEDGMENTS The author wishes to express his sincere appreciation to Dr. Lawrence E. Malvern who suggested the problem and under whose lofty inspiration, constant supervision, and unfailing interest this study was made possible. Thanks are also due to the members of his guidance com- mittee: Dr. Charles 0. Harris, Dr. George E. Hase, Dr. Richard Schlegel, and Dr. Charles P. Veils, as well as to the Applied Mechanics faculty who have contributed to the academic background necessary for the preparation of this paper. Finally, to his wife June, to whom this paper is dedicated, the writer is profoundly grateful for her patience, understanding, and moral support. TABLE OF CONTENTS LIST OF TABLES . LIST OF FIGURES LIST OF APPENDICES . LIST OF SYMBOLS CHAPTER I INTRODUCTION . l.l General Discussion . l.2 Background . CHAPTER ll FUNDAMENTALS 2.l Introduction . 2.2 Muskhelishvili Complex Variable Method . 2.3 Theory of the Electrically Conducting Sheet 2.h Theory of Conformal Mapping by the Conducting Sheet . . . . . . . . 2.5 Theory of Conformal Mapping by Numerical Means . 2.6 Determination of LJ (CZ/{MTC7I by Numerical Means . . . . . . . . . 2.7 Determination of U( (C/L/(C) by the Conducting Sheet . . . . . . CHAPTER III EXAMPLE - CONFORMAL MAPPING OF A SQUARE REGION . . . . . . . . . . . . . 3.l Introduction . 3.2 Conformal Mapping by the Conducting Paper 3.2.l Equipment . . 3.2.2 Construction of Models . 3.2.3 Test Procedure . 3.2.h Calculations . iii Page mVNW Ih ‘7 2| 26 27 30 30 3| 3| 32 35 35 3.3 Conformal Mapping by Numerical Means . 3.1+ Evaluation of U (Gibb/(C) on the Boundary of the Square Region . . . . 3.5 Conclusion . CHAPTER IV EXAMPLE - STRESS SOLUTION . h.l Statement of the Problem . h.2 Exact Solution by Superposition Integral h.3 Solution by Complex Variable Method h.h Finite-Difference Solution . h.5 Conclusion . CHAPTER V SUMMARY AND DISCUSSION . BIBLIOGRAPHY . Page ho AS #8 so 50 SI 53 58 65 67 7o 3.3 Conformal Mapping by Numerical Means . c7 /—-—, 3.4 Evaluation of («II )/LJ(CT) on the Boundary of the Square Region . . 3.5 Conclusion . CHAPTER IV EXAMPLE - STRESS SOLUTION . h.l Statement of the Problem . 4.2 Exact Solution by Superposition Integral h.3 Solution by Complex Variable Method h.h Finite-Difference Solution . h.5 Conclusion . CHAPTER V SUMMARY AND DISCUSSION . BIBLIOGRAPHY . Page ho as us 50 50 SI 53 58 65 67 7O I 3.3 3.h-i h.3-I h.h-l LIST OF TABLES Values of a and B at Grid Point of the Square Region from the Numerical Mapping Method . Fourier Coefficients for UICI/L/(C) of Square Region from Numerical and Conducting Paper Methods Fourier Coefficients for fl + if2. Traction on Square Region . Finite-Difference Solution to Biharmonic Equation for Square Region Page “3 57 6O .2-I .2-2 .2-3 .Z-A wuww H.3-l ALI-kl “-2 4.4-3 m-u tell-s A.I-I A.2-I A.2-2 LIST OF FIGURES Conducting Sheet . Geometric Representation of Conformal Mapping by Conducting Sheet Geometric Representation of Conformal Mapping by Numerical Method First Stage Model - Square Region Second Stage Model - Square Region . Analog Field Plotter with Pentagraph . Curvi linear Net of CI. ,B Conducting Paper Method Log r Versus (I -(3) - to Determine Constants A, C . Real and Imaginary Parts of MIC) T’IG‘) Versus 6 Boundary Traction for Square Region Finite-Difference Approximation at Boundary Stress Distribution in Square Region, Stress Distribution in Square Region, Stress Distribution in Square Region, Stress Distribution in Square Region, Schwarz-Christoffel Mapping - Polygonal Region . Schwarz-Christoffel Mapping - Square Region Schwarz-Christoffel Mapping - Approximating Square Region vi for Square Region by Page IS 22 33 33 34 36 39 ‘47 56 59 6| 62 63 61+ 7t. 75 77 LIST OF APPENDICES APPENDIX A CONFORMAL MAPPING OF POLYGONAL REGIONS . A.l Theory of Schwarz-Christoffel Transformation A.2 Conformal Mapping of a Square Region A.3 Proof That All Schwarz-Christoffel Interior Mappings Satisfy the Smirnoff Condition . APPENDIX B FINITE-DIFFERENCE SOLUTIONS TO PARTIAL DIFFERENTIAL EQUATIONS . B.I Discussion 8.2 Laplace and Biharmonic Equations Page 73 73 7h 78 82 82 83 9(2) G(I)EEu + iv LIST OF SYMBOLS Inner radius of circular ring in g -plane. Complex coefficients of power series for ¢(€) and W(C). Real constant for mapping function of unit circle, Eq. (2.h-8). Complex constants in Schwarz-Chrlstoffei mapping formula. Arbitrary constants of ¢(z), real and imaginary respectively. Coefficients of complex Fourier series for “(aye/(Ci ' Radius of inner boundary of region In z-plane. Length associated with distributed load on square region. Constant related to inner boundary of region in z-plane. Subscripts meaning exterior, interior, and boundary; used for calculation of stress function outside boundary of region under consideration. integrals of boundary tractions. Function which maps unit circle into C( , [3 -plane. Function which maps region in z-plane into C: ,1/3 -plane. Holomorphic function in conducting sheet, Eq. (2-3-5)- Hoiomorphic function for numerical mapping, Eq. (2.5-3). Grid Interval. Imaginary unit, \l-I. vlli Current density in conducting sheet. Total current function or flux function in conducting sheet. Integrals associated with Smirnoff condition, Appendix A.3. Total current flow In circuit, conducting sheet. Coefficients of complex power series for g I=(,,/"(z). Green's function of the first kind. Regular harmonic function - part of Green's function. Inner boundaries of regions in g and z-pianes respectively. Outer boundaries of regions in 4' and z-plenes respectively. Outward normal to boundary of region in z-plane. Constants related to coefficients of linear simul- taneous algebraic equations. Intensity of distributed load on portion of square region. Polar coordinates in z-plane. Resistance in ohms. Total resistance of circuit in ohms. Functions associated with Mid/m , Eqs. (2.6-3). Arc length on boundary contour of region in z-plane. Thickness of conducting sheet. Airy's stress function. Function satisfying Laplace's equation. Electrical potential in conducting sheet. Width of path for current flow in conducting sheet. Cartesian coordinates. P o 0‘ Cx' Cy' Txy ax’ay’Txy dim. (big) X(z) Wk). i/Né) wig) Cartesian coordinates. Boundary tractions. Complex variable In plane of arbitrary region. Angle in €-piane such that 0- 6k =a. Non-dimensionaiized variable measuring i. Non-dimensionallzed variable measuring V(x,y). Real constants appearing in Schwarz-Chrlstoffel formula. Angle in C-plane, such that 7= 3a. Arbitrary complex constant of UV(1)- Arbitrarily small angle associated with 6 in g -piane. Arbitrarily small length, so that p= l- E . Complex variable in plane of unit circle or annular region. Angle between line tangent to boundary and x-axls. Cartesian coordinates in g -plane. Polar coordinates in 4 -plane. Sheet resistivity - resistance of any square element of conducting sheet. Resistivity in ohms per unit length of a conductor with unit cross-sectional area. Value of C on boundary of unit circle, C‘= ei9 . Normal and shear stresses referred to x,y. Normal and shear stresses referred to §,9. Holomorphlc functions for representation of stresses. Holomorphlc function of the complex variable 2. Holomorphlc functions for representation of stresses. Mapping function. CHAPTER I lNTRODUCTION I.l General Discussion In recent years considerable progress has been made in solving plane elasticity problems approximately by methods such as the finite- difference method, utilizing high-speed digital computers, and by techniques involving photoelastic and electrical analogies. However, these solutions frequently lack generality or are deficient either in accuracy or convenience. The finite-difference method, in particular, leads to a system of linear equations not solvable by ordinary itera- tive schemes; the solution of the system by direct elimination requires either a very large computer or a very coarse mesh. This paper indi- cates how the Muskhelishvlli complex variable method with approximate conformal mapping transformations can be applied as a general method of analysis with a minimum of effort to give accurate results for simply connected regions. The complex variable method involves the solution of a boundary value problem by mapping the given simply connected region conformally onto the unit circle. The problem of the exact construction of a function effecting this mapping very often presents great or even insurmountable difficulties. Runge and Walsh (see, for example, Kantorovich and Krylov, I958) have demonstrated that if the boundary 0f the simply connected region is a simple curve, i.e., a curve Without multiple points, then a polynomial can be chosen so as to differ arbitrarily little from the transforming function everywhere in the region, Including even the boundary. However, even though the means for representing the transformation is known, Its structure will usually be so complex that use of it entails a great expenditure of labor. it is sufficient to recall the problem of determining the constants in the Schwarz-Christoffel formula for the function trans- forming a circle into a polygon to be convinced of this. Further. truncation of the polynomial after a reasonable number of terms may admit prohibitive errors, as is demonstrated by truncating the series for the square in Appendix A. Therefore, for the approximate solution of problems, the simpler method of effecting the transfenmation numerically or through electroanalogic methods ls of very great value. it permits a method of solution of elasticity problems suitable for use with a digital computer of moderate capacity, and such that each specific problem requires only easily accessible information an the geometry of the region and the applied traction. The present investigation applies these simple mapping techniques to a particular case -- a square region. We do not make use of the symmetry properties of this example, since we wanted to test the effort involved and efficacy of the method for a general region not possessing such symmetry. He then compare the solution based on these mappings with other methods for a particular loading of the square, selected because the exact solution is known and because It involves a region in which the stress varies rapldIy'wIth distance and therefore provides a critical test 0f the ability of an approximate solution to reproduce details of the stress distribution. We find remarkable agreement with the exact solution except, as expected, on the boundary itself (where we already know two of the three components of stress) and near the corners of the square. The solution is found to be considerably more accurate than that attainable with the same computer facilities by a finite- difference approximation to the biharmonic equation. The reason for this is that in the iteration method of numerical mapping we are able to use a grid of 8hl points while for the finite-difference solution by direct elimination we were limited to 36 points. The method presented here should be of great interest to stress analysts, since it provides a straightforward solution to plane elasticity problems for simply connected regions with remark- able accuracy. Before considering the mapplng example and the com- parison problem, we first present a brief review of plane elasticity and then summarize the fundamentals of the complex variable method, the conducting sheet, conformal mapping by the conducting sheet, and numerical conformal mapping. l.2 Background The solution to two-dimensional problems in elasticity reduces to the integration of the differential equations of equilibrium together with the compatibility equations and the boundary conditions. 111a British astronomer, G. B. Airy, in l862, was the first to note that the equilibrium equations imply the existence of a scalar function, now called the Airy stress function, such that the stresses are given by second derivatives of this function; the compatibility e(Nations demand that this function satisfy the biharmonlc equation. By using strain energy methods or taking for the stress function polynomials of various degrees or trigonometric series, a number of important problems of simple geometrical shape can be solved (see Timoshenko and Goodier, l95l). For more complex shapes the difficulties of obtaining analytic solutions become formidable, but these difficulties can be avoided somewhat by experimental methods such as the membrane, electrical, and photoelastlc methods. Brief reviews of the fundamentals of these methods are given by Mindlin (l939) and Higgins (l956). The membrane and electrical methods are based on the fact that Laplace's equation governs the variation of the sum of the principal stresses. These methods are applicable only to regions along whose boundaries the principal stresses can be previously determined. The most useful experimental procedure, the photoelastlc method, is based on the principle that when transparent materials are stressed, their optical properties undergo changes, which can be measured and related quanti- tatively to the state of stress. The accuracy attainable and the difficulty of simulating the applied loads limit the application of these methods. With the arrival of high-Speed and large-capacity computing rnachines, it became possible to find solutions to plane elasticity Problems by finite-difference methods. However, the system of linear Slinuitaneous equations obtained by the difference approximation to the biharmonlc equation cannot be solved by an iterative scheme, since tTiis technique produces divergent approximations. The system of linear 6questions must be solved by an exact means, such as triangulation of the augnented matrix, which Iimi ts drastically the number of variables (or consequently grid points) we can use. The MISTIC digital computer at Michigan State University, for instance, can solve a maximum of 39 linear simultaneous equations by the triangulation method; this results In a very coarse grid. The application of complex variable theory in plane elasticity was proposed by Kolossoff in I909, and later developed by Muskhelishvlli and other Russian investigators. A complete account of this work, including several different methods of applying the complex variable theory, is given by Muskhelishvlli (I953). The method used in the present work involves the mapping of the given region onto the unit circle; it provides a solution to the biharmonlc equation by reduction of the boundary value problem to the solution of certain functional equations in a complex domain. Savln (I956) (English translation now available, Pergamon Press, l96l) gives a summary of applications of this method to many important examples of stress concentration in the Infinite region exterior to a hole; in the cases reported the infinite series representing the exterior maps converge rapidly and may be approximated adequately by polynomials of low degree. For simply connected interior regions, however, it turns out that series repre- senting the mapping functions converge very slowly, so that solutions have been limited to a few regions whose exact mapping function is known. Thus, while an extremely powerful and general theory has been developed, it has been severely limited in its application to simply connected regions by the necessity of knowing the mapping function exactly. Parkus (I960) considers the application of the complex variable method to simply connected regions by approximate conformal mapping, using an electric tank to effect conformal transformations as in the work of Bradfleld, Hooker, and Southwell (I937). This combination of approximate mapping together with the complex variable method produces a powerful approach to the solution of plane elasticity Problems. The conducting sheet methods discussed In the present work are based on the same principles. However, questions arise as to the accuracy of these analog mapping methods. In the present investigation, therefore, a presumably more accurate numerical mapping technique has also been used and the results of the two mapping techniques have been compared with each other and the stress solution by numerical mapping has been compared with the known solution of a particular problem in plane elasticity. This concludes a brief review of plane elasticity. in the next chapter we shall begin with a summary of the Muskhelishvlli method of using complex function theory In elasticity and then present the fundamentals of the mapping techniques to be used. CHAPTER II FUNDAMENTALS 2.l Introduction This chapter concerns itself with the fundamentals on which the techniques of analysis proposed In the introduction are based. We will demonstrate that the solution to the boundary-value problem for plane stress or plane strain can be expressed in terms of two analytic functions of a complex variable. The components of stress can be obtained from these functions. Required in the solution of the boundary—value problem and the stress calculation is the conformal transformation of the given region onto the unit circle and the evaluation of LJ(€}{:%;5 e is the value where z .L‘AC) Is the mapping function and 6‘= el of C. on the unit circle. We will show that the point-to-polnt mapping can be achieved by numerical means or through use of the electricall conductln sheet. The function (vs/(7 ! ma also Y 9 um Y be obtained by the same processes. When this mapping and the calcula- tion of U%l(6‘) are available, the Muskhelishvlli complex variable method can be readily applied to any plane elasticity problem in the mapped region. 2.2 Muskhelishvili Complex Variable Method The stress components in plane strain, plane stress, or generalized plane stress problems can be expressed in terms of Airy's stress function U(x,y). (See, for example, Sokoinikoff, I956.) If body forces are absent, or are accounted for in a supple- mentary solution to be superposed, we have 2 2 2 c) u c) u c) u 6. a: . a” , =3 - 2.2-I x a y2 Cy _ (3 x2 Txy ()xay ( ) where U satisfies the biharmonlc equation vhu=o (2.2-2) in the region and the appropriate boundary conditions. When the traction components Xn(s) and Yn(s) are given in terms of the arc length 5 on the boundary, the boundary conditions on U assume the form 5 .52).”. z: -j Yn(5) ds 4» const. 3 HIS). Ox 50 (2.2-3) 5 3—9.. I xn(s) ds + const. E f2(5). Y s 0 where the constants are of little consequence since the stresses are determined by the second derivatives of U. The solution of the biharmonlc boundary-value problem is especially convenient by complex function theory. In the method of complex functions, the mathematical formula- tions of the boundaryovalue problems of plane strain and plane stress are identical; their solutions depend on the determination of two functions of a complex variable from certain functional equations. We shall proceed to derive these equations for a simply connected region bounded by a contour L2. For a more detailed development, reference is made to Muskhelishvili (I953). Goursat (see, for example, Muskhelishvili, l953) in l898 Imas the first to show that every real biharmonic function U(x,y) may be represented as zu - 2a, [291sz mm] = 295(2) + 255??) + Xiz) + 327;). (2.2-u) where¢(z) and X (z) are any two holomorphic functions of the complex variable 2 - x+iy. The boundary condition on U, Eqs. (2.2-3), now takes the form (95(2) + nghi + (Hulk +(—-y- + i 3” YIL2= fl+ if2, (2.2-5) where (U(z) '5 3% and L2 refers to the boundary contour. This boundary condition on the simply connected region suffices to completely deter- mine the two holomorphic functions ¢(z) and Wk) In the region apart from the arbitrary additive terms biz +7 and T’respectively, where b is a real constant and 7 , T’are complex constants. Eq. (2.2-5) completely determines the state ofstress of the body, since it is the second'derivatives of U which specify the components of stress. To remove all arbitrariness from the functions¢, 50 we require that gbm) - o, j[¢'(o)] .. 0, W0) a o, (2.2-6) “h." the stresses are given. The two-dimensional elastostatic problem ID with given boundary tractions then takes the form of determining two holomorphic functionsi¢b(z) and 97(2) In the given region subject to the boundary condition of Eq. (2.2-5) and the requirements of Eq. (2.2-6). As a first step in determining these functions, we consider the region bounded by L2 to be mapped conformally onto the circle |§|<| by the holomorphic function 2 = (,J(§). (2.2-7) The existence of such a holomorphic one-to-one transformation is assured by the Riemann mapping theorem (see, for example, Nehari , > >3 Ia I952). in the “a ~plane we let 5 pe where 9,9 are polar coordinates; on the boundary I4) =- I of the unit circle, we write 9"“ eieE G . Substitution of Eq. (2.2-7) i‘nto Eq. (2.2-5) results in the transformed boundary condition ¢(0)+ Z/L'—-:-:;—— ¢I+CVI SUI?) = f, + Ifz , (2.2-8) where we have, for simplicity, used the notation ¢(C) for ¢[U(C)] and yflc) for SUE/(£5 , and where f‘ + if2 has at each point of the circular boundary in the €-plane the same value that It has at the corresponding Image point of the boundary L2 in the z-plane, this value now being known as a function of 0“. Since ¢(z) and {/l(z) are holomorphic in the region bounded by L2, the functions¢(€) and Stag) will be holomorphic in the circular region and the problem has now been reduced to finding two functions $(4) and V(g) holomorphic II in the interior of the unit circle and satisfying the functional equation (Eq. 2.2-8) on the boundary. This equation is meaningful only if U(%.(C‘) is defined everywhere on the boundary; for this it is evidently necessary that the boundary value U(C‘) differ from zero everywhere. Conformallty only assures this in the interior, and there is some question about it at boundary points, especially when there are corners. Smirnoff (I932) has proved a sufficient condition that the boundary value approached by i,J’(C‘) is given by I../(§)| 5% , and U(G‘) ¢= 0 almost everywhere; (2.2-9a) €=°‘ the sufficient condition is that 277 i . ] IUIPPLBH d 9 remains bounded as (3 —--I. (2.2-9b) o This condition will certainly be satisfied if the boundary contour has Icontinuous derivatives up to the second order, i.e., the curvature changes continuously, and in this case IUIOLQII will actually iwave a positive lower bound, as is shown by Smirnoff (I932). For the square we will study later, it turns out that L/(Q—boo at a corner and U(Cym --> 0 , but it can happen that U’(€)—-p 0 b iit a corner, e.g., at a re-entrant angle of a polygon, according to time Schwarz-Christoffel formula, even though all Schwarz-Christoffel mappings satisfy the Smirnoff condition (2.2-9b). (See Appendix A.I). Since ¢IO and V(é) are holomorphic in the interior of the unit circle, they may be represented in power series l2 95¢) '00:,“ Ck . W(§)-:a'k g“ (2.2-Ia) which becanes 00 k 90 k 9D (0") '2 0.,0‘ . WWI-Z aLO‘ (2.2-ll) o - a when g - 0‘. The convergence properties of Eqs. (2.2-IO) on the circle boundary will be discussed when we present a method for solving for the coefficients 'k' a‘k. We assume the following expansions are possible: U(Gi— . + s cm". (2.2-I2) /L/IG‘) “ -oo and -+oo k f, + Ifz - Z Ak 0‘ . (2.2-l3) -oo The assumption. that (e/ki}<:fifi;) and f‘ + ifz are expandable L. in a Fourier series, lmposes some continuity restrictions both on the boundary shape and traction. it is sufficient (see, for example, Muskhelishvlli, l953) If the functions satisfy the Dirichlet conditions; then the Fourier series converge at all points on the Interval. This restriction does not practically limit our analysis. permitting, for exuuple, finite imp discontinuities In the boundary traction. Substituting Eqs. (2.2-ll), (2.l-l2), and (2.2-l3) into Eq. (2.2-B), we obtain k f E: kl k Eako‘ + E bxG kakO‘ + akO‘ / ’ ’ (2.2-III) Multiplying out the series and comparing coefficients of (Tm (m ‘3 I.2.3.?"), one finds M8 a... +' um. um”. = A... (m = I.2.3.---); (2.2-Is) k==l similarly, one obtains for 6"” (m = 0,l,2,'~-) o-o a,',, + k2 kak bwk- = A_,,, (m = O,l,2,---). (2.2-i6) II if the infinite system of Eqs. (2.2-l5) can be solved for ak (k = l,2,3---), the aL (k = O,l,2,---) are determined individually from Eqs. (2.2-l6). To completely determine the functions Qbhg) and V](§) we require, similarly to Eqs. (2.2-6), that mm) = 0, 925(0) = a0 .. o ,.i(a,) = o. (2.2-I7) Muskhelishvili (l953) shows the series for¢ (g) and SU(€) obtained In this way will satisfy the conditions of the problem provided the given functions f‘ and f2 are sufficiently regular, i.e., for example, their second derivatives with respect toe , satisfy the Dirichlet condition and the first derivatives exist everywhere (this requires continuous traction). With dug) and l/f(€) then evaluated by l4 Eq. (2.2-l0), the stresses can be obtained immediately from the relations I-I Re [¢I(z)] , G‘x +i3‘y = (2.2-l7) C‘y - cx + 2i Txy = 2 [3. Cp/(z) + W’(Z)] . or C§(+-C§ = h Re [iii g :L (2.2-l8) N G‘y-G‘x+2iTxy= %%(%)+% ; obtainable from Eq. (2.2-5). The solution of the plane elasticity problem can therefore be reduced to the finding of the holomorphic function (,j(§) mapping the unit circle of the C-plane onto the elastic body region of the z-plane. 2.3 Theory of the Electrically Conducting Sheet To effect by approximate means the conformal transformation required In the solution of the elasticity problem, we first consider the fundamentals of the conducting sheet and demonstrate the existence of a complex potential function, holomorphic in a simply connected region of the sheet. See, for example, Bradfleld, Hooker, and Southwell (I937). In Fig. (2.3-l) a sheet of conducting material of thickness t (either solid or liquid) is shown, whose boundary is subject to a given lS voltage distribution V(x,y). Consider an infinitesimal element Y i . ,V‘”) ”“33de [:de _.. dv ..... i dx x ix dx ix + 3: dx / '._./ (a) (b) Fig. 2.3-i. Conducting Sheet (dx dy t) of the sheet; the current leaving across its boundaries must equal the current entering across its boundaries (see Fig. 2.3-lb). Thus, the current balance equation (equation of continuity) is O‘X +9Ji=o , (2.3-l) c)”. dY where i and i are the current flows per unit length. The resistance x y of a finite element of the conducting sheet to electrical flow is R = (ZSL where R is the resistance in ohms, L the length of the element in the direction of the current flow, t the thickness, w the width, and (30 the resistivity in ohms per unit lenth of a conductor with unit cross-sectional area. if 0! 1.32 (2.3-2) defines the “sheet resistivity,” the resistance of a square element of 16 I the sheet (L = W) is then R = (3 . Since the current flow across each boundary is proportional to the voltage gradient at that boundary, we write X .17 Qv, iy=-_', 3.3.1, (2.3-3) (3 ox (3 OY where V(x,y) is the electrical potential and the negative sign is introduced because the current flows in the direction of decreasing potential gradient. Substituting Eqs. (2.3-3) into Eq. (2.3-l) and . I assuming (3 to be a constant, we obtain _:.)?_2_.\2.’.+.C.E))i‘2£=V2V—'=O, (2.3-ll) x v demonstrating that V(x,y) is a harmonic function. it follows from the theory of harmonic functions (see, for example, Churchill, l9h8) that Laplace's Equation (2.3-h) in a simply connected region implies the existence of a conjugate harmonic function cfi(x,y), unique except for an arbitrary additive constant, such that g( z) = 0’ + w (2.3-s) is holomorphic in the simply connected region under consideration. The contour curves V(x,yl = constant of the electrical lfibtential V(x,y) are called equipotential lines. A physical interpre- :tation of l(x,y) can be given by analogy to the two-dimensional motion ilf an incompressible fluid (see, for example, Churchill, l9h8). It follows from the Cauchy-Riemann equations and Eqs. (2.3-3) that l7 lx=+%-:7,iy=--—:-‘-, (2.3-6) where l(x,y) plays a role analogous to the stream function of the hydrodynamic theory and is designated here as either the total current function or flux function. Contour curves l(x,y) = constant will be called streamlines because of their similarity to streamlines in fluid flow. it follows directly from Eqs. (2.3-6) that the total current flow across an equipotential curve between any two points is given by the difference between the values of the flux function at the two points. That the holomorphic function g(z) has a definite value for «every 2 and that its existence may be detected by phySical means is the basis of the application of the conducting sheet to conformal transformations. 2.4 Theory of Conformal Mapping by the Conducting Sheet The theory of the conducting sheet will now be applied to effect the conformal transformation of an arbitrari iy-shaped doubly- connected region onto the interior of a ring bounded by two concentric (:iches. The conformal mapping of a simply connected region onto the uni t circle can be approximated by the mapping of a doubly connected region with a small circular inner boundary, since the equipotential lllmes near the center approximate circles. This follows from the '“lnerical mapping procedure we consider in the next section. Let us consider the three planes in Fig. (2.h-l). Fig. (2.h-la) shows an arbitrary doubly connected region in the z-plane, while l8 Fig. (2.h—lc) shows a ring bounded by concentric circles with inner radius a, in the C-plane. Suppose that the electrically conducting sheet occupies the regions between the closed electrodes L1, L2 in in the z-plane, and £1 ,,92 in the €'-plane.' Choose the electrical potential so that the voltage is zero at Li and also at,!', while the voltage is G at L2 and at I 2. Choose R for the electrical resistance of both circuits, i.e., between L‘, L2 and between lfi', fl 2. (b) Fig. 2.h-l. Geometric Representation of Conformal Mapping by Conducting Sheet ‘9 Then the total flow of current in each circuit will be i, where R l . v. if we divide both sides of Eq. (2.3-5) by \7. there results for the region in the z-plane Fix) 5 9(2) = 1' + i ¥; (2.1M) V Rl V or FM =C< + Ifl. (2.14-2) where C<'-§--;-. [32%. (2.4-3) so that (I and [B are non-dimensional quantities measuring l(x,y) and V(x,y) at a point whose coordinates are (x,y). By definition/2 wi ll have values ranging between zero on the boundary L2 and unity on the boundary L‘, while a must be cyclic; i.e., it must be multiple- velued in such a way that it increases by a fixed amount, the cyclic constant, for each circuit around the inner boundary. The difference between two values of l(x,y) along an equipotential line is the amount Of current flowing across the equipotential line between the two Points; therefore the cyclic constant of l is l, and the cyclic I constant of a Is .9... Similarly, for the annular region in the R I C-plane, if we select .53.. to be the same as in the z-plane, it R follows that f(€)= C(+l [2. (2.4.1,) 20 F(z) and f( g) are holomorphic in their respective cut regions, cut by C(-C( o , by the theory of the conducting sheet. The inverse Fol is holomorphic in thea , [3 region so that (A Q) .=. F-i [Hg )] (2.l+-S) is holomorphic since a holomorphic function of a holomorphic function is holomorphic. Thus, a point-to-point correspondence between regions in the z-and g-pianes, such that the values of C( and B at the corresponding points are the same, will be a conformal transformation. By exploring a sufficient number of contours of constant C( and a sufficient number of contours of constant [3 in the regions in the z-and é-planes, a curvilinear net is constructed by which the whole of the arbitrary region in the z-piane will be mapped onto the whole of the region of the circular ring. The contours of Ct and )8 in the z-plane must be found by an approximate method, but the annular region admits an exact analytical solution (see, for example, Streeter, l9’+8) of the form 'f(§)=O(+F/3 =-iAL09€, (2,1.-6) Where A is an arbitrary constant dependent on the size of the region. Sllmce g I F)°H3 C(+i,0 =A9-IALogp; (2.l+-7) I it follows, since the cyclic constant of (x is .9: , that A a (2.h-8) pl 2m ’ 2l andfor fl=lon£' a e {27’ "ET , (2.l-#9) determining the inner radius a. The contours for C( and [3 in the annular region can now be found analytically from Eqs. (2.l-b7) and (Lu-8). The conformal mapping of the arbitrary region onto the unit circle can therefore be reduced to the finding of the curvilinear net of constant a and constant B . While the physical conducting sheet presents a simplified experiment to effect the conformal trans- formation, a question arises on the accuracy that is achieved, especially when the mapping of a simply connected region is approxi- mated by that of a doubly connected region, using a small circular electrode for the inner boundary. Greater accuracy would be expected from the numerical mapping procedure described in the next section. 2.5 Theory of Conformal Happj ng by Numerical Means We now present a method whereby a simply connected region can be mapped conformally onto the unit circle by numerical means instead of by the experimental method of the preceding section. This has the advantage that the entire stress problem can be evaluated by numerical methods as well as the fact that we should expect an increase in the accuracy of the conformal mapping. Consider Fig. (2.5-l), similar to Fig. (2.h-l), where the inner boundaries of the regions in the z-and €~planes are each con- tracted to a point, and the curvilinear net of C( and B is 22 retained, where C: and B have the same definitions as in Section 2.14. z -(t/(§)_ “r y W Fig. 2.5-l. Geometric Representation of Conformal Transformation by Numerical Methods The analytic solution for the (X .[j is Eq. (2.4-6): net of the circular region C1+i}8 =-i’ALog(: , (2.5-i) With the boundary condition [3 == 0 on I 2 and where the inner boundary condition [3 = l is replaced by [3 = 00, at g = 0. A is an arbitrary constant which for simplicity is selected to be “hity. For the simply connected region, the (I . fl region in Fig. (2.5-lb) is the strip extending to infinity in the fl-direction. 23 by the Riemann napping theorem, see Nehari 0952), a holomorphic function ge (.j" (1) exists which maps the region inside I.2 conformally onto the unit circle in the g-piane in a one-to-one manner and such that “4(0) - O. This function can therefore be represented bY a series of the for. 4- ER" 2" in any circle i centered at z - 0 and lying entirely inside L2, where it. it 0 and (k' + k2 z. '0- k, 22 +--°) ¢ 0 in the circle because otherwise a second point would nap into g- 0 violating the one-to-one mapping assuaptlon. This last series also converges in the circle (its coefficients are dominated by those of the series expansion of the derivative dU’lldz) and is therefore bounded. Hence, in any circle centered at z - 0 and lying inside L: a + ifl - - i Log [10" + k2: + R32? +...)] (2.5-2) - - i Log 2 + 6(2) (2.5-3) share 6(2) is a holomorphic function. in the doubly connected region lying inside l.2 and outside the circle, C( + i [3 must be analytic and have the same cyclic properties as - i Log 2. Hence the representa- tion of Eq. (2.5-3) is valid throughout the region inside L2, mere 9(2) is a holmorphic function throughout the region, defined by con- tinuation of the function represented by -i Log (k, + k2 z + k, :2 4"") inside the circle. 6(1) can be represented as 6(1) - u(x,y) + l V(x,y). (2.5-h) 21+ The imaginary part of Eq. (2.5-3) Yields for z = reIQ fl 3 - L09 l' + V(xoY)- (2'5-5) Since 6(2) is holomorphic, V(x,y) is harmonic, v2 v 0, and the boundary condition v Log V x2+y2 on L2 (2.5-6) follows from Eq. (2.5-5) since fl = 0 on L2. Laplace's equation for V(x,y) with the boundary condition Eq. (2.5-6) present a Dirichlet problem which can be solved by various numerical schemes, see Appendix B. With v determined, we can evaluate ‘/3 throughout the arbitrary region in the z-plane from Eq. (2.5-S). The real part of Eq. (2.5-3) than yields (X = G) +U(X.yl (2.5-7) where G) is the polar angle in the z-plane and u(x,y) is the complex conjugate of v(x,y), determined to within an arbitrary constant by the Cauchy-Riemann equations once v has been determined. The variable CI has the limits (10$ C( $ C(o‘l'27T. We can now show the relation between the Green's function for C simply connected region and the conformal mapping of that region on the circle. lie define the Green's function K(z,zo) of a plane region in the z-plane with smooth boundary by the following three require- ments (see, for example. Bergnan and Schiffer, l953): 25 (a) K(z,zo) = Log I + k(z,zo) is a regular harmonic function of z in the region except at z 2 20. K(z,zo) is real valued. (b) Near z = 20 the function k(z,zo) is regular harmonic in z. (c) K(z,zo) is continuous in the closed region in the z-plane except at z = 20, and has the boundary value zero. it follows from Eqs. (2.5-5) and (2.5-6) that [I (the non- dimensionalized voltage in the conducting sheet) plays the role of a Green's function with 20 = 0. We write from Eqs.(2.2-7) and (2.5-i), with A =1, that i; a U-i(z) =e iioz+ i fl). (2.5-8) (:1 is the conjugate of _/3 ; this function maps the region in the z-plane in a one-to-one conformal manner onto the interior of the unit circle in the C -plane. For additional information on the relation between the Green's function and conformal mapping, see Kellogg, U929) . We see, then, that the problem of determining K(z,0 ) is equivalent to determining a conformal map onto the unit circle such that U(O) = 0. Thus when C1 and [2 have been determined as functions of X and y, the point-to-polnt correspondence with known points in the C -piane having the same values of CI and fl is a conformal 26 mapping of the arbitrary simply connected region onto the interior of the unit circle. This conformal transformation effected by either analogic or numerical methods and the complex variable solution to the boundary-value problem combine to form the method of analysis for plane elasticity problems whose potentialities are explored in this investigation. We will also need the boundary values 0f bjtzpéJKCT) and we now show two different ways to evaluate these boundary values. 2.6 Determination of ‘*/0:(/1;/1C7) by Numerical Means In the solution of the boundary-value problem of Section 2.2 for the complex functions @ . y] - '5'“ function Umfl/(L/(G‘) must be evaluated and expressed in terms of a Fourier series. This can be accomplished from the results of the numerical mapping as follows: if the curvature of the boundary contour changes continuously, or if U(G‘) satisfies the weaker condition, Eq. (2.2-9b) of Section 2.2, then (Jule) = - i 6‘6 2%). . (2.6-i) Since e'le = Cos e - i Sin 9 , there results fil‘ 3+3? [(xR + yS) + i (yR - xS) ]L2 (2.6-2) where R s (file—jfzcdse -‘g.’é.)12 Sine , (2.6-3) M II _ (%%_)!2Cos 9 +615)!( 2 Sin 6 . 27 and L2, f 2 refer to the boundaries of the region in the z-plane and of the unit circle respectively. Eqs. (2.6-2) and (2.6-3) can be used to calculate LACK/(6‘) at points on the boundary of the unit circle. The boundary values of U(G%/i(6.) may alternatively be obtained from measurements on the conducting sheet as will be shown in the next section. “T" 2.7 Determination of LJwC7)/tj(67i by the Conducting Sheet Here we will demonstrate in a manner similar to Parkus (l960) how the conducting sheet can be utilized to evaluate the function Lj(ci}/LJXCT) . This will provide a check for the numerical process discussed in the previous section. Consider the doubly connected region, as in Fig. (2.h-l), where CI and /C? represent current density and voltage distribution respectively. Then I , \/0 g—g-icui/jiw ui§i[%§+i£5——;] .. ' c) ~ “‘§’[5§*'%9l° (2.7-1) 3 3 Now along the boundary L2 where [J 0. 3Q+ifl2--_gflncds,l +i_%%5in,l, (2.7-2) clY c)“ where ,l is the angle between the tangent line and the x-axis and n is the outward normal. Eqs. (2.7-l), (2.7-2) yield f' %3- “1+ ifl) =+ U‘(€) glgi [Sink + i Cos/l ,] (2.7-3) 28 Oi’ L(o(+i3) =+ U’Uisifii Mg?) (2.7—b.) <1; / 5 ,3 n Similarly, at the boundary of the unit circle I2 of Fig. (2.li-l) 01+ i/j) — +i e"9(%€_)£. (2.7-5) 2 9.... “C Eq. (2.h-6) yields upon differentiation lg =8 -(_A...)I = -A , (2.7-6) do P 2 so that Eq. (2.7-5) becomes __ “1+ i/3) = - i e“9 A. (2.7-7) From Eqs. (2.7-h) and (2.7-7), we obtain Lj'(6‘) = .253. ei (A -9-%). (2.7-8) (6 ,i it follows from “(5‘): r eig and Eq. (2.7-8) that --—L=’“°‘ ”ii-ii: i(@+/l-9-“/2.). (z_7-9) To evaluate this equation at a point on the boundary for some value of a . which determines 6 , we need only measure r, @ , A , and '08 d [1 ; the constant A is determined by the method discussed in Section 2.“. The only measurement offering any difficulty at all is that of I O8 3 n , which can be accomplished with relatively good 29 accuracy by graphical differentiation of curves of measured voltage along interior normals. See Section 3.“ for a comparison of values obtained in this way with those obtained from the numerical mapping for the square region used as an example to illustrate the methods devel0ped in this chapter. CHAPTER Ill EXAMPLE - CONFORHAL MAPPING OF A SQUARE REGION 3.1 Introduction We will now show by example that the methods of conformal mapping by numerical and electroanalogic means are not difficult to use in practice, and we will also make some estimate of the accuracy of our methods. Let us take for the example a square region of dimensions two units by two units; this enables us to refer to the Schwarz-Christoffel transformation for comparison purposes. The curvilinear net of C1 and [3 contours is obtained by both the conducting paper and the numerical methods; a comparison shows the two plots to be indistinguishable except of course inside the inner boundary. Reference to Fig. (A.2-2) shows that the convergence of the Schwarz-Christoffel series on the boundary of the square is pro- iwibitively slow, indicating that the numerical or analog mapping iarocedure can be used to advantage for polygonal regions as well as for arbitrarily shaped regions. To solve the boundary-value problem, Eq. (2.2-IA), we need as find the function LJlC:><:7E;3 and expand it into a Fourier series. This is accomplished by both the numerical and conducting sheet methods as described in Sections 2.6, 2.7. A comparison of the results obtained by the two methods shows good agreement. 30 3| 3.2 Conformal Mapping_by the Conducting Paper The discussion of the application of the conducting paper technique to the conformal transformation of a square region will be subdivided into four parts: (a) equipment, (b) construction of models, (c) test procedure, and (d) calculations. These four areas will now be discussed briefly here. As a reference we suggest Karplus and Soroka (I959) and Bradfield, Cooper, Southwell (I937). 3.2.l Eguipment The conducting sheet paper is type L Teledeltos paper approximately 0.00h Inches thick, having a resistance value of approximately #000 ohms per square, which has been found by the manu- facturer to be very well suited to use in field plotting. it consists of carbon-graphite-impregnated stock; and the uniformity of resistance of sample areas has been found to run within the range of from :3 per cent to': 8 per cent maximum deviation from a mean value, with the paper showing the lower values when measured lengthwise of the roll as compared with the higher values obtained when measured crosswise. By making each equipotential line plot twice, with one plot rotated 90 degrees from the other relative to the plotting paper, and averag- ing the two positions found for each plotted point, we were able to locate points readily to within‘: 0.0I inches with the instrument Plotting circuit used. The control unit consists of an instrument used as a zero voltage or "null detector,” D-C power supply, and voltage divider; it is constructed by the Sunshine Scientific instrument company and ‘5 designated as catalog No. Zh-IB. The voltage divider consists of 32 ten turns of a helical-wound resistance element of rated linearity ofl: l/2 per cent; hence a fine and accurate adjustment of the tracing voltage level anywhere between zero and ICC per cent of the voltage is possible. The low resistance electrode areas on the region boundaries may be conveniently established on the conducting paper by means of silver conducting paint. Drying time for the painted areas is over- night. Since it is essential that the electrode areas be of very high conductivity, we used bare copper wire held in contact to the paper by the silver conducting paint as a safeguard against any voltage drop in segments of appreciable length along the boundary of the model. 3.2.2 Construction of the Models We desire to obtain superimposed plots of both the equipo- tential lines and streamlines on a single orthogonal map. The family of curves consisting of equipotential lines can be traced on the conducting paper using the control unit described above. The other fannly, the streamlines, could then be constructed graphically as the ”orthogonal trajectories” of the first ones. However, it is more convenient to construct the streamlines experimentally. This is accomplished by inverting the first model (first-stage model) used to construct the equipotential lines, so that the equipotential and streamlines are interchanged. The inverted model is called the second- Stage model. Mathematically, we now must redefine the complex poten- tial function of Eq. (2.li-2) by allowing F( 2) = [3+ iC(. Plotting equpotential lines on the second-stage model corresponds to plotting 33 2 Units l: =% -05 Units Lead".-“§4 ' --«*\e/~\.»-w~~. __ .05 Units ) if 05 .I. i . Units-T l 2 Units i i . h\\\\q\ Electrode ' Boundaries Lead Fig. 3.2-l. First-Stage Model for Square Region Lead 4”,,'Drawing Board x" / ’ 1 /A, . 1 3‘7 // fl/AV /_ i~'———-.I .05 Units .4. \\\\\\\\\\\\\ \ l Unit / Cardboard Backin V 9 Electrode -—-"""" L/,/’ Boundaries Lead Fig. 3.2-2. Second-Stage Model for Squar e Region 3h Control Unit Stylus . \ ~ '5 ‘ ./ r/. Graph Paper Fig. 3.2-3. Analog Field Plotter with Pantograph 35 the streamlines on the first-stage model. The inverted model is constructed by making its electrode boundaries coincide with two known streamlines of the first-stage model. The first-and second-stage models for the square are illus- trated in Figs. (3.2-l) and 3.2-2). The two known streamlines selected from Fig. (3.2-l) to be used as equipotential boundaries in Fig. (3.2-2) are the rays (3 = 0° and C) = “5° emanating from the origin; these must be streamlines because of the symmetry of the square region. The plotting of points for both families of curves is accomplished by the same test procedure applied to these two models. 3.2.3 Test Procedure The arrangement of the equipment is illustrated in Fig. (3.2-3). A one-to-one pantograph consisting of a double parallelogram linkage pivoted at the center of the plotting board was used to record the plotted lines.. The test is conducted as follows: By setting the potentiometer, the searching probe or tracing stylus is adjusted to the desired potential level of the equipotential line to be traced. 11: take a reading, the tracing stylus is drawn smoothly over the paper Starface, and at the point where the pointer of the null detector svvings through zero, the recording stylus is pressed into the record- lrtg paper sufficiently to register a permanent mark. 3.2.h Calculations Thea,/3 net for the square region is illustrated by Fig. (3.2-b), in which, because of symmetry, only l/8 of the net need be :6 3 ‘A Table 3.2-‘4. '/ (Continued) 37 38 shown. To allow for the anisotropy of the paper (lengthwise-crosswise), one-quarter of the square was tested; Fig. (3.2-h) represents averaged values. The values for the constants A and a of Eqs. (2.h-8) and (2.h-9) may be calculated in the following manner (see Bradfleld, Hooker, and Southwell, l937): For the square, the equipotential lines for r<:i/2 unit or J£32>0.20 are indistinguishable from circles; thus the streamlines for r<:l/2 unit are radial. For this central region we write in analogy with Eq. (2.li-6) that B a: - A Log i: s (3.2'i) i i where A ‘21-” g and d is anvarbitrary constant. Let d =C.e/A; then Eq. (3.2-l) becomes 3 ‘ - . Logr LogC+K(l U). (3.2 2) Hhenlfl3 - l, we observe that C is equal to the radius of the inner boundary of the square. If observed values of Log r and (l -/3 ) are plotted, we expect the points to fall in a straight line and the values of'A and C may be evaluated. From Fig. (3.2-5), we find that C = 0.0250 units and i. - 3.77. Considering the manner of its derivation, this value for C Is very satisfactory; the measured radius was 0.0250 units. The inner radius (a) of the circular annulus, to which the square is mapped, is determined from Eq. (2.h-9) to be 0.002h3 units. By example, we have demonstrated how by means of the electric current flow in a sheet of conducting paper, conformal transformations can be obtained experimentally. The alternative procedure of numerical vco < mucmumcou mo co_aec_Eeouoo acusmcooxm coo-e mc_uu:pcou ,.mi~.m .mwu 39 .u AMT: m.o m.o 0.0 m.o :.o m.o «.9 _.o p n pl - L h p b A a A a 1 q q - 110.:I ma.m- u m~o.o mo; 0.— NN.M LIO.MI o\0 o\ O\ tro.~i o\ o\ o\ o\ \O\ 1.10.... o \ if 0 J 501 hO conformal mapping will be applied to the square region in the next section. 3.3 Conformal Mapping byiNumerical Means in this section we will describe the method used to obtain the point-to-point mapping of the square into the unit circle by the numerical procedure of Section 2.h. The calculations were performed utilizing the electronic digital computer (MISTIC) located at Michigan State University. The results, values of(:( and /3 at each grid point, are listed in Table (3.3-l); only one-eighth of the square is shown because of its symmetry properties. A series of programs were written to accomplish the mapping; they will be discussed briefly here. Program A This program evaluates the boundary values of v, where v = Log sz + y2 on the external boundary of the given region. These values are to be used in the solution of the Dirichlet problem for v (x,y). Program 8 With the boundary values for v, a complete program from the MISTiC library (G-l) is used to solve the Dirichlet problem. The solution is by the Liebmann process (see, for example, Salvadore and Baron, I952). This method uses the finite difference approximation Q) n2 v2 - ® (@DCD + o 012); (3.3-!) to the Laplace operator the solutions are achieved by iteration (see Appendix B). hi Program 8 evaluates the harmonic function v at the lattice points of a square grid superimposed on the given region. Next we evaluate its complex conjugate U(x,y) by the Cauchy-Riemann equations. Programs C, D, E, F prepare for and solve the Cauchy- Riemann equations. PrOgram C 0v This program evaluates 5— at each lattice point. The method x used is Neville's (l93h) repeated root method of successive linear interpolations. Program 0 This program rearranges the output of Program C for use in Program E. Program E To minimize the error incurred by the Simpson's rule integra- tion of Program F, this program uses Neville's method to interpolate for values of Q1 between the lattice points. x Program F This program evaluates u, the conjugate harmonic function of v, at the lattice points by the Simpson's rule integration of Y u 6L 33—1. dy + const. The Simpson's rule approximation (Salvadore and Baron, I952) is A =-i31(ro + M, + 2r2 + M3 + + 2r",2 + Mm, + f") + 0 on“). (3.3-2) #2 With u and v evaluated at the lattice points, CI .‘[3 can be calculated by Eqs. (2.5-5). (2.5-7). This is accomplished by Programs 6, H respectively. Program 6 This program evaluates (I , from C( - @+ u, at the lattice points within the region, using the results of Program F for u. Program H This program calculates [3 , from [3 t - Log r + v, at the lattice points within the region, using the results of Program 8 for V. Program i This program calculates the values of x. y for points equally spaced on the boundary of the unit circle; these will be required later in Program J of Section 3.“. This is accomplished by inverse interpolation of the results of Program 6 ( C(-C((x), a - C((y) ) for equally spaced values of CI . Again Neville's process was used. With the values of CI and [3 known at each lattice point, 'see Table (3.3-l), we have a polnt-to-polnt conformal mapping between the square and the unit circle. The curvilinear net of the values of and D , normalized so that [3 - l at r - 0.025 units to compare it with Fig. (3.2-h) shows remarkable agreement with the results of the conducting paper plot; within plotting accuracy no differences were observable. (it should be noted that the conducting paper plot. Fig. (3.2-h). does not exist for r<:0.025 units) For this reason these curves are not repeated. ill/“6 lB/lh lZ/lh ll/l# l0/lh 9m. 8/lh 7/11. 6/114 5/”. am. am. Z/IH l/lh 0 00 l/lh .785398 2.367837 0 2/lh .785398 l.67h8l7 .h6369h l.909703 O 2.7lh399 2.02l229 #3 3/lh .785398 .26987l .58823h .h32l78 .32l935 .563072 0 .6l56h8 h/lh .785398 .983562 .6hhlh7 .l06023 .h6h386 .2l6787 .ZhShhl .297515 0 .327642 S/ih .785398 .763286 .676l23 .8606l3 .Sh2265 .95252h .382l23 .030729 .l98320 l.08hh85 0 l.l03806 6/lh .785398 .586ihh .697266 .665760 .59l687 .7h25hl .h6738h .8l2250 .32Q707} .869l77 .l66767 .906930 0 .9202l3 Table 3.3-l. Values of C1 and [3 at Grid Points in the Square Region from the Numerical Mapping Technique. 7/lh .785398 .hh0h63 .7l2789 .506h32 .626675 .5708l8 .526235 .63l280 .hll352 .68h6l8 .283l5h .726973 .lhhh89 .75hh2h 0 .763956 X lh/lh l3/lh l2/lh llllh lO/lh 9/lh 8/lh 7/lh 6/l4 S/lh h/lh 3/lh 2/lh l/lh 8/lh .785398 .319809 .725l89 .37hhhl .653725 .h28l7l .57052l .h79655 .h75h32 .527029 .36893l .567956 .252386 .599830 .l28250 .620l75 O .627l79 9/lh lO/lh .785398 .Ih0613 .785398 .7u54uu .220538 .l759ll .735812 .696209 .2651u6 .210942 .676252 .637576 .309248 .2h5280 .60645u .5695hl .352072 .278295 .5263h7 .h923on .392523 .309lhh .h36l9o .h063h7 .h2918l .336788 .336739 .312563 .h6o37h .360055 .229388 .212297 .uau337 .377759 .ll6236 .l0737l .h99h7h .388862 0 o .50h656 .3926h8 Table 3.3-l. ll/lh .785398 .07896l .75h555 .l05303 .7lh772 .l3l520 .6660lh .l57397 .608325 .l82600 .Shl875 .20665h .h67028 .22895h .38hh05 .2h8775 .29h9h5 .265332 .19993l .2778h6 .l00986 .285657 0 .28883l (Continued) lZ/lh .785398 .035077 .763h25 .0526ll .7326h9 .070097 .693076 .08hh69 .6hh76h .l0h506 .587858 .l2l0h5 .5226h0 .l36752 .h49573 .l5l23h .36935l .l6h038 .282928 . 171+678 .l9l538 .l82687 .096667 .l87669 0 .l8936l l3/lh .785398 .00877l .772207 .0l7538 .750250 .026293 .7l95h0 .0350l2 .680l26 .0h36h7 .632ll8 .052ll8 .575722 .060309 .5ll279 .06806h .h39306 .075193 .36053l .08l#75 .27S9lh .08668l .l86653 .095890 .09hl57 .0930l6 0 .093839 lh/lh .785398 0 .780009 0 .766889 0 .7h5009 o .7lh39h 0 .675108 0 .627278 0 .S7ll29 0 .507023 0 .h35h9h 0 .35728l O .2733H5 O .l8h87l 0 .0932hh O 0 0 #5 For the solution of the elasticity problem we need to have, in addition to the point-to-point mapping, the values of (“flax/(C) ‘ \ on the boundary, where z = U(C) is the mapping function and C8 6" ea. on the boundary of the unit circle. L] : lgjon the Boundary of the Square Region 3.“ Evaluation of in this section we evaluate the function U(%-{d—') on the boundary of the square region by the numerical and conducting paper methods described in Sections 2.6, 2.7. The Fourier coefficients of this function are required in the solution of the boundary-value problem, Eq. (2.2-in). To evaluate the function numerically we again utilized MISTIC to solve Eqs. (2.6-2), (2.6-3). A brief review will be given of the programs devised for this purpose. Program J This program evaluates the values of (935‘ ' (9.5) needed deiiz d 1" for Eqs. (2.6-3). It uses the results of Program I, which calculated the values of x, y for equally spaced intervals on the boundary of the unit circle. The derivatives are calculated by Neville's method of successive linear interpolation. Program K This program uses the results of Program J to evaluate Eqs. (2.6-2), (2.6-3), thus determining the values of LJ(€}<:/ejj at equally spaced intervals on the boundary of the unit circle. Q6 Fourier Coefficients from Conducting Paper Data Fourier Coefficients from Numerical Data \DCDNO‘UT-PUN—O 7s WWWWWWUWUWNNNNNNNNNN——-————I—-——n— \ooowoxmrwN—ommummkwn—ommwmmpwn-o bk b_k bk b-k +.93868l ~93'635 .5‘7737 +.522h70 -.D7l607 .055l59 .l35480 -.l3l654 +.023809 .02#03l .069823 +-°6'387 -.Dlll60 .009893 .ouuese - 0“#423 +.006038 .00fi965 .D3l9lS +.0293h8 -,oo3h7h .003387 .02h369 - 0232“' +.0020l2 ~00390‘ .0l9h6h +.0l9570 -,oo|101 .00l938 .Ol606l -.Ol7030 +.oooh95 .000710 .0l3583 +.0lh33l -.oooo7l ' -000602 .013709 -.0l2263 Table 3.h-l. Fourier Coefficients for LJ(C7) 6771;) of Square Region Obtained by Numerical and Conducting Paper Mappings. .co_mom ocoaom mo >cooc30m ecu mo mucom xcoewmoe_ one .eom ..u:.m .m_u menace: .mu_coE:z co» noncom co_c:ou oououeach nu mc_aaot .eu_coeez.iii Ill! cocoa mc_uu:ocou mu [1" m..- 9..» m.on m.o o.— m.— h8 Program L This is a complete program from the MlSTiC library which evaluates the Fourier coefficients of a given function. For the mathematical steps involved, see, for example, Chapter 8, Muskhelshvili (l953). The coefficients for L“/K;Z:7E5 are listed in Table (3.h-l). For comparison purposes, we have plotted the real and imagin- ary parts of WWW) obtained by: (a) numerical procedure. (b) conducting paper method, and (c) Fourier series from numerical results where the coefficients bk are truncated for k = :_39. Only one-half of the plot for Uta/U15) is shown, since the real and imaginary parts are even and odd functions respectively. The agree- ment between the results of the two methods is very good. The maximum deviation between the numerical and conducting paper results is three Per cent. The truncated Fourier series for (e/«i}€3%;5 shows excel- lent agreement with lts limit value except (not unexpectedly) near the value of 6 corresponding to a corner of the square. The Fourier series converges for all values ofe , since the function satisfies the Dirichlet condition, but in the corners the convergence is very slow. 3.5 Conclusion We have demonstrated by means of an example that the conducting paper conformal mapping is just as accurate as the numerical procedure for the point-to-point plotting and almost as accurate for the deter- mination of WWW on the boundary. ' The conducting paper method is simpler to use and more economical than an electrolytic tank, and the results given here show that when Q9 used carefully it is quite accurate. The numerical procedure is slightly more accurate and has the advantage of offering an almost completely automatic solution to the mapping problem. This advantage can only be obtained if automatic high-speed digital computing facilities are available as well as personnel trained in their use. The conducting paper plotting is much simpler and can be learned very quickly. For the solution of the elasticity problem presented in Chapter iv, a digital computer is almost a necessity in one step and very helpful throughout. The numerical mapping method is therefore especially suitable as a part of the total procedure, while the con- ducting paper then provides a convenient check on the numerical mapping. CHAPTER IV EXAMPLE - STRESS SOLUTION h.l Statement of the Problem in this chapter, we shall apply the complex variable method with numerical mapping to a specific elasticity problem. To demon- strate the effort involved and the accuracy of this method, we con- sider a problem that admits an exact solution. For further comparison, results are also obtained by the finite-difference approximation to the biharmonlc equation. Specifically, we consider a square region loaded along the center quarter of one side by a distributed load of intensity P0. The loading of the three other sides is taken to be identical with the stresses that would act there if the square were a portion of the semi- infinite plane. The stress distribution in and on the square can then be calculated by superposition techniques (see, for example, Timoshenko and Goodier, l95l). This problem was selected because the exact solution permits an estimate of the accuracy of the numerical methods in a problem where there is a region of rapidly varying stress (near the jump discontinuity in the boundary traction). Although this loading does not satisfy the sufficient conditions given by Muskhelishvili under which the solution of Eqs. (2.2-l5) and (2.2-l6) furnish series for $(C) and V(C) satisfying the conditions of the problem, comparison with the known exact solution of the present problem will show that in this case at least the conditions are satisfied. 50 Sl h.2 gxact Solution by Superposition integral ds "'il" Po +c -c -“’-" -‘-’_ ""1 Y 9 i A. ,i -—-i _ i fl (y-S) i (if) : I i i ii Fig. h.2-l Represented in Fig. b.2-l is the half plane with a uniformly- distributed load of intensity Po over the interval (-c, +c). The stress function for a single concentrated force on the half plane (Section 3.2, Timoshenko and Goodier, i95l) is 3.9. U = - ,. r19 Sinf} . (h.2-l) For the infinitesimal distributed load Po ds we write Po ds dU3-TreSine. (“.2-2) We obtain from Fig. (h.2-l) the following identities: r Sine - y - s, i. = r Sine, r2 1-sz + (y - s)2. and a - Tanml 1.3. (“.2-3) X The barred coordinate variables x,y are used to distinguish them from x,y coordinates to be introduced later in the square region, with origin at the center of the square. 52 Substituting Eq. (14.2-3) into Eq. (14.2-2) and integrating over -c_‘_s.‘.fl:, we obtain u - £274 [hi-dz + £11..." (77.15) - [(;+c)2 + .22] (“J-h) \ Tan"I 2?; + 2c!) . Eq. (16.240) represents the stress function over the half-plane, from which we may now determine the stresses at any point. We represent the stress in terms of the stress function U (see Section 2.2) as —_qu - du -_-$U ,- (5'5? ’ “‘07 ’ TJ‘Y‘W (”5) where the bar over the stresses denotes they are referred to (x,y). Substituting Eq. ((0.240) into Eq. (ill-5) there results _ p _ -_ _ .. - cy-'_O Tanl'YTE-Tanl-Y-Ec-+f +5 (“.2-6) 77 x x x +(y+C) ‘ Po -l "-c -l '+c R '+c 0x = 7?; [Tan (if?) - Tan (13-) - #32 (14.2-7) i '-c + 2§+<9-.)2] ' f .-.':2 .23... -_2.’__ - xy [2249-02 22““), . (m a) 53 For the problem specified in Section h.l, c equals L/8. Eqs. (4.2-6), (4.2-7), and b.2-8) enable us to calculate the stress anywhere in the square. The tractions on the edges 9 e L, y a f; 2 and R - L, evaluated from this exact solution for the semi-infinite plate, will be used, together with the prescribed traction on the edge R - O of the square, as boundary conditions for the Muskhelishvlli complex variable method of solution in the square. 4.3 Solution by the Muskhelishvili Complex Variable Method This solution involves the evaluation of f1 + ifz from the given boundary traction, the determination of<¢5, yyfrom Eqs. (2.2-l5), (2.2-l6), and the calculation of the stress distribution from Eq:(2.2-l7). The numerical mapping and evaluation of 5/6i}<:7€;) for the square region have been calculated earlier-——Section 3.“. Fig. (h.3-i) shows the coordinate system used here and a sketch of the boundary tractions according to the exact solution of Sec. “.2 . Programs M, N, L prepare for and calculate the Fourier coefficients of f] + ifz. Programs 0, P, Q prepare for and calculate 9b(z), yy(z) at the grid points in the given region. Then Programs R, S calculate the stresses at the grid points in the z-plane. Program M With the boundary traction given, this program evaluates f} + ifz, (see Eq. 2.2-3) at equally spaced intervals on the boundary of the square. The integration is performed by Simpson's rule. However, to express f. + ifz in a Fourier series, It must be evaluated over equally spaced intervals on the boundary of the unit circle. Sh Program N This program evaluates fl + if2 over equally spaced intervals on the boundary of the unit circle. From the conformal mapping we know CZ,- C: (s), 5 being the arc parameter on the boundary of the square; and from Program M we know f‘ - f2(s), f2 - f2(s). For a given value of C1 on the unit circle, its corresponding value of s is calculated from (3( - CI (5) by inverse interpolation; then f1 and f2 are calculated from f‘ - fl(s) and f2 = f2(s) by interpola- tion. This is repeated automatically for C1 equally spaced on the surface of the unit circle. Program L was then used to calculate the Fourier coefficients Ak, A-k of f‘ + ifz. They are tabulated in Table (b.3-l), truncated with k = 39. With the values for Ak' A-k and bk, b-k from Tables (4.3-l) and (3.h-l), we can proceed to solve Eqs. (2.2-l5) and (2.2-i6) for the power series coefficient of QD , 3%[ respectively. This operation is Performed by the use of Programs 0, P and Q. Program 0 This program, a standard library routine, solves the system of linear simultaneous Equations (2.2-l5) for ak by the “normal direct” method. The augmented matrix is reduced to triangular form and the roots of the equations are obtained by back substitution. The program is for real numbers only so that the complex equations must be separated into real and imaginary parts. 55 Program P This program evaluates the constants a'k from Eqs. (2.2-l6), one at a time. A series summation of 39 terms is required. We may now evaluate ¢ and ill. Program 9 With ak and a'k evaluated in Programs 0 and P respectively, we can now calculate (I), V] at the grid points in the square from Eqs. (2.2-l0). Series summation of 39 terms is involved for each value of ¢ , VJ. With ¢(z), i/f(z) known at each grid point, we are able to find the stress distribution from Eqs. (2.2-l6). Programs R and S are constructed for this purpose. Pragram R This program differentiates a function by using the Lagrange differentiation formula, Nielsen, 0956). Higher order derivatives are obtainable by repeated use of this program. Thus the values of ¢'(z), ¢u(z), 51/22) are evaluated at the lattice points in the z-plane. Program 5 This program evaluates the stresses at the lattice points in the z-plane from Eqs. (2.2-l6), utilizing the result of Program R. By means of Programs M through S, we are able to evaluate the stress at the lattice points in the z-plane. For a portion of the square, these stresses are plotted in Figs. (h.h-2) to h.§-5); these figures also show the exact solution and the finite-difference solution S6 .025'97 ,l0392u .l020l8 Po - Unit Load i ._ ________..+_ ________ _. .__ Q- X .l57520 (a) Normal Tractions Y .050l69 1079550 .050l69 (b) Shear Tractions Fig. 4.3-l. Boundary Tractions for the Square Region, Traction in (Unit Load)/(Unit Length)(Unit Depth). Square of Dimension 2 Units by 2 Units. 57 k A+k A-k D -.O75408 i —.l60794 +.078620 2 -.036Sl8 +.036660 3 -.022973 +.02l259 4 -.0l600l +.021679 S -.00337l +.OO92l9 6 -.0076h3 +.006429 7 -.00h375 +.00h706 8 -.001820 -.000097 9 -.003242 +,001276 i0 +.0008l7 -.00005h ii +.00098h -.00ll8h i2 +.0009li +.000050 i3 +.002587 -.00l606 lh +.000830 -.00l3hh i5 +.00l0l7 -.00086Q l6 +.00l090 -.00l67h l7 -.000l78 -.000h12 l8 +.000667 -.00029h l9 +.000i26 -.000250 20 -.000357 +.000752 2i +.000l32 +.000263 22 -.0008h9 +.000565 23 -.000750 +.00085# 24 -.00062l +.000333 25 -.00lih2 +.000858 26 -.000458 +.000683 27 -.000532 f;000hh2 28 -.00055# +.00077S 29 -.0000lh +.000229 30 -.000398 +.0002i5 3i -.000l6l +.000239 32 +.0000l9 -.000l9h 33 -.000282 +.000ll4 3h +.000l3h +.ODOOl8 35 -.000009 -.000060 36 -.000lh8 +.000293 37 +.000063 +.000072 38 -.000368 +.000239 39 -.000358 +.000hl9 Table h.3-l. Fourier Coefficients for f' + ifz. Traction on Square Region. 58 for comparison. Very close agreement is observed between the exact stress values and the values obtained by the complex variable methods except on the boundary itself and in the neighborhood of the corners of the square as expected. h.h Difference Solution The finite-difference approximation to the biharmonlc equation of Appendix B is used to evaluate the stresses within the square region. This enables us to make a comparison between the stress values obtained by the exact solution, by the complex variable method, and by the finite-difference method. The square region was subdivided into a square net of six by six since, as explained in Appendix B, the resulting linear simultaneous equations are limited to 39. Because of the symmetry of the example considered, a finer mesh could have been used for the half-square. But since no use was made of symmetry in the numerical conformal mapping solution, it was thought that a fairer comparison of the two methods would be obtained by not taking advantage of the symmetry in the finite-difference solution either. in applying the biharmonic equation in difference form at a mesh point in a row adjacent to the boundary, we need the values of the stress function U not only on the boundary but also at a fictitious point one mesh interval outside the boundary. The values of U and its normal derivative on the boundary are obtained by integration of the given boundary tractions along the boundary (see Appendix B). The value at the fictitious point e outside the boundary (see Fig. h.h-i) is then approximated by using the first central difference formula at point b 59 on the boundary as 519. ”e = Ui + 2h (ah )b , (“.“-i) where i is the interior point and h is the mesh interval. The boundary values of U for the square, loaded as in Fig. (“.3-l), are presented in Table (“.“-l). Fig. “.“-l. Finite-Difference Approximation at Boundary Solution of the system of simultaneous equations resulting from the application of the biharmonic operator, Eq. (B.2-“), was accomplished by Program 0, described in Section “.3. This gives us the values of U within the region, which are also tabulated in Table (“.“-l). With these values of U, the stresses were calculated from Eqs. (8.2-3). To obtain the stress at points coincident with the grid used in the complex variable approach, the stresses were inter- polated and are plotted in Figs. (“.“-2) to (“.“-S). it is quite apparent that particularly in a region of rapid stress variation the stness distribution calculated by the finite-difference method cannot hope to show the actual stress gradient with good accuracy, while the numerical mapping results are much closer to the exact solution. 60 .eo_u:_0m oucocowu_ououwc_u mo mu.:mox .momocucOLem c_ new .3 $0 o:_e> >ceoceom ;u_3 eo_mo¢ acmeum mo m_e: .co_uo:ou u_coEcoz_m on co_u:_0m eunucomm_nuou_c_u ._u:.: o_nnk X Aewmmmw..v My. Ammm:no nu ~m~m_mo.n~ flowed". a“ Ham—nw__uv AMmeMN_uv w_wa~.u Am_m_:m.-v A_:mmac.-v mmoxwo.- A_:mmco.-v Amsmaa_.-v xm~.m:..-v A_emnm_.-v A_o_o:~.-v smm m~.- Ansnmsm.-v A~mmmng.-v m.:m«..- ANwmmw..-v A~moqo_.-v Amooua..-v Aw.mqm~.-v Amsemmw.-v sum ~m.- Amemosn.-v $033.; an - - - - - . u - s... - 3a.- seas; seasons Echoes $3.3; a: are siesta sauces tars; 6| y l e D Y = 0 o—o-a-Q-o-oQ-e—o—o~5 C] {2 C) D - 0.5 _i x = l3/l“ units {2 \ o 9\g_é ‘C>-<>-E3-{3“-Eizsi5»-§}a=g3‘\€3\\ E] 23 \ E] g C) {3 \ -. . -— x = l2/l“ units a} \ o\ ('3 °~o - ‘—- -{3-—ES :r E3 E3 53’ 15:3 €5-{3\\ 6\ as E} \ - 0.5 — E} x - ll/l“ units ‘\\E5 (3 \\:3‘_§5 <> Exact Solution C3 Numerical Mapping O Fi ni te-Di fference Fig. “.“-2. Stress Distribution in Square Region - G‘x. Coordinates Refer to Fig. “.3-i. Stresses in Unit Load per Unit Length per Unit Depth. 62 Cl y B l . y = 0 E 5:65“ 6‘9 6\ D a\{} \D -0.5 '— x I lO/l“ units \ o\ op use W” \ Q\ a\ 0 a \ , 9\ -o.5 - x = 9/l“ units °\D 8 E -9 “9‘00~o\§ \ Ci §\ g’0. B \o . \o\ 3 x = 8/l“ units 5~ <> Exact Solution 1:] Numerical Mapping O Finite-Difference Fig. “.“-3. Stress Distribution in Square Region - 0x. Coordinates Refer to Fig. “.3-l. 63 E] C] Y ' i D y -= 0 otdz-o- _ __ .— g 8 ()5 B‘a\a \ x = i3/i“ units g}‘~g} C) \ 0.5-‘ (3 [3 [J E; —* 6—5-6: - c 83’ ‘5‘155“"E}-g5 x = l2/l“ units ‘9_ _. O E] D a 57-75:: - \ £3 g}‘Z;S-'S}“£5 x =3 il/“f units \b\§’g~.a\ a 0 *—€} 8... i D x =- lO/l“ units a x m 9/l“ units x = 8/1“ units <> Exact Solution [3 Numerical Mapping (D Finite-Difference Fig. “.“-“. Stress Distribution in Square Region -’ 6;” Coordinates Refer to Fig. “.3-i. 6“ ’5' o D Q -_ “'0" ”609:9;fi‘ - ’53 Cl E)\QO\ 6’6 9 x = 13/”. units G‘Q/ 0—086-‘9 2*.51‘6 {3 . C] x = l2/i“ units \\g}_’g5,/ W—o~—G ‘3 Q\ a ,., x = ll/l“ units O\{} -c” €§:: ——g§ ‘9‘ _ / 9 OQ‘Q‘Q\5 O 6 ,B/ x = 9/l“ units “§}"€3"£3 x = 8/l“ units <> Exact Solution D Numerical Mapping C) Finite-Difference Fig. “.“-5. Stress Distribution in Square Region - 7ky. Coordinates Refer to Fig. “.3-l. 65 The numerical mapping solution using an iteration method can use a grid with a maximum number of about 950 points on the MiSTlC computer as compared to a maximum of 39 available for the finite- dlfference solution with MlSTiC. “.5 Conclusion The numerical conformal mapping method for treating plane elasticity problems has been applied to a square region with a unit distributed load along the center quarter of one side and with the loading of the other three sides taken to be identical with the stresses that would act there if the square were a portion of the semi-infinite plane. This problem, which has an exact solution and includes a region of high stress gradient, was also solved by the finite-difference method. The results are plotted in Figs. (“.“-2) to (“.“-S) where a comparison in accuracy can be made between the exact stresses and those obtained by the complex variable and finite-difference methods. The finite-difference solution seems to produce "average” values of stress and, because of the coarse grid, cannot hope to demonstrate the stress gradient in the vicinity of the boundary load discontinuity. The complex variable solutions agree remarkably well with the exact stress except on the boundary itself (where two of the three components of stress are already known) and in the vicinity of the corners. Poor results on the boundary are due to the fact that the power series for 9b», y) converge slowly there. The results at the corners are, in addition, affected by the slow convergence of the Fourier series for U(GyJi-Gti . 66 By example, we have shown that the Muskhelshvili complex variable method with numerical mapping furnishes very good approxi- mations to the stress distribution, and that the entire procedure can be set up on an electronic digital computer. CHAPTER V SUMMARY AND DISCUSSION A method has been demonstrated for treating plane elasticity problems for simply connected regions. It consists of use of numerical conformal mapping methods to map the simply connected region onto the unit circle in order to apply the Muskhelishvili complex variable method. This solution consists of two parts: the mapping of the given region onto the unit circle by numerical means, and the solution of the boundary value problem on the unit circle by the method of undetermined coefficients. It is believed that this is the first time numerical mapping techniques have been used for the solution of plane elasticity problems by complex variable methods. This approach now makes the whole Muskhelishvili complex variable method susceptible to automatic solution on a digital computer, pro- vlding greater accuracy than automatic finite-difference equation solutions. To demonstrate that this method is not difficult to use in practice and to make some estimate of the accuracy, it was applied to an example involving a large stress gradient and for which the exact solution was obtainable. For further comparison, the stress distribution was also obtained by the finite-difference method whose accuracy was limited by the capacity of the computer to solve, by direct elimination, systems of linear algebraic equations, since automatic iterative solutions do not converge for the finite-difference approximation 67 68 to the biharmonic equation. A comparison of the results shows that the complex variable method with numerical mapping furnishes better results than the finite-difference solution and agrees remarkably well with the exact solution except on the boundary and in the neighbor- hood of the corners. The numerical procedure should be of great practical value to stress analysts, since it provides for the automatic solution of plane elasticity problems. The conducting paper techniques for the conformal mapping were used in this study mainly as a check on the numerical procedure. However, the accuracy achieved by this approach suggests a possible modification in our method of solution for some problems. Suppose the functions fl + if2 and MiG-ya; are continuous and have continuous derivatives of a high order; their Fourier series repre- sentation would then converge rapidly. if the conformal mapping and the function “(Ci/m are obtained by the electroanalogic method, the remaining numerical calculations to determine (P , 31/ and the stresses can then be accomplished easily without resorting to an electronic computer, since the number of simultaneous equations to be solved will be small. This then provides a method of solution for some elasticity problems not requiring the use of a digital computer. While for such problems one could also expect a coarse mesh to be satisfactory for finite-difference solutions of the equations, a coarse mesh presents difficulties with the boundary conditions on curved boundaries. The electroanalogic mapping method is, on the other hand, easily adapted to curved boundaries; and for curved boundaries it may be advisable to perform the mapping in this 69 way even when the boundary tractions are discontinuous and a digital computer is required to calculate the large number of series coef- ficients needed. For these reasons, further application of electro- analogic methods seems worthwhile. The Muskhelishvili method has also been adapted to the solution of the plate equation for laterally loaded thin plates (see, for example, Muskhelishvili. l953). The present numerical mapping methods may therefore find application to plate deflection problems as well as to plane elasticity. lO. ll. l2. l3. l“. BIBLIOGRAPHY Bergman, S. and Schiffer, M., Kernel Functions and Elliptic Differential Eguations in Mathematical Physics, Academic Press, Inc., New York, I953, p. 78? Bradfleld, K. N. E., Hooker, S. 6., and Southwell, R. V., l'Conformal Transformation with the Aid of an Electrical Tank,” Proc. Roy. Soc. Lond., Ser. A., vol. l59. I937, pp. 3IS-3“6. Churchill, R. V., introduction to Complex Variables and Applications, McGraw-Hill, New'York, I9“31 Faddeeva, V. N., Computational Methodgof Linear Al ebra, Dover Publications, New York, I959, pp. ll7-i53. Higgins, T. J., "Electroanalogic Methods,” Applied Mechanics Reviews, vol. 9, no. 2, I956, pp. “9-55. Kantorovich, L. B. and Krylov, V. i., Approximate Methods of Higher Analysis, Interscience Publishers, Inc., New York, l9'B. Kellogg, 0. 0., Foundations of Potential Theory, Frederick Ungar Publishing Company, New York, l929. pp. 365-370. Dover Publications, New York, i935. Karplus, w. J. and Soroka, w. w., Analog Methods, McGraw-Hill,. New York, I959, pp. 386-“29. Miller, K. 5., An Introduction to the Calculus of Finite Differences and DiffErenceEguatlons, Henry Holt and Company, New York, I960. Mindlin, R. D., ”A Review of the Photoelastic Method of Stress Analysis,” Journal of Applied Physics, vol. l0, I939, pp ‘ 273-2%0 Muskhelishvili, N. I., Some Basic Problems of the Mathematical Theory of Elasticity, P. Noordhdff, Groningen, l953. Nehari, 2., Conformal Mapping, McGrawhHiIl, New York, l952, pp. l73-l32. Neville, E. M., ”Iterative Interpolation," Journal Indian Mathematical Society, vol. 20, l93“, pp. 87-I20. Nielsen, K. L., Methods in Numerical Analysis, Macmillan Co., New York, l956, pp. 3“2-3“5. 7O l5. l6. l7. l8. I9. 20. 2|. 22. 7| Parkus, N., ”Elastic Thermal Stresses in Delta Wings, Part l," Report Sponsored by the United States Government, Contract No. 6l(052)-2l“; l960. Salvadori, M. G. and Baron, M. L., Numerical Methods in En ineerin , Prentice-Hall. Inc., New York, I952, pp. lg7-2i6. Savin, G. N., Stress Concentration on the Boundary of Holes, Government Publication, Moscow-Leningrad, l95i (in Russian). German Translation, Veb Veriag Technik, Berlin, I956. English Translation, Pergamon Press, New York, l96l. Smirnoff, V. I., ”Uber die Randerzuordnung bei Konformer Abblldung," Mathematlsche Annalen, vol. l0“, I932, pp. 3l3-323 (in German). Sokoinikoff, l. 5., Mathematical Theory of Ejasticity, McGraw- Hill, New York, l956. Southwell, R. V. et al, "Relaxation Methods Applied to Engineering Problems,” Trans. Roy. Soc. Lond.. Ser. A, vol. 239. I9“5. pp. “19-577. Streeter, V. L., Fluid Dynamics, McGraw-Hill, New York, l9“8. Timoshenko, S. and Goodier, J. N., Theory of Elasticity, McGraw-Hill, New York, i95l. APPENDICES 72 APPENDIX A CONFORMAL MAPPING OF POLYGONAL REGIONS In this appendix we will present the theory of the Schwarz- Chrlstoffel transformation for polygonal region and then apply it to the square region. We will demonstrate that the slow convergence of the mapping function series on the boundary of the unit circle prohibits its use in the calculation of the function U(CyLZ—GT). required in the solution of the stress problem by the Muskhelishvili complex variable method. We also demonstrate that the Schwarz- Christoffel interior mapping function satisfies the Smirnoff condition, Eq. (2.2-9b). A.I Theory of Schwarz-Christoffei Transformations In this section we shall treat the conformal mapping of a general polygon onto the unit circle. The Riemann mapping theorem assures us that any two simply connected domains can be mapped con- formally upon each other. While a simple solution of the general mapping is not known, the interior of a polygon can be mapped onto the unit circle by the Schwarz-Christoffel formula (see, for example, Churchill, l9“8): z =A (Hum-11' (l-82§)'u2' (hang-“n dg + a, (A.I-I) where A and B are arbitrary constants, and ak and [L k (k = I to n) are defined in Fig. (A.I-l). 73 7“ z-plane 32 5‘ / [32 a3 é: \ >( ,z / n; -piane 54/3 * I Fig. (A.l-I) Mapping of Polygonal Region onto Unit Circle From Fig. (A.l-l) we see that r—e BK'ILJT' 2.. l3..-= 2W"- Zia.- 2. «42.57. - l< LLK|ai-aj >20>0,i#j, (A.3-2) 79 where D is a sufficiently small number. Then near any corner, say ak for which Mk > 0 we have ié' aJ| > 0 for all j qt k. (A.3-3) Since “1 fllg-aj| <1 for [,LJ>0, (A.3-“) ifk Eq. (A.3-I) becomes ‘wiél ‘ < N ( g- akleuk (A.3-5) where N '5 11)." and n = 2-Mk. (A.3-6) We can write the Smirnoff condition, Eq. (2.2-9b), as 0,45 I. = J lu'meie )i d9 k (A.3-7) 9W5 _ 9 d < "L [IPCos e + I (DSlne - Cosa‘ - i Sinek‘]uk k N i2 . Algebraic manipulation of Eq. (A.3-7) yields 9,6 9 d '2 " [ (,2. I - 2 p c... WP—iik- (”'3’ 9i. 80 Now, if we let (9- ... Eq. (A.3-8) becomes =5 ,9-9k=c< ,and(3=l-€ MIX. (5 do: I = 2 . 2 / [62+“(I-€)Sin2g]%uk O 2 If we let Eq. (A.3-IO) becomes 26 l2‘2 d7 % 0 [52+h(l-€)5in27’] #k For (5 sufficiently small we can write that Sin 2V ;2-%; , and Eq. (A.3-l2) becomes 25 i2_ SE 2 dir [62+(‘-€)72] 241k g 2 07 '5 2 _ Ti'uk ii-efi‘f‘kfl‘k (i-EIH‘LI‘ i-uk O s (A.3‘9) (A.3-ID) (A.3-ll) (A.3-I2) (A.3-l3) (A.3-I“) 26 8i In Eq. (A.3-I“), we see, since I -/.,Lk> 0 , as (S --—i-0 that I2 is bounded; consequently, i' is also bounded as 5—- 0. Thus the Schwarz-Christoffel formula for interior maps satisfies the Smirnoff condition, Eq. (2.2-9b). APPENDIX B FINITE-DIFFERENCE SOLUTIONS TO PARTIAL DIFFERENTIAL EQUATIONS B.I Discussion The finite-difference method is one of the most commonly used methods for solving boundary-value problems in partial differential equations. The literature is extensive; for bibliography see Miller (I960). In this method the differential equation is replaced by an approximating difference equation, and the continuous region in which the solution is desired by a set of discrete points. This permits one to reduce the problem to the solution of a system of algebraic equations, which may involve hundreds of unknowns. The solution of a large number of linear equations by direct elimination methods ls prohibitively long, and some iterative technique has to be devised to solve such equations. Iterative methods, however, do not always converge and cannot be applied to all systems. In this appendix we discuss briefly the solutions to the Laplace and the biharmonic equations. The convergence test for the Seidel iterative method of solution, for the simultaneous linear algebraic equations resulting from the Laplace and biharmonic equations, show that the first is borderline while the second Is divergent. Conse- quently, the simultaneous equations from the biharmonic equation must be solved by direct elimination methods. Digital computers are limited in 82 83 the number of equations they can solve by these means; the MISTIC is limited to 39 equations. This imposes a severe limitation on the use of finite-difference methods to solve the biharmonic equation. 8.2 Laplace and Biharmonic Equations An important problem in potential theory is the Dirichlet problem in which we are asked to find a function harmonic in a given region and assuming preassigned values on the boundary. If the domain is divided Into a square net with n interior points, n finite-difference equations are obtained, in each of which the operator V2 can be represented schematically as (see Salvadore and Baron, I952): this... -a~—-@ .1 . , . h "2 V2 “ @’.‘ GD . + 0 (hzl . (B.2-l) I i-o—i where h is the mesh interval, and the difference equation at point 0 becomes VI + V2 + V3 + V“ ' “V0 ‘ 0 (8.2-2) where VI' V2, V3, V“ are the values of V at the four nearest neighbor mesh points, shown in the operator "star“I In Eq. (B.2-I). The numbers 8“ in the circles of the star are the coefficients of the appropriate terms in the difference equation and the + D (hz) indicates that the error in approximating the differential operator by this difference operator is of the order hz. We can now solve these equations together with the boundary conditions for the values of the function at each grid point in the domain. Faddeeva (I959) shows that a sufficient condition for a System to be solvable by Seidel iteration is that in each equation the sum (N) of the absolute values of all the coefficients except the largest be less than the absolute value (M) of this largest coefficient .31 <:l ). The smaller the quotient |aiin comparison to unity, the more rapid is the convergence. “-1 For Laplace's equation we note that ‘g‘ = 3'- , so that this is actually a borderline case. Experience shows that convergence does occur but at a very slow rate. The compatibility equation in terms of Airy's stress function U (see Section "2.2) reduces to the biharmonlc equation We... .23.... on... + A ... ...-..— . 8.2.3) 0x4 ' ‘Oxzdy” 83'" ( subject to the boundary conditions Eqs. (2.2-3) 3 g-E. - [Yn(s) ds + const. 5": fl (5) i e (B.2-“) s we + fxn(s) ds + const. 2 f2 (5) i 5e where the traction components X n(s), Y n(s) are given in terms of 85 the arc length 5 on the boundary. The stresses in terms of U are given by Eqs. (2.2-i). 2 2 2 0x 3.9.2., 6‘ ..(2—9.’ 7" ....QJ_ (3.2-5) Y "Y 8Y2 OX2 dx y If we divide the region into a square grid of n inner points, n finite- difference equations are obtained, in each of which the Operator can be represented schematically as (see Salvadore and Baron, I952) lh'hl Q)- oeo . “iv“ @@ .69 CD +003). (8.2-6) ®.@ Q) Application of the convergence theorem stated above shows that if” = 5% > I, so that Seidel iteration of the resulting simultaneous equation cannot be used. We must solve the equation by a direct elimination method. This imposes a severe limitation on the finite- difference solution of the biharmonic equation, since electronic digital computers are limited in the number of linear simultaneous equations they can solve by this method. The MISTIC can solve a maximum of 39 equations; this results in a square grid of six by six, which is certainly too coarse for many practical problems. A possible alternative 86 is to use relaxation methods in which convergence can be forced by the skill of the operator. See, for example, Southwell et al,(I9“5). The shortcoming of this approach is the prohibitive amount of effort required for each solution, since it is not readily adapted to auto- matic solution by a digital computer for systems not amenable to other iteration methods. 990i 023:: 0 an \ .1 3 i. LY s E R A R m L V .H s R E W N U E T A“ HIIZ iiliiiij