SIAM J. MATH. FOOL. \begin{matrix}c \\ ◯ \end{matrix} 1988 Society for Industrial and Applied Mathematics Vol. 1, No. 2, pp. 0–12, February 1988 002

THE COMPUTATION OF ONE-PARAMETER FAMILIES OF BIFURCATING ELASTIC SURFACES

FRANK E. BAGINSKI†

Abstract. We consider the problem of constructing the middle surface of a deformed elastic shell from its first and second fundamental forms, â_{αβ} and b̂_{αβ}. The undeformed shell is a spherical cap of radius R and thickness h with an angular width 2θ_{0} where 0 < θ_{0} < π/2. The cap is subjected to a constant uniform load λ and is simply supported at its edge. We seek to compute the one–parameter families of buckled states which branch from the unbuckled state of the shell. This is accomplished in two steps. First, a finite element method is used to solve the governing shell equations, a pair of fourth–order nonlinear partial differential equations. A solution of this system is a curvature potential w, a stress potential f, and the load λ. Using Liapunov-Schmidt reduction, it can be shown that solutions possessing a variety of symmetries bifurcate from the unbuckled state of the shell. In the work that is presented here, we will numerically continue these local branches. We parametrize solution branches in terms of a pseudo-arc-length parameter ρ (i.e., (λ, f, w) = (λ(ρ), f_{ρ}, w_{ρ})), enabling us to track them around turning points. The second step in our solution process is to solve numerically for the parametrization X̂_{ρ} corresponding to the middle surface of the buckled shell Ŝ_{ρ}. We do so by integrating the partial differential equations of Ŝ_{ρ}. The coefficients in these differential equations involve the first and second fundamental forms of the deformed shell Ŝ_{ρ} which can be computed from (λ(ρ), f_{ρ}, w_{ρ}). A number of bifurcation diagrams corresponding to the first three branch points of a spherical cap of size θ_{0} = 12.85^{∘} are presented. For this example, a secondary bifurcation point was found connecting two distinct nonaxisymmetric solution branches. Computer graphics are used to display images of various buckled surfaces which branch from the unbuckled state of the shell.

Key words. Bifurcation, finite elements, nonlinear elasticity, buckling, shells

AMS(MOS) subject classifications. 35B32, 65N30, 73C50, 73H05, 73K15

1. Introduction. In order to measure the stresses in an elastic shell, a local measure of the deformation is needed. For such a purpose, the deformation of the shell can be described by the stretching or bending of the middle surface of the shell (see [27, p. 56]). The most natural measures of the strains in an elastic shell are the differences in the fundamental forms of the deformed and undeformed middle surfaces (see [24]). However, this means that the actual displacements of the shell are known only implicitly through the equilibrium equations and the compatibility equations for the middle surface of the shell (or some approximate version, see [24]). This is not a serious restriction. In the work that we present here, we will use a result from differential geometry which states that a surface (in our case, the middle surface of the buckled shell) can be determined from its fundamental forms. In particular, this means that we do not need to compute the displacement vector directly.

† Department of Mathematics, George Washington University, Washington, DC 20052. The research described in this paper was supported by grants from the George Washington University, IBM Corporation, and the National Science Foundation (DMS-9005503). Part of this research was performed while the author was a visitor at the Center for Geometry, Analysis, Numerics, and Graphics (G.A.N.G.). G.A.N.G. is supported by the National Science Foundation under grants DMS-9011083, DMS-8802858, and DMS-9101903, and the U.S. Department of Energy under grant DE-FG02-86ER250125 of the Applied Mathematical Science Subprogram of the Office of Energy Research.

The model equations that we employ here are referred to as the John shell equations ([19, 20]). The shell equations of John occupy an intermediate position between shallow shell models (such as Marguerre’s equations [26]) and geometrically exact models (e.g., [1]). These equations are derived under the condition that the strains are small, but without a condition on the size of the displacements of the shell. The John model admits axisymmetric and nonaxisymmetric solutions and does not assume that the initial curvature of the shell is small.

There are two steps to our solution process. In the first step, we solve our model equations, which for a spherical cap are equivalent to a pair of fourth–order nonlinear partial differential equations. The model equations are discussed in Sections 2–3. These equations are solved using a finite element approach (see Section 4). A solution is a pair of potential functions from which we can determine the fundamental forms of the deformed middle surface, â_{αβ} and b̂_{αβ}. The second step in our solution process involves determining the surface Ŝ with â_{αβ} and b̂_{αβ} as its fundamental forms. Since our method of determining Ŝ requires only the knowledge of the first and second fundamental forms of Ŝ, this approach could be easily adapted to other shell models in which â_{αβ} and b̂_{αβ} are available.

For a spherical shell, the John equations can be reformulated in terms of a single equation defined on a certain Hilbert space W. In [5], the author showed that a variety of solutions in W bifurcate from the trivial solution of Eqn. (∗∗) (see Sec. 3). These solutions correspond to buckled states of a spherical cap. The operators in (∗∗) are compact and the nonlinear terms are o(||w||) on bounded λ intervals near the trivial solution w = 0. In the generic cases, the eigenvalues of the linearization of (∗∗) in W have multiplicity one. Thus, the global branching theorem of Rabinowitz [28] allows us to continue these local solution branches. However, this theorem does not give information regarding the finer structure of the solution branches. For example, the theorem does not say what type of branching takes place at a bifurcation point (i.e., critical, supercritical, or transcritical), nor does the theorem predict when a branch will “turn around.” The Liapunov–Schmidt method leads directly to criteria which determine the type of bifurcation that takes place (see [5]). However, the Liapunov–Schmidt method gives only a local description of the solution set near bifurcation points. Using a combination of the topological results of Rabinowitz, the local constructive method of Liapunov–Schmidt, and numerical computations, we are able to develop tools which enable us to give a more complete description of the global structure of the solution branches. In this paper, we will describe numerical methods that enable us to carry out this analysis.

Figure 1a. A spherical cap S Figure 1b. A deformed cap Ŝ with a uniform load p

The physical problem that we study is that of a simply–supported elastic spherical cap (denoted by S, see Figure 1a) which is subjected to a constant uniform load p. In this work, the middle surface of the buckled shell will be denoted by Ŝ and the unknown parametrization will be denoted by X̂ (see Figure 1b). The John model and its relation to other mathematical models are discussed in Section 2. A solution of the John equations can be expressed in terms of λ (a parameter proportional to p), a stress potential f, and a curvature potential w. The fundamental forms of the middle surface of the buckled shell can be computed from (λ, f, w) (see Eqns. (2.3)–(2.4) and (3.1)–(3.2)).

In Section 3, we introduce the weak formulation (Eqns. (∗)) of the John shell equations and discuss its relationship to classical solutions and related Hilbert space solutions. For a spherical cap, the John equations are invariant under the group of orthogonal linear transformations that keep the z–axis fixed. To avoid the problems caused by the equivariance of the model equations, we will consider solutions that are even about some plane ϕ = t^{∗} and choose t^{∗} = 0 for convenience.

In Section 4, we describe the numerical methods that are used to solve the John equations for the potential functions w, f. The goal is to replace the infinite dimensional formulation (∗) by an appropriate finite dimensional problem that we can solve (see Eqn. (4.11)). Our approach is fundamentally a Galerkin approximation. The construction of the finite element space W^{h} is discussed in Appendix B. Two linear eigenvalue problems are solved in Section 4 using our finite element method. These problems can be solved exactly via direct methods and they allow us to evaluate the effectiveness of our numerical method.

The approximation to the John equations leads to a nonlinear eigenvalue problem that is discussed in Section 4. We apply a pseudo–arc–length continuation method (see [9], [21]) to track bifurcation curves. In particular, we approximate solution branches of (∗) by solutions in the form (λ(ρ), f_{ρ}^{h}, w_{ρ}^{h}) where f_{ρ}^{h}, w_{ρ}^{h} ∈ W^{h},

lim_{ρ→0^{+}}(λ(ρ), f_{ρ}^{h}, w_{ρ}^{h}) = (λ_{k}^{m}, 0, 0),

ρ is a pseudo–arc–length parameter and λ_{k}^{m} corresponds to a bifurcation point (see (3.10)–(3.11)). For ρ ≠ 0, the corresponding middle surface of the buckled shell will be denoted by Ŝ_{ρ}.

It was not our intention to give an exhaustive analysis of the behavior of elastic spherical caps under a wide range of physical parameters. Our purpose was to demonstrate how a buckled shell could be constructed using our techniques. This is accomplished by presenting an example which illustrates the essential features of the solution process. For this reason, we consider a shell with physical parameters identical with one used in the experimental work of [25] (the example is referred to as Shell SS-62). Section 4 contains a number of bifurcation diagrams for Shell SS-62 along with other numerical results. We should point out that no single example could demonstrate the variety of solution sets that are possible. In particular, a small change in one of the parameters R, h, or θ_{0} could change the structure of the solution set. However, our techniques can be applied to these shells.

In Section 5, we describe how a “buckled surface” Ŝ_{ρ} is constructed from a solution of the John equations. A standard result in differential geometry states that if certain compatibility equations are satisfied, then a surface can be determined from its first and second fundamental forms. This is a constructive result and Ŝ_{ρ} can be determined by integrating a system of differential equations along the coordinate curves of Ŝ_{ρ} (see Eqns. (5.3)).

In Section 6, we apply the techniques that are described in Sections 2–5 to the spherical cap SS-62. Computer generated images of buckled states of spherical caps are presented. In Appendix A, we define some notation related to the eigenvalues of the Laplacian on a spherical cap.

2. The shell equations. In the following, S will denote the middle surface of a spherical cap of thickness h and radius R. We define

Ω_{θ_{0}} = {(θ, ϕ) | 0 ≤ θ ≤ θ_{0}, 0 ≤ ϕ ≤ 2π}

and

∂Ω_{θ_{0}} = {(θ, ϕ) | (θ, ϕ) ∈ Ω_{θ_{0}}, θ = θ_{0}}.

A point on the middle surface of the cap is given by

X(θ, ϕ) = R(cos ϕ sin θ, sin ϕ sin θ, cos θ) ∈ ℜ^{3}

for (θ, ϕ) ∈ Ω_{θ_{0}}. We will consider shells for which 0 < θ_{0} < π/2.

Various models have been used to study the buckling of spherical shells. Some models assume a priori that the shell is shallow. For a spherical cap, the nonlinear shallow shell model was developed by Marguerre [26]. This model has been used in a number of works, including [17] and [34]. Shallow shell equations, which are essentially perturbations of plate models (such as the von Kármán equations), have also been studied (see, e.g., [29]). Other models (see, e.g., [1], [6]) assume that the buckling is axisymmetric and although these provide accurate mathematical models, they do not admit asymmetric solutions.

In the work that is presented here, we will consider asymmetric as well as axisymmetric deformations of a spherical cap. A model which admits such behavior and is not based on a “shallow shell” assumption is a system of nonlinear partial differential equations called the “lowest order interior shell equations” of John (see [19–20]). These equations are accurate when the strains are small, and their nonlinear structure is rich enough to yield solutions possessing the same variety of symmetries as those which are observed experimentally (see [8], [10]).

Using a variational approach, Koiter derived a related set of shell equations called the Koiter–Sanders equations [24, pp. 34–36]. In the case of zero surface loads, the Koiter–Sanders equations agree completely with the lowest order interior shell equations of John as derived in [19, p. 260], and in the following, these will be referred to as the John shell equations. We will consider the John shell equations in an equivalent form developed by Koiter [24, p. 34–36],

