THE UNIVERSITY OF MICHIGAN INDUSTRY PROGRAM OF THE COLLEGE OF ENGINEERING DYNAMIC STABILITY OF A SHALLOW CYLINDRICAL SHELL Carl H. Popelar 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 October, 1965 IP-721

Doctoral Committee: Associate Professor Ivor K. McIvor, Chairman Professor Samuel K. Clark Professor Joe G. Eisley Assistant Professor Roger D. Low

ACKNOWLEDGMENTS It is a pleasure to express my appreciation and gratitude to my Chairman, Professor Ivor K. McIvor, for his excellent guidance and suggestions throughout the course of this work. I also wish to thank the other members of my committee for their time and interest. Appreciation is also expressed to the University of Michigan Computing Center and The Ford Foundation for the use of the computing facilities. Last but not least, I wish to thank my wife, Joyce, for her patience and understanding throughout my graduate program and for her help in preparing the preliminary drafts of this thesis. ii

TABLE OF CONTENTS Page ACKNOWLEDG MENTS....................................................... ii LIST OF FIGURES....................................................... iv N'Ot~tCLATURE........................................................... v INTRODUCTION.......................................................... 1 CHAPTER I. MATHEMATICAL FORMULATION OF THE PROBLEM................. 3 A. Strain-Displacement Relations....................... 3 B. Equations of Motion................................. CHAPTER II. LINEAR VIBRATIONS..............11.............-. 1i A. General Solutions of the Linear Equations........... 11 B. Normal Mode Shapes of a Simply Supported Shell......, 14 C. Response to an Arbitrary Initial Radial Velocity.... 22 D. Example of a Specific Initial Radial Velocity.,.,.. 26 CHAPTER III. SYMMETRIC DYNAMIC BUCKLING............................ 29 A. Derivation of Governing Equations......3.......... 30 B. Method of Solutions and Results.................... 3 CHAPTER IV. EXCITATION OF ASYMMETRIC MODES.......................... 37 A. Derivation of Governing Equations..................37 B. Initial Growth of the Asymmetric Modes.............. 40 APPENDIX A. ORTHOGONALITY OF THE RADIAL MODES.......................64 APPENDIX B. EVALUATION OF THE INTEGRALS cl THROUGH cg.......... 66 REFERENCES............................ 71

LIST OF FIGURES Figure Page 1. Coordinate System and Shell Geometry.......................... 52 2. Location of a Point Before and After Deformation.............. 95, 3. Symmetric Mode Shape Zo........................................5 4. Symmetric Mode Shape Z1. 4 5. Symmetric Mode Shape Z2...................................... 6. Symmetric Mode Shape Z3....................................... 56 7. Nondimensional Normal Force Versus Nondimensional Time........ 57 8. Nondimensional Central Deflection Versus Nondimensional Time.........................................................., 9. Maximum Response Parameter Versus Initial Velocity Parameter.. 59 10. Critical Velocity Parameter Versus Shape Parameter or Symmetric Buckling........................................... 60 11. Asymmetric Resonant Conditions............................... 61 12. First Zero of Fl, F2, F3 Versus Shape Parameter.............. 62 13. Critical Velocity Parameter Versus Shape Parameter........... 63 iv

NOMENCLATURE A, A1, --- A5 constants of integration A, B, C first order fundamental quantities of the undeformed shell A, B, C first order fundamental quantities of the deformed shell Bn, Cn constants in Chapter II, functions of T elsewhere Dn, En variational part of the solution E modulus of elasticity F arbitrary function of time F1, F2, F3 conditions of stability H rise of the shell L dimensionless Lagrangian M, N bending moment and normal force per unit length of the shell M dimensionless bending moment = 12(1 - v2) a M/Eh3 dimensionless normal force = (1 - v2) N/Eh S, Sn normalizing factors T kinetic energy per unit length of the shell T dimensionless kinetic energy of the shell = (1 - v2) T/ (Eah) To, Tn harmonics of the linear symmetric motion Tn harmonics of the linear asymmetric motion V potential energy per unit length of the shell V dimensionless potential energy = (1 - v2) V/(Eah) Xn complex form of the variational part of the solution a mean radius of the shell bn, Cn pourier coefficients c wave velocity = [E/(1 - v2) p] 2

Cl) -~- cg integrals of riode shapes and their derivatives ds, ds element of arc length in undeformed and deformed state respectively f arbitrary function of Q and t f(G) spatial variation of initial velocity fe(Q), fo(G) even and odd portion of f(G) respectively g(g) spatial variation of initial displacement h shell thickness j, k, 1, m, n integers p eigenvalue r radial coordinate s characteristic exponent in Chapter II and characteristic multiplier in Chapter IV t time V) w tangential and radial displacement of the midsurface respectively tangential and radial displacement anywhere in the shell respectively Vo, wo initial radial velocity and displacement respectively Yk;, fundamental system of solutions A j measure of the detuning = Aj - 2 wn Zo' Zn symmetric, radial mode shapes Zn asymmetric, radial mode shapes To) En symmetric, tangential mode shapes En asymmetric, tangential mode shapes Q period = 2T/(l - ) 2 thickness parameter = h2/12a2 vi

semi-opening angle response parameter 3mn Kronecker delta E initial velocity parameter = voit2/cP2 Ecr critical initial velocity parameter Er radial strain Erg shearing strain EG tangential strain dimensionless radial displacement = w/a 0 polar angle geometric parameter = [ \48 H/h]2 small parameter ALn circular frequency of the nth symmetric mode v Poisson's ratio p mass density per unit volume Cro tangential stress T dimensionless time = ct/a P, On) ~n phase angles dimensionless tangential displacement = v/a Awn circular frequency of the.. nth asymmetric mode vii

