We study the effect of a dimplelike geometric imperfection on the critical buckling load of spherical elastic shells under pressure loading. This investigation combines precision experiments, finite element modeling, and numerical solutions of a reduced shell theory, all of which are found to be in excellent quantitative agreement. In the experiments, the geometry and magnitude of the defect can be designed and precisely fabricated through a customizable rapid prototyping technique. Our primary focus is on predictively describing the imperfection sensitivity of the shell to provide a quantitative relation between its knockdown factor and the amplitude of the defect. In addition, we find that the buckling pressure becomes independent of the amplitude of the defect beyond a critical value. The level and onset of this plateau are quantified systematically and found to be affected by a single geometric parameter that depends on both the radius-to-thickness ratio of the shell and the angular width of the defect. To the best of our knowledge, this is the first time that experimental results on the knockdown factors of imperfect spherical shells have been accurately predicted, through both finite element modeling and shell theory solutions.

## Introduction

where *E* is the Young's modulus, *ν* is the Poisson's ratio, and $\eta =R/t$ is the dimensionless radius of the shell, of radius *R* and thickness *t*. For decades, this theoretical prediction was found to be in disagreement with the experimental results [5–10], and attempting to reconcile the two has been a cornerstone in structural mechanics [3]. Throughout this paper, we focus exclusively on spherical shells.

In Fig. 1, we plot a survey of historical experimental results from the literature for the knockdown factor $\kappa d=pmax/pc$, defined as the ratio between the experimental buckling load $pmax$ and $pc$, from Eq. (1), as a function of the dimensionless shell radius, *η*. In all these combined experimental studies, the dimensionless shell radius was varied in the range $76\u2264\eta \u22642834$, resulting in a widespread of knockdown factors: $0.05\u2264\kappa d\u22640.9$. The general trend of these data is that $\kappa d$ decreases for increasing values of *η*, albeit with a broad spread. Low-precision spherical shells produced by metal spinning [6,7] or plastic vacuum drawing [8] were found to buckle at relatively low values of the critical pressure ($0.17<\kappa d<0.8$), compared to the classic prediction of Eq. (1), presumably due to significant material and geometric imperfections imparted through the fabrication process. By contrast, high-precision shells fabricated by machining aluminum [9] tended to attain higher buckling pressures ($0.45<\kappa d<0.9$), but still with considerable scatter. Note that most of these experimental investigations were conducted with shallow spherical shell caps. Complete spherical shells were fabricated by electroforming [10,11], with a quality of surface finish that could be systematically improved through a chemical polishing treatment, thereby increasing the knockdown factor from $\kappa d=0.05$ to 0.86. Combined, all these findings have led to the now well-established recognition that the critical buckling load of a shell structure is highly sensitive to imperfections.

Von Kármán et al. [12–14] offered an explanation for the large discrepancies between theory and experiments by finding equilibrium states of the shell involving large deflections that can be maintained by a much lower applied load than $pc$, thereby proposing that the knockdown factors were connected with the elastic postbuckling behavior of shells. Subsequently, Tsien [5] assumed the existence of arbitrary disturbances and attributed the knockdown factors to the highly unstable postbuckling behavior of the shells, and compared his theory against the experiments.

In 1945, Koiter [15] made a groundbreaking contribution to the field by developing the general theory of stability for elastic systems subject to conservative loading. In this work, he introduced an asymptotic method to connect the initial postbuckling behavior with the sensitivity to imperfections. Following the English translation (from Dutch) in 1967 of Koiter's seminal work, there was an upsurge of research on the imperfection sensitivity of the buckling of slender structures, and his general theory was applied to a variety of shell structures and loading conditions [16]. In these investigations, the discrepancies between theory and experiments were attributed to variations of the shell thickness, nonuniformity of loading [17], boundary conditions [18], influence of prebuckling deformations [19], and deviations from the perfect shell geometry [20]. Focusing on cylindrical shells, Babcock [1] performed a direct comparison of the effect on buckling between different types of imperfection from these various studies [17–20] and concluded that the most important factor was the presence of geometric imperfections.

As noted by Babcock [1], fundamental experimental research to help advance the understanding of imperfection sensitivity has typically lagged significantly behind theoretical analysis and computation. As a result, the practical design of shell structures tends to be based on classical results, such as Eq. (1), together with empirical corrections [2]. Still, attempts to experimentally validate theories on imperfection sensitivity were done extensively for cylindrical shells [3] and, in fewer cases, for spherical shells [8,21]. In these experiments, the shape of the specimen was precisely measured prior to carrying out the buckling test. However, the experimental fabrication protocols typically impart randomness to the size and shape distributions of shell defects. As such, deterministic relationships have rarely been found between representative imperfection distributions and the load-bearing capacity of the shell. To circumnavigate this, statistical approaches have been developed to study the problem of shell buckling [3] but these have not yet been widely adopted for design purposes due to the lack of high-precision experimental information on the characteristic distributions of the imperfections and knowledge of their influence on buckling.