(2.1)  \begin{matrix}ε̃^{ασ}ε̃^{βμ}[γ_{αβ|σμ} + b̃_{αβ}ρ_{σμ} + \frac{1}{2}ρ_{αβ}ρ_{σμ}] = 0, \\ ε̃^{ασ}ε̃^{βμ}ρ_{σμ|β} = 0, \\ m^{αβ}|_{αβ} − (b̃_{αβ} + ρ_{αβ})n^{αβ} = −p, \\ n^{αβ}|_{β} = 0, \end{matrix}

together with the constitutive relations,

(2.2)  \begin{aligned}γ_{αβ} & = (Eh)^{−1}[(1 + ν)n_{αβ} − νã_{αβ}n_{γ}^{γ}], \\ m^{αβ} & = Eh^{3}γ^{2}[(1 − ν)ρ^{αβ} + νã^{αβ}ρ_{γ}^{γ}]. \end{aligned}

The unknowns in (2.1)–(2.2) are the strain tensor,

(2.3)  γ_{αβ} = \frac{1}{2}(â_{αβ} − ã_{αβ}),

and the change of curvature tensor,

(2.4)  ρ_{αβ} = b̂_{αβ} − b̃_{αβ}.

Here, ã_{αβ} and b̃_{αβ} are the first and second fundamental forms of the undeformed middle surface of the shell; ε̃^{αβ} denotes the usual alternating tensor (when R = 1, we write a_{αβ}, b_{αβ}, and ε^{αβ} to denote the corresponding tensors). The tensors â_{αβ} and b̂_{αβ} are the first and second fundamental forms of the deformed middle surface of the shell. In addition, m^{αβ} is the tensor of stress couples, n^{αβ} is the tensor of stress resultants, E is Young’s modulus, ν is Poisson’s ratio, γ^{2} = [12(1 − ν^{2})]^{−1}, and p is the normal load on the shell (p is positive when the pressure is directed inward). The symbol “_{|α}” denotes covariant differentiation with respect to the metric tensor ã_{αβ} (or a_{αβ}). The rules for raising and lowering of tensor indices follow the usual conventions.

The John shell equations have been used in [23] to model the buckling of complete spheres and in [5] to model the buckling of spherical caps. In [2], a system of equations was developed from the John equations by assuming that θ_{0} is small and ignoring terms of order θ_{0}^{2}. This system reduces to the Marguerre equations when simply supported boundary conditions are imposed (see [2, Eqns. (2.31)–(2.32)]). In [2, Theorem 2], we showed that the subcritical portion of the solution branch that bifurcates from the trivial solution at the critical buckling load λ_{c} (the smallest eigenvalue of a related linear problem) does turn around and gain stability. In particular, we showed the existence of nontrivial stable subcritical solutions.

In [5], we showed the existence of buckled states possessing circular, pear–shaped, elliptical, triangular, square–shaped, pentagonal and a variety of other symmetries. The solutions were obtained by applying the Liapunov–Schmidt method to the John shell equations. The results obtained were local in the sense that the solutions are “close” to the unbuckled state of the shell. The work that we present here is a natural extension of these local results. We will use the results of [5] to obtain a good initial approximation to a solution that branches from the unbuckled state of the shell and then numerically track solution branches away from the unbuckled state.

3. Formulation of the weak equations. By introducing a stress potential F and a change of curvature potential W and setting K = 1/R^{2}, we find that for a spherical shell, the change of curvature tensor and the tensor of stress resultants can be written as (see [33])

(3.1)  \begin{aligned}ρ_{αβ} & = W_{|αβ} + KW ã_{αβ}, \\ n^{αβ} & = ε̃^{ασ}ε̃^{βμ}F_{|σμ} + KF ã^{αβ}. \end{aligned}

Applying (3.1) to (2.1)–(2.4), one is lead to a pair of fourth–order nonlinear partial differential equations (see [23]), which after rescaling with

(3.2)  W = hγw̃, F = Eh^{3}γ^{2}f̃,

become

(3.3a)  Δ^{2}f̃+ 2Δf̃− α(Δw̃ + 2w̃) + \frac{1}{2}\{w̃, w̃\} = 0 on Ω_{θ_{0}},
(3.3b)  Δ^{2}w̃ + 2Δw̃ + α(Δf̃+ 2f̃) − {w̃, f̃} = −2αλ on Ω_{θ_{0}},

where

α = \frac{R}{hγ} , λ = \frac{pR^{3}}{2Eh^{3}γ^{2}} ,

{u, v} = (ε^{σκ}ε^{βμ}u|_{κμ}v_{,σ})|_{β} + ∇u · ∇v + uΔv + vΔu + 2uv.

The Laplace-Beltrami operator is given by Δw = a^{αβ}w|_{αβ} and ∇u · ∇v = a^{αβ}u_{,α}v_{,β}. We shall require boundary conditions (corresponding to a simply supported edge) in the form,

(3.4)  w̃ = Δw̃ = 0, f̃ = −λ, Δf̃ = 0 on ∂ Ω_{θ_{0}} .

Equations (3.3) together with the boundary conditions (3.4) will be referred to as the simply supported cap problem (SSC). A classical solution of SSC is a pair of functions, w̃ and f̃, that are four times continuously differentiable on Ω_{θ_{0}}, twice continuously differentiable on Ω̄_{θ_{0}} and satisfy (3.3) and (3.4) pointwise. We will consider solutions in the form,

(3.5)  (f̃, w̃) = (−λ + f, w_{0} + w),

where w_{0} is a constant and λ is as previously defined. Prior to buckling, the cap retains its spherical shape (i.e., w_{0} = w = f = 0). A trivial solution of SSC is a pair

(f̃, w̃) = (−λ, 0).

A trivial solution with λ ≠ 0 corresponds to a spherical shape that is in a state of stress. A nontrivial solution of SSC (corresponding to a buckled state of the shell) is a pair

(f̃, w̃) = (−λ + f, w),

where (f, w) ≠ (0, 0). It follows from (3.3) and (3.4) that (f, w) satisfies

(3.6a)  Δ^{2}f + 2Δf − α(Δw + 2w) + \frac{1}{2}\{w, w\} = 0 on Ω_{θ_{0}},
(3.6b)  Δ^{2}w + 2Δw + α(Δf + 2f) − {w, −λ + f} = 0 on Ω_{θ_{0}},
(3.6c)  w = Δw = 0 on ∂ Ω_{θ_{0}},
(3.6d)  f = Δf = 0 on ∂ Ω_{θ_{0}} .

Remark. Edge conditions of the form (3.6c)–(3.6d) are discussed for a plate model in [32, p. 71] and for spherical caps in [2], [5].

The weak formulation of SSC is obtained in the usual manner. Let ϕ ∈ Λ be a test function that is smooth on Ω_{θ_{0}} and vanishes on ∂ Ω_{θ_{0}}, i.e.,

Λ = C^{∞}(Ω_{θ_{0}}) ∩ {u = 0 on ∂ Ω_{θ_{0}}}.

One could use C^{2}(Ω_{θ_{0}}) in place of C^{∞}(Ω_{θ_{0}}) (see [32] for a further discussion on weak solutions). Multiplying (3.6a)–(3.6b) by ϕ, integrating by parts, and applying the boundary conditions, we are lead to the following set of equations,

(*)  \begin{aligned}ℓ(f; ϕ) + αa(w; ϕ) + \frac{1}{2}b(w, w; ϕ) & = 0, ϕ ∈ Λ, \\ ℓ(w; ϕ) − αa(f; ϕ) − λa(w; ϕ) − b(w, f; ϕ) & = 0, ϕ ∈ Λ. \end{aligned}

Remark. The boundary conditions Δf = Δw = 0 are natural and will be satisfied by a classical solution that is generated by a weak solution.

The operators ℓ(w; ·), a(w; ·) and b(w, f; ·) are defined by

(3.7)  ℓ(w; ϕ) = ∫ (ΔwΔϕ + 2ϕΔw)dS,

(3.8)  a(w; ϕ) = ∫ (∇w · ∇ϕ − 2wϕ)dS,

(3.9)  b(w, f; ϕ) = ∫ {w, f}ϕdS,

where all integrations are carried out over Ω_{θ_{0}} and dS = sin θdθdϕ. We say that (λ, f, w) is a weak solution of SSC if (f, w) satisfies (∗) for all ϕ ∈ Λ. In [5], we showed that the weak formulation of SSC leads to an equivalent Hilbert space formulation. The Hilbert space is denoted by H and is the completion of Λ relative to the norm defined by (3.7), i.e., ||w||^{2} = ℓ(w; w).

Problem (∗) can be reformulated in terms of an operator equation on H given by (see [5])

(**)  L_{λ}w + αQ(w) + C(w) = 0,

where

L_{λ}w = w − λAw + α^{2}A^{2}w,
Q(w) = \frac{1}{2}AB(w, w) + B(Aw, w),
C(w) = \frac{1}{2}B(w, B(w, w)).

The existence of the operators A and B are a consequence of the Reisz Representation Theorem and the definitions

(Aw, ϕ) = a(w; ϕ),
(B(w, f), ϕ) = b(w, f; ϕ).

It follows from their definitions that A : H → H and B : H × H → H.

Equations (3.6a)–(3.6b) and (∗∗) are invariant under the group of orthogonal linear transformations that keep the z–axis fixed. To avoid this problem, we will consider functions which are even about the plane ϕ = 0. A further discussion of this topic can be found in [2] and [5]. The resulting space of functions satisfying this condition is also a Hilbert space and will be denoted by W where

W = {w ∈ H |w = ∑_{m=0}^{∞} ∑_{k=1}^{∞} a_{k}^{m}u_{k}^{m}} .

The u_{k}^{m}’s are eigenfunctions of the Laplacian and the set {u_{k}^{m}} is complete set in W (see Appendix A). It follows from their definitions that A : W → W and B : W × W → W. In the following, we will consider solutions of (∗∗) in W. Notice that if (λ, w) is known, then f can be recovered from the first equation in (∗).

A branch point of the trivial solution is a pair (λ′, 0) such that in every neighborhood of (λ′, 0) there is a nontrivial solution. The branch points of SSC are in the form (λ_{k}^{m}, 0) where

λ_{k}^{m} = μ_{k}^{m} + \frac{α^{2}}{μ_{k}^{m}}

and μ_{k}^{m} is an eigenvalue of the Laplacian on a spherical cap domain (see Appendix A). The generic branch points are those for which α^{2} ≠ μ_{k}^{m}μ_{j}^{n} with (m, k) ≠ (n, j) and θ_{0} ∉ Θ_{0} (Θ_{0} is the set of values of θ_{0} for which μ_{k}^{m}(θ_{0}) = μ_{j}^{n}(θ_{0}) and (m, k) ≠ (n, j), see [3]–[4]).

In the generic cases, two types of solution branches were found. When a bifurcation point is in the form (λ, w) = (λ_{k}^{0}, 0), there is a single transcritical branch,

(3.10)  (λ, w_{k}^{0}(λ)) = (λ_{k}^{0} + αη, ηβ(η)u_{k}^{0} + U(ηβ(η), η)) , 0 < |η| < η^{∗},

where η = (λ − λ_{k}^{0})/α, β(η) and U are analytic in η, and ∥U∥ = o(|η|). When the bifurcation point is in the form (λ, w) = (λ_{k}^{m}, 0) with m ≠ 0, there are two subcritical branches,

(3.11) \begin{aligned}(λ, w_{k}^{m,±}(λ)) = (λ_{k}^{m} & + αη, |η|^{\frac{1}{2}} β^{±}(η)u_{k}^{m} \\ & + |η|v(β^{±}(η)u_{k}^{m}) + V (|η|^{\frac{1}{2}} β^{±}(η), η)), η ∈ I_{0}, m ≠ 0, \end{aligned}

