Nonlinear Dynamics of Spherical Shells Buckling under Step Pressure
1 Nonlinear Dynamics of Spherical Shells Buckling under Step Pressure
Jan Sieber a , John W. Hutchinson b and J. Michael T. Thompson c a CEMPS, University of Exeter, Exeter EX4 4QF, UK b SEAS, Harvard University, Cambridge MA 02138, USA c DAMPT, University of Cambridge, Cambridge CB3 0WA, UK
Abstract
Dynamic buckling is addressed for complete elastic spherical shells subject to a rapidly applied step in external pressure. Insights from the perspective of nonlinear dynamics reveal essential mathematical features of the buckling phenomena. To capture the strong buckling imperfection-sensitivity, initial geometric imperfections in the form of an axisymmetric dimple at each pole are introduced. Dynamic buckling under the step pressure is related to the quasi-static buckling pressure. Both loadings produce catastrophic collapse of the shell for conditions in which the pressure is prescribed. Damping plays an important role in dynamic buckling because of the time-dependent nonlinear interaction among modes, particularly the interaction between the spherically symmetric โbreathingโ mode and the buckling mode. In general, there is not a unique step pressure threshold separating responses associated with buckling from those that do not buckle. Instead there exists a cascade of buckling thresholds, dependent on the damping and level of imperfection, separating pressures for which buckling occurs from those for which it does not occur. For shells with small and moderately small imperfections the dynamic step buckling pressure can be substantially below the quasi-static buckling pressure.
Keywords: dynamic buckling, spherical shells, nonlinear dynamics, imperfection-sensitivity
1. Introduction
Together with the axially compressed cylindrical shell, the complete elastic spherical shell under spatially uniform external pressure is one of the two archetypal examples of shell buckling. For this reason, it is studied extensively as a proving ground for advanced shell-buckling theories, involving the notorious imperfection-sensitivity. Recently there has been a renewed interest in this problem, stimulated by recent developments in shell buckling on several fronts. These include efforts underway in China, Europe and the U.S. to update design criteria for shell buckling accounting for advances in computation and shell manufacturing since the criteria were established in the mid- . Also there have been recent advances in experimental and theoretical aspects including new experiments on spherical shells with precisely manufactured imperfections [1], accurate formulations and simulations [1,2], and proposals to assess buckling and imperfection-sensitivity by experimental probing techniques [3] executed in [4,5]. Work on spherical shell buckling so far has concentrated on the static behaviour under the slow increase of the spatially uniform loading, under both (dead) pressure control and (rigid) volume control. Particular attention has been given to the imperfection sensitivity and associated energy barriers against collapse under operating conditions. Our aim in this paper is to initiate some high-precision dynamical studies of imperfect spherical shells, following a common historical pattern by first addressing the problem of step-loading. In this, we examine the highly nonlinear dynamical response of imperfect shells when, at rest under zero load, the shell is subjected to a rapidly applied dead pressure of magnitude p which then remains constant over the time, ๐ก , which is effectively allowed to tend to infinity. This loading process is then repeated at a fine set of different p values, and a record is kept as to which values of p result in a buckling collapse inferred by the dynamical response undergoing a dramatic increase in magnitude. The study limits imperfections and deformations to be axisymmetric which nevertheless captures all the essential nonlinearity of spherical shells buckling under uniform pressure. In this sense, the spherical shell is ideal for an in-depth investigation of dynamic buckling of imperfection-sensitive structures. At very large deflections, non-axisymmetric departures from the axisymmetric response can occur, though in a range far beyond that relevant to the onset of buckling [6]. In the simplest scenario at a fixed imperfection magnitude, ๏ค , the set of increasing step-loads will exhibit no collapse until some D p p = but guaranteed collapse at D p p ๏พ . We can then identify ( ) D p ๏ค as the dynamic buckling load at this ๏ค , and its graph can be compared with the static imperfection sensitivity curve ( ) S p ๏ค . This simple scenario does however often break down due to the dynamical phasing as the energy barrier is approached. Unlike predictions based on simpler one- and two-degree of freedom nonlinear structural systems [7], collapse is not guaranteed for D p p ๏พ , but instead we observe a sequence of thresholds. Here, recent work on nonlinear dynamics provides insights into why there is not a unique threshold characterizing dynamic buckling of the shells under step pressure. We believe the concepts revealed for the spherical shell will carry over to step loading of other shell structures. Earlier work on dynamic buckling, some of which is reviewed in [7], has not revealed the complexity of dynamic buckling nor the insights offered by perspectives of nonlinear dynamics. In the rest of this Introduction we list the content of the sections and anticipate some of the findings. The equations governing the nonlinear behaviour of complete thin, elastic spherical shells are presented in Section 2 along with a brief outline of the numerical solution methods. A summary of relevant quasi-static results for the elastic buckling of spherical shells is presented in Section 3, including the effect of imperfections on reducing the buckling pressure and the energy barrier to buckling for shells subject to sub-failure pressures. The paper focuses on prescribed spatially uniform (dead) pressure loadings where the device applying the pressure is not influenced by the deformation of the shell. Under quasi-statically increased load, buckling takes place at the maximum pressure the shell can support and is followed by complete dynamic collapse. Dimple-shaped geometric imperfections will be considered which are both realistic and among the most deleterious to buckling. The shape and amplitude of the imperfections will be scaled such that our results for buckling will be essentially independent of the all-important shell radius to thickness ratio, / R h , for shells with / R h greater than about 25 to 50. Under step loading, the dynamic coupling between spherically symmetric breathing (vibration) mode and the incipient buckling mode plays a critical role, and thus, in Section 3, we also present selected results on these modes based on linearization about both stable unbuckled states and unstable buckled states. Dynamic responses for step loaded shells are shown in Section 4. Comparison between the dynamic and quasi-static buckling pressures is made dependent on the imperfection amplitude. For nearly perfect shells or those with relatively modest levels of imperfection, the lowest dynamic buckling pressure falls significantly below the quasi-static buckling pressure but, conversely, falls far above lower bounds based on energy barrier concepts. For larger imperfections, the step buckling load is only slightly below the quasi-static buckling pressure. We observe a significant delay between application of pressure step and occurrence of buckling. Especially at or just above the lowest threshold for step buckling, buckling occurs only after multiple oscillations of the breathing mode with a slow transfer of energy from the breathing mode to the buckling mode. For larger imperfections ( / 0.25 h ๏ค ๏พ ) imperfection sensitivity trends for the step buckling pressure are similar to those for quasi-static buckling (at slightly lower pressures) and not strongly dependent on damping. Because of the time-delay before buckling occurs, damping comes into play in the case of small imperfections and thus in the determination of the step buckling threshold. Damping in this paper arises from, and can be controlled by, the time discretization employed in the numerical simulations. Section 5 is devoted to a discussion of the role of damping on the determination of the step buckling pressures. While, for larger imperfections we will see imperfection sensitivity trends for the step buckling pressure which are similar to (at slightly lower pressures) those for quasi-static buckling and not strongly dependent on damping. Section 6 examines the cascade of step-buckling thresholds from the perspective of work in nonlinear dynamics. Each threshold pressure corresponds to a response where the pole of the shell performs one more large oscillation before crossing the buckling threshold. This implies that the surface forming the buckling threshold in phase space has a complicated geometry folding around the unbuckled state many times.
2. Notation, governing equations and dimensionless quantities
This paper considers thin spherical shells of radius R and thickness h . The shell material is isotropic and linearly elastic with Youngโs modulus E , Poissonโs ratio ๏ฎ , and uniform density ๏ฒ . Geometric imperfections, I w , in the location of the shell middle surface will be introduced. All numerical results in this paper are based on the small strain-moderate rotation theory (Sanders [8], Koiter [9]) for axisymmetric deformations of the shell. Specifically, we assume that motion is rotationally symmetric around the North-South pole axis with reflection symmetric about the equator. This theory, reviewed and employed for the symmetric case in [2], is accurate for thin shells, e.g., / 50 R h ๏ณ , if the largest deflections which occur at the poles do not exceed about R , which will always be true in the range of interest in this paper. This paper is only concerned with whether the shell buckles or not โ no attempt will be made to attain the collapse state. The onset of buckling can be determined from simulations of deflections within finite multiples of the thickness ( to h ). (a) Governing equations according to small strain-moderate rotation theory All displacements can be described as functions of the meridional angle ๐ โ [0, ๐/2] ( ๐ = 0 at the equator, ๐ = ๐/2 at the North pole, as shown in an inset in Figure 4). Rotational symmetry implies that each middle surface pointโs displacement on the shell ha s a single tangential component ( , ) u t ๏ฑ and an outward normal component ( , ) w t ๏ฑ in the radial direction. The middle surface strains ( , ) ๏ฑ ๏ท ๏ฅ ๏ฅ and bending strains ( , ) K K ๏ฑ ๏ท in small strain/moderate rotation theory for these conditions of symmetry are, with
W U ๏ช ๏ข= โ + as the linearized rotation, I U W W ๏ฑ ๏ฅ ๏ช ๏ช ๏ข๏ข= + + โ , tan U W ๏ท ๏ฅ ๏ฑ= โ + (2.1) ๏ฑ ๏ซ ๏ช๏ข= , tan ๏ท ๏ซ ๏ช ๏ฑ= โ (2.2) Here, ( ) ( ) / ๏ฑ๏ข = ๏ถ ๏ถ , and the dimensionless displacements and bending strains are ( , ) ( , ) / W U w u R = and ( , ) ( , ) R K K ๏ฑ ๏ท ๏ฑ ๏ท ๏ซ ๏ซ = . Subscript ฮธ refers to the meridional direction, while ฯ refers to the circumferential direction. An initial imperfection in the form of a normal stress-free displacement of the shell middle surface, ( ) I w ๏ฑ , from the perfect spherical shape has dimensionless form / I I
W w R = . The imperfection used in this paper is a set of identical inward dimples at each pole with shape specified by (at the upper pole with / 2 ๏ข ๏ฐ ๏ฑ= โ ) ( / ) I I w e ๏ข ๏ข ๏ค โ = โ with / 1 / I B R h ๏ข ๏ฎ= โ (2.3) and ๏ค as the imperfection amplitude. The radius I ๏ข is a measure of the width of the Gaussian shaped dimple. These imperfections are realistic with / h ๏ค usually not larger than about unity [1,2]. The scaling of I ๏ข ensures that the relation between the buckling pressure and the amplitude ๏ค is essentially independent of / R h for thin shells. For a given imperfection amplitude ๏ค , there is a value of B that produces the minimum buckling pressure. The choice B = used throughout this paper gives nearly the minimum buckling pressure in the range h ๏ค๏ฃ ๏ฃ (c.f., [1,2]). This choice also ensures that the imperfection is confined to the pole diminishing to zero at distances on the order of Rh from the pole. The underlying assumptions for this theory are small strains, ๏ฅ ๏ผ๏ผ , and moderate rotation, ๏ช ๏ผ๏ผ . The resultant membrane stresses, ( , ) N N ๏ฑ ๏ท , and bending moment quantities, ( , )
M M ๏ฑ ๏ท , which are work conjugate to the corresponding strains, have dimensionless forms ( , ) ( , ) / ( ) n n N N Eh ๏ฑ ๏ท ๏ฑ ๏ท = and ( , ) ( , ) / m m M M R D ๏ฑ ๏ท ๏ฑ ๏ท = (2.4) where / [12(1 )] D Eh ๏ฎ= โ is the bending stiffness. These dimensionless stresses are related to the dimensionless strain quantities by ( ) ( , ) , / (1 ) n n ๏ฑ ๏ท ๏ฑ ๏ท ๏ท ๏ฑ ๏ฅ ๏ฎ๏ฅ ๏ฅ ๏ฎ๏ฅ ๏ฎ= + + โ , ( ) ( , ) , m m ๏ฑ ๏ท ๏ฑ ๏ท ๏ท ๏ฑ ๏ซ ๏ฎ๏ซ ๏ซ ๏ฎ๏ซ= + + (2.5) In terms of dimensional quantities, the principle of virtual work for axisymmetric behaviour is ๏ป ๏ฝ๏ป ๏ฝ /22 0 /2 /22 20 0 r R M K M K N E N E dR p w d R f w f u d ๏ฐ ๏ฑ๏ฑ ๏ฑ๏ฑ ๏ท๏ท ๏ท๏ท ๏ฑ๏ฑ ๏ฑ๏ฑ ๏ท๏ท ๏ท๏ท๏ฐ ๏ฐ ๏ฑ ๏ฐ ๏ค ๏ค ๏ค ๏ค ๏ฑ ๏ฑ๏ฐ ๏ค ๏ฑ ๏ฑ ๏ฐ ๏ค ๏ค ๏ฑ ๏ฑ+ + + =โ + + ๏ฒ ๏ฒ ๏ฒ (2.6) where a positive pressure p acts inward. The DโAlembert โacceleration forcesโ are / r f h w t ๏ฒ= โ ๏ถ ๏ถ and / f h u t ๏ฑ ๏ฒ= โ ๏ถ ๏ถ . The equilibrium equations generated by this principle, expressed in terms of the dimensionless quantities, are ( ) ๏ ๏ tan ( ) r m m n n n p f ๏ฑ ๏ท ๏ฑ ๏ฑ ๏ท ๏ฑ ๏ก ๏ช๏ข๏ข๏ข ๏ขโ โ + + + = โ + (2.7) tan tan m m n n n f ๏ฑ ๏ท ๏ฑ ๏ฑ ๏ท ๏ฑ ๏ฑ ๏ก ๏ช ๏ฑ๏ฉ ๏น๏ข ๏ขโ โ โ โ + =๏ซ ๏ป (2.8) Here, ( ) ( ) , , , cos , , , n n m m n n m m ๏ฑ ๏ท ๏ฑ ๏ท ๏ฑ ๏ท ๏ฑ ๏ท ๏ฑ= , R t ๏ก ๏ฎ= โ , ( ) ( ) / ๏ฑ๏ข = ๏ถ ๏ถ and cos , ( , ) cos , r pR W Up f fD ๏ฑ ๏ฑ ๏ฑ ๏ด ๏ด๏ฆ ๏ถ๏ถ ๏ถ= = โ ๏ง ๏ท๏ถ ๏ถ๏จ ๏ธ (2.9) with dimensionless time D thR ๏ด ๏ฒ= (2.10) For solutions symmetric about the equator of the shell, the boundary conditions at the pole and at the equator require U = , ๏ช = and W ๏ข๏ข๏ข = . The dimensionless system is defined by the parameters ๏ฎ , / R h , / h ๏ค and I ๏ข . Most of our results will be essentially independent of / R h because the values of this parameter chosen are large enough such that the solutions approach the solution limit for large / R h . (b) Discretization in space Equations (2.7) โ (2.9) form a nonlinear system of partial differential equations in space (angle ๐ โ [0, ๐/2] ) and dimensionless time ๐ . For quasi-static equilibrium computations (figures 1 and 2), this system is solved with zero acceleration forces ( ๐ ๐ = ๐ ๐ = 0 ) and varying ๐ as a free parameter to obtain the curves of equilibria, as shown in figures 1 and 2 below. The functions ๐ and ๐ are approximated by continuous piecewise polynomials of degree , consisting of pieces. This piecewise polynomial collocation approximation in the angle ๐ is supported by embedded boundary-value solvers for ordinary differential equations (ODEs) such as the collocation toolbox of COCO [11], which was used for figures 1 and 2 below, and similar solvers in previous publications on quasi-static buckling problems [2, 11]. The mesh in ๐ is non-uniform: the length of a subinterval is approximately close to the equator and approximately โ3 near the pole, equi-distributing an error estimate. (c) Discretization in time For dynamic simulations with a time step size ๐ฅ ๐ก we keep the space ( ๐ ) mesh constant over time, and use the BDF-2 rule for approximating time derivatives. BDF-2 approximates the time derivative ๐ฆฬ (๐ก) of function ๐ฆ(๐ก) at time step ๐ก ๐ = ๐๐ฅ ๐ก by the finite backward difference BDF (๐ฆ)(๐ก ๐ ) , which depends on the values of ๐ฆ at the current and previous two time steps, ๐ฆฬ (๐ก ๐ ) โ BDF (๐ฆ)(๐ก ๐ ) = ๐ก [๐ ๐ฆ(๐ก ๐ ) + ๐ ๐ฆ(๐ก ๐โ1 ) + ๐ ๐ฆ(๐ก ๐โ2 )], (2.11) where (๐ , ๐ , ๐ ) = (1.5, โ2,0.5) for ๐ > 1 and (๐ , ๐ , ๐ ) = (1, โ1,0) for ๐ = 1 (see [21]). We used the overdot to denote derivative with respect to dimensionless time ๐ , and dropped the argument ๐ in (2.11). For the dynamic simulation we solve at each time step ๐ก ๐ = ๐๐ฅ ๐ก the nonlinear system (2.7) โ (2.9) the same way as for equilibrium computations to obtain the solutions ๐(๐ก ๐ ) and ๐(๐ก ๐ ) . We introduce the additional variables (๐ ๐ (๐ก ๐ ), ๐ ๐ (๐ก ๐ )) approximating the time derivatives (๐ฬ (๐ก ๐ ), ๐ฬ (๐ก ๐ )) . In the (now non-zero) acceleration forces (see eq.(2.9)) we replace the term ๐ ๐/๐ ๐ by BDF (๐ ๐ )(๐ก ๐ ) in ๐ ๐ , the term ๐ ๐/๐ ๐ by BDF (๐ ๐ (๐ก ๐ )) in ๐ ๐ (both in (2.9), entering (2.7) and (2.8)), and add the equations ๐ ๐ (๐ก ๐ ) = BDF (๐)(๐ก ๐ ) and ๐ ๐ (๐ก ๐ ) = BDF (๐)(๐ก ๐ ) . This results in a closed nonlinear system of equations (as many equations as variables) for the variables (๐(๐ก ๐ ), ๐(๐ก ๐ ), ๐ ๐ (๐ก ๐ ), ๐ ๐ (๐ก ๐ )) at every time step ๐ก ๐ = ๐๐ฅ ๐ก . This system has the same structure as the problem solved for equilibrium computations. It can, thus, be solved with the same solver. In fact, in the limit ๐ฅ ๐ก โ โ , this nonlinear system for the dynamic simulation approaches the nonlinear system for equilibria. For dynamic simulations with a pressure ramp the variable ๐ is a given function of time. The discretization introduces an error of order ( ) t ๏ . (d) Energy, work and volume We conclude this section by listing some fundamental quantities in dimensionless form which will be employed in the paper. Symmetry about the equator allows the expressions below for the full shell to be reduced to integration over the upper half of the shell. At any stage of deformation, the strain energy in the shell, SE , is [2] ( ) ( ) SE R dD h ๏ฐ ๏ฑ ๏ฑ ๏ท ๏ท ๏ฑ ๏ฑ ๏ท ๏ท๏ฐ ๏ซ ๏ฎ๏ซ๏ซ ๏ซ ๏ฅ ๏ฎ๏ฅ๏ฅ ๏ฅ ๏ฑ ๏ฑ๏ฐ โ ๏ฌ ๏ผ๏ฏ ๏ฏ๏ฆ ๏ถ= + + + + +๏ญ ๏ฝ๏ง ๏ท๏จ ๏ธ๏ฏ ๏ฏ๏ฎ ๏พ ๏ฒ (2.12) The linearized expression for the decrease of volume of the shell from its unstressed state, V ๏ , is sufficiently accurate for this study and is given by /23 /2 cos2 V W dR ๏ฐ๏ฐ ๏ฑ ๏ฑ๏ฐ โ ๏ = โ ๏ฒ (2.13) The kinetic energy of the shell, KE , is KE W U dD ๏ฐ๏ฐ ๏ฑ ๏ฑ๏ฐ ๏ด ๏ด โ ๏ฌ ๏ผ๏ถ ๏ถ๏ฏ ๏ฏ๏ฆ ๏ถ ๏ฆ ๏ถ= +๏ญ ๏ฝ๏ง ๏ท ๏ง ๏ท๏ถ ๏ถ๏จ ๏ธ ๏จ ๏ธ๏ฏ ๏ฏ๏ฎ ๏พ ๏ฒ (2.14) with the dimensionless time ๏ด defined in (2.10). The cumulative work done by the pressure, p W , for a shell that is unstressed at t = is /20 /2 p W Wd p dD ๏ด ๏ฐ๏ฐ ๏ด ๏ฑ๏ฐ ๏ด โ ๏ถ= โ ๏ถ ๏ฒ ๏ฒ (2.15) While the governing equations (2.7)-(2.9) conserve the energy balance p W SE KE = + , the time discretization introduces some loss of energy (damping) that increases with โ t. In our discussion of the results we will always specify the time step โ t and discuss the relation between damping and โ t in Section 5. Finally, we will present results using a second dimensionless time defined by ห / t T ๏ด = , where T is the period of the sinusoidal spherically symmetric vibration mode of the unpressurized shell, also known as the breathing mode, (correspondingly, ฯ =2ฯ/ T is the breathing frequency) and readily derived as T RE ๏ฐ ๏ฎ ๏ฒ๏ฐ๏ท โ= = (2.16)
3. Selected results for buckling under uniform pressure relevant to dynamic buckling
As background for the dynamic buckling, we present a brief overview of results for the buckling of a spherical shell subject to a quasi-statically applied uniform inward pressure.
Figure 1 reminds the reader of the axisymmetric buckling of a complete spherical shell under spatially uniform external pressure, presented here based on our formulation in Section 2 and using the current notation. In each graph the vertical axis displays the ratio of the pressure to the classical critical buckling pressure, C p , of the perfect shell from the linearized analysis. The horizontal axis in figure 1a is the inward deflection at the pole divided by the shell thickness h , while the horizontal axis in figure 1b is the volume decrease normalized by the volume decrease of the perfect shell at the classical buckling pressure, C V ๏ . The results based on moderate rotation theory in Fig. 1 agree closely with results computed independently using FEM modeling in [1]. The classical values are [1] C E h Rp ๏ฎ= โ and C R hV ๏ฐ ๏ฎ๏ฎโ๏ = โ (3.1) Both graphs display the static equilibrium states of imperfect shells with the imperfections indicated. The upper curve for the smallest imperfection, / 1 / 640 h ๏ค = , is a close approximation to the behaviour of the perfect shell. Under the slow quasi-static increase of the controlled (dead) pressure, the buckling pressures (the maximum pressures) are indicated by small black dots. At these limit points (called folds or saddle-node bifurcations in dynamics) the imperfect shell in a noise-free environment will jump dynamically (snap buckle) to a collapsed state outside the range of this theory. Figure. 1 Buckling behaviour of spherical shells with / 100
R h = , ๏ฎ= for various dimple imperfections having B = . (a) Pressure versus pole deflection. (b) Pressure versus change in volume with the energy barrier to buckling illustrated for a prescribed pressure, / 0.25 C p p = , for the shell with / 3 / 5 h ๏ค = . As the load increases slowly towards these limit points, a shell is in a more and more precarious meta-stable state protected by a diminishing energy barrier against small disturbances (e.g., noise) typically present in a physical experiment. The magnitude of this barrier Barrier E of a given shell at a prescribed (fixed) pressure p can be identified as the โtriangular - shaped โ area on the plot of ( ) p V ๏ in figure 1b. This area is simply the difference of the energy of the shell/loading system between the unstable buckled state on the falling segment of the curve (a dynamicist โs saddle) and the stable un-buckled state on the rising segment (a dynamicist โs node). This area is the difference in the strain energy of the shell in the two states less the work p V ๏ that would be performed by the external pressure. Energy barriers have been accurately calculated in [6, 12] for perfect and imperfect spherical shells under both (dead) controlled pressure and (rigid) controlled volume loading conditions. The barrier for prescribed pressure has been recomputed and presented in figure 2 using the same imperfection amplitudes employed in figure 1. For prescribed pressure, the quasi-static system energy is SE p V ๏น = โ ๏ . The results shown in figure 2 for various levels of imperfection are independent of / R h , to a very good approximation for thin shells with / 50
R h ๏ณ . The energy barrier vanishes at the buckling pressure and remains very small for pressures or volume changes somewhat below the buckling value. At pressures well below the buckling value the energy barrier increases dramatically and becomes relatively insensitive to the imperfection level. The energy barrier in figure 2 is normalized by / C C p V C h R ๏ where C ๏ฎ ๏ฎ๏ฉ ๏น= โ โ๏ซ ๏ป . Note that C C p V ๏ is the strain energy in the perfect shell at the classical buckling pressure. It follows from figure 2 that, because of the factor / Ch R , the energy barrier is a very small fraction of the energy stored in a thin shell (or, equivalently, of the work done on the shell by the pressure) . Moreover, the ratio of the energy barrier to the stored energy decreases for thinner shells directly in proportion to / h R . This is due to the fact that the deformation in the buckled state is localized at the pole in the form of a dimple whose size scales with Rh and thus decreases in size relative to the shell itself as / h R diminishes, as will be discussed more fully later. The implication of this will be discussed later with regard to dynamic buckling under step loading. In the limit for very thin shells, / 0 h R โ , the energy barriers for prescribed pressure and prescribed volume change are the same. The barrier in the case of prescribed volume change does depend somewhat on / h R , as discussed in [11]. Figure. 2 Energy barrier, barrier E , for shells loaded to prescribed dead pressure p . The barrier is presented for various levels of imperfection; these results are essentially independent of / R h for / R h greater than about 50. For these results, / 100
R h = , ๏ฎ= and B = . We conclude this section a few additional results relevant to dynamic buckling. The vibration frequency of the breathing mode of the perfect shell ( / 0 h ๏ค = ) undergoing small spherically symmetric oscillations ( ( ) 0, ( ) 0 w u ๏ฑ ๏ฑ๏ข = = ) and subject to no applied pressure introduced earlier in (2.8) is (equivalent to (2.16) E R ๏ท ๏ฎ ๏ฒ= โ (3.2) The buckling mode of the perfect shell associated with the lowest eigenvalue of the classical buckling problem, C p , is also an important reference mode. The lowest vibration frequency of the perfect shell vanishes at C p and the associated mode is the classical buckling mode. The normal deflection of the classical mode can be expressed in terms of a Legendre polynomial โ explicit expressions are given in [2,13]. The spectrum of frequencies of the modes linearized about applied pressures below C p is also revealing and relevant to the understanding of dynamic buckling. An illustration is presented in figure 3 which shows a selection of normalized complex frequencies, / ๏ท ๏ท , and associated modes for a shell with imperfection amplitude, / 1 / 4 h ๏ค = , and subject to pressure C p p = . Note from figure 1 that there are two equilibrium solutions associated with C p p = , one the stable unbuckled state (written here as โnodeโ) where the pole deflection is approximately w h = โ and the other the unstable equilibrium point at the saddle point of the system energy (written as โsaddleโ) where the pole deflection is approximately w h = โ . The results in figure 3 are obtained by linearizing the equations about these two equilibrium solutions. The time dependence of the linearized solution has the form t e ๏ท where ๏ท is the complex frequency for the respective mode. Note that with normalization used in figure 3a, the reference breathing mode for the perfect unpressured shell has / i ๏ท ๏ท = ๏ฑ . The spectrum of frequencies in the two states at C p p = are plotted in figure 3a and two of the most important associated modes shapes for each state are plotted in figure 3b. Since the shell has an imperfection there is no strictly spherically symmetric motion, but the mode identified as the breathing mode in the unbuckled state, which has / ๏ท ๏ท nearly equal to i ๏ฑ (pair of black dots in figure 3a identified by A). The associated normal deflection of mode A is plotted in the upper half of Fig. 3b deviates from the uniform normal deflection of the breathing mode of the perfect shell due to the imperfection. If one t racks this โbreathingโ mode through the equilibrium solutions to the buckled state at C p p = , one finds that the normalized frequencies hardly change from i ๏ฑ ( pair of black dots labeled D in the saddle spectrum), and the associated mode has even more distinct ๏ฑ -variations associated with the non-uniformity of the buckled state about which the linearized solution has been obtained. Next, we focus on the lowest vibration frequency of the unbuckled shell at C p p = for the mode identified as the โbucklingโ mode, which as seen by the second set of dots labeled C in Fig. 3a has / 0.35 i ๏ท ๏ท ๏ ๏ฑ . The normal deflection of the associated mode plotted in figure 3b is very similar to that of the classical buckling mode. Tracking this mode to the unstable buckled state at C p p = leads to the โbuckling modeโ with / 0.34 ๏ท ๏ท ๏ ๏ฑ in figure 3a (labelled B) corresponding to exponential growth/decay. An important feature regarding the buckling mode shape is that the undulations have localized to the polar region in the form of a dimple. Figure. 3 Linearized modes about the two equilibrium states at C p p = , the stable unbuckled state (node) and the unstable buckled state (saddle). The shell has imperfection, / 1 / 4 h ๏ค = and B = , with / 100 R h = and ๏ฎ= . The normalized spectrum of the frequencies, / ๏ท ๏ท , is given for the two states in a) and the radial displacements of two of the most important modes are presented in b). The frequencies and their associated modes are identified and discussed in the text. All modes have been scaled to have maximum modulus 1.
4. Buckling under step loading with spatially uniform pressure
As noted, attention in this paper is limited to the response of spherical shells under pressure loads that are spatially uniform. In subsequent sections we will consider dynamic buckling under time varying spatially uniform pressure, ( ) p t , in which the shell starts at rest in an unpressured state. From the stationary starting state, the pressure is ramped up rapidly and then held constant at its end value, which we also call p (without argument t ) in our discussions. We refer to this loading as step loading. For the limit in which ramping is instantaneous, we use the terminology โ instantaneous step loading โ . Figure 4 shows the dynamic response of several quantities of interest for a step loaded imperfect shell with / 0.25 h ๏ค = such that its quasi-static buckling pressure is ( ) 0.57 S C p p ๏ค = . The pressure is ramped up rapidly from 0 at / 1 t T = to C p p = at / 1.25 t T = (about 10% below the static buckling pressure) and held at C p for / 1.25 t T ๏พ . The plots in figure 4a show the time-variation of the work done by the pressure, / 2 p W D ๏ฐ , the strain energy in the shell, / 2 SE D ๏ฐ , and the kinetic energy, / 2 KE D ๏ฐ . In the simulation run for figure 4 the time step size โ t is chosen sufficiently small such that the energy balance, ( ) / 2 0 p SE KE W D ๏ฐ+ โ = , is satisfied to a high degree of accuracy (see Section 5 for a detailed discussion of damping).
Figure. 4 Time profiles of applied pressure, / C p p , together with the cumulative work of pressure, / 2 p W D ๏ฐ , strain energy, / 2 SE D ๏ฐ , kinetic energy, / 2 KE D ๏ฐ , volume change, / C V V ๏ ๏ , and pole deflection, / w h . These are for an imperfect shell with imperfection amplitude h ๏ค= , The other parameters are ๐ /โ = 100 , ๐ = 0.3 , ๐ต = 1.5 with time step size ๐ฅ ๐ก = 0.022 . The static buckling pressure is C p and the steady dynamic pressure in these simulations is C p . The inset shows the geometry underlying Eqs. (2.1-2.9). The lower plots in figure 4b show the prescribed time variation for the pressure and the variation of the volume change, / C V V ๏ ๏ , and the pole deflection, / w h . Note that for seven cycles the strain energy, work done by the pressure and the volume change are nearly sinusoidal with period T . During these cycles the kinetic energy has a period / 2 T as is typical for vibratory systems. Prior to / 7 t T ๏ the motion is dominantly a spherical symmetric โbreathingโ oscillation with ( , ) 0 w u ๏ข = . However, from the beginning, the deflection at the pole responds differently from most of the shell and, for this example, starting at roughly / 7 t T ๏ conditions at the pole give rise to localized snapping into a dimple buckle. In all cases in this paper, an inward pole deflection exceeding 3 or 4 time the thickness results in buckling of the spherical shell. Once the pole deflection reaches this magnitude the shell cannot resist dynamic collapse. The buckling behaviour is brought out more clearly in figure 5 which displays two representations of the normal deflection of the shell and one representation of the in-plane displacement as functions of both position ๏ฑ and time. By / 8 t T = it is evident that the deflection has taken the form of a dimple localized at the pole surrounded by spatial undulations decaying away from the pole. The amplitude of the dimple at the pole doubles between / 8 t T = and / 9 t T = . Under the fixed pressure the dimple grows unabated leading to complete collapse of the shell. The in-plane displacement is roughly two orders of magnitude smaller than the normal deflection, which is typical for spherical shell buckling. Two aspects of the dynamic process stand out from figures 4 and 5: the significant delay in the formation of a buckle until about 6 or 7 overall oscillations of the shell in this particular case, and the localization of the buckle at the pole as it develops. Further insight into the delay in buckling will emerge when results are presented shortly for the responses of the shell to a full range of step pressures. The localization helps to explain why the various energy variations of figure 4 are very large compared to the energy barrier to buckling, and it will be useful at this point to highlight these energy differences. Note that for the example in figure 4 the cyclic variations in the strain energy and kinetic of the shell prior to buckling are / 2 5 SE D ๏ฐ ๏ and / 2 1
KE D ๏ฐ ๏ . At a pressure / 0.512 C p p = , the energy barrier separating the static un-buckled and buckled states for an imperfect shell with / 0.25 h ๏ค = can be obtained from figure 2 as / [ / ] 0.00436 Barrier C C
E p V Ch R ๏ ๏ . The conversion between the two normalization factors in these dimensionless energy ratios is ( ) C C
Ch hp V DR R ๏ฐ๏ฎ๏ = โ (4.1) For the shell with / 100
R h = and ๏ฎ= , the conversion factor is h R ๏ฎโ = and thus / 2 0.000633 Barrier
E D ๏ฐ ๏ . The energy barrier is a tiny fraction of the energy variations taking place in the shell. Apart from one aspect mentioned later related to the level of imperfection, the barrier has essentially no quantitative relevance to the uniform pressure step loading because only a small region of the shell near the pole participates in the buckling process. Most of the shell undergoes breathing motion with ( , ) 0 w u ๏ข ๏ which accounts for nearly all of the energy variations seen in figure 4. The coupling between the breathing motion and the emerging dimple buckle at the pole requires seven cycles before buckling occurs. Figure 5. Normal deflection, ( , ) / w t h ๏ฑ , (with the in-plane displacement, ( , ) / u t h ๏ฑ , included on the right) plotted as a function of position, ๏ฑ in radians, and time, / t T , for the step loading example presented in figure 4. The dynamic responses at the pole of an imperfect shell with / 0.1 h ๏ค = and subject to various levels of step loading is shown in Fig. 6a revealing whether buckling takes place or not. The companion plot in figure 6b, perhaps the most important figure in this paper, summarizes the dynamic buckling behaviour under step loading over the range of imperfections from / 0 h ๏ค = to / 0.6 h ๏ค = . We systematically ran a sequence of simulations for a range of amplitudes ๏ค (from โ/640 to ) and a range of rapid uniform pressure ramps from up to a final pressure ๐ . Each dot in the right panel of figure 6 corresponds to one simulation: a dot at ( ๏ค , ๐) means that a simulation was run with pressure increasing from to ๐ between ๐ก/ T = 1 and and then held constant up to time ๐ก/ T = 26 . If the pole deflection dropped below โ4โ at a time ๐ก ๐ / T before the end time, we record the simulation as โbuckledโ, colo uring the dot red and indicate the delay (๐ก ๐ โ 1.5)/ T in the contours. Otherwise, the dot is coloured green indicating that the simulation did not show buckling. The left panel of figure 6 shows time profiles of the pole deflection ๐ค (๐ก) for a sequence of simulations for a fixed imperfection amplitude / 0.1 h ๏ค = . Time profiles in green colour did not buckle before ๐ก/ T = 26 , those in dark red did buckle. The transparent surface shows the static saddle equilibrium value for the pole deflection (with a small part of the node equilibrium surface close to the fold value S C p p = ). Blue ellipses indicate where the time profiles cross this saddle surface. All time profiles that cross the saddle surface wound up buckling. Figure 6: Parameter scan over end pressures ๐ for pressure ramp times of / 2 T and dynamic buckling summary for imperfection amplitudes in the range h ๏ค๏ฃ ๏ฃ . Each dynamic simulation is carried out from t T ๏ฃ ๏ฃ . Left panel: Sequence of time profiles for pole deflections ๐ค (๐ก) for fixed / 0.1 h ๏ค = and varying end pressure ๐ ( evenly spaced values between ๐ and C S p p = ). Right panel: Overview of simulation results in the (๐, ๏ค ) plane showing parameter values where buckling occured (red dots) or did not occur before ๐ก/ T = 26 (green dots). The static buckling pressure ( ) S p ๏ค is the black curve. The colour encodes the delay of the buckling after the end of the pressure ramp in units of ๐ . Parameters: ๐ /โ =100 , ๐ = 0.3 , ๐ต = 1.5 , ๐ฅ ๐ก = 0.087 . Figure 6 permits two observations. First, there is not a clear-cut buckling threshold for step loading. In the left panel a ramp up to a final pressure of ๐ leads to buckling, while a ramp up to a somewhat larger final pressure of ๐ does not give rise to buckling. Second, for imperfections that are not too small (for / 0.15 h ๏ค ๏ณ ), the difference between the static buckling pressure and the lowest step pressure causing buckling is uniformly small (about ๐ ), while for small imperfections (e.g., / 0.15 h ๏ค ๏ผ ) the difference is significantly larger, with reductions up to ๐ (when / 0 h ๏ค ๏ ). Especially for shells with small imperfections, the static buckling pressure is not an accurate estimate of the buckling pressure for a step pressure loading.
This finding is at odds with the dynamic buckling predictions for step loading based on 1-degree-of-freedom (DOF) models analyzed in [14,15]. These authors considered two types of imperfection-sensitive, 1-DOF models: one with unstable symmetric bifurcation behaviour (with a cubic nonlinearity) and the other with asymmetric bifurcation behaviour (with a quadratic nonlinearity). In each case, for every level of imperfection, it was possible to relate the instantaneous step buckling load, D ๏ฌ (generalizing p D ), to the static buckling load, S ๏ฌ (generalizing p S ), and the static buckling load of the perfect model, C ๏ฌ (generalizing p C ). For the symmetric model, this relation is D C DS C S ๏ฌ ๏ฌ ๏ฌ๏ฌ ๏ฌ ๏ฌ๏ฆ ๏ถโ= ๏ง ๏ทโ๏จ ๏ธ (4.2) while for the asymmetric model it is D C DS C S ๏ฌ ๏ฌ ๏ฌ๏ฌ ๏ฌ ๏ฌ๏ฆ ๏ถโ= ๏ง ๏ทโ๏จ ๏ธ (4.3) The relations between the step buckling load and the static buckling load from the 1-DOF models discussed above are in complete agreement with asymptotic results obtained by Thompson [16] derived from a general ( n - DOF) formulation for discrete elastic systems for which the perfect system has a unique buckling mode. The analysis is purely quasi-static and asymptotic for small imperfections, but did compare well with step-buckling experiments on structural frames of the type built by J. Roorda [24]. For both unstable symmetric and asymmetric bifurcations, Thompson determined: 1) the relation between the static buckling load S ๏ฌ (the max load) and the imperfection, and 2) the relation between the โ astatic load โ , N ๏ฌ , and the imperfection. For a given imperfection, the astatic load is that load ๏ฌ at which the work done by the fixed ๏ฌ equals the strain energy in the unstable buckled state. The asymptotic results for the astatic load coincide with asymptotic limits of (4.2) and (4.3) for small imperfections if N ๏ฌ is identified with D ๏ฌ . For instantaneous step loading of the 1-DOF models discussed above, it is straightforward to prove that the astatic load N ๏ฌ must be a lower bound to the dynamic buckling load D ๏ฌ . The fact that N D ๏ฌ ๏ฌ= for these models is due to the fact that, at the lowest step load for which the model buckles, the model attains the static unstable buckled state, momentarily coming to rest, such that the astatic condition is exactly satisfied. Applied to our spherical shell, the astatic condition is easily visualized as shown in figure 7. When the two areas, Area and Area , on the pressure-volume plot in figure 7 are equal, the astatic condition for N p is met, i.e., N B B p V SE ๏ = with B denoting the unstable static buckled state. At pressures below N p , the strain energy in state B exceeds N B p V ๏ , and vice versa. Thus, if an instantaneous step-load occurred to a pressure lower than N p , the loading system would not be able to provide enough work to reach the saddle point represented by B . Although we have not proved that the astatic pressure, N p , is a lower bound to the instantaneous step buckling pressure, it seems likely that this is the case. For N p to qualify as a strict lower bound (even an excessively conservative one) one would have to establish that there are no other mountain passes associated with other unidentified saddle equilibria. However, the detailed investigation of the thresholds later in Section 6 indicate that escape occurs indeed near the saddle used in computing the energy barrier presented in figure 2 and discussed further in connection with figure 7. Figure. 7 Illustration of equal area construction for determining the astatic pressure N p based on the condition N B B p V SE ๏ = . The case illustrated for / 1 / 640 h ๏ค = has / 0.188 N C p p = . We have computed the normalized astatic pressure / N C p p as a function of the imperfection amplitude for the shells considered in figure 1. We find that the astatic pressure is / 0.2
N C p p ๏ป with almost no dependence on the imperfection amplitude for imperfections in the range / 0.6 h ๏ค ๏ฃ . Compared with the step buckling pressures in figure 6, the astatic pressure is unrealistically low and of little predictive value, at least for the damping level associated with the results in figure 6. Hoff and Bruce [17] made an early use of the astatic load in their study of the buckling of shallow arches subject to step loading of a pressure distributed along its length. The shallow arch is like the spherical shell in that it undergoes dramatic changes in deflection when buckling occurs โ so called snap buckling. It differs in that the entire arch buckles while the buckling deflections of the sphere are localized near the pole. In the one specific example Hoff and Bruce considered, they found the dynamic step loading prediction agreed quite well with the astatic pressure, both giving an estimate of the dynamic buckling pressure that was about 20% below the static buckling pressure. This result is similar to what one might expect based on the 1-DOF models discussed above and on similar models in the book on dynamic stability [7]. The relation between the dynamic and static buckling loads for the spherical shell is at odds with the corresponding relations for the simple 1-DOF models in equations (4.2) and (4.3) as seen in figure 8. For the shell, the largest reductions of the dynamic step buckling pressure relative to the static buckling pressure occur for the shells with the smallest imperfections. By contrast, the dynamic buckling load of the simple models is only slightly reduced from the static buckling load when the imperfection is small. For the models, the largest relative reductions occur for the largest imperfections, while for the spherical shell the opposite is true. Figure. 8 Ratio of the dynamic buckling load under step loading to the static buckling load plotted as a function of the fractional reduction of the static buckling load from the static buckling load of the perfect structure. Predictions of the 1-degree of freedom models: symmetric model from (4.2) and asymmetric model from (4.3). The trend line from figure 6 for the buckling of the spherical shell under spatially uniform pressure is shown. The main factor at play in the different dynamic buckling behaviours of the simple models and the spherical shell is associated with the interaction among the different modes activated in the step-loaded shell. The initial response of the sphere is dominated by the oscillatory motion of the breathing mode which absorbs most of the work done by the pressure. When the step pressure is sufficiently large the nonlinear coupling between the breathing mode and the incipient dimple causes the dimple to grow and to snap buckle leading to full collapse of the shell (under the prescribed pressure considered here). The simple 1-DOF models discussed earlier do not encompass modal interaction. Dating from the early work of Goodier and McIvor [18] on the buckling of long cylindrical shells (effectively rings) under dynamic radial pressure there has been a large literature on the coupling of breathing and buckling modes, often leading to a Mathieu equation governing the early stages of the coupling. The ring buckling problem considered in [18] has this form but it is not imperfection-sensitive and the nonlinearity is such that snap buckling does not occur. Instead, in their problem the nonlinear mode interaction gives rise to a gradual amplification of the buckling mode. Tamura and Babcock [19] carried out an early nonlinear mode interaction analysis for step loading of a finite length, imperfect cylindrical shell under an axial step load. This structure/loading combination is imperfection-sensitive. The oscillation of the axial compressive stress (the breathing mode in this case) excited by the step load was treated approximately and coupled to two interacting buckling modes. The authors analyzed only one specific imperfect shell for which the dynamic buckling load associated with an abrupt increase in the shell deflection was found to be approximately 60% of the static buckling load. More recently, the dynamic buckling of conical shells under step loaded axial compression has been investigated [20]. This problem also has features in common with the spherical shell problem in that the structure/loading system is imperfection-sensitive and results in snap buckling once the buckling mode is sufficiently amplified. In plots of the ratio of the step buckling load to the static buckling load as a function of the imperfection level, results [20] show a trend similar to that in figure 8 for the spherical shell. In particular, they find that conical shells with small imperfections have ratios of dynamic to static buckling as low as about 0.6 and that this ratio increases for larger imperfections, similar to the trend in figure 8. The authors of [20] suggest that their numerical results apply to conditions where damping is negligible, and they do not identify the cascade of buckling thresholds of figure 6. We will return to these issues in the next section. The plots of the energy barrier for the spherical shell in figure 2 also shed some light on the trend for the dynamic buckling pressure for the spherical shell in Fig. 8. Note that for small imperfections the energy barrier remains very small for values of p as low as 60 to 70% of the static buckling load whereas for larger imperfections the energy barrier remains small for smaller reductions of p relative to the static buckling load. This is consistent with figure 8: a relatively perfect shell is more susceptible to buckling at pressures within a given fraction of its static buckling load than more imperfect shells loaded to the same fraction of their static buckling load. When snap buckling requires several oscillations of the breathing mode of the shell, as illustrated in figures 4-6, it is obvious that damping effects will influence the dynamic buckling load. Damping is present for these results associated with the numerical algorithm used in solving the dynamic equations. Section 5 which follows discusses some of the issues related to this algorithm and the role of damping in dynamic buckling.
5. Balance between damping and nonlinear coupling between modes
Even though the small strain moderate rotation theory does not include any dissipation, some damping is introduced by the numerical BDF-2 time stepper (2.11) for the simulation. As the spectrum of equilibria in figure 3 suggests, the introduction of damping is necessary to make dynamic simulations numerically feasible. Without damping small disturbances of equilibria will lead to near quasi-periodic behaviour composed of oscillations with arbitrarily high frequencies, where the range of frequencies is determined by the resolution of the space discretization. The importance of damping in regularizing the numerical analysis of dynamic structural systems features prominently in modern treatments of the subject [21]. In particular, the introduction of damping has a strong effect on the long-time behaviour of a conservative system such as (2.7) โ (2.9). To provide a good estimate of this effect on buckling thresholds, we recap briefly how much damping a time stepper based on the BDF-2 approximation introduces. We also illustrate this effect in figure 9. Figure. 9 (Left) Change of spectrum near stable equilibrium (node) from original conservative case (blue) by numerical discretization with small and larger time stepsize. (Right top) Volume oscillations and loss of energy
๐พ๐ธ + ๐๐ธ โ ๐ ๐ for small and larger stepsize. (Right bottom) Pole deflection for small and larger stepsize. Other parameters: ๐ฟ/โ = 0.25 , ๐/๐ ๐ = 0.512 . (a) Damping at linear level The amount of damping a numerical scheme introduces is well understood only for linear systems. In this case one may study the behaviour of the time stepper for each linear mode separately, since the damping depends only on its frequency ๐ . The damping at frequency ๐ is determined by inserting the numerical approximation (in our case the BDF-2 (2.11)) into the linear ODE ๐ฆฬ = i๐๐ฆ . The solution of ๐งฬ = i๐๐ง ( ๐ > 0 ) asymptotes to ๐ง(๐ก) = exp(๐ BDF ๐ก) , where BDF tt t d ๏ท ๏ท๏ฌ ๏ท ๏ท+ + ๏= + + ๏ = ๏ โ ๏ (5.1) Both the numerical growth rate ๐ and frequency shift ๐ฅ ๐ are negative for time step size ๐ฅ ๐ก > 0 (going to zero for ๐ฅ ๐ก โ 0 ) such that the time stepper introduces artificial numerical damping โ๐(๐, ๐ฅ ๐ก ) and a slow-down per period โ๐ฅ ๐ /(2๐) for a mode with frequency ๐ . The damping โ๐(๐, ๐ฅ ๐ก ) can be approximated over the range of frequencies shown in the spectra in figure 3 (to roughly twice the breathing frequency) by expanding the real part of (5.1) in ๐ = 0 : ( , ) 4 10 tt t d ๏ท๏ท ๏ท๏โ ๏ = + ๏ (5.2) an approximate quartic in the frequency ๐ . Figure 9 shows the effect and the amount of damping for two different step sizes. The smaller stepsize, ๐ฅ ๐ก = 0.022๐ , was used for the single example trajectories shown in figures 4 and 5, the larger stepsize, ๐ฅ ๐ก = 0.087๐ (four times the smaller stepsize) was used for the parameter study in Figure 6. Otherwise, all parameters are identical to figures 3, 4 and 5. The left panel of figure 9 shows that the numerical scheme introduces frequency dependent damping, suppressing high-frequency oscillations more strongly, according to the approximately quartic frequency-damping relationship (5.2). Thus, a single small-amplitude breathing oscillation around the stable equilibrium gets damped by less than % for ๐ฅ ๐ก = 0.022๐ but by % for ๐ฅ ๐ก = 0.087๐ (in one unit of time by our scaling). (b) Damping of shell motion after pressure ramp The top right of figure 9 shows that the damping factors derived for small-amplitude breathing oscillations carry over to the motion of the shell after the pressure ramp as in figure 4. The volume oscillations are dominated by breathing oscillations and the decay rate of these larger scale breathing oscillations matches the predictions from the linear approximation exp(๐ bdf ๐ก) (dotted curves in figure 9, top right panel). The dashed curves show the loss of energy along the trajectory, which is also times higher for the large stepsize ๐ฅ ๐ก = 0.087๐ . We also observe an additional loss of energy during the rapid ramp-up of the pressure in the time period from ๐ to , which is not directly predictable from linear theory. However, this loss of energy is, consistently, also higher for the larger stepsize. (c) Conclusion for calibration of damping Real shells and other numerical discretizations may have damping with frequency dependence different from the one shown in figure 1 and approximated by expansion (5.2). Figure 9 suggests that in these cases damping should be compared at the breathing frequency. In experiments the damping of the breathing vibration can be measured by applying small pressure load ramps far from buckling pressure. According to figure 9 this linear damping carries over to the motion with larger amplitude. Damping in structures has many possible sources, including air resistance, dissipation at joints and boundaries and material damping of various kinds. While the damping generated by the numerical discretization used in the present study may not represent all the physical sources of damping, it is clear from figure 9 and Equation (5.1) that the damping in results from figures 4-6 is comparable to damping in other numerical schemes and empirical data (after calibration at the breathing frequency). The bottom right panel of figure 9, showing the motion of the pole deflection for both stepsizes (pole deflection of the saddle is shown for reference), demonstrates that the larger damping for the larger stepsize causes the shell to avoid buckling (while it does buckle for lower damping at ๐ฅ ๐ก = 0.0022๐ as shown also in figures 4 and 5). We expect this to hold in general โ lower damping implies lower buckling threshold. One argument for this is given in the next section.
6. A cascade of buckling thresholds for non-zero damping
For the nonlinear dynamicist the study in this paper raises a number of interesting fundamental questions, and points to their relevance in practical applications. To examine some of the issues, let us focus on a conservative autonomous system, as is our spherical shell after the pressure has been step-loaded to a fixed value. Additionally, assume that there is only a single potential energy saddle and barrier-height V S that is preventing escape to a โremoteโ region of phase space (such a single saddle is not rigorously established for our shell buckling). This might be thought of as a well understood problem, but this is not the case, especially because our system has many degrees of freedom: strictly an infinite number, but we will treat the shell for simplicity as if it has a large number of mechanical degrees of freedom with a phase-space of N dimensions (twice the DOF number). Within these restrictions, there is a lower bound for both damped and undamped systems because a trajectory starting with total (kinetic plus potential) energy, E , cannot escape if E < V S . The question that remains is what happens if
E >V S and the situation is remarkably obscure. Even in the extreme case of no damping and with the elapsed time going to infinity, there is no guaranteed escape due to a number of complex blocking actions. These are still being explored in the multi-body problems of astronomy and chemical kinetics. The systematic parameter study in figure 6 shows that buckling under step loading can occur in ranges of pressures that are far below the static buckling pressure ๐ ๐ , but also far above lower bounds given by energy barrier. For some imperfection levels ๐ฟ multiple buckling thresholds are visible. This section investigates the thresholds in more detail, using the dynamic buckling pressure thresholds for ๐ฟ = โ/4 as an example. (a) Centre-stable manifold of the saddle Considering the spectrum of the linearization in the saddle at ๐ฟ = โ/4 and / 0.512 C p p = in figure 3, we see that the saddle has one stable eigenvalue and one unstable eigenvalue. Their respective eigenvectors correspond to directions in which trajectories exponentially converge to or diverge from the saddle. Without damping, the saddle appears as linearly neutral in all other directions: but with a little non-zero damping these other directions would be stable modes with trajectories converging to the saddle. This implies that close to the saddle the set of all initial conditions that do not diverge rapidly from the saddle forms a hypersurface, splitting the phase space near the saddle into two subsets (and the boundary hypersurface). One subset contains those initial conditions that buckle immediately. The other subset contains those initial conditions for which trajectories do not buckle immediately but instead oscillate around the node and either ultimately buckle or possibly converge to the nearby node if damping is present. The boundary is the set of all initial conditions whose difference to the saddle is spanned by the eigenvectors corresponding to the stable and neutral directions (directions that are neutral without damping become weakly stable with damping). Mathematically this boundary set is called the centre-stable manifold (CSM). Close to the saddle this CSM is approximately a hyperplane, a linear space of co-dimension 1, that is, of dimension one less than the entire space ( N โ in our simplified argument). Further away from the saddle the CSM is no longer a hyperplane but a (differentiable) curved N โ -dimensional hypersurface. It is known that CSMs of saddles can fold back on themselves dramatically, even in low-dimensional systems, allowing the system to become chaotic [22,23]. The CSM of the saddle depends on the parameter ๐ (as does the location of the saddle itself). The sketch in figure 10 shows the geometry that we have been discussing in a heuristic three-dimensional projection from a notional N dimensional phase space. The base plane shows the well-known 2D phase portrait of a one-degree-of-freedom system generated as a saddle and a node approach one another to give a saddle-node fold (or limit point). For easy visualization, the trajectory heading towards the node in this plane has been given a much higher damping level than we are currently discussing. The single axis normal to the base plane has to represent all the other phase dimensions. The first important sub-space of this third axis is the CSM which is illustrated as a green transparent upright surface (dimensionality, N โ Figure 10: Geometric arrangement of thresholds and buckling trajectories in the full phase space (including position and velocities). Somewhere between these two typical trajectories lies the special red trajectory which lies precisely in the green CSM. Initial conditions behind the surface generate immediate buckling, while those in front do not. However, the purple trajectory will make another round trip around the node. As the (green) CSM folds strongly further as it extends away from the saddle, after its next round trip the purple trajectory may lie behind those further folds of the green manifold, leading to multiple thresholds. The final divergence of the trajectories is intimately related to the second important sub-space, the so-called centre manifold drawn in orange. This centre manifold is a subset of the green CSM, but has been expanded in the 3D projection of figure 10. When there is no damping this multi-dimensional manifold (dimension N โ
2) contains a variety of periodic and quasi-periodic orbits. On the introduction of damping the orbits inside this manifold drift slowly downwards towards the saddle equilibrium point. All trajectories close to the threshold (the green CSM) are caught up in these circling motions before they are thrown off in opposing directions, including the blue, red and purple examples in the sketch. Since dynamic simulations create a small amount of damping, a simple criterion for the different sets is the long-time limit of the pole deflection ๐ค (๐ก) . Let us denote by ๐ค (๐) the value of the pole deflection ๐ค of the saddle equilibrium at pressure ๐ (transparent surface in left panel of figure 6). Then, (noting the pole deflections of interest are negative) a trajectory after pressure ramp to ๐ buckles, if ๐ค (๐ก) < ๐ค (๐) โ โ for large times ๐ก . 2. avoids buckling, if ๐ค (๐ก) โ ๐ค (๐) converges to a positive value for large times ๐ก (namely ๐ค (๐) โ ๐ค (๐) , where ๐ค (๐) is the pole deflection of the unbuckled stable equilibrium at pressure ๐ ), 3. is on the threshold (in the CSM ), if ๐ค (๐ก) โ ๐ค (๐) goes to zero for large times ๐ก . With the small (and physically necessary) damping created by the dynamic simulations, the convergence to ๐ค (๐) will be slow for threshold trajectories (being slower for smaller damping). After an initial exponential approach to the CSM, damping will introduce a drift to the saddle equilibrium, which is the point of lowest energy in its own CSM. For zero damping, the CSM contains periodic and quasiperiodic orbits. These orbits are all themselves of saddle type (thus unstable), and, hence, not visible in dynamic simulations. (b) Thresholds as connections to the centre-stable set of the saddle Despite the slow convergence, the above distinction provides a simple criterion for establishing thresholds more accurately than shown in figure 6. A pressure ramp to low ๐ leads to a trajectory of the non-buckling type, while for ramps to ๐ = ๐ ๐ the trajectory will buckle. Thus, for fixed end time ๐ก ๐ธ we may apply a bisection in ramp pressures ๐ to find a pressure ๐ that leads to a trajectory that has ๐ค (๐ก ๐ธ ) = ๐ค (๐) . Figure 11. (Left, top) Time profiles ๐ค (๐ก) for different pressure ramps, near ๐ = 0.5104๐ ๐ (green, purple) and near ๐ = 0.5112๐ ๐ (blue, red). (Left, bottom) distance of same trajectories from saddle as function of time. (Right) Same trajectories in (๐ค (๐ก), ๐คฬ (๐ก)) plane. Other parameters: ๐ฟ/โ = 0.25 , time stepsize ๐ฅ ๐ก = 0.087๐ . Figure 11 shows the bracketing trajectories for the result of the bisection for ๐ก ๐ธ = 4๐ (blue and red), and ๐ก ๐ธ = 6๐ (green and purple), for imperfection ๐ฟ = โ/4 and time step ๐ฅ ๐ก =0.087๐ . As the top left panel of figure 11 shows, the pole deflection ๐ค (๐ก) performs oscillations around the saddle value ๐ค (๐) (grey horizontal line) for considerable time before ๐ก ๐ธ (larger than ๐ ). During this time the trajectory is close to the saddle (as the bottom left panel shows). Hence, we can draw a first conclusion that the saddle computed for figures 1 and 3 indeed plays a key role in the buckling. However, the buckling threshold is not given by a trajectory that connects to the saddle, but rather a trajectory that connects into the CSM of the saddle (into a small amplitude periodic or quasi-periodic motion near the saddle). Panels on the right of figure 11 show the threshold-bracketing trajectories in their projection to the (๐ค (๐ก), ๐คฬ (๐ก)) plane. This projection also shows how the threshold trajectories make a small number of excursions where ๐ค (๐ก) is between ๐ค (๐) and before reaching the CSM. The difference between the two bracketing pairs is the time it takes before reaching the CSM. The red/blue pair brackets the threshold trajectory reaching the CSM before (during the initial time up to ), while the green/purple pair reaches the CSM only after time . Since the threshold pressures used in figure 11 are close to each other, the threshold trajectory for ๐ก ๐ธ =6๐ is nearly identical to the threshold trajectory for ๐ก ๐ธ = 4๐ . All trajectories shown in figure 11 only diverge from each other while spending time near the saddle as small oscillations that are part of the CSM: see the near-periodic orbits in the projected phase portraits on the right panels of figure 11, and how previously nearly identical trajectories diverge from these small oscillations. The diverging trajectories follow opposite directions in the unstable dimension (the outset [22]) of the small amplitude oscillation in the centre manifold. This can be seen by comparing the end pieces of the red versus the blue (left (๐ค (๐ก), ๐คฬ (๐ก)) projection in figure 11), or green versus purple trajectories (right (๐ค (๐ก), ๐คฬ (๐ก)) projection). (c) Time ordering of buckling threshold trajectories and pressures From these observations we expect that there is a discrete sequence of buckling thresholds when considering a range of step load pressures ๐ . The sequence is ordered by the time it takes for the threshold trajectory to get close to the CSM of the saddle. This order is not necessarily the same order as in the pressures ๐ . For example, between the two thresholds pressures shown in figure 11 there may be more threshold pressures, for which the trajectory reaches the passive set much later in time (especially for small damping). Thresholds that do not take a long time (such as the threshold given by the blue and purple trajectories in figure 11 for g t T = depend only moderately on the damping (that is, they have a well-defined limit for zero damping). However, the number of thresholds increases as the damping goes to zero, as additional thresholds may occur later and later in time. A rough estimate how many additional thresholds to expect for a particular damping level can be obtained by observing the amplitude and energy level of the small oscillations in the CSM that the first threshold trajectories converge. In figure 11 the small amplitude oscillations for the second (green/purple) pair near the saddle are already much smaller than the oscillations of the first (red/blue) pair. Thus, we expect at most one more threshold occurring after the two observed in figure 11 (in time ordering).
7. Conclusions
Accurate calculations for the buckling of elastic spherical shells under step pressure loading, based on small-strain/moderate rotation theory, have revealed nonlinear features of the dynamic buckling of imperfection-sensitive structures that seem not to have emerged in earlier studies. The most notable is the fact that there is not necessarily a clear threshold between pressure levels that cause buckling and those that do not result in buckling. Instead, particularly for a shell with a relatively small imperfection, there exists a cascade of buckling thresholds. The cascade of thresholds is sensitive to structural damping. For the spherical shell, and probably for other imperfection-sensitive shell structures as well, it appears that, the smaller the damping, the smaller the lowest pressure at which buckling occurs. For the spherical shells with the realistic damping levels employed in this paper, the lowest step buckling pressures were reduced by about 30% below the corresponding static buckling pressures for shells with relatively small imperfections (c.f., figure 8). For shells with larger imperfections, which buckle statically below about 60% of the classical buckling pressure of the perfect shell, the lowest step buckling pressures are reduced by less than 10% below the corresponding static buckling pressures. These dynamic step buckling trends for the spherical shell differ significantly from trends predicted using simple 1-DOF imperfection-sensitive models. The simple models suggest than nearly perfect structures will buckle under step loads only slightly below the corresponding static buckling load, and that the ratio of the step buckling load to the corresponding static load increases as the imperfection increases. We have also found that the lower bound (the astatic pressure) on step buckling pressure for the spherical shell based on overcoming the energy barrier associated with the saddle of the energy landscape lies far below the computed step buckling pressure, especially for shells with small imperfections. By contrast, the lowest step buckling load of the simple 1-DOF models coincides with the astatic load. Since the damping in our dynamic simulations was non-zero and the computed lowest step buckling pressure of the spherical shells depends on damping, it remains an open question as to whether shells with no damping might ultimately after long periods of dynamic oscillation buckle at pressures just above the astatic pressure. The oscillatory interaction between the so-called breathing mode and the modes contributing to buckling, first investigated for ring buckling in [18], appears to be ubiquitous. In some of the buckling literature this type of interaction is referred to as โparametric resonanceโ [7]. For spherical shells buckling at the lowest step pressures, this interaction amplifies the modes contributing to buckling interactions to the point where snap buckling takes over. At step pressures sufficiently above the lowest buckling threshold, snap buckling can occur almost immediately without the preliminary oscillatory interactions. Because the lowest buckling threshold depends strongly on damping, the damping in the simulations should be calibrated to the particular experimental situation studied. Our simulations suggest that most energy loss occurs at the breathing frequency such that matching the damping to observations at the breathing frequency is more important than the particular damping model. The lesson from the present study is that damping is an important consideration in the determination of the lowest step buckling threshold, because lowering the damping level adds additional thresholds that cause buckling with larger delay after the pressure step.
Funding Statement
J. Sieberโs research was supported by funding from the European Unionโs Horizon 2020 research and innovation programme under Grant Agreement number 643073, by the EPSRC Centre for Predictive Modelling in Healthcare (Grant Number EP/N014391/1) and by the EPSRC Fellowship EP/N023544/1.
References [1] Lee, A, Marthelot, J, Jimenez, FL, Hutchinson, JW, Reis, PM. 2016. The geometric role of precisely engineered imperfections on the critical buckling load of spherical elastic shells.
J. Appl. Mech., , pp. 111005-1-11. [2] Hutchinson, JW. 2016. Buckling of Spherical Shells Revisited. Proc. R. Soc. A.
Int. J. Bifur. and Chaos (2), 1630003 (25 pages) (doi: 10.1142/S0218127416300032). [4] Virot, E, Schneider, T, Rubinstein, SM. 2017. Stability landscape in shell buckling. Phys. Rev. Lett. , 224101. (doi: 10.1103/PhysRevLett.119.224101) . [5] Marthelot, J , Lรณpez Jimรฉnez, F , Lee, A , Hutchinson, JW , Reis, PM. 2017. Buckling of pressurized hemispherical shells subject to probing forces.
J. Appl. Mech . , 121005. [6] Hutchinson, JW, Thompson, JMT. 2017. Nonlinear buckling behaviour of spherical shells: barriers and symmetry-breaking dimples. Phil. Trans. R. Soc. A , 20160154 (doi.org/10.1098/rsta.2016.0154) [7] Simitses, GJ. 1990.
Dynamic Stability of Suddenly Loaded Structures.
Springer-Verlag, New York. [8] Sanders, JL. 1963. Nonlinear shell theories for thin shells . Quart. Appl. Math . , 21- 36. [9] Koiter, WT. 1966. On the nonlinear theory of thin elastic shells. Proc. Kon. Ned. Ak. Wet.
B69 , 1-54. [10] Dankowicz, H, Schilder,F. 2013.
Recipes for Continuation . Computer Science and Engineering. SIAM. [11] Thompson, JMT, Hutchinson, JW, Sieber, J. 2017. Probing shells against buckling: A nondestructive technique for laboratory testing.
Int. J. Bifur. Chaos , 1730048. [12] Hutchinson, JW, Thompson, JMT. 2018. Imperfections and energy barriers in shell buckling. Int. J. Sol. Struct. , 152-168. (doi.org/10.1016/j.ijsolstr.2018.01.030). [13] Koiter WT. 1969. The nonlinear buckling behavior of a complete spherical shell under uniform external pressure, Parts I, II, III & IV.
Proc. Kon. Ned. Ak. Wet.
B72 , 40 โ Proceedings of the XI International Congress on Applied Mechanics,
Munich, Germany, edited by Gortler, H, Springer-Verlag, 637-651. [15] Hutchinson, JW, Budiansky, B. 1966. Dynamic buckling estimates.
AIAA Journal , , 525-530. [16] Thompson, JMT. 1966. Dynamic buckling under step loading, Int. Conf. Dynamic Stability of Structures , Northwestern University, Oct. 1965, ed. G. Herrmann, Pergamon Press, Oxford, 215-236. [17] Hoff, NJ, Bruce, VC. 1954. Dynamic analysis of the buckling of laterally loaded flat arches.
J. Math. Phys. , 276 โ J. Appl. Mech. , 259-266. [19] Tamura, YS, Babcock, CD. 1975. Dynamic stability of cylindrical shells under step loading. J. Appl. Mech.
J. Appl. Mech . , 74-80. (doi:10.1115/1.2178836). [21] Bathe, K-J, Hoh, G. 2012. Insight into an implicit time integration scheme for structural dynamics. Comp. Struct. , 1-6. [22]
Thompson, JMT, Stewart, HB. 2002.
Nonlinear Dynamics and Chaos . 2nd ed. Chichester, UK: Wiley. [23] Krauskopf, B, Osinga, HM. 2003. Computing geodesic level sets on global (un) stable manifolds of vector fields.
SIAM J. Appl. Dyn. Sys. โ
69. [24] J. Roorda. 1965. Stability of structures with small imperfections.
J. Engng. Mech. Div. Am. Soc. Civ. Engrs, , 87-106., 87-106.