More recently, a rapid, versatile, and precision fabrication technique has been developed to manufacture thin elastic shells with controlled geometrical and mechanical properties [22]. This technique involves the coating of hemispherical molds with a polymer (elastomer) solution, which upon curing yields an elastic shell of nearly uniform thickness. Elastomeric shells allow for large elastic deformations to occur at operating pressures that are significantly lower than that for metallic shells, which significantly reduces the experimental complexity. While thermoplastic shells can be produced through injection, rotational, and blow molding, these techniques are typically geared for mass production and less suitable to a laboratory research setting, where reproducible, adaptable, and inexpensive prototyping tools are desirable. Still, to make the experimental technique developed in Ref. [22] relevant to the study of imperfection sensitivity, there is a need to adapt it to also fabricate shells that contain precisely designed defects of known geometric properties. Concurrently to these experimental developments, recent computational advances have yielded powerful and accurate numerical tools for large systems of highly nonlinear ordinary differential equations (ODEs) that can be readily ported to solve the shell buckling equations [23]. The time is therefore ripe to readdress the canonical mechanics problem of buckling of imperfect shells, with the goal of developing a predictive framework that relates the geometry of defects and the critical buckling conditions.

Here, we combine experiments and numerical analysis to revisit the buckling of spherical elastic shells under pressure loading, with an emphasis on determining the geometric role of precisely engineered imperfections on the buckling pressure. First, we develop a novel experimental technique to fabricate thin elastomeric shells containing a single “dimplelike” defect of known geometry, and measure their buckling strength under pressure loading (Sec. 2). In Fig. 1, we plot the knockdown factors of our shells ($20\u2264\eta \u2264108$ for nearly perfect shells and *η* = 108 for shells containing a geometric imperfection), on top of other experimental studies from the literature. We find that $\kappa d$ spans a wide range, but in a way that can be controlled, reproduced, and predicted. Finite element method (FEM) simulations (Sec. 3) are used to characterize the shape of these defects and analyze the buckling behavior of our imperfect shells, in excellent agreement with the experiments (Sec. 4). Moreover, a first-order shell theory is specialized to both perfect and imperfect spherical shells, and a set of nonlinear ODEs are derived to describe the mechanical response of our shells and solved numerically (Secs. 5 and 6). Excellent agreement is found across the triangle of experiments, FEM, and ODE simulations for both the critical buckling pressure as a function of the amplitude of the imperfection and the load–deformation behavior (Sec. 7). Finally, we find that beyond a critical defect amplitude, the buckling pressure becomes independent of the amplitude of the defect and quantitatively characterize this plateau (Sec. 8).

Overall, our results demonstrate that small deviations from the spherical geometry result in large reductions in the buckling pressure, in a way that can be accurately predicted by knowing the shape of the defect.

## Experiments

We have performed precision model experiments to investigate how the buckling strength of hemispherical elastic shells, under pressure loading, is affected by a geometric imperfection. In this section, we start by describing the rapid prototyping technique used to fabricate our elastomeric shells containing a well-defined dimplelike defect at their pole. The profile of this dimple defect is then characterized through digital image processing. Finally, we present the experimental apparatus used to pneumatically load the thin shells and measure the critical pressure at which buckling occurs.

### Fabrication of Precisely Imperfect Thin Elastic Shells.

Our thin elastic shells were manufactured by coating a spherical mold with a polymer solution, following a protocol similar to that introduced in a previous study [22], the basis of which is highlighted next. Two variations of this technique enable us to first fabricate flexible (elastic) molds, which are then used to produce thin elastic shells containing a single dimplelike defect.

The hemispherical elastic molds were fabricated by coating the surface of a rigid hemisphere (radius *R* = 24.85 mm, machined out of polyacetal by computer numerical control milling) with a polymer solution of vinylpolysiloxane (VPS, Elite Double 32, Zhermack, Italy), a silicone-based elastomer. VPS was mixed with a ratio of base to curing agent 1:1 in weight, for 10 s at 2000 rpm (clockwise), and then 10 s at 2200 rpm (counterclockwise) using a centrifugal mixer (ARE-310, Thinky USA Inc., Laguna Hills, CA). The VPS solution was poured onto the hemisphere and cured in approximately 20 min at room temperature (20 °C). Upon curing and peeling from the rigid hemisphere, a VPS shell of thickness *t* = 195 *μ*m was produced. Repeating the process multiple times enabled us to systematically increase the thickness of the shell, which once thick enough, itself became the flexible mold employed to fabricate the thin shells used in the experiments. Three different molds were fabricated with thicknesses, *t*_{mold}={585, 975, 1170} *μ*m, by repeating the coating three, five, and six times, respectively. The Young's modulus of cured VPS was measured to be *E* = 1.255 MPa, and its Poisson's ratio was assumed to be *ν* = 0.5.

The actual thin spherical specimens used in the experiments were fabricated following the same protocol described above, but using the thick elastic shells, themselves employed as molds. The VPS solution was poured onto the concave underside of the mold and turned upside down to drain the excess polymer and produce a thin lubrication film. The curing of this liquid film yielded a thin shell with *t* = 230 *μ*m. Note that this value of *t* was slightly higher than that reported above for a single coating step of the mold due to a slightly longer waiting time between the mixing of the polymer and pouring onto the mold [22], to allow sufficient time to prepare the apparatus and indent the pole (more below). The thin shells obtained this way had uncontrollable imperfections that were intrinsic to the fabrication process, for example, systematic variations of the shell thickness (6.6% standard deviation from pole to equator [22]), air bubbles, homogeneity of the polymer mixture, and surface roughness of the mold. Still, these imperfections were overshadowed by the single dimplelike defect that was precisely introduced in the shell fabrication protocol, as is described next.

