Computational stability investigations for a highly symmetric system: the pressurized spherical membrane

Thin membranes are notoriously sensitive to instabilities under mechanical loading, and need sophisticated analysis methods. Although analytical results are available for several special cases and assumptions, numerical approaches are normally needed for general descriptions of non-linear response and stability. The paper uses the case of a thin spherical hyper-elastic membrane subjected to internal gas over-pressure to investigate how stability conclusions are affected by chosen material models and kinematic discretizations. For spherical symmetry, group representation theory leads to linearized modes on the uniformly stretched sphere, with eigenvalues obtained from the mechanics of a thin membrane. A complete three-dimensional geometric description allows non-axisymmetric shear modes of the sphere, and such instabilities are shown to exist. When the symmetry of the continuous sphere is broken by discretized models, group representation theory gives predictions on the effects on the critical states. Numerical simulations of the pressurized sphere show and verify stability conclusions for sets of meshing strategies and hyper-elastic models.


Introduction
A large variety of thin pressurized membranes are common in biology and engineering [31,32,34,39]. The mechanical definition of a membrane implies that the structure will carry all loads acting on it by establishing normal force resultants. With the exception of cases where a plane membrane is subjected to tensile forces, this means that the membrane can only carry forces of distributed transversal intensities, i.e., pressure loads. With such loading situations, the response of a membrane to loading is fundamentally non-linear, when a combination of increasing stresses and increasing deformations carry an increasing pressure. This implies a need for smoothness and continuity in the deformed configuration of the membrane.
Rubber-like materials are commonly considered for membranes, and are then described by hyper-elastic material models. Due to computational convenience, the Mooney-Rivlin model, [35,43], is often used in a two-parameter form B Arne Nordmark nordmark@kth.se 1 KTH Engineering Mechanics, Royal Institute of Technology, Osquars backe 18, 100 44 Stockholm, Sweden related to the first two strain invariants [26]. Several studies have shown that the relation between the two constitutive parameters significantly affects the qualitative response of simulated membranes [14,40]. Many other material models for the simulation of membrane structures exist in literature, where a common one is given by Ogden [38]. Comparisons of models are given in, e.g., [8]. Parametric material models give large flexibility in describing materials, but can also cause unexpected results, e.g., non-intuitive instabilities [15].
Being thin structures, membranes are strongly affected by instabilities. Common views demand for static stability a minimum total potential energy at the equilibrium state [4,19,30,47]. Stability is, however, fundamentally dynamic, and can be evaluated from a Liapunov stability condition, i.e., a capacity to remain close to the static equilibrium situation after a minor disturbance [33]. A special aspect of membrane stability is wrinkling, when compressive forces create wavy transversal deformations with crest lines orthogonal to the principal compressive stress [18]. Relaxed strain energy forms allow the calculation of wrinkling as an average stress field, but lead to ill-defined problems under pressure loads [41].
For the stability of thin membranes, analytical results can be obtained only for simple geometries [23,37,45], but general treatments are also available [5,7,21]. Of particular interest has been the stability of pressurized spheres. This problem is frequently studied in literature [29,36,37,48,50, and many others]. An extensive treatment of this problem is given by Haughton and Ogden [24,25], where an axisymmetric analytical model and two hyperelastic material models provide detailed stability conclusions from the total potential energy in the deformed state. The first part of the present work can be seen as a review and an extension of this, with a more general kinematic description, a wider range of material models, and a stability analysis based on vibration properties [17]. As a motivation, the latter part of the work is also directed towards the representation of instabilities by numerical simulation models.
The stability analysis is in the present work restricted to the isolation of critical equilibrium states on the primary sequence, while secondary paths emanating from these are not a main objective. Both stable and unstable equilibria are thereby sought, as the treatment of bi-stable structures is considered as an interesting area for membranes [1,2,6].
The description of the geometry of a spherical surface is a well-known and frequently studied problem, not least in cartography, meteorology and oceanography [22,42,49, as examples]. When physical phenomena need be described on the spherical surface, the parameterized mapping of the surface causes problems [44], which can be handled by the Cubed sphere approach, with six local mappings together forming the sphere [46]. Another approach is to base the description of physics on the spherical harmonics functions [10]. The latter will be used in the present analysis, while the former is used as a tool to create one form of discretized mesh.
Nature shows many aspects of symmetry, and engineering frequently uses symmetric or repetitive designs, from both aesthetical and functional optimality reasons [11]. Symmetries of an object can be analyzed by group theory [9,12,13], where the elements of the symmetry group are the rigid transformations which bring a point on the sphere surface to another, leaving the center point unmoved, and keeping the spherical geometry [11].
In numerical simulations, symmetry is often introduced as mirror reflections in the main coordinate planes, or as rotational symmetry in axi-symmetric situations. A special aspect of symmetry is when an unsymmetric meshing of a fully symmetric domain can lead to surprising results [16]. As the present work utilizes flat triangular membrane elements to cover the complete sphere, this approximates the geometry by meshes of lower symmetry, which contain some, but not all, symmetry properties [28,51]. Classifying the resulting meshes with respect to their respective symmetry groups gives further insights into how a discretized mesh can significantly affect simulation results.
The analytical treatment of an isotropic hyper-elastic uniformly pressurized sphere is given in Sect. 2, discussing the combination of geometric properties and material models, and deriving stability conclusions from the vibration frequencies of the pressurized sphere. The relation between the spherical continuum and discretized element meshes is given in Sect. 3. Section 4 presents results from numerical experiments with finite element models, focussing on two common incompressible hyper-elastic material models, and on effects from the discretization. Section 5 draws a set of conclusions from the study. Appendices give further details.

Analysis of pressurized sphere
This section considers the linearized problem of small incremental motions for a spherical membrane of original radius a and uniform thickness t, when inflated to a radius λa by an internal constant over-pressure p; λ is then a uniform stretch. The membrane is assumed to be described by an isotropic hyper-elastic material model. For the dynamics, the inertia is assumed to come only from the membrane itself, with original uniform density ρ.
The motion around the equilibrium state is analyzed by finding all the eigenmodes of the membrane and the associated eigenvalues. The analysis is related to the primary equilibrium solution of uniform stretch, where critical solutions are found, with zero eigenvalues. At these, modes corresponding to vanishing eigenvalues exist, but secondary, non-spherical solutions are only very briefly mentioned here.
The discussion is related to two parts, where the first considers the conclusions that can be drawn from the geometry of the sphere alone, and the second considers the mechanics of a thin membrane formulation. These are brought together in a stability formulation, where stability of an equilibrium state is defined through the existence of vibration modes, i.e., positive eigenvalues for all modes of incremental displacements from the equilibrium state.