where η = (λ − λ_{k}^{m})/α, β^{±}(η) and V are analytic in η, ∥V ∥ = o(|η|) and v(u) is homogeneous of degree two (v(βu) = β^{2}v(u)); I_{0} is of the form {η | − η^{∗} < η < 0} or {η | 0 < η < η^{∗}}, depending on a certain condition related to the solution of the reduced branching equations (see [5, Theorem 2]). At the critical buckling load, we showed that I_{0} = {η | − η^{∗} < η < 0}. This agrees with our numerical computations presented in Section 4 where the critical buckling load is λ_{c} = λ_{1}^{1}. When λ_{k}^{m} is not the critical buckling load, the condition [5, Eqn. (5.34)] determines the form of I_{0}. We also found a pair of nonaxisymmetric solutions that branch subcritically from (λ_{1}^{2}, 0) where λ_{1}^{2} ≠ λ_{c}.

Solutions in the form (3.11) are asymmetric and can be characterized by their symmetries, (see [5, Section 4]). Solutions in the form (3.10) are axisymmetric. In [5], we referred to the eigenfunction u_{k}^{m} that appears in (3.10) or (3.11) as the fundamental buckling mode of the solution. Solutions possessing more than one fundamental buckling mode were shown to exist in [5], but these are nongeneric cases and will not be considered here.

The results presented in (3.10)–(3.11) are local in the sense that η^{∗} and ∥w∥ are small. However, under the same conditions under which solutions (3.10)–(3.11) are shown to exist, we can apply the global results of Rabinowitz to (∗∗) and extend these local branches (see [28] or [30, Theorems 13.9–13.10]). The key condition, which is satisfied in the generic cases, is that the null spaces corresponding to branch points are one–dimensional. The compactness of the operators and their behavior near w = 0 is discussed in [5]. While [30, Theorem 13.9–13.10] does give some information on the global behavior of solution branches, it does not yield any intermediate information such as when a subcritical solution will turn around. However, we do know that a subcritical solution must either intersect the trivial solution at another branch point or turn around and go off to infinity, for if (λ, w) is a solution, then

(3.12)  λ ≥ μ_{1}^{0} − \frac{α^{2}}{8μ_{1}^{0}}

(see [22]). In any case, if we are interested in obtaining more precise information on the behavior of bifurcating solutions, we must seek numerical solutions.

4. Numerical Solutions of Weak Equations. In this section, we describe the numerical methods that will be used to solve (∗). Our approximations are based on a Galerkin approach. In order to measure the effectiveness of our numerical method, we consider two linear eigenvalue problems that are related to (∗). We estimate the eigenvalues using our finite element approach and compare these results to the true eigenvalues which can be obtained by solving for the roots of certain Legendre functions. The last part of this section presents results related to the solutions of the nonlinear problem (∗). Residuals for the numerical solutions of (∗) along with a number of bifurcation diagrams are presented also.

The solutions that we construct will lie in the finite element space W^{h}, the linear span of the set

(4.1)  B^{h} = {ψ_{j}^{h}(θ, ϕ) | j = 1, 2, …, N}.

The set B^{h} and the vector space W^{h} are described in Appendix B. A member of W^{h} will be denoted by w^{h} where the h is a measure of the discretization, which in our applications is the gap between successive breakpoints (see Appendix B).

We define the operator P^{h} : W → W^{h} as follows,

(4.2)  w^{h} ≡ P^{h}w = ∑_{i=1}^{N} w_{i}ψ_{i}^{h},

where w^{h} ∈ W^{h} and the coefficients w_{i}, i = 1, 2, … N are obtained by solving the linear system,

(4.3)  ℓ(∑_{i=1}^{N} w_{i}ψ_{i}^{h}; ψ_{k}^{h}) = ℓ(w; ψ_{k}^{h}), k = 1, 2, …, N.

Relation (4.3) leads to a linear system,

(4.4)  ( \begin{matrix}L_{11} & L_{12} & … & L_{1N} \\ L_{21} & L_{22} & … & L_{2N} \\ ⋮ & ⋮ & & ⋮ \\ L_{N1} & L_{N2} & … & L_{NN} \end{matrix} ) ( w_{1} ; w_{2} ; ⋮ ; w_{N} ) = ( ℓ(w; ψ_{1}^{h}) ; ℓ(w; ψ_{2}^{h}) ; ⋮ ; ℓ(w; ψ_{N}^{h}) )

where

(4.5)  L_{mn} = ℓ(ψ_{n}^{h}; ψ_{m}^{h}), m, n = 1, 2, …, N.

For convenience and future reference, we define

(4.6)  A_{mn} = a(ψ_{n}^{h}; ψ_{m}^{h}), m, n = 1, 2, …, N.

The bilinear forms a(u; v) and ℓ(u; v) were defined in Section 3. We define also two N × N matrices, L and A by

(4.7)  L = [L_{nm}], A = [A_{nm}]

respectively. For the applications that are presented here, the matrices L and A were verified numerically to be positive definite. Linear Eigenvalue Problems.

Next, we consider two eigenvalue problems, restricting our attention to solutions that are even about ϕ = 0. Solutions to both these problems are related to the eigenvalues of the Laplacian (see (Δ) in Appendix A). In the following, μ_{k}^{m} will denote an eigenvalue of the Laplacian and u_{k}^{m}(θ, ϕ) will denote the corresponding eigenfunction.

We say that μ is an eigenvalue of (P_{1}) if there is a w ≠ 0 and a corresponding μ ≠ 0 which satisfy,

(P_{1})  \begin{aligned}(Δ + 2)(Δ + μ)w & = 0, Ω_{θ_{0}}, \\ Δw = w & = 0, ∂ Ω_{θ_{0}} . \end{aligned}

It is easy to show that the eigenvalues of problems (P_{1}) and (Δ) are identical.

The weak formulation of (P_{1}) is obtained by multiplying (P_{1}) by a test function ϕ ∈ Λ and then integrating over Ω_{θ_{0}}. In particular, we are lead to

ℓ(w; ϕ) = μa(w; ϕ) for all ϕ ∈ Λ.

The corresponding Galerkin approximation is obtained by replacing w ∈ W with w^{h} ∈ W^{h} and replacing Λ with W^{h} (or equivalently B^{h}). We are led to the discrete problem, i.e., find a w^{h} ∈ W^{h} such that

(4.8)  ℓ(w^{h}; ψ_{i}^{h}) = μ^{h}a(w^{h}; ψ_{i}^{h}) for all ψ_{i}^{h} ∈ B^{h} .

In the notation introduced above, we see that μ^{h} is an eigenvalue of (4.8) if and only if μ^{h} is an eigenvalue of the following matrix equation:

(P_{1}^{∗})  Ly⃗ = μ^{h}Ay⃗, for some y⃗ ∈ R^{N}.

The eigenvalues of the Laplacian can be computed by solving for the roots of certain Legendre functions (see (A.4)). Table 1 summarizes the results which are obtained when solutions of (P_{1}^{∗}) are used to approximate solutions of (P_{1}) for a cap with θ_{0} = 12.85^{∘}. Column 2 summarizes the results that were obtained when Bsplines of order k_{θ} = 6 with a uniform breakpoint sequence were used to estimate the eigenvalues of (P_{1}). Column 3 summarizes the results that were obtained when Bsplines of order k_{θ} = 6 with a nonuniform breakpoint sequence were used to estimate the eigenvalues of (P_{1}). The eigenvalue μ_{1}^{3} was not estimated since it was not needed for the examples presented in Section 5.

The results corresponding to nonaxisymmetric eigenvalues using uniform breakpoints are acceptable. However, the results corresponding to axisymmetric eigenvalues using uniform breakpoints are poor. This is most likely due to the large value of the second derivative (calculated with respect to θ) of the finite element nearest the pole θ = 0. By spreading out the breakpoints near the pole, we were able to obtain a better approximation to the axisymmetric eigenfunctions. While the smoothness of the finite elements is not changed by spacing the breakpoints in this way, the net effect is to control the value of the second derivative near θ = 0. The spacing of the breakpoints near the pole is not as crucial when approximating eigenvalues corresponding to nonaxisymmetric eigenfunctions since ∂^{2}u_{k}^{m}/∂θ^{2}(0^{+}, ϕ) = 0 for m ≥ 2. While ∂^{2}u_{k}^{1}/∂θ^{2}(0^{+}, ϕ) ≠ 0, the behavior of u_{k}^{1} is modeled accurately by our choice of finite elements (see Appendix B).

While it was possible to get a better approximation to μ_{k}^{0} by increasing the order of the finite elements, higher order splines do not realistically model the behavior of an elastic shell. For a set of finite elements of order greater than 6, the elements behave stiffly and very little deformation from the initial spherical shape was observed when the nonlinear buckling problem was solved. In summary, the numerics indicate that to obtain an accurate estimate of an eigenvalue of (P_{1}) corresponding to an axisymmetric eigenfunction, one can choose the order of the B-splines to be k_{θ} = 6 with a non-uniform spacing of the breakpoints. Nonaxisymmetric functions can be adequately approximated with a uniform breakpoint sequence.

Following the same approach that was used to estimate the eigenvalues of P_{1}, we can estimate the eigenvalues of the Laplacian for a spherical cap (see (Δ) in Appendix A). The results for a cap with θ_{0} = 12.85 are presented in Table 2. The results presented in Table 2 shows that our method performs well for a second-order problem, regardless of the type of breakpoint sequence that is chosen.

k_{θ} = 6^{1}  k_{θ} = 6^{2}  k_{θ} = 6^{1} k_{θ} = 6^{2}
μ_{k}^{m}  (μ_{k}^{m})^{h}  (μ_{k}^{m})^{h}  μ_{k}^{m}  (μ_{k}^{m})^{h}  (μ_{k}^{m})^{h}
μ_{1}^{0} = 114.64221 176.95396 114.69618  μ_{1}^{0} =114.64221 114.63603 114.64188
μ_{1}^{1} = 291.89360 291.89361 291.89373  μ_{1}^{1} = 291.89360 291.89361 291.98361
μ_{1}^{2} = 525.35963 525.35984 525.36039  μ_{1}^{2} = 525.35963 525.35983 525.35985
μ_{2}^{0} =605.46672 671.46645 605.44741  μ_{2}^{0} = 605.46672 605.26023 605.46498
μ_{1}^{3} = 811.96694 – – μ_{1}^{3} = 811.96694 – –
μ_{2}^{1} = 978.51482 978.51469 978.61548  μ_{2}^{1} = 978.51482 978.51469 978.51836
μ_{1}^{4} = 1149.8349 1149.8354 1149.8648  μ_{1}^{4} = 1149.8349 1149.8354 1149.8364
μ_{2}^{2} = 1409.5762 1409.5759 1409.7771  μ_{2}^{2} = 1409.5762 1409.5759 1409.5879
μ_{3}^{0} = 1488.4974 1566.0310 1488.7010  μ_{3}^{0} = 1488.4974 1487.6575 1488.5185

Table 1. Eigenvalues of (P_{1}) Table 2. Eigenvalues of (Δ)

1 – Using uniformally distributed breakpoints N = 187, θ_{0} = 12.85^{∘}. 2 – Using nonuniformally distributed breakpoints N = 187. θ_{0} = 12.85^{∘}.

Nonlinear problems

We approximate the solution of (∗) in the form (3.10) or (3.11) by (λ^{h}, w^{h}, f^{h}) where

(4.9)  (λ^{h}, w^{h}, f^{h}) = (λ^{h}, ∑_{i=1}^{N} w_{i}ψ_{i}^{h}, ∑_{i=1}^{N} f_{i}ψ_{i}^{h}).

