THE UNIVERSITY OF MICHIGAN COLLEGE OF ENGINEERING Department of Engineering Mechanics Technical Report DYNAMIC SNAP-THROUGH STABILITY OF THE NONLINEAR, ELASTIC, PLAIN STRAIN, FREE VIBRATION OF A SHALLOW CYLINDRICAL SHELL Lee J. Ovenshire Ivor K. McIvor ORA Project 03021 supported by: NATIONAL SCIENCE FOUNDATION GRANT NO. GK-10725 WASHINGTON, D.C. administered through: OFFICE OF RESEARCH ADMINISTRATION ANN ARBOR October 1969

This report was also a dissertation submitted by the first author in partial fulfillment of the requirements for the degree of Doctor of Philosophy in The University of Michigan, 1969. ii

TABLE OF CONTENTS Page LIST OF FIGURES v NOMENCLATURE vii ABSTRACT xv CHAPTER I. INTRODUCTION 1 II. MATHEMATICAL FORMULATION OF THE PROBLEM 2.1. General Comments 5 22. Description of the Shell Configuration 5 2.3. The Validity of the Geometrical Approximations 12 2.4. The Nonlinear Partial Differential Equations of Motion 17 2.5. Conditions Just After an Impulsive Pressure Loading 21 III. SUFFICIENCY CONDITIONS FOR DYNAMIC STABILITY 24 3.1. General Comments 24 3.2. The Critical Configurations 25 3.3. Balanced Support Critical Configurations 28 3.4. Simply Supported Symmetric Critical Configurations 50 3.5. Simply Supported Nonsymmetric Critical Configurations 6 3.6. The Potential of the Simply Supported Critical Configurations 41 IV. INFINITESIMAL VIBRATIONS 49 4.1. General Comments 4c 4.2. General Solution of the Linear Equations of Motion 49 4.3. Normal Mode Shapes for End Supports Elastically Restrained Against Rotation 55 4.5.1. The boundary conditions 4.3.2. The characteristic equation 55 4.3.35 The linear radial displacement mode shapes Ei 4.3.4. The tangential displacement, tension and bending moment 4.4. The Infinitesimal Motion for Supported Ends 75 4.5. The Geometrical Parameter X 92 iii

TABLE OF CONTENTS (Concluded) CHAPTER Page V. METHOD OF SOLUTION OF THE NONLINEAR PROBLEM 94 5.1. General Comments 94 5.2. Representation of Nonlinear Displacement and Deformation Quantities 94 5.3. The Functions-t and ~ and the Constant Array N 96 5.4. Derivation of the Equations Governing the Nonlinear Motion of the "Fourier Coefficients" 99 VI. THE NONLINEAR MOTION 103 6.1. General Comments 103 6.2. Measures of Response 104 6.3. Response Parameters for the Critical Configurations 105 6.4. The Nonlinear Response 107 6.5. Topological Interpretation of Results 119 VII. CONCLUSIONS 126 APPENDIX A. IDENTITIES FOR THE CRITICAL CONFIGURATIONS 128 APPENDIX B. ORTHOGONALITY OF THE RADIAL MODES 131 APPENDIX C. QUANTITIES NEARLY INVARIANT TO THE SEMI-OPENING ANGLE 1355 LIST OF REFERENCES 136 iv