Spherical symmetry
The uniformly stretched sphere maintains the full spherical symmetry group O(3) of the original shape. The infinitedimensional space of all small incremental displacements from this state can be split into a direct sum of finitedimensional subspaces, each invariant under the spherical symmetry. For example, the 1-dimensional space of equal and purely radial displacement is clearly such a subspace. Group representation theory for O (3) classifies these subspaces as each corresponding to one of a list of known irreducible representations of the group, cf. "Appendix A".
On the other hand, the subspace of eigenmodes belonging to a particular eigenvalue is clearly also an invariant subspace under spherical symmetry. The end result is that the space of small incremental displacements can be split into a direct sum of invariant eigenmode subspaces, each of which corresponds to a particular irreducible representation of O(3). These eigenmode spaces have finite dimensions, and thus eigenvalue multiplicities, given by the representation. For the stretched sphere, the degree of symmetry is so high that the mode shapes are almost entirely determined by symmetry alone.
When considering all possible modes of the sphere, and not only axi-symmetric ones, a small incremental displacement vector for a point of the stretched sphere will be represented by a scalar outward normal component u together with a 2D vector v in the tangent space of the sphere. A basis for scalar functions is given by the set of real-valued scalar spherical harmonics Y m , where = 0, 1, 2, 3, . . . and m = − , − + 1, . . . , − 1, , cf. "Appendix A". If further ∇ (S) denotes the surface gradient of a scalar function on the stretched sphere, and * is the operation of rotating a tangent vector clockwise a right angle about the outward normal, then a basis for vector-valued functions on the stretched sphere is given by ∇ (S) Y m and * ∇ (S) Y m , where now = 1, 2, 3, . . ., as a zero value for gives a zero gradient.
This basis is convenient since it corresponds to the irreducible representations of the group O(3). The representations will here be denoted D g , D u for = 0, 1, 2, 3, . . .. In particular, for any given even , the set of 2 + 1 scalar functions Y m transforms as the representation D g , as does the set of vector functions ∇ (S) Y m (for > 0), whereas the set * ∇ (S) Y m (for > 0) transforms as D u . For odd , the first two sets instead transform as D u , while the third set transforms as D g . Since each eigenspace can be taken to correspond to one of the irreducible representations, the basis above is already close to an eigenmode basis, the exception being where a scalar and a vector set both transform according to the same representation. To summarise: The D 0g representation gives a radial mode where the mode amplitude α 0 is an arbitary scaling, so this mode is unique. While this mode shape is fully determined by symmetry, the corresponding eigenvalue depends on the physics of the problem. The representation D 0u would give u = 0, v = 0, which means there are no modes for this representation.
For ≥ 1, we first get modes corresponding to the representations D u for even , and D g for odd , of the form -a (2 + 1)-fold shear mode for each , again with an arbitrary scaling γ 1 and with a corresponding physics-dependent two different (2 + 1)-fold coupled modes (for i = 2, 3).
Here, physics determines not only the two eigenvalues for each , but also the ratio of the amplitude constants α i and β i . The relation between representations and mode types is summarised in Table 1.

Isotropic hyper-elastic membranes
The eigenproblem for solving the coefficients α i , β i , and γ i in Eqs. (2)-(3) needs the consideration of the physics of a particular problem on the sphere. This section gives the general equations for an isotropic hyper-elastic membrane, and specializes them for a uniformly stretched sphere. An isotropic hyper-elastic material is characterized by a strain energy density W e , i.e., per material volume. This is a function of invariants of the right Cauchy-Green tensor C. For a membrane, two coordinates X α , such as the two spherical angles θ and ϕ for the sphere, are introduced to describe points in its surface, while the third direction is taken to be its outwards normal. We will use the convention that lower case component indices run over the two in-plane directions, whereas an upper case index runs over all three directions. Using local plane stress, the covariant components of C then have the property C 3α = C α3 = 0, and C 33 can be solved for in terms of the in-plane components C αβ from either incompressibility or plane stress, cf. the material model examples below. Thus, W e can be expressed in terms of C αβ and by isotropy, from two invariants of the in-plane part C 2 contain-ing the components C αβ The two invariants allow the definition of the strain energy density as The general equation can be specialized for a spherical membrane of initial thickness t isotropically stretched by λ, for which J 1 = J 2 = λ 2 . This is expanded up to second order in small incremental displacements. As further shown in "Appendix B", the first order terms in both J 1 and J 2 are equal, but the second order terms are not. This means that to find the second order terms in variations in W e , it is enough to know the values of the fundamental material functions with t the initial thickness, and the arguments showing that the strain energy density derivatives should be evaluated for the arguments (J 1 = λ 2 , J 2 = λ 2 ). The description based on three functions was also used in [25], using the linearly transformed forms In Eq. (6), the expressions e, f , and g have the same physical dimension, force per length. For the unstretched material with λ = 1, f = −e for a material in equilibrium, and e and g must both be positive for a stable material model. It is also obvious that the expressions in Eq. (6) are easily obtained for any hyper-elastic material model which can be written in the form of Eq. (5).
The three functions in Eq. (6) express the dependence on a particular material model, given the particular stretch state existing in the pressurized sphere. It is not surprising that the derived functions are identical to the ones for a bi-axially isotropically stretched plane square given in [15], but with other expansions of J 1 , J 2 .
A few common material models are elaborated below.

Equation of motion for small deviations of the inflated sphere
Starting from the spherical geometry and symmetry, and introducing the mechanics of a uniformly stretched membrane, the equations of linearized motion around an equilibrium state are now presented. The uniform tension needed to keep the membrane in equilibrium at radius λa is related to the internal over-pressure p. As shown in "Appendix B", this can be related to the tension defined as the integral of the Cauchy normal stress over the current thickness as with the fundamental material functions from Eq. (6). The linearized equations of motion for the small incremental displacements described by u and v are derived in "Appendix B" as where (S) is the surface Laplacian, and the uniform mass per stretched area is The equations of motion are converted to an eigenvalue problem by replacing the second time derivatives through multiplication by −κ (if κ is positive, this corresponds to a harmonic time variation with angular frequency ω = √ κ), yielding one scalar and one vector equation

Eigensolutions to equations of motion
The eigenmodes from Eqs. (1), (2) or (3) can be used to compute the corresponding eigenvalues. Using the properties of the spherical harmonics Y m , cf. "Appendix A", the eigenvalues are dependent on the index , but not on m. As in "Appendix A", the short-hand notation q = Y m , r = (λa)∇ (S) q, s = * r is used in this section. First, considering the radial mode = 0 gives the single eigenmode u = α 0 q, v = 0 from Eq. (1). After some algebra, cf. "Appendix A", two equations are obtained, according to when dropping the λ arguments to g, h and ψ, and noting by the index on κ that the single eigenvalue belongs to = 0. Equation (15) gives, by introduction of Eq. (13), the eigenvalue corresponding to this mode. The eigenvalue is thereby related to the current stretch and the fundamental material functions in Eqs. (6) and (11), i.e., the material model and its parameters. Next, the case ≥ 1 with the first form of eigenmodes from Eqs. (2), is considered. Again, after some algebra, the eigenvalue equations in Eq. (14) then become -where a simplified notation is introduced. This gives the final expression for the eigenvalue now dependent on the mode order . As these eigenvalues κ 1 do not depend on m, they have multiplicities 2 + 1.
Since the corresponding eigenmodes have zero divergence, they were called pure shear modes above and correspond to purely circumferential incremental displacement in an axisymmetric model. As a special case, when = 1, i.e., θ = 0, the triple eigenvalue κ 1 1 is vanishing, corresponding to three modes of rigid rotation.
Finally, the coupled modes u = α 2,3 q, v = β 2,3 r from Eq. (3) are considered. Equation (14) then gives the two coupled equations where indices on κ, α, and β include both the order and the 2, 3 notation, which indicates that two new eigenvalues are obtained, each of multiplicity 2 + 1. They require the solution to the 2 × 2 eigenvalue problem which can be ensured to give real eigenvalues. For = 1, θ 1 = 0, the triple eigenvalue corresponds to three modes of rigid translation. Another triple mode is found from Eq. (21) for = 1 as The modes in Eq. (23) are commonly described as pearshaped, as they primarily consist of an axial deformation where one pole moves outwards, and the opposite one inwards.
For ≥ 2, the eigenvalues and eigenmodes can be computed by solving a second order algebraic equation from Eq. (21). It may, however, be easier to discuss the eigenvalues in terms of the trace and determinant of Eq. (21), which can be written By considering Eqs. (16), (19) for > 1, (23) for = 1, and Eqs. (24) for > 1, it is noted that all eigenvalues at a specific uniform stretch value λ can be easily evaluated through the order number , and the three fundamental material functions coming from the hyper-elastic model. This, and in particular when any of these eigenvalues vanishes, will be used for stability investigations below.