We can choose the basis functions in such a way that a solution in the form (4.9) will be sufficiently smooth to apply the results of Section 5. To simplify the notation, we will drop the superscript notation on ψ_{i}^{h} and write ψ_{i} to denote an element of B^{h}. We will use a superscript ^{h} when referring to an approximate solution (λ^{h}, w^{h}, f^{h}) in the following sections.

Upon substituting (4.9) into the weak equations (∗), we obtain (4.10)

\begin{matrix}∑_{n=1}^{N} ℓ(ψ_{n}; ϕ)f_{n} + α ∑_{n=1}^{N} a(ψ_{n}; ϕ)w_{n} + \frac{1}{2} ∑_{n=1}^{N} ∑_{q=1}^{N} b(ψ_{n}, ψ_{q}; ϕ)w_{n}w_{q} = 0, \\ ∑_{n=1}^{N} ℓ(ψ_{n}; ϕ)w_{n} − α ∑_{n=1}^{N} a(ψ_{n}; ϕ)f_{n} − λ ∑_{n=1}^{N} a(ψ_{n}; ϕ)w_{n} − ∑_{n=1}^{N} ∑_{q=1}^{N} b(ψ_{n}, ψ_{q}; ϕ)w_{n}f_{q} = 0, \end{matrix}

where the unknowns are (f_{1}, …, f_{N}, w_{1}, …, w_{N}). Setting ϕ = ψ_{m}, for m = 1, … N in (4.10), we are lead to the following system of 2N quadratic equations, (4.11)

\begin{aligned}G_{m}(λ, f^{h}, w^{h}) & = ∑_{n=1}^{N} L_{mn}f_{n} + α ∑_{n=1}^{N} A_{mn}w_{n} + \frac{1}{2} ∑_{n=1}^{N} ∑_{q=1}^{N} B_{qnm}w_{n}w_{q} = 0, m = 1, 2, … N, \\ G_{N+m}(λ, f^{h}, w^{h}) & = ∑_{n=1}^{N} L_{mn}w_{n} − α ∑_{n=1}^{N} A_{mn}f_{n} \\ & − λ ∑_{n=1}^{N} A_{mn}w_{n} − ∑_{n=1}^{N} ∑_{q=1}^{N} B_{qnm}w_{n}f_{q} = 0, m = 1, 2, … N, \end{aligned}

where B_{qnm} = b(ψ_{n}, ψ_{q}; ψ_{m}), and A_{mn}, L_{mn} are as previously defined.

In order to compute a solution of the nonlinear system (4.11), we first need to compute an initial guess. In particular, let β^{∗}u_{k}^{m} denote the lowest order term of w_{k}^{0}(λ), w_{k}^{m,+}(λ), or w_{k}^{m,−}(λ) (see (3.10)–(3.11)). Because o(||U||) = o(||V ||) = o(||v||) = o(|η|) near η = 0, a good initial guess for a solution of (4.11) is given by w_{0}^{h}, where

(4.12)  w_{0}^{h} = P^{h}[β^{∗}u_{k}^{m}] = ∑_{i=1}^{N} w_{i}^{0}ψ_{i}.

The coefficients w_{i}^{0} are found by solving

(4.13)  Lw⃗_{0} = b⃗,

where L = [L_{ij}], w⃗_{0} = (w_{1}^{0}, w_{2}^{0}, …, w_{N}^{0})^{t}, b⃗ = (b_{1}, b_{2}, …, b_{N})^{t}, and b_{i} = ℓ(β^{∗}u_{k}^{m}; ψ_{i}), i = 1, …, N. To determine f when η is small, we ignore the higher order terms in (3.6a) − (3.6b) and use the resulting approximation, Δf + αw = 0. Since u_{k}^{m} is an eigenfunction of the Laplacian and (Δ + 2)^{−1} exists on Λ (see [5]), we see that

(4.14)  f = −\frac{αβ^{∗}}{μ_{k}^{m}} u_{k}^{m}.

Hence, for small η ≠ 0, we find that an initial guess for a solution is given by

(4.15)  (λ_{0}^{h}, w_{0}^{h}, f_{0}^{h}) = (λ_{k}^{m} + αη, ∑_{i=1}^{N} w_{i}^{0}ψ_{i}, ∑_{i=1}^{N} f_{i}^{0}ψ_{i}),

where f_{i}^{0} = −(α/μ_{k}^{m})w_{i}^{0}.

Remark. While it is not possible to determine the scalar β^{∗} exactly when η ≠ 0 by only looking at the linear terms in (3.6a)–(3.6b), this was not a serious problem. In practice, our method was sufficiently robust so that for relatively large values of β^{∗}, the iterations converged to the correct solution. If |β^{∗}| was too small, the iterations converged to the zero solution. It was not difficult to find an adequate initial guess w_{0}^{h} for the solution of the nonlinear problem.

In the remaining sections of this paper, our numerical computations will be applied to a spherical cap with specifications based on Shell SS-62 in [25],

R = 3.0 in., E = 10, 800, 000 psi, ν = 0.3, h = 0.03 in., θ_{0} = 12.85^{∘}.

Using these values, we find that γ = 0.3026 and α = 330.5. Our boundary conditions are not the same as those used in the experiments in [25], so it is not possible to compare our results directly.

In [5], we showed that λ corresponds to a branch point if

λ = μ + \frac{α^{2}}{μ}

for some μ = μ_{k}^{m} where μ_{k}^{m} is an eigenvalue of the Laplacian. We define λ_{k}^{m} = μ_{k}^{m} + α^{2}/μ_{k}^{m}. From Figure 2, we see that the smallest λ_{k}^{m} for SS-62 is

λ_{c} = λ_{1}^{1} = μ_{1}^{1} + \frac{α^{2}}{μ_{1}^{1}} = 666.0024,

and the next two bifurcation points occur when

λ_{1}^{2} = μ_{1}^{2} + \frac{α^{2}}{μ_{1}^{2}} = 733.2172, λ_{2}^{0} = μ_{2}^{0} + \frac{α^{2}}{μ_{2}^{0}} = 785.8234.

Figure 2. λ = μ + α^{2}/μ

In order to measure how well w_{0}^{h} approximates β^{∗}u_{k}^{m}, we consider the quantity,

\frac{|ℓ(β^{∗}u_{k}^{m} − w_{0}^{h}; w_{0}^{h})|}{|β^{∗}| · ℓ(w_{0}^{h}; w_{0}^{h})^{1/2}},

for the eigenfunctions corresponding to the first three branch points of SS-62 (see Table 3). We divide by |β^{∗}| in the previous expression, because β^{∗} is initially small in our applications. The eigenfunction u_{k}^{m} was normalized so that a(u_{k}^{m}; u_{k}^{m}) = 1.

μ_{k}^{m}  \frac{|ℓ(β^{∗}u_{k}^{m}−w_{0}^{h};w_{0}^{h})|}{|β^{∗}|·ℓ(w_{0}^{h};w_{0}^{h})^{1/2}}
μ_{1}^{1} = 291.8936 0.5506 × 10^{−7}
μ_{1}^{2} = 525.3596 0.1159 × 10^{−7}
μ_{2}^{0} = 605.4667 0.22369 × 10^{−12}

Table 3. Approximating an eigenfunction The solution (4.15) is used as the initial guess for a Newton’s method solu-

tion of (4.11) with λ = λ_{0}^{h} fixed. The actual method used to solve (4.11) was a Levenberg–Marquardt algorithm, a variation of Newton’s method (see the IMSL subroutine DNEQNJ [18, p. 776]). The solution obtained in this way with initial guess (4.15) will be denoted by (λ_{1}^{h}, w_{1}^{h}, f_{1}^{h}). Note, λ_{0}^{h} = λ_{1}^{h}.

To measure the error in our approximation of the nonlinear problem, we consider the sum of the squares of the residuals of each equation in (4.11). We define

(4.16)  r(λ^{h}, f^{h}, w^{h}) = ∑_{i=1}^{2N} |G_{i}(λ^{h}, f^{h}, w^{h})|^{2}.

Table 4 contains the values of r for a number of solutions that bifurcate from the trivial solution at one of the first three branch points of the Shell SS–62. A negative value of β^{∗} was used in the calculations that are presented in Table 4. The third column of Table 4 gives r(λ_{0}^{h}, f_{0}^{h}, w_{0}^{h}). Note, these values could be lowered by choosing β^{∗} smaller. We include them in Table 4 for completeness. The fifth column of Table 4 gives r(λ_{1}^{h}, f_{1}^{h}, w_{1}^{h}).

In order to compute (λ_{2}^{h}, f_{2}^{h}, w_{2}^{h}), we use (λ_{2}^{h}, f_{1}^{h}, w_{1}^{h}) as an initial guess for the Newton iterations. The solution branch emanating from the trivial solution could be tracked by repeating this process. However, as discussed in Section 1, we know that subcritical branches must eventually turn around. Clearly, this method will break down near these turning points. The results for λ = λ_{1}^{h} and λ = λ_{2}^{h} will be used as the initial starting points for a continuation method.

λ_{k}^{m}  λ_{0}^{h}  |r(λ_{0}^{h}, x_{0}^{h})|  λ_{1}^{h}  |r(λ_{1}^{h}, x_{1}^{h})|  λ_{2}^{h}  |r(λ_{2}^{h}, x_{2}^{h})|
λ_{1}^{1} = 666.002 665.3728 0.118 × 10^{−1} 665.3728 0.59568 × 10^{−6} 665.3687 0.7781 × 10^{−9}
λ_{1}^{2} = 732.585 732.5843 0.167 × 10^{−1} 732.5843 0.42892 × 10^{−7} 732.5825 0.4298 × 10^{−12}
λ_{2}^{0} = 785.8234 785.5071 0.7835 × 10^{−4} 785.5071 0.64638 × 10^{−10} 785.5061 0.2684 × 10^{−21}

Table 4. Residuals of Nonlinear Problems

The dotted lines in Figures 3a–3b are the graphs of w_{0}^{h}(θ, 0) for (m, k) = (1, 1) and (2, 1). Figure 3c presents a graph of w_{0}^{h}(θ, 0) for (m, k) = (0, 2). The relative error,

(4.17)  max_{0≤θ≤θ_{0}} \frac{|u_{k}^{m}(θ, ϕ) − P^{h}u_{k}^{m}(θ, ϕ)|}{|u_{k}^{m}(θ, ϕ)|}, 0 < θ < θ_{0}, ϕ = 0,

was of order 10^{−5} for (m, k) = (0, 2) and of order 10^{−7} for (m, k) = (1, 1) and (2, 1). The solid curves in Figures 3a–3c represent a slice of the solution w_{1}^{h}(θ, 0) for 0 < θ < θ_{0} and (m, k) = (1, 1), (2, 1), (0, 2).

Figure 3a. ±w_{0}^{h}(θ, 0) and ±w_{1}^{h}(θ, 0), 0 < θ < θ_{0} with uniform breakpoints and (m, k) = (1, 1)

Figure 3b. ±w_{0}^{h}(θ, 0)) and ±w_{1}^{h}(θ, 0), 0 < θ < θ_{0} with uniform breakpoints and (m, k) = (2, 1)

Figure 3c. w_{0}^{h}(θ, 0) and w_{1}^{h}(θ, 0) 0 < θ < θ_{0} with nonuniform breakpoints and (m, k) = (0, 2)

A continuation method

To simplify the notation, we will drop the arrow over vectors in ℜ^{N}. For convenience, we introduce the following notation. Let G(λ, x) denote the left hand side of equations (4.11), where x = (f_{1}, f_{2}, …, f_{N}, w_{1}, w_{2}, …, w_{N}). In this notation, we can express equations (4.11) in the form:

(4.18)  G(λ, x) = 0,