In Fig. 2, we present a series of photographs, along with corresponding schematic diagrams, that illustrate the fabrication protocol of our *imperfect* thin hemispherical shells, containing a precisely engineered defect. After filling the mold with VPS and draining the excess liquid, the pole of the mold was indented by a flat plate attached to an universal testing machine (5943, Instron, Norwood, MA). We assume that the mold indentation results in the same displacement of the shell pole from its perfect spherical geometry, such that it defines the amplitude, *δ*, of the defect (this is validated in Sec. 3.1 through FEM simulations). To set *δ*, we programed the Instron to move the indentation tip at a constant velocity (0.3 mm/min) until a specific load was detected by a 10 N load cell, corresponding to the targeted amplitude (based on the linear load–displacement relation), and then fixed this position. The defect amplitude *δ* was therefore defined as the distance between the position where the onset of a nonzero load was first detected and the position at which the target load was reached. While holding the indentation constant, the VPS solution cured inside of the deformed mold. Upon curing and peeling from the mold, the final shells had thickness, *t* = 230 ± 11 *μ*m (uncertainty is standard deviation of ten fabricated shells), resulting in a radius-to-thickness ratio of *η* = 108. Moreover, this procedure of deforming the mold through indentation allowed us to produce shells with a single “dimplelike” defect at its pole, whose amplitude could be varied in the range 0 < *δ* (*μ*m) < 542. A localized thicker band (2 mm thickness) at the equator due to the accumulation of excess polymer ensured that the boundary conditions there were clamped.

### Experimental Profile of the Dimplelike Defect.

Whereas the fabrication technique presented above enables us to control the amplitude of the defect, *δ* (through the depth of the indentation), the precise shape of the dimple is self-selected by the elastic properties, and hence the deformation, of the mold. In particular, we are interested in characterizing the defect by the radial deviation from a spherical shape, *w _{I}*, as a function of the zenith angle,

*β*. Experimentally, we have measured this $wI(\beta )$ profile through digital imaging (Nikon D3200 camera, with a Micro-NIKKOR 60 mm lens) and then extracted the shell contour by an edge detection algorithm (image processing toolbox, matlab). A circle was fit to the region away from the pole, where the effect of the indentation is negligible, corresponding to the profile of the defect-free spherical shell. The difference between this circle and the digitized profile defines $wI(\beta )$. Two representative examples of experimental imperfection profiles are provided in Fig. 3(b), for two shells fabricated using molds with $tmold=585$ and 1170

*μ*m, both at the same defect amplitude

*δ*= 207

*μ*m. This profiles exhibit an inward, axisymmetric, and dimplelike deflection at the vicinity of the pole (for $\beta \u2009\u2272\u200920\u2009deg$), beyond which the shell remains spherical [ $wI(\beta \u2009\u2273\u200920\u2009deg)\u22480$ ]. We have also done FEM simulations to corroborate these findings, the details of which will be presented in Sec. 3.1.

### Measuring the Critical Buckling Pressure.

The experimental critical buckling pressure, $pmax$, was measured for each shell using the following procedure. The shell was mounted onto an acrylic plate with a hole at its center and connected to both a syringe pump (NE-1000, New Era Pump Systems, Inc., Farmingdale, NY) and a pressure sensor (MPXV7002, NXP Semiconductors, The Netherlands). The air inside the shell was extracted at the imposed constant flow rate of 0.1 ml/min, while monitoring its internal pressure at an acquisition rate of 1 Hz using the pressure sensor. The internal pressure decreased gradually with time, until a minimum value was reached, at which the shell buckled. The maximum pressure differential between the outside (atmospheric pressure) and the inside of the shell was defined as the critical buckling pressure, $pmax$.

### Experimental Procedure and Range of Parameters.

We proceed by describing the experimental procedure used to measure $pmax$ for a collection of shells containing precisely designed geometric imperfections, of different amplitude and width. First, 60 shells were fabricated following the protocol detailed in Sec. 2.1, using the three elastic molds with *t*_{mold}={585, 975, 1170} *μ*m (to change the width of the defect), and systematically varied the mold indentation depth (to obtain defect amplitudes in the range 0 ≤ *δ* (*μ*m) ≤ 542). Throughout, the radius and thickness of the shell were kept fixed at *R* = 24.85 mm and *t* = 230 *μ*m, such that $\eta =R/t=108$. For each shell, three identical experimental runs were conducted; each experimental data point represents the average of these measurements and its error bars represent the standard deviation, although these are typically smaller than the symbols size (e.g., as in Fig. 4).

## Finite Element Simulations

The FEM simulations were performed using the commercial package abaqus/standard. The model was simplified to be two-dimensional by assuming rotational symmetry. This reduced the computational cost by a factor of $\u224820$, compared to an equivalent model using a three-dimensional description of the structure using shell elements. The shells were treated as incompressible neo-Hookean solids, and reduced hybrid axisymmetric elements CAX4RH were employed. A convergence study was performed, which led to the selection of a regular mesh with 1000 elements in the zenith direction and an equivalent mesh size in the radial direction (between 6 and 30 elements, depending on the shell thickness). All analyses considered a nonlinear geometry.

Two different sets of FEM simulations were performed for the following purposes: (i) to characterize the shape of the shells obtained through the fabrication process and (ii) to calculate the buckling load and postbuckling response of the shells under external pressure, for shells with a variety of defect geometries.

### FEM of the Profile of the Imperfect Shells.

The goal of this first set of FEM simulations was to model the fabrication procedure and determine the shape of the engineered defect, for different levels of indentation of the flexible molds. Each mold was modeled as a flexible shell (thicknesses *t*_{mold}={585, 975, 1170} *μ*m), and the indentation plate was modeled as a rigid surface using RAX2 elements. A frictionless general contact was defined between all free surfaces. The indentation loading was modeled by imposing the vertical displacement of the plate, which resulted in the deformation of the mold. At the end of the simulation, the position of the inner surface of the mold was extracted and assumed to be equal to the outer surface of the fabricated shell. The defect is defined as the radial displacement *w _{I}* as a function of the zenith angle,