Wave-speed interpretation
As a further interpretation of the above expressions, it is seen that asymptotically for large , and thus large θ , the two equations (24) provide the solutions while Eq. (19) is still valid and can be expressed as The expressions for the eigenvalues are thereby all valid in high and arbitrary m, i.e., in the limit of short wave length, where curvature is unimportant. The three expressions correspond to transverse in-plane, transverse out-of plane, and longitudinal in-plane short waves with wave speeds (on the deformed sphere) respectively, and depending on the material model. The form of the equations is chosen to emphasize the similarity to common wave speed expressions for solids as the square root of a quotient of stiffness to density.

Stability
Rather than investigating the definiteness of the total potential energy, stability of an equilibrium state for a pressurized sphere is judged by the existence of vibration frequencies for linearized small motions around the state, considering the mass distribution of the model [17]. Stability thereby requires all eigenvalues to correspond to time harmonic modes, i.e., all κ i = ω 2 positive, except the two sets of triple zero eigenvalues κ 1 1 and κ 1 2 , which are related to rigid body motions. The signs of the eigenvalues from Eqs. (16), (19), (23) and (24) only depend on the stretch and the material model through the values of e(λ), ψ(λ), and g(λ). With e, ψ and g considered as independent quantities, the signs of the eigenvalues can be computed for any point in (e, ψ, g) space. This will give regions in (e, ψ, g) space corresponding to stability, and surfaces corresponding to critical stability with any vanishing eigenvalue. Once these regions and surfaces have been established, the increasing stretch λ of a sphere of any particular material will correspond to a curve in (e, ψ, g) space. These curves will pass through the regions and surfaces, thus giving the stability intervals and critical values for that material, in the particular situation of uniform stretch.
First, Eq. (27) is used to show that, with large enough , there are an infinite number of negative eigenvalues whenever any of e, ψ, or e + g is negative. Thus, any stable region must lie within a set A given by e ≥ 0, ψ ≥ 0, g ≥ −e. Within this set A, critical states are sought, for all .
The single mode for = 0 is critical according to Eq. (16) when The non-rigid-body mode for order = 1 is critical according to Eq. (23) when The critical surfaces for the 2,3 modes with ≥ 2 is given by the vanishing determinant in Eq. (24), which for increasing is As both e > 0 and ψ > 0 in the interior of the set A, the values g (e, ψ) form a decreasing sequence with for fixed e and ψ. As becomes large, g (e, ψ) converges to g ∞ (e, ψ) = −e. Thus, for given e > 0 and ψ > 0, a relation holds: The relation shows that the critical surfaces g = g (e, ψ) never intersect in the interior of A, but lie monotonously stacked below each other in (e, ψ, g) space for increasing , and bound by the simple surfaces g = ψ/2 and g = −e, as shown in Fig. 1. It can also be shown from Eq. (19) that no surfaces of type 1 exist in the interior of A.
The conclusion from Fig. 1 is that only the region above g = g 0 (e, ψ) is stable. For any point with e > 0, ψ = 0, and g > 0, e.g., the undeformed sphere, Eqs. (19) and (24) immediately show that all non-rigid-body eigenvalues are  (32), but also the critical surfaces e = 0, and ψ = 0. The quantity μt in the axis labels is the product of linearized shear stiffness and thickness stable The surface g = g 0 does, however, not belong to the stable region, nor does the boundary e = 0, as all modes of type 1 are critical there. This also shows that the more stable side of each of the surfaces g = g is the upper side, in the sense that the number of negative eigenvalue is increased by 2 + 1 on the lower side.
To conclude, the stable region consists of points where while there are n =0 (2 + 1) = (n + 1) 2 negative eigenvalues, all others being positive, in the region ψ > 0, e > 0, and g n+1 (e, ψ) < g < g n (e, ψ), (34) formed between the g n and g n+1 surfaces. Finally, an infinite number of eigenvalues are negative for This result agrees with the work in [25], except in the numbers of negative eigenvalues in the sub-regions, which are now higher due to the more comprehensive kinematic description. As a final remark, it is noted that these stability conclusions are valid for the case of given pressure. Other situations, for example controlling the pressure implicitly by prescribing the gas amount or gas volume are possible. As long as these indirect ways of specifying the pressure retain full spherical symmetry, the only eigenvalue that could be influenced is κ 0 , which would result in a different-or possibly absentlimit state surface g 0 . For example, for fixed volume, the surface g 0 disappears, and the stability region is bounded by g > g 1 = 0 instead.

Response to increasing stretch
The above treatment, and Fig. 1, have assumed that e, g and ψ are independent variables. Attention is now given to the curve in (e, ψ, g) space corresponding to a particular material when stretch λ increases from λ = 1 at a state with e(1), g(1) > 0, ψ(1) = 0. This belongs to the stable region, but is on the boundary, which defines the demands for an initially stable material. As the stretch λ increases from 1, the value ψ(λ) becomes positive, so the curve initially points into the stable region. Continuing on the curve with increasing λ, giving e(λ), g(λ), and ψ(λ), stability can be lost through three different mechanisms.
The first is related to g(λ) passing through g 0 (e(λ), ψ(λ)) from above, when the single critical mode is the spherically symmetric one, according to Eq. (1). This corresponds to a limit state, where the pressure p has a local maximum with respect to the stretch λ, cf. [25]. If g(λ) further decreases through g 1 (e(λ), ψ(λ)) = 0, a bifurcation with a 3-fold zero eigenvalue occurs for = 1. If g(λ) continues to move towards g ∞ (e(λ), ψ(λ)) = −e(λ), bifurcations occur, with an increasing number of unstable modes for = 2, 3, . . ., with increasing multiplicity. It is, however, fully possible that the function g(λ) could reduce instability by starting to pass back through the g (e(λ), ψ(λ)) values from below for increasing . It is even possible to regain stability at further increased stretch, cf. the Ogden model below.
A second mechanism for loss of stability is related to e(λ) becoming negative. As this, through Eq. (27) 1 , is related to shear modes according to Eq. (2), this makes all in-plane shear modes become unstable at the same value of the stretch λ, and the material thereby collapses in shear.
A third mechanism is related to ψ(λ) becoming negative. As this, through Eq. (27) 2 , is related to coupled modes according to Eq. (3), this makes all short wavelength out-ofplane modes unstable through a wrinkling behavior, which is also a collapse. Figure 2 shows how the e(λ), ψ(λ), and g(λ) functions vary during a uniform stretch variation λ for a few specific material models discussed below. The models show qualitatively very different responses, as clearly seen from the crossings with the critical surfaces and the vanishing e and ψ values. Note how the curve 'Case 1' will cross all critical surfaces before simultaneously reaching the g = −e and ψ = 0 surfaces.

Common hyper-elastic material models
The discussion above has considered a general hyper-elastic material model, where the strain energy density is written as a function of the two invariants of the in-plane strain tensor, Eq. (5). Stability conclusions can thereby be drawn for the pressurized sphere under uniform stretch, for three common material models.

Saint Venant-Kirchhoff model
The Saint Venant-Kirchhoff material model extends the linear elastic material model to the non-linear Green-Lagrange strains with capital letters denoting 3D indices (A, B = 1, 2, 3), and δ A B the Kronecker delta function. Introducing local planestress assumptions and formulating deviatoric strains gives the strain invariants with lower case letters denoting 2D indices, and overbars 2D deviatoric strainsĒ The strain energy density is then with G the shear modulus and K 2D = (9G K 3D )/(4G + 3K 3D ) the 2D plane stress version of the bulk modulus K 3D .
An example Since the constants G, K 2D > 0, the sphere is stable for all λ ≥ 1 as