where G : ℜ^{2N+1} ↦ ℜ^{2N}. There is a natural identification between a solution (λ, x) ∈ ℜ^{2N+1} that solves (4.18) and a solution (λ, f^{h}, w^{h}) ∈ ℜ × W^{h} × W^{h} that solves (4.11). The components of x are the coefficients of the ψ_{i}^{h}’s in (4.9). The problem of computing solutions of (∗) which bifurcate from the trivial solution at (λ_{k}^{m}, 0) can be approximated by the problem of determining the curve (λ, x(λ)) which bifurcates from the trivial solution of (4.18). While the formulation in (4.18) is adequate for tracking solutions near a branch point, it does not work as well near a turning point. To avoid these difficulties, we will incorporate a pseudo–arc–length Euler– Newton continuation method by Keller (see [9], [21]) and determine (λ(ρ), x(ρ)). For the convenience of the reader, we present a brief summary of this method. We will follow the notation used in [9].

Suppose that (λ(ρ_{0}), x(ρ_{0})) is a solution of (4.18). To determine (λ(ρ), x(ρ)) where ρ − ρ_{0} is small, we first compute a tangent to the solution branch,

(λ̇_{0}, ẋ_{0}) = (λ̇(ρ_{0}), ẋ(ρ_{0})),

where the dot denotes differentiation with respect to the pseudo–arc–length parameter ρ. In particular, (λ̇_{0}, ẋ_{0}) satisfies

G_{x}ẋ_{0} + λ̇_{0}G_{λ} = 0,
||ẋ||_{2}^{2} + |λ̇|^{2} = 1,

where x = (f, w), f ∈ ℜ^{N}, w ∈ ℜ^{N},

||x||_{2}^{2} = ||f||^{2} + ||w||^{2},
||f||^{2} = ∑_{i,j} L_{ij}f_{i}f_{j},
||w||^{2} = ∑_{i,j} L_{ij}w_{i}w_{j}.

The L_{ij}’s are defined in (4.5) and 1 ≤ i, j ≤ N. For our applications, [L_{ij}] is positive definite and the discrete norm || · || on ℜ^{N} is well defined.

We then proceed to compute an Euler step in the direction of the tangent,

(4.19)  \begin{matrix}x^{0} = x_{0} + (ρ − ρ_{0})ẋ_{0}, \\ λ^{0} = λ_{0} + (ρ − ρ_{0})λ̇_{0}. \end{matrix}

Finally, we return to the solution branch by using (4.19) as an initial guess in order to determine the solution of

(4.20a)  G(λ(ρ), x(ρ)) = 0,
(4.20b)  N(λ(ρ), x(ρ), ρ) = 0,

where

N(λ(ρ), x(ρ), ρ) = ||x(ρ) − x(ρ_{0})||_{2}^{2} + |λ(ρ) − λ(ρ_{0})|^{2} − |ρ − ρ_{0}|^{2}.

The definition of N forces the solution (λ(ρ), x(ρ)) to lie on a sphere of radius |ρ − ρ_{0}| that is centered at (λ(ρ_{0}), x(ρ_{0})). This process is repeated, generating a sequence of solutions,

(λ(ρ_{n}), x(ρ_{n})) ∈ ℜ^{2N+1},

or equivalently,

(λ(ρ_{n}), f_{ρ_{n}}^{h}(θ, ϕ), w_{ρ_{n}}^{h}(θ, ϕ)) ∈ ℜ × W^{h} × W^{h} .

Equations (4.20) were solved also using a Levenberg–Marquardt algorithm (see [18, p. 776]).

Remark. A choice for ±λ̇_{0} is made at each step to ensure that the algorithm will track solution curves for increasing values of ρ (see [9]). Thus, two solutions (corresponding to different values of λ^{h}) must be known before the pseudo–arc–length method can be employed. This can be accomplished by solving (4.18) near a bifurcation point for two different values of λ (see Table 4 and the results for λ = λ_{1}^{h} and λ = λ_{2}^{h}).

The energy associated with a solution in the form (3.10)–(3.11) is given by

E(λ, w, f) = \frac{1}{2}[(w, w) − λ(Aw, w) + (f, f)]

where (·, ·) is the inner product on W. We define the energy of (λ^{h}, w^{h}, f^{h}) to be

E^{h}(λ^{h}, w^{h}, f^{h}) = \frac{1}{2}[ℓ(w^{h}; w^{h}) − λ^{h}a(w^{h}; f^{h}) + ℓ(f^{h}; f^{h})].

Next, we present solution branches which bifurcate from the trivial solution for the spherical cap SS–62. We will use a nonuniform breakpoint sequence to generate a finite element space for axisymmetric solutions and a uniform breakpoint sequence to generate a finite element space for nonaxisymmetric solutions (see Appendix B). As we tracked the bifurcation curves, the residuals r(λ^{h}(ρ_{n}), f_{ρ_{n}}^{h}, w_{ρ_{n}}^{h}) remained about the same order of magnitude as given in Column 7 of Table 4. Exceptions to this occurred near turning points. However, it was found that the residuals could be lowered near turning points by reducing the stepsize ρ − ρ_{0}.

Before we present the bifurcation diagrams, we first define some notation. Two solution branches bifurcate from each of the branch points (λ_{1}^{1}, 0), (λ_{1}^{2}, 0), and (λ_{2}^{0}, 0). Given a solution of (4.11), in the form (λ^{h}(ρ), f_{ρ}^{h}(θ, ϕ), w_{ρ}^{h}(θ, ϕ)), we define

C_{λ_{k}^{m}}^{+} = {(λ^{h}(ρ), f_{ρ}^{h}(θ, ϕ), w_{ρ}^{h}(θ, ϕ)) | lim_{ρ→0} λ^{h}(ρ) = λ_{k}^{m}, lim_{ρ→0} sign[w_{ρ}^{h}(θ, ϕ)/u_{k}^{m}(θ, ϕ)] = +1} .

The second branch, C_{λ_{k}^{m}}^{−} is defined in an analogous fashion except that

lim_{ρ→0} sign[w_{ρ}^{h}(θ, ϕ)/u_{k}^{m}(θ, ϕ)] = −1.

Note, when (m, k) = (0, 2), the above definitions imply that C_{λ_{2}^{0}}^{−} is the subcritical branch and C_{λ_{2}^{0}}^{+} is the supercritical branch. By construction, all branches start at (λ_{k}^{m}, 0, 0).

The branches C_{λ_{k}^{m}}^{σ} for σ = + or σ = − are curves in ℜ × W^{h} × W^{h}. We will look at a number of two-dimensional curves constructed from these branches. To facilitate this, we define four operators:

S_{a}[C_{λ_{k}^{m}}^{σ}] = (λ(ρ), w_{ρ}^{h}(0, 0)),
S_{b}[C_{λ_{k}^{m}}^{σ}] = (λ(ρ), w_{ρ}^{h}(\frac{1}{2}θ_{0}, 0)),
S_{c}[C_{λ_{k}^{m}}^{σ}] = (λ(ρ), E^{h}(λ(ρ), f_{ρ}^{h}, w_{ρ}^{h})),
S_{d}[C_{λ_{k}^{m}}^{σ}] = (λ(ρ), ||w_{ρ}^{h}||),

for σ = + or σ = −. We will let C_{λ_{k}^{m}} = C_{λ_{k}^{m}}^{+} ∪ C_{λ_{k}^{m}}^{−}.

In Figure 4a, we plot the initial portions of S_{a}[C_{λ_{1}^{1}}^{+}] and S_{a}[C_{λ_{1}^{1}}^{−}]. Note, the traces of these two curves coincide. In Figure 4b, we plot initial portions of S_{b}[C_{λ_{1}^{1}}^{−}] (the solid curve) and S_{b}[C_{λ_{1}^{1}}^{+}] (the dashed curve). The branches C_{λ_{1}^{1}}^{−} and C_{λ_{1}^{1}}^{+} meet at a certain distinguished point, denoted by Q^{∗} in Figure 4b. Had we continued both branches for a sufficiently large enough ρ, we would see that C_{λ_{1}^{1}}^{−} and C_{λ_{1}^{1}}^{+} are closed curves with the same trace. It turns out that Q^{∗} also lies on the branch C_{λ_{1}^{2}}^{+}. Figure 4c presents the curves S_{c}[C_{λ_{1}^{1}}^{−}] and S_{c}[C_{λ_{1}^{1}}^{+}]. Figure 4d presents the curves S_{d}[C_{λ_{1}^{1}}^{−}] and S_{d}[C_{λ_{1}^{1}}^{+}].

Figure 4a. S_{a}[C_{λ_{1}^{1}}^{±}] – (λ^{h}(ρ), w_{ρ}^{h}(0, 0)) Figure 4b. S_{b}[C_{λ_{1}^{1}}^{±}] – (λ^{h}(ρ), w_{ρ}^{h}(θ_{0}/2, 0))

For every solution on a branch C_{λ_{k}^{m}}^{+} or C_{λ_{k}^{m}}^{−}, we can consider a profile of the corresponding solution w_{ρ}^{h}(θ, ϕ) for 0 < θ < θ_{0}, ϕ = 0 and ϕ = π (or equivalently −θ_{0} < θ < θ_{0} and ϕ = 0). In Figure 5a, we present the profiles of a number of different solutions on the branch C_{λ_{1}^{1}}^{−}. Figure 5b contains similar profiles for the branch C_{λ_{1}^{1}}^{+}.

Figure 5a. C_{λ_{1}^{1}}^{−} solution profiles

Figures 6-7 relate to the branch point (λ_{1}^{2}, 0).

Figure 6a. S_{a}[C_{λ_{1}^{2}}^{±}] – (λ^{h}(ρ), w_{ρ}^{h}(0, 0)) Figure 6b. S_{b}[C_{λ_{1}^{2}}^{±}] – (λ^{h}(ρ), w_{ρ}^{h}(θ_{0}/2, 0))

Figure 7a. C_{λ_{1}^{2}}^{−} solution profiles

Figures 8–9 are related to the branch point (λ_{2}^{0}, 0), with the difference that the solid curves represent curves branching subcritically and dashed curves represent curves which branch supercritically.

Figure 8a. S_{a}[C_{λ_{2}^{0}}^{±}] – (λ^{h}(ρ), w_{ρ}^{h}(0, 0)) Figure 8b. S_{b}[C_{λ_{2}^{0}}^{±}] – (λ^{h}(ρ), w_{ρ}^{h}(θ_{0}/2, 0))

Figure 8a′ presents a magnification of Figure 8a near the bifurcation point. Figure 8d′ presents a magnification of Figure 8d near the bifurcation point.

Figure 8a′. Magnification of (λ^{h}(ρ), w_{ρ}^{h}(0, 0)) near the branch point λ_{2}^{0} − η^{∗} < λ^{h}(ρ) < λ_{2}^{0} + η^{∗}

Figure 8c′. Magnification of (λ^{h}(ρ), E^{h}(ρ)) near the branch point for λ^{h}(ρ) < λ_{2}^{0}

Figure 9a. C_{λ_{2}^{0}}^{−} solution profiles

Figures 10 combine the diagrams that were presented in Figures 4, 6, and 8. Figures 10 are scaled differently than Figures 4, 6, and 8. Figure 10c′ is a magnification of Figure 10c. The curves presented in Figures 10 are dashed or solid to distinguish branches C_{λ_{k}^{m}}^{±} for (m, k) = (1, 1), (2, 1), (0, 2).

Figure 10a. S_{a}[C_{λ_{k}^{m}}^{±}] for (m, k) = (1, 1), (2, 1), (0, 2)

Figure 10b. S_{b}[C_{λ_{k}^{m}}^{±}] for (m, k) = (1, 1), (2, 1), (0, 2)