*β*. The amplitude of the defect,

*δ*, is equal to the deflection at the pole, $wI(0)$.

Our simulations show that the width of the defect, defined as the zenith angle at which the deflection *w _{I}* becomes negligible, increases with both the thickness of the mold and the amplitude,

*δ*. Figure 3(a) shows the profiles of shells with

*t*

_{mold}= 585

*μ*m and 30 ≤

*δ*(

*μ*m) ≤ 300. The defect is highly localized near the pole (

*β*= 0), and the small variation of the profile of the shell for increased values of

*δ*can be seen in the zoomed inset of Fig. 3(a).

The shape of different defects can be more easily compared when $wI(\beta )$ is normalized by *δ*. In Fig. 3(b), we compare the defect profiles obtained from FEM and experiments (see Sec. 2.2), finding excellent agreement. The results used in this comparison correspond to shells with the same amplitude, *δ* = 207 *μ*m, fabricated using two molds of thickness, $tmold=585$ and 1170 *μ*m. The clear difference between the profiles obtained with both molds demonstrates that the overall shape of the defect (e.g., its width) can be controlled by varying the thickness of the mold.

Given the good agreement between FEM and experiments, for the remainder of this paper, the reported defect amplitudes and the corresponding profiles will be computed from FEM from the corresponding experimental parameters, given the laborious procedure that would be required to systematically extract the same quantities from the experiments.

### FEM of the Imperfection Sensitivity.

A second set of simulations was then performed to investigate the effect of the geometry of the imperfections on the buckling load of depressurized shells. In this case, the loading was modeled as live pressure, applied on the outer surface of the shells. We found that using the BUCKLE analysis in abaqus significantly overestimates the buckling pressure, even with an existent defect. The reason is that this is a linearized buckling analysis that does not take into account the deformation that takes place in the principal solution, prior to the instability. In order to account for the nonlinear geometry, and given that the collapse of the shells is unstable [24], the simulations employed the Riks method [25] to simultaneously solve for loads and displacements, with the progress of the analysis measured by the arc length of the load–displacement. The buckling pressure was then defined as the maximum pressure attained in the analysis.

*t*= 230

*μ*m, and the geometric imperfections were directly introduced in the mesh. Two approaches to define the shape of the defect were followed. First, for direct comparison with the experiments, the profile of the shells was directly taken from the complete modeling of the full fabrication process (detailed in Sec. 3.1). In this set of simulations, the geometry of the defects changed for every value of the thickness of the mold and the applied indentation. The results from these simulations are shown and compared with the experimental results in Figs. 3 and 4. Second, in order to more thoroughly decouple the effect of the amplitude and the width of the defect, we chose the simpler defect profile of a Gaussian dimple

where *β*_{0} controls the width of the defect. This simple parameterization allowed us to perform a systematic study of the effect of the dimple geometry on the buckling pressure of the shells, presented in Secs. 7 and 8.

where *α* is the edge angle of a shallow spherical shell measured from the axis of symmetry. Kaplan and Fung [6] showed that the nonlinear buckling behavior of a shallow spherical shell is set by *λ*, and subsequent studies have tended to present the buckling pressure as a function of this geometric quantity [6–9]. In the results presented in Sec. 8, we will use a definition of *λ* that replaces *α* by the angular width of the imperfection *β*_{0} from Eq. (2), thereby assuming that the nonlinear deformation occurs only in the region of the shell containing the dimplelike defect. This is similar to the approach followed in the classic numerical analysis of Koga and Hoff [26].

## Comparison Between Experiments and FEM

We now follow the methodologies presented in Secs. 2.4 and 3.2 to compare the experimental and FEM results. In Fig. 4, we plot the knockdown factor, $\kappa d=pmax/pc$ (normalized critical buckling pressure), as a function of the dimensionless amplitude, $\delta \xaf=\delta /t$ (normalized by the shell thickness), of a single dimplelike defect. Three datasets are presented for shells fabricated from molds with *t*_{mold}={585, 975, 1170} *μ*m, resulting in defects of increasingly larger angular width, as characterized in Secs. 2.2 and 3.1.

Focusing first on the experiments, for a shell without an engineered defect ($\delta \xaf=0$), we find a knockdown factor of $\kappa d=0.69\xb10.06$, due to the uncontrollable imperfections that are intrinsic to the fabrication and experimental procedures. These include variations in the shell thickness from the pole to the equator [22], small air bubbles trapped in the elastomer during fabrication, and self-weight, all of which are not taken into account in the classic prediction of Eq. (1). With the presence of a defect, beyond $\delta \xaf>0$, the knockdown factor varies widely in the range $0.15<\kappa d<0.75$, but in a way that can be robustly and reproducibly set by systematically varying the geometry of the defect. The $\kappa d(\delta \xaf)$ data first decrease sharply for $0<\delta \xaf\u2009\u2272\u20091.5$, but, eventually, reach a plateau when $\delta \xaf\u2009\u2273\u20091.5$ at $\kappa d\u223c0.2$. For $\delta \xaf\u2009\u2272\u20091.5$, shells with wider defects (e.g., obtained by using molds with *t*_{mold} = 1170 *μ*m) have knockdown factors that are slightly higher than narrower defects (e.g., *t*_{mold} = 585 *μ*m), but this trend is inverted for $\delta \xaf\u2009\u2273\u20091.5$, even if the differences between the three datasets are relatively small.