Two-parameter incompressible Mooney-Rivlin model
A two-parameter Mooney-Rivlin model utilizes the incompressibility relation I 3 = 1, and defines a strain energy density as with two constants and the strain invariants, identified with the ones in Eq. (4), using that incompressibility gives C 3 3 = 1/(J 2 ) 2 . The strain energy density for the incompressible Mooney-Rivlin material model is thereby and the fundamental material functions become For a uniform stretch of λ = 1, this gives e(1) = 2t(c 1 + c 2 ), f (1) = −e(1), and g(1) = 3e(1), which identifies μ = 2(c 1 + c 2 ) as the linearized shear modulus at unit stretch.
As will be further shown by numerical examples below, the two-parameter Mooney-Rivlin material model can show three qualitatively different behaviors for positive c 1 , in addition to the sudden collapse for c 1 < 0: a monotonously increasing pressure under stretch, a pair of pressure limit states, or one pressure limit state followed by an infinite sequence of bifurcations. To what extent these instabilities can be reproduced by discretized models will be further discussed below.

Incompressible Ogden model
The strain energy density for the incompressible N -term Ogden material model is expressed using the three principal stretches λ J , [38], and 2N constitutive parameters μ i , α i .
A one-term restricted form of the model in Eq. (46) was studied in both [25,37] with conclusions in Table 2. Results are numerically confirmed below.
In two-term form, with four parameters, the Ogden material model can give a more complex behavior than the Mooney-Rivlin material model, and is then strongly dependent on the α i exponents. This will be further shown by numerical examples below. One example is Case 4 in Fig. 2.  (46), μ 2 = μ 3 = 0 All modes are stable at large stretch −3/2 < α 1 < −1 or 2 < α 1 < 3: Only the = 0 mode is unstable at large stretch 1 < α 1 < 2: A finite number = 0, 1, . . . , max modes are unstable at large stretch −1 < α 1 < 0 or 0 < α 1 < 1: Modes are unstable for all at large stretch

Symmetries of structures and models
Previous work [16,17] has shown how symmetry properties of a discretization can significantly affect the calculated response to loading. In discretized modelling, symmetry aspects will essentially appear in two contexts, where the first relates to a complete structure represented by a mesh which has, or does not have, particular symmetry properties. The second aspect is related to the modeling of a minimal part of the full structure, where sets of conditions are introduced on the subdomain boundaries [16]. For example, using a mesh with icosahedral symmetry, only a fundamental domain of 1/120 of the sphere need be modelled, but with a significant amount of bookkeeping needed to reconstruct the solutions and eigenmodes on the full sphere. The present work only considered modeling of the full sphere, and this section will discuss how the symmetry groups of a discretized mesh will affect the description of the critical eigenvectors of order on the sphere.

Splitting of representations
The representations D g and D u for = 0, 1, 2, . . . in Table 1 are irreducible for the full spherical symmetry group O(3). As noted above, the symmetry of a discretized mesh must be that of a finite subgroup G of O(3). While D g and D u are still representations for G, they are not necessarily irreducible, since the subset of group elements belonging to G also implies a subset of linear operators when testing for invariance. Like any reducible representation, it will split into a direct sum of irreducible representations of the group G.
"Appendix C" reviews how to use published character tables to compute the splitting. The conventional names for the irreducible representations follow a pattern, where A and B are one-dimensional representations, E twodimensional, T three-dimensional, G four-dimensional, and H five-dimensional. Numerical indices denote different representations of the same dimension, whereas the indices (·) g and (·) u denote symmetry and anti-symmetry with respect to inversions.
As an example, "Appendix C" shows that a mesh with tetrahedral symmetry will split the 11-dimensional representation D 5u into E+T 1 +2T 2 with dimensions 2+3+2·3 = 11. This shows that, for any reasonably fine mesh, a group of four close eigenvalues of multiplicities 2, 3, 3, 3 (in some order) will be observed instead of one eigenvalue of multiplicity 11. The same reasoning shows that the zero-crossings of the eigenvalues will not coincide.
For critical eigenvalues and modes, it is of special interest to identify whether a mode corresponds to the trivial representation 1 where each group element is associated with the one-dimensional identity operator. The eigenmode then has full symmetry, the same as the equilibrium solution. In general, a zero eigenvalue then shows a limit state on the primary branch. On the other hand, critical eigenmodes corresponding to non-trivial representations imply a symmetry-breaking bifurcation, where one or more secondary branches of lower symmetry can cross the primary branch.
Although several response aspects converge as expected with fineness of the meshing, some aspects like the capacity to reproduce bifurcations of higher orders are independent of the mesh refinement, and are purely effects from the mesh symmetry properties.

Meshes from regular geometries
A basic method to create a mesh on the sphere starts from a simple regular geometry, and refines it systematically. Only a few realistic instances exist, and the symmetry properties of three regular geometries are discussed below.

Icosahedron-based meshes I h
The best modelling of the sphere based on a triangle mesh, i.e., the one with highest symmetry is using the icosahedron symmetry group I h . The analysis of a mesh with I h symmetry results in the representation splittings in Table 3. Any mesh based on icosahedral symmetry will then represent the criticalities up to = 2 as isolated states with correct multiplicity, while the states for = 3, 4, 5 will be split into 2, 2, 3 closely spaced but separate bifurcations of lower orders, respectively. The presence of the trivial representation A g in the entry for = 6, corresponding to eigenvalues κ 6 2 and κ 6 3 , shows that the I h models will here give a limit state-instead of a bifurcation-and thereby fail to follow the correct sequence of uniform stretch.
The trivial representation A g is highlighted Table 4 Splitting of out-of-plane and in-plane modes by a T d (tetrahedral) mesh, cf. Table 1 D g D u The trivial representation A 1 is highlighted

Tetrahedron-based meshes T d
Conclusions can be reached similarly for discretizations based on the tetrahedron symmetry group T d , Table 4. A tetrahedron-based mesh will therefore show vanishing eigenvalues of at most multiplicity 3. A conclusion is also that any mesh based on tetrahedral symmetry will split the = 2 bifurcation into two bifurcations of multiplicities 2 and 3, respectively, and will fail to represent the = 3 bifurcation state correctly.

Octahedron-based meshes O h
The analysis of meshes with O h symmetry is similar, and shows that any such mesh will split the = 2 and = 3 bifurcations into isolated bifurcations of multiplicities 2 + 3 and 1+3+3, respectively, and will fail to represent the = 4 bifurcation state correctly. It is noted that a cube, which could be another possibility to create the mesh, is also of symmetry group O h .

Axis-based meshes
Symmetry aspects are in engineering practice most commonly observed through mirror planes. This corresponds to the creation of a well-defined subdomain, followed by a rep-etition of this through simple geometric transformations. In the present context, it is rather natural to introduce an equator mirror plane, and to define as the basic subdomain one sector, limited by mirror planes of constant longitude, of the half-sphere. This introduces an axial direction, around which the initial mesh is rotated. Consequently, the meshes are not equal in the three global axis directions, and are of lower symmetry. With n sectors around the circumference, the symmetry will be from the dihedral group D nh . In the present work, a basic mesh was created for n = 8, i.e., from one eighth of the half-sphere, introducing eight vertices in the mirror plane, together with one apex point on each side. The symmetry group for these meshes was thereby D 8h .

Dihedral meshes D 8h
An analysis of a mesh with D 8h symmetry will give results of the same form as above, and indicate that any mesh based on the D 8h symmetry splits already the = 1 bifurcation into two bifurcations of multiplicities 1 and 2, respectively, and will fail to represent the = 2 bifurcation state correctly.