Figure 10c. S_{c}[C_{λ_{k}^{m}}^{±}] for (m, k) = (1, 1), (2, 1), (0, 2)

Figure 10d. S_{d}[C_{λ_{k}^{m}}^{±}] for (m, k) = (1, 1), (2, 1), (0, 2)

Figure 10c′. Local view of S_{c}[C_{λ_{k}^{m}}^{±}] for (m, k) = (1, 1), (2, 1), (0, 2)

For small ρ, solutions on C_{λ_{1}^{1}}^{±} appear to preserve the symmetry of the eigenfunction u_{1}^{1}(θ, ϕ), i.e., w_{ρ}^{h}(θ, ϕ) = −w_{ρ}^{h}(θ, −ϕ). However, this symmetry is not preserved globally by the numerical solutions for large ρ. In fact, w_{ρ^{∗}}^{h}(θ, ϕ) = w_{ρ^{∗}}^{h}(θ, −ϕ) for some ρ = ρ^{∗}. This property characterizes the solution at the previously introduced point Q^{∗}. Q^{∗} is a secondary bifurcation point and is the unique point that lies on both branches C_{λ_{1}^{1}} and C_{λ_{1}^{2}} (this was verified numerically by comparing all the components of a solution on C_{λ_{1}^{1}} with a solution on C_{λ_{1}^{2}}^{+}).

5. Constructing a “buckled” surface. In the following section, it will be convenient to define (u^{1}, u^{2}) = (θ, ϕ). Vectors in ℜ^{3} will be denoted by bold face. As in Section 2, we will denote the parametrization of the middle surface of the undeformed spherical cap of radius R by X : Ω_{θ_{0}} ⊂ ℜ^{2} → S ⊂ ℜ^{3}. A parametrization for the middle surface of the deformed shell Ŝ will be denoted by X̂ : Ω_{θ_{0}} ⊂ ℜ^{2} → Ŝ ⊂ ℜ^{3}. Partial derivatives of X̂ with respect to u^{i} will be abbreviated X̂_{i} = ∂ X̂ /∂u^{i}. In the following, all surfaces under consideration are assumed to have a parametrization X̂ ∈ C^{3}(Ω_{θ_{0}}). For the purposes of this section, we will assume that we have computed a solution of (4.11) in the form

(λ^{h}, f^{h}, w^{h}) = (λ^{h}(ρ_{n}), f_{ρ_{n}}^{h}(θ, ϕ), w_{ρ_{n}}^{h}(θ, ϕ))

where ρ_{n} > 0 is fixed.

A standard result in differential geometry concerns the existence of a surface Ŝ with prescribed first and second fundamental forms, â_{ij} and b̂_{ij}. In general, â_{ij} and b̂_{ij} by themselves are not sufficient to guarantee the existence of such a surface. However, if three compatibility equations are met, then a surface Ŝ does exist with the fundamental forms â_{ij} and b̂_{ij}. These compatibility equations are commonly referred to as the Mainardi–Codazzi equations (see [31, Eqns. (6.22)]) and Gauss’s equation (see [31, Eqn. (6.26)]). In this section, we outline the method by which we determine the parametrization X̂(u^{1}, u^{2}) from the forms â_{ij} and b̂_{ij}. We first state a result concerning the existence of Ŝ:

Theorem ([31, p. 146]). Suppose that â_{ij} are twice continuously differentiable, b̂_{ij} have continuous first derivatives, and â_{11} â_{22} − â_{12} â_{21} > 0 in a neighborhood of a point (u_{0}, v_{0}). Furthermore, suppose that the vectors

(5.1)  X̂_{i}(u_{0}, v_{0}) = \frac{∂ X̂}{∂u^{i}} |_{(u_{0},v_{0})}, i = 1, 2,

are given along with the corresponding normal

(5.2)  N̂(u_{0}, v_{0}) = \frac{X̂_{1} × X̂_{2}}{|X̂_{1} × X̂_{2}|}|_{(u_{0},v_{0})}.

Then there exists a parametrization X̂(u^{1}, u^{2}) which is valid for all (u^{1}, u^{2}) in neighborhood of (u_{0}, v_{0}) with coordinate vectors {X̂_{1}(u^{1}, u^{2}), X̂_{2}(u^{1}, u^{2}), N̂(u^{1}, u^{2})} that satisfy the initial conditions (5.1)–(5.2).

The parametrization X̂ is unique up to a translation and orthogonal rotation. In practice, an initial position vector X̂_{0} is specified so that in our applications Ŝ is unique (see the remark following Eqn. (5.4)).

Remark. The first equation of (2.1) is an approximation to Gauss’s equation, while the second and third equations approximate the Mainardi–Codazzi equations. Although the fundamental forms of the deformed shell do not satisfy the exact compatibility equations, they do satisfy a set of equations which approximate the compatibility equations. In fact, this is the best that one could expect, since the constitutive laws are already approximations to begin with. For a further discussion of the compatibility equations for the middle surface of an elastic shell, the reader is referred to [24].

Using the methods outlined in Sections 3–4, we can approximate the solutions of the John equations f and w (and thus approximate F and W using (3.2)). Using equations (3.1)–(3.2) and the constitutive equations (2.2), we are able to compute â_{ij}(u^{1}, u^{2}) and b̂_{ij}(u^{1}, u^{2}) for all (u^{1}, u^{2}) ∈ Ω_{θ_{0}} in terms of the potentials W and F and the load λ. By construction, the elements of the basis B^{h} are four times continuously differentiable so that â_{ij} and b̂_{ij} have the desired amount of smoothness (see Appendix B). Thus, it is possible to construct the buckled surface Ŝ.

The differential equations of a surface Ŝ are 15 scalar differential equations involving X̂_{j} and N̂,

(5.3)  \begin{aligned}\frac{∂ X̂_{j}}{∂u^{i}} & = Γ̂_{ij}^{r} X̂_{r} + b̂_{ij} N̂, i, j = 1, 2, \\ \frac{∂ N̂}{∂u^{i}} & = b̂_{i}^{r} X̂_{r}, i = 1, 2, \end{aligned}

where repeated indices on the right hand side of (5.3) are summed. The eight Christoffel symbols Γ̂_{jk}^{i} can be computed directly from the metric tensor â_{ij} and its derivatives of order one (see [31, p. 135]). If the compatibility conditions are met, it is then valid to integrate along the coordinate curves to determine the parametrization for Ŝ. For example, suppose we hold the variable u^{1} = u_{0} fixed and vary u^{2}. To compute the position vector X̂ and the triad {X̂_{1}, X̂_{2}, N̂}, we consider the following twelve nonlinear differential equations:

(5.4)

\frac{∂ X̂_{1}}{∂u^{2}} = Γ̂_{12}^{r} X̂_{r} + b̂_{12} N̂,
\frac{∂ X̂_{2}}{∂u^{2}} = Γ̂_{22}^{r} X̂_{r} + b̂_{22} N̂,
\frac{∂ N̂}{∂u^{2}} = b̂_{2}^{r} X̂_{r},
\frac{∂ X̂}{∂u^{2}} = X̂_{2},

with the “initial conditions” (5.1)–(5.2) and X̂(u_{0}, v_{0}) = X̂_{0}. In general, a solution obtained by numerically integrating (5.4) along its coordinate curves will be a discrete set in the form

(5.5)  {(X̂(u_{0}, u_{k}^{2}), X̂_{1}(u_{0}, u_{k}^{2}), X̂_{2}(u_{0}, u_{k}^{2}), N̂(u_{0}, u_{k}^{2})) , k = 1, 2, …, n}

where 0 < u_{1}^{2} < u_{2}^{2} < ⋯ < u_{n}^{2} < 2π for some integer n. By reversing the roles of u^{1} and u^{2} one can integrate along the coordinate curves u^{2} = v_{0} in a similar fashion. A fifth–order integration scheme was used to integrate (5.3) along the coordinate curves.

Since we can carry out this procedure for any point on a solution branch C_{λ_{k}^{m}}^{±}

parametrized by ρ, we can think of the corresponding surfaces Ŝ_{ρ} as forming a one– parameter family of bifurcating elastic surfaces.

Remark 1. It should be emphasized here that the result concerning the existence of Ŝ is a local one (i.e., it is valid only in a neighborhood of (u_{0}, v_{0})). Classical results regarding the global existence of a surface typically require that the surface be closed, regular and convex (see, e.g., [31, p. 358]). Recent work regarding boundary value problems for surfaces require a condition on the Gauss curvature K (see, e.g., [13] or [16]). These results do not apply to the middle surface of the deformed cap. In practice, the integration of (5.3) posed no difficulties for a cap of size θ_{0} = 12.85^{∘}.

Remark 2. To fix an orientation and an origin for the surface Ŝ, we assumed that near (u_{0}, v_{0}) ≈ (0, 0), the parametrization had the form,

(5.6)  X̂(u^{1}, u^{2}) ≈ X(u^{1}, u^{2}) + W (u^{1}, u^{2}) N(u^{1}, u^{2}),

where N is the unit normal on Ŝ. The vectors, X̂(u_{0}, v_{0}), X̂_{1}(u_{0}, v_{0}), X̂_{2}(u_{0}, v_{0}), and N̂(u_{0}, v_{0}) were computed from (5.6). We do not assume that Ŝ has a global parametrization in the form (5.6). Eqn. (5.6) was chosen only to fix the initial position vector X̂(u_{0}, v_{0}) and orientation of the initial triad. Any other choice could have been used to begin the integration of (5.3).

6. Buckled states of a spherical cap. For each point on a solution branch, there corresponds a buckled state of the shell. We selected a number of points on each of the branches C_{λ_{1}^{1}}^{±}, C_{λ_{1}^{2}}^{±} and C_{λ_{2}^{0}}^{±} and constructed the corresponding buckled shells. In Figure 10b, the selected points lying on a branch in the form C_{λ_{k}^{m}}^{−} are indicated with an “o”, while points corresponding to C_{λ_{k}^{m}}^{+} are indicated with an “x.”

The actual deformation corresponding to C_{λ_{1}^{1}}^{±} was very small (see Figure 11a). To bring out the asymmetries of these states, we stretched out the z–coordinate (see Figures 11b–11e). Figures 11a, 12a–12d, 13a–13d, 14a–14d, 15a–15d are drawn to scale.

Figure 11a. Buckled state on C_{λ_{1}^{1}}^{+}

Figure 11b. Scaled buckled state on C_{λ_{1}^{1}}^{−} Figure 11c. Scaled buckled state near Q^{∗} on C_{λ_{1}^{1}}^{+}

Figure 11d. Scaled buckled state on C_{λ_{1}^{1}}^{−} Figure 11e. Scaled buckled state near Q^{∗} on C_{λ_{1}^{1}}^{+} Figure 12a. Buckled state on C_{λ_{1}^{2}}^{−} Figure 13a. Buckled state on C_{λ_{1}^{2}}^{+}

Figure 12b. Buckled state on C_{λ_{1}^{2}}^{−} Figure 13b. Buckled state on C_{λ_{1}^{2}}^{+}

Figure 12c. Buckled state on C_{λ_{1}^{2}}^{−} Figure 13c. Buckled state on C_{λ_{1}^{2}}^{+} Figure 14a. Buckled state on C_{λ_{2}^{0}}^{−} Figure 15a. Buckled state on C_{λ_{2}^{0}}^{+}

Figure 14b. Buckled state on C_{λ_{2}^{0}}^{−} Figure 15b. Buckled state on C_{λ_{2}^{0}}^{+}

Figure 14c. Buckled state on C_{λ_{2}^{0}}^{−} Figure 15c. Buckled state on C_{λ_{2}^{0}}^{+}

