THE UNIVERSITY OF MICHIGAN INDUSTRY PROGRAM OF THE COLLEGE OF ENGINEERING THE AXISYMMETRIC RESPONSE OF A CLOSED SPHERICAL SHELL TO A NEARLY UNIFORM. RADIAL IMPULSE David A. Sonstegard A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy in the University of Michigan Department of Engineering Mechanics 1965 September, 1965 IP-716

Doctoral Committee: Associate Professor Ivor K. McIvor) Chairman Professor Emeritus Ruel V. Churchill Professor Samuel K. Clark Professor Donald T. Greenwood Assistant Professor Roger D. Low

ACKNOWLEDGEMENTS This area of research was suggested by Professor Ivor Ko McIvor, my committee chairman, I gratefully acknowledge not only his invaluable instruction and guidance but also his personal enthusiasm and stimulus. I express my sincere appreciation to all committee members for their interest, suggestions, and efforts, I am indebted to the Alfred Po Sloan Foundation, the Consumers Power Company, the National Science Foundation, and the Ford Foundation for financial aid given to me during my graduate studies, David Ansel Sonstegard ii

TABLE OF CONTENTS Page ACKNOWLEDGEMENTS.. o o o o. a. ~ o o o o o o ao o o o ii TABLE OF CONTENTS 0 0 0 0 0 0. 0 o o o iii LIST OF TABLES o ~ o o a o ~ o. o o o o o o o. o o o o iv LIST OF FIGURES o o o o o o o o o o o o o o o...... v NOMENCLATURE o o o o o o o o o o o o a o o o o o a o o o vi CHAPTER lo INTRODUCTION o o o o... o 1 CHAPTER 2, A NONLINEAR FORMULATION. o o o o.. o 4 2.1, Representation of Deformation.. 0..... 4 2o2o Midsurface Quantities. o. o. o o o o o 5 2,3. Strain-Midsurface Quantities-Stress Relations... 8 2o4, Equations of Motion for Small Displacements....... 10 CHAPTER 3. A LINEAR STUDY o o o o o o o. o..... 16 3.1, The Linear Solution. o o o o o o o o. a o o o o.. 16 352. Mode Classification o.......o. o 20 5o3o Response to an Arbitrary, Symmetric Radial Velocity ~ 27 CHAPTER 4. RESPONSE TO A NEARLY UNIFORM RADIAL IMPULSE o. 34 4,01 Formulation.... o. o o o o..... 34 4, 2 The Approximation o o..... o o... o. 36 4,o3 Stability of the Breathing Mode.o.. o... 38 4, 4 Long-term Behavior............. 47 4o,5 Conclusions o.. o. o o o o........ 53 APPENDIX...... 0 0 0 0 60 BIBLIOGRAPHY o, o o o o o o.... o.. 62 iii

LIST OF TABLES Table Page 1o Angular Frequencies o o o O o. o o o o o o o 0... 21 2o Mode Indicators o o o o o o o o o o 0 0 o 24 iv

LIST OF FIGURES Figure Page 1. Undeformed and Deformed Configurations........... 4 2. View on a Meridional Cut.............. 8 3. View on a Meridional Cut................. 9 4. Plot of Amplitude Ratios..... o...... 31 5. Mathieu Stability Diagram................. 55 60 Stability Diagram for Vmax............ 56 7. Stability Diagram for.6 Vmax.......... 57 8. Long-term Behavior for Vmax.............. 58 9. Long-term Behavior for.6 Vmax........... 59 v

NOMENCLATURE Ao, Anmn Anc Constants of integration An Amplitude of critical deviation after complete energy transfer E Young s Modulus Eb Sum of maximum strains minus sum of midsurface strains Em Sum of midsurface strains L Lagrangian Mm Momentum along polar axis T Kinetic energy Ti Energy imparted by impulsive pressure. V Potential energy Vnax Maximum initial dimensionless velocity for ~y/E = 5(10)-3 VN Potential energy of Nth deviation after complete energy transfer a Radius of spherical shell aoan Coefficients of Legendre polynomial expansion of., bn Coefficients of associated Legendre polynomial expansion of 4. c Wave velocity, [E/p ( —v2)]2 dso0, dson Midsurface elemental lengths before deformation ds, dsn Midsurface elemental lengths after deformation h Shell thickness i, j, k Unit vectors along x, y, X vi

ry G 0 Spherical coordinates t Time t9 tng n Unit midsurface vectors u Meridional displacement u Meridional displacement with respect to moving reference Vo Uniform radial velocity w Radial displacement, a-r o~ w Radial displacement with respect to moving reference x, y, e Cartesian coordinates z Distance from midsurface Lo, ALnc 9 Lanm Lnc, Anm Amplitudes of motion due to arbitrary radial velocity distribution Tip Period of fundamental mode Qns A.n Parameters in Mathieu equation a2 TThickness parameter, h2/12a2 09o anc, anm Constants of integration PffiB~ ~Angle between midsurface normal and radial ray YO Midsurface shear strain 5ncy 5nm Amplitude ratios EAdd E0O Midsurface normal strains en Dimensionless perturbation quantity ce En Normal strains JtO a2 TrNon-dimensional displacements vii

%an I 14Midsurface curvatures v Poisson's ratio,.3 Polar angles P Mass density ao Stress due to breathing mode Oy Yield stress ax, o Normal stresses T Dimensionless time T Dimensionless time, UQT o,} rnci9 Mnm Angular frequencies of breathing, composite, and membrane modes respectively t2c Angular frequency of fundamental mode viii

CHAPTER 1 INTRODUCTION Early studies of the dynamics of spherical shells were attempts to construct a theory of bells. Lord Rayleigh(l)considered an inextensional representation of the motion, ie. he assumed that no stretching of the midsurface occurs during deformation. Love(2) took exception with the inextensional representation for two reasons: first, a closed shell can not undergo deformation without midsurface stretching; second, the restriction of inextensionality does not allow satisfaction of boundary conditions for a freely vibrating open shell, namely a bell. Love's formulation of the problem(3) constitutes the classical bending theory of shells now known as Love's first approximation. This formulation includes both flexural and extensional effects. In actually solving the problem, however, Love assumed the vibrating shell to have negligible resistance to bending. He considered only extensional axisymmetric motion. Lamb(4) also used an extensional formulation in his study of the radial motion of closed spherical shells, By the extensional approach Love found two sets of modes of vibration. Both sets contain an infinite number of modes. For the first set, called the upper set, the frequency spectrum is unbounded. All upper mode frequencies are higher than those of the second set, called the lower set, which has a bounded frequency spectrum. The boundedness of the lower set spectrum means that intervals between natural frequencies must become arbitrarily small, a physically untenable conclusion. -1