The experimental results presented above corroborate the seminal numerical predictions by both Hutchinson [27] for defect shaped with the critical buckling mode at onset, and by Koga and Hoff [26] for axisymmetric dimplelike defects. Note, however, that the defect shape considered by Hutchinson was different from ours, and Koga and Hoff overestimated the effect of the dimpled defects [23]. Moreover, the maximum defect amplitude considered by both of these previous studies was $\delta \xaf=0.75$ [27] and 0.5 [26], such that they did not observe the development of the plateau, whereas we were able to fabricate shells up to $\delta \xaf=2.36$.

In Fig. 4, we superpose numerical FEM results onto the experiments, for identical parameter values, and find remarkable quantitative agreement. Specifically, the FEM data show the presence of a clear plateau at high values of $\delta \xaf$, as well as the crossing and subsequent inversion in the relative buckling strength for shells with different angular widths, when $\delta \xaf\u2009\u2273\u20091.5$. For the parameters explored, the level of this plateau lies in the range $0.17<pplateau/pc<0.20$, such that the buckling pressure has a lower bound at these values. In Sec. 8, we will further explore the FEM simulations to systematically quantify the level and onset (in $\delta \xaf$) of the plateau, as functions of the defect geometry.

To the best of our knowledge, this is the first time that experimental results are reported showing a direct relationship between the critical buckling pressure of spherical shells and the systematically varied geometric properties of an imperfection. Moreover, for a given defect geometry, we are able to accurately predict the associated knockdown factors through FEM. Our results are in stark contrast to the broad spread in the experimental data extracted from the literature shown in Fig. 1, as well as the inability for the classic theories, e.g., Eq. (1), to predict them.

We proceed by supporting this comparison between FEM and experiments with an analytical description based on a first-order shell theory. Specializing this theory for shells containing a single dimplelike defect yields a set of nonlinear ODEs that will then be solved numerically and compared directly with FEM (as in Sec. 7).

## Formulation of the Shell Theory

We now formulate shell buckling equations using a small strains and moderate rotations theory [23]. By focusing on the maximum pressure that the shell can support, we shall demonstrate that middle surface strains remain “very small” and rotations remain “moderately small.” In nonlinear shell theory, this translates into middle surface strains *ϵ* satisfying $|\u03f5|\u226a1$ and rotations $\phi $ satisfying $\phi 2\u226a1$. Rotations about the middle surface tangents are the most important, while rotation about the normal to the shell middle surface turns out to be small in the spherical shell buckling problem. Nevertheless, the equations employed accommodate moderate rotations about the normal. Our analysis indicated that there is essentially no difference between dead and live pressure for the behavior of interest in the current study. Accurate equations for first-order shell theory with small strains and moderate rotations were given by Sanders [28], Koiter [29,30], and Budiansky [31]. These are specialized below for initially perfect spherical shells followed by the introduction of initial imperfections (as in Sec. 6) that resemble those of our experimentally fabricated shells (presented in Sec. 2.1).

*r*as the distance from the origin,

*ω*as the circumferential angle, and $\theta =\pi /2\u2212\beta $ as the meridional angle ranging from 0 at the equator to $\pi /2$ at the pole. The radius of the undeformed middle surface of the shell is

*R*. A material point at $(\omega ,\theta ,R)$ on the middle surface of the undeformed shell is located on the deformed shell at

where $(i\omega ,i\theta ,ir)$ are unit vectors normal and tangent to the undeformed middle surface associated with the respective coordinates. For general deflections, the displacements $(u\omega ,u\theta ,w)$ are functions of *ω* and *θ*. For axisymmetric deflections, $u\omega =0$, while the other two displacements are functions only of *θ*.

*w*from the perfect spherical shape are considered with $(u\omega ,u\theta )I=0$, but imperfections in the form of thickness variations or residual stresses will not be investigated. In addition, our attention is limited to

_{I}*axisymmetric imperfections*such that

*w*is a function of

_{I}*θ*, but not of

*ω*. Assuming that

*w*itself produces small middle surface strains and moderate rotations (a condition easily met in all our examples), $E\alpha \beta I$ denotes the strains in Eq. (7) arising from

_{I}*w*. The total strains due to $(u\omega ,u\theta ,wI+w)$, where

_{I}*w*is additional to

*w*, are denoted by $E\alpha \beta I+U$, and the strains that give rise to stress arising from displacements additional to

_{I}*w*are $E\alpha \beta =E\alpha \beta I+U\u2212E\alpha \beta I$

_{I}Given that the bending strains are linear in the displacements and their gradients, the same process reveals that Eq. (8) still holds for the relationship between the bending strains and the additional displacements, with no influence of *w _{I}*. From here on, the additional displacements $(u\omega ,u\theta ,w)$ will simply be referred to as “the displacements.” An imperfection contribution also arises for live pressure loading which will be introduced shortly.

*S*denoting the spherical reference surface specified by

*r*=

*R*and the Euler angles $(\omega ,\theta )$, the strain energy in the shell is

*p*on the shell is the negative of the work done by the pressure. For

*dead pressure*(per unit original area of the middle surface and acting radially), we have

*live pressure*(per area of the deformed middle surface and acting normal to the deformed middle surface), the potential energy is the negative of the pressure times the change of volume $\Delta V$ within the middle surface. The exact expression for $\Delta V$ is a cubic function of the displacements and their gradients when expressed as an integral over the reference spherical hemispherical surface [23]. It is worth recording that Koiter [30] has given an expression for $\Delta V$ which appears to include errors or misprints. For axisymmetric displacements and