Systematic mesh refinement
With either of the two approaches above, a minimal basic mesh is created, consisting of a set of triangles. This mesh introduces the symmetry properties which will also be valid for recursive refinements to computational meshes. The basic meshes used here were created as a set of identical surfaces, but it is noted that meshes created from regular geometries consisted of equilateral triangular surfaces, while the axisbased mesh consisted of isosceles triangles.
A requirement for kept symmetry is a systematic procedure for the refinement from the basic mesh. Aiming at a general method for any domain, a recursive process divided each existing triangle into four, by introducing nodes at the midpoints of existing triangle edges. All new nodes were then radially moved to the unit sphere, before going to the next recursive step. This is essentially similar to the approaches in [3,42]. It is noted that, although all vertices are correctly placed on the sphere, the recursive introduction of new vertices will lead to a hierarchy of nodes Two steps of refinement of an equilateral triangle in Fig. 3 shows how new generations of nodes are created, with different positions in relation to the symmetry of the initial mesh; nodes with the same notation are symmetrically identical in the uniform stretch. Figure 4 shows how the initial regular tetrahedron with four flat surfaces and four vertices was refined to a T d mesh with 1024 surfaces and 514 vertices, keeping the initial symmetry.
As discussed in [27], this method for meshing on the sphere is not perfect, as the sphere surface will be divided into elements of unequal sizes, with larger elements in the central  regions of initial triangles. The effect is particularly visible for the tetrahedron-based meshes in Fig. 4. For the present purposes, and primarily focussing on the icosahedron-based meshes, the effects from non-uniform meshing are limited, as long as symmetry is kept. Another non-iterative method for mesh generation is similar to the method shown in [46]. A cube is then regularly but non-uniformly divided into six side meshes. When the cube sides are mapped onto a circumscribed sphere, this mesh becomes slightly more uniform in element sizes.
A practical aspect of mesh refinement is that very high precision must be used, in order to keep the intended symmetry; a relative inaccuracy in nodal coordinates of ≈ 1 · 10 −6 was in [17] shown to destroy symmetry in results.

Meshes for numerical tests
Simulations of the gas-pressurized spherical membrane were primarily performed with meshes of triangular elements from four different symmetry classes. These meshes are below denoted using the symmetry group name followed by the number of elements within parentheses. The four classes were: -A regular icosahedron mesh I h (20) (16) mesh was created according to the description above, and using 10 nodes. The used D 8h (4096) mesh with 2050 nodes needed 4 refinement steps, which gave 81 different nodal classes, due to the isosceles initial triangles.
In addition, a variation of the cubed sphere mesh from [46] was created. The mesh was created from an internally symmetric mesh of 2N × N (N even) triangular elements on each side of a cube, with the nodes then moved radially to correct radius. In particular, a C(3072) mesh was created with N = 16 giving 1538 nodes of 45 different nodal classes. The resulting cube-based mesh maintains the O h symmetry of the cube.

Numerical verification
The previous sections have given an analytical treatment of the stability of gas-pressurized hyper-elastic spheres, and discussed how different discretizations for the continuous sphere will affect the obtained results. The present section will demonstrate how the analytical results and predictions are reproduced in element-based discretizations. Results are presented for parameter variations within two hyper-elastic material models discussed above, and for different meshes.

Problem setting, and basic formulation
The case studied was a thin sphere of midplane radius a = 50 mm and thickness t = 0.05 mm pressurized by an internal over-pressure p. Incompressible material models were considered, with linearized shear stiffness for the unstressed membrane μ = 422.5 kPa. The material models were included in the formulation of a faceted triangular thin membrane element with linear interpolation of translations and local plane stress conditions [14]. No handling of wrinkling was included, but all presented results were verified to be wrinkling-free.
The numerical tests used meshes described in Sect. 3.5 above, with the I h (5120) mesh as default, and with one initial node at (0, 0, a). The radial displacement u a of this node is used to quantify stretch as λ = (a + u a )/a in results.
Mainly solutions on the primary solution path for the discretized problem were obtained, i.e., the path where solutions maintain the full symmetry of the mesh. Results are represented by the relation between stretch λ and over-pressure p. The non-linear equilibrium problem was set up as described in [17]. Stability of equilibrium states was evaluated from the Jacobian matrix, with no inertia associated to constraint equations. Bisection was used to high precision, in order to isolate states giving a change of stability. Evaluation of secondary paths was only occasionally performed, even if this could be done in the present setting.

Mooney-Rivlin models
Simulations were performed for Mooney-Rivlin models with different values for k = c 2 /c 1 in μ = 2(c 1 + c 2 ). Results are given in Fig. 5, as pressure-stretch curves for six values 0.18 ≤ k ≤ 0.23, and a critical sequence for a larger range.
The critical sequence verifies that the double limit states only exist for k ≤ 0.21446. To given accuracy, the same limiting value for k was obtained with the coarser I h (1280) and I h (320) models, but the stretch at this limiting value was slightly different for the three meshes. Negative constants For a constitutive constant k < 0, also bifurcations appeared for the Mooney-Rivlin models [17]. The first bifurcations were well represented by the I h (5120) model. For c 1 > 0, c 2 /c 1 = −1/2, the limit state was found at λ = 1.1441, and the two first bifurcations at λ = 1.1606 and λ = 1.1871, respectively, cf. Sect. 2.5.2. Equilibrium sequences obtained for a few different negative k constants are shown in Fig. 6. A magnified view on the case k = −0.2 is also given in Fig. 6b. The figure shows how the = 5 bifurcation was split into three distinct bifurcation states for this model, while the = 6 bifurcation was incorrectly represented. Higher order bifurcations > 6, although existing in the analytical solution, were not reached (but could be evaluated for a sequence restarted at larger stretch). The pressure-stretch relation shows how the solution deviates from the correct relation for the continuous sphere. This deviation occurs

Ogden models
Simulations were also performed with an Ogden model for different constitutive parameters. Initially, these were based on the 'well-known' six parameters discussed above, aiming at a qualitative description of the influences from the constants μ i , (i = 1, 2, 3), but keeping the constants α i . Results are given in Fig. 8, showing also cases where μ 2 = 0, μ 3 = 0 or μ 2 = μ 3 = 0. The results in Fig. 8 show qualitatively that the two additional terms gave similar effects on the overall response when compared to Case 4 with only μ 1 non-zero, as the incompressibility condition makes either the in-plane or the transversal stretches above 1, with the other below 1. The examples below only considered at most two terms in the Ogden formulation.
It is noted in Fig. 8 that bifurcations appeared for all cases. For Case 1, the two limit states were found at λ = 1.3733 and λ = 4.3179, respectively, while bifurcations appeared for λ = 1.7739 and λ = 2.5383, which agrees perfectly with the analysis in Sect. 2.5.3, and well with numerical results in [48].
One-term Ogden models Tests were then performed for a one-term Ogden model, where only the exponent α 1 was varied, with μ 1 = 2μ/α 1 . Results are given in Fig. 9, and verify that large positive and negative exponents gave monotonous responses, but that values closer to zero gave more complex responses, cf. Sect. 2.4.1. Further experiments showed that values 0 < |α 1 | < 1 gave an infinite number of bifurcation states for the continuous case, but that only a finite set of these were shown by the discretization, while 1 < α 1 < 2 gave a limited number of bifurcations. For α 1 = 1.5, only the = 1 bifurcation was found, while α 1 = 1.3 in Fig. 8, Case 4, gave the = 1, 2 cases, and a further case with α 1 = 1.1 was correct up to = 3 (split into one 4-fold and one 3-fold bifurcation). The α 1 = ±0.5 cases both showed the = 5 bifurcations, before the I h (5120) model was unable to follow the continuous solution, cf. Fig. 7.

Two-term Ogden models
More complex behaviors can result for cases with more than one term in the Ogden material model, and interest was focussed on two-term models.
Based on Fig. 9, parameters were chosen according to α 1 = 1.05, α 2 = −2, with μ 2 = −k μ μ 1 , μ = (μ 1 α 1 + μ 2 α 2 )/2. Some results are given in Fig. 10, showing that even for a very small positive value for k μ (i.e., a negative μ 2 together with the negative α 2 ), the response eventually regained stability. While the case k μ = 0 (not drawn) passed through the limit pressure state = 0, and then successively up to = 5 (but no higher) for increasing stretch, a case with k μ = 10 −5 passed = 5, and then successively went back through lower values until a stable solution was again reached at λ = 62.09, from which stretch only positive eigenvalues were present. For k μ = 0.1, only the two limit pressure states were found, and for all k μ > 0.17694, the sphere was always stable, as seen from a critical equilibrium sequence, similar to the one in Fig. 5. An evaluation for the = 1 bifurcation showed that this 3-fold critical state only exists for k μ ≤ 0.03053.