-2An analytical study by Silbiger(5) and experimental work by Baker(6) were based on the extensional analysis of Love, Naghdi and Kalnin:s(7) applied classical bending theory to obtain a solution for torsionless axisymmetric motion. Also included in this paper is a study of asymmetric motion based on extensional theoryo Kalnims(8), using classical bending theory, investigated the flexural vs membrane strain energy content of axisymmetric modes of vibration. Two mode sets were foundo For one set the deformations were predominately extensional, For the other set, however, deformations were neither predominately extensional nor predominately inextensional. In the limiting case of infinitely large radius-to-thickness ratio this latter mode set was shown to degenerate to the lower set found by an extensional analysis. The paradoxical boundedness of the lower set frequency spectrum as predicted by an extensional analysis was thus shown to be caused by the effects of bendingo In 1883 Lord Rayleigh(9) reported an experiment in which a stretched string was attached to one prong of a tuning fork. The fork vibrated in the direction of the string. Lateral vibrations of the string of frequency f were produced by fork vibrations of frequency 2fo This phenomenon is termed heteroparametric excitation (the prefix'hetero" being commonly omitted), Parametric effects occur frequently in mechanics. Numerous examples are presented by Minorsky(1O). Bolotin's book on dynamic stability(ll)is devoted entirely to heteroparametric excitation.

-3Autoparametric excitation is a similar phenomenon. An example of this behavior is a mass suspended on a spring as reported by Gorelik and Witt (see(l0) po 506)o The system is allowed two degrees of freedom, spring elongation and plane pendulum motiono For certain elongation modependulum mode frequency relations nonlinear coupling between the modes alters the initially excited motion markedlyo Energy originally imparted to elongation motion is transferred to pendulum motion. Both hetero-and auto-parametric excitation are generally governed by differential equations with periodic coefficients. Two distinctly different physical phenomena occur however. For heteroparametric excitation an external energy source, such as Lord Rayleigh's tuning fork, produces the parametric effecto For autoparametric behavior no outside energy is supplied to the initially excited system, only an energy redistribution occurso The energy exchange is due to nonlinear coupling between distinct modes of motiono Autoparametric behavior has been reported for the cylindrical shell by Goodier and McIvor(l2). Here the energy exchange is effected by nonlinear coupling between extensional and inextensional mode sets. Analogously, two mode sets exist for the closed spherical shell. However, one mode set is not of a distinct extensional or inextensional charactero The possibility of autoparametric excitation is, therefore, of particular interest and gave impetus to this investigation.

CHAPTER 2 A NONLINEAR FORMULATION In this chapter a Lagrangian representation of spherical shell deformation is presented. The formulation is restricted to torsionless axisymmetric motion, Nonlinear equations of motion are derived for a second order theory which is restricted to small deformation of thin shells. 2.1 Representation of Deformation The undeformed configuration is the closed spherical shell of radius a. A point P on the undeformed midsurface is located by angles C and i, Figure la. During deformation P goes to P* on the deformed midsurface and is located by spherical coordinates r,, and 0, Figure lb. s. I y y Xu 1 U r YC ri Firx e la x. lbFigure 1, Undeformed and Deformed Configurations -4

-5Motion of point P is described in the Lagrangian sense. The deformed state coordinates are taken as functions of initial position and time, ioe. r = r ( t ) t) (21o.la) a = (, r t) (2.1.1b) 0 0 ( t3 t) (2.1.1c) 202 Midsurface Quantities In the Lagrangian representation of Section 2.o unit vectors on the midsurface are expressed as t.. _ 1 [rSinoCose - reSin0Sine [r2 + r2&2Sin20 + r222] + r Cos0Cos~] i + [rSin0Sin~ + r@Sin0Cose + r0Cos0Sin~] J + [rCosO - r OSin0] ) (2,2.la) t = 1. r'Sin~CosG - rE'Sin~SinG [r'2 + r28'2Sin2~+r2t2] 2 + ro'CosoCose] i + [r'SinjSine + rG'SinOCose + ro'CosoSinE] j + [r'Coso - re'Sine] k' (2.2.lb) where i, j, k are unit vectors along x, y, a respectively, the dot denotes a/S and the prime denotes a/ao. From differential geometry the normal strain ec and the curvature KoS of the midsurface along t are Eo2 = [r2in2 + rSin2 ]2 1 (2.2.2a) a

