.. 4““ a I 9 . . .«fifl‘ =22+. fiv f . :rxzotf. fivfifin '- 3. .1. H»... u»; 13.3...2n 3.... a ‘ «earnifnnc .l‘k. t9 ~x.l,..\é+.su.¢%. xanxxik.‘ r r! x sax , 11:: ‘i l :3; .. p 303‘ It: t. y i i 1.3 7 2. «RS». :r... IT. + i... w 51$ ,._ , ' '1 9.. 3 " V ' £4 , . 931‘ ‘d EX... «Sit-3: 15'! “by .32).): ‘ .(yi : 71,, .3: 1er ,a ..\ ,.....|rv v.1. rtv‘rrsi. ; $12.} r 17;” .1 You via-u: 2.1!. I, , Ezhfl whqfikfihuz . «1" r»r.,wzw:.~ .2, .411» 1,12. :7: ..»w. ‘iiaawb‘vq.§v.1< a... ..v.yn.'zy.: , 411‘ V , t4 haw «air. Vt Vtrx‘mw. V I QJu4QJ5'fVofithVV1-ia. \ c. a $51.... $.33...» huh m5. 8. , .u t. rxmwurgfi :3 L. :1 v. . l , xv. 9.. I.“ ‘s....m . a x. s. ht 2ft ‘ ‘ E3. 4., 3&0: t, . Y. . r p ‘ y? C.» . 131%.. $1“ ban». in: firui... ¢. 1.9623373. #2. vx .a 2:335 , 3.1. . 1:“ 5‘13 2:. 23:! 91.5.1»? x .1. V .,e...».,i..r.i.:. .3112: .. I}? a, 59’: ELK-13.7. It!) ¥1IVk $31; 92,; :5 , .. r? «to. . ’- .’!sb.n1.. ~I . I b I 1‘: MMC|iv¢111Iru’v 31.1.... I]. a )1, .. 14...}..3d . ‘ n...x.:zi-_ ~ rl%r 1.916 ~ 3.3;: ,7? iii. .23 L;:....;~r..\vfia!..u2 .v... 1:21.; 3.3335“! 1...! 3.1.3.3. iii its... v A 1 I 33.3.: v. {a .; ‘Vft‘wlét ’vx IF: . ’«kv\ £231 .11: . a v . u: .11.!» ».»u‘». any? Ir: 7 k . 4...»..vnvu.‘ L. .61 l, , ‘4 .1 . i: . .2’1Vvvatliint1 31‘...qmlv‘ vafloiv‘ci! 1.. 3:41.! .. . vs: .1...‘ (:31: h...n..u..fi$f.. 1.). Jfign}fl«:i1$£€i This is to certify that the thesis entitled ONE-DIMENSIONAL ELASTIC WAVE PROPAGATION IN A BAR WITH THERMALLY-INDUCED LONGITUDINAL INHOMOGENEITY ’ presented by Rajinder K. Khetarpal has been accepted towards fulfillment of the requirements for Ph.D. degree in Mechanics Department of Metallurgy, Mechanics and Materials Science \ / .- ' ”f f/ , , :‘\, h. - u, - B I k .--'7,.!'/ __;‘ \ f L/t CC‘L‘VL—Lr’trfiv‘cii‘: L J ' Rte/tap 4%. t Major professor Lawrence E. Malvern '\ Date March 283 1966 1| 0-169 ABSTRACT ONE-DIMENSIONAL ELASTIC WAVE PROPAGATION IN A BAR WITH THERMALLY-INDUCED LONGITUDINAL INHOMOGENEITY by Rajinder K. Khetarpal A temperature gradient along the length of an elastic bar gives rise to a variation in the modulus of elasticity E and the density p along the length of the bar. Hence the longitudinal elastic wave speed c = (ET/T; becomes a function of the axial coordinate. When the thermal gradient and the dependence of E and p on temperature are known, the problem becomes one of wave propagation along an inhomogeneous elastic bar with known values of the variable wave Speed C(x). The change in density due to a change in temperature has a small effect on the wave Speed for the cases considered, the primary effect being produced by the change in the modulus, which may be reduced as much as 20 percent in steel when the temperature is raised to 12000 F. Numerical integration by the method of characteristics was programmed and used on the Control Data Corporation 3600 computer. The numerical solutions for two cases were compared with the experimental data obtained: (1) for a stainless steel bar with a temperature gradient at the middle, and (Z) for a stainless steel bar with the loaded end hot and the other cold; both bars were loaded at one end by a flat-tOpped stress pulse. Experimental records of the transmitted pulses in both cases agreed very closely with the numerical solutions. For the bar heated at the middle a reflected RAJINDER K. KI-IETAR PAL pulse of strain with a magnitude of 5% of the incident pulse was observed experimentally as compared with 7% predicted by the numerical solution. The transmitted strain pulse, after passing through the gradient of rising temperature, showed at the point of highest temperature a continuous slow rise in amplitude after the initial jump, while at a second room-temperature gage beyond the hot region the amplitude was nearly the same as in the incident pulse. The transmitted pulse in the bar with one end hot was slightly different and had a higher magnitude than the incident pulse recorded on the striker bar. An analytical study was also made of the periodic longitudinal vibrations in a free-free bar with one end hot and the other end cold, excited at the hot end in the first case by a sinusoidally varying displacement and in the second case by a sinusoidally varying stress. An explicit solution is possible when the elastic modulus E is a linear, exponential, or power function of the position in the bar. For a numerical solution, programs were written to solve the periodic vibration problem by finite-difference methods for an arbitrary temperature distribution. For a 20-inch-long Type 303 stainless steel bar, the effect on the periodic longitudinal vibration was studied for several temperature distributions varying from room temperature to 1200017‘ for excitation frequencies of 5000; 7500; 10, 000; 12, 500; 15, 000; 17, 500; and 20, 000 cps. The numerical solutions agreed with the explicit solutions up to six significant figures in the cases where explicit solutions were available. ONE-DIMENSIONAL ELASTIC WAVE PROPAGATION IN A BAR WITH THERMALLY-INDUCED LONGITUDINAL INHOMOGENEITY by R ajinder K. Khetarpal A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY in Mechanics Department of Metallurgy, Mechanics and Materials Science 1966 ACKNOWLEDGMENTS Iwish to express my sincere appreciation to Professor Lawrence E. Malvern for his suggestion of this problem, and for his invaluable aid, guidance, and counsel throughout this research. For his aid in the planning and setting up of the apparatus and his help to me in so many other ways I shall always be grateful. Sincere thanks are also due to Dr. Terry Triffet, Dr. George E. Mase, and Dr. Norman L. Hills for their services on my guidance committee. I wish to also thank Mr. F] T. Bromley and Mr. Donald Childs for their help in mounting the strain gages and thermocouples; and special thanks are given to Mr. Robert B. Engle for his help in planning the apparatus and in the experimental work. I am grateful to Mrs. Barbara Judge for her careful typing of the manuscript. This project was supported by the National Science Foundation under Grant No. G-24898 through the Division of Engineering Research at Michigan State University. For my parents in India whose personal sacrifice and trust has made my education possible I will always have the greatest love and respect. I am also thankful to my perfect wife for her help and assistance throughout this project. ii 3% .- .ufdfiuumfiw t. In- 1,5154 an... Joining. a filflh ._ . . TABLE OF CONTENTS LIST OF FIGURES. . . . . . ................ NOMENCLATURE O O O C ..... O O O O O O O O O O 0 O 0 CHAPTER I. 1.1. 102. CHAPTER II. CHAPTER III. 3.1. 3.2. 3.3 3.4. 3.5. CHAPTER IV. 4.1. 4.2. 4.3. 4. 4. INTRODUCTION............... Purpose and Objectives ......... . . . Background .................. DERIVATION OF THE WAVE EQUATION. . PULSE PROPAGATION. . . . . . . . . . . Solution by the Method of Characteristics . . Numerical Procedure . . . . . . . . . . . . . (a) Leading Wave Front. . . . . . . . . (b) Impacted End . . . . . . . ....... (c) General Interior Point. . . . . . . . . Computation of c and the First Characteristic Calculation Procedure. . . . . . . . . . . . . (a) Impact at the Hot End of the Bar . . . . (b) Hot Region in the Middle of the Bar . . . Description of the Problem ..... . . . . . PERIODIC VIBRATIONS . . . . . . . . . . Introduction . . . . . . . ........... Analytical Solutions . . . . . . . . . . . . . . (a) E(x)=E0+kx.............. i. Displacement Wave Equation .. . . . ii. Stress Wave Equation . . . . . . . . (b)E(x)=Eekx i. Dispfacement Wave Equation . . . . ii. Stress Wave nEquation . . ...... (c)E(x)= EO(x/kL)r1 ............ i. Displacement Wave Equation . . . . ii. SpecialCasefornzZ . . . . . . . . (d) E(x)=E aconstant. . . . . . . . . . . Numerical Solution. . . . . . . . . . . . . . . (a) Formulation of Algebraic Equations. . . (b) Solution of the Algebraic Equations . . . i. Direct Solution. . . . . . . . . . . . ii. Iteration............... Description of the Problem . . . . ...... iii l3 13 15 15 17 18 2.0 22 Z3 Z4 Z4 Z9 29 3O 30 3O 34 37 37 41 44 44 50 53 57 57 63 64 7O 71 CHAPTER V. EXPERIMENTAL APPARATUS . . . . . . . 5.1. General Description . ...... . ...... 5.2. Specimen and Striker Bar . . . . . . . . . . . 5. 3. Hyge Shock Tester. . . . . . . . ...... 5. 4. Furnace and Temperature Measurement . . . 5. 5. Strain Measuring and Recording Equipment. . CHAPTER VI. EXPERIMENTAL PROCEDURE ...... 6212 Determination of E, p and c ......... 6. 2. Electronic Calibration of the Strain Gages . . 6.3. Test Procedure . . .............. 6.4. Reduction of Data ............ . . . CHAPTER VII. RESULTS. . . ............... 7.1. Pulse Propagation ...... . ..... . . . (a) Oscilloscope Records . . . . . . . . . . (b) Discussion of the Test Results ..... i. Results from the 6-Foot-Long . Bar Heated in the Middle . . . . . . 11. Results from the 4-Foot-Long Bar Impacted at the Hot End. . . .. . . . (c) NumericalResults. . . . . . . . . . . . 1. Numerical Solution Results from the 6— Foot- -Long Bar Heated in the Middle. . . .. ..... ii. Results from the 4- Foot- Long Bar Struck at the Hot End. . . . . . iii. Results of Lindholm's Problem and the Modified Lindholm Problem. 7.2. Periodic Vibrations . . . .......... (a) The Results of Periodic Vibration Solutions ...... . . . ........ (b) Discussion of the Results .. ........ CHAPTER VIII. SUMMARY AND CONCLUSIONS ..... BIBLIOGRAPHY ................. . . . . . . FIGURES . . ......................... APPENDIX A. INTERPOLATION ....... . ...... A.l.Introduction... A. 2. Aitken' 3 Method of Interpolation . . . . . . . (a) Linear Interpolation . . . . . ...... (b) Aitken's Repeated Process ..... . . (c) Programming.............. APPENDIX B. COMPUTER PROGRAMS ...... . . . . B.1. Pulse Propagation. . . . . . . . . ...... B. 2. Longitudinal Periodic Vibrations . . . . . . . B. 3. Fortran Programs .............. iv Page 73 73 73 76 77 79 83 83 84 86 87 88 88 88 88 88 89 9O 9O 91 93 95 95 96 98 104 106 149 149 150 150 151 154 156 156 162 166 to 148 .., ,1 W I11. , is. unfailinxmt. .. .1. I 1!. H.414! .1: .1 1; ul .1 I I ”In 1 U: 4 I“. G. .151. a. . . . . f H 1 . a. J... . .a.‘ ... . .. a... .u. I . n. . 4...... MU... {.1 wagilfibh . ,1- J. n. In: 11?: . u w n... a F4: ,1. ....,1 3.10 3.11 3.12 3.13 5. 3(a) 5. 3(b) 5. 6(a) 5. 6(b) LIST OF FIGURES The Characteristics in the XT-Plane. Temperature Distribution for the 6-Foot-Long Bar Heated at the Middle . Temperature Distribution for the 4—Foot-Long Bar Heated at the End. Elastic Modulus Versus Temperature Wave Speed Versus Distance for the 6-Foot Bar Heated at the Middle . Wave Speed Versus Distance for the 4-Foot Bar Heated at the End. Leading Wave Front Characteristic for the 6-Foot Bar Heated at the Middle . Leading Wave Front Characteristic for the 4-Foot Bar Heated at the End. Types of Temperature Distribution Considered Characteristic Field for the Bar with a Hot End . Characteristic Field for a Bar Heated in the Middle . Incident Pulse for the 6- Foot- L—ong Bar Heated intheMiddle Incident Pulse for the 4- Foot— -Long Bar Heated attheEnd. Schematic Drawing of the Apparatus . A General View of the Test Set-Up . Details of the Bars for the First Experiment. Details of the Bars for the Second Experiment. Hyge Shock Tester (after Chiddister, 1961) Close- -Up of the Furnace and the High- Temperature Gage. . WheatstoneBridge................ Potentiometer Circuit . ......... . ..... Page 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 121 122 123 124 124 .10 .11 .12 .13 .14 .15 .16 .17 .18 .19 .20 .21 Oscillosc0pe Records from the 6-Foot-Long Bar HeatedintheMiddle .. . . . . . . . . . Oscilloscope Records from the 4-Foot-Long Bar Impacted at the Hot End . Pulse at Various Points of the 6-Foot Bar Heated at the Middl e Stress Pulse at Various Points of the 6-Foot Bar Heated at the Middle . The Transmitted Pulses in the 4-Foot-Long Bar Impacted at the Hot End . . . . . . . ..... . The Transmitted Stress Pulses in the 4-Foot-Long Bar Impacted at the Hot End . Reflected Pulse in the Modified Lindholm Problem. Stress Amplitude for 5000 cps. Displacement Amplitude for 5000 cps. . Stress Amplitude for 10, 000 cps. Displacement Amplitude for 10, 000 cps. Stress Amplitude for 15, 000 cps. Displacement Amplitude for 15, 000 cps. Stress Amplitude for 20, 000 cps. Displacement Amplitude for 20, 000 cps. Stress Amplitude for 5000 cps. Displacement Amplitude for 5000 cps. . Stress Amplitude for 10, 000 cps.. Displacement Amplitude for 10, 000 cps. Stress Amplitude for 15, 000 cps.. Displacement Amplitude for 15, 000 cps. Page 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 .22 .23 .24 . 1 StressAmplitude for 20, 000 cps. Displacement Amplitude for 20, 000 cps. . Temperature versus Distance for the Assumed Distributions of E Range of Application of Interpolation Formula . vii Page 146 147 148 154 NOMENCLATURE Elastic wave speed, also a defined constant in Chapter 4. Elastic wave speed at room temperature Non-dimensional elastic wave speed, 2 c(x)/cO Non-dimensional elastic modulus, = E(x)/EO Constants in Chapter 4 Elastic modulus Elastic modulus at room temperature Bessel's function of the first kind and order p A constant, with units of (sec. )-1 in Chapter 3 Also other defined constants in Chapter 4 Length of the bar A defined constant Amplitude of the stress end condition Non-dimensional stress, = cr/EO Time Non-dimensional time, = t/k Longitudinal displacement, considered positive when in the negative x-direction Amplitude of the displacement end condition Particle velocity, considered positive when in the negative x-direction Non—dimensionalized particle velocity, = v/cO Initial coordinate of cross section Non-dimensional distance, = kx/cO Non-dimensional strain, = 6 viii Bessel's function of the second kind and order p Strain, considered positive in compression Density Stress, considered positive when compressive Circular frequency ix CHAPTER I INTR ODUC TION 1. 1. Purpose and Objectives It is well known that the mechanical properties of a material change with the temperature. Among other changes in mechanical pr0perties due to a change in temperature, the changes in the modulus of elasticity and the density affect the prOpagation of a wave. The purpose of the present investigation is to study the effect of a thermally-induced inhomogeneity on the pr0pagation of a one-dimensional elastic disturbance. The problems investigated in this study have the following possible applications: (a) Interpretation of pressure-bar records from hot or radioactive environments (see for example Chiddister*, 1961), (b) Finding the critical frequencies and solving resonance problems of bars subjected to varying temperatures, (c) Calibration of high-temperature strain gages, and (d) Obtaining the value of the elastic modulus at elevated temperatures. A temperature gradient along the length of an elastic bar gives rise to a variation in the elastic modulus E and the density p along the length of the bar. Hence the longitudinal elastic wave speed c =N} 3% in the bar becomes a function of the axial coordinate. \'a '|‘ Surnames followed by dates refer to Bibliography. When the thermal gradient and the dependence of E and p on temperature are known, the problem becomes one of the wave pr0pagation along an inhomogeneous elastic bar with known values of the variable wave speed C(x). The change in the density due to a change in temperature is quite small. Based on the value of the coefficient of expansion for Type 303 steel given in the Metals Handbook (1948), the difference in the values of the density at 75°F and 1200°F is 3. 89% and the resulting difference in the values of the elastic wave speed at IZOOOF is only 1. 98% (the values are given in Sec. 6.1). Therefore the density can be considered constant, and the primary effect is produced by the change in the elastic modulus, which may be reduced as much as 20% in steel when the temperature is raised to IZOOOF. The inhomogeneity in this study is therefore prescribed as a variation in magnitude of the elastic modulus only. The present investigation had the following objectives: (a) To study the effect of the thermally-induced inhomogeneity on the propagation of a pulse in an elastic slender bar with an arbitrary temperature distribution along its length. (b) To study the effect of the thermally induced inhomogeneity on the periodic vibrations "of an elastic slender bar with an arbitrary temperature distribution along its length and various boundary conditions. (c) To set up an experiment to study the effect of the thermally-induced inhomogeneity on the pr0pagation of a' pulse produced by an impact in an elastic slender bar with a temperature distribution along its length, and compare it with the outcome of the objective (a). Numerical integration by the method of characteristics was programmed and used on the Control Data Corporation 3600 computer. The program was used (11 five Specific problems, calculating results: 1. to compare with the experimental pulse-propagation study on the bar heated at its center; to compare with the experimental pulse-prOpagation study on the bar heated at one end; to verify the results of Chiddister (1961), which he obtained by approximating an experimentally-obtained thermal gradient by a series of steps; to verify the results of Lindholm (1963), which he obtained by a series method for the Special case x n E 2 E0 (1:13) ; and to solve a modified version of Lindholm's problem with variation according to a cubic power law in one-fourth the bar, while the rest of the bar had a uniform elastic modulus. Two experimental pulse propagation studies were made on Type 303 stainless steel test bars: 1. Bar Heated at the Middle. A 6-foot-long bar was heated at the middle to a maximum temperature of 12000F by a 5-inch-long coaxial furnace and impacted longitudinally at one end by a 4-foot-1ong striker bar. The pulse produced in the test bar was recorded by three strain gages, two at room temperature and one high-temperature gage at the middle of the bar. The first room-temperature gage recorded the incident pulse and the reflection from the thermal gradient, while the second recorded the transmitted pulse in the region beyond the heated middle portion of the bar. 2. Bar Heated at the Impact End. A 4—foot-1ong bar was heated at one end by the same coaxial furnace and impacted as in the first experiment. In this second experiment one of the room- temperature gages was on the striker bar. Results in both experiments agreed very closely with the numerical solutions. In the case of the bar heated at the middle, the maximum reflected pulse was 7% of the incident pulse according to the calculations and 5% from the experiment. The initial jump and the amplitude of the transmitted strain pulse recorded at the second room-temperature gage were nearly the same as the incident pulse. The amplitude of the pulse at the point of highest temperature in the bar heated at the middle showed a continuous slow rise after the initial jump. The transmitted pulse in the bar struck at the hot end was slightly different and had a higher magni- tude than the pulse recorded on the striker bar. An analytical study was also made of the periodic longitudinal vibrations in a free-free bar with one end hot and the other end cold, excited at the hot end in the first case by a sinusoidally varying displacement and in the second case by a sinusoidally varying stress. An explicit solution is possible when the thermally-induced inhomo- geneity is one such that the elastic modulus E is a function of the position in the bar of the form E : E0 + kx Datta (1956) or E : E0 ekx Sur (1961) x n or E : Eo (kl-3) Lindholm (1963) These problems have been solved by converting the wave equation to a form of Bessel's equation. The solution of the displacement wave equation has been obtained by the same method as that used by the authors cited. In addition, the explicit solutions of the stress wave equation with the same type of inhomogeneities as in the case of the displacement wave equation have been obtained in the present investigation. In the. present investigation programs were written to solve the periodic vibration problem numerically by finite differences for an arbitrary temperature distribution. These pr0grams were used to determine how certain thermal gradients affect the periodic longitudinal vibration of a 20-inch-1ong Type 303 stainless steel bar. For the excitation frequencies of 5000, 7500, 10, 000, 12, 500, 15, 000, 17, 500, and 20, 000 cps, explicit solutions of the stress and displacement were obtained for the thermally-induced inhomogeneity in which the elastic modulus varies with the distance as Eo + kx, Eo ekX and E0 (Exi)2 ; and numerical solutions of the stress and displacement were obtained for the same three cases as a check and in addition for variation as Eo (fig/2 and for atemperature distribution measured in the laboratory under conditions similar to those of the pulse-propagation experiments. The value of k in each case was such that one end of the bar was at 750F and the other end at 12000F. With these choices of the constants the exponential and power function variations of E(x) do not depart much from the linear variation in the 20-inch bar, with the result that those solutions are also very close to the solution for linear variation. The numerical solutions agreed with the explicit solutions up to six significant figures in the cases where explicit solutions were obtained. The solutions were also obtained for the same bar with a uniform temperature of 750F in one case and 12000F in the other. Plotting all the above solutions, it was found that the solutions, for the cases of varying temperatures in the bar, were always between the solutions obtained for the uniform temperatures of 750F and 1200OF. The presence and the location of the nodes, the critical frequencies, and the amplitude, were all affected by the inhomogeneity in the bar. 1. 2. Background A considerable amount of work has been done on inhomogeneity (see Olszak, 1959). Sternberg and Chakaravorty (1959) have considered the problem of the prOpagation of a shear shock wave from a circular hole in a semi-infinite plate, where the hole is subjected to suddenly rising uniform shearing tractions and the shear modulus of the material of the plate is pr0portional to an arbitrary power of the radial distance from the center of the hole. Cristescu (1959) has considered the prepagation of waves in elastic-plastic non-homogeneous thin and semi-infinite rods. In the elastic region he took 0' : E(x)€ . Perzyna (1959) also considered the propagation of elastic-plastic waves in a non-homogeneous bar for a general type of non- homogeneity. He analyzed the problem with o~ : f(€, x). Both Cristescu and Perzyna formulated the problem by the method of characteristics in the same manner as will be done in this investigation, but they did not include any numerical solutions. Datta (1956) has studied the propagation of sinusoidal and impulsive disturbances in a bar having linear variation of the elastic parameters; he obtained the solution for the impulsive loading by Laplace transform techniques. Ghosh (1961) has studied the problem of extensional vibration of a bar having linear variation of the elastic parameters, and excited by the impact of an elastic load; he used Operational methods. Sur (1961) has studied the pr0pagation of sinusoidal and impulsive disturbances in a bar having exponential variation of the elastic parameters, and he too obtained the solution for the impulsive loading by Laplace transform techniques. Chiddister (1961) studied the problem of pulse propagation in a bar with thermally-induced longitudinal inhomogeneity by approximating the thermal gradient by five step-discontinuities separating regions of constant temperature and superposing the waves calculated by simple reflection theory. Chiddister used his approximate solution to interpret the pres sure-bar records from the hot environments where his specimen was located. As was mentioned above, this is one of the important applications of the present investigation. Lindholm (1963) solved the problem of an elastic disturbance propagating in a nonhomogeneous bar of finite length by using the . 51x11 .‘flwgll .. . a. 3.4.... £11 1% . . r .. . . .111 .U (.11.. lduljill . 31.37.?” t 4 principle of virtual work. He prescribed the nonhomogeneity as a modulus of elasticity continuously varying with position in the bar, given by E : Eo(fi)n, and solved the pulse pr0pagation problem by a superposition of sinusoidal solutions at different frequencies. In his investigation of the pulse propagation he found an apparent absence of reflections. An explanation of this apparent difference from the results of the present investigation is given in Sec.7. 1(c. iii). CHAPTER II DERIVATION OF THE WAVE EQUATION Three different forms of the governing equations of wave propagation will be derived here: a set of two first-order partial differential equations for the unknown stress and particle velocity; a second—order wave equation for the disPIacement; and a second- order wave equation for the stress (see Timoshenko, 1955). These governing equations of motion for the one-dimensional theory of longitudinal wave propagation in a bar of uniform cross- section, with the elastic modulus varying along the length of the bar, are obtained by assuming that plane sections remain plane and that across these‘sections the stress is uniform, and that the displace- ments are small. At a time t, let u(x, t) be the displacement of the cross- section initially at a distance x from the left-hand end of the bar. At the section under consideration the strain E and the particle velocity v are then _ 22 e — 8x (2.1) V : its (2.2) In this discussion the stress and strain are considered positive when they are compressive, and the disPlacement and the particle velocity are considered positive when they are in the negative x direction. 10 By differentiating Eq. (2.1) with respect to time and Eq. (2. 2) with respect to x, and eliminating u, we get the equation of continuity 0.9 m e .531 (2.3) Q) q Q) < — = P(X)— (2.4) E(x) E = 0' or, 86 80" E(x) 5'? = '5'? (Z. 5) Substituting Eq. (2. 5) in Eq. (2. 3), 8 8 x) at q < 1H Q) E X A We therefore have a system of two linear first-order partial differential equations. The system is written as 1 E 2.1 — 0 1 E(x) 8t - 8x — 1 (2.6) 80 8v _ 5;; 'MX)??? -0 or, 11 From equations (2.4) and (2. 2) 2 8 8 5% : p (x) ——g (2. 8) at so that, 8 811 azu 5;{E(X)_8-£} = P(X)g—tz— (2-9) Eq. (2. 9) is the one-dimensional wave equation. Perzyna (1959) gave a similar derivation for the wave equation, Eq. (14) in his paper. Lindholm (1963) begins with the same wave equation. From equations (2. 3) and (2. 5) 80' _ 8v 5‘? - W) 5:; so that 1 820 2 82V (2 10) E(x) atZ atax Differentiating Eq. (2. 4) with respect to x yields a 1 ae 82v 532 (p(x) 5;) : axat (2.11) From equations (2.10) and (2.11) a 1 Be 820 E(x) 53: (pee a—x') = M2 (2°12) If the variation in density is negligible, we get the following stres 5 wave equation 12 2 2 8 8 E(x)——92- = p g (2.13) 8x 8t 2 2 or, c2(x) 8 g = 8 : (2°14) 3x at where, c(x) = 1251‘) (2.15) c(x) is the elastic bar -wave speed. The system of two first-order equations will be used in the third chapter to study pulse propagation by the method of charac- teristics, while the two second-order wave equations will be used in the fourth chapter to study the sinusoidal stress and displacement standing waves. CHAPTER III PULSE PR OPAGATION 3. 1. Solution by the Method of Characteristics We have obtained a system of two linear first-order differ- ential equations, Eq. (2. 6), for the pulse pr0pagation in a bar. 1 __1 9.1 21—0 Ex) 8t - 8x — f (3.1) 30' 3v _ '53: - P1X) _—t - 0 J Since the system is hyperbolic, it is suitable for numerical solution by the method of characteristics (Courant and Hilbert, 1962). i 1 For Eq. (3.1), we have 2( §_u__ +b .3.) z 0 1 aki ax ki a 1 i_ 0' _ O -l — ——E(x) 0 u ‘ ' aki ‘ ’ bki ‘ u 1 0 O -p(x) And the characteristic curves are defined by dx = + c(x) dt (3. 2) dx - c(x) dt where c(x) = 12(5):) is the speed of propagation of the elastic wave front. The interior differential equations holding along the charac- teristics correSponding to Eq. (3. 2) are 13 14 I l 0 d0 - p(x) c(x) dv (3.3) (10 + p(x) c(x) dv = 0 Both Cristescu (1959) and Perzyna (1959) give the equations defining the characteristics, and also interior differential equations relating v and 6 for the elastic-plastic problem, but the papers cited included no actual numerical solutions. For most of the metals, e. g. steel, the coefficient of expansion per unit volume is quite small. And because density is mass per unit volume, the change in density due to the temperature change is also small (see Sec. 6.1) and will be neglected in this investigation. The system of Eqs. (3. 2) and (3.3) can be treated by a numerical procedure. The variables are first non-dimensionalized by the following transformation. __<:_ _1<_x S_E X—C o o _ _E(X) Y—e D.E o V_V_ ) "C (3:21):— o c o Tzkt where E0 and c0 are the elastic modulus and the wave propagation 2 speed at room temperature; (E0 = p cO ); k has the units of sec. - , and its value is chosen for convenience. The transformed characteristic curves and the interior differential equations along them can be written as 15 dS-CdV ‘l 0 along the curve dX = CdT ! l .4 t((3 ) dS + C dV 0 along the curve dX = - CdT For numerical calculations these equations are treated as finite -difference equations, with dS and dV replaced by A S and A V, respectively and C replaced by an average value on the characteristic segment. The XT —plane is subdivided by a mesh of characteristic curves at finite intervals of A T and A X. The characteristics are plotted with a constant interval of A T, and A X is allowed to vary. As shown in Fig. 3.1, for X > 0 there are two characteristics passing through each mesh point. Therefore at any mesh point P of the XT-plane'the solution can be obtained by solving two difference equations along the apprOpriate characteristics, if we already have the solution for the two neighboring points A and B. For the points along the edges of the bar, one of the unknowns is prescribed, so that one equation will be sufficient to give the solution. Further details of the numerical solution are given in the next section. 3. 2. Numerical Procedure We write the interior differential equations, Eq. (3. 4) as difference equations, with the average value of C along a segment approximated by the arithmetic mean of its two end values, and obtain 16 (0 +0) (SP'SA)' 132 A (VP'VA):O \: 1, (3.5) (c +C) (SP-SB)+ P2 B (VP-VB):O f These equations are solved simultaneously to give the values of the dimensionless stress S and velocity V at the point P of Figure 3.1 in terms of the values of S and V at A and B. Proceeding point-by-point in the mesh and using Eq. (3. 5), we obtain the values of S and V for all the points in the field of characteristics. The characteristic field can be constructed by solving the characteristic difference equations, based on Eq. (3. 2), which can be done Without solving the pr0pagation problem. The calculations with Eq. (3. 5) then give the solution at the mesh points. There are three possible types of points to be considered for a wave pr0pagation in a semi -infinite bar: Leading Wave Front; Impacted End; and General Interior Point. For a finite bar the points at the other end also require special treatment. 3. 2(a). Leading Wave Front Along the leading wave front the values of the stress, strain, and the velocity are zero for a continuously rising pulse. For a discontinuous jump, the values after the jump are obtained by equating impulse to the change of momentum and by continuity to give the jump c onditions 17 AS:-CAV AV=-CAY (am AS=C2AY 3. 2(b). Impacted End Only one equation is available for the impacted end. Assuming that one boundary condition is known at this end, the solution can be obtained by considering the propagation along dX : - CdT. Two possible cases are (a) The Stress Boundary Condition is Known: Let the stress at the point X: 0 be SO(T). Then Eq. (3. 5) can be written as (0 +0) (SP-SB)+ P 2 B (VP-VB) = o (3xn The solution of Eq. (3. 7) can be written as SP : so YP : :59 (3.8) P . 2 1 VPZVB'EETC—Bwo'sB) j (b) The Strain Boundary Condition is Known: Let the strain boundary condition be Y(0,'I) : YO(T). The equations for this boundary condition are S 18 YP = YO SP = DPYP )(3.9) (c + c ) 9 P B l (SP-SB)+ 2 (VP-VB) _ 0/: The solution of Eq. (3. 9) can be written as \ SP 2 DP Yo YP = YO >(3.10) VP:VB'C 2c: (DPY 'SB) P + B 0 1 (c) General Interior Point For a general interior point the Eq. (3. 5) can be written as (0 +0) \ (SP-SA)-—P—Z-—A (VP-VA) = 0 ‘ f )(3.11) (cp+cB) (SP-SEN 2 (VP-VB) = o) The solution of Eq. (3.11) can be written as _—_ cp+cB s + cl:,+cA -s 2CP+CA+CB A ch+cA+cB B (3.12) + ((280 Cf; (Cf; ()SB) (VB ' VA) P A B YP =‘SD'E P v 1 P : 2CP+CA+CB { (C P+CA) VA + (CP+CB) VB ‘ZSA+ZSB} 19 The equations, Eq. (3.8), Eq. (3.10) and Eq. (3.12), have been programmed in Fortran for the CDC-3600 computer. The program is given in Appendix B._ Two of the possible boundary conditions at the other end of the bar (opposite from the impact end) are free end and fixed end. At a free. end, stress is equal to zero, so that at all the mesh points falling at that end of the bar the non-dimensional stress S will be zero. At a fixed end, particle velocity is equal to zero, so that at all the mesh points falling at the end of the bar the non-dimensional particle velocity V will be zero._ In both cases only one of the Eq. 3 (3. 5) is used to obtain a solution. The program has been used to solve five different problems: 1. Pulse pr0pagation in a 6-foot stainless steel bar with the experimental temperature distribution in the middle. The maximum temperature was 12000F. The measured temperature distribution is given in Fig. 3. 2. Temperature was measured at 24 points 0. 8 in. apart on one side of the temperature distribution. For the other side, the temperature distribution was taken to be a Inirror image of this side. 2. Pulse pr0pagation in a 4-foot stainless steel bar with the experimental temperature distribution such that the impacted end was at 12OOOF. The measured temperature distribution is given in Fig. 3. 3. Temperature was measured at 20 points one inch apart. 3. Verification of the results obtained by Chiddister (1961). 20 4. Verification of the results obtained by Lindholm (1963). 5. Solution of a modified version of Lindholm's problem with variation according to a cubic power law in one -fourth of the bar, while the rest of the bar had a uniform elastic mo dulus . Each mesh point was designated as (i, j) where i is the number of the characteristic with positive slope through the mesh point; for example, i=1 identifies the leading wave front X = CT; and j is the number of mesh point on the characteristic (see Fig. 3.1). In the first problem the final mesh had 900 horizontal intervals, spaced so that the time At taken by the wave to travel across each interval was the same, namely At equal to 0. 4104 microseconds. In the second problem the final mesh had 960 horizontal intervals each covered-by the wave in time At equal to O. 2564 microseconds. For the third and fourth problems 100 horizontal intervals were used and for the fifth problem 160 intervals. 3. 3. Computation of c and the First Characteristic The elastic wave speed c(x) is equal toJE. Garafalo (1960) found that the elastic modulus of 18-8 stainless steil could be represented as a linear function of temperature. The slope of his curve of modulus versus temperature was used to construct the linear plot given in Fig. 3. 4 for the Type 303', 18-8 stainless steel Specimen bar used in the present study. The room temperature value was determined experimentally and used with the slope from Garafalo's paper to construct the plot. At 750F the experimentally determined value was 21 E0 = 29.0x106 psi., while from Garafalo's paper SlOpe = - 6'. 47 x 103 psi/OF. From the values of the temperature measured at various points asmentioned in Section 3. 2, the values of the elastic modulus were obtained from the linear relation illustrated in Fig. 3. 4 at those points. These in turn were used to get the elastic wave speed versus distance curves shown in Fig. 3. 5 and Fig. 3. 6. The procedure of Sections 3.1 and 3. 2 leads to the values of S, Y, and V at mesh points of the characteristic net. To make this into a solution in terms of X and T, it is necessary to determine the characteristic net, i. e. to find values of X and T at each mesh point. For this we find a curve between X and T, which will be the leading-wave front characteristic curve. Other characteristic curves of the same family are parallel to it at intervals of ZAT above it. To obtain the leading wave front we start from the equation of the characteristic dXZCdT (3.13) _ X 1_ or T—fo C dX Eq. (3. 13) is integrated by Simpson's rule, using the numerical data cited above in Figures 3. 5 and 3. 6 (non-dimensionalized). The interval from 0 to X is divided into 2m equal intervals by points 15‘ aw... 22 X0, X1, X2, . . . , sz, and the integral is evaluated by the following numerical Simpson's rule. (See for example Milne (1949), page 121.) X f 2 X0 m h ' de—§{yo+y2m+4(yl+y3+... JrYZm-l)+ 2 (y2 +y4+ ... +y2m_2)} In Fig. 3. 7 and Fig. 3. 8 curves between dimensional values of x and t are shown for the leading wave front. These were obtained by integrating the dimensionalized form of Eq. (3. 13). The reason for the equal x intervals here is that the experimental data on temperature was taken at points spaced at equal distances. The mesh-point values for equal time intervals were subsequently obtained by interpolation. At first glance these curves may appear straight, but a closer look will disclose definite deviations from linearity in the heated region, around 30 inches from the end in Fig. 3. 7, and near the hot end in Fig. 3. 8. 3. 4. Calculation Procedure Three types of temperature distribution conditions are of interest for pulse pr0pagation in a bar subjected to a longitudinal impact. 1. One end hot and the other cold. Impact at the hot end. 2. One end hot and the other cold. Impact at the cold end. 3. Both ends cold. Hot region in the middle of the bar. The numerical and experimental procedures for the first two cases are similar, so only the first and the third cases are considered here. 23 3. 4(a). Impact at the Hot End of the Bar Fig. 3. 9(a) Assuming that the shape and the magnitude of the incident pulse are known, the field of characteristics is calculated first, choosing a suitable time interval. The time interval was considered sufficiently small if a solution obtained by choosing the grid still smaller only changes the sixth digit of the solution. Along the leading wave front the values of the stress, strain, and the velocity are zero for a continuously rising pulse. The procedure is illustrated by the schematic drawing of Fig. 3.10. The SlOpe of the characteristics varies in the region where the temperature varies, and is constant where the temperature is uniform. For a finite bar, for which we have to consider the reflections from the ends, the total time taken for the wave front to travel from one end to the other end of the bar is calculated from the leading wave -front XT-curve. Then a time interval of the grid is calculated by dividing the total time by the number of divisions estimated to be large enough for sufficient accuracy. The grid Spacing can also be manipulated so that the grid point falls on the bar at a point where a strain gage is mounted in the experimental measurements, in order to obtain directly the calculated values without plotting extra curves. Following the procedure mentioned in Sec. 3. 2 and solving Eq. (3. 8), Eq. (3. 10), and Eq. (3. 12), the value for non-dimensional stress and velocity can be calculated at all the points of the grid. 24 3. 4(b). Hot Region in the Middle of the Bar This is the case in which a traveling pulse passes through a temperature gradient like that of Fig. 3. 9(b). Again assuming the incident pulse is known, the field of characteristics is calculated first with a suitable time interval. The field of characteristics is shown schematically in Fig. (3.11). Again characteristics have varying slope in the region of varying temperature and become straight where temperature is uniform. Grid spacing is again calculated for convenience, accuracy, and position of the gage station in experimental measurements. Again solving Eq. (3.8), Eq. (3. 10), and Eq. (3.12), and introducing the boundary conditions, the non-dimensional value of the stress and velocity can be calculated at all the points of the grid. 1 Plotting the value of stress or strain on all the grid points on a line X equal to a constant X the stress versus time or 1’ strain versus time curve can be obtained for the point on the bar at distance X1 from the impacted end. 3. 5. Description of the Problem Two types of pulse propagation problems were solved by the numerical methods. The solutions, obtained by the method of characteristics, are given in Sec. 7. The two problems solved to compare with the experiments are 1. Bar Heated at the Middle. A 6-foot-1ong bar was heated at the middle to a maximum temperature of IZOOOF and impacted longitudinally at one end by a 4-foot-long striker bar. The 25 temperature distribution in the bar is given in Fig. 3. 2. The solutions of the pulse propagation were obtained and compared at the location of the strain gages in the experimental set-up described in Sec. 5. 2, i. e. two room-temperature gages located one foot from each end of the bar and one high-temperature gage located at the middle of» the bar. 2. Bar Heated at Impacted End. A 4-foot-long bar was heated at one end to a maximum temperature of 12000F and impacted as in the first case. The temperature distribution in the bar is given in Fig. 3. 3. Again the solutions were obtained' and compared at the location of the strain gages in the experimental set-up described in Sec. 5. 2. In the first problem the temperature was measured at 22 points spaced 0. 8 inches apart, starting from the center of the bar, and in the second problem the temperature was measured at 16 points spaced one inch apart starting from the hot end. The description of the experimental measurements is given in Sec. 5. 4. With these measured values of temperature the values of the elastic modulus were obtained from the linear relation of Fig. 3. 4 and the wave speed c(x) = (Mp calculated as plotted in Fig. 3. 5. For the first problem, the incident pulse was obtained from the output (Fig. 7. l) of the first room-temperature strain gage, located one foot from the impacted end (Fig. 5. 3) before the reflections from the thermal gradient reached it. The pulse being flat-t0pped when reflections are absent (Fig. 7. 1a), the incident pulse was assumed to be 486 microseconds long and flat-t0pped with an amplitude given by 26 the first room-temperature gage. The incident pulse obtained from the oscilloscope record of Fig. 7.1b is shown in Fig. 3.12; it was used in the solution of the problem as input data for the computer program described in Appendix B. This strain pulse at the room— temperature gage station was converted to stress by multiplication with the room-temperature elastic modulus. Because the recording station is so far from the impact end, three-dimensional effects are believed insignificant in this problem. For the second problem, the incident pulse was obtained from the output (Fig. 7. 2) of the room-temperature gage located on the striker bar (Fig. 5. 3) 25 inches from the impact end. The incident pulse was assumed to be 486 microseconds long and flat- t0pped with an amplitude obtained from the gage on the striker bar. The incident pulse obtained from the oscillosc0pe record of Fig. 7. 2b is shown in Fig. 3.13. It was used in the numerical solution of the problem as input data for the computer program described in Appendix B. This strain pulse in the room-temperature striker bar is assumed convertible to stress at the impact end of the Specimen by multiplication with the room-temperature elastic modulus. The assumption neglects three-dimensional effects in the vicinity of the impact ends of the two bars. Since three—dimensional effects are known to be present at the impact end (see Bell, 1960), some error may be introduced by the simplifying assumption used. The incident pulses shown in Fig. 3.12 and Fig. 3.13 are obtained by drawing a smooth curve through the experimentally measured points. As can be seen in Figs. 7.1 and 7. 2 there was 27 noise on the signal. Experimental points were measured only at points clearly on the pulse. The Fortran programs for numerical solution of both of the above mentioned problems have been described in Appendix B and programs are given at the end of the Appendix. The results obtained from the numerical solution of these problems have been given and discussed in Sec. 7. 3. The programs were also used for three other problems: 3. Verification of Chiddister's (1961LResults. This was the first calculation made with the program written for Problem 1 above in order to check the program. The results of this calculation showed good agreement with Chiddister's results, obtained by approximating the thermal gradient with five step discontinuities. These results are not reproduced in this dissertation. 4. Verification of Lindholm's (1963) Results. Lindholm calculated the propagation of a half-sine -wave pulse along a'bar with E = EO(EXE)n . The results are discussed in Sec. 7.1 (c.iii). 5. Modified Lindholm Problem. To demonstrate that reflections would be obtained from a steeper gradient varying as a power law, the following problem was solved. With the same values of maximum and minimum E, the length of the part of the bar where the inhomogeneity existed was reduced to one-fourth the length of the bar in Lindholm's third problem, with the variation following the same cubic law as in his problem, and E was taken uniform along the other three-fourths of the bar. The length of the pulse was again taken approximately one-fourth of the length of the bar. This gave a 28 c(x) of the form 3 103/ 2] c(x) = c0[1+.296(x-Z- in the region where E was varying. The results are given in Sec. 7.1 (c. iii). CHAPTER IV PER IO DIC VIBR ATIONS 4. 1 . Introduction The solution of the wave equation for the periodic vibrations of a bar with a thermally-induced longitudinal inhomogeneity can be obtained numerically for a general kind of inhomogeneity, and analytically for some Specific kinds of inhomogeneity. Analytic solutions of the wave equation in a bar with a thermal gradient have apparently been obtained only for cases where the inhomogeneity in the modulus of elasticity caused by the thermal gradient is of the following types: 1. E(x) 2 EO +kx (Datta, 1956) kx 2. E(x) = Eoe - (Sur, 1961) x n 3. E(x) : E0015? (Lindholm, 1963) and the variation in density is assumed negligible. The analytic solution method is discussed in Sec. 4. 2. For a general case of inhomogeneity we have to use numerical methods for the solution of the wave equation because of the absence of analytic solutions. The numerical method is described in Sec. 4. 3. Both the analytic solutions and numerical solutions are obtained for two types of boundary conditions: first a free-free bar with a source of disturbance producing periodic displacement at the hot end; 29 30 and second a free-free bar with a source of disturbance producing periodic stress at the hot end. The solutions in each case are obtained both for the displacement amplitude and the stress amplitude of the periodic vibrations. The analytic and numerical results obtained for a 20-inch-long Type 303 stainless steel bar with various types of inhomogeneity are discussed in Section 7. 4. 4. 2. Analytic Solutions 4. 2(a). E(x) 2 EO + kx 4. 2(a. 1). Displacement Wave Equation From Eq. (2. 9), neglecting variations in p, the differential equation is 2 8 au 5 u -— (E(X) —) = P —— 8x 8x at2 Substituting the value of E(x) weiget a an azu ax {(E0+kx)ax} =9 atz (4-1) To obtain a periodic solution of Eq. (4.1) let u(x, t) = u(x)eiwt The resulting ordinary differential equation for u(x) is 2 (Eo+kx)9—%+k%—:+p02=o (4.2) dx We transform this equation to a Bessel's equation, Eq. (4. 4) below, by the following substitutions. First let - ,. m .. if“. ...-..merw . 3,4: . iiinqlvawa. _ 31 z = E + kx o to obtain 2 2 zd121+gg+p£i§ =0 (4.3) dz k Then substitute s‘Zk‘2 2 N} - or S = 22% 4pw _<_1_§_ :w'x/fi-‘p Z-l/Z nd dzs _ mfg? -3/2 dz k a 2 ‘ ' 2k Z _ dz (02 to transform Eq. (4. 3), after division by ——23 , to k 2 d u 1 du __ 2+;gg+u—O (4.4) ds which is a Bessel's equation with general solution u = AJO(S) + BYO(S) (4. 5) Hence the complete solution as given by Datta (1956) is u(x, t) -— (“022%. “2) + BYO(§-“’—k—P— “Zn e1“t (4. 6) where z=E +kx 0 A and B are determined from the Boundary Conditions. Case 1 Let there be a source of disturbance producing periodic displacement ert at x = L, and let the end x = 0 be 32 free. The, boundary conditions are __ iwt 11]sz Ue au _ -8—X]X:O _ O From Eq. (4. 6), 22-_ _ k [A{Zb)'\/—we J- (ZU‘VEZ)} 8x _ ZJX k l k (4. 7) + B{ Zwk'xlp Y1(Zw'\/kpz)}]eiwt Now, if we denote z at x = 0 by 20’ and z at x = L by ZL’ the second boundary condition yields Y1 (—2w—1<£_Q M) Z - 4o A B 2%,ng ( 8) J1 (__k. ) so that 20.)sz _ 1 20)sz o u—B Bum/p?0 {-JO( k )Y1(——k ) J1‘ k ) (4. 8) 2&N pz . 20W 2 z 0 not + YO( k ) Jl( k ) e With 20W p zL 20W p 20 20W p zL 20W p zoI D1:Y0( k )J1( k )-Jo( k )Yl( the first boundary condition yields . . . 1 .....tv ...: .. .m ( . .475... .. . ..v . nu , .. . 1 . 2 .7. 11.11 a . . .. in}... 33 ZwN/p z \ A = - y— Y (———3) D1 1 k and > (4. 9) B : ‘2 J (20%) p 20) D 1 k ) Eq. (4. 9) is the same as given by Datta (1956). Substituting these values of A and B in Eq. (4. 6), we obtain the value of u(x, t). The stress 0 (x, t) is given by _ _ QB _ 92 0'(x,t)—E(x)€— (E0+kx) 8x — kz Bz so that 0' (x, t) = -m/ pz’ { A31 (29—?— “kz) + BY Z—w—E— 1(2)} eiwt (4.10) 1 ( where A and B are given by Eq. (4. 9). Case 2 Let there be a source of disturbance producing a periodic stress Pelwt at x = L, and let the end x = 0 be free. The boundary conditions are (E + kx) 33] = Pelwt 0 8x x: L 8 0%] = 0 x=0 Again the second boundary condition yields Eq. (4. 8). Using the first boundary condition and writing Zw'x/ sz 20W pzo 20%}sz ZLN p 20 Dz Z Y1(“T‘)J1(_‘1’<—”" J1(_T<—_ Y1(‘—E—) 34 we get whence and Eq. (4.11) is the same as that given by Datta (1956). (4.11) Substituting these values of A and B in Eq. (4. 6) we obtain the value of u(x, t), and in Eq. (4.10) we obtain the value of 0’ (x, t). 4. 2(a. ii) Stress Wave Equation The differential equation is Eq. (2.13) when variations in p are neglected. 2 2 8 0' a 0' E(x) 2 : P 3x at Substituting the value of E(x) we get 2 2 (E0+kx) “Ba—:- 2 p g—Zg 8x at To obtain a periodic solution of Eq. (4. 12) let 0‘ (x, t) = 0' (x) ei‘”t (4.12) 35 The resulting ordinary differential equation for 0 (x) is d20' 2 (Eo+kx) --—Z +pw 0‘ :0 (4.13) dx We transform this equation to a Bessel's equation, Eq. (4.16) below, by the following substitutions. First let z = E + kx o to obtain 2 d20' 2 zk +pw0' = O (4.14) dz Then substitute 2 2 z : S k or s _ 263V pz k 4pw 2 d 0 1 d0 _ ds brings Eq. (4.15) into a convenient form of the Bessel's equation 2 + m|-' CLIQ. m (D 1 2 +(1-—-2-)e=0 (4.16) ds 3 whose general solution is 9 = A J1(s) + B Y1(s) s 0 that 3e 0 = S{A 11(5) + B Y1(s)} (4.17) Hence the complete solution of the wave equation Eq. (4. 12) is 0'(x,t) — EEL‘M {AJ1(ZL_L M) k (4.18) + B Y1 (Zw‘VkEZ)}elwt A and B are determined from the boundary conditions. Example Let there be a source of disturbance producing a periodic stress Pemic at x = L, and let the end x = 0 be free. The boundary conditions are _ iwt 0—]sz — P U]x=0 - 0 The secondary boundary condition yields Zw'x/pzo Y (——) l k A = B 2 , wquo so that 5(x t) : depz B {Y (EL “2) J (11:12) ’ k Zwsz; 1 k 1 k J1 ( k ) 2(1)sz . 2 «I _ Y (__2) J" (Jo—fl.” em’t (4.19) 1 k 1 k With 20.)sz 20)sz 2(0sz 200sz D=Y(———L-)J( O)-Y( °)J( L) 2 1 k l k 1 k 1 k 37 the first boundary condition yields 201sz D P - L 2 B R J (Zak/7723) l k whence ZwN/‘p-El' o k Y1(_ k "" 1 A = - P 21.1.)“ pZL D2 and 5 (4. 20) 2w'\/p z J (...—2.) B _ P k 1 k _ —_./' D 2(1) sz 2 / Substituting the values of A and B in Eq. (4. 18) we obtain the value of 0'(x, t). 4. 2(b). E(x) = E0 ekx 4. 2(b. 1). Displacement Wave Equation From Eq. (2. 9). neglecting variations in p , the differential equation is 2 u 09 8 Bu $115311“ 5;) = P 9’11 Substituting the value of E(x) we get 2 8 kx 8u _ 8 u E(Eoe 5; “P Ta (4'2“ 1'. Let -kx z = e Then Eq. (4.21) becomes 2 2 2 3 u _ 8 u EOkZ—T—p—Z (4.22.) 38 To obtain a periodic solution of Eq. (4. 22) let u(z,t) = u(z)eiwt Then Eq. (4. 22) yields an ordinary differential equation for u(z) 2 2 d u w z ——7 + _L_2 u = 0 (4.23) dz Eok To transform this into Bessel's equation, Eq. (4. 25) below, let SZkZE - —&”E Z — 01' S - -k-— E 4pm 0 §_2_L dz - k E z 0 L02 After division by 79—— , Eq. (4. 23) becomes k E o dzu 1 du EZ-S.d§+u:o (4.24) The change of dependent variable 11 = 30 then transforms Eq. (4. 24) into dZG 1 d0 1 gg+sfi+ujflezo (AM This Bessel's equation has the general solution 0 = AJ1(S) + B Yl(s) so that s {A J1(S) + B Y1(s)} (4. Z6) {3 ll 39 Hence the complete solution of Eq. (4. 21), as given by Sur (1961), is with z = e-kx u(x,t)= H/ O{AJl(—k—w (4%? O)+BY1( —k—“’ E:e)} t (4.27) A and B are determined from the boundary conditions. Case 1 Let there be a source of disturbance producing periodic displacement ert at x = L, and let the end x = 0 be free. The boundary conditions are _ iwt u]x=L — Ue Bu _ .23—x1x20— 0 Now, from Eq. (4. 27) Z Bu _ 20) z 21.) 2 21.) z iwt a; — ' T P—EO {AH—k LEO) + B YJ—k‘ LEO” 6 (4°28) Hence from the second boundary condition . 2w pZO Yo (‘1.— (if; ) A=-B (4.29) 20) pz J0(k O) 0 Writing (note 20:1) pz _ 20) j o _ _ Se- ? ? z]x=o-Z land Z]x:L ZL and (20) pzo 2w sz 2w pzo —1 (4.32) pz J (2w 0) k 0 k E0 B: 'P 2 D / 20) p 4 Eq. (4. 32) is the same as obtained by Sur (1961). Substituting the values of A and B in Eq. (4. 27) we obtain the value of u(x, t), and in Eq. (4. 3,1) we obtain the value of 0 (x, t). 4. 2(b. ii). Stress Wave Equation The differential equation, neglecting variation in p, is again from Eq. (2.13) 820 820' E(X) = P — 8x2 atz Substituting the value of E(x) we get E e —— =p9——- (4.33) To obtain a periodic solution of Eq. (4. 33) let 0' (x, t) = 0" (x) em)t Eq. (4. 33) yields an ordinary differential equation for 0' (x). .....s_.tfi:fl.....f...n.....s is... ,..Immense...”......z...££.._.i}..-......,3.i.t.1:...1...1_§3.sisal: . 4.543111: 4i 2 E 6 (LE 4,pr5 2 0 (4.34) o 2 dx This can also be brought into a standard form of Bessel's equation, Eq. (4. 36) below, by the following substitutions. First let -kx z = e to obtain dzcr do 8 Z Z 2 + d— + w 2 0" : O (4.35) dz 2 E k 0 Then substitute E s‘2 k'2 0 Z I 4 p w or S : Z_w P_Z d_s : Q _P_ k E ’ dz k E z o o (02 to transform Eq. (4.35), after dividing by —Z—P— , to k E o dzo 1 do 2+———+0'=O (4.36) 3 ds ds whose general solution is 0' : A Jo(s) + B Yo(s) (4. 37) Hence the complete solution of Eq. (4. 33) is o~ A = - B M 2w 0 Jo T E > 0 so that B iwt a (x, t) — Jog) { -J0(s) Yo(so)+ Yo(s) Joe-:3} e (4. 39) Again with ._ 2w IDZL 2w 920 2w pzL 2w F’Zo D4--Jo<‘k—/E)Yo(—k' EHYJ—k— _E ”on? E) O O O O the first boundary condition yields 2w pZo \ Yd? To ) A = - P D4 P (4. 40) 2w 20 510‘? E > B = + P D / 44 Using the values of A and B in Eq. (4. 38) we obtain the value of 0' (x, t). 4. 2(6). E(X) — ka( Ln) Here E0 is the elastic modulus at x = kL, n may assume any real value, and k is a constant whose value determines the total amount of variation of the modulus in a bar of given length. In the examples treated the bar will lie between x = kL and x = (k+1) L (see Lindholm, 1963). 4. 213. i). DiSplacement Wave Equation From Eq. (2. 9), neglecting variation in p, the differential equation is 2 8 Bu 5 u — (E(x) —) 2 p —— 8x 3x atZ Substituting the value of E(x) we obtain L{E (_X_)n 22} — 2.22 (4 41) 8x 0 kL 3x _ P atZ ° 2 n-l Z x 3 u n x Bu _ 8 u Eo(kL) ex 2 + Eo kL(kL) ax ‘ p at2 For a periodic solution let u(x, t) = u(x) eiwt to obtain an ordinary differential equation for u(x) 2 d u n x n-l du Z _ EO-k-( er1) —2 + E0 —kL(—kL) —dx+ pm u _ o (4.42) 45 For n )é 2 we transform this to a Bessel's equation, Eq. (4. 44), by the following subsitutions. The case n : Z is treated separately, beginning with Eq. (4. 52). First let z = x/kL to obtain 2 2 2 2 du+39‘—1-+4‘3—————ka Luzo (4.43) z dz E n dz 0 z (71 "1) Then let u = z 6 whence -1- 1- 213‘- — -n z(_z£)9 + z(—Z—n_) 51-6— dz dz and .l-n 2 -3-n -1-n (——) 2 d u : (l-nH-l-n) z(_—Z—)9+ (1_n) z(—-z—)_(_i‘§.+z Z d 9 2 Z 2 dz 2 dz dz l-n Substituting these in equation (4. 43) and dividing by z( 2 ) we obtain d29+1_d_9__{(1-n)2_l___1_9w2k2L2 )6—0 2 dz 2 Z n E — dz 2 z o Nowlet Z-n 2 __ y=(—Z_n)ka EL z( z) 0 4x_ _P_ -n/‘2 dz _ ka -E o 2 -Z—n d - :31 _P_ ( > 2 — 2 LokL E z 2 46 kaZLZ -n to obtain, after collecting terms and dividing by E9 Z . O 2 2 E é—%+_l_%g.+{1_.(l_'2n.)_ z(n'2)__ZTOZ._}G:O dy Y V w k L p which reduces to 2 2 e _ d__z+1_g_9_+{1_(12_2 —1_Z_}e:o (4.44) dY Y y m Y since yZ : 42 kaZLZ fig z(Z-n) (2.-n) 0 Eq. (4. 44) is a Bessel's equation of order p = ) %-:—:1i) with the general SOIUthl’l 9 : A .I y + B Y y whence l-n u = 2(T) { A Jp(y) + B Yp (y) (4.45) Hence the complete solution of the wave equation Eq. (4. 41), for all values of n except 2, as given by Lindholm (1963), is (L?) u(x, t) = z { A Jp(y) + B Yp(y)} eiwt (4. 46) where __}L Z ‘ kL l—n p = [—2 nl and ~< l NA I N 4|, 6 ... L“ AK 0 to) WA r34 ”I 47 A and B are determined from the boundary conditions. Case 1 Let there be a source of disturbance producing periodic displacement Uelwt at x = (k+l)L, and let the end x : kL be free. The boundary conditions are _ iwt ‘11 x:(k+1)L ' Ue 33.] z 0 8x x=kL Now, from Eq. (4.46) for n< 1 or n> 2 93— (122n)w—L{AJ ()+BY ()} 1“” (447) ax ‘ Z E 1 V 1 V e ' 0 -Zén -2-n andfor lin< Z 93- (1—22n)w -P—{AJ ()+BY ()} 1“” ax ‘ "z E0 1 V 1 Y e Z-n Z-n From the second boundary condition: for 1 <_ n < 2 Y21n(ykL) A = - B J - ( ) 1 VkL Z-n whfle for n< l and n> 2 (4.48) Y__1_ (ykL) A z _ B Z-n J (YkL) A Z-n 48 With D5 2 { Yp(y(k+l)L) JerL) ' Jp(y(k+l)L) Yr(ykL)} the other boundary condition yields l-n D UelLOt : (11:1) 2 B{ J 5(y ) } 816”: r kL so that l-n Y (y ) _ k r kL A ‘ mm) 2 { D - 5 and (4.49) l-n —— J (Y ) _ k 2 r kL B-Ufifi’ { D5 } where r : —1-— for 1 < n< 2 Z-n — rz—l— for n2 Z-n By substituting these values of A and B in Eq. (4.46), we can obtain the value of u(x, t). The stress 0‘ (x, t) is given by 0'(X,t) = E(x)€ = E (— Therefore from Eq. (4. 47) 0‘ (x, t) =7L(n-2) E0 zl/Z w / EP— { A Jr(y) + B Yr(Y)} eiwt (4. 50) o . .. .. .. a..-....w%_w§i.u§fi§3r¥f¢at£.cti¢w.411141.4u1351133mfiflua: 49 where r:-Zl— forlin<2 -n rz-l— for n2 Z-n A and B are given by Eq. (4.49). Case 2 Let there be a source of disturbance producing wt periodic stress Pe1 at x = (k+l)L, and let the end x : kL be free. The boundary conditions are E(x) g—u] = Peiwt X x:(k+1)L 2.11] : 0 8x x=kL Again, from the second boundary condition, with D6 2 Yr(y(k+l)L) JerL) ' Jr(y(k+1)L) Yr(YkL) we obtain 1 { r ykL) } ) wk] Eop D6 A = r(Z-n)U (k—1:—1)]7/2 (4. 51) J (Y ) k 1/2 1 r kL B = r(n-Z) U (—) ———___, {— where rz—l— for1<_n<2 N I lr—I :3 r=- forn2 50 Substituting these values of A and B into Eq. (4. 46) we obtain the value of u(x, t) of 0' (x, t). 4. Z(c.ii). Special Case for n = 2 For n : 2, the differential equation, Eq. Eo 28u kzlg (X 2 '5; For a periodic solution let u = u(x)e1wt. where Q. {3 with general 3 olution : X'Z/2{ A cos(gz) + B sin(gZ)} where N/c-124 g: and into Eq. (4. 50) we obtain the value (4. 41), is (4. 52) Then u(x) satisfies (4. 53) (4. 54) The complete solution of Eq. (4. 52), as given by Lindholm (1963), is (with z = loge x) u(x, t) = x-1/2{ A cos(gz) + B sin(gz)} e lwt (4.55) 51 A and B are determined from the boundary conditions. Let there be a source of disturbance producing Case 1 periodic displacement Uelwt at x = (k+l)L, and let the end x = kL be free. The boundary conditions are _ iwt u] x:(k+l)L — Ue 8 a—E] = 0 x=kL From Eq. (4. 55) 8 -3 2 . 5% :{x / g[-A31n(gz) + B cos(gz)] - 1? x-3/Z [A cos(gz) + B sin(gz)]}e1wt (4. 56) From the second boundary condition 1 . g cos(gzo) - -2- Sin(gzo) A = B (4. 57) gsin(gz )+ 1— cos(gz ) o 2 o where Z0 2 loge kL Using the other boundary condition, and writing 2L 2 loge(k+l)L . l YIO = g s1n(gzo)+ :- cos (gzo) _ . 1_ Y1L — g Sln (gzL) + 2 cos (gzL) 1 . YZO g cos (gzo) - E Sln (gzo) 52 Y2 = cos(z)-1—sin(z) L g g1. 2 g1. we obtain Y2 \ _ 1/2 0 A _ U{ (k+l)L} YZO cos(gzL) + Ylo sin(gZL) )(4. 58) Yl _ 1/2 0 B _ U{ (k+1)L} YZO cos(gzL) + YIO sin (gzL)J Substituting these values of A and B into the Eq. (4. 55) we obtain the value of u(x, t). The stress 0" (x, t) is given by 2 — _ X 2.2 0' (X, t) — E(x) E — EO (bl-(I) 8x or from Eq. (4. 56) E X1/2 1 (T (x, t) = O [-A { g sin(gz) + — cos(gz)} 2 2 2 k L + B { g cos(gz) -1§ sin(gz)}] elwt (4. 59) where A and B are given by Eq. (4.58). Case 2 Let there be a source of disturbance producing periodic stress Pelwt at x : (k+l)L, and let the end x = kL be free. The boundary conditions are E(x) g—u] = Pelwt X x=(k+l)L 21 5‘37] = o x=kL Again from the second boundary condition we obtain equation (4. 57). Using the other boundary condition, and writing 53 Y1O = g sin(gzo) +é- cos (gzo) YlL = g sin(gzL) + 13 cos (gzL) YZO = g cos(gzo) - 1Esin(gzo) YZL : g cos(gzL) - 1gsin(gzL) we obtain A : kZLZ { YZO } p EO{(k+l)L}1]‘Z YIOYZL-YZOYIL (4. 60) B : kZLZ { Y10 } p Eo{ (k+1)L} 1/2 YIOYZL - YZOYIL /‘ Substituting these-values of A and B in Eq. (4. 55) we obtain the value of u(x, t) and in Eq. (4. 59) we obtain the value of 0' (x, t). 4. 2(d). EFL) : E a Constant The differential equation is 2 2 E 8 u : 8 u (4. 61) 2 2 8x 8t icot 41(1),): u(x)?— For a periodic solution let u Knew/t so that u(x) satisfies (12 2 —‘—1-+P—ou=o (4.62) 2 E dx with general 5 olution u=Asin(wa-]%)+Bcos(wx ’EE.) The complete solution of Eq. (4. 61) is u(x,t):{Asin(wx«/g‘)+Bcos(1.1x/-}%)}emt (4.63) 54 A and B are determined from the boundary conditions. Case 1 Let there be a source of disturbance producing iwt periodic displacement Ue at x : L, and let the end x = 0 be free. The boundary conditions are u]X:L _ Ueiwt 22;] — 0 8X x20 — From the second boundary condition A = O, and u(x, t) = B cos (LOX /%) eiwt From the other boundary condition The complete solution is cos (to XJE) cos (wL J? iwt u(x, t) = U The stress a (x, t) is given by 0'(X,t) : E— whence sin (wx /-E) . x, t) = - UwVpE ' E €31“)t U( cos (toL 1%) (4.64) (4.65) (4. 66) 55 Case 2 Let there be a source of disturbance producing periodic stress Pelwt at x : L, and let the end x = 0 be free. The boundary c onditions are From the second boundary condition we get A = 0, while the first boundary condition yields B : — P 1 uNpE sin (wLE) cos (wa-EI) . u(x, t) = - P l E emt (4. 67) “PE sin(wL LIE) and sin (to x £) 0' (X, t) — P E sin (1.1L / .123 In both the cases of constant E (4. 68) 1. u(x, t) = 0 when cos (1.) x fig ) = O, or displacement nodes are JEP n=0.1.2,3,... (4.69) when these positions fall inside the bar. located at Elli- x=(n+lz) Tr 56 Z. 0' (x, t) = 0 when sin (to x (fi) = 0, or stress nodes are located at x=n1T,1-(EF n=0,1,2,3,... (4-70) when these positions fall inside the bar. 3. For the displacement boundary condition, Case 1, resonance occurs when cox (“LN/"1333) = O; and for the stress boundary condition, Case 2, resonance occurs when sin (01L (3%) = O; or when Casel 1 W E w 1‘ (n+3) "I: (F n=0,l,2,... (4.71) CaseZ _ n_w E _ L0 — L p n—1,Z,3,... For a 20-inch bar of Type 303 stainless steel some of the resonance frequencies are: For Case 1. (Displacement B. C.) i. Bar at room temperature (75oF) 2,470; 7,410; 12,350; 17, 290; 22, 230; cps. ii. Bar at 1200°F 2, 138; 6,416; 10, 690; 14, 966; 19, 242; cps. For Case 2. (Stress B.C.) 1. Bar at room temperature (750F) 4, 940; 9, 880; 14, 820; 19, 760; 24, 700;. . . . cps. ii. Bar at lZOOOF 4, 275; 8, 551; 12, 826; 17,101; 21,376; cps. ..i. . . J1... ajj__§¥§w¥.§mhir . wififliionigfilflfiiflnmih. - 1.41: . 1.1.. .91.]... .... ...... 57 4. 3. Numerical Solution (Case of General Variation of Temperature) 4. 3E). Formulation of Algebraic Equations From Eq. (2. 9), neglecting variation in p , the displacement wave equation is Z a au _ 8 u a (E(X) 5;) - P _2 (4.72) at For a periodic solution let u(x, t) = u(x) elwt so that u(x) satisfies the ordinary differential equation d du 2 _ a; (E(x) a?) + p (.0 u — O dzu 1 dE(x) du 902 (J + E(x) dx 3; + E(x) u :0 (4°73) Divide the length of the bar into a mesh of n equal intervals of length h : L/n. The central-difference quotients approximating the derivatives are d3 z 1+1 1-1 dx 2h dZu z u1+1 + u1-1 - 2 ul de h2 —- 1 < j < n+1 n _ _ Also x}. = (j-1)h = (j-l) Eq. (4. 73) is thus approximated at each mesh point by ui+1 ’ Z U‘i + ui.1 + 1 E(XiH) “E(Xi—1)} (ui+1 ”Ii—1‘ I hZ E(Xi) 2h 2h 2 L : + E(xi) ui O >‘In order to agree with the computer program language the numbering of the points varies from 1 to n+1 instead of from O to n. 58 or - 4 E(xi) ui-l . 2 2 p w h + (-————-—E(Xi) — 2) ui {E( . ) - 1*3( - )} +[1+ x1“ x1'1 ] u. = 0 (4.74) 4 E(xi) We write Eq. (4. 74) at each interior mesh point and also at an end point where a stress boundary condition is imposed. The stress boundary furnishes a condition on 3% , which permits expressing the value of u at one fictitious mesh point outside the bar in terms of the value at an adjacent interior mesh point. Solution of the simultaneous system of linear equations then gives the value of u at each mesh point. We consider two different cases of boundary conditions. Case 1 Let there be a source of disturbance producing a Lot periodic displacement Ue1 at the end x = L, and let the end x = 0 be free. The boundary conditions are _ iwt u] x=L — Ue u - u 32] : 0 = Z O x x=0 2h Now, writing equations for points 1 to n in Eq. (4.74) and using the 59 boundary conditions, we obtain for i=1, wzh2 (BE—(gr-Z)ul+2u2=0 i=2, {E(x3) — E(xl)} w2h2 [1 ' 4 E(xz) ] u1 H E(xz) ’ 2’ u2 +[l+ 4E(x2) u3 --0 1:3, { E(x4) — E(x2)} 02112 [1 ’ 4 E(x3) u2 + ( E(x3) ' Z) “3 ”1+ 4E(x3) “4:0 n’ { } E(x +1) - E(x 1) thZ [1‘ 4E(xn) n Jun-1+(E(x ) -Z)un {E(x ) - E(x )} ”1+ I4+iz(xn) M ] U=0 We have n equations and n unknowns. These equations can be written in the matrix form AU'rK where A is the tri-diagonal matrix with ai,i+1 The rest of the f. ._ 1r : 6O 12 O ’ ’ ’ ’ " " 0 a 22 a23 ” ” ' ‘ ” 0 l 32 a33 \ \ \ a43 \ \ \ \ \ (4.75) \ \ \ I \ \ \ \ ‘ \ ‘ an-1,n.2 an.1,n-1 n-l,n _ — - O 3'n n-l ann _) thZ = 1%}? for i=1, 2, ., n i :[1 {E(x-+1) - min” 4 E(xi) ] for i=2, 3, , n — 2 4 E(Xi) for i=2, 3, , n-l aij are zero. To 0 where 0 {E(x ) - E(x )} ° n _ 4 E(x ) n k n.) 61 Now to obtain the value of stress from these values of displacement we use the relation 0'(x) = E(x) %% so that u. -u. ) 0- (xi) = E(xi> ( 1“ 2h 1'1 ) (4.76) Equation (4. 76) gives us the value of stress at any interior mesh point. For the stress at the end x : L, we can use a backward differ ence formula 1 un+1 — un a (an) = §{E(xn) + E(XDHH <—h——> (4.77) Case 2 Let there be a source of disturbance producing a periodic stress Pewt at the end x = L, and let the end x : 0 be free. The boundary conditions are 8 un+2 ' un iwt itot E(x) 5%] = E(xn+l) (T) e = P8 x=L 23] — uz - uo — 0 8x x=0 _ 2h _ Again writing the equations for points 1 to n+1 and using the boundary conditions we obtain the same equations as in the first case for 121 to n-1, and in addition we get for _ : “3......’inimgfiiiTa;...quaiafiwfiflahfifihfimigu1¥nH.115 .\ 62 1:n { E(x ) - E(x )} Z Z [1 ' nllmx ) l un-l (%'h_) ' 2) un n Xn + 1+{12(an) — E(xn_1)} u _ 0 4E(x ) n+1 — n i=n+l pwzh2 Zun+(E(X ) -2) un+1 n+1 2h P { E(an) ' E(1%)} _ +E(x )[1+ 2E(x ) *0 n+1 n+1 This again can be written in the matrix form A U : K where the tri-diagonal n+1 by n+1 matrix A ('— - — — — — _ — — — all a12 0 O O a‘21“‘“22 a23"-—“—‘-0 \ I 0 6132 6‘33 \ \ l O O a \ \ \ l _ 43 \ \ A — . . \ . . I (4.78) l | \ \ \ I I ) \ \ \ l ' \ \ \ ' l ' an,n-l an’n an, n+1 L0 O h _ _ _ _ _O an+1,n an+l,n+l has elements _ BALL _ aii — (E(Xi) Z) 1—1, 2, ., n+1 8L1,1.1 2“ ‘ 4E(x.) ] 1:2; 3’ ’ n 1 63 an+l,n — 2 a12 ’ 2 E . -E . a : [1+{ (Xlfl) (xl'1)}] i=2 3 n i,i+1 4E(Xi) ’ The restofthe aij are zero. p—ul - ‘0 '-1 u 0 2 'where u3 O ZhP { E(xn+1) ' E(Xn)} U= K= k+1=-B(——) 1+ 21:6 > n xn+1 n+1 u 0 n un+1 kn+1 Again stress at any point of the mesh can be obtained by using Eq. (4. 76) or Eq. (4. 77), depending upon the point. 4. 3(b). Solution of the Algebraic Equations We have to solve the equation AU=K where 63 an+l,n _ 2 a12 * 2 ai,i+1 : [1+ 4E(xi) 122’ 3' ‘ ’ n The rest of the aij are zero. u 0 2 where u3 O _ ZhP { E(Xn+l) ' E(Xn)} U- K: kn+1 E(x )[1+ 2E(X ) n+1 n+1 u 0 n un+1 kn+1 Again stress at any point of the mesh can be obtained by using Eq. (4. 76) or Eq. (4. 77), depending upon the point. 4. 3(b). Solution of the Algebraic Equations We have to solve the equation A U = K where all 5‘12 0 0 ’ ' ’’’’ ” ’ 2121.422 a23-—-—--*— 0 3L32 a33 \ \ \\ \ A :: I \ \ \ I \ \ ' I \ \ \ I \\ \ \ I \ \ 1 I ‘ ‘ | . m-l, m-Z m-l, m-l O 0 - - - - O a m, m-l ul O 112 O u 0 u: 3 K: tum tkm 4. 3(b. i). Direct Solution The matrix equation can conveniently be solved by the direct method (see Hammings, 1962). We can write the tridiagonal matrix A as the product of a lower triangular matrix L and an upper triangular matrix Z. the original matrix equation LZUzk, Then 65 by use of ZU H < becomes LV K The resolution into two triangular matrix factors is made as follows . For the tridiagonal matrix A it is possible to choose the triangular matrices so that all elements are zero except those on the diagonal and one adjacent parallel row. Fl “ 11 1 21 22 <:> £32 £33 \ \ \ \ fl ‘ ‘ \ \ \ \j m, rn-l m, m ha. 11 a21 1- By matrix multiplication 111 Z11 2 a11 [11 Z12 2 a12 F Z 11 z12 Z22 z23 z33 23 33 a34 \\\\ \ \ and the rest of zij 66 :Oforj>2. If we arbitrarily choose all diagonal elements of L to be unity then Now 01‘ Similarly yielding 01‘ N ll (1: {suff’hikf undone-1f... ..1... .. II , .I...» t . .\ . M _ s . . , . .. . . _ s5” .1. ( ..- .- ,. . (hung—«4.3, .. .... .. .... ,. . .. , .. 67 01‘ Z33 2 a33 ' 132 Z23 In general I - : rr 1 Zr, r+1 ar, r+1 Z11 a11 _ ar+l, r z _ a z r+1,r_ z rr rr r,r-1 r-1,r rr a'r, r-l a - — a : rr Zr-I, r—l All other elements of the two triangular matrices are zero. Now to find V components of the equation L‘V == K 11 /, v1 k1 21 £22 I \) V2 k2 1 ‘\’/ k 32 33 V3 3 \ Z \ \ \ \ A \ \ o J \ \ v k m m m, m-l mm we proceed as follows: From Eq. (4. 80) we obtain (4.79) r—1,r (4. 80) so that v1 2 k1 : .. 2 V2 k2 21 V1 V3 2 k3 ‘ £32 V2 or : - I Vr kr r, r-l Vr-—1 Further we have Z U = V " '1 (- — (- n Z11 Z12 ' _ u1 V1 Z22 Z23 O uz V2 . \ \ \ \ \ \ u3 v3 \ \ . z \ . \ zm-l,m . z u v mm m m L. ..J .... _. .... .- from which 2 ‘ = v mm m m zm-l, m-l um-l + zm-l,m um — Vm-l So in general — < < .. er ur zr,r+l ur+1 v for l_r_ m1 z u = v for r =m rr r r 68 (4. 81) (4. 82) ) (4.83) IIIJI‘ in , .- . 1. Wyn .. , 44.11113: .0. ...» u 1 _. Ir ., . ,_ -.....l..n|.. . . . . I 4. n. p , . . , I . . I .L-11k _ n - ...) .1 ... ... . ....a .J. .... inf: . u . . Inigo. Jim-.1 -. 1 K ,,._ .... . 69 From Eq. (4. 83) m u : m 2 mm (4.84) and u = (v - z u ) -1-— r r r,r+l r+1 z rr The above is the direct method of solution as given by Hammings (1962). For application of this method, we had for Eq. (4. 81) kr = O for r )4 m and km )é 0 Therefore vr = O for r [g m V = k for r = m r r Hence Eq. (4. 83) and Eq. (4. 84) reduce to =0 forl:r vomnm mbus o.m .Mfim “monofldv m¢m mmB mo 92m Bow Mme 20mm m02w ooo.oww ooo.oon ooo.oo~ (was / WI) (133.18 mm 112 0.86“: on» an condom Adm vacate on» you oapmmnopomhmg pdofim 283 mnflvwoq fin .wflm amonocuv m¢m Ema mo 92m Mme 20mm moz Hwhoflow ¢ N .m .me 121 anosdhonkm Uncoom 0:» non 000m on» Mo mafiduom Anvm.m .wgm hum. noawoonm .n 0 0m0¢ ham xuhum 0m0¢ honufifimawfi I. 09500.,000809uamum 096 nouwuhe u 41 , // :mmeoe am. 0 \ , T... / :0 3 Tai 0Hmm¢ mmabm 20mm Dmmbm .ufim A0235; moz.9mum hocosvobfl $6 . omflm {Samoa n Tmoham .o .m. 444 (ax) ssaums 114-5 .000 ooomme you m0sufiags< pacemowammfin 0m.~ .wum A monondv WoZdBmHQ 0m , 0P ON“..- .ml prnofiahcnxm _ 0mg! 0.mn wxwp.emnm .m .m m.mno ooo.me honosvohm _ ‘0. It.“.ll.‘.hll‘ lllllllllllll .4 l o.“ 0.? . (a Lgtx) muawaovudsxa 114.6 ON .mao 000.0N you mcsuaage< unoppm mm.~ .000 Amonondv moz Y(M)=Z(M+L-l) DO 905 K=1.26 DO 900 11:6“! 1:205 00 900 J=1oll J1=J+I~1 P(IQJ)=((P(I~191)*(Y(J1)-X(K)))-(P(I-19J+l)*(Y(l-1)-X(K))))/(YIJ1) 16'? DO 804 K1=269430.10 K2=K1+10 DO 802 M=lo5 P(IoM)=C(M+L-l) 802 Y(M)=Z(M+L-1) DO 805 K=K1oK2 DO 800 1:2o5 1126-! DO 800 J=10Il JI=J+I-1 800 P(IoJ)=(-G(K))))/(Y(J1) C-Y(I-l)) 605 X(K):P(5¢l) L=2 ‘ DO 304 K1=29o430o10 K2=K1+9 00 302 Malos P(1.M)=Z(M+L-1) 302 Y(M):T(M+L—1) DO 305 K=K1.K2 DO 300 I=2o5 1126-! DO 300 J=loII J1=J+I‘I 300 P(IoJ)=((P(I-1.1)*(Y(J1)-G(K))1—(P(I-I.J+1)*(Y(I-1)-G(K))))/(Y(J1) C-Y(I-1)) 305 X(K)=P(5o1) 304 L=L+l L=43 DO 733 M=1.5 P(IoM)=Z(M+L-I) 733 Y