Discretization
Results above were all obtained with the I h (5120) model, which is considered to be an accurate representation of the continuous sphere. Even this discretization, however, has shortcomings.
The main aspects coming from the discretization of the continuous sphere into a finite number of flat linearly inter-  Fig. 6, showing how the bifurcations for ≥ 3 were distorted into either a splitting of a high order bifurcation (with many simultaneously vanishing eigenvalues) into several lower order ones, or an inability to represent a bifurcation existing for the continuous sphere. These aspects were studied for different discretizations, with a focus on the symmetry properties and the fineness of the meshes.

Fineness of I h model
Referring to Fig. 6, showing result from an I h (5120) model with a Mooney-Rivlin material with k = c 2 /c 1 = −0.2, icosahedron meshes of different finenesses were analysed, with respect to critical states. Results are given in Table 5, as the number of vanishing eigenvalues and the corresponding stretches λ. The four finer meshes show the same orders and splittings of the bifurcations for ≤ 5, and the same failure to represent the bifurcation at = 6, where analytical results show that the critical modes should be of multiplic-

Symmetry properties of the mesh
The case with a two-term Ogden model for k μ = 10 −5 shown in Fig. 10 was used to compare the effects of meshes with different symmetry properties. For the I h (5120) mesh, critical states up to = 5, and then down to = 0 were found, Fig. 11. No other meshes were able to reproduce the critical states up to order = 5, i.e., the two 11-fold bifurcations. Not only did the models with lower symmetry fail to reach the = 5 bifurcation, but they also split the lower order critical states in different ways, cf. Sect. 3.1. Table 6 shows how the critical states for increasing were represented as one or more critical states of lower multiplicities.
It is noted from the table that the = 4 bifurcation was split into another sequence, 4 + 5 rather than 5 + 4, compared to the Mooney-Rivlin model reported in Table 5, showing that the splitting sequence order is both problem and mesh dependent. It is also noted that the cubed sphere mesh C(3072), as expected, showed a similar behavior as the octahedral O h mesh, even if the deviation from the purely spherical stretch took another direction, Fig. 11. Figure 12 shows how the original symmetry of the meshes was kept during stretch, even when the results in Fig. 11 seem to indicate that the fundamental paths were left by the simulations. The subfigures show, for the icosahedral and the tetrahedral meshes, how all nodal radii of the models developed with increasing stretch. The figures verify that the meshes consisted of 30 and 102 nodal classes, respectively. The response curves in Fig. 12 are therefore primary solution sequences for the discretized meshes, and no symmetrybreaking bifurcations were introduced in the solutions. The imperfect spherical stretch is thus due to the inability of the discretizations to represent all displacement modes of the sphere.
From the figures can be noted the significantly different scales of the radial values, and also that the stretch of the sphere was outwards between the initial icosahedral vertices, while it was inwards for the tetrahedral mesh. This, however, might be different for other meshes with the same symmetry. The T d model simulation reported in Fig. 12b returned to an almost perfect spherical shape after a very severe deviation, where equatorial radius was less than 10% of the polar for stretches around 10.
As noted above, the symmetry of a mesh can be very sensitive to disturbances. A practical consequence relates to the effects of rounding nodal coordinates. In the I h (20) mesh of the unit sphere, the nodal coordinates consisted of six different values (except 0 and 1), with either sign. In the simulations presented above, these values were given with machine precision. Rounding all these values to a lower number of decimal places would keep a certain regularity in nodal placement, but break the I h symmetry. It was by simulations verified that a rounding of all unit sphere coordinate values to six decimals before refinement was enough to destroy the correct instability behavior.

Secondary paths
The present work has focussed on the presence, isolation and identification of critical states along the primary solution sequence for the pressurized sphere, noting the difference between analytical and discretized models. Obviously, the bifurcation states will be origins for secondary branches. For bifurcations of multiplicities higher than one, several secondary paths will normally exist, and group theory can be used to provide a list of such branches [20]. As an example, the = 1 bifurcation, with the axi-symmetric 'pear-shaped' mode, gives a 3-dimensional critical eigenspace of zero eigenvalues, expressed by a basis of three arbitrary but orthogonal eigenvectors. While an infinite number of rotation axes exist for the continuous sphere, the discretized models will only be able to yield pear-shaped modes along certain symmetry axes of the mesh. As all I h meshes have 6 axes through initial icosahedron vertices, 15 axes through midpoints of the initial edges and 10 axes through initial side centers, each of these 31 symmetry axes was verified to define a secondary sequence, which could be easily followed between the bifurcations. Due to the nodal hierarchy, the three types of sequences were marginally different in their representations of the spherical mode. The T d meshes similarly give four axes through a vertex and a side center, and three axes through two edge midpoints, giving a total of 7 secondary paths passing through the = 1 bifurcation, which is correctly reproduced. It is also noted that these conclusions are completely independent of the lev- Table 6 Splitting of critical states for meshes of different symmetry properties Order The numbers indicate the multiplicities of isolated critical states The higher order bifurcations lead to even more complex secondary path situations.

Conclusions
The paper has discussed the numerical evaluation of stability for hyper-elastic membranes, with the pressurized sphere as an example, and focussed on the dependence on the material model and the discretization of the complete geometry.
The material description has developed a general framework for hyper-elastic material models, where most common models can be introduced, the only demand being that a strain energy density can be formulated from first and second invariants of in-plane strain. The derivation shows that three fundamental material functions describe the material response for the case of uniform stretch of the sphere. From these, considering the case of infinitely small stretch, conditions on the material parameters for initial stability can be formulated. Stability conditions for the stretched equilibrium can also be expressed in these functions, for the particular material model. For common material models, the possible instability responses are thereby considerably more complex than previously found, in particular for cases where wider ranges of material parameters can be allowed within the requirements of initial stability. For instance, two-parameter Mooney-Rivlin models with one negative constitutive constant can show an infinite number of bifurcation states in the response to pressurization, whereas a two-term Ogden material model can give a response with a high but limited number of bifurcations followed by a return to stable conditions at larger stretches. The developed expressions allow a straight-forward evaluation of critical states.
For the continuous modelling of the sphere by spherical harmonics functions, an infinite series of such functions can completely represent the incremental deformations from the deformed sphere. The eigenmodes, representing either vibration or instability modes, are constructed by vector analysis from the spherical harmonics functions, introducing the fundamental material functions.
Symmetry in representation of the spherical membrane is a key aspect in discretization. The main aspects deciding what instabilities can be correctly identified by a discretized element model of the complete sphere are described through group theoretical concepts. One main conclusion is that meshes derived from a basic icosahedral geometry will give the best possibilities to represent the instabilities of a pressurized sphere, but that some limitations exist also for these.
Numerical experiments were used to verify analytical results, and in particular to show how variations in material and geometry are represented by the modelling. Performed simulations thereby accurately verified the analytical predictions of the critical sequences for the Mooney-Rivlin and Ogden material models, showing also the higher order bifurcations existing for certain sets of constitutive parameters. The numerical simulations allowed the prediction of parameter ranges where specific instabilities exist. Based on triangular flat elements on the spherical surface, simulations also showed that symmetry properties of the used mesh will affect the obtained results in a manner predicted by the group theoretical analyses. The main conclusion is that symmetry aspects of a studied problem must be carefully and very accurately considered if a full instability investigation of a structure is the aim. In particular is shown how simplified considerations of symmetries can be misleading in this context. The results also clearly show that a lack of symmetry in the mesh can not in all aspects be compensated for by a finer mesh.
Acknowledgements Open access funding provided by Royal Institute of Technology. Financial support for this project from the Swedish Research Council (Dnr: 2014-5686) is gratefully acknowledged.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