LIST OF FIGURES Figure Page 1. Cross section of the shell in its undeformed state. 6 2. Shell cross section for general plane strain deformation. 7 3. Circular arc configurations. 13 4. Cross section with external loads. 19 5. Plots of Y and, functions. 34 2 6. Plots of cos q and (qX). 40 7. Plot of V versus \ for both nonsymmetric and symmetric critical configurations. 45 11`1 4 2 8. Plots of 8V(X/z) versus X /4. 47 9. Comparison of Hsu's stability sufficiency curve with present curve V N. 48 10. Plots of the characteristic function A for symmetric modes for P =.25 and X = 3, 6, and 15 for the simple support case; also, plots of the tangent function. 60 11-16. Plots of ratios, Zn(e)/Zn(0), of the n'th symmetric modes to their value at 0 = 0 for 3 =.05 with the shape parameter X as the parameter (n = 1,2,...,6). 67-72 17-20. Plots of Tn/A2n versus the dimensionless angle r = E/P with X as the parameter for the simple support case for the n'th mode and a semi-opening angle P =.05. 74-77 21-26. Plots of a dimensionless midsurface strain Nc/vo (= -~AVGC/vo) versus T for P =.05 for X values of 3, 7, 8.28, 10, 13.35, and 15 (ten modes were retained). 50-35 27-32. Plots of a dimensionless central deflection (0, T)c/vo versus T for P =.05 for X values of 3, 7, 8.28, 10, 13.35, and 15 (ten modes were retained). 36-91 335. Plot of X versus P with 1/h as the parameter. 9. v

LIST OF FIGURES (Concluded) Figure Page 34. Response curves for X = 4,. =.05 and c = 2.1 for: (1) k = 0 and (2) k = 0.1. 109 35. Response curves for X = 4, 3 =.05, and k = 0 with six symmetric modes for: (1) c = 1.2 and (2) c = 2.22. 110 36. Response curves for X = 4, 3 =.05, and E = 6.44 with six symmetric modes for triangular and lowest radial displacement mode shaped impulses. 111 37. Response curves for x = 5, and P =.05 for: (1) c = 1.75 and k = 0; (2) c = 1.9 and k = 0; and (3) c = 1.85 and k =.1. 112 38. Response curves for X = 6, k = 0, and six symmetric modes for: (1) P =.05, C = 1.1; (2) P =.05, = 1.2; (3) =.01, e = 1.2; and (4) e =.1, E = 1.2. 11 39. Response curve for X = 6, P =.05, k = 0, e = 1.4 (E = 12.75), and six symmetric modes. 114 40. Response curves for X = 6, B =.05 and e = 1.55 for: (1) k = O. four symmetric modes; (2) k = 0, six symmetric modes; and (3) k =.1, six symmetric and six asymmetric modes. 115 41. Response curve for X = 7, P =.05, k = 0, and c = 1.0 (E = 12.3) with a six-symmetric mode representation. 116 42. Response curves for X = 10, P =.05, and k = 0 for: (1) e =.81, eight symmetric modes, and (2) c = 1.2, nine symmetric modes. 117 43. Plots of first swing 5MAX versus c with X as the parameter. 120 44. Plots of the critical initial velocity parameter CCR versus the geometrical parameter X for immediate dynamic snapthrough instability. 121 45. Plots of energy versus X for a stability sufficiency condition and for dynamic snap-through conditions with a cosine impulse. 122 vi

NOMENCLATURE ~~(')> (1)0:-( ) () Quantity pertaining to the n'th "symmetric-like" mode ( ) Quantity pertaining to the n'th "asymmetric-like" mode A,.*.*,A Coefficients of the linear radial displacement mode shapes with the mode number unspecified A ln5.. A Coefficients of the n'th linear radial displacement mode shape B1..,B 4 Coefficients for the critical configurations (z = B1 sin qF + B2 cos qF + B sin rF + B4 cos rr + N/a) Bln...,Bn, Coefficients used in calculating, (f): (Bl..,B ) = (A,-v A, A,wA ) B ln.B 4n vn ln n 2n n n'n 4n C (1) Abbreviated notation for cos r in the critical configuration analysis (2) Abbreviated notation for cosh w in the analysis of infinitesimal vibrations (3) Center of curvature for S and S C* Center of curvature of midsurface when deformed into a circular arc C Abbreviated notation for cosh w n n E Young's modulus of elasticity E Total nondimensional energy, T+V E Total invariant nondimensional energy, T+V F(Y) Right side of the configuration vector equation Y=F F Right side of the nonlinear system Y =F, n=1,2,... G(q) = 0 Characteristic equation for simply supported symmetric critical configurations vii

H(0) Midsurface rise of the undeformed shell at 0 H Maximum undeformed midsurface rise, H(O) o H* Midsurface rise of the deformed shell I(0) Dimensionless impulse (=z(0,0)= O p dT) e I..(m,n,r) Functions used to calculate Q (r) 1j mn K Elastic bending stiffness when the two ends are balanced, i.e., K = K = K K K 2 Left- and right-hand elastic bending stiffness constants L Span of the arch midsurface M Moment resultant per unit length of cylinder M Dimensionless bending moment (7M = bending part of strain) M External nondimensional bending moment applied to an edge N Tensile force resultant per unit length of cylinder N Midsurface strain (and a dimensionless tension) N External nondimensional tension applied to an edge e N Array used in calculating N in the nonlinear analysis (=n ) m m N Nonlinear part of N O Material point on a cross section S of the midsurface o 0* Spatial point that 0 occupies at time t P A material point on Sp radially inward from 0 by a distance z P* Spatial point that P occupies at time t Q Material point on S a differential arc distance dl from P o P QA* Spatial point that Q displaces into Q Nondimensional shear resultant externally applied to the shell at an edge viii

R* Position vector CP* that locates P* after deformation R* Position vector CO* that locates O* after deformation o S (1) Abbreviated notation for sin r in the critical configuration analysis (2) Abbreviated notation for sinh w in the analysis of infinitesimal vibrations (3) Part of the exponential representation Z = Ae S Abbreviated notation for sinh w n n S Cross section of the undeformed midsurface, with radius a 0o S* Curve that S deforms into o o S Undeformed circular arc cross section concentric to the midsurface with radius (a-z) S* Curve that S deforms into P P T (1) Abbreviated notation for tan r in the critical configuration analysis (2) Abbreviated notation for tanh w in the analysis of infinitesimal vibrations (3) Kinetic energy of the cylinder per unit length of cylinder (4) Time varying modal "Fourier" coefficient (mode number is unspecified) T Nondimensional kinetic energy, T(l-v )/(Eah) f" - 2 T Invariant nondimensional kinetic energy, T/(a2 ) V Strain energy per unit length of cylinder 2 V Nondimensional strain potential energy, (l-v )V/(Eah) V Invariant nondimensional strain potential energy, V/(pa ) X or X* Distance from 0 or 0* to the undeformed plane of symmetry Y(T) Configuration vector whose components Y consist of Y i=T n 2i i and Y2i_=Ti i=l,2..,N; n=-,2,...,2N Z,Z n Linear radial displacement mode shape ix

a Radius of the undeformed midsurface a* Radius of midsurface when deformed into a circular cylinder Coefficients for the critical configurations (a = (~1 sin qF + ~2 cos qF + 3 sin rF + I4 cos rF + 1)N/a ) c (1) Abbreviated notation for cos q in the analysis of the critical configurations (2) Abbreviated notation for cos v in the analysis of infinitesimal vibrations (3) Wave speed, E/(p(l-v ) c Abbreviated notation for cos v n n dli Differential arc length PQ along Sp dlp Deformed arc length P*Q* of dp along SP f Shape of asymmetric part of initial velocity distribution used in several examples ( - 3 r(l-F )/2) f Function appearing in infinitesimal vibration analysis g Function appearing in h and in an approximation to G for the analysis of critical configurations (= 2q /3-4q /K ) h The shell thickness h(q) Function used to calculate B for nonsymmetric critical configurations (h = g-l/cos q4 1B1 q /P ) 1, Unit horizontal and vertical base vectors (J points upwards) n,n. Coefficient for midsurface strain for infinitesimal vibration of a linear mode (N = n sin (jT + 0)) n Unit normal to the undeformed midsurface S at 0 (or to Sp at P) A n* Unit normal to S* at P* p (1) Inward radial pressure (2) Part of the separation constant in the derivation of linear modes p Position vector OP x

2 bPe ~ Nondimensional inward radial pressure, (1-v )ap/(Eh) p* Position vector O*P* / — -i q Eigenvalue for a critical configuration, VA + -1 q*(K) Larger of two values of q that satisfy'= 1 (used to determine A in the critical configuration analysis) r Quantity used in the critical configuration analysis, /q = A - VA - s (1) Abbreviated notation for sin q in the critical configuration analysis (2) Abbreviated notation for sin v in the analysis of infinitesimal vibrations s Abbreviated notation for sin v n n t (1) Real time (2) Abbreviated notation for tan q in the critical configuration analysis (3) Abbreviated notation for tan v in the analysis of infinitesimal vibrations t Unit tangent vector to the undeformed midsurface S at 0 (or to S at P) Unit tangent vector to S* at P* P u Quantity in the infinitesimal vibration analysis, B p = (v2_P2) /2 v (1) Tangential component of 5, i.e., v = 6 (2) Eigenvalue for the linear modes of vibration ( p+i ) v,v Eigenvalues v for "asymmetric-like" and "symmetric-like" linear modes of vibration, respectively v Radial component of midsurface velocity at the center o w (1) Inward radial component of 5 i.e., w = n (2) Quantity in the infinitesimal vibration analysis, Sp-il = Sv2: x,y Cartesian components of the midsurface displacement, i.e., _ A A 6 iX - Jy o xi

z Distance from midsurface to generic point P b (q) Difference between slopes of / and X used in a proof r (rE) Subscripted functions in the nonlinear analysis, / Z' Zd -1 m n Part of approximate characteristic equation for simply supported critical configurations (= qg) C ~ Lim = 2q5/3 -*woo r (v) Right side of the characteristic equation t = f(v) for "asymmetric-like" linear modes V (v) Right side of the characteristic equation t = (v) for "symmetric-like" linear modes 2 1'4 Array appearing in the nonlinear analysis, a m m dF ru _-1 n m JH ~ Part of the approximate characteristic equation for simply supported critical configurations, q sec2 q -tan q'q Nonlinear part of differential equations governing the motion of the "fourier" coefficients T q L y(qu) Part of the equation cos q = whose roots determine where B is positive for nonsymmetric critical configurations 1'2 (r) Function used in the nonlinear analysis IT19 O2 Portions of expression for V r Normalized angle coordinate, 5/5 A A constant coefficient for the differential equation governing the shape of the critical configurations T9.n Tangential displacement mode shapes ~a ~ A slenderness ratio, h/(aX2 ) pB0 ~ Semi-opening angle before deformation pB"A ~ Semi-opening angle for shell when deformed into a circular arc 7 Dimensionless distance from the midsurface, z/a xii

1 1 1 1 6(T) Response parameter, 2 / a: d' /( / H dF) 2-1 -1 MAX Maximum value of 6 attained in the time interval T = 0 to MAX the next zero of 6 6 Displacement vector 00 for 0 C (1) Lagrangian strain along S (2) Impulse parameter, v n2/(c P2) (3) Dimensionless time duration of impulse ER Value of c for MAX 1.5 during first swing CR MAX E Midsurface strain, N Dimensionless inward radial component of midsurface displacement, 6 ~ n/a o 1 1 CAVG Linear part of N (= - 2 f ~dr) 0 Clockwise angle between plane of symmetry and the ray CP before deformation 0* Clockwise angle between plane of symmetry and the ray C*P* after deformation into a circular arc ~~X Geometrical shape parameter, i/H 0 41/4 (H /h)l/2 N Value of X for which a new pair of symmetric critical configurations arise as X is increased XN Value of X where the N'th nonsymmetric critical configuration first arises xN Value of X for which the highest symmetric critical configuration has the root q = Nr gJ~ ~ Angular dimensionless frequency for infinitesimal vibration of T ( = (v2-_2)/%2) v Poisson's ratio ~a ~ Tensile stress in tangential direction T Dimensionless time, ct/a Phase angle in T = T* sin (U.T + ) xiii

Phase angle in T - T* sin (u T- + 4 ) n n n n 9~f ~ Dimensionless tangential component of midsurface displacement, t * S /a 0 w Angle of counterclockwise rotation of a radial section xiv

ABSTRACT The dynamic snap-through stability of the nonlinear, elastic, plane strain, free vibration of a shallow circular cylindrical shell is studied using equations of motion based upon thin-shell theory. Critical configurations, i.e., solutions to the equations of motion that are independent of time, are derived for general linear elastic bending restraints against rotation at the boundaries. The strain potential energy for configurations with no bending restraints is presented. Based on the energy of the critical configurations, a stability sufficiency condition is determined. If the total energy of the shell never exceeds a certain level, the shell cannot snap through. It is found that stability sufficiency conditions based upon thin-shell theory and shallow-beam theory are qualitatively the same. A stability sufficiency condition does not, of course, determine the magnitude of initial conditions at which instability occurs for a given spatial distribution. Here it is necessary to integrate the equations of motion directly. In the present investigation, the conditions of dynamic snap-through are presented for cosine impulsive loading of a previously undeformed and motionless, simply supported shell. As a preliminary to the latter investigation the linear mode shapes and the infinitesimal motion are determined for the case of supported edges with general linear elastic restraints against rotation. Plots of the characteristic function and linear mode shapes and of central deflection and tension versus time are presented for the simple support case. A beating effect is noted for certain geometries. xv

Nonlinear ordinary differential equations are obtained through Galerkin's method by representing the displacements as a sum of linear mode shapes, treating the time varying coefficients as generalized coordinates, applying Hamilton's Principle, and integrating out the spatial dependence. For a given set of initial conditions the time varying coefficients are integrated numerically. The required number of modes of representation is examined and found to depend upon both the size of the impulse and a geometrical parameter. It is found that previous investigations of the shallow arch did not use enough modes. Numbers of modes adequate for the convergence of a response parameter b during its first half cycle are presented. More modes are probably needed for longer periods of time. Several results of the present investigation are found to be virtually independent of small semi-opening angle P. These include the critical configurations, the infinitesimal vibrations and the nonlinear solution for the first two oscillations of t. However when P is changed appreciably, important quantitative differences in the nonlinear solution appear after longer periods of time. This suggests that investigations of long-time, large-amplitude motion of other shallow snapping elements may also need geometrical parameters that characterize ratios of rise to span in some way. Dynamic snap-through is observed only for geometries in which nontrivial critical configurations exist. It is found that a delayed snap-through phenomenon sometimes occurs for energy levels insufficient for immediate snapthrough. In this investigation the energy levels of snap-through are always observed to be above the second lowest of the nontrivial critical configuraxvi

tions. A heuristic argument in terms of the topology of the potential energy surface is presented to give both a general interpretation of results and an implication of truncating the number of terms in the modal representation. This argument suggests that the appropriate number of modes retained should be adequate to represent the critical configurations. xvii

CHAPTER I INTRODUCTION Dynamic stability of structures has received considerable attention in recent literature. There are a great variety of dynamic stability problems even when material behavior is restricted to the elastic range. The motion of the structure may become unstable under time-varying loads which may be periodic or aperiodicc Also the motion of a freely oscillating structure may be unstable in the sense that small changes in initial conditions give rise to substantially different motions. A number of authors have investigated problems in which the loading is applied throughout the motion of interest. Loadings considered include those (1-3)* which are suddenly applied and maintained and those which have periodic (4) variations about a mean.( The initial value or free vibration problem has also received much recent attention. A prominent example is the impulsively loaded structure. If a high intensity loading is applied over a sufficiently short time, an easily determinable velocity distribution is imparted to the structure during a process involving negligible displacement. The conditions just after completion of the impulsive process serve as initial conditions for the free vibration prolem. The ensuing motion may be unstablem in the sense that small changes in initial conditions result in significantly different motion. Numbers raised in parentheses designate references. 1

2 In a recent paper Hsu considers dynamic instability "in the large" for free motion problems in terms of trajectories in functional configuration space. To obtain an analytical stability criterion is in general difficult. Hsu presents sufficiency conditions for stability. The motion is stable if the energy imparted to the structure is less than the potential energy associated with the first critical or static equilibrium configuration encountered by an expanding level energy surface in the functional configuration space. Hsu illustrates his approach with the sinusoidal arch problem using curved (6) beam equations derived by Fung and Kaplan. Several of Hsu's ideas are used in the present study. Chapter III is devoted to obtaining a sufficiency condition for a simply supported circular cylindrical shell undergoing plane motion. The nonlinear equations of motion are derived in Chapter II and are based upon thin shell theory. It is found that using the different equations gives quantitatively different results. This approach obtains an energy level for initial conditions of any shape below which the motion is dynamically stable. A sufficiency condition does not, of course, determine the magnitude of initial conditions at which instability occurs for a given spatial distribution. Here it is necessary to integrate the equations of motion directly. The usual approach is to represent the displacements in terms of one or more spatial functions whose magnitudes are time varying coefficients serving as generalized coordinates. A system of ordinary differential equations in these coefficients is obtained by a modified Galerkin procedure. * (7) See, for example, Kantorovich and Krylov.

5 Several authors have used this approach in treating the dynamic snapthrough of shallow structures. The various investigations differ primarily in two respects: (1) the derivation of the partial differential equations of motion, and (2) the spatial functions chosen for the representation. To the author's knowledge, with the exception of a study by McIvor and (8) Popelar, all previous studies of the shallow arch are based upon curved (6) beam equations of the type developed by Fung and Kaplan. The partial differential equations derived by McIvor and Popelar are a slightly simplified (9) form of those derived by Goodier and McIvor for the complete circular cylindrical shell undergoing plane motion. These equations are used in the present investigation. Several different sets of representative spatial functions have been used. One group of investigators has selected a single displacement function (one degree of freedom model) which approximates the buckled configuration and sat(10) isfies boundary conditions. This was done by Humphreys and Bodner ) in their analysis of the symmetric buckling of spherical caps and cylindrical (11) panels and later by Fulton for shallow conical shells. Another group has used the linear mode shapes of an equivalent flat element. Hoff and Bruce) used the lowest symmetric and the lowest asymmetric mode shapes of a straight beam to represent the deformation of a sinusoidal arch. To study the spher(12) ical cap, Budiansky and Roth represented displacements as a sum of the five lowest axisymmetric mode shapes of a circular plate. Once the representative spatial functions are chosen, the conditions of dynamic stability are usually obtained by direct numerical integration for various initial conditions of a specified spatial distribution. One or more

4 measures of overall shell deformation are usually computed. If the measures vary markedly for small changes in initial conditions the motion is considered unstable. In most investigations the spatial distribution of initial conditions is fixed. Usually a comparison is made only of the maximum values of the deformation measure attained for the various magnitudes of initial conditions. Fitting into these generalizations are the investigations of Simitses, (1) (2,3) (1) (12) Lock, ( Hoff and Bruce, and Budiansky and Roth among others. This type of analysis is made in the present investigation for the shallow cylindrical shell undergoing plane motion. Using the partial differential equations of motion derived in Chapter II, the linear modes of vibration are obtained in Chapter IV for general elastically restrained bending conditions. In Chapter V the nonlinear displacements are represented by an expansion in the linear mode shapes. Treating the coefficients as generalized coordinates, a system of ordinary differential equations is derived. A detailed investigation of snap-through buckling is presented in Chapter VI for a range of geometries. The initial conditions correspond to an impulsive inward radial pressure of cosine spatial distribution. This numerical investigation is limited to the case of simply supported boundary conditions.

CHAPTER II MATHEMATICAL FORMULATION OF THE PROBLEM 2.1. GENERAL COMMENTS In this chapter a set of nonlinear partial differential equations governing the dynamic behavior of a shallow, circular cylindrical shell is derived. (14) The derivation follows that in Popelar and is given here for completeness. Since the major interest in this dissertation is the motion following impulsive pressure, the appropriate initial conditions for this case are also obtained. The analysis is based on the following assumptions: (a) The external pressure loading is applied radially inward and does not vary along the straight line generator. (b) The deformation is that of plane strain. (c) Material lines normal to the midsurface before deformation remain so after deformation and do not change their length. (d) Damping may be neglected. (e) The normal stress acting on internal surfaces parallel to the midsurface is neglected compared to the other stresses. (f) The material of the shell is elastic, homogeneous and isotropic and remains so during subsequent motion. (g) The shell is shallow. 2e 2. DESCRIPTION OF THE SHELL CONFIGURATION In this section the necessary geometric relations characterizing the shell configuration are derived. A number of geometric approximations are 5

6 introduced in the derivation. The resulting error and limitations of the theory are discussed. The undeformed geometry is shown in Figure 1. Figure 1. Cross section of the shell in its undeformed state. The quantities shown are: a The midsurface radius h The shell thickness P The semi-opening angle H(O) The midsurface rise at 0 Figure 1. Cross section of the shell in its undeformed state. The quantities shown are: a The midsurface radius h The shell thickness P The semi-opening angle H(e) The midsurface rise at e H The maximum midsurface rise o L The span of the arch A general plane strain deformation is shown in Figure 2.

7 S & O Figure 2. Shell cross section for general plane strain deformation. The quantities shown are: S A cross section of the undeformed midsurface, a circular arc of radius a Sp A circular section concentric to SO with radius (a-z) o A material point on S at angle 9 from the vertical P A material point on S radially inward from 0 p p p 5A AA~~ P,n Unit tangent and normal vectors to Sp at P or to S at 0 ~~~~~~~* * 0,P Spertical, P at time t ~~~~~~~~~* * S,S Curves that S and S have deformed into at time t A dl* _ * * Sp A circular section concentric to So with radius (a-z) 0 A material point on S0 at angle 9 from the vertical P A material point on Sp radially inward from 0 C The center of curvature of Sp and So p OP nn Unit tangent and normal vectors to Sp at P or to S. at 0 0,P Spatial points that 0 and P occupy at time t 0X XJ P 0P

8 A A * * t*,n* Unit tangent and normal vectors to S at P P A A t n Unit tangent and normal vectors to S at 0 o o R Position vector CP for P - *- A A w,v The radial and tangential components of 5, i.e., 6 = wn + vt O 0 Q A material point on S a differential arc distance dip from P Q The spatial point that Q displaces into dip The deformed arc length P Q of dip The Lagrangian strain along Sp is dl - dl _l P (2. 1) 2dlP The required squares of differential lengths may be written as *2 _*t, * dip (= e(R de ) d) (2.2a) and 2 2 2 dip = (a-z) 2d (2.2b) where the prime denotes differentiation with respect to 8. From Figure 2 we see that _* A - - R = na + oo + p (2.3) P_* AX(A* p = zn (2. 4)

9 It is convenient to introduce the following dimensionless quantities: 7 = z/a 5 = w/a (2.5) r = v/a Substituting Equations (2.4) and (2.5) into (2,3) yields -* A A*(2 R = a[an(]-l) + t + 7n ]. (2.6) The derivative of this with respect to e is R = a[n( + -+1) + 7n ] (2.7) A* A* Under Assumption (c) on page 5, the unit vectors t and n are parallal to..XA~~~~~~-.-^ab A-e- a-t thtisr t and n for every 7. Thus the vector R /a can be evaluated at the midsuro o face to give a vector T = (-+)+ n('g ) 0 * * -that is tangent to S at P. We will see later that the following inequalities are valid for the deformations anticipated: 1 ~ It'- >> ~ c'+ l _- ~ ~~~~A* A*According to these inequalities, T is nearly a unit vector and t and n may be approximated by t = t - A A2.8a) t = t - wn (2.8a)

10 and IA* " A A n = n + wt (2.8b) where - (='+*). (2.9) The quantity (e is the infinitesimal counterclockwise rotation of the radial section. By substituting Equations (2.8) and (2.9) into (2.7) and rearranging terms, we get R /a - n(l-7)W +( l +;+')] Using this and Equation (2.2a), we get dl 2/(ad)2 = (1-) 22 + [-( "+f)+(-y) (l+'"+') ]2 (2.10) Combining Equations (2.10) and (2.2b) with (2.1) and expanding l/(l-y) yields 2 25 2 E = () -1)/2 + [(, +')- (, + )( 7+ +... ) +] /2 The second order displacement quantities arising from the square bracket are negligible. By dropping these, neglecting r in comparison to i', and retaining only terms through first degree in 7, we obtain the strain displacement relation E = N + 7M (2.11) The ratio,r/~' is of order CL(3 ) for the type of deformation anticipated. This is discussed in Section 2.2.5.

11 where N = (,'-) + () 2/2 (2.12) and M = -(f+0). (2. 15) From the form of Equation (2.11), it is evident that the midsurface strain is N whereas the bending strain is yM. We now show that N and M may be interpreted as dimensionless tension and moment resultants, respectively. The tension and moment resultants per unit axial length, N and M, are given by h/2 N $/2. c dz (2.14a) -h/2 and h/2 M = /2a zdz (2.14b) where a is the tensile stress in the tangential direction. Throughout this thesis the constitutive relation is assumed to be o = eE/(l-v2) (2.15) Using Equation (2.15), Equations (2.14) yield the relations h/(2a) N(1-v )/(Eh) = (a/h) _h (a dy = N (2.16a)

12 and h/(2a) M12a(l-v2)/(Eh) = 3(2a/h) h/() 7 d M. (2.16b) -h 2a)= Equations (2.16) may be regarded as defining relations for dimensionless tension and moment resultants per unit length N and M, respectively. 2.3. THE VALIDITY OF THE GEOMETRICAL APPROXIMATIONS In the previous section a number of geometric approximations were introduced based on the relative magnitude of displacement and displacement gradients. In this section the validity of these approximations is considered. Since the actual nonlinear response is not known, it is not possible to explicitely determine the errors involved in these approximations. We can, however, get an estimate by considering hypothetical, but reasonable configurations. If we neglect tangential inertia, it will be shown that the midsurface strain is spatially uniform. Thus here we only consider configurations having uniform midsurface strain. In particular, let us consider a family of hypothetical shell configurations having circular arc sections each with uniform midsurface strain as shown in Figure 3. The uniform midsurface strain condition leads to aE a E aG _ a~ ** a5 or * * = D,

15 Undeformed midsurface The quantities shown are: a The radius of the undeformed midsurface O A midsurface material point with cartesian coordinates (X,H) and polar coordinates (a,9) H The rise of point 0 n The unit normal vector to the undeformed midsurface () quantities corresponding to () for the deformed midsurface H The midsurface displacement vector 00 () ~atte o rr s o d n o( o h eomdmduf

lh Using this and some elementary trigonometry, we find that sin sine 1 y =a cos 9 - cos - - cos s + cos ( (2. 18a) sin P sin P and sin x = a s sin sin sin (2. 18b) Lsin where x = X - X. The midsurface displacement vector may be expressed by A A A A = ix - jy = na t + ta r The components in the two coordinate systems are related by a = - x sin 0 + y cos 9 and as = x cos 9 + y sin. Substituting Equations (2.18) into these and dividing by a, we get sin [ sin ( )sin sin 5 = " sin 9 -- - sin 9 + cos 9rcos n -, cos ( 9) + * cos - cos 0 (2. lc ~L sin / sin

15 and = cos 0 Csin sin in sin [os 8 sin sin sin + sin cos - cos O. (2.20) Lsin P3 / sin These last two equations are exact displacement expressions for a circular arc shaped configuration with uniform midsurface strain. They form a basis for estimating the relative magnitudes of displacement and displacement gradient quantities. The quantity 1 may be regarded as the parameter for the family of configurations. The counterclockwise rotation is = E = =(1- */13). (2.21) Thus the normal and tangent unit vectors to the deformed midsurface are A* A n = n cos (1-P /P)9 + t sin (l-P1/P) (2.22a) and A* A * at = -n sin (1-P /P)e + t cos (1-* /P)). (2.22b) The midsurface strain is C = a 13 /(a3) - 1 = 1 sin 3/(3 sin 1) - 1. (2.23) To examine the relative magnitude of the various quantities, they were

16 * evaluated for two different values of: (1) 3 =0 corresponding to the dead center position, and (2) P = -f corresponding to the completely snapped-through position. * For 3 = 0, we find, by using L'Hospital's rule whenever appropriate, that = 1 - (O sin P sin O)/P - cos 0 cos, r = (O cos a sin P)/P - sin 0 cos P o) = 9e (2.24) and = (sin P)/P 1 * For 1 = -P P we get = 2 cos E (cos 9 - cos ), = 2 sin 0 (cos 9 - cos P), D = 20? (2.25) and C = 0. 0 Assuming that P is small, we find that in both cases a attains its maximum value at 9 3/-7. After expanding in a power series in P, we find that at e = P/^|T5 1'/\ - 6/2 for D = 0 and'/'| - 2/0P for * = -

17 With: = 0.1 and: = -I, for example, j'/* = -200 and a can be neglected in comparison with' with an error of about one half of one percent. All the approximations of the previous section were examined for the two configurations in a similar manner. In all cases the error of the approximations is ( 2). This indicates that the analysis subsequent to Section (2.2) is valid 2 to order P. Thus the analysis is restricted to semi-opening angles less than about 0. 1 to insure that the errors are within a few percent. The assumption of elastic material behavior may also place a further restriction on P. For example, with mild steel the strain must remain less than about 2 x 10. If the shell attains the dead center configuration, the midsurface strain given by Equation (2.24.) is approximately 2 c = P2/6 0 For the steel shell P must be less than about.11 radian for no midsurface yield. The actual dynamic response will not necessarily bring the shell into the dead center configuration. Thus it may be possible for steel shells with larger semi-opening angles to deflect below the dead center position without plastic deformation via dynamical configurations with smaller midsurface strains. Nevertheless, if P exceeds the value 76YE, where Y is the yield stress, then the strain should be examined for possible yielding. 2.4. THE NONLINEAR PARTIAL DIFFERENTIAL EQUATIONS OF MOTION In this section the equations of motion and the natural boundary conditions are obtained from Hamilton's principle. The kinetic energy per unit length along a generator is

18 P h/2 -* -* T - * Rt (a-z)dz de 2 -3 -h/2 6t 6t If the rotary inertia is neglected, we can replace R with the corresponding _ * __* midsurface position vector, R = CO. After introducing the dimensionless variables we obtain a nondimensional kinetic energy 2 P T - E(ah 1 - [i2 *2 (2.26) Eah -2 - in which the dot denotes differentiation with respect to T, where r - ct/a (2.27) and 2 E c2.- (2.28) p(l-v2) The strain energy per unit length of cylinder is P h/2 V - h/2 to(a-z)dz dO. (2.29) When the radial stress is neglected, Hooke's law reduces to = cE/(l-v2). (2.15) Substituting Equations (2.15) and (2.11) into (2,29), carrying out the integration and dividing by Eah/(l-v ) gives a nondimensional strain energy. It is V (l-v ) 1 2.) V ____ 1-ta (222) dO (2.o0) E V- -~ (2 2~

19 where a = h/(ai ) *. (2.51) The external loading is shown in Figure 4. p(, t) e e Mg~ i~ ~~e ~(",t) NQ (-(ttt) e Nj-3,t) e Figure 4. Cross section with external loads. The external bending moment, tension and radial shear force resultants per unit length along a generator are M, N and Q, respectively. The inward e e e surface pressure is p. We shall take the direction of the differential pressure force acting on an element of the shell to be the same as that of the current midsurface normal; its magnitude is considered as p multiplied by the undeformed midsurface differential area. These approximations are consistant with approximations already introduced. It is convenient to introduce the following dimensionless external loads: = (l-v2)N /(Eh) (2. 2a) e e M = 12a(l-v )M /(Eh) (2. 2b) e e bp = a(l-v )p/(Eh) (2. 2c) e = (l-v2)ge/(Eh). 2d)

20 The equations of motion can be obtained through Hamilton's principle. They are 2T - - -'- (M +M) - N- (N_) p (2. S) and _- 9 1= - f p.(2.34) (N-N )b = 0 (2.35) (N i'C M+Q )~s = 0 (2.36) and (Me-M)b[L' ] = 0. (2.37) The quantities in the parentheses always vanish at the ends. They represent force boundary conditions when there are no restrictions on the variations. When a variation quantity is required to vanish at the end, it replaces the corresponding vanishing parenthesis as the boundary condition; there the external loading is an unprescribed constraint force. Henceforth we shall restrict ourselves to boundary conditions where the ends are fixed against translation and have elastic restraints against rotation. The fixed end conditions imply that; and r are zero at e = ~P. The restraints against end rotations considered here have the form given by

21'K2/ -M = O at a = (2. 8a) and'K1/ + M = 0 at e = -P (2. 8b) where K1 and K2 are elastic constants. The factor 1/P has been included for later convenience. To complete the statement of the forced vibration problem, initial conditions must be specified on ~ and 4, 2.5. CONDITIONS JUST AFTER AN IMPULSIVE PRESSURE LOADING In this dissertation we take the external surface loading as an impulsive pressure applied to the shell at T = 0. If the duration of the impulse is e, then a dimensionless impulse is 1(0) = f pe (,T)dr. (2.59) 0 e The geometrical configuration of the shell changes very little during the impulse for sufficiently small ~. Thus we take (0e,~) = 5(9,0) (2.40a) and E (ua ) = ( (9,0) a (2. 40b) Equations (2 33) and (2.34) can be integrated to obtain an impulse-momentum

22 relation. Since the force quantities which depend upon deformation alone are always finite, they will make a negligible contribution for sufficiently small ~. Discarding these terms and taking the geometry to be constant during the impulse gives ~(a,~) = S(e,0) + I(e) (2.41a) and;(e,e) = 4r(e,o) -'(e,o0)I(). (2. 41b) In most problems of interest, the shell is undeformed and motionless prior to impulsive loading. Letting E approach zero, we take the conditions just after the impulse to be at T = 0. Thus the initial conditions at the onset of free vibration are (e9,0) = 0 (2.42a) (e9,0) = 1(e) (2.42b) (e,0o) = 0 (2.43a) v(e,o) =. (2.43b) For these initial conditions the tangential velocity component V may be * * neglected in comparison to ~. If it is neglected in the kinetic energy expression (2.26), then the tangential inertial term does not appear in Equation (2.34). Whereas Equations (2.33) and (2.34) constitute a fourth-order system * (15) This approximation was introduced by Reissner. Recently it was examined in some detail for the linear vibration of a cylindrical shell of finite length of Lovell. (16)

25 with respect to time and hence require four initial conditions, by neglecting tangential inertia we reduce the system to second-order in time. With this we may no longer impose Equations (2.43); * and { must assume whatever values are necessary so that N remains spatially constant.

CHAPTER III SUFFICIENCY CONDITIONS FOR DYNAMIC STABILITY 3.1 GENERAL COMMENTS In a recent paper Hsu considers dynamic stability "in the large" for free motion problems in terms of trajectories in functional configuration space. For dissipative systems he discusses asymptotic stability regions; within each region all motions asymptotically approach a statically stable equilibrium configuration. Hsu presents sufficiency conditions for stability. He considers the motion to be stable if the energy imparted to the structure is less than the potential energy associated with the first critical or static equilibrium configuration encountered by an expanding level energy surface in functional configuration space. For a conservative system Hsu tacitly considers orbital stability and points out that the above sufficiency condition applies there also. Hsu asserts that the sufficiency condition is also necessary for conservative problems. From continuity considerations it is evident that the trajectories within Hsu's asymptotic stability regions are also orbitally stable. Hsu illustrates his approach with the sinusoidal (6) arch problem using curved beam equations derived by Fung and Kaplan. Several of Hsu's ideas are used in the present study in which we treat the circular arch using thin shell theory. In this chapter we develop conditions which are sufficient to insure orbital stability. For our analysis the critical configurations and their strain potential are determined as a function of geometry. 24

25 3.2 THE CRITICAL CONFIGURATIONS The critical configurations are deformation fields which satisfy the equations of motion and are independent of timee For the shallow cylindrical shell the differential equations governing the critical configurations are obtained by setting the time derivatives to zero in Equations (2.33) and (2.34) For no radial pressure we have a2 ( 2"T)v + 2f" + ) N' (1 + ) 1) and - N' = 0 (3.2) From Equation (5.2) we conclude that N is a constant for each critical configuration. Thus Equations (3.1) and (3.2) imply that'+v + 2 A t + N/ (33) where A = 1 - N/(2a2). (3.4) The general solution to Equation (3.3) may be written 2 B sinqF + B cosqF + B sinrF + B4 cosrF + N/a (535) 1 2 1 4 In the footnote on page 22, we discussed an approximation in which the tangential inertia term is neglected. Here, however, there is no approximation involved.

26 where q = P + A2 - (3.6a) r = A- - (3.6b) and r = e/ (3.7) For supports fixed against translation we have = 0 at r = + 1 (3.8) t = 0 at r = ~ 1. (359) For the elastic restraints against rotation considered here, Equations (2.38) may be written' K1/P - - " = 0 at 1 = -1 (3.10a) and' K2/p + + 0 = at r - +1 (3.10b) After substituting Equation (3.5) into Equations (3.8) and (3.10), we can express these two sets of boundary conditions in matrix form by s c S C b 1 1 r (511) 2 2 2 2 A -(q s-K2qc) -(q c+K qs) -(r S-K rC) -(r C+K rS) b 0 -(q s-Klqc) (q c+Klqs) -(r S-KlrC) (r C+K rS) b -s c -S C b 1

27 where A 2 b. -aB./N for i = 1,2,3,4 (3.12) 1 i and we have used the following abbreviated notation: s = sin q c cos q (3.13) S = sin r C - cos r A Equations (3.11) may be solved for the b.. The determinant D of the square array in Equation (3.11), after some algebra, may be expressed as -D/(2cC2) = [(q -r )+K2(qt-rT)][T(q2t-K )-t(r2T-Kr)] (3.14) + [(q -r )+Kl(qt=rT)] [T(q t-K2)-t(r2T-K r)] where we have introduced the additional abbreviated notation: t = tan q (3.15) T = tan r. A After some algebra, the bi coefficients may be expressed as A 2 2 cb (q +K2qt)(r +K rT) -D 1 +K 1 2 1 -D~~~~~~~ _ 2 ~~~~~~~~~~(3.16a) 2c2C2 T ~ 2 2 (5.1a) 2c2C2 T -(q2+K qt)(r +K rT) 1 2 D A (r2+K2rT)[t(r2T-KLr) -T( q t-K lq)] 2c2 cb = 2 (3.16b) 22c 2 2 +(r +KlrT)[t(r2T-K2r)-T(q t-K q)] 1 ~~~~22

28 -D Cb3 -(q2+K qt)(r2+K rT) 2 2 -D 1 2c2C t 2,2 (3.16c) +(q +Klqt)(r +K2rT) -D A (q +K2qt) [T(q t-Klq)-t(r T-K1r)] 2cC ~ Cb2 2 2 (3.16d) 2c C 4 +(q2+K qt) [T(q t-Kq)-t(r T-K2r)] The dimensionless midsurface strain (i.e., the dimensionless tension) is related to geometric variables through Equation (2.12). It is N ='- + (')2/2. (5.17) From Equation (3.2) N is independent of 0. Thus Equation (3.17) can be integrated with respect to r to give 1 1 1 1 2 N = r= ~ r=- ~ - t dr) + r (). (3-.8) 2 F=1 ~=-l 4 By imposing the boundary conditions (3.9) on this we get 1 1 1 219) N - f / dr + (') dr (3.19) 2 -1 -1 Equations (3.14), (3.16), and (3.19) define the five coefficients B1,...,B4 and N that appear in Equation (3.5) for 5. 3.3 BALANCED SUPPORT CRITICAL CONFIGURATIONS If the two elastic restraints against rotation are equal, i.e., K =K =K (3.20) then Equations (3.14) and (3.16) can be rewritten as

29 D 4DD2 (3.21) A A Db Db 0 (3.22) 1 3 Db2 = C(r + KrT) (3.23a) 12 A 2 Db4 = -c(q + Kqt) (3.23b) where D = cC [(r q2) + K(rT- qt)] (3.24) and D cC [T(q t + Kq) - t(r T + Kr)] 2 2 cC (t[(q - r )T + Kr] - TKq) (3.25) In Equation (3.11) let us subtract the fourth and second rows from the first and third, respectively; then add the first and fourth rows. After dividing each row by two we get s 0 S 0 b1 0 2 2 A (-q s+K qc) Kqs (-r s+KrC) KDrS b2 0 A=A < 3(5.26) iKqc g(q2c+qs) KrC (r2C+KArS) b A 0 c 0 C b, 1 where K^ = (KI + K2) ( 7) A 2 1 2~~~~~~~~~~~~~~3-7

30 and K - (K - K (3.27) D 2 1 2 For balanced restraints (i.e., K K K) System (3.26) divides into two sets: s S b 0 p p ^A — / (3-28) (q s-Kqc) (r s-KrC) b O and c C b 1. (3.29) (q2c+Kqs) (r C+KrS) b4, 0 O The determinants of the square arrays in Equations (3.28) and (3.29) are the same as D2 given by Equation (3.25) and D given by Equation (3.24), respectively. Two possibilities exist for satisfying Equations (3.22) and (3.23): A A either (1) b1 and b3 are both zero, or (2) D2 vanishes. Equation (3.28) reA A lates b3 to b1 for the second possibility: A b A S$(. b_ 3.4 SIMPLY SUPPORTED SYMMETRIC CRITICAL CONFIGURATIONS A For balanced supports the possibility in which b and 3 both vanish corresponds to symmetric critical configurations. Here, Equation (3.19) serves as the characteristic equation. For the simple support case (i.e., K = K = 0) the symmetric critical 1- 2

31 configurations may be written 2-i~~v 2 ~cosrr 2 q)] FL( 1 cos qr) cosr s = [ILq 2 ( cos q r2 ) q where we have used identities (A.1) and (A.6) from Appendix A. We will see 2 later that q > 1 and that r < 3. Hence an approximation for ~ may be obtained by using the relations 2 ( << 1 (5.32) and r <<. (3.33) Thus a good approximation for ~ is = 12 (1 cos () - 1 )] (3.34) 2cos q 2 The term p (1-r )/2 approximately represents the dead center position. The 2 2 remaining term, P (1-cos qr/cos q)/q, represents deviations from the dead center position; the deviations are smaller in magnitude and have more ripples for larger values of q. To obtain the characteristic equation for the simple support case, we substitute Equation (3.31) into Equation (3.19). After using Definition (A.7) and Identity (A.8) and rearranging, we get the characteristic equation The geometrical parameter X introduced by Equation (A.7) is discussed in Section 14.5 on page 92.

32 G(q) = 0 (3 35) where 2 2 S 4 4 G = 4cC[1 + (q) ](Cs-cq[C+q2 Sr ] - C() [-() ] 2 4cc[cq(S)-sC] + [l-() ]qC2(l sin 2q) -r. [1-(p)4] 1,sin 2r- (3) 2 3 2r - e qxpadin (3.36) By expanding C, S, and sin 2r and using the relations (3.52) and (3-33), we obtain the approximation G ~ sc + q (1-gc 2 (3-37) where g - 2q/3 - 4(q/x) (3.38) By setting the right side of Equation (3.37) to zero we can approximately express the characteristic equation as Jo(q) ='(q) (3539) where ~(q) - qg = q[2q2/3 - 4(q/x) ] (3.40) A careful numerical inves ation varified that this approximation is valid over the entire range of interest.

33 and (qs) = q sec q - tan q. (3.41) Sketches of and, are indicated in Figure 5. Inspecting Equation (3.40) we see that for a given value of q,. monotonically increases with K. Hence the Y curves never intersect for q > 0 and are all bounded by the limiting curve y = 2 q3/3. (3.42) To see if Equation (3.39) is satisfied in the open interval 0 < q < t/2, consider the difference of slopes dj& d_ 2 tan q oYU(q) = = - (3.43) dq dq qcosq We have the following inequalities: sin q < q < tan q for 0 < q < n/2 and sin q cos q < — < 1 for 0 < q < q Hence sin q tan q. 1 < < for 0 < q <-. q cos q q cos q 2 Thus r(q) > 0 for 0 < q <'.

34 20T y =m = 4q5/15 1 1| / (Locus of Maxima of Figure 5i Plots of yand,21f(Locus of 16 Minima of 12-89 Figure Plots Minima ofns. 22.89 Legend: Broken Line -.-. y = Solid Line y = Figure 5. Plots of and ~functions.

35 and it follows that Equation (3.39) has no roots on the interval 0 < q <. 2 The local minimums of occur at (Ni, Ni) where N = 1,23,... Because Scontinously varies with \, there are values of X, kN, for which Equation (.539) is satisfied at these points. They are 4 6 Nit) AN=....(3.44) N 2 ~N (Nit) 3/2 The slopes of at these points are dq (Ni) =2 [1 | 2 (Ni)2 < 0 for N = 1,2,.... (.45) From continuity considerations we can conclude that as X increases special values of k, N', are encountered where the 5 curve just touches a new branch of the j curve and gives rise to a new pair of roots. Because of the Inequality (3.45) we conclude that %N < AN (3.46) A lower bound to the 1 may be obtained as follows. The maximums of? occur at 2 q - *- ( (3.47) 4-0 By using this to eliminate X from Equation (5.40) we obtain the locus of maximum of ~: y~4, -^4 -~ y (3'48)

36 This intersects the line y = q at q = 1.935 >- 2 2 2 From Figure 5 we see that a lower bound for 1 is X such that = n. We get 1 X =6 1 [ 6 2.68 The upper bound to X1 is readily calculated from Equation (3.44). We get;_ ~ 22.895. Hence, 2.68 < 1 < 2.895. (5349) 3.5 SIMPLY SUPPORTED NONSYMMETRIC CRITICAL CONFIGURATIONS For balanced supports (i.e., K = K2 = K) the possibility in which D (given by Equation (3.25)) vanishes corresponds to nonsymmetric critical configurations. Here the vanishing of determinant D2 serves as the characteristic equation; once we have an eigenvalue q, Equations (3.19) and (3.30) and A A Identity (A.6) enable us to calculate bl, b and N. The characteristic equation is

37 ~t TKq (q2-r2)T + Kr (3-50) For completely rigid supports, the characteristic equation is T t = -q -q. (3.51) r For the simple support case of interest here, we set K to zero and get the * roots q = n; n = 1,2,3,... (3.52) The simply supported nonsymmetric critical configurations may be written 2 2 2' {-s2 (b sin qr + 1 cos ) 1 1- ) cos rr + vC ( ( 353) where 4- 1b ~~ 2' (3. 4) As with Equation (5531) we have used Identities (A.1) and (A.6) from Appendix A. Equation (3.53) may be used for the symmetric configurations by setting bl (or B 1) equal to zero. Equations (3.53) and (3.54) are well approximated by The root q = 0 corresponds to no deformation. We also have roots q - -ni (n = 1,2,...) which correspond to configurations whose shape is the mirror image of configurations with q - +nn.

38 2 1 cos qP q (- [b1 sin qr + (1 - c )] q c + 1 ( ( 2) (55) 2 where B1 b1 2s 1T' 1 (3.56) After substituting Equation (3.53) into Equation (3.19) and solving for ~2 b we get 2 [1 + 2T) [ 2 _!< (i) [L (q) ] m r2 q g b -4 I - I n 1 2 r2 sin 2r l T 21 12r 42 r (3 ~ -~ (3.57). (2r)2 _ [1 - ()4) By expanding T, C, and sin 2r and using the relations (3.32) and (3.33) we obtain the approximation 2 b ~ h(q) (3. 58) where h(q) c2 + g(q). (5.59) The function g is given by Equation (3.38). A careful numerical investigation varified that this approximation is valid for the entire range of interest.

59 After examining Equation (3.59) we see that h is negative for sufficently small q and for all q greater than some constant that depends upon X. If h = 0 has no roots, there are no nonsymmetric critical configurations. For each odd pair of two successive roots of h = 0 that are end points of an interval containing integral multiples of t, the multiples of j are eigenvalues q for nonsymmetric critical configurations. Thus we next study the roots of h = 0. The roots of h = 0 may be obtained from the equation cos q = V(q,x) (3.60) where X (q)=_ 1 (361a) 2 q - 4(-) 4 3 X or /= 2( qj (3.61b) In Figure 6 we have sketched the left and right-hand sides of Equation (3.60). An examination of Equation (3.6la) reveals that a is a monotonically decreasing function of K; hence the curves never cross each other. Further, the positive branches are always above the curve y = Ml(q,oo) 2 * (3.62) 2q

40 y-y=2/q,oo) 5 V- y=3/q (Locus of Minima for2y) y=/\\, (=q,X) Y C __ I ^q y=cos2q 1 2 3 4 5 6 7 8 9 q Figure 6. Plots of cos q and (q,.). This latter function has no intersections with the curve 2 y = cos q (3.63) in the interval 0 = q < J (3.64) 2 ( 3 e 64) Hence no members of the family y ='(q) (3.65) intersect the curve (3.63) on the interval (3.64). The second root of= 1 is given by

41 q = q() (1 + (.66) The function q*(X) is a monotonically increasing, unbounded function of k. As X is continuously increased from values of X,, are encountered that satisfy qn (s) = n, n = 1,2,... (.67) When X is increased to X, the n'th nonsymmetric critical configuration first arises. After a little algebra, we find that X is given by n =4 (nit) -/ (3.68) n (nit4 3/2 This is identical to the expression for X, Equation (3.44). Since the expressions for X and X are identical, the nonsymmetric and the even numbered n n symmetric critical configurations have the same eigenvalue, q = nt, for X = x. n It is of interest to note that as X is increased, a pair of symmetric critical configurations always arise slightly before a single nonsymmetric configuration; we can see this from Inequality (3.46). This result agrees with the results of Hsu for the sinusoidal arch using curved beam equations. 3.6 THE POTENTIAL OF THE SIMPLY SUPPORTED CRITICAL CONFIGURATIONS It is convenient to introduce a second dimensionless form of the strain potential as V "~~~~~~~= -^~~~j~~ (5.69)

42 where V is given by Equation (2.30). After imposing Equation (3.2) on Equation (2.30) and substituting V into (3.69) we get ~"~2 1 1 v N2 1 f+ dr. (3.70) V ~2 -1 After substituting Equation (A.7) and Identity (A.8) into V, we have 4 1 V = [1 - ()24 (q + J M2 dr (3571) q 2 -1 The bending part of V is given by VB f M dr (35.72) B -2 where M = -(W" +:). (2.13) Substituting Equation (3.53) into (3.72) we get, after some algebra, 2 4 V _ q 1 ( sin 2 q V - ("B 2- 2 2 1 2q [L+ (f)] + +r + (3573a) 2c 1 2 where 1 r1 sin 2r [3 r 1 c2(r)2 2 2r L2 qj q + - [1 + cos 2r] [1 + (5,75b) 2 q

43 and 3;= 2 1 q 4~ 4- - [ P2] t + 2 2q {4(5)a::- 2 (3.73c) Limiting processes, such as we have already done several times, lead to the approximations W - 1 (3.74a) and rt^~~~~~ —24~ _*~ i.(3.74b) 2 2q Thus a good approximation for VB is 1~ 2 sin 2q 1 t +_ V =-b (1 - ) + - + (3.75) B 2 1 2q 2e 2q The total potential is well approximated by 4 b V () + (1 sin 2q 1 3t V +- (1- ) + +1. 6) X 2 2q 2c - 2q For the symmetric configurations, b is zero and V simplifies to 4 V (a) + + 1. (3.77) s X 2- 2q The approximate expression for the characteristic equation, Equation (3.37), may be rewritten as 2 4,ti-22 -! —.^ * (3.78) ( 2q ) =2E 3 ) _ q Using this in Equation (3.78) we have

44 2 4 V 1 + - ( t (3.79) s 3 X q For nonsymmetric critical configurations, by using Approximation (3.58) and q = nt, we obtain the accurate approximation 2 4 (V )= + (nt)- (n) (3.80) n It is interesting to note that Equation (3.79) gives the same result for q = nrt. Since we already established that when X = X, the eigenvalues are q = nit for both the nonsymmetric and the even numbered symmetric critical configurations, we now see that the potentials are also the same for this X. Plots of V versus X using the exact expressions (3.73) with q also obtained from exact expressions are shown in Figure 7. The values of \ where each nonsymmetric curve begins are given by Equation (3.68). If we eliminate nit between Equations (3.68) and (3.80) we obtain the locus of nonsymmetric critical configuration beginning points. The locus is 4 n V= 5/4 + 72 (1 + 4 (3.81) n Although it is not obvious from Figure 7, the V curves for the even numbered ** symmetric configurations cross this locus. Turning back to the original purpose for the critical configuration analysis, the second lowest V curve constitutes a stability sufficiency *The calculations were also made using the approximate expressions for V and for q; the agreement was excellent. It was not determined if the V curves for the odd numbered symmetric configurations also cross this locus.

/' 60 CW 6b/Legend: / - 1> |2n'th Symmetric Configuration: - - -- // n=5 - 50 - n'th Nonsymmetric Configuration: - - - n (2n-l)'th Symmetric Configuration: P Beginning of n'th Nonsymmetric Configuration at = /=/ 40 (2n'th Symmetric Configuration Curves also Pass Through This Point): + _ Locus of'+' Points (Equation: Y 5/4+A 1+i -367 /72:/ 30 CO.H 20t 4-) -p n=3 0 0 10n=2 c | n=l1 0 1 2 3 4 5 6 7 8 9 10 Geometrical Parameter X Figure 7. Plot of V versus. for both nonsymmetric and symmetric critical configurations.

46 boundary; if the total energy of the shell is less than this level the shell cannot exhibit snap-through instability in the sense discussed in Section 3.1. In Figure 8 we have redrawn the V versus X curves shown in Figure 7 in k4 2 the form 8V(X/t) versus X /4. These last two quantities are directly comparable to ones used by Hsu for the sinusoidal arch. The stability sufficiency curve given by Hsu is also shown in Figure 8. Hsu's curve begins at 2 /4 = 2 (i.e., X 2.828) with a symmetric critical configuration; however as X /4 is increased beyond j (i.e., as A is increased beyond 2.91), his curve is associated with a nonsymmetric critical configuration. For the present analysis, the stability sufficiency curve begins at X = k (see Inequality (3.49)) with the second symmetric critical configuration; as X is increased beyond X = 2.895, the curve corresponds to the lowest nonsymmetric configuration. In terms of V, Hsu's sufficiency curve is below our present curve by a constant value, (1- 2 /5- it /24) 0.25. In Figure 9 we have shown the percentage that Hsu's sufficiency curve is below ours.

47 600 Legend: 2n'th Symmetric Critical Configuration: - - n'th Nonsymmetric Critical Configuration: 500 -(2n-l)'th Symmetric Critical Configuration: - Beginning of n'th Nonsymmetric Configuration at \=. (2n'th Symmetric Configuration Curves Also Pass Through These Points): + 400 - // // |I^~~ ^ —Stability Sufficiency Curve |-1^ ^< sof Present Analysis 300 / 200 Hsu's Stability Sufficiency Curve n=2 100 n= 0 1 2 24 / /4 4 2 Figure 8. Plots of 8V(%/Tr) versus /4.

48 10 9 S 8 H a7 0 o a) cd.4 0, 5 co 2 U -P a) 1 - 0 1 2 3 4 5 6 7 8 9 10 Geometrical Parameter X Figure 9. Comparison of Hsu's stability sufficiency curve with present curve V1N

CHAPTER IV INFINITESIMAL VIBRATIONS 4.1. GENERAL COMMENTS In Chapter II we derived the nonlinear equations governing plane motion of a shallow cylindrical shell. If the motion is sufficiently small, the second-order terms are not significant and the equations governing the motion are linear. Although this thesis is primarily concerned with nonlinear motion, the linear solution is obtained to effect approximate solutions of the nonlinear differential equations. 4. 2. GENERAL SOLUTION OF THE LINEAR EQUATIONS OF MOTION The equations of motion for infinitesimal free vibrations can be obtained by neglecting the nonlinear terms in Equations (2.33) and (2.34) and setting p equal to zero. After neglecting the transverse inertial term the equations become: - a (M +M) - N = 0 (4. la) and N = 0 (4.lb) where M = (~ +~) (4. c) and 49

50 N = r - (. (4.ld) By differentiating (4. la) with respect to e and using (4. b) and (4.1c) we find that 2 v!' = a (v+2 + ). (4.2) We seek a separable solution of the form C = Z(e)T(T). (4.3) Substituting this into (4.2) leads to Zv + 2Z"tZ'(1-p2) = 0 (4.4a) and T + apT = 0 (4. 4b) 2 2 where (-a p ) is the separation constant. The solutions of (4.4a) and (4.4b) are Z = A e9 (4.5a) and T = sin(apT+0) (4 5b) where ~ is a constant and s satisfies s5 + 2s + (l-p2)s = 0. (4.5c)

51 The roots of this last equation are s = O, ~i-T7i, ~\F-,TI. (4.6) A solution to Equation (4.4a) is Z = A1 sinp 9 + A2 cos-p+i 9 + A sin p+ A4 cosh ip-11 + A (4.7) Substituting Equations (4.7) and (4.5b) into (4.3) gives the following solution to Equation (4.2): (,T) = LA1 sin p+ 9 + A2 cos p+ A 9 + A sinhp-1 9 + A, cosh p-l 9 + A5 sin (apT+~). (4.8) It is convenient to introduce the following quantities: X = D/~7 (4.9a) v = v (4.9b) u - \V= vok - (4.9c) w = P = 7-2' (4.9d) u2 v2 P2 t == ap = - = (4.9e) x2 where the signs of the radicals are taken to be positive. The quantities i, u, w, and p may be considered as functions of v. Using Equations (4.9) in (4.8) we have

52 = Z(O)sin(-T++) (4.10a) where Z = A sin vF + A2 cos vF + A3 sinh wF + A4 cosh wF + A. (4.10b) By substituting Equation (4.10a) into Equation (4.1a), solving for N and using 4 = P4/2, we get N(T) = n sin ([T+1) (4.11a) where 2 2 n = A vw (4.11b) 5 Thus N is independent of 0 which is consistent with Equation (4.lb). If we solve Equation (4.1d) for f' and use Equations (4.10) and (4.11) we get'(ST) = [A sin vT + A cos vr + A sinh wI 2 2 + A cosh wr + A (1 - v sin ( IT+) Integrating this yields V(e,9T) = P(e) sin (1rT+~) + F(T) (4.12a) where

53 A A T= - (cos v - cos vF) + - (sin vF + sin v) v v A A - L (cosh wF - cosh w) -~ - (sinh wl' + sinh w) w w + A (1 - v (r+i) (4.12b) For each of a countable infinity of eigenvalues v, the value of v, the ratios of constant coefficients, A:A2:AA 4 A A and the arbitrary function F(T) can be determined from six boundary conditions. For each eigenvalue the magnitudes of both Q and the set of coefficients can be determined from two initial conditions. 4.3. NORMAL MODE SHAPES FOR END SUPPORTS ELASTICALLY RESTRAINED AGAINST ROTATION 4.3.1. The Boundary Conditions For infinitesimal vibrations the boundary conditions for supports elastically restrained against rotation are 1 K2 + + + = 0 at e = (4.1a) 1 at 8; ~B (4.l1b) K - ~ = 0 at 0 =. (4.1b) (D(,T) = 0 and 5(-P,T) = 0 (4.14a) 4(D,T) = 0 and (-P,T) = 0. (4.14b) Imposing the latter of conditions (4.14b) on Equation (4.12a), we see that F(T) = O. (4.15)

54 The first of Equations (4.14b) and Equations (4-15) and (4.12) imply that A sin v A sinh w A f - + + 1 = 0 (4.16) v w vw where / 2 2 f - vw V w. (4.17) Imposing the Boundary Conditions (4.14a) on Equation (4.10a) leads to A1 sin v + A2 cos v + A3 sinh w + Al cosh w + A = 0 (4.18a) and -A1 sin v + A2 cos v - A3 sinh w + A4 cosh w + A5 0. (4.18b) In view of Equations (4.14a), Equations (4.13) imply that 2 2 -A (v sin v - K v cos v) - A2(v cos v + K v sin v) 2 2 +A (w sinh w + K w cosh w) + A4 (w cosh w + K w sinh w) = (4.19a) and 2 2 -A (v sin v - K v cos v) + A (v cos v + K v sin v) 2 2 +A (w sinh w + K w cosh w) - A4(w cosh -w + K w sinh w) = 0.(4.19b) 1this chapter the following abbreviated notation is introduced: In this chapter the following abbreviated notation is introduced:

55 s = sin v, c = cos v, t = tan v S = sinh w, C = cosh w, T = tanh w (4.20) With Equations (4.20) the Boundary Conditions (4o16), (4.18), and (4.19) can be expressed in matrix notation as s c S C 11 A v2 2 2 2 -(v s-K2vc) -(v c+K2vs) (W S+KwC) (W C+K2wS) 0 A 0 = ws 0 vS f A +(vs-Kvc) -(v2c+K vs) -(2S+K wC) (w C+KlwS) 0 A -s C -S C 1 A (4.21) 4.3.2. The Characteristic Equation For Equation (4.21) to be satisfied, the determinant of the square array must vanish. This condition gives the characteristic equation which may be expressed as [-(2u T+K2w)t + K2vT]-[(w3 + K (2u2T-vfl) t + (v3T-u 2fl) - KLTwf] +[-(2u2T+K w)t + K vT]-[(w3 + K2(2u2T-vfl) t + (v3T-2u fl) - K2Twfl] -. (4.22) For brevity we introduce the quantities 2 a. = 2u T + Kiw 5. = w3 + K (2uT- vfl) 1 1 i= KvT - 2u fl - K. Twf 5i = KivT i -= 1 or 2. (4.23)

56 Using these, Equation (4.22) can be expressed as At - 2Bt - C = 0 (4.24) where A 2= 21 + 12 (4.25a) 2B (= 2( 2 + 1 2) - (2Y + 12) (4.25b) C = 271 + 12 (4.25c) Solving Equation (4.24)for t gives t = m(v) (4.26a) and t = (v) (4.26b) where 2 _f^ B + 2 + AC (v) = B - B +C (4.27a) and ~, (B- B 2 + AC7b).(v) A (4.27b) These results simplify for several special cases. When the two elastic restraints are the same, Equations (4.27) reduce to

57 (v) = KvT (4.28a) 2u T + Kw and -(v5T-2u2f ) + KTwf X(v) = 1 (4.28b) w + K(2u2T-vfl) where KK = K2. (4.29) For this case the characteristic Equation (4.22) reduces to the product of two factors [-(2uT+Kw)t + KvT][w3 + K(2u T-vf )Jt + v3T-2u f) - KTwfl] = 0.(4..0) Equation (4.28a) is associated with the vanishing of the first factor whereas Equation (4.28b) is associated with the vanishing of the second. We shall see later that Equation (4.28a) leads to asymmetric modes while Equation (4.28b) gives rise to symmetric modes. Thus in the more general case the set of roots given by Equation (4.26a) can be considered to give rise to "asymmetric like" modes, i.e., modes that become asymmetric when K1 approaches K2; similarly, the set of roots given by Equation (4.26b) gives rise to "symmetric like" modes, i.e., modes that become symmetric as K1 approaches K2. For the case of rigid end supports, the equations for the two sets of roots are obtained by dividing the numerator and the denominator of the right

58 side of Equations (4.28) by K and letting 1/K approach zero. They are vT (v) = T (4.31a) w and Twf (v) = 2 -l (4.31b) 2u T - vf For the simple support case, we can set K to zero in Equation (4.28) and obtain f(v) = O (4.32) and -v3T + 2u2f 4(v) = 1.(4.33a) w Equation (4.32) implies that v = nT; n = 0,1,2,... (4.33b) The zero root corresponds to a trivial solution. In numerical work that follows, attention is focused on the simple support case. To obtain roots associated with Equation (4.33a) we need to examine the We will be able to see this from Equations (4.41a), (4.42), and (4.43) and the identity v2 + w2 = 2u2.

59 behavior of T. Figure 10 shows plots of 7 versus v for K = 3, 7, and 15 and P =.25. Although difficult to detect in Figure 10, there are two roots near the origin: v =, i.e., u = 0 (4.34a) and v = e- i.e., w = 0. (4. 4b) These two roots lead to trivial solutions. The function C is essentially independent of the semi-opening angle P and hence X is the only significant geometrical parameter. This parameter will be discussed in Section 4.5. In general, the function r increases to a maximum at v approximately equal to.688X. It then decreases very sharply. As can be seen from Figure 10, the roots of Equation (4.26b) will be nearly odd multiples of </2 except for one or two roots in the vicinity of v = X or v = j/2. In general we will denote the nontrivial roots of Equation (4.26a) as v and the nontrivial roots of Equation (4.26b) as v. For the simple support case, it is clear from the above discussion that (2n-l) < v < (2n+l) 2 (4.35a) 2 n 2 (2n-l) - < v < (2n+l) n = 1,2... (4.55b) 2 n 2 ( This is also probably true for the general case when K1 and K2 are not equal to zero.

60 7t'* n = 15 h = ~ II s- k = 6 y 0 0 T/2 i 3T/2 2r 55</2 Figure 10. Plots of the characteristic function I for symmetric modes for 5 =.25 and. = 3, 6, and 15 for the simple support case; also, plots of the tangent function.

61 4.35.5 The Linear Radial Displacement Mode Shapes For each eigenvalue obtained from either Equation (4.26a) or (4.26b), the ratio of coefficients A:A:A3:A4 AA can be determined through Equation (4.21). When an eigenvalue and the corresponding ratios of coefficients are introduced into Equation (4.10b), we obtain the linear radial displacement mode shape as Z = A sin v f + A cos v -i A sinh w r + A4n cosh w r + A5 n in n 2n n 5n n 4n n 5n (4.36) where v represents either v or v. We will denote the mode corresponding to n n n v as Z and the mode corresponding to v as Z. A similar notation will be n n n n used for the constants A in For the special case where K = K = K we will find that Z is a sym-t p c w 1 K 2' n metric function while Z is an asymmetric function, i.e., A = A = (437a) ln 3n and A2 A = A O, n = 1,2,., (4.37b) 2n 4n 5n Also we will find that A = 0, n = 1,2,... (4.38) 3n Anticipating these results, we solve for A2 5 A, A, and A in terms of Y 2n' 5nn 4n' 5n Al but solve for Aln, An, A4, and A5 in terms of A2 For the tilde set lnwe getl2n we get

62 A T - [(w K(+vK2t) - (w +wK2T)(v +vKlt)] (4.39a) A2 2 s [(w2 wK1T) (v +vKt) (w+K2 T(v +vKlt) (39b) 2 4 c r [2uTt - K(vT-wt)](v2+vK2t (49c) C L+[2u Tt - K(vT-wt)' (v +vKl A 2 2 5 p2u2 t K(wT+vt) -2Tt+ (VTwt) (4.39d) — _ c2 AS "[[2u + K (wT+vt) [-2u Tt + K2(vT-wt)] 2 2 where [2uTt K2(VT-wt)] (w2wT) (4.59e) L[2u2Tt - K(vT-wt)](w wK2T) For the barred set we get A2 T 2u Tt +K(wt-vT)]'(w2iwK T) (4.40a) AL =~ -[2u Tt + K (wt-vT) (w2+wKT A1 1 A-A sL (2TwKT) (+vK (4+.40b) Ar = s L(.W 2T).(vK

635 F-2 2 4 c [2u Tt + K2(wt-vT)].(v +tvK t) _~~~~~~~'=~~~~~~~ ~ ^ 2(4.40c) A1 VS +[2u2Tt + K (wt-vT)] (v2t+vK2t) and A5_ -c [2u + K2(wT+vt)].[(2u2Tt + Kl(wt-vT) U - = "T 2? (4.40d) A1,T +[2u + K1(wT+vt)] [(2u Tt + K2(wt-vT)] where - = (v +vKt)(w +wK1T) + (v +vKlt) (w +wKT). (4.40e) For the special case where K1 = K = K, Equations (4.39) and (4.40) simplify to A A 1 3 A 0 (4.41a) A A 2 2 A 2 A4 c(v +Kvt ) (4.41b) A C(w +KwT) A 2 5 -c[2u + K(wT+vt)] %~ O 2 (4.4lc) A2 (w + KwT) and

64 A2 A A A -=- - = = 0. (4.42) A1 A A1 A Use was made of the characteristic equation in establishing that A2/A = 0 = A4/A1 in Equation (4.42). From Equations (4041a) and (4.42) it is now evident that Z is symmetric while Z is asymmetric for the special case of K K = K. For this reason let us term Z and Z as "symmetric like" and "asym2 n n metric like" mode shapes for the more general case. For the case of simple supports, i.e., for K = K2 = 0 we get further simplifications: A 2 4 cv /, i _ 4-:= - (4.43a) 2 A2 Cu and A 2 A_ -c2u A5_^~~~ c-~~~u~ -^. ~(4.43b) 2 A w These equations are used in the numerical work to follow. Finally we select the coefficients A2 and A1 so as to normalize the mode shapes, i.e., so that 1 2 f1 Z dA = 1 (4.44) where r = e/:. By substituting Equation (4.10b) into Equation (4.44) and in

65 tegrating, we get 1 = <l- + + + A -s A1 sin 2 21 + sin 2v + A2 Fsinh 2w 3 2v 2 2Uv 2w 2 sinh 2w A22 + A A 2 4 2w A1A3 2 Lw s L u - wSc+vCs 4s 4S + AA4 2 [w + -AA + A5. (4.45) U 22v 4+ A Equation (4.45) can be used to obtain A2 and A1 for normalizing Z and Z, respectively For the special case where K = K2 = K, we find by using Equations (4.45), (4.41a), and (4.42) that / 2 F 0 j v \ 2 v- 1 1sin 2v 4 sinh 2w o 1 + + 2v 2w wSc+vCs 4s A 4S -- 2- - - ^ v v (4.46a) 2 v w L u /A A and 1 sin 2v -2 = -2v s v = v (4.46b) A 1 For simple supports, i.e., K = K2 = 0 Equation (4.46b) gives A1 = 1. (4.47) Equation (4.36) gives the normalized inward radial mode shapes when v n

66 are given by either Equations (4.26a) or (4.26b). The quantities A2 and A 2n In satisfy the normalizing condition, Equation (4.45), with the ratios of coefficients given by Equations (4.39) and (4.40). In Appendix B on page 131, it is proved that these normal modes are orthogonal to each other. Thus we have the useful orthogonality relations 1 Z Z dr =; p,q = 1,2,... (4.48) 1 p q pq where Z represents either a "symmetric like" or "asymmetric like" mode shape and 5 is the Kronecker delta. Figures 11 through 16 show plots of pq Z (e)/Z (0) versus r = e/P with X as the parameter for the simple support case for P =.05 for the first six symmetric modes. The corresponding plots for the asymmetric modes are ordinary sinusoids. 4.3.4. The Tangential Displacement, Tension and Bending Moment When the eigenvalues and amplitude ratios are introduced into Equation (4.12b), we obtain the tangential displacement mode shapes. Using Equation (4.16), Equation (4.12b) becomes A A Alm A2m /( =r) (cos v - cos y ) + -- sin v r m v m m v m m m + - (cosh w F - cosh w ) + - sinh w r w m m w m m m v'vI + A m Fmr (4.49) where vn again represents either vn or vn. Accurate plots were also made for f =.25 but no differences could be detected by eye.

67 1.0 ~.5 -.5 -1.5 r = e/p Figure 11. Plots of first symmetric mode shapes normalized to the center, l(e)/iZl(O), versus F = O/n for O =.05 with the shape parameter A. as the parameter.

68 1.5 \z7 1. k- 3.5oJ -.5 -1.- N =10 -1.5....... 0.2.4 c6.8 1.0 r = e/p Figure 12. Plots of second symmetric mode shapes normalized to the center, Z2(O)/Z2(0), versus r = ~/P for P =.05 with the shape parameter X as the parameter.

69 1.5 1 —- =10 1. oN \ \\ = 155 -l.0-1 00.2.4 6.8 1.0 nr - e/: Figure 13. Plots of third symmetric mode shapes normalized to the center, Z3(0)/Z3(O), versus r = e/B for --.05 with the shape parameter A as the parameter.

70 2 4 61.0 1.5 0 0 -1.5 0.2.4.6.8 1.0 r = e/P Figure 14. Plots of fourth symmetric mode shapes normalized to the center, Z4(e)/4(0), versus r =e/) for 3 =.05 with the shape parameter as the parameter.

71 1 5 J —- 0 J | l - cD -1.5 0.2.4.6.8 1.0 Figure 15. Plots of fifth symmetric mode shapes normalized to the center, Z5(e)/5(0), versus r = e/3 for B =.05 with the shape parameter X as the parameter.

72 1.5 1.O \o0 x = 15 1 I I y.hr = 10 0.2.4.6.8 1.0 Figure 16. Plots of sixth symmetric mode shapes normalized to the center, Z 6()/Z6(0), versus r = e/p for 3 =.05 withthe shape parameter \ as the parameter,

75 Unlike the inward radial displacement mode shapes Z, the tangential mode shapes are neither orthogonal nor normalized. Figures 17 through 20 show plots of T /A2 versus r = Q/3 with K as the parameter for the simple support case for 3 =.05 for the first four symmetric modes. The constant tension associated with these mode shapes can be obtained from Equation (4.11b). We have 2 2 v w n n n A= -AA (4.50) n 5n 4' Finally, we introduce mode shapes for the dimensionless moment resultants: m. = "(Z + Z ) (4.51) n n n In Equations (4.50) and (4.51) we will again denote these quantities with a bar when referring to asymmetric like modes and a tilde when referring to symmetric like modes. 4,4. THE INFINITESIMAL MOTION FOR SUPPORTED ENDS In the previous section we found a countable infinity of modal solutions to the equations of motion, Equations (4.1). The general solution is a linear combination of these solutions. Thus = E [T Z Z (4.52) n=l n n n n * Accurate plots were also made for: =.25 but no differences could be detected by eye.

74 5 -3.1.5 - -.1 /10 l-k=15 -.0.2.4 6.8 1.0 r = e/B Figure 17. Plots of l/iA21 versus r = e/P with. as the parameter for the simple support case for the first symmetric mode and p =.05.

75. e3 C" x=l0 - 0 2..6 v.8 1. 0.2.4.6.8 1.0 r =e/ Figure 18. Plots of W2/A22 versus r = 9e/ with X as the parameter for the simple support case for the second symmetric mode and =-.05.

76.4 3 - = 10.1 -,2 -.3 -.4 ~0.2.4 6.8 1.0 rY - e/ Figure 19. Plots of {3/A23 versus F = e/: with X as the parameter for the simple support case for the third symmetric mode and.05.

77 = \-k=7.1 -.3 0.2 4 6.8 1.0 r = e/B Figure 20. Plots of Y4/AL24 versus r = 9/ with X as the parameter for the simple support case for the fourth symmetric mode and a =.05.

78 For brevity, we write g = Z T Z (4. 53) n n where T = T sin (l T + 0 ). (4.54) n n n n The summation will sometimes denote the right side of Equation (4.52) or may also denote, T Z alone or Z T Z alone, depending on the context. The quannn nn * tities T, 5 nand n are constants, Similarly, we have n n?n = T T (4.55a) n n N = ~ T n (4.55b) n n M = Z T m. (4.55c) n n From Equation (4.53) we have at T = 0 = m T )Z (4.56a) and g= ZT (O)z. (4.56b) m m Using Equations (4.56) and the orthogonality condition (4.48), T (O) = T sin = f1 N(,0)Z dF. (4.57a) m m mand and

79 T (O) = z T cos 0 = 1 (eo)z mdr. (4.57b) m mm m -.1 m When the loading is a pressure impulsively applied to a previously motionless and undeformed shell, the initial conditions are given by Equations (2.42)e Introducing these into Equation (4.57) gives =0 (4.58a) m and T = i IZ d. (4.58b) m m 1 m where I is the nondimensional impulse. Equations (4.53), (4.54), (4,55), and (4.58) completely define the infinitesimal motion for the type of initial conditions considered here. As a specific example, numerical calculations were made for the dimensionless impulse distribution I(e) (9,0) = ^ (eo) = Zos (4.59) where v is the initial velocity of the shell at 9 = 0. Figures 21 through 26 show plots of cN/v (= -cCAVG /v) versus T; Figures 27 through 52 show plots of ci(O,T)/v versus T. For each series of plots, the semi-opening angle: is O05 and the X values are 3, 7, 8.25, 10, 13.35 and 15. The infinite series were truncated after the first ten symmetric modes for these plots. The dimensionless midsurface strain (i.e., the dimensionless tension) shows a beating effect for certain values of X. The effect is quite pronounced

1.0 0 ^o.I ii i f "".2!> 0' 0 C)H 2 cr.2 U) 0 4o 80 120 16o 200 24c) Dimensionless Time T Figu~re 21. Plot of dimensionless midsurface strain Nc/ o( rAG/o vru ouin -^ntsm vbai f. a (e oe e rtnAe ri- i -1.0 —--- - ------------— I 0 140 80 120 i60 200 2140 Dimensionless Time T Figure 21. Plot of dimensionless midsurface strain Nc/v0 (= -tAVGC/vo) versus T during infinitesimal vibration for B =.05 and %. = 5 (ten modes were retained).

C> o I I 1.0 0 -10,, 0 4o 80 120 160 200 240 Dimensionless Time T Figure 22. Plot of dimensionless midsurface strain Nc/vo (= -~AVGC/vo) versus X during infinitesimal vibration for ~ =.05 and X = 7 (ten modes were retained). -rp c 0 4 Dimensionless Time ~ Figure 22. Plot of dimensionless midsurface strain Nc/v0 (= -rAVGC/vo) versus T during infinitesimal vibration for P =.05 and \ = 7 (ten modes were retained).

0.6 0 C)'~.2 U)0 C2 2ii2 co.c.'6 I,o|1' 0 -1. 01 Dimensionless Time T Figure 215. Plot of dimensionless midsurface strain Nc/v0 ( C/v) versus T du ingr| J I' w during ini vibra 1 or | in =82 t me re a d *r-4 0 40 80 120 l60 200 240 Bimnzensi-onless Time T Fig;ure 25~ Plot of dimensionless midsurface strain Nc/vo (= -^r'9/^ ~) versus T during infinitesimal vibration for P.05 and \ = 8.28 (ten mnodes were retained) o

o 1.0 | I 0) ( I II 0I~E~~V 2 00 | C) -1. I 0 Dimensionless Time T Figure 24. Plot of dimensionless midsurface strain Nc/vo (= -tAVGC/Vo) versus during infinitesimal vibration for f =.05 and X = 10 (ten modes were retained)~ W.l ed -l O —-------— t —--- ---— * —-— i —--- 0l 408 2 6 0 4 o ~ ~ ~ ~ ~ ~~~imninesTm

0 o 6.6 1. 10 - - 4C3 -1.0 U -) =C/ 0 40 80 120 160 200 240 Dimensionless Time T Figure 25. Plot of dimensionless midsurface strain Nc/vo (= -AVGc/vO) versus T during infinitesimal vibration for 3 =.05 and X = 135.5 (ten modes were retained).

1.0 0 0 f.6 Cdd II A 0:o J ( I t 11 h. h I~~~~~~~~~~~~~~~~~~~~~~~ cto <u ~ i i i' i i I i'1 0d i I! 1I t I' 1 { I In ^ i j 1 s. * 1 t I i i I:: l CH - 2 I I ro 0 *H -1.0 i 0 40 80 120 16o 200 240 Dimensionless Time T Figure 26. Plot of dimensionless midsurface strain Nc/v (= -tAVGC/vo) versus T during infinitesimal vibration for 5 =.05 and \ = 15 (ten modes were retained).

2,0 o 1.2 H O 4 0 u 1.2 -r2 Dimensionless Time T Figure 27. Plot of a dimensionless central deflection [(O,)c/v versus T during Hinfinitesimal vibration for I = 05 and \ = 3 (ten modes were retained). infinitesimal vibration for 8 =.05 and X = 3 (ten modes were retained).

2.0 1" L, l,, I I I,, 0 1.2 0 0 -. Q) 0 4,0 fl 0 - 4-, \ -r-1.2 -2.0 0 40 80 120 160 200 240 Dimensionless Time r Figure 28. Plot of a dimensionless central deflection r(O,^)c/vo versus T during infinitesimal vibration for 5 =.05 and X = 7 (ten modes were retained).

0 3.0 0 0 0 \ -H CQ 0 -p o G 0 1 0 CQ4 V -5.0 0 40 80 120 l6o 200 24O Dimensionless Time T Figure 29. Plot of a dimensionless central deflection (OT)c/vo versus T during infinitesimal vibration for P =.05 and \ = 8.28 (ten modes were retained).

4 C 0 -p O 0 H 0 no.rq H cd 0 C D A 0 -2U'rA 0; 4o 80 12 ii 20 240 Dimensionless Time Tr Figure 50. Plot of dimensionless central deflection *(Om)c/v0 versus T during infinitesimal vibration for D =.05 and %. = 10 (ten modes were retained).; ll 0 40 80 120 l60 200 240 Dimnensionless Ti~me T Figure 50. Plot of dimensionless central deflection C(0^r)c/vo versus T during infinitesimal vibration for B =.05 and \ = 10 (ten modes were retained).

5 0 o o 1 qH 0 /qo -p o O O -3 -5 0 40 80 120 160 200 240 Dimensionless Time T Figure 31. Plot of a dimensionless central deflection t(0, T)c/v0 versus T during infinitesimal vibration for 3 =,05 and X = 135.55 (ten modes were retained).

0 0o 2 a 0 A I o ~ —---- a_ H G3 GH a -2. 0 40 80 120 160 200 240 0.* 0 40 80 120 160 200 240 Dimensionless Time T Figure 32. Plot of a dimensionless central deflection t(O,T)c/vo versus T during infinitesimal vibration for 5 =.05 and \ = 15 (ten modes were retained).

92 for X = 8.28 and for X = 15. For small values of X the dimensionless central deflection consists mainly of the fundamental mode; but for higher values of X the presence of higher harmonics is evident. 4 5. THE GEOMETRICAL PARAMETER X Previously we found that the modal solution is virtually independent of the semi-opening angle: and essentially depends only on X. To appreciate the geometrical significance of X let us first note that the rise of the midsurface is H(e) = a[cos 0 - cos. The power series expansion of H = H(O) is 2 P Ho a0 2 - 4: *J In view of Equations (4.9a) and (2.31), this implies that K approximately equals H 48/h. Thus 2 is essentially proportional to the ratio of the midsurface rise to the thickness. The square of the geometrical parameter may be expressed as 2 = k = P7J (4.60) where 1 = 2a: is the undeformed midsurface arc length. Figure 33 shows a plot of X versus D with (1/h) as the parameter. Also indicated are two forbidden regions: (A) 1/h < 10 in which the shell cannot be regarded as thin, and The ratio h/a is sometimes used as a measure of thinness. Since h/a may be expressed as 23h/l, this ratio is always small when h/l is small for the semiopening angles allowed in our analysis.

93 (B) P < 0.1 in which geometric approximations based on a small semi-opening angle are not valid. There is no theoretical upper bound for \. However values of X corresponding to 1/h ratios much larger than, say, 1000 hold little practical interest. 15K - /h = 1000 -\ / ( / / ^Region A: l/h < 10 10- 1\ /h = 500/ Thin shell theory is not valid. Region B: e >.10 -~/ /~ / ~/ Geometrical approximations 5D / rl/h = 100 based on small f are not ~- ~/ ~/ \ ovalid..05 P.10 Figure 33o Plot of. versus P with 1/h as the parameter.

CHAPTER V METHOD OF SOLUTION OF THE NONLINEAR PROBLEM 5.1 GENERAL COMMENTS In Chapter II, we derived the nonlinear partial differential equations governing plane motion of a shallow cylindrical shell. To obtain an exact solution does not appear to be feasible. In this chapter the nonlinear problem is reduced to a system of ordinary differential equations derived from a Galerkin procedure. The method requires that the solution be represented in terms of a series of orthogonal spatial functions which satisfy the boundary conditions. Here we choose the linear vibrational modes for these functions. The coefficients in the series are functions of time and may be viewed as generalized coordinates. The method yields an approximation whose error is orthogonal to the spatial functions. 5.2 REPRESENTATION OF NONLINEAR DISPLACEMENT AND DEFORMATION QUANTITIES For the nonlinear analysis let us represent the radial displacement by a "modal Fourier series" analogous to Equation (4.53). Thus = Z- TZ (5.1) m m where the T are time varying "Fourier" coefficients. The coefficients can be regarded as generalized coordinates. We assume that tangential inertia may also be neglected in the nonlinear analysis. Then from Equation (2.34) with p set to zero, 94.

95 N = 0. (5.2) With this the tangential displacement H may be obtained by solving Equation (2.12) for 4' and integrating with respect to 0. Imposing the latter of Equations (4.14b) and using Equation (5.2) gives 1 r 1 2 4r = (r + 1) N + f [ -')] d (5.5) -l 2 Substituting Equation (5.1) into (5.3) yields 1 1 i = (r ) + 1) N + T [-27(F) - T () (5.4a) m=l m=l n mn where ~(r) = - J Z dr (5.4b) -1 and r -(r) = Z' Z'dr (5.4c) The nondimensional tension may be obtained by imposing the first of Equations (4.14b) on Equation (5.4a) and solving for N. Thus N = -AVG. (5.5a) where 1 00 AG = - 2 dF - f T N (5.5b) -1 m=l and

96 -1 00 N 1 / ( )2 dF T T (5.5c) 4 M rm mn n -1 m,n=l in which N =? (1) - f Z drF (5.5d) -1 and 1 ~ -= (1) = Z'Z'd F (5.5e) mn mn - m n The arrays N and ~ can be calculated directly from the functions 9 (r) m mn m and m (r) once the integrals of Equations (5.4b) and (5.4c) are derived. mn These functions are given in the next section. The dimensionless moment resultant can be obtained by substituting Equation (5.1) into (2.15). The result is 00 M = Z T m (5.6a) m m m=l where m - (z" + Z ). (5.6b) m m Thus once we obtain the T (T), we can calculate 5, M, N, and i using Equations (5.1), (5.6), (5.5), and (5.4), respectively. 5.~5 TE FUNCTIONS? D AND D THE CONSTANT ARRAY N.5 3 THE FUNCTIONS -/m A mn m To express the functions 7 and nd defined by Equations (5.4b) and (5.4c) in a concise manner, it is convenient to introduce the abbreviated notation:

97 c = cos v, s sin v F, C = Cosh w, S =Sinh w r (5.7) n n n n n n n When no confusion will result we will drop subscripts for brevity. Introducing Z from Equation (4.36) into (5.4b) gives 5m as m m AlC A2s A C A4S F n- ( I- -+ + k+ A5+ -l (9.8) ) 2 ( v v w w The array obtained from Equation (5.8) is A s A S N - ( (1) = (2 + + AA ) m m v w 5m m m =l By adding Identity (4.16) to this, we get 2 2 =V w N.m A.. (5.9) m X4 5m (9 Thus N reduces to the convenient expression for n given by Equation (4.50); m m 00 this is to be expected since Z T N is the linear part of N. m=l m m Introducing the mode shapes into Equation (5.4c) gives 4 (F) = B. I..(m,n,fr) B (5.10a) -rru-i.m ij jn i, j =l where B l v A wA B- vA B - w A = w A (5.10b) ln n ln 2n n 2n' 3n n 2n 4n n 4n and Jc dc dF fc s dr Jc c dr fc S dr m n mn mn mn js c dr fsJ s'' fs c d. s f s s d,,'m n m n mmn I..(m,n,r)= m n (5.1Oc) n J J'C c dr' cc s dr c c dr Jc S dr mn mn mn m n fs cdr d r /S s d S C dr fSSdF fmcdF mn mn imnn

98 The limits are understood to be from -1 to F. The functions I. (m,n,r) can j13 be obtained by direct integration. For m + n the matrix for I. (m,n,I) is v s c -v c s I v s s +v c C r w s c +'v Cs w C c +v s S F mmn nmn n mn m n n m n n m mn n n m m n 22 22 22 2 2 2 v-v v -v v+w v +w m n - n n -1 m n -1 m n -1 vs c -vs c r w s S -vc c ws C -vc S r I1(nm,) nmn mn n m n nmmn n mn m m n 122 2 22 22 v -v v +t-w v +w m m n m n -1 l(nm I(nmr),w S C — w S C w -w S S w -w w -w m n -1 n m -1 w C S -wS C I(n,m,r) 124 (n,m,r) I4(n,m,r) n n 2 2 w -w m n -1 (5.la) For m - n the matrix for I i(m,n,^) is ij Fr+ sin 2vFr1 sin vr wSc+vCs wCc+vs [2 4v 2v v2+w2 v+w rF sin 2vr wsS-vcC wsC-vcS IF_ -d 12 L2 4v J v+w 1 v+ 1 Sinh 2wF + ]' Sinh wI 15 235 4w -1 2 I14 I2 I Sinh2wr r r 1r4 X I24 I4 L 4w 2 i-1 (5.11b)

99 5.4 DERIVATION OF THE EQUATIONS GOVERNING THE NONLINEAR MOTION OF THE "FOURIER COEFFICIENTS" After neglecting the transverse velocity in the kinetic energy expression, Equation (2.26), we get T = /5 2de. (5.12a) 2 Since N is independent of 9 we can replace Equation (2.30) with 2 2 V = p + 2. f Md8. (5.12b) The Lagrangian equations of motion may be written q q where the T are generalized coordinates. Substituting Equations (5.12) into q. (5-13) we get 1d OT N 2 d. Jd (7 ) dr + 2N a N +2 S M d r -= O. (5.14) -1 q q 1 q Substituting Equations (5.1), (5.4), and (5.5) into this and using the implied summation convention and using Ni = ni, we have i i' T /1 Z Z dr + 2N(n [T ~ + T ) m 1 m q 4 n qn mmq 2 -1 + a T 1 m m d' = 0. (5.15) m m q where the e are defined by Equation (3.5e). Using the orthogonality rela-. qn tion, Equation (4.48), and noting that ~ is symmetric, we can simplify nEquation (15) Equation (5.15) to

100 T + 2N (n +- T ) + T = (516a) q q 2 m mq m mq where 1 = c 1 m m dr. (516b) mq -1 m q For infinitesimal vibrations, Equation (5.16a) reduces to T + T (2n n +' m ) = 0. (5-17) q m mq mq Since this must be equivalent to Equation (4.4b) we can conclude that (i2 ) - 2(n2) for m =q (5.18a) - 2n n for m I q. (5.18b) mn q which can also be verified by direct calculation. Substituting Equations (5.18) into (5.16a) and using Equation (5.5c) we get T + ( T) + = 0 (5.19a) q q:-q where 2Nn + NT ~. (5.19b) q q m mq The quantities q, N, n, N and mq are given by Equations (4.9e), (5.5c), (5.9), (5.5a), and (5.5e). The last term of Equation (5.19a), q' is a nonlinear coupling term which becomes negligible for infinitesimal Xq oscillations. After writing out the coupling term in detail, Equation (5.19a) becomes

101 of 2 1 T + G( T) + - T (n ~ + n ~ + n ~ q q 2 m q mn n mq m nq + - T ~ ~ )T =. (5.19c) 2 p pq mn n The initial conditions for the generalized coordinates are T (0) = 0 (5.20a) and 1 1 T (0) = J (,O)Z dF f= IZ ddr (5.20b) -1 -1 for an initially undeformed shell. The system given by Equations (5.19) and (5.20) governs the nonlinear motion. No closed form solution was found for the highly nonlinear system of equations given by Equation (5.19c)e If the system is truncated to find N unknowns, however, it is readily integrated numerically using the Runge Kuta method. To do this, it is convenient to convert the given system to the system of first order equations. Y -; p = 1,2,...,2N (5.21a) P P where Y - T (5.21b) 2q q Y(2-1) T (5.21c) F2q -. T - (5.21d) F(2-1) 2 q (5=21e) q = 1q2,...

102 Equation (5.21a) can be expressed by the matrix equation Y = F(Y) (5.22) where Y and F are vectors having 2N dimensions whose elements are given by Equations (5.21). Letting y be the vector with elements 0 (y,) = T(O) (5.23a) Yo2q q and (yo)2q = Tq(0), (5.2b) the initial condition can be expressed as Y(O) = y. The solution to the system defined by Equations (5.21) and (5.23) is discussed in the next chapter.

CHAPTER VI THE NONLINEAR MOTION 6.1. GENERAL COMMENTS For the freely oscillating structure with prescribed boundary conditions, the dynamic stability problem is to determine the critical initial velocity and displacement fields associated with instability. A technically important case, one of primary interest in this dissertation, arises from impulsive loading of structures with a prescribed displacement field. Here only the critical initial velocity fields are sought. In Chapter III we obtained a sufficiency condition for stability. Violation of this condition does not, however, imply that snap-through will occur. In this chapter we seek the actual critical initial velocity fields associated with dynamic snap-through instability. We consider instability here as a discontinuity in an appropriate response parameter with respect to changes in initial conditions. To obtain all possible critical fields does not appear feasible. For the numerical work in this chapter (except as noted) we restrict the initial conditions to zero displacements and cosine shaped initial velocity distributions with a slight additional asymmetric component. The initial conditions being considered may be expressed as (9,) = (0) -2 cos( + kf) (6.la) and 105

104 (0o0,) = o (6.lb) where 2 V Tt E = - (6.1c) cP is an initial velocity parameter and f =_ L r (i - r2). (6.id) 2 The function f is an asymmetric cubic having a maximum of unity at r = 1/1 and k is a constant. In most of the subsequent numerical work, the initial conditions are symmetric, i.e., k = 0. To examine the sensitivity of the nonlinear response to slight deviations in shape from the pure cosine initial velocity distribution, numerical results are also presented for initial velocity distributions having the form of (1) the lowest symmetric radial displacement mode shape, (2) a triangle, and (3) a cosine plus a small asymmetric component (k = 0.1). 6.2. MEASURES OF RESPONSE We have considerable latitude in choosing a measure of response. A significant discontinuity in the motion with respect to initial conditions will be attended by a discontinuity of many continuous functions of the elements of the vector Y. Since the elements of Y are dimensionless, we could consider the length of the configuration vector Y, IYI =7/^ YA X f1 r +;1 dr 2,1/2 (6.2) m=l _ -1 -1 _

105 as a measure of response. We recognize the second integral in Equation (6.2) as 2- 2 - T (E - V) where E is the nondimensional total energy. Thus - 1 2 2 - - 1/2 Y I= [lf dr +-(E - )] (6.3) -1 and hence Y depends upon current deformation and not explicitly upon the velocity. The response parameter which is adopted for presenting our results is defined as 1 1 2 f a~ dF if Hdr 2 -l It represents the ratio of the average radial displacement to the average rise of the shell. This measure does not directly reflect the asymmetric component of the current shell configuration. The asymmetric component of [ only affects 6 indirectly through its effect upon the symmetric component of 5. Nevertheless, it is a good measure of the overall shell configuration when the motion is essentially symmetric and has been used by a number of investigators. 6.3. RESPONSE PARAMETERS FOR THE CRITICAL CONFIGURATIONS For interpreting our results, it is of interest to make comparisons between actual configurations encountered and the critical configurations. For simply supported critical configurations we have

106 5 - ) 1 (6.5a) 7(1 q q2 or v + ( 4 (6.5b) q - 1 + - (6.6) 2 3 q 2 where q is either nir for the n'th nonsymmetric critical configuration or is a root of Equation (3.35). The mode shape "Fourier" series expansion of Equation (5334) for symmetric critical configurations is 0o Zt BZ (6.8) n n n=1 where the coefficients are in v sinh w B ~3Ar{ 2 + n A n n 2 v n w n l / n n sin (vn qsinin (v + q) 2a + n_ _n n_ - -- + + +(w sinh w cos q v - q q 2 2 n nn n n w +q n \ sn q 1 + q cosh w sin q) + 2 sin q 1 2- [2v cos v 2 1 n 2 A J q cos q (v - 2) sin v 1- - [ [(w + 2) sinh w - 2w cosh w ]v w n nr r6.9)

107 6.4. THE NONLINEAR RESPONSE To ascertain the nonlinear response, the system of Equations (5.20) was numerically integrated on a digital computer using the Runge Kuta method. Along with each solution point Y, various quantities were calculated including N, Y(0,T), T, V, E, and the response parameter 5 defined by Equation (6.4). As a check on the numerical accuracy, the total energy E was monitored at each time point. It never varied from its first solution point value by more than 0.2% and the variation was usually considerably smaller. The time step size and the number of degrees of freedom were selected by trial and error. It is found that a time step size of one-twentieth of the infinitesimal vibration half period corresponding to the highest linear mode shape retained is sufficiently small. Steps of about this size were used while determining the number of modes required. The number of modes required varies both with the shell geometry and with the size (and possibly the shape) of the impulse. For very small vibrations, the modes are essentially decoupled and the smallest modal expansion that contains most of the energy of the initial conditions is adequate. However, for the large amplitudes considered here this expansion may not have enough terms. For example with k = 6.0, f3 =.05, and c = 1.55 for a cosine impulse, fourand six-mode expansions contain all but.00347% and.00231%, respectively, of the theoretical input energy; yet the time response of 5 for the two expansions begins to differ noticeably during the second oscillation (see Figure 40). Direct numerical integration reveals that the number of modes needed tends to increase as the geometrical parameter X is increased. The number of modes selected for subsequent numerical work was such that using an additional mode

108 changed the maximum value of E, FMAX' by no more than one percent. The time interval used here was from T = 0 to the next zero of 6. The impulse parameter C for this determination was such that b exceeded unity. A two mode representation was found to be inadequate for all geometries considered. The following numbers of symmetric modes give acceptable convergence for E for the first half cycles: four for 3 < 5, five for X < 4, six for \X 7, and eight for X - 10. More modes are probably needed for longer periods of time. Figures 54 through 42 show plots of the time responses of the parameter 5 for geometrical parameters \ of 4, 5, 6, 7, and 10. Except for Figure 36, these curves are for initial conditions given by Equations (6.1). Of these, all are for pure cosine impulses (k = O) except for three in which k = 0.1: one curve each in Figures 34, 37, and 40. With the exception of two curves shown in Figure 38, all plots are for semi-opening angles 5 of.05. Figure 42 illustrates that for a given % and dimensionless impulse 1(E), the response is virtually independent of B for a few oscillations when P is small. As large amplitudes develop, however, variations in the motion appear for different values of P. Moreover, it was found that the numerical solution was completely different for f ~ 0.2. This is probably due more to the geometrical approximations introduced into our analysis than to a physical phenomenon. The qualitative invariance to small P suggests that conditions for immediate dynamic snap-through for one particular value of B are valid for a range of D. On the other hand, if an accurate description of the long term motion is wanted, the correct semi-opening angle must be used. In the present discussion we are seeking to characterize dynamic snapthrough of the shell. As we have previously noted, such behavior should be

2.4 i s 3135) 1.6 / k=0, E=6. 04 and k=0.1, E=6.08 and Sixs _s _.8g \Six Symmetric Symmetric and Six.8 - / \\ Modes Asyetric Modes 0o -.8. II... I I,i,. II...I 0 10 20 50 40 50 60 70 80 90 100 Dimensionless Time T Figure 54. Response curves for X = 4, 5 =.05, and e = 2.1.

2.4 I v is 3.135) 1.6 - =2.22 (E=6.44) / \ 2.1 (E=6.04) o. 8 cCd.8 a) 0 10 20 30 40 50 60 70 80 90 100 Dimensionless Time v Figure 35. Response curves for A = 4 and P =.05 with a six-symmetric mode representation.

2.4 \ Triangular 1.6 X. // \\ /Shape Impulse \ ( i // \\ / \ /~\ /.-Impulse with 2s 0)h |D~~ // \ / hA / he Shape of the....\ \ \ / Lowest Symmetric --- O-.8/ 0 -.8 I I 0 10 20 30 40 50 60 70 80 90 100 Dimensionless Time ~ Figure 36. Response curves for \ = 4, P =.05, and E = 6.44 with a six-mode representation for two different shaped impulses.

2.4 51s(V1s = 3.193) 1.6 / \ _ IN \ / _______=1.85, k=0.1, Six Symmeho /^ N 7 /'\ / tric and Six Asymmetric.8- 2s(V =6.150) Modes (E=11.093) 0) 0 0 - e=1.75, k=O (E=9.825) c e=1.9, k=O (E=11.581) -.8 I.II 0 10 20 30 40 50 60 70 80 90 100 Dimensionless Time T Figure 37. Response curves for X = 5 and B =.05 with a six symmetric and six asymmetric mode representation.

2.4 2 4ls(V=3.214) ~ 5=0.1 and e=1.2 I is is - ~~~1.6 0^ p=.01 and e=1.2 (E=9.36) / 1.6 - / I *T T2S(V2s=6. 465) < o.8 - =.O05 and c=1.2 (E=9.36) P~ 0 0- ^ P=.05 and E=1.1 (E=7.88) 0.8 10 20 30 I 50I I I 80 I90 0 10 20 350 40 50 60 70 80 90 100 Dimensionless Time T Figure 38. Response curves for X = 6 for a six-symmetric mode representation with k = 0.

2.4 2.4 SLs(V l= 3.214) 1.6 k I t- sq12(v, = 6.465) 4c co 0 \ -.8.8 I. I I I,I I I 0 10 20 30 40 50 60 70 80 90 100 Dimensionless Time T Figure 39. Response curve for \ = 6, X =.05, k = 0, and c = 1.4 (E = 12.75) with a sixsymmetric mode representation.

2.4 1 ls(V =5.214) 1.6 / n 6 -- 60 -, k=O and Four Symmetric Modes 2s 6_6..8 k=O.1 and Six Symmetri ^f^^\ \^ /^^y\ / and Six Asymmetric Modes (/ / / (m=15.8); k= and Six Symm c M s k=O and Six Symmetric Modes (E=15.6) -.8 I I I I I I I 0 10 20 30 40 50 60 70 80 90 100 Dimensionless Time r Figure 40. Response curves for A = 6, P =.05, and e = 1.55.

2.4 1.6 ) (Vl=3.225) $.8 0 -.8 I IV K;. To /, 0 10 20 50 40 50 6 0 70 80 90 100 Dimensionless Time T Figure 41. Response curve for A = 7, f =.05, k = 0, and e = 1 (E = 12.3) with a six-symmetric mode representation.

2.4 F- i(V. =,.250) __. ^~~~~s is___________ 1.6 / -~c=1.2 and Nine Symmetric Modes (E=7.59) ^ /k~ I/~ Y~ —~ f(V=6.725) 4-)~~~~~~~~~~~~~~~s c.8 0 0 — H 0 E =0.81 and Eight Symmetric Modes (E=5.57) -.8 I - I I I a a I I I o 10 20 50 40 50 60 70 80 90 100 Dimensionless Time T Figure 42. Response curves for % = 10, P =.05, and k = 0.

118 indicated by a discontinuous jump in our response parameter at some critical initial conditions. From a numerical viewpoint it is difficult to find such discontinuities precisely. As a feasible alternative, we determine the conditions for which the response parameter changes from relatively small amplitudes to relatively large amplitudes for small changes in initial conditions. The interpretation here is somewhat arbitrary. For definiteness we will consider that the shell has snapped-through whenever t exceeds 2 if the motion is sym2s metric or exceeds 1N if the motion is nonsymmetric. The constants 2s and 1N are values of F corresponding to the second lowest symmetric and the lowest nonsymmetric critical configurations, respectively. If the initial energy imparted to the shell is sufficiently high, the shell snaps-through immediately (i.e., without intervening oscillations). This type of behavior is illustrated by Figure 356, Figure 37 with C = 1.85 and k = 0.1, Figure 40 with k = 0.1, and Figure 42 with c = 1.2. A delayed snap-through phenomenon is illustrated by Figure 34 with k = 0.1, Figure 35 with c = 2.22, Figure 38 with c = 1.2, Figure 39, Figure 40 with k = 0 and Figure 41. The modal convergence of & was established only for the first cycle of oscillation; the accuracy of these curves is not known over the full intervals of time shown. The shape of the impulse may have an important effect upon dynamic snapthrough. Adding an asymmetric component seems to enhance the possibility of snap-through. This is illustrated by Figures 34, 37, and 40. From Figure 37 we see that for \ = 5, the shell does not snap-through with E = 11.581 when k = 0, whereas it immediately snaps through with E = 11.095 when k = 0.1. By comparing Figure 35 (E = 6.44) with Figure 36 we note the response to three different

119 shaped symmetric impulses having the same energy. The general character of the responses is somewhat the same. However, snap-through was delayed for the cosine distribution, whereas it was immediate for a triangular distribution and for a distribution in the shape of the lowest symmetric mode. To study immediate dynamic snap-through, the maximum value of E(= MAX) in the time interval T = 0 to the next zero of ~ is plotted versus c with X as the parameter in Figure 43. For x' 2.8, the curves have a nearly vertical part, denoting the onset of snap-through. For definiteness, the critical value has been selected as the value of (=CR ) corresponding to & = 1.5. A plot of CR CR versus. is shown in Figure 44. CR It is difficult to obtain critical conditions for the phenomenon of delayed snap-through. Extensive numerical computation is required, since there is no a priori condition for determining the time interval that should be examined. It is not clear that delayed snap-through will occur for all geometries. The cases of delayed snap-through that were observed, of course, have values of c below the cCR curve of Figure 44. In Figure 45, plots of the invariant form of energy versus X are shown for our stability sufficiency condition as well as for various dynamic snap-through conditions. All cases of delayed dynamic snap-through fall between curves 1 and 2. 6.5. TOPOLOGICAL INTERPRETATION OF RESULTS It is expected that the critical configurations have an important role in determining the type of motion that the shell may exhibit. We have already used the critical configurations for determining a sufficiency condition which insures that dynamic snap-through cannot occur. The numerical investigation

120 2.5 2.0 1.5 co 0 <-XX o ( ~-D UX ^ Cl KC i / c/u 1.0 o 5 0 1 2 3 4 5 Initial Velocity Parameter E Figure 43. Plots of first swing AJ~X versus c with \ as the parameter.

121 300 2.5 0, 2.0 0 o \.5 4-P 1 0 C) 4-, -)\ 0 2 6 8 10 Geometrical Parameter k Figure 44. Plot of the critical initial velocity parameter cCR versus the geometrical parameter X. for immediate dynamic snap — Ge ometric al arameter X through instability.

122 1. Immediate Dynamic Snap-Through with a Cosine 35 Shaped Impulse. 2. Invariant Strain Potential Energy ViN for the Lowest Nonsymmetric Critical Configuration (A Stability Sufficiency Condition) 50 5. Delayed Snap-Through (See Figure /5) 4. Delayed Snap-Through (See Figure 358) 2. Delayed Snap-Through (See Figure 41) CQ / 0 15 10 4 Xp^~~~~ 3 2 0 I I I I I 0 2 4 6 8 10 Geometrical Parameter A Figure 45. Energy versus \ for a stability sufficiency condition and for dynamic snap-through conditions with a cosine impulse.

123 has shown, however, that a wide variety of dynamic behavior may be exhibited. It also appears that the sufficiency condition is quite conservative, at least for some distributions of impulsive loading. The following heuristic argument is an attempt to give a general interpretation to our results. As a motion trajectory is traced out in functional configuration space, a corresponding surface path is traced out on the potential energy surface. The general features of this potential energy surface are determined by the critical configurations. In a generalized sense, stable static equilibrium configurations are the low points of valleys in the surface. Unstable configurations will occur either at the high point of a hill or as a pass or saddle point in the ridges separating different valleys. The undeformed configuration is, of course, a stable equilibrium point located at the bottom of a valley. For impulsive loading of a previously undeformed structure, all surface paths originate from this point. These paths can only lead out of this valley through passes or over ridges. If the tangent to the surface path has just the right direction at the origin (i.e., the right shaped initial conditions), the path may exit the valley directly through a pass with the energy level of the critical configuration which is located at the pass. On the other hand, if the direction of the tangent is held constant (i.e., the shape of the initial conditions is prescribed), the surface path during the early stage of the motion may not come near any critical point. Thus immediate snap-through can occur only if the initial kinetic energy is sufficient for the surface path to reach the top of a ridge. The total energy of this trajectory will be higher than the lowest pass. A surface path with energy insufficient for immediate exit will be a com

124 plicated curve traced out over the valley. In some circumstances it may encounter a lower ridge or even a pass and then exit (i.e., exhibit delayed snapthrough). In the present investigation the static stability of the critical configurations was not examined. For the case of simple supports the lowest symmetric critical configuration is statically stable. It is likely that all other critical configurations are statically unstable. The lowest energy for which a surface path could leave the valley is the energy associated with the lowest passe This pass is above the lowest (nontrivial) critical configuration. To determine the specific critical configuration located at the lowest pass involves a fairly extensive investigation. In the present problem, however, it is probable that this point corresponds to the critical configuration of second lowest energy: the second symmetric critical configuration for \ < X < X1 or the lowest nonsymmetric critical configuration for k1 < X. The rather high energy levels associated with the actual dynamic snap-through found here are probably due to the particular form of initial conditions. We would expect that distributions with large asymnetric components would snapthrough at lower values. Finally we consider the implication of truncating the number of terms in our modal representation. In our numerical work we are attempting to span functional displacement subspace with the N-dimensional Fourier coefficient vector T = (T1,T2,...,TN)

125 We are, in effect, constraining the motion to paths for which all coefficients with index greater than N vanish. It is difficult to say how constrained and unconstrained surface paths of impending snap-through would differ. It is clear, however, that the constrained surface paths cannot come near critical points whose Fourier coefficients with index greater than N make a significant contribution. Thus in numerical work it is necessary to retain enough modes to accurately represent all of the critical configurations. As a check on the number of modes used for % = 5, the Fourier coefficients for the first three symmetric critical configurations were computed. The convergence of this representation is illustrated in the following table: MODAL FOURIER COEFFICIENTS FOR EXPANSION OF EQUATION (3.34) FOR % = 5 AND F =.05 Mode Symmetric Critical Configuration Number First Second Third 1.764597.531247.315983 2 -.671470 -.291530 -.505843 3.313105 x 101.149214 x 101.170316 x 101 4 -.540739 x 102 -.264157 x 10 -.299848 x 10 5.148549 x 10 2.730080 x 10-3.828655 x 10 3 6 -.753254 x 10- -.380660 x 10-3 -.434175 x 10 3 The first six terms were used to compute the strain potential energy. The results agreed with those given by the analytical expression, Equation (5.72), to within 1.16%, 0.170%, and 0.153%, respectively. We conclude that the number of modes used was adequate for representing the critical configurations.

CHAPTER VII CONCLUSIONS The dynamic stability of plane strain motion arising from prescribed initial conditions of a shallow circular cylindrical shell has been analyzed. Critical configurations based upon thin-shell theory have been determined; a stability sufficiency condition has been obtained in the form of the strain potential energy of the appropriate critical configuration. Quantitative differences are found between these results and Hsu's which are based upon a shallow beam theory and apply to the sinusoidal arch. As a preliminary to finding the nonlinear solution, an analytical solution was obtained for infinitesimal vibrations of the shell with general elastic restraint bending conditions; detailed results were presented for the simple support case. Using a Galerkin technique, a system of ordinary differential equations has been derived that governs large amplitude motion for general elastic restraint bending conditions. For the simple support case, the nonlinear response has been examined and conditions of immediate dynamic snap-through instability have been obtained for initial conditions corresponding to cosine impulses. It is found that a delayed snap-through phenomenon sometimes occurs for energy levels insufficient for immediate snap-through. The energy levels for impending, delayed snap-through have been observed to be between that of the lowest nonsymmetric critical configuration and that of impending immediate snapthrough. 126

127 The required number of modes of representation has been examined and found to depend upon both the size of the impulse and the geometrical parameter \(\4 / h ). It is found that previous investigations of the shallow arch did not use enough modes. Numbers of modes adequate for the convergence of a response parameter & during its first half cycle have been presented. More modes are probably needed for longer periods of time. For several examples of delayed snap presented in this investigation, the modal convergence of t was not established for the full period of integration. Several results of the present investigation are found to be virtually independent of small semi-opening angle P. These include the critical configurations, the infinitesimal vibrations, and the nonlinear solution for the first two oscillations of &. However, when p is changed appreciably, important quantitative differences in the nonlinear solutions appear after longer periods of time. This would seem to suggest that investigations of long time, large amplitude motion of other shallo snapping elements may also need geometrical parameters that characterize ratios of rise to span in some way.

APPENDIX A IDENTITIES FOR THE CRITICAL CONFIGURATIONS Here we develop some identities that are useful for studying the critical configurations. If we multiply Equations (3.6a) and (3.6b) together we get the identity qr = 2 ~ (A.1) By squaring Equations (3.6a) and (3.6b) and taking their sum and difference we get 2 2 q + r 2qr A (A.2) and 2 2 2- 2 ~2 q r = 2qr. (A.3) From (A.3) we have 2 22 (q - r ) = (A- 1)(A +1) (A.4) (2qr)2 Solving (A.2) for A and substituting A into (A.4) we get 2 22 (q= - (A - )(q + r) 2qr Solving this for (A - 1), we conclude from Equation (3.4) that 128

129 __.._. q...-r) _ 2 - =A - 1 =(q- ) (A.5) 22 2qr 22 Using (A.1) in (A.5) we conclude that 2 2 (A.6) 2?2 2(-) Let us define a geometrical parameter \ by* h = P'T4. (A.7) Then, 2 ()4 ()()4 (A.8) Solving (A.6) for the midsurface strain N and using (A.8) we have N = -( ) (A.9) From physical considerations it is clear that From physical considerations it is clear that N < O (A.10) for the critical configurations. From (A.5) this implies that A > 1. (A.l1) *The geometrical parameter i is discussed in Section 4.5 on page 92.

130 Using this inequality, we conclude from Equation (3.6) that q and r are both real and that q > r. (A.12)

APPENDIX B ORTHOGONALITY OF THE RADIAL MODES The orthogonality of the linear radial modes is proved in this appendix. Equation (4.la) may be rewritten N = + t + 2 a2 "1 + +) (B.1) By substituting Equation (4.3) followed by (4,4b) into this we have 2 tv 2 2 N = T(-T)oa[Z + 2Z" + (1 - p /a)Z] (B.2) In view of Equation (4.lb), this can be rewritten N = r T(T) (B.3) where ~ is a constant. By using Equation (4.3) in Equation (4.1c), we can write M = m(9) T(T) (B.4) where m = -(Z" + Z) (B.5) Solving Equation (4.1d) for 4' and using Equations (4.3) and(B.3), we may write *' = (r + Z)T. (B.6) This may be written 131

132 l' = C'(o)T (B.7) where C' = + Z. (B.8) By integrating Equation (B.37), we get * = CT + R B. (B.9) where R. is a constant corresponding to rigid body rotation. We shall assume r vanishes at 0 = + 3; hence we may set R to zero and write ~ = C(O) T(r) (B.10) where C = (e + )i + f Zda (B.11) -P By substituting Equations (4.3), (B 35) and (B 4) into Equations (4.1a) followed by (4.4b) and then canceling a T we get -pZ - m" + Z" - Z - 1 = O (B.12) a After rearranging we have for the m'th mode 2 ^m (1 - pm)Zm - m"+ Z = (B.1) m m m 2 For brevity let us introduce the notation <[ ]> = f [ ] de (B.14) -5

1355 [ ]}) = [ ],= - [ ] =-1 (B.15) Then multiplying Equation (1,.13) by Z and integrating we have 0 = (1 - p )<Z Z > - <m" > <ZZ > + <ZZ > - < > (B16) m m n m n m n 2 m n. Interchanging subscripts and subtracting we get 2 2 (pZ P )< > - Z Z > - <m"Z > + <ZTZ > - <Z Z > n mmn mn nm mn nm m n 2 >+ — < > = (B.17) 2 mn 2 n From Equation (4.lb) we have n = o. (B.18) This implies that - <f1C > = 0 (B.19) m n Integrating Equation (B.19) by parts, we get - < C> = (mCn + < c>. (B.20) m n m n m n Using Equation (B.3), this becomes <- TmC > = - (C n + <rq > + <r Z >. (B.21) m n m n m n m n (.1 Interchanging subscripts and subtracting, we get - <T'C > + <rlC > = O = - { C] + (n C } + < rl Zn > - < Z > mn nm mn n mn nm (B. 22)

134 Appending Equation (B.22) multiplied by 1/a2 to Equation (B.17) and the rearranging we get 2 2 (P2 -2) <z z > Z -< "Z + <Z "Z > n m m n n m nm m n m 1 1 + 7- nC 7 [nrC. (B.23) 2 [~mCn 2 n a a Integrating the first two integrals on the right-hand side by parts twice, Equation (B.23) reduces to 2 2 (P2 p2) <Z Z > = (m'Z - (mIZ ) - (mZ + (m Zt P m n m mn nZm mn n m + - Cnm 2 n 7nCm n (B.24) a a. Since t and f vanish at both ends, this reduces to 2 2 (P2 P) <Z Z > = (m Z' - m Z') (B.25) n m m n n m m n For the elastic bending restraints given by Equations (4.13a) and (4.13b), the right side of Equation (B.25) identically vanishes and orthogonality is proved.

APPENDIX C QUANTITIES NEARLY INVARIANT TO THE SEMI-OPENING ANGLE 5 In Chapter II we found that the shapes of the critical configurations are relatively insensitive to the semi-opening angle P. The same thing was found to be true of the linear mode shapes in Chapter IV. In Chapter V we found that the system of ordinary differential Equations (5.19) does not depend explicitly upon 5; it has only a mild, implicit dependence through Relations (4.9c), (4.9d), and (4.9e). Thus it seems evident that there are several quantities that are nearly invariant to the semi-opening angle Q. Quantities found to be nearly invariant to P are = V r /( c ), an initial velocity parameter. 2- 2 P N/a, an invariant midsurface strain.* V = V/(Pa ), an invariant potential energy.** T = T/( )T an invariant kinetic energy. T /(P ), an invariant "Fourier" coefficient. Z, a linear radial displacement mode shape. m Al and A2, asymmetric and symmetric mode normalizing coefficients. lm 2m IY1/a, an invariant configuration vector-length. 5, a response parameter defined by Equation (6.4). *See, for example, Equation (A.6) for critical configurations. **For example, cf. Equations (3.70) and (3.80). 1355

LIST OF REFERENCES 1o Hoff, N. J. and Bruce, J. C., "Dynamic Analysis of the Buckling of Laterally Flat Arches," J. Math. and Phys., Vol. 32, No. 4 (January 1954), pp. 276-288. 2. Lock, M. H., "Snapping of Shallow Sinusoidal Arch Under Step Pressure Load," AIAA Journal, Vol. 4, No. 7 (July 1966), pp. 1249-56.. Lock, M. H., "A Study of Buckling and Snapping Under Dynamic Load," Air Force Report No. SAMSO-TR-68-100 (December 1967). 4. Bolotin, V. V., Dynamic Stability of Elastic Systems, Holden-Day, Inc., San Francisco, 1964. 5. Hsu, C. S., "On Dynamic Stability of Elastic Bodies with Prescribed Initial Conditions," Int. J. Eng. Sci., Vol. 4 (1966), pp. 1-21. 6. Fung, Y. C. and Kaplan, A., "Buckling of Low Arches or Curved Beams of Small Curvature," NACA TN 2840, California Institute of Technology (1952). 7. Kantorovich, L. V. and Krylov, V. I., Approximate Methods of Higher Analysis, Interscience Publishers, New York, 1958. 8. McIvor, I. K. and Popelar, C. H., "Dynamic Stability of a Shallow Cylindrical Shell," J. of the Eng. Mech. Div., ASCE, Vol. 93, No. EM3, Proc. Paper 5289 (June 1967), pp. 109-127. 9. Goodier, J. N. and McIvor, I. K., "The Elastic Cylindrical Shell Under Nearly Uniform Radial Impulse," J. Appl. Mech., Vol. 31 (1964), pp. 259-266. 10. Humphreys, J. S. and Bodner, S. R., "Dynamic Buckling of Shells Under Impulsive Loading," J. of the Eng. Mech. Div., ASCE, Vol. 88, No. EM2 (April 1962), pp. 17-36. 136

137 11. Fulton, R. E., "Dynamic Axisymmetric Buckling of Shallow Conical Shells Subjected to Impulsive Loads," J. Appl. Mech., Vol. 32, No. 1 (March 1965), PP. 129-134. 12. Budiansky, B. and Roth, R. S., "Axisymmetric Dynamic Buckling of Clamped Spherical Shells," NASA TN D-1510 (December 1962), pp. 597-606. 13. Simitses, G. J., "Dynamic Snap-Through Buckling of Shallow Spherical Caps," AIAA/ASME 7th Structures and Materials Conference, Cocoa Beach, Florida (1966). 14. Popelar, C. H., "Dynamic Stability of a Shallow Cylindrical Shell," Doctoral Dissertation, The University of Michigan, 1965. 15. Reissner, E., "On Transverse Vibrations of Thin, Shallow Elastic Shells," Quarterly of Appl. Math., Vol. 13, No. 2 (1955), pp. 169-176. 16. Lovell, E. G., "The Cylindrical Shell of Finite Length Under a Nearly Uniform Radial Impulse," Doctoral Dissertation, The University of Michigan, 1967.