7. Concluding remarks. The numerical work of [6] and [7] showed that the global structure of numerical solution branches can be complicated. In the work of [6] on rectangular plates and the work of [7] on hemispheres and spheres, the authors discovered smooth transitions between solution branches of different modes. Our numerical computations indicate that a similar transition takes place at Q^{∗} (see Figure 10b). The point Q^{∗} is a secondary bifurcation point that connects the branches C_{λ_{1}^{1}}^{±} and C_{λ_{1}^{2}}^{+}. This behavior persists if the number of finite elements is changed. The numerical solutions in [6] were based on finite differences and used the von Kármán equations to model the buckling of a plate. While our work dealt with spherical caps, used a different model, and a different numerical scheme, similar global behavior of numerical solution branches was observed. It appears that a numerical solution can lose certain symmetries as one moves along a branch (see, e.g., C_{λ_{1}^{1}}). The transition from C_{λ_{1}^{1}}^{±} to C_{λ_{1}^{2}}^{+} takes place in a smooth fashion. At Q^{∗}, a critical point appears at the north pole of the solution potential w_{ρ^{∗}}^{h} on C_{λ_{1}^{1}}, where originally there was none. To prevent this type of behavior, one could consider (4.11) on certain subspaces of W^{h} where a certain symmetry can be enforced. For solutions branching from (λ_{1}^{1}, 0), one might consider the subspace,

M^{h} = \{w^{h} | w^{h} ∈ W^{h}, w^{h}(θ, \frac{1}{2}π − ϕ) = −w^{h}(θ, \frac{1}{2}π + ϕ)\} .

Solutions in M^{h} have nodal lines fixed at ϕ = ± \frac{1}{2} π. This type of approach has been used in [11], where numerical solutions with prescribed symmetries and nodal sets were found in certain invariant subspaces. The analytical work of [14] and [15] utilizes the preservation of nodal structure to obtain global characterizations of solution branches of certain nonlinear partial differential equations on planar domains. In order for such an approach to work, the domain must possess certain symmetries to begin with. If one seeks solutions in an invariant subspace, one may not be able to detect certain secondary bifurcations which may be present in the original system. Had we restricted attention to the subspace M^{h}, then we would not have observed the transition from a C_{λ_{1}^{1}}–solution to a C_{λ_{1}^{2}}–solution. From (3.11), we see that

(7.1)  |η|^{− \frac{1}{2}} w_{η}^{±}(θ, ϕ) = β^{±}(η)u_{k}^{m} + |η|^{\frac{1}{2}} v(β^{±}(η)u_{k}^{m}) + o(|η|^{\frac{1}{2}} ), m ≠ 0.

The form of v(β^{±}(η)u_{k}^{m}) and [5, Eqn. (A.7)] suggests that w_{η}^{±}(0, 0) ≠ 0 when η ≠ 0. Therefore, the nodal structure of w_{η}^{±} must change near θ = 0.

Although the John equations were derived under the assumption that the strains are small, many solutions that are presented in Section 6 could not be classified as having small strains. Nevertheless, even when the strains were not small, it was possible to integrate (5.3) and obtain a reasonable solution.

Our numerical computations were carried out on an IBM 3090VF model J. Using 187 elements, it took approximately 0.23 minutes to take one step along a nonaxisymmetric solution branch following the continuation method described in Section 4. The image presented in Figure 13d, required computing X̂(θ, ϕ) at 2100 grid points, \{(θ_{i}, ϕ_{j}) | θ_{i} = \frac{θ_{0}i}{30}, i = 1, …, 30, ϕ_{j} = \frac{2πj}{70} , j = 1, …, 70\} and took 3.27 minutes of cpu. Note, there is no correlation between the number of grid lines presented in the figures in Section 6 and the number of finite elements used in Section 4. The total time could be reduced by decreasing the number of grid points, since an integration along a typical coordinate curve ϕ = ϕ_{j} only requires about 0.0467 minutes of cpu. For an axisymmetric solution, the computational time was reduced since only 38 elements are used. On the average, it took approximately 0.008 minutes of cpu to take one step along an axisymmetric solution branch. It took 1.9 minutes to complete the computations needed to construct the surface presented in Figure 15d. This time could be reduced by a factor of 70, since only one integration along a curve ϕ = ϕ_{j} is needed for an axisymmetric solution.

The three-dimensional images presented in Section 6 were produced at the Center for Geometry Analysis Numerics and Graphics, University of Massachusetts, Amherst using software developed by James T. Hoffman.

Acknowledgement. The author would like to thank Professor David Hoffman for making available the facilities at the Center for Geometry Analysis Numerics and Graphics (G.A.N.G.) in support of this work. The author is grateful to David Oliver and Ed Thayer for their technical assistance at G.A.N.G. The author also would like to thank Augustin Dubrulle, IBM Scientific Research Center, Palo Alto for a number of helpful suggestions for optimizing the code which was developed for this research.

Appendix A. The eigenvalues of the Laplacian on a spherical cap. On a spherical domain, the Laplacian (Δ) is given by

(A.1)  Δu = u_{θθ} + \frac{1}{sin^{2} θ} [u_{ϕϕ} + sin θ cos θu_{θ}].

If (μ, u) satisfies

(Δ)  \begin{matrix}Δu + μu = 0 on Ω_{θ_{0}}, \\ u = 0 on ∂ Ω_{θ_{0}}, \end{matrix}

with u ≠ 0, then μ is called an eigenvalue of the Laplacian and u is called an eigenfunction. The eigenfunctions of (Δ) are:

(A.2a)  u_{k}^{m}(θ, ϕ) = c_{k}^{m}P_{ν_{k}^{m}}^{m}(cos θ) cos mϕ, m = 0, 1, …, k = 1, 2, …,
(A.2b)  v_{k}^{m}(θ, ϕ) = d_{k}^{m}P_{ν_{k}^{m}}^{m}(cos θ) sin mϕ, m = 1, 2, …, k = 1, 2, …,

where P_{ν}^{m} is a Legendre function of the first kind of degree ν and order m. The eigenvalue corresponding to u = u_{k}^{m} is given by μ = μ_{k}^{m} where

(A.3)  μ_{k}^{m} = ν_{k}^{m}(ν_{k}^{m} + 1),

and

(A.4)  P_{ν_{k}^{m}}^{m}(cos θ_{0}) = 0.

The constants c_{k}^{m} and d_{k}^{m} are normalizing coefficients chosen so that a(u_{k}^{m}; u_{k}^{m}) = a(v_{k}^{m}; v_{k}^{m}) = 1 (see Eqn. (3.8)).

In general, the μ_{k}^{m}’s need not be distinct. For example, when θ_{0} ˙=35.8^{∘}, we find that μ_{1}^{6} = μ_{3}^{1}. However, the set

Θ_{0} = {θ_{0} | 0 < θ_{0} < π, μ_{k}^{m}(θ_{0}) = μ_{j}^{n}(θ_{0}), (m, k) ≠ (n, j)},

has measure zero (see [3]–[4]). Thus, for almost every θ_{0}, μ_{k}^{m}(θ_{0}) has multiplicity one when m = 0 and multiplicity two when m ≠ 0. When we restrict our attention to functions which are even about ϕ = 0, we see that the eigenfunctions are given by (A.2a). The set {u_{k}^{m}} is complete in W (See [5]).

Appendix B. Constructing a basis for the solution space. The numerical solutions w^{h}, f^{h} that are presented in Sections 4–6 were constructed to lie in W^{h}, the linear space of functions spanned by the elements of B^{h}. In this section, we will describe how the set of finite elements

B^{h} = {ψ_{j}^{h}(θ, ϕ) | j = 1, 2, …, N}.

is constructed. We will follow the terminology that is defined in [12] and [18] and the reader is referred to these references for a more complete exposition on splines.

Our choice of B^{h} is motivated by the theorem in Section 5 that we wish to apply. Since the first and second fundamental forms of the deformed shell need to be sufficiently smooth, our numerical solutions also must possess a certain amount of smoothness. In particular, we must have â_{ij} ∈ C^{2}(Ω_{θ_{0}}) (and therefore, ρ_{ij} ∈ C^{2}(Ω_{θ_{0}})) and b̂_{ij} ∈ C^{1}(Ω_{θ_{0}}). Using the constitutive relations, we see that f^{h} must be at least four times continuously differentiable on Ω_{θ_{0}} and w^{h} must be at least three times continuously differentiable on Ω_{θ_{0}}. Solutions in the form (3.10)–(3.11) were shown to be analytic in [5]. We will require f, w ∈ C^{4}(Ω_{θ_{0}}) and build this property into the space W^{h}.

In the classical formulation, f and w must satisfy the boundary conditions (3.6c)– (3.6d). We will be working with the weak equations (∗) where the boundary conditions require that a test function ϕ vanishes on ∂ Ω_{θ_{0}}. Classical solutions generated by weak solutions will satisfy the appropriate boundary conditions. In the following, we will show how the continuity requirements and boundary conditions are built into the space W^{h}.

For ψ_{j}^{h}(θ, ϕ) ∈ B^{h}, we require that the following conditions must hold: (C_{1}) ∂^{n}ψ_{j}^{h}/∂θ^{n} are continuous for θ ∈ (0, θ_{0}), n ≤ 4. (C_{2}) ∂^{n}ψ_{j}^{h}/∂ϕ^{n} are continuous for ϕ ∈ (0, 2π) and periodic with period 2π for n ≤ 4. (C_{3}) ψ_{j}^{h}(θ_{0}, ϕ) = 0, 0 < ϕ < 2π. (C_{4}) If w^{h} = w^{h}(θ) ∈ W^{h}, then

\frac{dw^{h}}{dθ} (0^{+}) = 0.

Remark. Note that (C_{4}) is a necessary condition for an axisymmetric function in C^{1}(Ω_{θ_{0}}).

We begin by defining some terminology and notation:

(1) P_{k,ξ} is the linear space of piecewise polynomial functions of order k with breakpoint sequence ξ = (ξ_{1}, …, ξ_{ℓ+1}) (see [12, p. 87]).

(2) P_{k,ξ,ν} is the linear subspace of P_{k,ξ} which satisfy continuity requirements specified by ν (see [12, p. 100]).

(3) B_{i} = B_{i,k,t} is the i^{th} B–spline of order k with knot sequence t. The knot sequence t consists of elements of the breakpoint sequence ξ and is used to specify continuity requirements at the internal breakpoints ξ_{j}, j = 2, 3, … ℓ. In particular, if the j^{th} derivative is required to be continuous at such a ξ_{j}, then ξ_{j} must occur k − j + 1 times in the corresponding knot sequence t. If M denotes the length of t, then {B_{i} | i = 1, …, M − k} is a basis of dimension m = M − k for the space P_{k,ξ,ν}. Note, as a basis for the space of piecewise polynomials, the B_{i}’s are linearly independent.

For our applications, we set the order k = k_{θ} = 6 and require four continuous derivatives at interior breakpoints so that ξ_{2}, …, ξ_{ℓ} each appears at most once in t. We choose the length of our knot sequence to be 46, ξ_{1} = 0.0 and ξ_{ℓ+1} = θ_{0}. In particular, t_{i} = ξ_{1} = 0 for i = 1, …, 6 and t_{i} = ξ_{ℓ+1} = θ_{0} for i = 40, …, 46. In our work, we used the ”not-a-knot” rule to generate t from the breakpoint sequence ξ (see [12]).

We considered two types of breakpoint sequences. In the first type, the breakpoints are distributed uniformally on (0, θ_{0}). In the second, breakpoints near the pole were spaced further apart than those near the edge θ = θ_{0} (this does not alter the smoothness of the corresponding B_{i}’s). It was found that by choosing the breakpoints in the second way, it was possible to get a better approximation to an axisymmetric function. For our test case, we find θ_{0} = 0.224 radians, M = 46, and the dimension of the corresponding space P_{k,ξ,ν} is m = 46 − 6 = 40. With equally spaced knots, we find