A.1 Group representations
A brief review of how representations naturally come out of the eigenvalue problem follows. Each transformation from the group O(3) maps the stretched sphere onto itself, but also maps any small incremental displacement to a new small incremental displacement, and this latter mapping is linear. Each element of the group is thereby associated with a linear mapping on the space of small incremental displacements. Further, two elements of the group successively applied will correspond to composition of the linear operators. To conclude, each group element g is associated with a linear operator L g , and composition is preserved in that L g•h = L g L h . This composition-preserving association between a group element and a linear operator is termed a representation of the group.
This representation is reducible in the sense that there are subspaces of the tangent space that are invariant under all of the L g . One example is the one-dimensional subspace of all radial and equal displacements given by a small uniform change in radius, which is invariant under full spherical symmetry.
The infinite-dimensional full space of small incremental displacements will split into a direct sum of finitedimensional invariant subspaces. On each such subspace of dimension n, the association of a group element g with the restriction of the linear operator L g is an irreducible representation of the group. Given a set of n basis vectors for the subspace, each linear operator will be represented by a square n × n matrix, and the basis can be taken such that the matrix is orthogonal for every operator. In such an irreducible representation, group compositions thus correspond to matrix multiplication: M g•h = M g M h , with no further splitting then possible. This splitting of the small incremental displacement space into invariant finite-dimensional subspaces and corresponding representations also splits the fully symmetric eigenproblem, as any eigenmode must be transformed into a new eigenmode with the same eigenvalue under the group operations. Thus, any invariant subspace must correspond to a single eigenvalue, which will have a multiplicity at least equal to the dimension of the subspace. In general, one can expect different subspaces to have different eigenvalues, with multiplicities precisely the dimensions of the subspaces. Knowing the irreducible representations for a group thus gives direct information about the multiplicity of eigenvalues and eigenmode shapes. Complete lists of irreducible representations are known and published for all groups considered here, cf. "Appendix C".

A.2 Representations and eigenmodes for the spherical membrane
Since  (49) where i = 1, 2, 3 denotes the three solutions to a physicsand -dependent 3 × 3 eigenvalue problem for the amplitude constants α, β, and γ . Also, by varying m, 2 + 1 basis eigenmodes with the same eigenvalue exist for given ( , i).
Without the consideration of physics, only the mode shape corresponding to = 0 is fully determined. For higher values of , the amplitudes given by α, β, and γ needs the actual physics to be determined, and will give three possibilities for each from a 3 × 3 matrix eigenvalue problem.
These three can be separated into two groups of one and two, respectively, by next considering the inversion symmetry group C i . Here there are two irreducible representations.
In the representation named A g , inversion leaves the shape unchanged, while in the representation named A u , inversion causes a sign change. The spherical harmonics can be expressed as homogeneous polynomials of degree in the variables x, y, and z. Since inversion means changing the sign of all coordinates, it is easy to see that for even , both Y m and ∇ (S) Y m are unchanged by the inversion. The * operator on the other hand will change from a rotation clockwise by a right angle to a rotation counter-clockwise. Since rotating a vector counter-clockwise by a right angle is the same as rotating the vector clockwise followed by a multiplication by −1, * ∇ (S) Y m changes sign under inversion if is even. For odd , the sign changes follow the opposite pattern. This means that the eigenmodes in Eq. (49) are obtained by taking either γ = 0 or α = β = 0.
As the modes in Eq. (2) are fully described for any , the symmetry considerations have reduced the eigenmode problem to a set of 2 × 2 problems for the amplitudes α i and β i in Eq. (3), giving two possibilities i = 2, 3 for each . Formulation and computation of the eigenvalues for the then known mode shapes requires the actual equations of motion, i.e., the physics of the problem on the sphere geometry.

A.3 Scalar spherical harmonics
The scalar spherical harmonic functions are the basis for the symmetry consideration above. The Laplace spherical harmonics functions Y m constitute a complete set of orthogonal functions defined on the surface of a sphere. They may be organized by angular frequency, and are basis functions for representations of functions on the sphere, somewhat similar to the functions in a Fourier series. The shapes of the used functions can be roughly visualized as a two-dimensional product of trigonometric functions with crests and node lines, corresponding to a mesh with −|m|+1 half-waves between the north and south poles, and 2|m| half-waves along the latitudes [10]. As an example of how each sperical harmonic function generates mode shapes for small incremental displacements of the sphere, we can consider the case of > 0, m = 0, where Y 0 is independent of the longitude angle. Then the shear mode displacements are along curves of constant latitude and radius, whereas the coupled mode displacements lie in planes of constant longitude, and all displacements are axially symmetric.

A.3.1 Function and operator expressions
In the further description below, a number of differential operator expressions on the spherical harmonics functions are used. The operators are all defined on the stretched sphere. In particular, all operators can be constructed from the surface gradient ∇ (S) and the quarter turn rotation operator * , in combination with the surface divergence operator ∇ (S) •.
As these are homogeneous and isotropic operators on the stretched sphere, so are all operators constructed from them.
For the operator expressions used, it can first be shown that for a general scalar function d and Further, for a general vector function e * * e = −e, The spherical harmonics functions also have the fundamental property From these basic expressions, it is straight-forward to obtain a set of relations for scalar and vector functions based on the spherical harmonics. For a scalar function one obtains, based on Eqs. (51) and (54), -the operator gives a scalar multiplication for the argument. For vectors of the format and and from Eqs. (58), (57) and (50) Finally, Eqs. (53), (59), (60) and (57) give while Eqs. (53), (59), (60) and (58) give -again with results which are scalar multiplications of the arguments, with the scalar dependent on the mode order and the stretch value λ.

B Derivation of the equations of motion
This appendix derives the equations of linearized small motion around a uniformly stretched pressurized spherical membrane, using the fundamental material functions for a hyper-elastic material model.

B.1 Small displacements of a curved surface
We first review some basic notions concerning the differential geometry of surfaces. Consider a surface given by R(X ), where X α (α = 1, 2) is an arbitrary coordinate system for the surface. Let N(X ) be a unit normal vector to the surface. The normal direction should be consistently chosen so it varies continuously. The shape (up to rigid motions) of the surface is completely described by two tensors. The first fundamental form, or metric tensor, is where ,α denotes partial derivation with respect to X α . The second fundamental form is where the second expression shows that L αβ is symmetric. While the metric deals with properties of the 2D space itself, the L tensor deals with how the space is imbedded in 3D. As with any space with a metric tensor, the metric tensor and its inverse g αβ can be used to raise and lower tensor indices, Christoffel symbols Γ α βγ can be computed and thus the covariant derivative ;α can be defined. Now, let us at each point X take a vector v α (X ) in the tangent plane of the surface and a scalar u(X ), and use these to define a finite displacement of the surface in two steps. First, we follow a geodesic starting at X and with initial tangent vector v for a unit value of the affine parameter, to reach a point in the surface with coordinates x α (X , v(X )). For this point we move along the normal vector N(x) with coefficient u(X ) (note that this is u as a function of the original point X , not the new point x). Thus, the displaced point is The right Cauchy-Green tensor is and C is needed up to second order in v and u to obtain the linearized equations of motion about a stationary solution.
The result is, with terms of third and higher orders neglected, where R αi jβ is the Riemann curvature tensor, and everything is evaluated at the point X . Raising one index to C α β = g αi C iβ , we can compute the two invariants where R αβ = R i αβi is the Ricci tensor, which may alternatively be computed as R αβ = −g αβ Det(L i j ). Alternatively, we can express J 2 as where both terms involving only v are divergences.