*live pressure*

This result is applicable to either a full spherical shell or any shell segment, such as the hemisphere considered here, which is constrained, and $u\theta $ vanishes on the boundary.

*w*using the process described above for the strains, where

_{I}*w*becomes additional to

*w*. Because it is linear in

_{I}*w*, PE for dead pressure remains as Eq. (12). For live pressure, the process using Eq. (13) gives a lengthy expression which is abbreviated here as

As an aside, it is worth noting that Donnell–Mushtari–Vlasov (DMV) shell theory also generates accurate solutions for the problems considered in this paper. The equations for that theory for a spherical shell are immediately obtained as follows: one omits $\phi r2$ in the strains in Eq. (9), and the in-plane displacements, $u\omega $ and $u\theta $, are also neglected in the rotations in Eq. (6). In addition, the expression for dead pressure given by Eq. (12) is usually assumed for DMV.

## Axisymmetric Deformations of Clamped Hemispheres Containing Axisymmetric Imperfections: A Set of Nonlinear ODEs

*w*, and

*w*are functions of

_{I}*θ*and $u\omega =0$. Hemispherical shells $(0\u2264\theta \u2264\pi /2)$ clamped at the equator and subject to uniform inward pressure

*p*are considered. Dimensionless displacements are defined as $U=u\theta /R,\u2009W=w/R$ and $WI=wI/R$. Let $d()/d\theta =()\u2032$. Then, with

where $(n\u0302\omega ,n\u0302\theta )=[\alpha \u0302/(Et)]cos\u2009\theta (N\omega \omega ,N\theta \theta ),\u2009(m\xaf\omega ,m\xaf\theta )=(R/D)cos\u2009\theta (M\omega \omega ,M\theta \theta ),\u2009p\u0302=(R3/D)cos\u2009\theta p$, and $\alpha \u0302=12(R/t)2$. There are additional terms in these equations for live pressure multiplying $p\u0302$ which have not been shown. The clamping boundary conditions at the equator require $U(0)=W(0)=W\u2032(0)=0$.

*U*and

*W*or, equivalently, in terms of $\phi $ and

*W*, with $U=W\u2032+\phi $. The most highly differentiated terms are $\phi \u2034$ and $W\u2034$, thereby yielding a sixth-order, nonlinear ODE system. In all the problems considered in this paper, the axisymmetric behavior is such that the inward deflection at the pole, $\u2212W(\pi /2)$, increases monotonically, while the dimensionless pressure, $p\u0303=R3p/D$, increases in the early stages and then usually attains a limit point after which it decreases. For this reason, it is effective to treat $p\u0303$ as an unknown, to introduce an extra ODE, $dp\u0303/d\theta =0$, and to prescribe $\u2212W(\pi /2)$ as the “load parameter.” This augmented system can be reduced to seven first-order ODEs in the standard form

with $m\xaf\omega =\u2212sin\u2009\theta \phi +\nu \u2009cos\u2009\theta \phi \u2032,\u2009m\xaf\omega \u2032=\nu \u2009cos\u2009\theta \phi \u2033\u2212(1+\nu )sin\u2009\theta \phi \u2032\u2212cos\u2009\theta \phi ,\u2009m\xaf\theta \u2032=cos\u2009\theta (\phi \u2033\u2212\nu \phi )\u2212(1+\nu )sin\u2009\theta \phi \u2032,\u2009n\u0302\omega =\alpha \u0302\u2009cos\u2009\theta (E\omega \omega +\nu E\theta \theta )$, and $n\u0302\theta =\alpha \u0302\u2009cos\u2009\theta (E\theta \theta +\nu E\omega \omega )$, where $E\omega \omega $ and $E\theta \theta $ are given by Eqs. (16) and (17) using $U=\phi +W\u2032$. The derivative, $n\u0302\u2032\theta $, is directly computed in terms of $\phi $, *W*, and their derivatives.

*W*= 0. The functions $\phi $ and

*W*are analytic at the pole, with $\phi $ being odd and

*W*even about the pole such that $\phi =\phi \u2033=W\u2032=W\u2034=0$ at $\theta =\pi /2$. At the pole, $f2=0,\u2009f3=\phi \u2032,\u2009f4=0,\u2009f5=W\u2033,\u2009f6=0$, and $f7=0$. A somewhat lengthy expansion about the pole provides the following expression for $\phi \u2034$ at $\theta =\pi /2$:

Solving Eq. (20) using a modern nonlinear ODE solver provides highly accurate results. In particular, the buckling pressure, i.e., the maximum pressure attained at the limit point, can be calculated accurately and efficiently. We have used the ODE solver routine DBVPFD in IMSL [32], which incorporates Newton iteration to satisfy the nonlinear equations and an automatic mesh refinement to meet accuracy tolerances. As already noted, the inward pole deflection serves as the loading parameter, and it is increased in incremental steps using a converged solution at one step as the starting guess for the next step. The solution process is fast and robust. As will be illustrated, the solutions can be readily obtained at deflections well past the limit point, beyond the onset of buckling. For the problems that we shall consider, our simulations have shown that there is virtually no difference between predictions for dead and live pressure. The results reported throughout have been computed assuming live pressure.