\begin{matrix}t = ( & 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.017, 0.023, 0.028, 0.034, 0.040, 0.046, 0.051, 0.057, 0.063, 0.069, \\ 0.074, 0.080, 0.086, 0.092, 0.097, 0.104, 0.109, 0.115, 0.121, 0.127, 0.132, 0.138, 0.144, 0.150, 0.155, \\ 0.161, 0.167, 0.173, 0.178, 0.184, 0.190, 0.196, 0.201, 0.207, 0.224, 0.224, 0.224, 0.224, 0.224, 0.224). \end{matrix}

For a non-uniformally distributed breakpoint sequence, we find

\begin{matrix}t = ( & 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.105, 0.109, 0.112, 0.115, 0.119, 0.122, 0.125, 0.128, 0.132, 0.135, \\ 0.138, 0.142, 0.145, 0.148, 0.152, 0.155, 0.158, 0.161, 0.165, 0.168, 0.171, 0.175, 0.178, 0.181, 0.185, \\ 0.188, 0.191, 0.195, 0.198, 0.201, 0.204, 0.208, 0.211, 0.214, 0.224, 0.224, 0.224, 0.224, 0.224, 0.224). \end{matrix}

After imposing boundary conditions, we disregard B_{m}(θ) because B_{m}(θ_{0}) ≠ 0. We define m_{θ} = m − 1, so that in our applications m_{θ} = 39. By construction, we see that B_{i}(θ_{0}) = 0 for i = 1, …, m_{θ}. We are led to the following set of functions,

{B_{1}(θ) + B_{2}(θ), B_{2}(θ), …, B_{m_{θ}}(θ)} .

Note, by definition (B_{1} + B_{2})(0) = 1, (B_{1} + B_{2})′(0) = 0, B_{2}(0) = 0, and B′_{2}(0) ≠ 0. If we are interested in axisymmetric solutions, then B_{2}(θ) is excluded from the above set. Figure B.1 presents a graph of a typical interior spline where k_{θ} = 6. In Figures B.2 and B.3, we plot a number of splines for k_{θ} = 6. In Figures B.1, B.2, and B.3, the breakpoints are uniformally distributed.

We are now able to define our basis functions. Let n_{i} denote nonnegative integers, n_{0} = 0, n_{1} = 1, n_{i} < n_{i+1}. Define

(B.1)

ψ_{1,0}^{h}(θ, ϕ) = B_{1}(θ) + B_{2}(θ),
ψ_{i,0}^{h}(θ, ϕ) = B_{i}(θ),  i = 3, 4, …, m_{θ},
ψ_{i,1}^{h}(θ, ϕ) = B_{i}(θ) cos ϕ, i = 2, 3, …, m_{θ},
ψ_{i,n_{2}}^{h}(θ, ϕ) = B_{i}(θ) cos n_{2}ϕ, i = 3, 4, …, m_{θ},
⋮
ψ_{i,n_{ϕ}}^{h}(θ, ϕ) = B_{i}(θ) cos n_{ϕ}ϕ, i = 3, 4, …, m_{θ}.

Remark. The labeling of the ψ_{i,j}^{h}’s is not important here and to simplify the notation we will refer to them with a single subscript in Sections 1–7. At this point we see that the quantities L_{nm}, A_{nm}, B_{nmq} used in Sections 3–4 are well-defined and can be computed numerically. one-parameter families of bifurcating elastic surfaces 37

Figure B.1 A typical interior spline of order 6

Figure B.2. B_{1}(θ) + B_{2}(θ), B_{3}(θ), … B_{7}(θ) Figure B.3. B_{33}(θ), B_{34}(θ), … B_{39}(θ)

Since the quantity Δψ_{i,j}^{h} occurs in the calculation of the matrices [L_{nm}], [A_{nm}], and [B_{nmq}], it is worthwhile to note that Δψ_{i,j}^{h} remains bounded as θ ↓ 0. Using the properties of the functions B_{i}, it is easy to see that the ψ_{i,j}^{h}’s are well behaved at θ = 0^{+}, except possibly ψ_{2,1}^{h}. Using the definition of ψ_{2,1}^{h}, it follows that

(B.1)  Δψ_{2,1}^{h} = B′′_{2}(θ) cos ϕ + \frac{1}{sin^{2} θ} [−B_{2}(θ) + sin θ cos θB′_{2}(θ)] cos ϕ.

By the properties of B_{2}, we see that for small θ,

−B_{2}(θ) + sin θ cos θB′_{2}(θ) = O(θ^{2})

and Δψ_{2,1}^{h} remains bounded as θ ↓ 0.

If u_{k}^{m} is an eigenfunction of the Laplacian, then

\begin{matrix}\frac{∂u_{k}^{m}}{∂θ} |_{(0^{+},ϕ)} & = 0, m ≠ 1, \\ \frac{∂u_{k}^{1}}{∂θ} |_{(0^{+},ϕ)} & ≠ 0. \end{matrix}

Using elements from B^{h}, one can construct functions which behave the same way near the north pole. In previous work of this type, it is not uncommon to require that the second derivatives of the finite elements vanish at θ_{0} = 0 (see [17] and [34]). Our work differs in this respect. By construction ∂^{2}ψ_{1,0}^{h}/∂θ^{2} and ∂^{2}ψ_{2,1}^{h}/∂θ^{2} do not vanish at θ = 0.

In the nonaxisymmetric cases, we chose

k_{θ} = 6, m_{θ} = 39, n_{ϕ} = 6, n_{0} = 0, n_{1} = 1, n_{2} = 2, n_{3} = 4, n_{4} = 6.

One can show that for the examples considered here, solutions do not possess the symmetry of the dihedral group D_{3} so that the function cos 3ϕ is not needed for B^{h}. Moreover, frequencies higher than n_{ϕ} = 6 were not needed, since the corresponding amplitudes were negligible. For nonaxisymmetric solutions, the number of elements in B^{h} is N = 187. In the axisymmetric cases, we used

k_{θ} = 6, m_{θ} = 39, n_{ϕ} = 1 n_{0} = 0,

yielding a total of N = 38 functions in B^{h}. one-parameter families of bifurcating elastic surfaces 39

References

[1] S. S. Antman & Ke-Gang Shih, Qualitative properties of large buckled states of spherical shells, Arch. Rational Mech. Anal., 93 (1986), pp. 357–384.

[2] F. Baginski, Axisymmetric and nonaxisymmetric buckled states of a shallow spherical cap, Q. Appl. Math., 46 No. 2 (1988), pp. 331–351.

[3] F. Baginski, Ordering the zeroes of the Legendre functions P_{ν}^{m}(z_{0}) when considered as a function of ν, JMAA, 147 (1990), pp. 296–308.

[4] F. Baginski, Comparison Theorems for the ν–zeroes of Legendre functions P_{ν}^{m}(z_{0}) when −1 < z_{0} < 1, Proc. of AMS, 111(2) (1991), pp. 395-402.

[5] F. Baginski, The buckling of elastic spherical caps, J. of Elasticity, 25 (1991), pp. 159–192.

[6] L. Bauer, H. B. Keller, & E. L. Reiss, Axisymmetric buckling of hollow spheres and hemispheres, Comm. Pure Appl. Math, 23 (1970), pp. 529–568.

[7] L. Bauer & E. L. Reiss, Nonlinear buckling of rectangular plates, SIAM J., 13 (1965), pp. 603– 626.

[8] L. Berke & R. L. Carlson, Experimental studies of the postbuckling behavior of complete spherical shells, Experimental Mech., 8 (1968), pp. 548-553.

[9] J. H. Bolstad & H. B. Keller, A multigrid continuation method for elliptic problems with folds, SIAM J. Sci. Stat. Comput., 7 (4) (1986), pp. 1081–1104.

[10] R. L. Carlson, R. L. Sendelbeck & N. J. Hoff, Experimental studies of the buckling of complete spherical shells, Experimental Mech., 7 (1967), pp. 281-288.

[11] Y.S. Choi & P.J. McKenna, A mountain pass method for the numerical solution of semilinear elliptic problems, Nonlinear Analysis Theory, Methods & Applications (to appear).

[12] Carl de Boor, A practical guide to splines, Springer–Verlag, New York, 1978.

[13] B. Guan & J. Spruck, Boundary value problems on S^{n} for surfaces of constant Gauss curvature, Annals of Math. (to appear).

[14] T. Healey & H. Kielhoffer, Symmetry and nodal properties in global bifurcation analysis of quasi–linear elliptic equations, Technical Report (1990), Cornell University.

[15] T. Healey & P. Negron-Marrero, On the global preservation of nodal structure of the von Kármán equations for rectangular plates, preprint.

[16] D. Hoffman, H. Rosenberg, & J. Spruck, Boundary value problems for surfaces of constant Gauss curvature, Comm. Pure Appl. Math., 45 (1992), pp. 1051-1062.

[17] Nai-Chien Huang, Unsymmetrical buckling of thin shallow spherical shells, J. of Appl. Mech., 31 (1964), pp. 447–457.

[18], IMSL Stat/Library User’s Manual, Version 1.1, International Mathematical and Statistical

Libraries, Houston, 1989.

[19] F. John, Estimates for the derivatives of the stresses in a thin shell and interior equations, Comm. Pure Appl. Math., 18 (1965), pp. 235–267.

[20] F. John, Refined interior shell equations for thin elastic shells, Comm. Pure Appl. Math., 24 (1971), pp. 583–615.

[21] H.B. Keller, Numerical solutions of bifurcation and nonlinear eigenvalue problems, in Applications of Bifurcation Theory, P. Rabinowitz, ed., Academic Press, New York, 1977.

[22] G. H. Knightly & D. Sather, A priori bounds in nonlinear shell theory, in Nonlinear Phenomena in Natural Sciences, Academic Press, New York, 1982.

[23] G. H. Knightly & D. Sather, Buckled states of a spherical shell under uniform external pressure, Arch. Rational Mech. Anal., 72 (1980), pp. 315–380.

[24] W. T. Koiter, On the nonlinear theory of thin elastic shells, Proc. Kon. Nederl. Akad. Wet. Amsterdam, B69 (1966), pp. 1–54.

[25] M. A. Krenzke & T. J. Kiernan, Elastic stability of near perfect shallow spherical shells, AIAA, 1 No. 12 (1963), pp. 2855-2857.

[26] K. Marguerre, Zur Theorie der gekrümmten Plate grosser Formänderung, Proc. 5th Int. Cong. Appl. Mechs., Cambridge, Mass., 1938, pp. 93–101.

[27] F. I. Niordson, Shell Theory, Elsevier Science Publishers B. V., Amsterdam, 1985.

[28] P. H. Rabinowitz, Some global results for nonlinear eigenvalue problems, J. Funct. Anal., 7 (1971), pp. 487–513.

[29] P. H. Rabinowitz, A note on topological degree for potential operators, JMAA, 51 (1975), pp. 483–492.

[30] J. Smoller, Shock Waves and Reaction–Diffusion Equations, Springer–Verlag, New York, 1982.

[31] J. J. Stoker, Differential Geometry, Wiley, New York, 1969.

40 frank e. baginski

[32] G. Strang & G. Fix, An Analysis of the Finite Element Method, Prentice–Hall, 1973.

[33] C. Truesdell, invariant and complete stress functions for general continua, Arch. Rational Mech. Anal., 4 (1959/60), pp. 1-29.

[34] H. J. Weinitschke, On asymmetric buckling of shallow spherical shells, J. Math. and Phys., 44 (1965), pp. 141-163.