-6Ko~ =.... 1 Ir42Sin2o + 4r4p2 [2+ r2&2Sin2~ + r22]3/2 L + r [- 4r2r' (2Sin2O + (2) + 4r3 ( 8 GSin2O + e2~ SinoCos0 +'] 222 2 42 4k 2 422 + r2[ (8 Sin2 + 4Sin2 + 3 4Sin Z + -2+ 4e2~2 + 4e2~2Sin20 + 4 e ieSinoCoso - 2&2~Sin0Coso + 4i4 ) - 2 r r(6 eSin2 in2Co + 2 SnCos +') + 2(e2Sin2 + ~2)] + r3[2r (&35Sin40 + e4Sin35Coso + Q20o Sin2o + e 2Sin26+~2 3Sin~Coso + i3 ) - 2 r(4Sin41 + 2e2~2Sin20 +' 4)] + r4[66Sin40 + 4(3 2Sin20 + 2Sin20Cos20 - 2(Sin35Cos0) + 2 a sini3Cos + 2(Si2 2 + 3 4+ 4Cos2 - 4 02,Sin0Cos0) + ~^2s in2o + 2 e e ~(2 2Sin0Cos0 - 0Sin20) + *06] 2 (2.2.2b) Similarly, along t~ the normal strain and curvature are con 1 [r'2 + r202Sin20 + r22'2] - 1 (2o2o2c) a Sin t K011 Kg with dots replaced by primes (2o022d)

-7The midsurface shear strain yo is 7o =..r r' + r2e e'Sin20 + r2".' (22.o2e) [r2 + r2~2Sin2o + r202]2 [r'2 + r2e'2Sin2o + r2 f,2]2 One of the difficulties encountered in even a small displacement approximation of (2.2o2) is the non-orthogonality of the ~-r curvilinear coordinate system. Orthogonality results from the introduction of axisymmetry, whereby r' =:' = O (2.2.3a) 8 - =0 (2.2o3b) The latter requirement precludes circumferential or torsional displacements. Imposing conditions (2.2.3) on (2.2.2) yields the following midsurface strain and curvature expressions: 1 e = 1 [r2 + r2~2] 2 -1 (2.2.4a) a - r Sin 0 1 (2.2.4b) a Sin > K - 1 4 r 4 + r [45r 4 - 4r2r }2] [r2 + r2"2] 5/2 + r2[r2 (2 + 4~4) -2 r'r 0' + r22] + r3N[2r 3 - 2' 4] + r46} 2 (2.2.4c) KOT= 1 (2.2. 4d) a Sin 0 7o = 0 (2.2.4e)

-82.3 Strain-Midsurface Quantities-Stress Relations An element as viewed on a meridional cut before and after deformation is shown in Figure 2. I/ z / / // 0' Figure 2. View on a Meridional Cut The angle between the surface normal vector i? and the radial ray may be expressed as p = tan- r (2.3.1) r Normals to the undeformed midsurface are assumed to be unstretched and to remain normal to the midsurface after deformation. Midsurface elemental lengths along te before and after deformation are respectively dsog = adS (2.53.2a) ds5 = (1 + ~oe) add (2.3.2b) The normal strain along tS a distance z from the midsurface follows as

-9E = (1 + o) (1 + z Ko) - (1 + Z/a) (2.53) 1 + /a The right side of equation (2.353) is expanded in powers of z/a. The shell is assumed sufficiently thin to allow neglect of z2/a2 with respect to unity. The expansion then yields e = Co + z(0Ko - 1) (1 + eoe) (1 - z) (2.5.4) a a A deformed element as again viewed on a meridional cut is shown B A/ / in Figure 3. / KSin(0-P) CO Cm / i a/ Figure 3. View on a Meridional Cut/ Figure 3. View on a Meridional Cut

-10 - Midsurface element lengths along t before and after deformation are respectively dso = a Sin d dr (2.3.5a) ds~ = (1 - E~o) a Sin g dr (2.3.5b) The midsurface strain along t a distance z from the midsurface is En = [1 + z Ko Sin (0 - )] [1 + Co] - (1 +Z/a) (2.5.6) 1 + Z/a Again neglecting z2/a2 with respect to unity, expansion of the right side of (203,6) in powers of z/a gives E = Eoq + z [so0 Sin (0 - ) - 1] (1 + e0o) (1 - z/a) (253~7) a The continuum is assumed to be isotropic and homogeneous. Neglecting stress normal to the midsurface, Hooke's Law yields biaxial stress-strain relationso They are =_ E (C + v E~) (2.3.8a) 1 - V2 a = (E E + v ES) (235.8b) 1 - v2 where E is Young's Modulus and v is Poisson's Ratio, 2~4 Equations of Motion for Small Displacements Equations of motion may be derived by an energy formulation or by consideration of the equilibrium of a shell element. Since the energy formulation is convenient in the subsequent analysis it is to be used here,

-11Assuming negligible rotary inertia the kinetic energy T is 2 2 T = [ ( r I + 2 f ]dV (2.4.1) volume where p is the mass density0 A non-dimensional time TA a non-dimensional radial displacement, and an angular meridional displacement 4 are defined as T = E 2 t (2.4.2a) a P (l - v2) = 1 (a - r) (2.4,2b) a' = 0- ~ (2.4.2c) Substitution of T,,, and / into (2.4.1) followed by an integration over the circumference and thickness yields T = Ha2Eh F ~_2 + (1-20) I 2 Sin d 5 (2.4.3) 1-v2! -r 0 L() Displacements and their derivatives are assumed small;therefore, termsof order greater than three are neglected in (20 453) For free vibrations of the assumed conservative system the potential energy V is the strain energy stored during deformation, namely v = 1 rio + el dV (2.4.4) J volume From substitution of stress-strain relations (2.358) there results V = E A [c + 2 v + ] dV (2.4o5) 2(1 - v2) volume

-12The midsurface quantities and the angle P may be expressed in terms of f and ~ by substitution of (2o4o2) into (20o24) and (2.3o1). The results are o1 = -; - A 1+ 3 (2.4.6a) 2 Eol =- + Cot - - C or Cot i (2a4.6b), -=1 (1 +; + + i2 + 1 I2 + 2 - _ - 2 ) (2.4.6c) a 2 Ko= 1 (1 + C - V Cot i + +2 + *2 Csc2 C - 4 r Cot) (2.4.6d) a Sin t P3 = ~- +: - h (2.4.6e) Consistent with the retention of third order terms in (24o.3), terms of order greater than two are neglected in (2.406). The strain energy (2.4o5) may now be expressed in terms of the midsurface quantities by substitution of (2.036) and (2o357), and subsequently in terms of g and 4 by substitution of (2.4.6)o Terms up to the third order are again retained after integration over the circumference and thicknesso These terms occur in two forms, the first with coefficient h and the second with coefficient h3/12a2. To be consistent with the thin shell assumption of Section 253, third order terms of the latter form are neglected. The resulting potential energy expression is

-135Ir V = n a2E h, 2 + I2Cot25 + 2(1 + v) ( 2 1 - v2 J 0 - 5 i Cot - t Jf + a2i Cot 5 + 22,) - 2 f (2 + 2~ _ i2+ f * 2 - l3Cot 5 - 2 ~ f2Cot2S + 2v(i f Cot 1 - ~ 2i -_ 2 ~ if f Cot 5 + ~ i f2 ~ 5 I2 + 1'2i Cot ) + a2[82+ i2Cot2t + 2 +2 Cot + 2 o + + 2 Cot2 + 2v(i i Cot | + [ i Cot C +' i Cot + C o Cot ) Sin ~ d ~ (2.4.7) where a2 9 a thickness parameter, is defined as a h2 12a2 The Lagrangian L is defined as the difference between the kinetic and potential energies, io.e L = T - V (2.4o8) Hamilton's Principle states that the motion realized in nature is that particular one for which the time integral of L assumes a stationary value. The principle may be expressed r t2 tl(2.4.9) J1

-14where ti and t2 are arbitrary but fixed times, and. denotes the customary variational operation. The Lagrangian is expressed as (2.4.3) minus (2.4.7). In order to satisfy (20 49), t and. must satisfy the two nonlinear partial differential equations (see, for example, Hildebrand(13) p. 137)' + 4 Cot! - K(v + Cot2!) - (1 + v) - 2 2 - 2 G t +:: - 2 t: Cot! + 2 52Co-t — (1- v) + ~ 2Cot 5(3 - v) -; (1 - v - 2 Cot2) - 2 v V Cot t + 2 5 i (1 + v) + a2[" + + Cot 5 - t(v + Cot2S)j +' + Cot 5 - ~(v + Cot2))] = (1 - 2 ) a - 2) 2 _. (2.4.10a) aT2 aT aT (l + v) ( + t Cot - 2 2 - t - 2 - 2 2 2 f Cot - 5 - t Cot 5 + I ) + 82 +'1 + 2Cot 5 + v(2 t t Cot t +' r Cot t - + t) - Q[* + 2tCot | - (1 + v + Cot2) + Cot g(2 - v + Cot2a) +"+'C Cot -'(l+ + + Cot2) + Cot (2 -v + Cot2])] *^ + -)- (2.4.10b) aT2 6T

-15Satisfaction of (2.4.9) also imposes natural boundary conditions on ~ and ro These conditions are, at 0 = 0 and H, 6:<[: t - G: + v(* r^ Cot - ] Sin t + o2[(4 + V) (1 - v) Cos22 Csc + - (4 +t) Sin g + vQ + i) Csc i - (t -') Cos i] = 0 (2.4.11a) 6: + v( + v( + ) Cot ~] Sin = 0 (2.4.11b) 5 j - + 2- 2 6 A 2 2+ +( _ - + 52+ ~ Cot 5 -~ ~2 - 2 Cot )+ 2[ + + 2v([ Cot + + Cot ~]} Sin > = 0 (2.4.11c)

CHAPTER 3 A LINEAR STUDY The linear motion of a closed spherical shell has received the attention of numerous investigators. A brief survey is contained in Chapter 1o Available results5 however, are not in a form convenient for the present nonlinear study. A linear analysis is therefore presented in this chapter. 301 The Linear Solution For sufficiently small deformations the second order terms in (204.10) may be neglected. The resulting expressions constitute the classical bending formulation of torsionless axisymmetric motion. The two linear partial differential equations governing free motions are' + ~ Cot * - t(v + Cot2 ~) - ~(1 + v) + a2F. + t Cot ( - t(v + Cot2) +' + Cot - ~ (v + Cot2 )] = 2 (3.11a) 2 ~T2 (i+v) ( C + V Cot 2 - 2 ) -2 [ + 2 * Cot t - R(+y-;+ +Cot2 ) + 4 Cot ~(2 - v + Cot2~) + f + 2' Cot 5 -' (1 + v + Cot2 ) + i Cot ~(2 - v + Cot2~)] = a (35. lb) rT2 -16

-17The displacement t is expanded in a series of Legendre polynomials of the first kind. Legendre polynomials of the second kind are singular at the poles and are therefore omitted from the expansion. The expansion is 00 - = Z a (T) Pn (Cos M) (3.1.2) n=o The displacement 4 is expanded in a series of associated Legendre polynomials of the first order, first kind. Again those of the second kind are omitted. Then1 00 4 = Zbn (T) Pn (Cos ) (3.1.3) n=l The assumed forms (3.1.2) and (3.1.3) satisfy the natural boundary conditions (2.4.11). Substitution of (l3.12) and (3.1.3) into (3.1.1) and repeated use of the Legendre relation Pn(Cos:) + Pn(Cos 5) Cot 5 + n(n + 1) Pn(Cos o) = O (3.1.4) yield D2ao + 2(1 + v) ao = 0 (3.1.5a) DT2 1 Pn(Cos ) is more commonly expressed with the change of variable x = Cos S whereby Pn(Cos ) =- (1 - x2) d [Pn(x)] = P(x) See, for example, Magnus and Oberhettinger(14) p. 53.

-18bn[1 - v - n(n + 1)] (1 + a2) - an{ 1 + v - (2[l - v - n(n+l)]} = D2bn n > 1 (31.5b) DT2 - bn(l + v)n(n + 1) - 2 an (1 + v) - O(an + bn) [n2(n + 1)2 (1 - v) n(n + 1)] D2an n > 1 (3.1.5c) DT2 The solution of (3.1o5) is ao = Ao Sin (woT + ao) (3.1.6a) an = Anm Sin (cunmT + canm) + Anc Sin (uncT + anc) n > 1 (3.1.6b) bn = SnmAnm Sin (conmT + canm) + bncAncSin (WLncT + cnc) n > 1 (3.1.6c) where Ao, Anm, Ancy, ao, om, and ancare arbitrary constants. The amplitude ratios are 6n = 1 + v - c2(1 - v) + o2 n(n + 1) (3.1.7a) [1 - - n(n + 1)] (l+12) +

-19bnc = 1 + v - c2 (1 - v) + a2n(n + 1) (3,1.7b) [1 - v -n(n + 1)] (l+a2) + oic Angular fre';ency Co0 is )Wo - [2(1 + v)]2 (3.1.8) Angular frequencies anm and onc are the real, positive roots of the frequency equation2 o4 - o2n [1 + 3v - C2(l - v) + n(n + 1) + v o2n(n + 1) + a2n2(n + 1)2] - 2(1 - 2) (1 + C2) + n(n + 1) (1 - v2) + ac2n(n + 1) (5 - v2) - 4 a2n2(n + 1)2 + c2n3(n + 1)3 = 0 (35.19) Solving for twice the square of the frequency yields 24 = 1 + 3v - a2(l - v) + n(n+l) (l+v a2) + 2n2(n+1)2 + [9 + 6v + v2+ 2 a2(3 + v) (1 - v) + 2n(n+l) (2v2+ 3v - 1) + n2(n + 1)2+ 2 a2n(n+l) (5v2+ 2v - 11) + 2 a2n(n+l) (9+4v) - 2 c2n3(n+l)3- 2 a4n(n+l) v(l - v) + a4n2(n+l)2(v2+ 2v - 2) + 2v a4n3(n + 1)3 + a4n4(n + 1)4]2 (31.10) 2 Equations (351.5) are noted to be analogous to a linear, two degree of freedom system; therefore, such roots are known to exist.

-20The positive sign in (3.1.10) gives the frequency cam whereas the negative sign gives oc. Substitution of (3.1.6) into (3.1.2) and (3.1.3) yields the solution of the linear problem, namely 00 f(S, T) = Ao Sin(cnoT + co) + Z [Anm Sin(onmT + Anm) n=l + Anc Sin(cnC T + anc)] Pn(CoS ) (3.1.11a) v(L, T) = Z [nmAnmSin(chm T + nm) n=l + 5ncAncSin(Dnc T + anc)] Pn(Cos I) (3.1.1lb) 3.2 Mode Classification The n=0 mode, pure radial motion with amplitude Ao, is called the breathing mode. Introduction of time by substitution of (2.4.2) into (3.1.8) gives the actual breathing mode angular frequency as O, [1 i) 2 E 12 (3.2.1) actual a L p(l - V) J a result first derived by Lamb(4). The distinguishing characteristics of the two mode sets are most readily studied by a quantitative presentation. The frequencies anm and Cnc are given in Table 1 for v =.3, 1< n < 10, and two values of a/h.

-21Table 1 Angular Frequencies a nm nc a/h n \ 20 200 20 200 1 1.975 1.974 0 0 2 2.723 2.722.703.701 3 3.636 3.635.841.830 4 4.598 4.596.916.881 5 5.578 5.575.988.906 6 6.565 6.562 1.081.921 7 7.556 7.552 1.206.931 8 8.550 8.546 1.369.939 9 9.546 9.541 1.572.946 10 10.542 10.536 1.813.954 For all values of a/h the lowest mode frequency is C2c. This mode is called the fundamental mode. A striking difference between the wnm and wnc mode sets is shown in Table 1. Angular frequencies for <nm modes are relatively independent of a/h; however, onc frequencies are sensitive to changes in a/h. This significant effect suggests a membrane-bending, or extensionalinextensional, mode classification.

-22For a membrane, or extensional, representation of deformation, first presented by Love(3), no consideration is given to the thickness of the shell. All resistance to deformation is assumed to be provided by stretching of the midsurfaceo A convenient measure of membrane behavior is the sum of the midsurface strains Em, Em = Eo~ + ~Eo (3.2.2) For a bending, or inextensional, representation of deformation, first presented by Rayleigh(l), the shell midsurface is assumed to be rigid3o Eb, a measure of bending behavior commensurate with Em, is the sum of the maximum strains occurring in the shell less the sum of the midsurface strains. From (2.3.4) and (2o357) the maximum strains occur at z = h/2 and their sum is i + En = co + o + h 1 - h) {(1 + - 1 {!max + (1 + co )[Ko Sin(-P) - 1 } (32.53) Hence, Eb is given as Eb = h 1- h {(1 + Eo) [ 1 ] + (1+oe )[ Sin(-) - ]} 2 2a a a (352.4) 3 Midsurface stretching must accompany all deformations of a closed shell (see Love(2), p. 542). Rayleigh was, therefore, able to include only open shells in his inextensional analysis. His limiting case of the opening angle - Hn was physically untenableo

-23The order of magnitude of Eb is not changed by the neglect of Eoe, E', and h/2a with respect to unity. With this, E - h. { - 1 + [ ~ Sin(0-) - 1] } (3.2.5) 2 a L a Em and Eb are obtained in terms of the displacements ~ and 4 by substitution from (2.o53). The linearized expressions are Em = - 2 + + * Cot e + ~ (3.2.6a) Eb = h (2 5 + 5 Cot e + ~) (3.2.6b) 2a The ratio IEm/Ebl is an indicator of mode behavior. A large ratio indicates predominately membrane behavior; a small ratio predominately bending behavior. For every n > 1 solution (3.1.11) yields the following indicator for modes associated with either anm or Wnc: Em =2 + n(n + 1) 5n o 2a (3.2.7) I E2 - n(n + l) h in which 5n is 5n for an om mode and 5nc for an cnc mode. Numerical values of |Em/EbI for the am and Onc mode sets are given in Table 2 for two values of a/h.

-24Table 2 Mode Indicators _EM/Eb_ a/h = 20 a/h 200 n |unm W nc | )nm | nc 1 00 -- 2 56.4 3.80 564 38.0 3 4005 2.16 405 21.6 4 36.0 1.32 360 13.2 5 34.1.88 341 8.8 6 33.2.64 332 6.4 7 32.5.48 325 4.8 8 32.0.536 320 3.6 9 31.8.28 318 2.8 10 31.7.24 317 2.4 ~a a o o 0 o o i 53008 0 308 0 For a/h = 20 the relatively large values of IEm/Ebl indicate that the cnm mode behavior is predominately membrane for all n. However, the indicator values for the same a/h show cnc modes are both membrane and bending in nature. Their behavior is changed from predominately membrane to predominately bending in the finite range of n's given.

-25From (35.27) the mode indicators are directly proportional to the radius-to-thickness ratio. The value a/h = 20 is about the minimum for athin shell theoryo All indicator values may be made arbitrarily large by increasing a/h, a point emphasized in Table 2 by comparison of a/h = 20 with a/h = 200. The relatively large values of jEm/Ebi for u>m modes at a/h = 20 and their subsequent increase as a/h increases substantiate a membrane mode classification for this seto The corresponding increase of hnc mode indicators causes the transition from membrane to bending behavior to occur at higher values of no Conversely, for any n an cnc mode may be made predominately membrane in behavior by choice of a sufficiently large value of a/h. Classification of the cn set as bending modes, as their dependence on a/h suggests, would be equivocal. Therefore, the Cwnc modes are classified here as composite modes of vibration.4 An extensional treatment of the problem is obtained by setting a2 0. Frequency equation (351,10) then becomes 4 Kalnins(8) studied the modes of vibration by use of strain energy coefficients defined as ratios of bending strain energy to total strain energyo This approach is not unlike that of IEm/Ebi mode indicators as presentedo Modes here classified as composite modes were classified as bending modes and were associated with predominately bending strain energyo This association is dependent on a/h and n as is indicated by (35o27) and Table 29 and, therefore, must be qualifiedo See McIvor and Sonstegard(15)o

-262 = 1 + 3v + n(n + 1) + [9 + 6v + v2 + 2n(n + 1) (2v2+ 3v - 1) + n2(n+l)2] 2 (3.2.8) It follows that as n - o5: cnnm -~00 (352. 9a) wnc - [1 - v2] 2 (3.2b9b) Membrane mode frequencies increase without bound, as expectedo However, composite modes constitute an infinite set in a bounded frequency spectrum, an untenable result.6 The increasingly membrane behavior of composite modes as a/h increases was noted above. The degeneracy of this set to the paradoxical extensional analysis set istherefore, not surprising and further substantiates a composite, and not a bending, classification.7 5 n cannot increase without bound in keeping with the neglect of shear deformation and rotary inertia in this formulationo Consideration of this case is of interest, however, in providing a better understanding of mode characteristics. 6 This is the lower mode set of paradoxical behavior referred to in the Introduct ion. 7 Association of degenerate " bending " modes with the paradoxical extensional modes was made by Kalnins(8).

-27The relative independence from a/h of membrane mode frequencies, Table 1, and the results obtained from the extensional formulation indicate that (3.2.8) provides an adequate means for their determination. However, composite mode frequencies may be determined by an extensional analysis only for sufficiently large values of a/h and/or small values of n. Both factors must be considered before use of (302.8) is justified. 353 Response to an Arbitrary, Symmetric Radial Velocity A closed spherical shell is suddenly enveloped by a pressure impulse symmetric about a polar axis. For a loading time of sufficiently short duration the result of the impulse is an arbitrary, symmetric initial radial velocity distribution. This velocity distribution is expressed in a series of Legendre polynomials of the first kind as 00 W = Z vn Pn (Cos 5) (353-1) aT n=o where w = a - r Expression (33..1) imparts a rigid body translation along the negative polar axis. This motion is eliminated by use of a reference frame with its origin at the shell inertial center moving with an appropriate velocity. This matter is considered in detail in the Appendix, The velocity distribution with respect to the inertial center contributes only to the shell deformation. It is

-28W = Vo + v -1 Vn Sn P1 (Cos E) Sa~t' 4 n=1,3... 00 + z v Pn (Cos ) (335.2a) n-=2,3,,.. 00 au = 1 E vn Sn P1 (os ) (3.2b);t 4 n=1,3,... where vw/at is the radial velocity distribution, au/3t is the meridional velocity distribution and Sn = (2n).' 4 + 1 n 4 22n(n')2 4-n2 1 2n-1 4-(n-2)2 + 1-3 n(n- 1) 4 + 1-2 (2n-1) (2n-3) 4-(n-4)2 The brackets contain n + 1 terms. Motion resulting from the impulse is determined by imposing (3.3.2) as initial conditions on solution (351.11). The problem is thus formulated as 00 i(ST) = AoSin(coT + Go) + Z [AnmSin(cnMT + anm) n=l AncSin(cncT + nc))] Pn (Cos ) (51.11a)

-2900 ~(5, T) = Z [5nm Anm Sin(n T + cnm) n=l + 6nc Anc Sin(ncl T + UCnc)] Pn(Cos 5) (35.11b) ~(~, 0) = 0 (3.35.a) ~(S, 0) = 0 (3,3.3b) B (~, 0) =1 aw (33.3.c) aT c aT _ (i, 0) =1 a (33..3d) T7d ac dT where E 2 = I p(l-v2) The arbitrary constants in (3.1.11) are now found to be ao = nm = cn = 0 (533.4a) A0 = AO = vO (3.3.4b) 00 Alm = V1 Z Vn Sn (5.3.4c) n=1,3,5... c 1Im Ale = 0 (3.3.4d)

-30Anm = nc Vn n >2 (3.3.4e) ~nm(bnc - 5nm) c Anc nm n n > 2 (5.3.4f) cnc(Snm - bnc) c The solution is of the form:(~,T) = Ao Sin (oT + AlmSin Wolm Pl(Cos ~) 00 + Z [Anm Sin nOrm + ncSin WncT] Pn(Cos ~) (353.5a) n=2 (,yTr) = zlm Sin wlmT P1(Cos 5) 0 _ __ + Z [ AnmSinmT + AncSinfncT] Pn(Cos 5) (353.5b) n=2 am and Anc are the radial amplitudes of the nth membrane mode and the nth composite mode respectively. Similarly Anm and Qnc are the meridional amplitudes of the nth membrane mode and nth composite mode respectively. For all n > 2 the membrane mode radial and meridional amplitudes, Anm and A., are compared to the composite mode radial amplitude nc. Then | sm = 5nc Wnc (353.6a) nc nmn nm

-31_nm| = 5nc nc (3.3.6b) Xnc ^'nm These two amplitude ratios are plotted in Figure 4 for 2 < n < 10, v =.3, and a/h = 20. Values shown are negligibly changed for a/h = 200..12.10 Ratio \ | nm.o6 \nc.04.02 2 3 4 5 6 7 8 9 10 n Figure 4. Plot of Amplitude Ratios

-32Figure 4 indicates that for all appropriate a/h and n > 2 the motion is predominately due to vibration in composite modes. Hence, it is assumed that the membrane contributions to: and 4I for n > 2 are negligible and that the response to velocity distribution (553.2) may be expressed as 00 t(~T) = Vg Sin LOT + [vl- 1 Z VnSn] Sin ClmT Pl(Cos ) c — __ n=1,3,. o Co c lm 00 + Vn bnm Sin WncT Pn(Cos ) (3.3.7a) n= —2 "nc bnm - bnc 00 lr(ST) = I vnSn Sin wlm Pl(Cos ) n=l,3,.. 00 + Z v 5nm 5nc Sin ancT Pn(Cos ) (3.3.7b) n=2 -- nt-Si ccnc nm - nc Expressions (3.3.7) do not satify initial conditions (353.3) but an approximate set which is (~, o0) = 0 (3.3.8a) (~, o0) o0 (53.38b) 00 Z. vnSnPn(Cos a: (e, 0) = vo + vi - Z vnS P,(Cos n) aT c n=1,3,5.. c

-33+ Z Vn nm Pn(Cos ) (3.3.8c) n=2 C 5nm -5nc a (e, 0) = 1 Z vnSn Pl(Cos A) aT 4 n=l,3,... C 00 + Z vn 5nm 5nc Pn(Cos 5) (3.3.8d) n=2 - c 5nm-5nc Since wnm frequencies are higher than the fundamental frequency w2c (see Table 1) the maximum amplitude of the fundamental mode is attained after those of the membrane modes. Hence, the displacement comparison of Figure 4 is not initially valid and the entire solution (5.3.6) must be retained. The period of the fundamental mode Tf is Tf = 2n (3.3.9) cW2c All composite modes will have attained their maximum amplitudes within a time II/2ca2c; thereafter, the transition to predominately composite mode deformation will have occurred and approximations (33..7) and (553.8) may be used. Since bending effects in the shell are predominately due to composite mode motion an accurate determination of flexural stress may be obtained from the approximate solution. The neglect of membrane mode contributions, however, prohibits determination of membrane stress from this solution.

CHAPTER 4 RESPONSE TO A NEARLY UNIFORM RADIAL IMPULSE The basic response of a closed spherical shell to a uniform radial impulse is the breathing mode of motion. Small non-uniformities in the impulse may cause this response to be dynamically unstable. The results of Chapter 3 indicate that composite mode deformations will predominate any deviation from the basic response. Breathing mode stability is therefore studied with respect to perturbations of composite modes of motion. Stability criteria are thus determined by the character of nonlinear coupling between linear modes of motion, i.e. by the possibility of autoparametric excitation. 401 Formulation When an uniform impulsive pressure of short duration suddenly envelopes a closed spherical shell, a radial velocity vo is imparted to every shell element. In terms of the non-dimensional quantities previously introduced the initial conditions are (, 0) = 0 (4.1,la) 4(5, O) = (4.1.b) a (e, 0) = o (4.1.1c) aT c at (e, 0) = o (4.1.1d) aT -34

-35where 1 2 l p(l - v)2 Governing equations of motion (2.4.10) and initial conditions (4.1.1) are satisfied by a harmonic breathing mode,:(SJ T) = v0 Sin M0T (4.1,2a) Co0 (Sg v) =.0 (4.1.2b) where {o = [2(1 + v)]2 If the impulsive pressure is slightly non-uniform, this basic response may be dynamically unstable. This possibility is to be studied by applying an axisymmetric perturbation to the initial velocity, giving initial conditions t(S, 0) = O (4.1.3a) t(9, o) = 0 (4.1.3b) 00 a (e, 0) = vo [1 + Z n Pn(CoS e) (4.1.3c) - - n=l'T ("C

-366a ( 0, o) = o (4.l.3d) aT where n << 1 The basic response will be termed stable if the amplitudes of the deviation from the breathing mode remain of the order of magnitude of eno Conversely, if these amplitudes grow to magnitudes of higher order, the basic response will be termed unstable. 4o2 The Approximation An exact solution of nonlinear equations (2.4.10) satisfying initial conditions (401o3) has not been found. An approximation is feasible by utilizing results of Chapter 35 The perturbation in the uniform initial velocity, equations (4.ol,3), constitutes initial conditions for the response to an arbitrary, symmetric radial velocity distribution, equations (35302). From Section 5o3 the corresponding linear response is predominately composite mode motiono This result is extended to the slightly larger motions anticipated here, The perturbation of the basic response is limited to a superposition of composite mode motion on the breathing mode. The investigation is thus a study of the stability of the breathing mode with respect to composite mode deviation only.

-37An approximate set of initial conditions is now formulated. Conditions (4oo13) are expressed in the form of (35308) as 5(, 0) = 0 (4 2. la) r(S, 0) = 0 (4,2.Ib) (^ 0) V2 1 + (+l+ T) PJ(Cos ~) aT c 00 + Z En Pn(CoS )(4.2.lc) n=2 d3 (5, 0) = v Pi(Cos,) _r c. + EnbncPn(cos ) (4.2, d) The quantity E represents the effect of having chosen a reference system so as to remove rigid body motion imparted by (4.1.3). T may be determined by the method. outlined in the Appendix. The contribution of C to the response will be of an order of magnitude of a perturbation and is neglected. The approximation implies a perturbation free of a first mode (membrane) component; hence, E1 is also neglected. The resulting initial conditions are (, 0) = 0 (4.2.2a) 4(~, 0) = 0 (4.2.2b)

-38(~, 0) = Vo 1 + E EnPn(Cos ) (4.2.2c) aT C _ 00 5a (e, 0) = vo Z en 5nc Pn(Cos 4) (4,2.2d) --- ~ — n=2 EAT c The assumed solution consists of the breathing mode together with composite modes. From equations (3.357) the appropriate representation is 5(5, T) = ao(T) + an(T) Pn(COs 5) (4.2.3a) n=2 00 ~(r, T) = z an(T) 8nc Pn(Cos 5) (4.2.3b) n==2 4.3 Stability of the Breathing Mode The approximation of Section 4,2 allows a study of the stability of the breathing mode using the energy formulation of Section 2.4. The Lagrangian L is (2.4.3) minus (2.4.7). The neglect of 2~ with respect to unity gives II L = n a2E h 2 \+ 2 2- 2_ 2Cot2 1 - V2 \ oT FT I-v dS/ v T - 2(1 + v) (2_ t Cot e - 4 + 52* Cot + + 2~) + 22_ 2f + 2 2_ t2+ Cot + 2 ~ r2Cot2t - 2v(* 4r Cot e - 4 424 _ 2 CtJr Cot e + 2,2 -_ 2 2

-39+ ~ ~2c Cot 2) - a2[r2+ ~2Cot25 + 2 + 2r + 2 ~ Cot2S + +2 + 2Cot25 + 2v(4 r Cot + + ~ Cot t + /r Cot g + ~' Cot ~] Sin 5 d 4 (4.3.1) Representation (4.2.3) is now substituted into the Lagrangian. Initially an(n > 2) is of an order of magnitude of a perturbation, i.e. an << ao. Lagrangian terms of order three will therefore be retained only if ao occurs to at least the first power. After integration the Lagrangian valid for initial motion may be expressed L =-lla2E h 2 Da2 +2 ( Dr [l + n 1-v2 DTr 2n+l 2 2 =2 Z 5ncan n(n + l) [n(n + l) - l + v] _ 2n+l 00 2 an (1 + 5nc)2 n(n + 1) [n(n + 1) - 1 + v] n=-2 2n+l 00 + 8(1 + v) ao Z 6nc a2n n(n + 1) 2n+l

+ 4 ao Z nc an n(n + 1) [n(n + 1) - 1 + v] n —22n + 1 + 2(1 + v) ao0 a2 n(n + l) (1 - bnc) (43.2) n=2 -- - 2n+l Quantities ao and an (n > 2) are identified as generalized coordinates. Motion is, therefore, governed by Lagrange's equations expressed as D roL' -L =0 (4.3.3a) DTi a TDaOr ao DT D aL - aL =0 n>2 (4.3.5b) DT Danj aan DDT Introducing L from (4.35.2) into (4.5.53) yields the equations governing initial motion. They are D2a + 2(l + v) ao - an n(n + 1) 1 + v DT2 n2 2n+ 2 + 2 5nc(l + v) + 52c[n(n + 1) - 3-v = 0 (4.3.4a) 2 J D2an [1 + 2nc n(n + 1)] + {2(1 + v)[1 + +nc n(n + 1)] DT2

-41+ n(n + 1) [n(n +1) - 1 + v] [82n + c62(1 + 5n)2] an - (1 - nc) n(n + 1) (1 + v) + 2 2nc n(n + 1) [n(n+l) - 1+v] + 4 6nc(l + v) n(n + 1) aoan = 0 n > 2 (4.3.4b) Since a2 << a initially, equation (4.3.4a) may be rewritten as D2ao + 2(1 + v) ao = 0 (4.5.5) DT2 the solution of which is the breathing mode motion ao = vo Sin LO (4.3.6) c'Wo The constants of integration have been evaluated for initial conditions following from (4.2.2) and (4.2.3a). Solution (4,3.6) is now substituted into (4.3.4b). A change of independent variable T = COT gives D2a + (n n Sin) an = n 2 (4.5.7) D72

-42where n = _1..... F2(l+v)[l + 5nc n(n + 1)] ([D+o 5c n(n+l)] L + n(n+l) [n(n+l) - 1 + v] [52 + Co(l + 5nc)2] (4.3.8a) n = VO 1 ( (1 - 6nc) n(n + 1) (1 + v) cw' 1+52c n(n+l) L + 2 52 n(n + 1) [n(n + 1) - l+v] + 4 bnc (1 + v) n(n + i (4.3.8b) Equation (4,-38a) may also be expressed as On = nc (4.3.9) thereby emphasizing the character of (4.3.4b). The linear portion of the equations represents motion in composite modes of vibration. The nonlinear, or ao an terms represent coupling of the breathing mode with composite modes, Equations (4.3.7) are Mathieu equations. Thus deviations from the breathing mode are Mathieu functions and their growth characteristics

-43are determined from the well known Mathieu stability diagram shown in Figure 5 (see for example Stoker(l6) p. 202)1. Whether or not a particular deviation amplitude an remains of an order of magnitude of a perturbation is readily deduced by location of the point (sn, [n) on Figure 5. If the point falls in an unshaded region or on a boundary curve other than ~n = 0, the amplitude an will increase exponentially and the breathing mode is unstable. Such a deviation will be called a critical deviation2. Otherwise an will not change order of magnitude and the breathing mode is stable. The two bounding curves for a particular unstable region intersect the positive Q axis at Q = k2/4 (k = 1,2,3,...). In the discussion to follow a region will be denoted by this intersection, such as region.25. 1 Breathing mode stability for the cylindrical shell was also found to be governed by Mathieu equations by Goodier and McIvor(l). Their study includes analog computor results which demonstrate the stable and unstable growth characteristics of Mathieu functions. These results will not be reproduced. 2 Floquet theory, see Stoker(16) p. 193, predicts that a critical deviation will grow without bound, an impossible consequence for the fixed, finite energy input of the initial pulse. Terms in (45.34) predominating after growth starts have been neglected. Their effect is examined in Section 4.4,

-44The membrane stress generated by the breathing mode is co = E vo Sin C0T (4.3.10) - v ccuo For the material to remain elastic, ao must remain less than the yield stress. Thus for a material of yield stress ay the maximum radial velocity imparted by the initial pulse may be expressed as (vo mx y Wo (1 - v) (4.3.11) max About the largest ay/E value of practical interest is 5(10)3. The maximum value of vo/c for this case is denoted Vmax and is Vmax = 5.646 (10)'3 (4.3.12) A pertinent portion of the Mathieu stability diagram for Vmax is shown in Figure 6. The a/h curves have significance only at intersections with the dashed, nearly horizontal lines denoting integer values of n. Those intersections occurring in unshaded regions or on their boundaries indicate deviation growth and breathing mode instability. For example, for a shell of a/h = 100 intersections occur in region 1.0 at n = 21 and 22. Deviation amplitudes a21 and a22 therefore increase exponentially and the breathing mode is unstable. For a/h = 50 no intersections occur in an unshaded region or on a boundary

-45and the breathing mode is stable. By interpolation between curves, al2 will cause instability at about a/h = 355 It is notable that no instability can occur for n < 9, a/h > 20. In Figure 6 critical deviationsoccur in region.25 for sufficiently large a/h. At a/h = 150 deviations 11 through 17 will cause breathing mode instability (other critical deviations occur outside region.25, a matter to be discussed later). More and more critical deviations are included in region.25 as a/h is increased. Referring now to Figure 5, regions 1.0,2.25,... become increasingly narrow strips for decreasing Ap. Should a point (On, n) fall within such a strip an unstable solution is indicated. The corresponding deviation is not classified as critical, however, because an unjustifiably precise knowledge of the system constants is presumed in order to locate (tn,.n)o Furthermore, the inevitable presence of damping cuts off portions of the unstable regions which border on the f axis, see Bolotin(12), ~ 9. From Figure 6 it is clear, however, that critical deviations may occur in regions 1.0,2.25,o.. for sufficiently high n and large values of a/h. Therefore, as a/h increases consideration must be given to occurance of critical deviations in more and more regions. Lubkin and Stoker(17)showed that in the presence of an arbitrarily small amount of damping a value Qo exists for which the region Q > Qo, < Q is entirely stable. For large n,

-46tn Vo 1 (4.3.12) ~n cabo Q2 n2 Hence, for all a/h (expect the extensional case) some n exists above which critical deviations will not occur. It is interesting to note that the Mathieu diagram, Figure 5, is symmetric about the Q axis. Consequently the stability criteria for an inward initial velocity remain unchanged for an outward initial velocity, The stability diagram for.6 Vmax, or, what is equivalent, for the maximum value of vo/c when oy/E = 3(10)-, is shown in Figure 7. The effect of this reduced initial velocity is apparent from comparison of Figures 6 and 7. Since On is independent of vo/c, the ordinate In of all points is simply multiplied by the factor.6. The stability diagram for any other vo/c may be determined accordingly. Referring to examples given on p.44 for Figure 6, reduction of Vmax by.6 leaves only a21 as a critical deviation at a/h = 100. a/h = 50 still gives a stable breathing mode. al2 will no longer give instability for any value of a/h; in fact, no critical deviations can now occur for n < 14. As before, consideration of region.25 and of the portion of regions 2,25,... with significant width is necessary for correspondingly higher a/h and n.

-47Two striking characteristics of stability behavior are noted in Figures 6 and 70 The first is that relatively high order modes are excited. At Vmax, ay/E = 5(10) 3, no modes of order less than 10 cause instability. For smaller values of Oy/E this lower limit is increased. The second is that for thinner shells instability will almost always occur. 4.4 Long-term Behavior The Mathieu equations governing breathing mode stability are valid only for initial motion. Only those third order terms containing ao to at least the first power were retained in Lagrangian expression (4.302). As indicated on p. 43, terms must be added to this expression after deviation growth starts or a paradoxical energy buildup will occur. The ratio of meridional to radial displacements 5nc is much less than unity for the composite modes which may be strongly excited. Significant third order contributions to the Lagrangian (4.3.1) are thus assumed to come solely from the term (1 + v) 5 ~2 since all other such terms contain 4 and/or t. Using representation (4.2.2), the Lagrangian is now given by (4.3.2) plus the term +1 2 a2E h <(1 + v) z anPn(x) [ amPm(x d (4.4.1) 1.V2 \n -1 -1 -~~~~~~~~~~~~~~_

-48Subsequent substitution into Lagrange's equations gives D2ao + -0oao l+v Z n(n + 1) a2 = 0 (4.4.2a) ----- n=2 DT2 2 2n + 1 D2an+ n + 2 n @o n ] an DT2 DT vo - (l+v) f amPm(x)] aP(x) dx = 0 can [ m=2 [qZ~2 aqpl(x) ian -1 L J L - J (4.4.2b) The last term in each of equations (404,2) is significant only for those n denoting critical deviations. In (4,4.2a) this term gives a reduction of breathing mode amplitudes as deviations grow. This term, containing the integral, in (4.4o2b) will give two forms of contributions. The first form will contain products of the a's with equal subscripts; the second will contain products with unequal subscripts. The latter form represents interactions between critical deviations. The relatively high orders of Legendre polynomials involved prohibit the practical evaluation of the integral term in (.4.42b). The effect of this term has been determined at Vmax for one critical deviation, al0, and at 1o4 Vm for two critical deviations, a9 and a 10 The max 9 0 numerical solutions of (4. 42) showed the integral term to be of negligible consequence in these cases,, 3 Since the solutions are qualitatively similar to other cases to be presented later they are not pursued in detail,

-49This observation lends credence to an approximation which facilitates the study of long-term behavior for higher values of n. It is assumed that the second order terms which would result from the retention of the integral expression are negligible with respect to the linear term 2 a The resulting equations governing long-term motion, 00 D2ao + 2oao - l+v Z n(n + l) a2 =0 (4.4.3a)..... n=2 ----—. DT2 2 2n + 1 D2an + 02o [n - c)o tn ao an =0 (4o4.3b) DT2 V have been solved numerically for several cases of instability. A representative long-term behavior is shown in Figure 8 for a/h = 100, vo/c = Vmax o Critical deviations 21 and 22 grow in amplitude, as was predicted in Figure 6, and are seen to enter a cyclic energy exchange with the breathing mode. The results for a reduced impulse of.6 Vmax is shown in Figure 9. As predicted in Figure 7 only deviation 21 is critical and again a cyclic energy exchange with the breathing mode is observed. An interesting feature of the energy exchanges which was observed in all solutions obtained is demonstrated in Figures 6 through 9. From Figure 8, the time required for energy exchange to occur with a22 is less

-50than that with a21 l By comparison of Figures 8 and 9, the time required for exchange with a21 is significantly increased for the reduced impulse.6 Vmax~ Noting the corresponding positions of these deviations on the stability diagrams, Figures 6 and 7, it is observed that the time required for significant amplitude growth is less for increased An and for increased distance from the boundaries of the instability regions. The extent to which energy is transferred from the breathing mode to critical deviations is of considerable interest. The total energy Ti imparted to the shell by the impulsive pressure is Ti = 21 a2p h v (4.4.4) The potential, or strain, energy stored in the shell at any time T is V = 211 a2E h 2^ a2 - v2 L 00 + 1 (Qna2 co naoan) (4.45) n=-2 2n+l In (4.4.5) the terms corresponding to the second order terms neglected in (4.43.) have been deleted, Should all the energy be contained in critical deviation N, the potential energy at some time will be VN = 21a2Eh o2 QN aN (4.4.6) 1-v2 2N+1

-51Equating VN to Ti gives the deviation amplitude AN for complete energy exchange. It is 2 AN = N + 1 v (4.4.7) L N cab For N = 22 and Vmax, A22 = 6.25 (vo (4.4.8) For N = 21 and.6 Vax, A21 = 6.49 vo \ (44.9) From the maximum deviation amplitudes attained in Figures 8 and 9 it is concluded that the energy exchanges are essentially complete. In Figure 8 the different time intervals for deviation growths result in a significant energy content in but one deviation at some instant. In both Figures 8 and 9 the amplitude a0 of the breathing mode is seemingly significant at times of essentially complete energy exchange. Its energy content is, however, small. The total energy was calculated as a check at each step of the numerical analysis. The energy content of the system remained constant and thus the terms neglected in the approximation do not significantly contribute to the total energy.

-52The sum of the maximum flexural stresses in a case of complete energy transfer is (a + ar) = E h [N(N + 1) - 2] AN (44.10) flex 1-v 2a By substitution from (4.4.7), (ao + =4 E h N(N+1) - 2 2N+l 2 vO (4.4.11) flex 1Ov 2a L N J cOo flex ^-v 2a -L In terms of ao, the stress produced by the basic breathing mode, (4.4.11) becomes ( + a, - h [N(N + 1) - 2 2N + 1 2 CO (4.4.12) flex 2a- ON For the two examples of essentially complete energy exchange given above it follows that i(a + a)l 15.75 ao (4.4.13a) L =22 [c + +ra)flexj = 14.38 ao (4.4.13b) N=-21 Hence, the stresses caused by breathing mode instability may be considerably in excess of those predicted by consideration of the breathing mode alone.

-53The neglect of higher order terms in (4.4.2) limits numerical consideration of long-term behavior. The limitation may be observed in the potential energy expression (4.4.5). For sufficiently small in and large An this expression may become negative, a physically untenable result. Such was the case in attempted numerical solutions for Vmaxg a/h > 150, i.e. as modes grow the contributions of the integral term in (4.4.2b) become significant. However, for all these attempted solutions the deviation growth was initially as predicted by the Mathieu stability diagram and the time required for significant growth was again characterized by the location of the deviation with respect to An and the bordering curves. 4.5 Conclusions The closed spherical shell under nearly uniform axisymmetric impulsive loading may undergo autoparametric excitation. The. parametric effect is generated by nonlinear coupling between the breathing mode and certain composite modes of motion, i.e. by the interaction of the membrane stress with the flexural curvature. The breathing mode seldom excites relatively low order composite modes. High order modes are almost always excited for thinner shells. In those cases of long-term behavior for which numerical solutions have been obtained an essentially complete, cyclic energy exchange occurs between the breathing mode and the parametrically excited modes. Resulting displacements and stresses are in excess of those predicted by a consideration of the breathing mode response only.

-54A thorough numerical study of long-term behaviors for very thin shells (a/h > 150) has not been made. This limitation is due to extensive computation required in the evaluation of an integral expression involving relatively high order Legendre polynomials. Initial behavior for these cases is not affected by the limitation; therefore, the same general conclusions are implied although quantitative results are not available. Further pursuit of long-term behavior for very thin shells does not seem appropriate for this axisymmetric formulation. A realistic appraisal of the physical situation for the relatively large values of n entailed for such shells indicates that asymmetric effects would be significanto

.55-X 7 6 ~~F1i e MtiuSalt..

-56a/h = oo o o.500 - a/h-m =X1000 4 -47,_~ _ 9_; l 18~46 4J4 200 43 2.9 <1~~~~~~~~~~~~~~42..... 41.....339 ______ 0 36 3 34 0 31 - - 32 30 Figure 6. Stability Diagram forV- 28 - 20 Figure 6. Stability Diagram for Vmax

-57a/h = o 1000 500 - - -rlE~l\\> / _ 60 59 -58 -57 -- - 56 55 K- = 54 3 -- - - - 1 -6 — 53 - -- =-452 F51 " 5 49 28 47 45 - 4 200 0 1 2 3

,-*1 I — To 3.88 a. 7o/ f 50 100 J~50 200 250 6 -H T21 3.86 2 a21 vo/wao 50 -oo / 0g 207 ^ 250-2 -6~-I~ j*~T22 35.64 4 a22 6 150 /| /5 2 /%5 -2 / N -i' ^l ~f' -4 -6 Figure 8. Long-term Behavior for Vm, a/h = 100, e =.05/n5

u/;o-' = "U'00o = 1 q/ LIN _-Q 1H L — i oo0 oOt | 00oo 00o 00T N.................................oo.

-60APPENDIX In Section 3o3 an impulsive pressure acting on a spherical shell is assumed to impart a velocity distribution aw/at symmetric about the polar axis. This distribution can be represented as 00 aw = z vn Pn (Cos ) (Al) at n=O The Legendre polynomials of the first kind may be expressed in cosine expansions (see, for example, Magnus and Oberhettinger(14), p. 50 )o The expansion is Pn(Cos () = (2n). Cos n + 1; n Cos (n-2) 22n(n.)2 1 2n-1 + 13 n(n-l) Cos(n-4) r +o.. 1 (A2) 1?2 (2n-1)(2n-3) The brackets contain n+l terms. The momentum along the polar axis Mz is M, =/ p h aw d(Area) (A3) J at Area Substitution of (Al) and (A2) into (A3) and integrating yields 00 MZ = -J p a2 Z VnSn (A4) n odd where Sn = (2n)' 4 + n _ 4 22n(n')2 4-n2 1 2n-1 4-(n-2)2 + 1'3 n(n-l) 4 +.. 1*2 (2n-1)(2n-3) 4-(n-4)2

-61Again the brackets include n+l terms. Expression (A4) corresponds to imparting to the center of mass of the shell a velocity V along the polar axis, 00 = -t E VnSn (A5) n odd The distribution aw/at is now considered relative to a shell whose center of mass is moving with velocity Vo This relative distribution will contribute only to shell deformation. It is 00 00 a = vO + [v - VnSn] P(Cos ) + Z VnPn(Cos ) (A6a) 6t n odd n=2 00 au = ( S vnSnPl(Cos 5) (A6b) at n odd where w and u are radial and meridional displacements respectively.

-62BIBLIOGRAPHY lo Lord Rayleigho " On the Infinitesimal Bending of Surfaces of Revolution", Proceedings of the London Mathematical Society, XIII (1888), po4 2, Love, AoEoHo The Mathematical Theory of Elasticity, Dover Publications, New York, 1944. 3o Love, AoEoHo " The Small Free Vibrations and Deformation of a Thin Elastic Shell", Philosophical Transactions of the Royal Society, 179A (1888), p. 491. 4. Lamb, Ho " On the Vibrations of a Spherical Shell", Proceedings of the London Mathematical Society", 14 (1882), po 50~ 5o Silbiger, Ao " Free and Forced Vibrations of a Spherical Shell", ONR Report U-106-48, December, 1960. 6. Baker, WoEo " Axisymmetric Modes of Vibration of Thin Spherical Shell", The Journal of the Acoustical Society of America, 33 (1961), p. 1749. 7. Naghdi, PoM. and Kalnins, Ao t" On Vibrations of Elastic Spherical Shells", Journal of Applied Mechanics, 29 (1962), p, 65. 8~ Kalnins, Ao " Effect of Bending on Vibrations of Spherical Shells", The Journal of the Acoustical Society of America, 36 (1964), p. 74, 9. Lord Rayleigh "' On Maintained Vibrations", Philosophical Magazine, 15 (1883) Series 5, p. 229. 10o Minorsky, No Nonlinear Oscillations, Do Van Nostrand Company, Inc., Princeton, 1962, 11o Bolotin, VoVo The Dynamic Stability of Elastic Systems, HoldenDay Inc., San Francisco, 1964. 12. Goodier, J.No and McIvor, IoKo " The Elastic Cylindrical Shell Under Nearly Uniform Radial Impulse", Journal of Applied Mechanics, 31 (1964), p. 259. 13o Hildebrand, FoBo Methods of Applied Mathematics, Prentice-Hall, Inc., Englewood Cliffs, 1956,

-63 - 14. Magnus, W. and Oberhettinger, Fo Formulas and Theorems for the Functions of Mathematical Physics, Chelsea Publishing Company, New York, 1949. 15. McIvor, IoKo and Sonstegard, DoA. "Discussion of'Effect of Bending on Vibrations of Spherical Shells' "' The Journal of the Acoustical Society of America, 37 (1965), po 931. 16. Stoker, JoJo Nonlinear Vibrations, Interscience Publishers, Inc., New York, 1950 17. Lubkin, So and Stoker, JoJo " Stability of Columns and Strings Under Periodically Varying Forces", Quarterly of Applied Mathematics, 1 (1943), po 215.