Thus far, we have exclusively considered axisymmetric imperfections, and both the FEM and ODE analyses assume axisymmetry. It is conceivable that nonaxisymmetric bifurcations could occur for this system. Nevertheless, a recent analysis [23] employing the moderate rotation theory presented in Sec. 5 found no evidence for such bifurcations, neither for perfect shells (even for large pole deflections up to $w/t=10$) nor for shells containing a dimple imperfection (before the maximum pressure of the axisymmetric state). On the other hand, the previous experimental and FEM studies with thin elastic shells under a variety of loading conditions [33–36] have found that an originally axisymmetric buckled shell may develop nonaxisymmetric localized angular structures, referred to as s-cones, in the advanced postbuckling regime. This mechanical behavior is, however, beyond the scope of our work, and the remainder of this paper focuses entirely on an axisymmetric analysis and response. We therefore leave a more detailed investigation of the possibility of nonaxisymmetric imperfections and/or response for a future study.

## Comparison Between ODE and FEM Results

We proceed by comparing the results for the mechanical response of the shells obtained from both the ODE solution and FEM, which also serves as a joint verification of the two frameworks. Figure 5 shows the effect that imperfections with different amplitudes have on the prebuckling and the postbuckling behavior, with a focus on the evolution of the internal pressure during deformation. For the remainder of this paper, no results are provided for very small values of defect amplitude: $\delta \xaf<0.15$ for FEM and $\delta \xaf<0.2$ for ODE. The reason is that such deflections are too small to be accurately and reliably taken into account by the respective numerical algorithms.

In Fig. 5(a), we plot the normalized pressure, $p\xaf=p/pc$, versus the normalized volume change, $V\xaf=\Delta V/V0$, where *V*_{0} is the total volume change of the perfect shell immediately prior to the onset of buckling, for both the ODE solution (solid lines) and FEM (dashed lines). For a nearly perfect shell with a small imperfection $\delta \xaf=0.03$ (black lines), the pressure first increases linearly with increasing $V\xaf$ and reaches a maximum value, $p\xafmax=pmax/pc$, just before buckling occurs. Past this point, $p\xaf$ decreases with decreasing $V\xaf$, closely following the upward branch, and then turns around to eventually decrease with increasing $V\xaf$. If $V\xaf$ is imposed and increased monotonically, then the shell becomes unstable and undergoes snap-buckling almost immediately after the maximum normalized pressure $p\xafmax$ is attained. If the shell were perfect, there would be a pressure drop when $V\xaf=1$ to the lower branch at $p\xaf\u22480.2$, after which $p\xaf$ would continue to decrease for increasing $V\xaf$. For shells containing defects with higher values of $\delta \xaf$, the volume change required for buckling decreases and the peak pressure is progressively lower. Thus, even though increasing $\delta \xaf$ weakens the shell, buckling is less catastrophic. When the imperfection amplitude is sufficiently large (e.g., $\delta \xaf\u22651$), the pressure decreases smoothly after the maximum value is attained, without a pressure jump. It is important to highlight that in all of these results, the FEM and ODE data (dashed and solid lines, respectively) are nearly indistinguishable, which serves as a joint verification of both numerical approaches.

In Fig. 5(b), the same data presented in Fig. 5(a) are now replotted as a function of the normalized displacement of the shell pole, $w\xaf=w/t$, to obtain the load–displacement behavior. For all curves (different values of the defect amplitude, $\delta \xaf$), $p\xaf$ initially increases sharply with $w\xaf$ in the early stages of deformation (linear response), until a maximum pressure is reached at $w\xaf\u22481$, after which the pressure decreases. With increasing $\delta \xaf$, the value of $p\xafmax$ decreases, and the postbuckling decrease of $p\xaf$ with $w\xaf$ becomes less abrupt. Note that all the $p\xaf(w\xaf)$ curves for the different values of $\delta \xaf$ explored approach one another in the later stages of deformation ($w\xaf>5$). Again, an excellent agreement is found between the ODE solutions (solid lines) and FEM (dashed lines) results, with at most 0.9% relative difference in $p\xafmax$ between the two.

## Parametric Exploration of the Knockdown Factor

Having characterized the load-bearing capacity of the imperfect shells with defect of different amplitudes, $\delta \xaf$, we now return to characterize the knockdown factor. First, we use a single geometric parameter to characterize the imperfect shells, and then focus on the plateau observed for $\delta \xaf\u2009\u2273\u20091.5$ (first reported in Sec. 4). In particular, we focus on the dependence of the level and onset of this plateau on both the angular width of the defect, *β*_{0} (Eq. (2)), and the radius-to-thickness ratio of the shell, $\eta =R/t$. Given the excellent agreement found in Sec. 4 between the experiments and FEM (validation), as well as between the FEM and ODE solutions in Sec. 7 (verification), we center this discussion exclusively on the FEM and ODE results.

### Characterization of the Imperfect Shell by a Single Geometric Parameter.

Following the works of Kaplan and Fung [6] and Koga and Hoff [26], we report our results with respect to the geometric parameter *λ* introduced in Eq. (3), but with $\alpha =\beta 0$, which considers the combined effect of *η* and *β*_{0} for a dimensionless characterization of the defect geometry. We performed FEM simulations and ODE calculations for two sets of shells with *η*={100, 200} containing defects in the range $1\u2264\lambda \u22645\u2009(2.34\u2009deg\u2264\beta 0\u226416.54\u2009deg)$. In Fig. 6, their corresponding knockdown factors are plotted versus the imperfection amplitude, $\delta \xaf$.