INTRODUCTION A problem that has received considerable attention in elastic stability is the snap-through buckling of slightly curved elements. If the applied loading causes the curvature to decrease, the element may "jump" or snap-through to a new equilibrium configuration. The new equilibrium or buckled configuration is characterized by having a curvature opposite to the undeformed shape. The classic example is the static loading of the bottom of an oil can. Upon unloading, the element may snap-through to its original shape, as in the case of a properly designed oil can, or it may remain in the deformed shape. The dynamic analogue of this problem has also received attention. (1>1 Studies include loadings which are suddenly applied and maintained(l) and periodic variations about some mean load.(2) A third type of loading is the impulsive loading. In this case the element is given an initial velocity without any appreciable displacement. For a given distribution the critical velocity is that initial velocity necessary to cause the element to snap-through. All previous analyses of impulsive loading employ an energy method of solution after assuming the spatial form of the displacement functions. The essential difference in these studies lies in the choice of the spatial form. One group of investigators selects a displacement function which approximates the buckled configuration and satisfies the boundary conditions. Humphreys and Bodner(3) did this in their analysis 1Numbers raised and in parentheses designate References. -1

-2of the symmetric buckling of spherical caps and cylindrical panels. This resulted in a single-degree-of-freedom model. Fulton(4) applied the same method to shallow conical shells. Another group prefers to use the mode shapes of an equivalent flat element. Hoff and Bruce used the beam mode shapes to simulate the symmetric and asymmetric motion of their two-degree-of-freedom model of a sinusoidal arch. The five lowest axisymmetric modes (five-degree-of-freedom model) of a circular plate were selected by Budiansky and Roth(5) to study the spherical cap. It seems reasonable to use the vibrational mode shapes for an investigation of the initiation of dynamic buckling. When buckling commences, the deflected shape may not resemble the final buckled configurationo Experimental evidence(6) indicates that the shape during buckling may be asymmetric, whereas the final configuration is symmetrical. In the present study we consider a simply supported, shallow cylindrical shell (see Figure 1). An energy method of solution is used. The assumed displacement function is composed of the linear, vibrational mode shapes of the shallow cylindrical shell. Following a symmetric initial radial velocity, the shell may snap-through or undergo large amplitude oscillations. However in the presence of slight asymmetric irregularities in the initial velocity, the asymmetric modes may grow due to nonlinear coupling between the normal force and the change of curvature. In the present investigation criteria are developed for symmetric dynamic snapthrough buckling and for the initial growth of the asymmetric modes.

CHAPTER I MATHEMATICAL FORMULATION OF THE PROBLEM In this chapter we derive the differential equations and boundary conditions governing the dynamical behavior of a shallow, cylindrical shell. The analysis is based on the following assumptions: a. The loading is normal to the generators and does not vary in the axial direction, i.e. our analysis is restricted to the case of plane strain; b. Elements normal to the midsurface before deformation remain so after deformation and do not change their length. The normal stress acting on surfaces parallel to the midsurface is neglected compared to the other stresses; c. The shell is undergoing undamped, free vibration produced by an initial radial disturbance from the static equilibrium position; d. The material of the shell is elastic, homogeneous and isotropic and remains so during the subsequent motion. A. Strain-Displacement Relations With assumption "a" the strain in the direction of the generators of the shell is zero. Thus the generators will remain straight. The deformation of the shell will be the same in every plane normal to the generators. We let the deformations of the shell be described by the radial displacement w and the tangential displacement 7. In general w and v will be functions of the polar coordinates r and 8 and time t. The radial strain

=4Er and the shearing strain Erg are r = _ - (1.la) 8r E +- a r ~iJ: ii (1.lb) erg = 2 [rg..( r;] (lb) Assumption "b" requires that Er and Erg vanish; therefore Er = - (1.2a) ErG r g] o 0. (1.2b) Integrating (1.2a) leads to w = w(g,t) (1.3) where w is the midsurface radial displacement. Introducing (1.3) into (1.2b) and integrating with respect to r yields Tv =- + rr(Gt), (1.4) in which f is an arbitrary function of G and t. With v as the midsurface tangential displacement we have v = v(G,t) = v(a,G,t) (1.5) where a is the radius of the midsurface of the shell. Setting r = a in (1.4) gives

-5and then (1.4) becomes 6wv = - + 6a (1.6) Equations (1.3) and (1.6) give the displacements at a general point in terms of the midsurface displacements. Taking x1 and x2 as the Cartesian coordinates before deformation of a point lying on a curve whose element of arc length is ds, we have ds2 dx + dx2 = A dr2 + 2Bdrdg + CdQ2 where A,B and C are the first order fundamental quantities of the undeformed state. We let xl and x2 be the Cartesian coordinates of the same point after deformation. The element of length in the deformed state ds is ds2= dx2 + d2 = A dr2 + 2BdrdG + CdQ2 (17) in which A, B and C are the first order fundamental quantities of the deformed state. The Cartesian coordinates in terms of the polar coordinates and displacements are (See Figure 2.) xl = (r-7) sin Q + v cos 0 X2 = (r-w) cos 0 - v sin 0. For finite deformations the tangential strain Eg may be defined as C - C. 2C In the undeformed state A = 1, B = O, C = r2. If we consider an element of arc length in the deformed state such that r is held constant, then dr = 0.

-6After substituting dxl and dx2 into (1.7), we obtain 2 =2 2 + (w + 2 + 2r( ) + r2 d Thus C = - + 2 r + r Therefore co =r -V W +- v 6 + 27 + V 9 r V w[1 2r 9 2w 7 rv The quanitity r a(7 is the tangential strain associated with infinitesimal deformation and may be neglected compared to unity. Then, e@ becomes r ar- Iv-) (1.8)+ 2 EQ + Substituting (1.3) and (1.6) into (1.8) gives = + 12 |v + 3 I + 1_ _ a (19) EQ =a 79- r 2a2 6g a r a a r ~' We introduce into (1.9) the change of variables r = a (1 - - Caz

-7where z is a coordinate measured radially inward from the midsurface. Expanding l/(l-z/a) and retaining terms through second degree in z yields the strain-displacement relation 1a (- ( +w). (1.20) B. Equations of Motion The equations of motion may be derived from Hamilton's variational principle. To use this principle we form the Lagrangian L = T-V, where rI and V are the kinetic and potential energies per unit length of the shell respectively. L is a function of v and w and their derivatives. The principle states that the motion in the interval t1 to t2 renders the integral t2 I = L dt stationary with respect to all variations of the motion vanishing at t and t2. If rotary inertia is neglected, the kinetic energy is 2T P [ 7] [ 1- a a dzdQ -h where p is the mass density per unit volume and h is the thickness of the shell. If the shell is shallow and if the disturbance from equilibrium is produced by an initial radial displacement and/or velocity, the tangential velocity will be small compared to the subsequent radial velocity(3) Therefore the tangential velocity may be neglected in comparison with the

radial velocity2. With this approximation T is 2 T -pah f( 2aw d 2 -_ If we neglect shear deformations the potential energy in the absence of external and body forces is h V 1 I i Eg a (1- Z)dzdg. (1.21) 2 a h 2 With assumption "d" the tangential stress cr is related to Eg by Hooke's iaw, = E EG (1.22) l-v2 E being the modulus of elasticity and v being Poisson's ratio3. If (1.22) and (1.20) are substituted into (1.21) and the integration with respect to z performed, V becomes V 2 {[ -W)- v+ ]+ (2 + 2 Neglecting a7/6t in T is equivalent to neglecting the tangential inertia. This approximation has the limitation that we can not satisfy initial conditions on both the radial and tangential motion. In the subsequent discussion we will satisfy only the initial conditions on the radial motion. 3 The results are equally valid for a circular arch if v is taken to be zero.

-9It is convenient to express T and V in nondimensional forms by defining -w (1.23a) a V (1.23b) a 1 T = E, t (1.23c) - = h2 (1.23d) 12a2 - = T (1-v2) (1.23e) E ah E ah Then, 2 d 0 (1.24a) V - ) + + 2 2 )- ]+2+ ) } dO (1.24b) 2 2 in which the dot denotes differentiation with respect to T and the prime indicates differentiation with respect to 0. The nondimensional Lagrangian L is I + r[(1' )2 + (+ - )22 )2. (1.25) 2 n LW2 The equations of motion are obtained from (1.25) through HamiltonTc principle. They are

-10+ 2 + 2"+ A) - (4I- ) + 2 _ [('- 5) ( + ) + ( ) + = o (1.26a) 2 [(,_'- ) + ~ (, + <')2]' _(,'-_)(,+')) _ ~ (,+5')3 o. (1.26b) The boundary conditions are [ a2 (t +( ) - (t 0-)(-++') - (~+')3] {J } i (1.27a) -1 [~"+ +] {:} = 0 (1.27b) -1 [(47- ) + ~2 (2 + /) 2] {4} = 0. (1.27c) These conditions are satisfied by the vanishing of either the expression in the brackets or the term in the braces at G = + P. The initial conditions are, (O,0) = wo g (G) (1.28a) a ~ (G,0) = f (o) (1.28b) C where c = [ v)2p (1.29) is a wave velocity. The actual initial displacement and velocity distributions are wog (@) and vof (G) respectively. The mathematical description of the problem is completely given by the equations of motion (1.26), the boundary conditions (1.27) and the initial conditions (1.28).

CHAPTER II LINEAR VIBRATIONS In Chapter I the nonlinear equations were derived for plane motion of a shallow, circular, cylindrical shell. For sufficiently small initial disturbances the nonlinear terms in (1.26) and (1.27) may be neglected. The corresponding motion is referred to as linear vibrations and is studied in the present chapter. Love(7N formulated the problem of the incomplete ring. He assumed inextensional motion (no stretching of the midsurface), but retained the tangential inertia. Archer(0) numerically studied the inextensional motion of rings with opening angles from t to 2v. In this and the remaining chapters discussion is limited to shells simply supportedi along the edges Q = + f. The method of solution could be applied, at least in principle, for other boundary conditionso For this case approximate4 analytical solutions are obtained for the radial and tangential mode shapes. A. General Solution of the Linear Equations of Motion The equations of motion for the linear vibrations are obtained by linearizing equations (1.26). They are (r':-.) = + c 2( /V+25" +:) (2. la) 4They are approximate in the sense that only approximate roots are found to the associated frequency equation. -11

(t' )= o (2. lb) The introduction of (2.1a) into (2.lb) yields + f2(0V+2:"' ++) = o. (2.2) Nontrivial solutions of the form = Ae S sin(cppT+0) (2.3) are found to satisfy (2.2) if s5 + 2s3 + (1-p2)s = 0. (2.4) The roots of this equation are s= o, ~ i, ~ i ll Since (2.1) are linear, the general solution is the sum of the linearly independent solutions obtained by substituting the above roots for s in (2.3). We obtain = [A1 sin pJl @ + A2 cos p+l G + A3 sinh @p-l 0 + A4 cosh p-l 1 + A5] sin (apT + 0). (2.5) The introduction of (2.5) into the right-hand side of (2.la) leads to (*'t) = C 2(1_p2) A5 sin (OCpT + 0). (2.6) Equation (2.6) states that (4r'-t) is a function of r only which follows, of course, directly from (2.lb). If we consider the linearized version of (1.20) at z = 0, (t'-5) is the extensional strain of

the midsurfaceo Thus the midsurface strain is spatially constant and is periodic in time. This is a result of neglecting the tangential velocity in (l124a). Whenever A5=0 or p=~l, then (v'-0) = 0, and there is no extension or stretching of the midsurface. This motion is said to be inextensionalo Solving (2.6) for $' and integrating leads to = iCdg + 2(1-p2) A5 G sin (cOpT + M) + F(T) (2.7) where 5 is given by (2o5) and F(T) is an arbitrary function of T The constants of integration Al through A5 and 0, the function F(T) and the eigenvalue p are determined from six boundary conditions and two initial conditions. The linearized boundary conditions are [' +~ ]{~}_ = O (2.8a) ~, {+;s'} _ = o (2.8b) 0t'-:]{t} _: = a * (2.8c) The two initial conditions are ~(g,o) ='E" g(G) (2.9a) 5(,'o) =.f() (2.9b) Once the boundary and initial conditions are specified, ~ and ~ are uniquely determined(9)o

Bo Normal Mode Shapes of a Simply Supported Shell If the shell has the same type of supports at @ = i i, it is convenient to consider the motion to be composed of symmetric and asymmetric motion with respect to Q = 0. For a shell with both edges (Q = ~ B) simply supported the boundary conditions (2.8) reduce to C(PIT) = 0 (2.10a) / P(,T) + (nT) = 0 (2.10b) 4r(a,T) = 0. 10c) If the symmetric motion is considered first, ~ must be even in @ and * odd in G. This requires A1 = A3 = F(T) = 0. Substituting (2o5) and (207) into (2.10) yields the following set of homogeneous, linear, algebraic equations in A2, A4 and A5: cospl T A2 + cosh p-l i A4 + A5 = O (2.11a) -(p+l) cos il A2 + (p-l) cosh fp-i A4 = 0 (2.11lb) sinp+__ sinh A4 + (a2_a2p2+1)p A5O-.. (2.11c) For a nontrivial solution to exist the determinant formed by the coefficients of A2, A4 and A5 must vanish. Expanding the determinant yields the frequency equation tan f71 g + (p+" l, tanh fi - 2p 1t (2ca2p-2+1)p = 0. (2.12) 1)5~/2 ~p-l

Provided that /4T2 > T/2, graphical and numerical studies of (2.12) show that the roots of (2.12) are given with a good degree of accuracy by p+l = (2n+l); n = 1,2,3,... (2.13a) 4r1 (2.13b) p = 0, i 1. (2.13c) Equation (2.12) does not admit repeated roots. Due to the approximations for the roots we restrict (2.13a) and (2.13b) from being identical. It is convenient to define a geometrical parameter X as x= /c (2.14) The physical significance of X becomes apparent from a consideration of the ratio of the rise of the shell H to its thickness. This ratio is - = a (1-cos 2) (.15) h h For shallow shells B is small, and cos B can be approximated by (1-2/2). Noting that a/h = q1}, we may write (2.15) as H _ h a / t or X2 =4- H (2.16)

-16Thus k2 is proportional to the rise to thickness ratio. To require X > rt/2 is equivalent to requiring H/h > 0.36. This certainly does not impose a very strong restriction on the cases of interest. For H/h < 0.36 elementary plate vibration theory could be used without introducing any appreciable error. From (2.13a) p = (2n+l1)22 _1 4p2 For shallow shells B will be less than t/12. Therefore p > 323, and we may neglect unity compared to p. The same argument applies to the root 1/ a. Thus unity will be neglected compared with p, and y2 will be neglected compared with unity in the following discussion. Introducing (2.13a) into (2.11a) and (2.11c) gives A = _ (.-!)n2 A2n (2.17a) 5n (2n+l)a42-(2n+1)t+2 A4n =1 2 A2n (2.17b) [(2n+l)-An-(2n+1)A+2] cosh 2 n 2 where tanh (2n+1)r/2 has been approximated by unity and 2 n= +) (2.18) 16k is the square of the circular frequency of the nth symmetric mode. Introducing (2.13a) and (2.17) into (2.5) gives

-17A0[_(2n+1)iT9 (1)l2 cs (2n+)T)n9...- coI 2(2n+l)xIc 2n [(2n+l)jli2-(2n+l)j+2] cosh ~2n+i) 2 +(2n-_l)n~tn(2n2l~ ~2 sin(PnT +,n) n 1,2.... (2.9a) (2'n'-'('2'n+' l).j++2 (2.92a) Substituting (2.19) into (2.7) and performing the integration yields n = A2n [2nl) sin (2n+l)irg (_")n4p n An(2n-1)t -2 [(2n+l)int2 (2n2n+l)+2](2n+1), cosh (2n+-)l sinh (2n+)lQt (+ 1)n2(1- 2)9 sinh (2ni)- 2- sin(4nT+On); n=1,2,... (2.19b) (2n+l)1 4n-(2n+l)t+2 The introduction of (2.13T) into (2.11a) and (2.11c) leads to A A |sin X cos A50 =A2oj0flh -cos x] sin X 40 - 20 sinh X If we replace A5 and A4 in (2.5) by the above expressions, the solution corresponding to root (2.13b) is o = A20 [cos XQ sin X cosh X + sin X cos X sin(T+O0). (2.20a) -~- sinh X - tanh' Substituting (2.20a) into (2.7) gives -sih ~ sinh 1 ( to = A20 S [in. X-sin X inh XQ ]sin(T+0o). (2.20b) The roots (2.13c) lead tsin solution. The roots (2.13c) lead to the trivial solution.

Due to the linearity the radial displacement for symmetric motion may be represented by superimposing (2.19a) for all n and (2.20a). Likewise the tangential displacement is given by superimposing (2.l9b) and (2.20b). We have 00 00 = Z Bn Sn = Z Bn Zn (0) Tn(T) (2.21a) n=o n=o ~~~00 ~00' = Z Bn Vn = Z Bn Yn(Q) Tn(T), (2.21b) n=o n=o where To(T) = sin (T+0o) (2. 2a) Tn(T) = sin (tn T + On) n = 1,2,.... (2.22b) Zn(Q) is the orthonormal, symmetric, radial mode shape obtained by choosing A2n and A2,) in (2.19a) and (2.20a) so that f Zm Zn dQ = mn, (2.23) where imn is the Kronecker delta. The tangential mode shape is Tn and in general is neither orthogonal nor normal. The symmetric mode shapes are z --- [cos XQ sin X cosh Q + sinX os (2.24a) o 1/2L sinhX f tanh X J =- _ [ Q sin (2.24b) =, sin s inh (2.24b) 0 X D sinh X 5See Appendix A for proof of orthogonality of the radial mode shapes.

-19Zn = S cos (2n+l)eg + An cosh (2n+1)' + An (2.24c) 1 —~ 2P cosh (2n+l)i 2c 2 Sn 28 (2n+l)> 2PAn (2n+l)An (2n+t n(2nrl)n si g (2n+l)n cosh (2n+l)iSinh 2P + An(1-tn) ]; n = 1,2,3,..., (2.24d) in which An = (-1)2 (2.25a) (2n+1)7 _ I(2n+ 1 )+2 S 1 F( osX + sin X sin X + s si n X cos ]2 (2.25b) tanh X X sin X tanh X l+ (1)n4A ~2 6An2 -~ S 1 + ()n4An + 2An 2 ]n (2.25c) (2n+l>) ( 2n+) In these expressions [cosh (2n+1)A/2]-2 has been neglected compared to unity. Figures 3 through 6 show the symmetric radial mode shapes for various values of X For sin X =cOS X, tanh X A50 vanishes and the modes Zo and To are inextensional. The values of X satisfying this condition are given approximately by = (4m+l); 1,2,... 4; m = 1,...

Since A50 only vanishes for the above values of X and A5n tends to zero for large n only, the symmetrical motion is neither purely inextensional nor purely extensional. For the asymmetric case A2 = A4 = A5 = 0. Since A5 vanishes, we conclude that the asymmetric motion is inextensional. This result is to be expected from the form of the asymmetric motion. Applying the boundary conditions (2.10a) and (2.10b) yields sin jl A1 + sinh pT-1 i A3 = 0 - (p+l) siny1p+l Al + (p-l) sinh\p-1 P A3 = 0. For a nontrivial solution to exist the determinant of the coefficients of Al and A3 must vanish. The roots of the resulting frequency equation are (2.13c) and J1 =p; n= 1,2,.... (2.26) Proceeding as in the symmetric case, we have n = A sln sin sin[(22) T + n]; n = 1,2,... (2.27a) x 2 Substituting (2.27a) into (2.7) and performing the integration leads to =n Al cosn$i] sin [(n2,+2_-2) n7 1 c s Xs in + n + F (2.27b) *n nt P 2 I

Introducing (2.27b) into (2.8c) yields A. g + ~n ~22 PI(T) A l n (_l) sin ( T 7n Therefore ALnp [ nQ [(n2 It tn = [(-1)n cos ] sin [(n 22) ~ n]. (2.27c) The radial and tangential asymmetric displacements are given by respectively summing (2.27a) and (2.27c) over all n. They are oo 5 = Z C n = Cn Zn (Q) C n (T) (2.28a) n=l n=l 00 -= Z Cnn Zn = Cn -n (Q) Tn(T) (2.28b) n=l n=l where Tn(T) = sin [(.n22. T+~n1 (2.29) X2'7] (2.29) Zn () is the orthonormal, asymmetric, radial mode shape obtained by choosing Aln so that Zn satisfies (2.23) when Zn is replaced by Zn ~ -n(@) is the corresponding tangential mode shape. They are 1si nn (2.30a) n() W= — (/-l)n - cos n]; n = 1,2,3,.(2.30b) n nit

-22C. Response to an Arbitrary Initial Radial Velocity If the shell is given an arbitrary, initial radial velocity, then g(G) in (2.9a) is identically zero. The initial velocity distribution f(Q) may be written as f(@) = fe(Q) + fo(Q) (2..L3a) where fe() -= [f(o) + f(-Q)] (2.13b) 2 is an even function of @, and f,0() = -[f(Q) - f(-Q)] (2.31c) is an odd function of G. If f(Q) is piecewise continuous, then f(@) may be represented as an infinite series of the radial mode shapes. The representations of fe(g) and fo(Q) are 00 fe(g) = Z bn Zn (2.32a) n=o 00 fo(g) = Z Cn Zn (2.32b) n=l where bn = f fe(g) Zn dQ (2.33a) Cn = f fo(g) Zn dO. (2.33b)

With (2.32) and (2.33) the initial velocity representation is i(Qo) = [ bn Zn + Cn Z. (2.n4a) n=o n=l In addition we have:(Qj) = o. (2.34b) The radial displacement is given by (2.21a). Equation (2.34b) is satisfied by setting On = =n = 0 in (2.21a). Taking the partial derivative with respect to T and evaluating the results at T = 0 gives 00 00 (n2.g2_p21 O [(Q,0) = Bo Zo n Bn Z+ Z n Bn Zn + 2-) C (235) n=l n=l X Equating (2.35) to (2.34), multiplying both sides by Zo dG and integrating from G = -I to 5 yields Bo = bo0V (2.36a) c In a similar manner bnvo Bn = (2.36b) Cn = (n ); n=l, 2, 3,. (2. Substituting (2.36) into (2.21) and (2.28) and adding gives the response to an arbitrary, initial, radial velocity as

-24 - c [ n=l An + 2 Zn sin -n22-~f2) (2,37a) n=l n 2- Xn =- b= c o To sin T + n nsin n T cl.00 n=l In oo X2cn 2_2 2 + 22 2 n sin 2. (237b) n=l n 72 The tangential normal force and the bending moment per unit length acting on the radial cut, Q = constant, are defined as 2 h N Cr " az (2.38a) h 2 M = Z ag z. (2.38b) Even though plane sections remain plane, the stress does not vary linearly due to the initial curvature of the shell. For thin shells (a/h > 20) the distribution is nearly linear.(l0) ThLs, the assumption that the stress distribution is linear is reasonable for thin shells and consistent with Love's first approximation. If (1.22), in which ~Q is given by (1.20) after deleting the second and last term, is substituted into (2.38), we obtain

N(1-v2) (7') (239a) Eh l E2Ma(-v2) = (/ +) (2.j39b) where N and M are convenient dimensionless stress resultants. Substituting (2.37) into (2.39) gives V 0bos sin X bnAnSnln 1 N - c[l/2 tanhX - cos )sin T +n T (2.40a) vo[bosX2 XQ X2 sin X XQ sin M = - - + cosh + cos X- s in T Sc 1I2 C0 2 i s2 SinhX ~ tanh + + z b ( 2n-+)22 2 o (2n+l)it + An(2n+l)osh n=l 4 2P 4p2cosh(2n+l)~T 2P cT + Z cns (n)it- ) - Ansinn 5/2 sin sin T + (2.40b) n=l,8/- s i where unity has been neglected compared to (2n+l) 2 and 1 4p a As can be seen from (2.40a), the tangential normal force is spatially constant. Hoff and Bruce(l) assumed this in their analysis. Another point of interest is a beating that may occur between the first term and the mth term of the infinite series for N. It occurs whenever X is such that the first and mth term have nearly the same frequency, i.e. whenever (2m+1)2i2 4X2

-2La - If this reinforcement is significant, the normal force developed for a given initial velocity will be substantially larger for some values of X. Shells with such values of X2 could possibly be more susceptible to dynamic instability. If sin X cos X =, the tanh X first term in N will vanish as noted in Section B. D. Example of a Specific Initial Radial Velocity As an example we consider the initial radial velocity 5(O,0) v~ cos It (2.41) c 2P This initial velocity satisfies the boundary conditions. Such a distribution could reasonably be expected for blast loadings in which the duration of the loading is less than the fundamental period of oscillation. For this distribution f(@) is an even function f(G) = fe(G) = cos A. (2.42) 2P Since fo(g) = 0, it follows from (2.23b) that cn = O. Substitution of (2.42) into (2.33a) yields upon integration 162S/2[cosX + 2+4si)tan h X ] (2.43a) +4AnSngl/2 (2n+1)2 b' n =1,2.... (2.43b) nrr (2n+1)2+1

-7The introduction of (2.24) and (2.43) into (2.37) leads to Vo 16X2S2 cos X sin X G sin X G sin I = 2-4- + 1 cos -- - cosh + {16- (k2+422)tanh 2[ kg sin (n2+4k2)tanh c 242 +4X2)tanh X J sinh X tanh X 16X2 An Sn cos X sin T +. z nn+ Cos (2n+l)~Q 1 n= (2n+)2+1 2 An (2n+l)ir h (2n+ + An sin n T (2.44a) cos ~r(2n+)i co 2 ] n 2 Vo 16xpS2 cos X sin x ] * sin X ~l...... sin X si ksinh - Sin T c \ S [T'-4k2+ (r2+4k2)tanh sinh i]sin 4f= ~2 05k sik 2k sinh + 16X2 1 An S 2p sin (2n+l)rG _ 2_ An 5t3 n=l (2n+1)2+1 L (2n+) 2 (2n+l)Acosh(2n+l)/2 sinh (2n+1)tg + An(l-n)Q] sin 1n T (2.44b) Making the same substitutions into (2.40) gives -_Vo 16X2S2 [cos A + sin + ] sin -cos c it LT2- 2 ( kt2+4X2)th k tanh J tanh X 22 4 + An Sn (2n+l) sin n (2.45a) n=l x2 (2n+l)2+l

- = vo 16X2S2 cos + sin X in i X2 sin X o G c T S it 2-4x2 (.2+4X2) tanhX 2 sinh X i sin X + Cos X sin 1 i 16X2 An(S2 2n+1)222 (2n+l ( b) 4p2cosh(2n+l)t/2 2f n n Figure 7 is a plot of nondimensional normal force as a function of the nondimensional time for three values of.2. For -2 = 68.5 and m -- 2, Em = 0.9. From the discussion in Section C a beating should occur between the first and second term of the infinite series for N. This is clearly depicted on Figure 7. The maximum value of N is approximately seven times greater than the maximum values for other values of. The value X2 = 100 falls approximately in the middle of the range of values of X2 which are of interest. k2 = 178 is near the upper end of this range, and it is also a value which causes the first term in N to vanish. Figure 8 is a plot of nondimensional central deflection versus nondimensional time. There is no evidence of beating of the central deflection. For this example the first two terms in (2.44a) essentially determine the central deflection. This is clearly indicated in Figure 8 since the central deflection appears to be composed of harmonics whose circular frequencies are unity and 91T2/4X2

CHAPTER III SYMMETRIC DYNAMIC BUCKLING In this chapter we consider the large amplitude oscillations of the shell described in Chapter II. This motion may lead to what is referred to as "dynamic snap-through buckling." The objective is to determine configurations for which symmetric "dynamic snap-through buckling" is possible. In addition the initial velocity necessary to cause shells of these configurations to snap through is determined. An exact analysis would require the study of the stability of the nonlinear partial differential equations (1.26) with the appropriate boundary and initial conditions (1.27) and (1.28). However, there are many formidable mathematical obstacles. Here an energy method is used to study approximately this large amplitude motion. Such methods have been used in dynamic stability problems by a number of authors. Pertinent to the present study are the investigations of Goodier and McIvor,(ll) McIvor,(12) and Humphreys and Bodner.(3) The use of this method by Russian authors is described by Bolotin.(2) The particular method employed is an extension of that used by Goodier and McIvor. Briefly, the method consists of assuming that the motion can be represented by an infinite series of products of arbitrary functions of time and the linear mode shapes of Chapter II. This assumed form will satisfy the boundary condition (1.27). By specifying appropriate initial conditions on the arbitrary functions the initial conditions (1.28) may be satisfied. The arbitrary functions are considered as generalized coordinates. The series is truncated after two terms and the associated nonlinear Lagrangian equations are integrated numerically.

A. Derivation of Governing Equations The linear mode shapes obtained in Chapter II are geometrically admissible. Hence we assume for ~ and 4 the forms oo oo = Z Bn(T) Zn(@) + Z Cn(T) Zn(@) (3.la) n=o n —i 00 00 t = 2 Bn(T) n(9) + 2 Cn(T) En(@). (3.lb) n=o n=l Bn and Cn are arbitrary functions of T as opposed to constant coefficients in Chapter II. The first series in (3.1) represents the symmetric motion whereas the second corresponds to the asymmetric motion. With this choice of ~ and r linearization of the nonlinear differential equations gives the correct linear equations for harmonic motion. This seems to be a realistic extension of the linear problem as noted by Goodier and McIvor. The assumed solution, like the linear solution, yields a normal force that is spatially constant. As previously noted Hoff and Bruce also assumed the normal force to be constant. If the shell is given an initial radial velocity, then (3.la) must satisfy (2.34). Therefore Bn(O) = 0 (3.2a) Bn(O) = bn; n = 0,1,2,3... (3.2b) Cn(0) = 0 (3.2c) Cn() cn; n = 1,2,3,.... (3.2d) In the sequel we may neglect 4 compared to j' and ~ compared to t". For the assumed form of ~ and 4 this is equivalent to neglecting unity compared to n2i2/,2, (2n+1)2X2//2 and i/o. This is consistent

with earlier approximations. With these approximations (1.25) reduces to = 1 {2 [fit) 1(f)2]2 _ 2(,)2}. (3.3) _ For the symmetric motion Cn is set equal to zero. In Chapter II the central deflection was found to be closely approximated by two terms which are equivalent to the first two terms of (3.1a). It is reasonable to expect that they would predominate in the somewhat larger amplitude of the nonlinear motion. The mode shapes are so complex that an extensive numerical integration technique would be required if more than two terms were included. Therefore, in the representations of 5 and 4 only the first two terms will be retained. This represents the next refinement above the single-degree-of-freedom model. It is convenient to introduce the change of variables 5/2 - Bn = P Bn; n = 0,1,2,.... (3.4) Substituting the first two terms of (3.1) into (3.3) and using (3.4) yields L =2 {B + B2 1 - + cl + c B1 + c3BOBl + c4 -4 *-3- -2-2 — 3 -4 - C5Bo - C6BOB1 - C7BOB1 - C8BOB1 - c9B1} (3.) in which6 cl =-P 5/ (To-Zo)(zo) dG (3.6a) 6See Appendix B for the evaluation of the integrals.

2 5/2 2 C2 =- f [(iP-zl)(z') + 2(' -zo)Z'Zo] dG (3.6b) -030 0 lo 5/2 2 c3 -5/ / j [(o-Zo)(Zl) + 2(Ti-Z1) Ziz;] d@ (3.6c) c4 _5/2 j ({i-z )(zi) d) (3.6d) p5 p 4 C5 =- f (z~) dO (3.6e) -0 3 C6 - (Z) zt dO (3.6r) -13 C8 - fJ (z)(Z)3 d (3.6h) The differential equations governing Bo and B1 are obtained from Lagrange's equation of motion d (l| | = O (3.7) dT I I qi

-353by treating Bo and B1 as the generalized coordinates qi' Substituting (3.5) into (3.7) gives Bo- B1o-2 3 3 -2Bo + Bo 2- c.lB0 - c2BoBj - 2 c3B1 + 2c5Bo + 3 c6BOBl + c7Bo- +1 2 C8Bl o (3.8a) ~ 3o B1 + jBl 2 c4B1 - c3BOB1 12 - 2o c2B + c6Bo + c7oB1 c8BB1 + +2cgB - = (3.8b) where the dot now indicates the ordinary derivative with respect to T For the initial velocity (2.41), Equations (3.2) and (2.43) yield the initial conditions Bo(0) = B1(0) - 0 (3.9a) 16eX2SS cos X sin X (3.9b) j3 [2-4X2 (gt2+4X2)tanh ]J Bl(0) = 36 eAbS (3.9c) where the initial velocity parameter e is defined as = v2 (3.10) The symmetric motion is governed by the differential equations (3.8) and initial conditions (3.9). B. Method of Solution and Results Qualitatively, we say that the shell has snapped through if upon a small change of the initial velocity the deformation "Jumps" from

moderate to large. This motion is unstable. We may apply the classical method of Poincare(l3) to study the stability of the motion. The method requires finding the singular or equilibrium points associated with (3.8). Dynamic snap-through buckling may occur when unstable equilibrium points exist. Finding these singular points requires the solution of two simultaneous, nonlinear, algebraic equations. Since the coefficients of (3.8) are algebraically involved functions of X, this would have to be done numerically. Rather than solve these equations numerically, Equations (3.8) with initial conditions (3.9) were integrated numerically. This has the advantage that not only does it exhibit the instability but also yields the response in the case of no unstable equilibrium points. A comparison can be made of the response of a shell for which instability is possible with the response of a shell which is geometrically very similar yet does not exhibit snap-through buckling. It is convenient to introduce a response parameter(5) (w)ave (H)ave which is the ratio of the average radial deflection (w)ave to the average rise of the shell (H)ave. From (3.la) this is 5 = 3B sin cos + - [(3-2) A1-2 ] (3.11) tanh X Figure 9 is a plot of 6max, the first maximum value of 8, as a function of C for representative values of X. Considering X = 2.5, we see that for small values of E, bmax is nearly proportional to E indicating an almost linear response. As E increases this relation becomes nonlinear. Near E = 1.16 a small change of the initial velocity

produces a jump in the average displacement from less than the average rise to greater than the average rise; i.e. it jumps from a moderate to a large value, We conclude that the shell has snapped through. This value of E is taken as the critical initial velocity Ecr This curve resembles the load-deflection curve of a shallow shell with a dead static loading. (14) Similar curves were obtained by Budiansky and Roth.(5) The same phenomenon is observed for X = 1.80, but no jumps occur for X < 1.75. At the other extreme a jump occurs for X = 2.85, but none for X > 2.90. The conclusion is that symmetric dynamic snapthrough buckling can occur if 1.80 < X < 2085. This result is valid for arbitrary, initial radial velocities, provided the motion can be represented by the first two terms of (3.1), since the singular or equilibrium points are independent of the initial velocity. In contrast to these results are those obtained by Humphreys and Bodner. They found symmetric dynamic snap-through buckling could occur for X > 2.78 ~ The stability criterion is clearly sensitive to the assumed form of the solution since the formulations were similar. For initial velocities slightly less than critical the shell is undergoing essentially linear vibrations. Therefore, using the linear mode shapes would seem to be indicated. For the statically analogous problem Hoff and Bruce found that symmetric snap-through buckling due to a sinusoidal static loading was possible for 2 < X < 2.83. This range is the same order of magnitude obtained for the dynamic problem. Figure 10 shows the relationship between Ocr and X2 for symmetric dynamic snap-through buckling for 1.80 < X < 2.85. These results are about one percent lower than those obtained by applying

Poincare' s method when only the first term of (3.1) is used. For this range of X we conclude that the single-degree-of-freedom model is sufficient. In the region of common validity Humphreys and Bodner obtained critical velocities about twice those presented in Figure 10. Again this can probably be attributed to the assumed form of displacements and to a lesser degree to the fact that they assumed a uniform initial velocity. For values of X not in the set 1.80 < X < 2.85, the maximum response increases monotonically with E, and no jumps appear. From an experimental viewpoint, it may be difficult to determine if snapthrough buckling occurs for values of X near either end of the set. For instance if we consider the curves for X = 1.80 and X = 1.75, there is not a great deal of difference between them. Without sensitive instrumentation it would be difficult to decide whether a shell with X = 1.75 snaps through or just experiences large scale deformations.

CHAPTER IV EXCITATION OF ASYMMETRIC MODES Empirical evidence(6) indicates that the final configuration of symmetrically loaded spherical caps is symmetric. There is also evidence that the initial deformation may be partially asymmetrical and by means of a complex mechanism end in a symmetrical configuration. Theoretical studies (e.g., see References 15 and 16) of symmetrical static loading of shallow shells indicate that asymmetric deformations may initiate the instability. There is no reason to believe that this same phenomenon does not occur in the dynamic case. This chapter is concerned with the stability of the shell when asymmetric modes are excited by a small asymmetric perturbation of the symmetric initial velocity. A criterion is developed for the initial growth of the asymmetric modes. No effort is made to determine their behavior with increasing deformation. This would necessitate including the fourth degree terms in the Lagrangian which are neglected in the present analysis. It is likely that the asymmetric modes grow to a point where they become unstable so that the motion is again symmetrical. A. Derivation of Governing Equations The initial growth of the asymmetric motion is produced by the nonlinear coupling between the normal force and the change in curvature. If the phase of the normal force relative to the curvature change is such that the normal force tends to accentuate the curvature change, -37-.5

-38then the asymmetric motion grows and the shell is said to be dynamically uhstable. To include this nonlinear coupling in the initial motion it is necessary to retain only the third degree terms in (3.3). If we proceed as in Chapter III, the Lagrangian L with terms through third degree is computed for the displacement representations (3.1). The Lagrangian is then substituted in (3.7). Taking qi = Bo and Bn in (3.7) yields 3 2 sin X cos sin X Bo + Bo Sg tn cos B - 2S cos emBoBm 2 tanh X - g m AmSmmBoBm - ( tanh - cos kl Z dkmBkBm k~m 0000 00 - Z Z emAkSkkB kBm - tanh cos fmB2 1 ( s in X _ cos m2Cm = (4. la) 2 tanh X m=l +2 n ngBo S - sin Xcos eB-2AnSn- n emBoBm'~ + IBnBn - cs [ Bn +[nn 2 tanhX m=l sin X sin X - S tanh - cos ml dmnBoBm - S tanh - cos ) fnBBn mAn 0c2AS 20 00 00 2o - en Z AmSmpmBoBm An2Sn n E dkmBkBm - Z Z dmnAkSk0kBkBm m=l k=l m=l k=l m=l kkm mAn 0o oo 00 -fn Z AmSmBmBn - AnSn2 m2C2 0; m=l m=l m=l n = 1, 2, 3,..., (4.lb)

-39 - in which SkSm(2k+l) (2m+l) [ (-1)kAm(2m+l) (-1)m Ak(2k+l) dkm = dmk = ~ 5/2 i(2k+1)2 + (2m+1)2 (2m+1)2 + (2k+1)2 AkAm +(2m+l)+ (2k+l) (4.2a) S Sm(2m+l)t (-l)mcos X (-l)m sin X Am cos em =1- + tanh -(vkm s in X- c kos) pL52 kLm-l kim+l tanh X + m+l + I-lm 1 Tm sin X h sin X (4. 2b) SmX ktm r (-14 Am 2 Am -m =/ [1 )m (4.2c) f = I5/2 + (2m+l) (2n+l(4.2c) __ ___ S2mkMm -1in2 A (4.2d) g3= —1+ sin h - osin S 3sin s 9 52 [1 3 snX tanh X _cos X sinh2 (4.2d) where tanh(2n+l)ct/2 has been replaced by unity and [cosh(2n+l)rt/2]-2 neglected compared to unity. Setting qi = Cn gives 2 -n2ic2 sin X c + AmSmBm n = l, 2,..., (43) where 2 n44 (4. 4) The dot again indicates the ordinary derivative with respect to T The initial displacement is taken to be zero. The initial velocity distribution can be expanded in a series of the orthonormal radial modes. This representation is a (,o) vo bnZn(O) (4.5a) aT c n=o

In addition we consider a small asymmetric perturbation of this velocity which can be represented by 3T c nl where the an's are small compared to unity. The initial conditions on Bn and Cn become Bn(0) =0 (4.6a): Vo — bn; n = 0,1,2,3,.... (4.6b) Cn(O) 0= (4.6c) Cn - a~nn; n =1,2,13,... (4.6d) The stability of the shell is governed by the differential equations (4.1) and (4.3) with initial conditions (4.6). B. Initial Growth of the Asmetric Modes Initially the Cn's are small since the initial disturbance is small. Therefore, their squares may be neglected in (4.1). Further if the initial velocity (4.6b) is not too large, the unperturbed motion may be represented by the linear theory of Chapter II. Figure 9 shows that this is a good approximation if e is less than 0.4. Thus, Bo and Bn are taken as Bo = bo sin T (4.7a) Vo bn) Bn =- - sin nT; n - 1,2,3,.... ( ) c.n

-41Substitution of (4.7) into (4.3) leads to Hill's equation for the asymmetric motion,. 2O O00 Cn + Cn = en2[bo sin T + Z bm sin kmT] Cn; n = 1,2,3,... (4.8) m=l where b S sin X - bo = 1/2 tanh (4 9a) =bmAmSm-m; m = 1,2,3.... (4.9b) The term in the brackets of (4.8) is proportional to the normal force obtained in Chapter II. The growth of the asymmetric modes is governed by the stability of (4.8). The stability criterion may be obtained by applying the method of Struble. (17 18)7 It is convenient to write (4.8) as Cn = Xn (4.10a) 00 Xn + nCn_ = en2[bo sin T + Z b-m sin PmT] Cn; n = 1,2,3,... (4.10b) m=l Struble considers the solution to be composed of two parts. The first part is taken in the form of a variation of parameters solution, and the second part consists of a perturbation solution. Thus, a solution is assumed to have the form 0 Cn = Dn(T) sin WcnT + En(T) cos conT + Z Ej fnj(r) (4.11a) j=l 70ther methods(l9,20) usually require a priori knowledge of the period of the bracketed term in (4.8). Here the period is in general unknown.

00 Xn = %nDn(T) cos %nT -w nEn(T) sin CnT + Z Ej fnj(T). (4.llb) j=1 Substituting (4.11) into (4.10), retaining terms of zeroeth and first order in e and introducing suitable trigometric identities yields Dn sin anT + En cos )nT = 0 (4.12a) 2 En2Dn onDn COS )nsT - LnEn sin ounT + e(fnl+Cinfnl) = 2 D{o os (1-tnn)nT - CoS(l+W)7T] + 7. b-lm[COS(Gm-n)T - coS( pm+Th)T]} m=l 2 00 En En + 2 O[sin(l+)T] + sin(ml-c+)T] + Z bm[sin(Lm+un)T m=l + sin(tm-1n)T]}. (4.12b) Struble's method consists of associating all the members of the righthand side of (4.12b) except those which lead to very large or infinite magnitudes with the perturbation solution. The remaining terms are associated with the variation of parameters part of the solution. If cu1 is not near the circular frequency of any harmonic on the righthand side of (4.12b), Dn and En are constants and n2Dn- cos(l-wn)T Cos(l+cwn)T C - [ cos (lm-con)T fnl -2 o _ b(1. <n) -( - a-(l+~n) + m, (mn) cos2(Im+Cin)T n2 |b [sin(l+con)T sin(l-un)T n ~(m+Ln)21} 2 \ [ (l+on) +c-(1 —cn)2 + z bm [sin(km+cn)T +sin( mk-Mn)T] m=l (2-((km+m) 2nm ) 2 Clearly, terms with large amplitudes occur when the system is at "resonance", or when

and/or Cn = + (Pm + Ln) n Of all the possible sign combinations, the only physically possible conditions (an > 0) are C 2n =1 (4.13a) and/or -n 2 (4.13b) Significant deviation from the linear solution is possible if the circular frequency of the asymmetric mode is close to one half of the circular frequency of a harmonic in the normal force. This is a condition for which the normal force may accentuate the curvature change associated with the asymmetric motion. An analogous phenomenon is reported in Reference 11. Other examples of the 1:2 frequency ratio are cited by McLachlan(21) and Den Hartog. (22) With (2.18) and (4.4) the resonant conditions (4.13) may be written as n2c2 - 1 (4.14a) x~' 2 and/or n2 (2m+l)2 (4.14b) 8 There are no integral values of n and m which satisfy (4.14b). Thus for a given n, m is taken to be the integer "closest to" 12n - 1/2. Whether or not there are integral values of n satisfying (4.14a) depends on the value of X. In either case we take n to be the integer "closest to" X/ f2c. For this integer there is a second integer j which nearly satisfies (4.14b). Since (4.14a) cannot exist

-4) - alone, there are two cases of interest, the distinct case and the nondistinct case. The first case refers to condition (4.14b) and the latter occurs when conditions (4.14) exist simultaneously. Figure 11 is a plot of the resonant conditions. The variational part of the solution containing the resonant terms of the distinct case is Dn sin WnT + En cos WnT = O (4.15a) cnDn cos XnT - TnEn sin 2En Dnbj cOS(olj - Wn) + 2 en2Enbj sin(pij - Wn)T (4.15b) It is convenient to introduce a measure of the detuning; i.e. the deviation of the actual circular frequencies from the resonant conditions, as Aj = pj - 2n. (4.16) Introducing (4.16) into (4.15) and solving for the first derivatives yields Dn = [Dn cos(Aj + OnT) + En sin(Aj + nn)T] cos LnT (4.17a) 2wn En =2j E - 2n [ Dn cos(ej + +n)T + En sin(.j + Wn)T] sin TnT. (4.17b) Applying the averaging method of Kryloff-Bogoliubov, (23) we average the right-hand sides of (4.17) over the period 21t/wn of the fundamental harmonic. The results are Dn - [D cos + En sin (4.18a) Dn = E [Dn cos AjT + En sin AjT] (4.18a)

-45En -= n2 [-Dn sin AjT + En COS AjT]. (4.18b) Equations (4.18) may be written in complex form as =n =4 Xn e T (4.19) where Xn = Dn + iEn and Xn is the complex conjugate of Xn. The solution of (4.19) is Xn = Xon exp(p + 2 ) T in which Xon is a complex constant and p is a root of p2 E- -+ = 0. (4.20) 16 an 4 For the motion to be bounded and hence stable the real part of p must be negative. The stability criterion is F1 > 0 (4.21) where 2 2 j E2 n4b,2j F1 -= J. (4.22) A point of interest is that the same criterion within 0(e) would have been obtained if Struble's method had been applied to the Mathieu equation 2 2Cn + (2n - En bj sin pijT) Cn = 0 o Proceeding as in the distinct case and applying the averaging method to the variational part of the solution, we have for the nondistinct case En2DEn 2 Dn Dn (bo cos AT + bj cos AjT) + n (bo sin AT + bj sin AjT) 4cn 4kcun (4.23a)

En= n2D n (bo Asi n A + bj sin AjT) - en2En( os b 4aun 4acco (4.23b) where A = (1-2n). (4.24) To study the stability of (4.23), it is advantageous to introduce the change of variables yl(T) = Dn cos 2- T + En sin T (4.25a) j AI Y2(T) = -Dn sin + En cos. (4.25b) Equations (4.23), written in matrix form, are F n2b: Yl f'n2 Ys1-j)T sin(l-CIj) T yl 4o2nn 22 +en b An2 4 Y2 - 2 4 Y2 sin(l-[jj)T -cos(l-[j)T y2 (4.26) This is a set of linear differential equations with periodic coefficients of period = ( (4.27) Floquet theory(24) asserts that the matrix equation y = p(T)y, (4.28) where y is a column matrix and p(T) is a periodic square matrix of period Q, has a solution with the property y(T + Q) = sy(T) s is a characteristic multiplier. y may be written as Y = Y()d

where y is a square matrix of the fundamental system of solutions Yk2 and d is a constant column matrix. y(T) satisfy the initial conditions y(O) = I, (4.29) I being the identity matrix. Therefore (y(O) - sI) d = 0 The condition for a nontrivial solution is Y(Q2) - sIt = O. (4.30) The solution is asymptotically stable if the absolute values of the characteristic multipliers are less than unity. The difficulty with this analysis is the characteristic equation (4.30) depends upon knowing the solution which in general is difficult to find. Since the fundamental system of solutions of (4.26) has not been found, Malkin's procedure(25) is used to obtain an approximate characteristic equation. In essence the procedure obtains an asymptotic expansion for the fundamental system of solutions. The matrix p is assumed to depend analytically upon some small parameter p. so that it may be expanded in a power series in. P = Po + APl + k2P2 +'(431) (4.31) where Po, Pl P2,'.' are continuous periodic matrices of period. The fundamental system of solution is assumed to be represented by a power series expansion in p. Y = Yo + [Y1 + 2 Y2 + (4.32) Equation (4.29) is satisfied if Yo(O) = I Y1(o) 2(0) ='... = 0o. (4.33)

Substituting (4.31) and (4.32) into (4.28) and equating coefficients of like powers of p yields the following recursive system: Yo = PoYo Yl = Poyl + Plyo (4 34) Yi = Z PmYi-m m=o If the solution to the first equation of (4.34) can be found, e.g. if po is a constant matrix, then solutions to the successive equations can be effected. Equation (4.32) is substituted into (4.30) which yields the characteristic equation. If kt is taken as en bo then a comparison of (4.28) and (4.31) with (4.26) yields En 2b.j A. 4=n 2 Po _j _ ~n2 2 cos(l-4j)7 sin(l-:j)T Pl = sin( 1-j) -cos (1-kj) p2 = P3 =... =. The integration of (4.34) is straight forward since po is a constant matrix. It is necessary to retain terms of first order in p because only terms of this order were kept in Struble's method. After considerable

-49algebraic manipulation Yk (2) are found to be Yll(2) = cos [n2b en2bI (1-1j ] | 2 s (4.35a) (l-Cj) I48c 4any (l-(j) 4Y (l-j) Y12) 4= C)n7+ (1 - ~J) 2 _47 sin (4.35b) Y() = _Y2 1 +E n2b+oj E (-1j s j)n sin 2 - (4-.35c) 22 (l-+)'Y 1nj2 (l-#J) -4J (l-)(4.35) 2 -2 1 sin 1/2 where Y = t4-: nj (4.36) For stability we require that y be real, which is equivalent to the stability criterion of the distinct case. It is convenient to introduce the change of variables s = 1 + _. (4.37) 1 - This maps the interior of the unit circle onto the left-half of the complex plane. The condition Isl < 1 is equivalent to Re r < 0 and Hurwitz's criterion may be applied to the resulting polynomial. Introducing (4.35) and (4.37) into (4.30) yields F2r2 + 2F4r + F3 = O (4.38) where F2n4b2o (1-1lj -j)2._,, 2ity E n b 2 2 c(4.39a) F2 = 2 + 2 cos ___ -sin2 (49a (1_-j) 4~ [(l_-j)2 _ 472]2 (_l(j) F3 = 2 - 2 cos 2cn bo sin2 2Tyr (4.39b) F4 = e2n-os A(1-1), _ 2 n F4 0 k, - 4yj sin2 27(4.39c) F4 21 = ~ )2 -4212 s n2 F4~~~~~~ (}_-~j

Since F4 is non-negative, q will have a negative real part if and only if F2 and F3 are positive. The stability criterion is F1 > 0 (4.21) F2 > 0 (4.40a) F3 > 0. (4.40b) Instability is implied if any one of the inequalities are not satisfied. For a given initial velocity distribution and X, the critical value of E is the smallest of the zeros of Fl, F2 and F3. For the distinct case we consider n = 1 and from Figure 11 find the appropriate integer j. After the quantities in F1 are determined the zero is obtained numerically. The procedure is repeated for n = 2,3,.... We would supposedly have to consider all n, but bj tends to zero very rapidly so only the first few n's are necessary. For the nondistinct case we enter Figure 11 with X and determine the integer n and the corresponding j. The first zeros of F2 and F3 are determined numerically. For the initial velocity (2.41) Figure 12 shows a plot of the first zeros of F1, F2 and F3 versus X2 for various values of n The curves are similar to those for the asymmetric buckling of clamped shallow spherical shells subjected to a uniform static pressure.(15) Similar results are anticipated for other initial velocity distributions. Some curves for F2 and F3 appear to end abruptly. This is due to either y becoming imaginary, such as (21) near -2 = 24, or n switching to the next integer, e.g. (32) and (33) at -2 = 124. If n were allowed to vary continuously, many of the discontinuities would be eliminated. It is possible that some of the irregularities might be eliminated by retaining higher order terms in Malkin's method.

Figure 13 is the locus of the smallest values of e, i.e. Ecr, taken from Figure 12. If C is greater than Ecr, the asymmetric modes will grow. For X > 4.25, Ecr < 0.31; hence the results are consistent with assuming linear symmetric motion and using perturbation methods. However for X < 4.25, Ecr > 1 and the results are no longer valid. The reasons are the symmetric motion is no longer linear as assumed in the analysis, and the perturbation methods are based upon e being small. The three peaks are near the values of X for which bo vanishes. This is equivalent to the condition for modes Zo and To to be inextensional (see Section B, Chapter II). For X > 4.25 we conclude that dynamic instability is initiated by the growth of the asymmetric modes. This is also likely for the range 2.85 < X < 4.25 but this is not substantiated by the present analysis. An analysis using the nonlinear symmetric motion would probably verify this. Figure 15 must not be misconstrued as representing the initial velocity which produces a symmetric buckled configuration via asymmetric motion. An analysis considering finite symmetric and asymmetric motion would have to be considered to study this. What it does represent is the initial velocity necessary to excite the asymmetric modes.

-52Figure 1. Coordinate System and Shell Geometry. K2 K1(2 XI XI Before Deformation After Deformation Figure 2. Location of a Point Before and After Deformation.

1.2 1.0= X=0 X=2.5 I 0 0.0.4o.20.i 0.3 0.5.7 0.9 z,,(o) -0.2 - -0.4-0.6-0.8- -=4 -I.2 - -I.4 Figure 3. Symmetric Mode Shape Zo.

1.2 1.0 X 5 0.8 0.60.40.2zl(O ) m I.O e z:(O) 0.[ 0.2 0.5 0.6 0.7 0.8 0.9 -0.2 - -0.4 - X 7 -0.6 - = - X=:5 -X=:o -0.8 - -I.0 -:2.5 X=3 -1.4 - X=4 -Figure Symetric Mode Shape Z.6 Figure 4. Symmetric Mode Shape Z1.

1.2 1.0 0.80.6X=8 0.40.2z (8) o o.I 0.3 0.4 0.5 0.7 0.8 0.9 Z,(O) -0.2 -0.4 - \ - km -0.6- - -0.8 - -I.0-1.2 X=15 -1.4 - -I.6 - -1.8 Figure 5. Symmetric Mode Shape Z2.

1.2 1.0 0.8 0.6 0.4 0.2 z3 (8) 1.0 Zto) 0.1 0.2 0.3 0. 0.5 0.6 0.8 0.9 -0.2 - -0.4 0-0.6-, -0.8-1.0 - X=8 -1.2 -1.4 Figure 6. Symmetric Mode Shape Z3.

1.4 1.2 1.0 0.8 0.6 0.4 0.2 Nc 9 V0 30 N_-0.2 -0.4 -0.6 -0.8 -1.0 ------ X =68.5 -1.2- -- - =I00 -I,4 - — X-=178 Figure 7. Nondimensional Normal Force Versus Nondimensional Time

( / - N. _____C -t-- 22\ 3, L L / \'O 1 20 40 -3- A=68.5 -4 - =100oo -5 -.- 2=178 Figure 8. Nondimensional Central Deflection Versus Nondimensional Time

-592.0 I.8 1.6 1.4 1.2 8max 1.0 max 0.80.6 - 0.4 0.2 0 0 1.0 2.0 3.0 Figure 9. Maximum Response Parameter Versus Initial Velocity Parameter.

1.4 1.2 Ecr 1.0 o 0.8 I I I I 2 3 4 5 6 7 8 9 10 2 X Figure 10. Critical Velocity Parameter Versus Shape Parameter for Symmetric Buckling.

12 1 1 1 8 2j+i n 6 4 x n= 2 o ~~~II I IIII 0 2 4 6 8 10 12 14 16 jx Figure 11. Asymmetric Resonant Conditions.

0.3 (21) (11) (22) (32) (12) (22) (32) (12) (The first digit in parenthesis is the subscript of the respective function and the second digit is n) (3') 0.2 (33) (13) 4E o) (23) (ii) (31) 0.1 Q- ) /1 (33) 50 I00 150 200 Figure 2~1)~2 x Figure 12. First Zero of F1, F2, F3 Versus Shape Parameter.

0.3 0.2 Ecr OI~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~O 0.1 0 50 100 150 200 2 x Figure 13. Critical Velocity Parameter Versus Shape Parameter.

APPENDIX A ORTHOGONALITY OF THE RADIAL MODES The orthogonality of the linear radial mode shapes is shown. No such property can be shown to exist for the tangential mode shapes due to neglecting the tangential inertia. If the standard separation of variables method is applied to (2.1) by assuming, = R(G) T(T) C(G) T(T)) then C -R = (A. la) a2(R'+ 2R) + 2_ -2 + l) R - C O 0 (A.lb) where R, C and P are the radial mode shapes, tangential mode shapes and eigenvalue, respectively. Let Pm and Pn be two distinct eigenvalues, then (A.1) may be written as Cm - Rn = 0 (A.2a) 1/ / n C2(Rm + 2Rm) + (c2- G2p2+ 1) Rm - Cm 0. (A.2b) Multiplying (A.2a) by Cn and (A.2b) by Rn and subtracting the results yields C C + R C + a R2R+2RR) (a C Rn. (A.3a) -Cm Cn + Rm Cn + (2(Rmn+2FmR n) + (a2 2 p2m+l) RmRn - Cm Rn = 0. (A.3a)

Interchanging m and n gives - CnCm+ RnCm + a2RnRm+ 2RnRm) + ( n2-2pn+)1) RmRn - CnRm = 0. (A.3b) Subtracting (A.3b) from (A.3a) and integrating by parts from -B to P yields 2 2 2, 2 (EP -p m) RmRndg + a2 [(Rm+ Rm) Rn - (Rn+ Rn) Rm [(R Rm) Rn] (R+ Rn) R] - (Cm Rm) Cn (Cn Rn) Cm _ (A.4) If the boundary conditions (2.8) are introduced into(A.4), it reduces to (P2 _ p2) Rm RndG = 0. Since Pm and Pn are assumed to be distinct, then Rm Rn dg = O; m l n.

APPENDIX B EVALUATION OF THE INTEGRALS cl THROUGH C9 The values of the integrals (3.6) are C = S3 [ sin _ cos x][1+ ( sin2X - cos sinn i tanh X X tanh X s inh2J 2 sinX s 1 sin. + A1 C2 = 6 S2S1 tanh tanh k+1 ( Tsin k -cos k) + sin - + 2S A 1 +3 sin2 ) tarhta nh k -cos X sin ) -_ sinh2 k sinh2J 2 2 cos 1 sink Al c3 = 6 S S1 A1..+ (co7s _sin + - cos x) L 1-1 ~- 1+1 tanh X i 1+1 21 c S4%K4 [ ( 4 sin2 % + sin4 % - sin X eos X - 3 sin5\ cos )< 4 [ 4\ sinh X sinh4 7| 20x 10k _in_ 78 + 71 sin2 _ 126 sin2X%+ sin3 78 cos X + 63 sin X 20%tnh178+Yls20inx h- - I + 20Xtanh X sinh2 \ I 20Ox sinh2 k tanh3 I -66 -

c6 =- - 6 s3s 3 V-1 (1-I1) cos3 - (3 - C11) cos X + sin x 06 -1 6 S S~Qu 1ij (iv) (9-4) 2(1+j1) tanh X [ sin fX cos + 2( + 2) sin2X cos x + 2)2+4 tanh X 1 1 ( 1 - r [sin k cos 2?-2(~ -2)sn 1 sin2k cos 2 tanh J 2(~1-1) sinh2k + 1 [ (7 -1) ( 2 + 1 ) sin2k cos X + 4 sin,3 4(-1-1)2 + 16 sinh 2 tanh X 1' +1) 2 + s sin2k cos k - 4(41+1)+ 16 1sinh2X tanh X + sin5. [ p3 +1 _ 13 it 2S SlA1X<[6Jlsinx ( ln+9)(1+) tanh3X sinh X tanh %. (i1+. (9+ i 1) L - 3 (3 + o) cos X + (1 + 11K) \1 sin3X + 3(1 + 11) cos3x] + 3 S3Sj1Aj1.i |2 sin %L(X1 tanh X -1)_ (1 + tanh X) 4 sin 2 cos X 2 tanh? +1)2+ 4[ + (T +1) sin +) cos 2x] + 2tac X (+ ( )+ l-1) sinc.- s 2 1 - JJ +2~-, 3 S3S+ -4 3 S"S1A1l%~N. f sin2? 2 1 2 (sn +2)7 ta n - ( 2+anh X) sih sin. J1cl s X 2 (L1 - tanh ks inh2 [

+ sin22 t[ + 1 [ (NJ1-2) sin X- cos i (12) 2+1 tanh X sinh X 2 sin2 X (Nflsi - cos )} + S3S1A1\;J 1 [1 (ll+l) sinh23 2 H1-9 + 4 3 sinn3 ssinnh2 ( tanh3X sinh2X tanh jX i] _, ~ [ (N71 tanh X -1) sin3k ] sinh2% tanh c 1ss3K x _ (11-2) sin X cos. + (311+ 4) sin2 sin x? 2 2 2(Ki1-1 ) 2(G1+1) tanh X _ k sin- 11 X sin2% 1 1(2N/jT+l) sin cos _sin2 1] 2 sinh2 (21+1)2+1 tanh X (1-21)2+ [(1-2N/+1) sin k cos _ sin2 I] 3 SsA 111 1 1 ~tanh'2 {1 (N/'-1+2)2+4 [Icos 2 X + 2 (N/K+2) sin cos X ( 2)2+1[ 2X - 2 (NJ1-2) sin X cos ] 4 [ (tanh \ - /717) sin X cos X ( 2I -1)tanh X

-69+(13/2 tanh -12) s+ 2 in2 + 2 S 2 223 (N1+2)2+ 1 tanh X sinh2 X + (N/j -2) [2- 2 + 1 ] sin2 X si + 3 lAlP23Ili ( 1-2 )2+ I1 tanh % sinh 2 J / sinh 2 (Ni,-2) +I l-J 1 sin X cosh 2/;1 cos 2 + 2 [1 + 1 1[(2N+1) sin2 e2 - ( ( 1 + 21) (1 1+ 1) ( 21+ 1)12 tanh IJL - sin X cos ] + 2 1 (2 -1) sin2\ - sin X cos ] 51~1 Cil/2I 2+ 1 (2 1i-_1)2+1 tanh %J L + 2([-1) L'~ sinh2l tanh \ 2 \ sinh2X'8 = S S3 3 13/2 2 (cos + (7?11+1) sin h 8 = S S~X3~i tanh k 1(L -1) (91l-1) (1+1) (9111+ 1) + 3 A1 2(1 sin X - cos x) - (2N1 +1) cos x - ",sin x 2 1 + 1 (2 N1lJ+ l)2 + 41 + (2 N/71-1) cos X + \/1 sin X + 2(N/ltanh X -1) sin X (2 + 1) + tanh ) sin ( 1-1 tanh ) sin (1+1)(1 +tanh) si ( -)( - nh h () sin X (N/Z1+ 1)2 + 41l tanh k (41- 1) + 41ll tanh kj

A2- cos + 2 sin T-) cos 2-(sin - 3A2 (4/V1-) cos x= + 2Iiiisnx - (i+1) cos X-2 2iisin. 1- ( 1)2+ l + 1)2+ 4 1 cos % + (2N/ij+l)(1 + tanh N) sin X (2Ng1-1)(1 - tanh X) sin X ( 1-t1) (2N/1+1)2 + C1 tanh X (2f1'1) + 41 tanh xJ A3 F3 /lsin X cos X 3;V1 tanh k -1 sin x 9 11 + 1 9 1L - 1 tanh X 4 4 2.2 4 7 s1% ~1 x 32 Al + 3A 1- 1+ A1 4 15 5 3 In the above tanh (3T/2) has been replaced by unity and [cosh 3.t/2] has been dropped compared to unity.

REFERENCES 1. Hoff, N.J. and Bruce, V.C., "Dynamic Analysis of the Buckling of Laterally Loaded Flat Arches", J. Math. and Phys., Vol. 32, No. 4, (Jan.) 1954), p. 276-288. 2. Bolotin, V.V., Dynamic Stability of Elastic Systems, Holden-Day, Inc., San Francisco, 1964. 3. Humphreys, J.S. and Bodner, S.R., "Dynamic Buckling of Shells under Impulsive Loading", Proc. ASCE, Vol. 88, No. EM2, (April,1962), p. 17-36. 4, Fulton, R.E., "Dynamic Axisymmetric Buckling of Shallow Conical Shells Subjected to Impulsive Loads", J. Appl. Mech, Vol.32, No. 1, (March, 1965), p.129-134. 5. Budiansky, B. and Roth, R.S., "Axisymmetric Dynamic Buckling of Clamped Spherical Shells'" NASA TN D-1510, (December, 1962), p. 597-606. 6. Burns, J.J., Popelar, C.H., and Foral, R.F., "Experimental Buckling of Thin Shells under Transient Pressure Pulse Loading", Presented to American Rocket Society Launch Vehicles: Structures and Material Conference, (April 3-5, 1962), p. 2425-62. 7. Love, A. E.H., A Treatise on the Mathematical Theory of Elasticity, Dover Publications, Inc., New York, 1944. 8. Archer, R4R., "Small Vibrations of Thin Incomplete Circular Rings", Inter. J. Mech. Sci. No. 1, (Jan., 1960), p. 45-56. 9. Sokolnikoff, I.S., Mathematical Theory of Elasticity, McGraw-Hill Book Company, Inc., New York, (1956), p. 86. 10. Timoshenko, S., Strength of Materials Part II, D. Van Nostrand Company, Inc., New York, (1930), p. 429. 11. Goodier, N.J. and McIvor, I.K., "The Elastic Cylindrical Shell under Nearly Uniform Radial Impulse", J ppl Mech., Vol.31, 1964. 12. McIvor, I.K., "Dynamic Stability of Axially Vibrations Columns", Proc. ASCE, Vol. 90, No. EM6,(Dec., 1964), p. 191-210. 13. Stoker, J.J., Nonlinear Vibrations, Interscience Publishers, Inc., New York, (1963), p. 36. 14. Timoshenko, S., Theory of Elastic Stability, McGraw-Hill Book Company, Inc., New York, (195), p. 234. -71

-7215. Weinitschke, H.J., "Asymmetric Buckling of Clamped Shallow Spherical Shells",?NASA TND-1510, (1962), p. 481-491. 16. Schreyer, H.L., "On the Theory of Elastic Stability", Ph.D. Thesis, University of Michigan, 1965. 17. Struble, R.A., Nonlinear Differential Equations, McGraw Hill Book Company, Inc., New York, 1962. 18. Struble, R.A. and Fletcher, J.E. "General Perturbational Solution of the Mathieu Equation", J. Soc. Ind. and Appl. Math., Vol. 10, (1962), p. 314-328. 19. Series of Papers by Valcev, K.G., "On Hill's Method in the Theory of Linear Differential Equations with Periodic Coefficients", J. of Appl. Math. and Mech. (PPM), Vol. 24, (1960), p. 877-902, 1493-1505, Vol. 25, (1961), p. 460-466. 20. Hsu, C.S., "On the Parametric Excitation of a Dynamic System Having Multiple Degrees of Freedom", J. Appl. Mech., Vol. 20, (1963), p. 367. 21. McLachlan, N.W., Theory and Applications of Mathieu Functions, Dover Publications, Inc., New York, (1964), p. 3. 22. Den Hartog, J.P., Mechanical Vibrations, McGraw-Hill Book Company, Inc., New York, (1956), p. 339. 23. Bogoliubov, N.N. and Mitropolsky, J., Asymptotic Methods in the Theory of Nonlinear Oscillations, Hindustan Publishing Corp., Delhi - 6, 1961. (Translation from Russian) 24. Minorsky, N., Nonlinear Oscillations, D. Van Nostrand Company, Inc., New York, (1962), p. 124-127. 25. Malkin, I.G., Some Problems in the Theory of Nonlinear Oscillations, (A translation of Nekotorye Zadachi Teorii Nelineinykh Kolebani, Moscow 1956) United States Atomic Energy Comission, AEC-tr-3766, p. 220-225.