B.2 The elastic energy for small displacements of a an unstressed membrane
Consider an elastic membrane where the undeformed reference configuration is without any tension or strain. The elastic energy per unit material volume for an isotropic material under plane stress conditions is a function W e (J 1 , J 2 ) of J 1 and J 2 . The energy density per unit material area is t W e , where t is the membrane thickness. Using the definitions of the fundamental material functions we have for deviations from the undeformed material (where λ = 1), up to first order Since the material is in equilibrium in the undeformed configuration, the variation of the first order terms must vanish, so f (1) = −e(1). Using this, the energy can be compactly expressed to second order using the linear strain tensor and the deviatoric linear strain tensor which gives, to second order, Thus, e(1) can be interpreted as the 2D shear modulus and g(1) as the 2D bulk modulus integrated over thickness, giving the contributions from area preserving shear deformation, and shear-free area deformations, respectively.

B.3 The elastic energy for small displacements of a uniformly stretched surface
Now we will try to relate the expressions derived for small displacements of an undeformed membrane, to the corresponding expressions for small displacements of an already deformed membrane, where stretches are different from 1. The same coordinates can still be used for the deformed membrane. We will use superscript (M) (for material) to denote quantities computed in the reference configuration, and superscript (S) (for spatial) or no superscript for the deformed configuration. The values of the Cauchy-Green tensor only depends on the coordinates, so C For the special case of a surface deformed using a uniform isotropic stretch λ, the difference between material and spatial quantities is a simple scaling involving λ. For example g (M)αβ = λ 2 g (S)αβ , which gives The elastic energy per unit material area is, to second order in u, v, where in the last term we have used that J (S) 2 = J (S) 1 to first order, and J (S) 1 = 1 to zeroth order. If we instead compute the elastic energy per unit spatial area and ignore the constant term, we get While we must have e(1) = − f (1) for any material, so the material is in equilibrium when unstressed and unstrained, there is no reason why e and f should be equal for λ = 1. Thus, the first order part of the elastic energy is not zero, and some loading is required to keep the membrane in equilibrium. We will use an internal over-pressure to supply the needed loading.

B.4 Over-pressure loading on a membrane
For a closed membrane subjected to a constant over-pressure p from the inside, we can define a pressure energy − p(V − V 0 ), where V is the enclosed volume, and V 0 a reference volume. To relate a small displacement u, v to a volume change, we note that for any v, r(X , v(X ), ξ ) defines a 3D coordinate system (X , ξ) for points near the surface. We can then find the volume change per unit area by integrating the Jacobi determinant of the coordinate system from ξ = 0 to ξ = u(X ), to get to second order in v, u. The pressure energy per unit area is thus The same expression can actually be used for both a closed membrane, and an orientable membrane with fixed boundaries, since the surface can be closed by adding suitable fixed portions.
Returning to the uniformly stretched membrane, we have the total potential energy per unit spatial area Φ = Φ e + Φ p and variations of the terms linear in u, v must vanish in equilibrium when integrated over the surface. To linear order, this gives Now, any term which is a divergence can be converted to an integral over the boundary, and the contribution to the variation will be zero since u, v are fixed on the (possibly non-existing) boundary. Thus, terms that are divergences can be ignored. Here, the term involving v i ;i is of that type, since e, f are constant. Then, the u terms give for the overpressure Since p is constant, so must L i i (twice the mean curvature) be. We may then measure the loading using the constant surface tension per unit spatial length, instead of the pressure p. If f is replaced by ψ − e and terms in the energy coming from e, g, and ψ are collected, we get Ignoring terms that are pure divergences (of which the second expression for J 2 contains several), we conclude that Φ = e(λ)¯ i j¯ j i + g(λ) i i j j /2 +ψ(λ) g i j u ;i u ; j /2 − L i j L j i u 2 /2 .
The expression for Φ has formally the structure of small deviations of a linear, isotropic 2D material in equilibrium with 2D shear and bulk moduli e(λ) and g(λ), cf. Eq. (72), with an extra stiffness coming from the tension ψ load, but it should be remembered that the ψ contribution includes some terms from the elastic energy.
In vector notation, using ∇ (S) for spatial surface gradient, (S) for the spatial surface Laplacian, and v for a vector tangent to the surface, this can be written

B.6 Kinetic energy contribution
We assume that all kinetic energy comes from the membrane itself, not from the motion of the fluid providing the pressure. If the constant density in reference configuration is ρ, the kinetic energy per material unit area is For the uniformly stretched membrane, the kinetic energy density per spatial unit area is using the mass density per spatial unit area h = tρ/λ 2 . Variation, and using partial integration in time to remove the time derivative from the variation of u and v, gives or, in vector notation, Combining the contributions to the variation from kinetic and potential energy, the equations of motion for small displacements of a uniformly stretched sphere are found as (89)

C Using character tables to compute representation splitting caused by lowering symmetry
When symmetry is lowered by restricting to a subgroup, any irreducible representation of the supergroup is typically a reducible representation of the subgroup, i.e., it can be split into a direct sum of irreducible representations of the subgroup. For the spherical membrane, the symmetry group O(3) of the continuous model is lowered to that of a subgroup like T d by the introduction of a discretized mesh. Consequently, an irreducible representation like the 11-dimensional D 5u of O(3) will be split into a direct sum of the irreducible representations of T d : the one-dimensional A 1 , A 2 , the two-dimensional E, and the three-dimensional T 1 , T 2 . This splitting can be easily computed using published character tables for the two groups [12]. The character χ α for a representation α of a group is a complex valued function over the elements of the group. Specifically, if the representation α maps the group element g to the orthogonal matrix R(g), then χ α (g) = Tr (R(g)) , which can be shown to be independent of which particular orthogonal matrix was used. Further, all group elements in the same conjugate class have the same value of the character. Thus, the character tables typically would be rectangular, with rows corresponding to different representations and columns to different conjugate classes of the group. Any conjugate class of a subgroup thereby corresponds to a unique conjugate class of the supergroup, but not the other way around.
The formula for computation of the splitting is given as where N is the number of elements in the subgroup, α is a possibly reducible representation, β an irreducible representation for the subgroup, and the overbar denotes complex conjugate. The number I (α, β) will be a non-negative integer, denoting the number of times that β occurs in the splitting of α.
As an example, the number of times the representation T 2 occurs when splitting the representation D 5u is computed. The representation D 5u is reducible for T d , so we will typically not find it tabulated. To compute the character values of D 5u , we note that O(3) is the direct product of the rotation group SO(3) and the inversion group C i , which means that the character values of D 5u are the products of the corresponding values of the rotation representation D 5 and the inversion representation A u . Since SO(3) is a continuous group, the character values of its representations can not be tabulated in a finite table, but a formula can be found as χ D l (g) = l m=−l cos(mφ).
(92) Here, φ is the angle of the rotation corresponding to the element g. The representation A u has the character values 1 for the identity element and −1 for the inversion. Thus, to compute the character values of D 5u for elements in T d , the corresponding rotation angle is found, and also whether it is an inversion or not. The character values of the rotation and inversion parts are multiplied. The conjugate classes of T d can be denoted: E (the identity element, no rotation, no inversion), C 3 (rotation by φ = 2π/3, no inversion, 8 different axes), C 2 (rotation by φ = π , no inversion, 3 different axes), S 4 (rotoreflection by a quarter turn, φ = π/2, inversion, 6 different axes), and σ d (reflection, φ = π , inversion, 6 different axes), giving a total of 24 elements. Taking an element g belonging to the class S 4 as an example gives Similarly, the other conjugate classes can be computed to give the first row of Table 7.
For the irreducible representations, the character tables are published in several sources. The rows for A 1 and T 2 are included in the table as examples. From this, the number of times T 2 occurs when D 5u splits are easily calculated using Eq. (91) as I (D 5u , T 2 ) = 1 24 [1 · 11 · 3 + 8 · (−1) · 0 + 3 · (−1) · (−1) showing that the splitting of D 5u contains T 2 twice. Using the trivial A 1 row instead, gives I (D 5u , A 1 ) = 0, so A 1 does not occur in the splitting. Checking the other irreducible representations verifies that D 5u is split into E, T 1 , and T 2 (twice). As a small check, it can be seen that the dimensions match: 2 + 3 + 2 · 3 = 11.