In Fig. 6(a), we plot $\kappa d=p\xafmax$ versus $\delta \xaf$ for $1.5\u2264\lambda \u22645$ (in steps of 0.5). For each value of *λ*, the ODE solutions (solid lines) and the FEM (dashed and dotted lines) all collapse onto grouped curves. This indicates that for fixed *λ*, the buckling behavior is not affected by different values of *η*. Moreover, these results demonstrate that the single geometric parameter *λ* characterizing the defects governs the imperfection sensitivity of our imperfect shells. All curves exhibit an initial decrease of $\kappa d$ with $\delta \xaf$, followed by a plateau. As the geometric parameter of the defect is increased, the plateau appears at higher values of $\delta \xaf$, and with a level that decreases monotonically with *λ*. For example, the small defect with $\lambda =1.5$ has an initially sharp decrease of $\kappa d(\delta \xaf)$ and the plateau is attained at $\delta \xaf\u22481$, whereas the largest defect considered (*λ* = 5) exhibits a slower initial decay and the plateau is only reached after $\delta \xaf\u22483$.

In Fig. 6(b), we focus exclusively on FEM and present the same data as in Fig. 6(a), but with a higher density of data in the range $1\u2264\lambda \u22645$, in steps of 0.25. At each value of $\delta \xaf$, there is a critical $\lambda c$ that corresponds to the lowest buckling pressure, which is plotted in Fig. 6(c). The stepwise nature of these data stems from the discrete increase of *λ* in steps of 0.25, and a more continuous curve would have been obtained for a finer variation of this parameter. Koga and Hoff [26] also studied the critical conditions that minimize $\kappa d$, for dimplelike defects with amplitudes in the range $0.1\u2264\delta \xaf\u22640.5$, and reported a value of $\lambda c=4$. By contrast, we observe that $\lambda c$ increases monotonically with $\delta \xaf$, within the range of parameters studied, from $\lambda c=1.875$ at $\delta \xaf=0.15$ up to $\lambda c=5$ at $\delta \xaf=3$. This discrepancy is likely due to the rudimentary (but pioneering) computational tools available at the time of their study, as pointed out by Hutchinson [23].

The empirical description of Eq. (23), together with the data in Fig. 6(c), provides a design guideline for the shape that a defect should have in order for a shell to buckle at the lowest possible pressure. Whereas traditional applications in structural mechanics would typically seek to maximize $\kappa d$, these findings could be useful for the more recent movement of utilizing buckling as a mechanism for functionality [37,38].

### Buckling Plateau for Large Amplitude Defects.

where *ξ* is a threshold whose sensitivity is evaluated by choosing different values, $\xi ={0.005,\u20090.01,\u20090.025,\u20090.05}$. In Fig. 7(a), we plot $\u2329\kappa d\u232aplateau$ versus *λ* and find a monotonic decrease, from $\u2329\kappa d\u232aplateau=0.45$ at *λ* = 1 down to $\u2329\kappa d\u232aplateau=0.15$ when *λ* = 5. The level of the plateau is insensitive to the chosen values of *ξ* (with a variation of at most 0.35% across the four cases).

Figure 7(b) plots the onset of the plateau, $\delta \xafplateau$, as a function of *λ*. For small defects in the range $\lambda <2,\u2009\delta \xafplateau$ is approximately constant, but with a value that depends on the choice of *ξ*. As *λ* is increased, $\delta \xafplateau$ also increases but the curves with different values for the thresholds converge. Overall, we conclude that the plateau starts when the amplitude of the imperfection is at least larger than the shell thickness ($\delta \xafplateau\u2009\u2273\u20091$).

## Conclusions

We have reported results from experiments on the critical buckling load of spherical elastic shells under pressure loading, with an emphasis on how their knockdown factors are affected by an engineered dimplelike imperfection. A fabrication method was developed to produce elastomeric spherical shells containing a single defect with geometric properties that could be accurately controlled and systematically varied. Precision experiments were performed to measure the critical pressure for the onset of buckling for these shells. The experimental results showed a direct relationship between the critical buckling pressure and the geometry of the imperfection (amplitude and angular width). In addition, FEM simulations and ODE numerical analyses were conducted, showing excellent quantitative agreement with each other and with experiments. To the best of our knowledge, this is the first time that experimental results on the knockdown factors of imperfect spherical shells have been accurately predicted.

Our study is well aligned with efforts currently underway by NASA and others interested in large shell structures to replace the old empirical knockdown factors employed in design codes by an approach that (i) first, measures the topographic distributions of imperfections, (ii) then, predicts buckling loads from the measured data, and (iii) finally, determines the corresponding knockdown factors [39,40]. In contrast to a statistical approach that starts from measuring uncontrollable imperfections, here we have precisely and systematically controlled a single imperfection and were able to predict the associated knockdown factors. We also found a buckling plateau for large amplitudes of the imperfection and presented the results of FEM simulations and ODE solutions to characterize it. Both the level of the plateau and its onset are functions of a single geometric parameter set by the angular width of a defect and the radius-to-thickness ratio of the shell. Existing experimental data collected from the literature on the buckling of spherical shells (Fig. 1) provide an indication that the plateau may be connected to the apparent lower limit of the ensemble of historic buckling data. This suggests that replacing the current empirical lower limit curves [2] by a deterministic framework may be a goal worth pursuing.

We hope that our results will instigate a resurgence of interest on the mechanics of thin spherical shells and motivate future explorations on the effect of other types of imperfections on their buckling behavior. Shell buckling, in addition to its canonical status in structural mechanics, continues to be an industrially relevant problem. Furthermore, it is also of interest for the life sciences, such as in the contexts of viruses [41], capsules [42], and pollen grains [43]. This is therefore an area of mechanics research that is as relevant as ever and deserves further attention.

## Acknowledgment

This work was supported by the National Science Foundation (CAREER CMMI-1351449).