Constrained stability of conservative static equilibrium

Modelling of static structural stability problems is considered. Focus is set on problems where passive physical constraints affect the response to applied forces, and where more than one free parameter describes the setting. The existence of vibration frequencies at equilibrium states is investigated, as an indication of stability. The relevant Jacobian matrix is developed, with an emphasis on the necessity to formulate the constraint equations from an energy form in a conservative problem. The corresponding mass matrix is introduced, with zero mass contribution from constraint equations. Three different forms of the relevant Jacobians are considered, and alternative methods for the eigenvalue extraction given. Stability is discussed in a context of generalized equilibrium problems, where auxiliary parameters and equations can be included in a continuation setting. Examples show the formulation, implementation and interpretation of stability.


Introduction
Most structures react non-linearly to mechanical loading, even if significant non-linear effects may only appear at supra-operational loading situations. Non-linear response, and the related instability effects have often been seen as problems and failure mechanisms, but the phenomena are increasingly studied for morphing structures, in which the existence of multi-stable equilibrium solutions is the main advantage, and can lead to new engineering solutions. So does Reis discuss several situations where buckling is rather an asset than a problem [38]. Rafsanjani, Akbarzadeh and Pasini show how a meta-material with selected properties can be designed from geometrically non-linear effects in a non-homogeneous material [37]. Hamouche and co-workers derive closed-form expressions for bi-or tri-stability of thin shallow shells, and show how the shape of a structure can be controlled by active materials [25,26]. Emam and Inman give an extensive review of the morphing and energy harvesting potential in bi-stable composite laminates, where morphing is an "interesting feature of modern structures that enables them to change shape according to environmental or operational conditions" [14]. Arena  concept and design an adaptive inlet regulator, which snaps between open and closed configurations using morphing of a bi-stable component [3]. Similarly, Abdullah et al. show the bifurcation as the response of an activated bi-layer membrane, with functional implications [1]. Bertoldi discusses how architected cellular materials can be used for creating auxetic materials, for controlling the propagation of elastic waves, and for creating efficient energy-absorbers [7], and with co-authors discuss mechanical meta-materials in general terms [8]. Haghpanah et al. develop a designed material with very high energy dissipation potential through unit cells including elements with non-convex strain energy [24]. All the mentioned applications are related to sophisticated analyses of large configuration-changing displacements, instabilities and post-critical responses of structures.
Several types of non-linearities are inherently timedependent or at least depending on a sequential ordering of phenomena. They are then stated without inertia effects but with a fictitious time scale. Many important classes of non-linearities are, however, time-independent, and static equilibrium solutions can be sought as relevant representations. This is, for instance, the case for many geometrically non-linear settings, including conservative displacement-dependent forces.
Even if the mentioned problems are thereby by definition strictly static when a particular loading situation is considered, engineering viewpoints often lead to a treatment of 123 so-called quasi-static settings, where a variable load intensity is parameterized by a fictitious time, aimed to mimic an experimental situation [9]. These views contain no ingredient of real time, and solutions are just a set of independent static equilibrium states. The view on the response as continuous from an unstressed state up to an interesting load level, is deeply rooted in engineering simulations of mechanical response.
In the present work, it has been considered most appropriate to use the term sequence for the set of calculated individual equilibrium states, not to over-emphasize an assumed continuity between the states. This notation is thereby also relevant for evaluations of the variations in structural response, when the structure is parameterized in geometry or properties [12]. A particular setting is then when forces on a structure are multi-parametric, for instance, when the coupling between gas pressure and electrical activation can create unusual instability situations [34].
Stability is fundamentally dynamic in nature, but can under many practically relevant conditions be evaluated for static equilibrium states [6]. For a structure under acting forces, stability is a capacity to remain close to the static equilibrium situation after a minor disturbance. This view is related to a Liapunov stability condition, where a dynamic process starting near the equilibrium state stays close to it [35]. In practice, this means that any damping existing in the system will eventually bring the structure back to the static equilibrium. This view on stability is needed in order to give a basis for common stability investigations where otherwise "the precise notion of stability, always tacitly assumed essentially static in nature, is in fact, left undefined." [11, p. 5] The axiomatic definition of stability from the potential energy is discussed by Godoy, who notes that this has been in many cases an issue of "faith", but that the definition has been "of great value to improve our understanding of the buckling and postbuckling of structures" [21, p. 73].
The concept of stability thereby refers to one particular equilibrium state, and is a property of this state. Again, the typical engineering viewpoints will focus interest on how and when stability will be lost when following a parameterized load regime. Common methods thereby study critical states, by seeking the lowest force for which a deflected equilibrium exists for the structure. The view on stability of equilibrium thereby implicitly uses a simplified formulation of kinetic energy and inertia, analysing only the potential energy at equilibrium under fixed loads [43]. Although the reasoning is strictly not valid in continuous cases [11], theories for discrete systems lead to a conclusion that a sufficient condition for static stability is a minimum for the total potential energy under the considered forces. This typically demands a positive definiteness of the tangent stiffness, i.e., the second differential of the strain energy, as this matrix is established from all considered displacement components.
The present work focusses on the treatment of structures, where, in addition to prescribed loads, a set of physical constraints are defined. As opposed to active constraints, which are based on some external control strategy and addition of energy to the system [10,22,32], passive constraints are enforced by the considered system itself. Such systems are thereby conservative, and a total potential energy can be formulated, with the constraint considered.
The introduction of constraints inherently leads to the existence of auxiliary parameters in the static equilibrium setting. This implies that the positive definiteness of the tangent stiffness matrix is no longer conclusive, as the passive constraints introduce additions and modifications. The setting demands a clear definition of which parameters are fixed in the stability conclusion, and which are allowed to vary in, e.g., the fulfilment of constraints. This paper discusses in Sect. 2 the basic formulation for a multi-parametric discretized system, affected by some passive physical constraints. Section 3 discusses the stability of the established system, and the resulting eigenproblem, while Sect. 4 briefly discusses the used parameterizations, the solution method and the algorithmic implementation. Section 5 gives two numerical examples of the setting, whereafter Sect. 6 draws some conclusions from the presented work.

Mechanical modelling
Sophisticated simulation algorithms are needed for the evaluation and interpretation of non-linear static equilibrium solutions. The present work assumes that the structure and its forces are defined in discretized form, e.g., in a finite element context. The setting is thereby based on a displacement formulation, without any tuning modifications, e.g., for improving convergence. The formulation is assumed to be not overly sensitive to the numerical implementation, e.g., the scaling of the problem. Further, a static non-linear equilibrium problem is considered, where the system and the forces can be described through more than one free parameter, and the stability of an equilibrium state related to this context. A basic assumption is also that the problems considered are conservative, in the sense that they can be formulated through a total potential energy.

Basic problem setting
For the discretized conservative problems considered, the basic equilibrium problem is preferably stated using a total potential energy function where u consists of N u discrete displacement values, and consists of N Λ problem parameters, which can be structural parameters and parameters describing independent load effects. The splitting in the final part of Eq. (1) emphasizes that the potential energy consists of the internal strain energy W p , and a potential energy W f from external loads. From this, demanding stationarity of potential energy with respect to displacements u, the basic equilibrium relation is obtained as a set of N u residual equations, being the difference between internal and external forces corresponding to the terms in Eq. (1). A sub-index following a comma denotes a derivative, and a column vector form is used. In this form, N u can be a high number, whereas N Λ is typically low; the present work primarily focusses on cases with N Λ ≥ 2.
In the present setting, internal forces p can be affected by the parameters, whereas external forces f can be displacementdependent. The formulation in Eq. (2) thereby takes a more general form than the most common setting of the equilibrium problem, where just contains one load factor λ, p is independent of and f of u. In either situation, a large number of necessary structural parameters are implicitly considered as fixed, and are hard-coded in the simulation model.

Physical constraints
In many relevant problems, the static equilibrium solutions should fulfil additional physical constraints. These could, for instance, be a prescribed displacement relation, or a specified volume of liquid in a liquid-pressurized closed membrane, when pressure is used as the load parameter [46]. In examples below, physical constraints are also used to restrain rigid body movements with an objective method.
Keeping the energy-based formulation from Eq. (1), the basic energy expression can thereby be augmented by penalty terms considering a set of constraint functions g , accompanied by a corresponding set of Lagrange multipliers μ . The functions may be based on the displacements, but also a set of auxiliary parameters x . Even if the Lagrange multipliers are introduced as a technical tool to enforce the needed constraints, they are many times relevant results from the simulation. Also other types of constraints can, however, be valid for the equilibrium problem. These are expressed by additions to the total potential energy function of the form and dependent on a set of constraint parameters c and the auxiliary parameters x , but are independent of displacements, or they would be included in the term W f . Combining all constraining parameters, μ and c , into a set c leads to a splitting of into where c are constraint-enforcing, and x are auxiliary parameters. With the terminology introduced, this means that the total constrained potential energy is It is noted that these ΔW terms are considered, as they represent two typical modelling situations, even if more general expressions W c (u, ) could be relevant. Such expressions would not fundamentally modify the discussion below.
Differentiating the energy expression in Eq. (6) with respect to the parameters c , and demanding them to vanish, gives a set of augmenting constraint equations where functions g c depend on the discrete displacement components u, but in general also on the parameters . When Lagrange multipliers are used to enforce constraints, the equations set the penalty functions g = 0, but constraints of the form in Eq. (4) need be explicitly determined. Energy terms of the form in Eq. (4) will not give any contribution to the constrained residual. The added penalty term in Eq. (3) will, however, affect it, giving from Eq. (6) including additions to the residual in Eq. (2) of the form μ T g ,u . As the penalty terms will give force-like contributions from the Lagrange multipliers, the resulting constrained equilibrium equations are with internal forces independent of the constraint-enforcing parameters.
As the equilibrium equations combined with the physical constraints define the problem at hand, the setting demands solutions to an augmented equilibrium set [16] being N u + N c equations in N u + N Λ variables. The difference N x = N Λ − N c , which is the dimension of x , defines 123 the dimension of the solution space to the constrained equilibrium problem. Unless at least one parameter is included in x , the solution to the constrained static equilibrium problem is of dimension zero, and consists of one or more isolated single states.
When introducing N x > 0 parameters in x , and solving for one particular equilibrium state, N x selector equations are needed to make the state unique. Varying the parameters through these equations can be used to define subsets of equilibrium states fulfilling the equations in Eq. (10). In a solution algorithm, one of the functions demanded to vanish is a sequence selector, which is a generalization of common "arc-length" expressions. Other selector functions are further discussed below.
As described below, the N x ≥ 0 auxiliary parameters in Eq. (5) do not affect stability conclusions. Stability properties are, however, often evaluated as dependent on the parameters, which can be of diverse types, not necessarily related to a load case.

Differential response
A constrained static non-linear equilibrium solution demands fulfilment of Eq. (9), connecting the displacements u to the parameters . Variations to the residual functions are thereby (11) where parameters in are as in Eq. (5), and δ denotes a variation.
Similarly, the variations to the physical constraint functions in Eq. (7) are For the further considerations, the notation is introduced, which combines the state and the constraintenforcing parameters. For this form, the Jacobian matrix for the constrained equilibrium problem in Eq. (10) is obtained as again with an obvious matrix form of the second derivative.
The top left submatrix is the tangent stiffness matrix which, if energy is stated in the form of Eq. (6), may contain contributions from W p , W f and ΔW (if the penalty functions g are not linear in displacements), but not from ΔW c . It is noted that a particular constraint can be stated in several different forms g c , but as discussed in [46], physical constraints can -and should, for computational reasonsbe formulated to give a symmetric Jacobian. This also comes automatically when the problem formulation starts from a total energy expression, when the Jacobian is symmetric in the sense that with the negative sign coming from the definition of external forces in Eq. (8).
When the parameter space is as in Eq. (5), the total differential expression of the constrained equilibrium functions in Eq. (10) also identifies the differential expressions related to the auxiliary parameters x . In order to obtain a formulation similar to the form of common one-parameter equilibrium problems, a load-like (N u which in general consists of several columns and can contain diverse components. This leads to an incremental constrained equilibrium expression With N x = 1, this is a differential expression along the equilibrium sequence, i.e., a rate expression in fictitious time For any N x , this also defines the tangent space of the equilibrium manifold at a solution [z T , T x ] T , which is also the null space of the differential to the system in Eq. (10) through vectors fulfilling with t y organized similarly as a collection of all problem variables.

Stability
Stability is a fundamental issue in mechanical equilibrium problems, and is here restricted to the consideration of static stability. The stability property thereby refers to one static equilibrium state for a structure, and in essence relates to the capacity of the structure to remain close to an equilibrium state after a perturbation. The treatment below will discuss this stability property from two different starting points.

Stability without constraints
Common views demand for static stability a minimum total potential energy at the equilibrium state [20,29,42], which is a sufficient condition for a discrete setting. For an unconstrained problem with no g c , and thereby no c , this corresponds to positive definiteness of the tangent stiffness matrix, which is now evaluated at a particular equilibrium state, with fixed values for all parameters in = [ x ], and other problem parameters hard-coded. To be precise, this only applies to parameters which are present in the potential W c . Additional convenience parameters, e.g., Γ in Eq. (76) are not fixed, but do not affect stability.
The positive definiteness of K is thereby defined from all eigenvalues κ K i being positive in the one-matrix eigenproblem As an alternative view on the static stability criterion, the eigenproblem can be related to an un-damped free vibration equation, in the form of a two-matrix eigenproblem Given that all κ i are positive, the solution to this setting are modes i and frequencies √ κ i of natural vibrations, for 1 ≤ i ≤ N u , in small linearized vibrations around the equilibrium state. With any κ i ≤ 0, non-vibration responses exist. Stability is thereby judged by the sign spectrum of the eigenvalues κ i .

Comments and interpretation
The two seemingly very different views on the stability of a static equilibrium state above reflect the two basic criteria used. The positive definiteness of the potential energy and the tangent stiffness matrix, corresponding to the eigenvalue problem in Eq. (22), are related to a diagonalization of the current stiffness of the structure into its principal directions, as described by the eigenvectors K i . The eigenvalues thereby express the stiffnesses against movements along these directions, and positive eigenvalues show that a positive force is needed to displace the structure. As the eigenvectors span the incremental displacement space around the current equilibrium, no additional displacement can occur without an acting force. A zero eigenvalue also shows that the structure can be displaced, at least in some directions, without additions to the external force.
The vibration viewpoint, as expressed by the eigenvalue problem in Eq. (23), focusses on whether the structure will be able to sustain small linearized vibrations around the current equilibrium state, if an initial disturbing displacement is introduced, and the structure released from this state. The positive values for κ i , i.e., the real values for √ κ i , correspond to frequencies in a vibration-type response, which will continue without any external action. With a zero eigenvalue, introduction of a small incremental displacement in the corresponding eigenvector direction will not lead to vibration, but to another equilibrium state. It is emphasized that the vibration problem is evaluated at an equilibrium state, with internal forces affecting the tangent stiffness matrix. It is noted that either view is implicitly based on the observation that only the sign spectrum of the eigenvalues of the tangent stiffness matrix is of importance for the stability conclusion, not the precise values.
The criteria for stability are evaluated for one particular equilibrium state, where displacements u are corresponding to some external forces as expressed by a fixed set of parameters in x and hard-coded parameters in the model. This means that not just load parameters are fixed in the stability evaluation, when one equilibrium state is considered. This aspect is of fundamental importance for the interpretation of stability, when a multi-parametric setting is introduced.
For the common engineering setting, where a fixed structural model is considered when affected by some load case described by a single load parameter x = [λ], the nonlinear response evaluation implicitly sees the displacements as functions of this single parameter. The stability conclusions for the evaluated equilibrium states will thereby also be parameterized. As most simulations of this type start from an unloaded stable initial state, and focusses on finding the first critical equilibrium state for the model, the view is typically that the used load parameter is the one causing the instability, and expressions of the form that stability given with respect to this load are common. For the particular setting, this is also a fully valid view, as the critical stability implies that incremental displacements can occur without increments to this load parameter.
In the more general setting used here, with several parameters in x and possibly parameters of different kinds, the view on an instability-causing load is no longer valid. The basic definition of stability for one particular equilibrium state -with all parameters fixed-must be used, disregarding any parameterization which might have been used when finding this state. The notion of a particular parameter caus-ing the instability is thereby in general not relevant. Along a parameterized sequence of equilibrium states, the stability conclusions relate to certain parametric regions, which are delimited by critical states, where the stability properties change. For instance, Fig. 1c related to a numerical example below shows how a Mises truss will show different stability properties under a fixed vertical force, when the height of the truss is parameterized.

Stability with constraints
The stability of a constrained static equilibrium solution can be similarly evaluated, based on the Jacobian matrix in Eq. (14). As the constraint equations in Eq. (7) have no attached masses, the relevant undamped free vibration equation is written as where according to Eqs. (14), (21) and (16) and the relevant diagonal mass matrix is described in schematic form as with a zero matrix of indicated size. Given that allκ i > 0, the constrained natural vibration modes˜ i thereby describe vibration modes expressed by a vector corresponding to z in Eq. (13).
The demand for stability of the constrained equilibrium state is again that allκ i > 0. A critical state, representing a transition in stability is characterized byκ i = 0. At a critical state with one or moreκ i = 0, Eq. (24) shows that linearized solutions exist whereKδz = 0, with δz any linear combination of the critical eigenvectors. These are non-trivial solutions to Eq. (18) with x fixed.

Comments and interpretation
It is noted in a comparison to the unconstrained stability criteria in the previous section, that no expression corresponding to the one-matrix eigenproblem in Eq. (22) is obvious for the constrained case, due to the massless constraint equations. Only the vibration-type two matrix eigenproblem will thus be considered.
The displacement components u thereby vary together with the constraint-enforcing parameters c needed to maintain the physical constraint during the vibration cycle, with the same angular frequency √κ i .
Major interest in a mechanical investigation of a structure is focussed on the critical states, where the stability properties change under parametric variations. With more than one component in x , this is a generalization of the common formulation without constraints and only one force parameter. With an arbitrary number N c of constraints and N x = 1 auxiliary parameters in the problem formulation, similar conclusions regarding limit and bifurcation states can be drawn as in a one-parametric unconstrained setting. For N x > 1 a more general interpretation must be introduced.

Constrained free vibrations
The stability investigation for a constrained conservative equilibrium problem focusses on the eigenvalues of the problem in Eq. (24), based on matricesK andM. For simplified notation, the treatment below will see the Jacobian matrix as with c = g c,u and b = g c, c . As the evaluation is performed at a specific equilibrium state, the matrices are constant. It is in the treatment below assumed that the matrix c has, or has been reduced to have, linearly independent rows, but also that the matrix b is symmetric, which comes automatically from the energy form. The treatment of the constrained problem is dependent on properties of the block b: if it is zero, singular or has full rank, where zeroes typically come from the penalty terms in Eq. (3). The treatment will be divided into three cases below. Utilizing that only the signs of eigenvaluesκ i are interesting for the stability conclusions, and that this sign spectrum is not affected by the contents of the mass matrix M, as long as it is positive definite, the treatment below will also set the mass matrix in Eqs. (23) and (26) to M = I N u . This corresponds to the implicit assumption when the static stability criterion is a demand for a positive definite tangent stiffness matrix K.
With the present objectives, it is also assumed that the number of constraints N c is rather small, say in the order of ten, even if the number of degrees of freedom N u is several thousands.

Treatment of eigensolutions
In a systematic investigation of static structural stability, the evaluation of eigensolutions related to the relevant Jacobian matrix is a main objective. Here, the eigenvectors are of interest, as they can give a basis for the initiation of secondary equilibrium sequences, but a main focus is set on the eigenvalue sign distribution, i.e., the numbers of negative, (approximately) zero, and positive eigenvalues. With regard to precise values, the evaluation of eigenvalues close to zero are of main interest, while the clearly positive and negative ones are generally of lower importance. This agrees with the notion of eigenvalues being the principal stiffnesses of the considered system.
The efficiency of numerical methods for extraction of matrix eigensolutions is highly dependent on the properties of involved matrices [5,45]. Methods are developed for particular settings, where symmetry, sparseness and positive definiteness of the involved matrices are main aspects, as is an overall knowledge about the expected eigenvalue spectrum. Possibilities to handle a one-matrix setting may also in many cases simplify treatment, as do options to define a required accuracy in the eigenvalues.
The present setting, as described by Eq. (24), is characterized by symmetric matrices, normally sparse, a semi-definite diagonal mass matrix, and often of clusters of identical or very close eigenvalues, sometimes also situations where several eigenvalues are simultaneously zero due to symmetries in the model. Normally, only a limited number of eigenvalues are requested, sometimes only their signs, sometimes their exact values and sometimes also the eigenvectors, depending on the phase of the algorithm.
Even if it is conceivable to develop algorithms for handling of this exact setting, it has in the present work been chosen to set the problem in forms which can be handled by builtin functions in Matlab, 1 considered as industry standard for basic linear algebra tools. The basic Eq. (24) with matrices from Eqs. (25) and (26) is re-written with the objective to allow reliable and efficient handling of the separate cases by the available functions.
When evaluating eigenvalues without too high demands on precision, a strategy based on a Sturm sequence [41] can preferably be used, by an LDL-factorization of the Jacobian matrix shifted by a multiplier of the relevant mass matrix. This method is in the present implementation used to decide the numbers of negative eigenvalues below a certain tolerance − γ and of zero ones within 0 ± γ , but also in non-critical situations to give a coarse estimate to the eigenvalue closest to zero, which can be used for predictions of approaching critical states. A systematic usage of the same method in a bisection strategy can give also the eigenvalues, but is inefficient for an accurate determination.
gives a two-matrix eigenproblem system according to Due to the semi-definiteness of the matrixM, the system will demand eigenvectors with displacement parts fulfilling the orthogonality where the second term vanishes due to Eq. (30). Thereby, ψ i does not affect the equation, and the eigenvalue could be computed from if the eigenvector were known. This Rayleigh quotient expressions shows that the eigenvalues of the constrained problem are limited by the minimum and maximum eigenvalues of K [41,45], with consequences in constrained dynamics [4]. It is, however, not obvious how to practically find the eigenvectors φ i under the constraint in Eq. (30). If these could be obtained, the components of ψ i would be the constraint-enforcing parameters needed to keep the vibrations orthogonal to the constraints, and could be solved from the over-determined Eq. (29) 1 as The full eigenvectors, corresponding to eigenvaluesκ i would then be obtained by suitable normalization of˜ i in Eq. (28).
As another approach, the orthogonality condition can be introduced by a projection matrix which makes any N u -dimensional vector orthogonal to the columns of c T . As this matrix has N u − N c unit eigenvalues and N c zero ones, the two-matrix eigenvalue problem will give N u eigensolutions (φ i ,κ i ) for the constrained system, and, after a re-projection of the obtained eigenvectors, The obtained eigenvectors will include N c spurious ones, characterized by (almost) zero components in φ i after the transformation in Eq. (36). They can therefore be easily removed, leaving just the relevant eigensolutions, which are then treated by Eqs. (33) and (28).
Although formally correct and useful for small problems, the method becomes inefficient for practical problems, due to the non-sparsity of the matrix and consequently the matrices in Eq. (35). As the mass matrix can also not be assumed to be diagonal after the transformation, a two-matrix eigenproblem need be solved.
Two more realistic methods for solving the constrained eigenproblem are based on the discussion above. In addition, it is also conceivable to develop methods similar to power or subspace iterations which orthogonalize the iterative vectors with respect to the constraints.

Fictitious small masses related to constraints
A numerically feasible method for solving the eigenvalue problem described by Eqs. (24)-(27) is based on an introduction of small fictitious masses for the constraints. Replacing the mass matrix in Eq. (26) by a matrix of the form with ε a small value, the matrixM ε thereby is the square of a diagonal matrix This leads to a two-matrix eigenproblem of size (N u + N c ), which can be solved with standard methods as the matrices are symmetric and the mass matrix now positive-definite. The system can also be re-formulated through a transformation of Eq. (24) into a one-matrix problem of the same size, and the form which in essence has multipliedK andM ε by m −1 ε from both sides. Here, the lower-right block of the Jacobian matrix is zero, since b = 0 here. This eigenproblem can also be solved with standard methods, utilizing the sparsity of the matrix. The relevant eigenvectors from either of Eqs. (39) or (40) are transformed into the final result by˜ and a normalization. Either of the two above formulations will give a set of 2 N c spurious eigenvalues of large magnitude. As these correspond to eigenvectors with low values in all displacement components, the irrelevant eigensolutions are easily discarded. This removal of some obtained eigenvectors will reduce the complete set of (N u + N c ) eigenvectors to a number of (N u − N c ). If only a subset of eigensolutions is needed, the forms above allow this, if suitable basic eigensolution algorithms are utilized.
Reduced basis solution A potentially efficient method for solving a constrained eigenvalue problem of the considered type is based on a projection of the eigenvalue problem in Eq. (24). Similar to the projection matrix in Eq. (34), an N u − by − (N u − N c ) matrix P of linearly independent basis vectors fulfilling cP = 0 is then created. The eigenvalue problem for the displacement components, is then a two-matrix eigenvalue problem of size (N u − N c ), which implicitly removes the irrelevant eigensolutions.
Although the eigenvaluesκ i are immediately obtained, the corresponding eigenvectors are re-constructed in two steps, where the first back-projects the displacement components, according to and the second is identical to Eq. (33), giving ψ i and leading to the full eigensolution (κ i ,˜ i ) to Eq. (28). The efficiency of the presented method is completely dependent on the possibilities to find a sparse matrix P with columns spanning the orthogonal complement to the constraint matrix c T . As further discussed below, the creation of the matrix is based on a trade-off between the sparsity of the resulting matrices in Eq. (42) and the well-conditioning of the basis vectors.

The case b = 0, with non-singular b
The situation when the constraint equations explicitly involve the constraint-enforcing parameters is fundamentally different from the case above, as it is not just a restriction on the displacements. The treatment of this situation is dependent on the rank of the matrix b.
When the matrix b is non-singular, there exist N u eigensolutions to Eq. (24), in the (N u + N c )-dimensional space considered for z. Then, a formal treatment performs a partinversion of the Jacobian matrix in Eqs. (24)- (26), giving first a one-matrix eigenproblem of size N u for the displacement components This is followed by the introduction of the constraint part The whole eigenvector for the problem in Eq. (24) is then obtained as which can be suitably normalized. It is noted that no elimination of irrelevant eigensolutions is needed here. As the part-inversion in Eq. (44) will lead to a filled although symmetric matrix, the method is only useful for small problems.
The methods from above, with small fictitious masses ε 2 associated to each of the constraints, leading to Eqs. (39) or (40), but now with a non-zero b matrix, can also be used. The obtained N c eigenvectors to the full system with only small displacement components are easily discarded.
Also, similar to the procedure described by Eq. (42), a linearly independent basis matrixP of size (N u + N c )-by-N u can be established. This is orthogonal to [c, b] T , and reduces the Jacobian matrixK to size N u . This gives the two-matrix eigenproblem which implicitly removes the irrelevant eigensolutions. The final constrained eigenvectors, corresponding to the calculated eigenvalues are then obtained as

The case b = 0, with singular b
The situation when b is non-zero but singular will need special treatment, the reason seen from the formal expressions in Eqs. (44), (45). This case is the most general one, and the methods described are useful for also the special cases discussed above, but then somewhat less efficient. For the general setting, treatment of the full problem in Eqs. (24)-(26) is possible, but this becomes demanding if the algorithm used does not allow extraction of a limited number of eigensolutions to a two-matrix setting with a semi-definite mass matrix.
The methods above with fictitious small masses attached to the constraints, i.e, settings according to either of Eqs. (39) or (40) are still possible. When followed by a transformation according to Eq. (41), and a subsequent elimination of eigenvectors with very small displacement components, the eigensolutions are obtained. The number of relevant eigensolutions will in this case be N u − N s , where N r is the rank and N s = N c − N r the rank deficiency of the matrix bwith N s = 0 for regular b, and N s = N c for b = 0.
As a general method for treatment of this case, a singular value decomposition of the symmetric small matrix b can be used to obtain where the pseudo-singular values σ * are ordered in descending order of magnitude in the diagonal of matrix σ * . The slight modification is needed as the matrix b is not necessarily positive definite, and allowed as only the rank of the matrix is sought. Based on Eq. (49), the orthogonal matrix v is then used in an orthogonal transformation matrix which, when operating on the full Jacobian matrixK in Eq. (27), gives a stiffness-like matrix together with the singular constraints, according tõ associating a rank-containing part of the constraint functions to the tangent stiffness matrix, and pushing a set c of N s rows to the end of the set. As the mass matrixM is unaffected by this transformation, the transformed two-matrix eigenproblem is where the form is now agreeing with the form discussed in Sect. 3.4.1 above, as the lower-right block is zero. Analogously to Eq. (42), a transformation matrixP can be found, but now spanning the orthogonal complement to (c ) T . With the matrixP of size (N u + N r )-by- (N u + 2N r − N c ), the problem is defined as whereM appends N r rows and columns of zeroes toM. The eigenvectors are then back-transformed bŷ A procedure similar to the one in Eq. (33), but now as calculates the remaining N s components inˆ i = [φ T i ,ψ T i ] T , which, after one further transformation, and a subsequent normalization, gives the final eigenvectors to the constrained system. It is, however, noted that the eigenvalues were obtained already by Eqs. (52) or (53).

Creating an orthogonal linearly independent basis
The success of the methods involving a reduced basis described by Eqs. (42), (47) or (53) is completely dependent on the possibility to efficiently create a simple transformation matrix, which does not introduce extensive fill-in in the operating matrices. With Eq. (42) as prototype, the task is to create a matrix P which spans the vector space orthogonal to the columns of matrix c T in Eq. (27), i.e., fulfills cP = 0, where c is N c -by-N u and P should be N u -by- (N u − N c ). It is noted that the columns of P need not be orthogonal, but just not too close to parallel. A practical implementation of this method starts with a matrix P (0) = I N u , and evaluates the matrix product cP (0) . The matrix P (0) is then successively updated to P (i) , (i = 1, . . . , N c ) by orthogonalizing its columns to the columns of c T , one by one.
The most straight-forward method for going from matrix P (i−1) to P (i) , when all columns of P (i−1) are orthogonal to columns 1, . . . , i − 1 of c T , identifies all non-zero components in row i of the product cP (i−1) , i.e., indices (i) . Based on the coefficients where the notation emphasizes that row i of the matrix product is considered-the indexed columns of P (i−1) are updated according to where the colon index notation refers to the entire column. The non-indexed columns are just transferred to the new matrix (keeping the order of columns), as they already fulfil orthogonality to column i of c T . The updating is performed for all but the final index in (i) ; the last of the columns is discarded.
Due to the transformations in Eqs. (57), (58), the whole row i in cP (i) is vanishing. After handling all the N c rows of c, the matrix P = P (N c ) is the sought matrix, with remaining columns orthogonal to the matrix c T .
The method described above leads to a very low fill-in, and is easily described in efficient algorithmic form. When the number of non-zero components in the matrix c is high, the matrix P, however, easily becomes ill-conditioned in the sense that its columns become close to parallel. In these cases, a slightly more elaborate procedure is preferred. In this form, the index vector (i) is sorted such that it gives the non-zero components of (cP ) in ascending order of magnitude.
This re-ordering will in Eq. (57) lead to all |α (i) k | ≤ 1, which improves the orthogonality of the columns in P (i) , but at the cost of -sometimes significantly-increased fill-in. Also the refined algorithm is easily algorithmically described, with only the c matrix as input.
Practical experiences from experimentation with the larger examples below, with well-filled c matrices, indicate that the refined procedure is necessary for obtaining eigensolutions of sufficient precision. For other problems, with more point-wise constraints, the basic method has been deemed sufficiently accurate, and also overall more efficient.

Parameterized equilibrium sequences
The above discussion focusses on the stability properties of a constrained equilibrium state defined by z and with a fixed set of auxiliary parameters x . This section discusses the introduction of such parameters in a multi-parametric structural model, and the solution of parameterized equilibrium sequences.

Auxiliary parameters
When analyzing a multi-parametric equilibrium problem, parameters might be hard-coded in the problem formulation or variables in the simulation, i.e., introduced by the set x of size N x . This is a generalization of the engineering approach, which evaluates the response of a structure, and in particular its stability, as dependent on a one-parameterization of a load case.
The introduced parameters can be used to describe any kind of properties of the structural model or the load cases considered. Introduction of N x auxiliary parameters x will give an N x -dimensional solution space, which can be reduced by a set of auxiliary equations where an (N x − 1)-dimensional set of conditions leads to response curves. The functions give conditions for relations between variables and parameters, without introducing any constraints on the solution, and are often simple. Since they are for each studied equilibrium state constant, they never affect stability investigations. The following will discuss different auxiliary equations, and their relation to introduced parameters.

Sequence parameterization
A common engineering setting of static equilibrium analyses is to parameterize a main case by a load parameter , with all other parameters hard-coded. Even if other approaches have been suggested, e.g., based on asymptotic methods [30], the most common approach is to solve displacements u j for a set of λ j , even if technically so-called displacement or arc-length control [13,19,31,40] is often used. Such strategies parameterize an equilibrium sequence in the combined space (u, λ) [15], and can be described as one selector function g x appended to the residual expressions. Stability is judged based on the tangent stiffness matrix at each state, and thereby related to the load level λ.
Similarly, with N c constraint equations in Eq. (10), z according to Eq. (13), and an (N c +1)-dimensional parameter set = [ T c , λ] T , the first N c parameters are constraintenforcing, whereas the parameter λ allows a parameterization of the equilibrium sequence by any selector function involving the variables [z T , λ] T [15]. In a sequence-following method, where each step finds one new static equilibrium [z T i+1 , λ i+1 ] T for (i = 0, . . .), the selector function is updated step-wise, demanding a vanishing with y according to Eq. (20), and introducing a fictitious time increment Δτ i for each new solution state. The superscript * shows that a function of this form is always introduced automatically as part of the sequence-following algorithm.

Parameterized structural model
Geometric or material parameters for the considered structural model can be introduced by x . A sequence of static equilibrium solutions for a hard-coded load case is then obtained as function of the parameter; examples are shown below. Stability can then be evaluated for this load as dependent on the parameter.

Convenience transformation
With several parameters in x , any kind of auxiliary equations can be used to select specific parameterized equilibrium states [16,23]. Addition of further parameters in x and related functions in g x are then used for transformations between parameters in the simulation. This idea can be used for, e.g., summation of some displacement components, intended for post-processing of results, or for convenience in interpretation. Below is shown how gas pressure and gas amount in a closed membrane can be connected by a condition on the parameters, which allows different views on the loading.

Criticality function
Critical equilibrium states when increasing a load intensity, or when traversing a generalized equilibrium sequence, are commonly characterized by zero eigenvalues of the tangent stiffness matrix, as in Eq. (23), even if other similar definitions are sometimes used. When using physical constraint equations as in Eq. (10), the criticality of the equilibrium state can be similarly judged based on the vanishing of a function the lowest magnitude eigenvalue of the Jacobian matrix in Eqs. (24), (25). This is used in finding one critical equilibrium state for a hard-coded structure and one-parameter loading, replacing the arc-length function. It can also be used in a two-parameter setting to obtain a sequence of critical equilibrium states for a variable parameter. With more parameters and selector functions, this is thereby a more systematic approach, but similar in idea, to the specialized one used by Li and Healey in order to find stability boundaries in parameter space [33].

Critical sequence conditions
As further presented and discussed elsewhere, following of critical equilibrium sequences is not always a robust procedure if based on Eq. (61), and other methods are needed. As one alternative, auxiliary conditions can keep the solutions on a specific sequence. This is relevant for bifurcation states, where the model symmetry is broken on the secondary sequences. When the sought critical equilibrium states are situated on the primary sequence, conditions of the form are introduced, expressing the desired symmetry by a set of functions, and accompanied by forces f s . The inclusion of a set of such auxiliary conditions will decrease the condition number of the iteration matrix, but will not affect the equilibrium state, as the forces f s will converge to zero. The difference between these symmetry-enforcing conditions and the constraints introduced by Eq. (3) is obvious. It is also clear that the conditions in Eq. (62) must be chosen in relation to the specific symmetry conditions of a particular critical situation.
As another approach, in order o avoid iterating at the exact critical state, the critical selector condition in Eq. (61) can be replaced by the scalar functioñ where κ c are the N 0 lowest magnitude eigenvalues, when a sequence of N 0 -fold criticality is followed. The value κ ε is a numerical tolerance for what is considered as a vanishing eigenvalue. A sequence followed with this selector will thus trace a border line for the zone of critical equilibrium states. With a suitable value for κ ε , critical states very close to the exact ones will be solved, known also to be on the more stable side of the critical border.

Implementation
The above setting is implemented in a further development of the basic sequence-following algorithm [16]. The main difference is the systematic distinction among the added equations, as representing either physical constraints or selector conditions in the multi-parametric solution space. Following [46], the non-linear set of equations to solve for one equilibrium state is stated as According to Eqs. (20) and (13), the variable vector y contains the displacements u, the constraint-enforcing parameters c and the auxiliary parameters x . The expressions are the residual forces r c , the physical constraints g c and the selector functions g x . The problem specification thereby must provide N u +N Λ −1 function values from N u +N Λ variables in a problem-specific code segment. The major part then comes from general element functions and from the penalty functions, but the constraint functions in g c , and all but the final function in g x need be defined for a particular case.
Similarly, the corresponding differential expressions must be provided for a specific problem setting, according to where r c,u ≡ K is mainly formulated from general element expressions. The algorithm traverses a sequence of equilibrium states with a continuation algorithm [2,16,39], assuming an initial equilibrium solution [z T 0 , λ 0 ] T is known. When solving for a new equilibrium solution y i+1 = y i + Δy i , from a previous solution, the prediction to the increment in the new step is cf. Eq. (19) [16]. The step length Δs i is chosen before the prediction step. This choice affects the sequence selector function in subsequent iterations in the increment.
In an iterative formulation for Newton corrections j = 0, . . ., the differential G ,y is also evaluated at y j i+1 = y i +Δy j i . Iterations according to Eq. (67) are performed to a low tolerance, i.e., until whereafter solution y i+1 is stored, and its properties evaluated.
The algorithm isolates and identifies special equilibrium states along the traversed sequence. The special states include solutions where the number of negative eigenvalues of the relevant Jacobian matrix-K in Eq. (21), orK in Eq. (25)-changes, but also solutions where any of the problem parameters or their tangent components vanish. The implementation utilizes the continuity of the current problem class in an accurate isolation of special states on the sequence.

A demonstration problem
The concepts above are first demonstrated for a Mises truss, Fig. 1a. Neglecting the possibility for local buckling in this demonstration example, which would demand a choice of finite deflection beam formulation, the problem was modelled as a truss, by analytical treatment of the non-linear Green-Lagrange strains and linearly elastic second Piola-Kirchhoff stresses with the unstressed state as reference. For L = 1 m, a scaled axial stiffness E A = 1 N, left two displacement components and the corresponding explicitly applied forces (u i , f i ) (i = 1, 2) and a geometrical parameter h as basic definition of the problem. The basic potential energy for the system was Equilibrium sequences were obtained for a single parameter in x together with a sequence selector function g * (y) automatically inserted. Results for hard-coded h = 1.5 m, f 1 = 0, and parameterized x = [ f 2 ] are shown in Fig. 1b, and for hard-coded f 1 = 0, f 2 = − 0.15 N with Fig. 1c. For both settings, primary solutions with u 1 = 0 were obtained without any conditions or constraints. Dashed lines in subfigures (b) and (c) emphasize that the figures are different sections through the same manifold. Isolated critical states are marked, and the numbers of negative eigenvalues for the tangent stiffness matrix K are shown.
As a perturbed, but not constrained, problem, two parameters x = [ f 1 , f 2 ] T were introduced, with a hard-coded h = 2.5 m. Introducing the selector functions equilibrium states with displacement u 1 = X = 0.2 m were selected, cf. Fig. 2a.
In general, f 1 = 0 for this sequence. Special solutions were found for f 2 = 0, for f 1 = 0, and from vanishing tangent components. A singular tangent stiffness matrix was found at four equilibrium states, of which two were close to, but not at, the states with extremum f 2 , cf. Fig. 2b. The special states were found by bisection to high accuracy without any particular formulation. Similar situations were found when the critical state appeared close to states with f 1 = 0.
In order to find the parametric variation of critical states with height h, three parameters x = [ f 1 , f 2 , h] T , were introduced. Demanding criticality together with a maintained symmetry, i.e., a condition u 1 = 0, three selector functions were introduced, containing the lowest magnitude eigenvalue of the tangent stiffness matrix. Starting from isolated critical states reported in Fig. 1b, critical sequences in Fig. 2c were obtained, verifying that a bifurcation is reached before the first limit state if h ≥ √ 3. Introducing a physical constraint on u 1 , the problem was formulated through a Lagrange multiplier μ 1 as with f 1 = 0 hard-coded. The problem was solved for X = 0.2 m, introducing c = [μ 1 ], g c (u) = X − u 1 , and a corresponding load-like vector h = − g T c,u = [1, 0] T , giving a symmetric Jacobian. With x = [ f 2 ], and one sequence selector function g * (y), the total external forces on the structure are f( ) = [μ 1 , f 2 ] T . The solution sequence obtained was identical to the one in Fig. 2a, but only two critical states were found, coinciding with extremum f 2 .
For a displacement constraint of the form u 2 = Xu 1 , with X a constant and f 1 = 0 hard-coded, the potential energy form was written , and the sequence selector function g * (y), solution sequences were obtained for hardcoded X . The total force acting on the structure was f( ) = It is noted that for neither of the two constrained settings, the parameter in x affected the tangent stiffness matrix, as the penalty functions g were linear in the displacement components. With, e.g., a circular constraint function g = (u 2 1 + u 2 2 − ( * ) 2 ), an addition to the basic tangent stiffness K = W ,uu would appear.

A large-scale example
The inflation by gas of a spherical membrane of unstressed radius R = 50 mm, and thickness t = 0.05 mm was studied. The material was assumed as an incompressible isotropic Mooney-Rivlin model, with the constants related to the two first strain invariants c 1 = 0.25 MPa, c 2 = k c 1 .
The membrane was pressurized by an internal overpressure p, in relation to an external, ambient, pressure p 0 . The included amount of gas was an alternative parameter, and measured by Γ = ( p + p 0 )V , with V ≡ V (u) the enclosed volume at a state u. The instabilities existing for this problem are extensively studied in literature [28,36,44].
It is noted that with present data the over-pressure is significantly lower than ambient pressure, aiming at the representation of balloon-like membranes on earth. A study of the dependence on ambient pressure is also discussed below.

General settings
The basic total potential energy for this problem was with the strain energy W p dependent on geometric and material parameters hard-coded in the expressions. Rigid body motions were constrained by an additional energy term with μ R containing six Lagrange multipliers. The expression considers displacements u k and initial positions X k of all nodes in the model. In all converged solutions μ R ≈ 0 to high precision. The needed constraining functions g c (u), below denoted as g R , and their corresponding contributions to the residual equations were obtained from Eq. (75).
Recognizing that the symmetry properties of a discretized mesh can significantly affect simulation results [27,47], different computational meshes were used. All were based on flat, triangular membrane elements with assumed local plane stress conditions [17]. Main results are given for a 5120 element mesh based on a regular refinement of an icosahedron mesh, and fulfilling I h symmetry in Schoenflies notation. 2 The refinement means that the final mesh will use nodes of 30 different positions with respect to symmetry. In particular, the nodes at the axis crossings with the sphere are different, and give slightly different radial expansions. Results are here given for a node initially at Z = R, being a vertex of the icosahedron.
Two other models were based on 4096 elements: a regular subdivision of a tetrahedron, fulfilling T d symmetry, and based on an 8-sector and one mirror mesh, fulfilling D 8h symmetry. Low resolution element meshes of these symmetry classes are shown in Fig. 3.
In all simulations reported below, a fictitious mass setting according to Eq. (39) was used for eigen-analysis, with ε = 10 −5 .

Response to over-pressure forces
For an otherwise un-constrained setting, parameters were introduced as c = [μ R ] and x = [p, Γ ] T . Selector functions were introduced as defining the gas amount Γ as an auxiliary parameter, from a hard-coded p 0 . The Jacobian of the constrained system in Eq. (14) was symmetric from Eq. (75), as f , c = − g I c,u , while g c, c was zero. The matrix h contained a first column − V ,u and one zero column, while Γ gave a contribution to the differential g x, x .
With the I h model, and ambient pressure p 0 = 100 kPa, the relations between radial expansion and over-pressure are given in Fig. 4a, for a set of k. Local maximum pressure states were isolated for low k, with accompanying local minima if k > 0. Figure 4b shows the relation between over-pressure p and enclosed amount of gas Γ . Also the gas amount can show a limit state under inflation for low ambient pressure and k < 0, Fig. 4c. This limit amount situation was here noted only from a change of sign for the Γ tangent component.
Critical states, for fixed p, are marked in the figures. These are primarily limit states, but bifurcation states, of different multiplicities, were found for k < 0. Compared to the analytical solution, the mesh introduced an imperfection. This is seen for the case k = − 0.05, which gave a limit state as well as a turning state in the considered radial expansion measure, and returned-in the used projection-parallel to itself, as can be vaguely seen from Fig. 4c.

Symmetry of model
An analytical treatment of the gas-pressurized sphere (presented elsewhere) shows that bifurcations will occur from the primary equilibrium sequence for certain hyper-elastic These bifurcation states will be more or less well represented by a chosen mesh. Multiple bifurcations are thereby split into several lower order bifurcations, or converted into limit states by mesh un-symmetries acting as imperfections [18]. Such limit states indicate that the primary equilibrium sequence for the discretized case deviates from the one valid for the continuous case, however, without a bifurcation. Figure 5 shows the computed response curves for a case with k = − 0.05 and p 0 = 100 kPa, and the three different meshes. The figure shows at which bifurcation state for the continuous case the respective models are no longer able to yield an essentially spherical solution. Figure 5a shows the pressure-expansion relation, which is magnified in Fig. 5b, while Fig. 5c shows the deviation from spherical expansion, evaluated from displaced nodal coordinates. The figures show, in addition to sequences started from zero over-pressure, a sequence for the I h model restarted for larger expansions on the primary continuous sequence. The observation is that the continuous case will be possible to represent with a set of unconnected equilibrium sequences.

Gas amount as main load parameter
Simulations were performed with gas amount as load parameter, based on a potential function As the constraint term is of the ΔW c form in Eq. (6), the corresponding seven constraint functions were , and p 0 a hard-coded ambient pressure. It is noted that the second expression in Eq. (78) is the same as the first expression in Eq. (76), but reformulated to give a symmetric Jacobian matrix. The moving of the function from a condition in g x to a constraint in g c is also of significance. With over-pressure p now included in c , the Jacobian included W c, pp = − Γ /(p 0 + p) 2 = 0 as a diagonal term of the b submatrix in the schematic Eq. (27), and a zero row was appended to the c matrix. The effective load vector h only contained W c, pΓ = 1/( p 0 + p). With N x = 1, a sequence selector function g * (y) was automatically added by the algorithm.
The equilibrium sequences [u T , p, Γ ] T , were the same as above, but stability conclusions differed, as the present case studied stability with fixed Γ , while over-pressure p was allowed to vary. Figure 6a shows that the present case did not show a critical state at the maximum over-pressure.
The magnified Fig. 6b shows that the maximum amount state was now noted by a singular Jacobian. The same bifurcation states along the sequence were isolated for both forms. Figure 6b

Response to ambient pressure
Simulations were performed for a spherical membrane, which was inflated by an amount of gas, closed, and then subjected to variable ambient pressure p 0 . The potential of this setting was the same as in Eq. (77), but now with Γ hard-coded, c = [μ T R , p] T , and x = [p 0 ], with constraint functions as in Eq. (78). The setting also leads to the same g c,u , f , c and g c, c as above, but an effective load vector Sequence-following added a single selector function g * (y), and did not need any symmetry condition.
A situation similar to those in Fig. 4a-b was first studied for p 0 = 100 kPa and k = 0.18, giving three equilibrium states for p = 0.85 kPa: stable at u R = 22.76 mm, Γ = 0.1624 kNm; unstable at u R = 45.95 mm, Γ = 0.3723 kNm; stable at u R = 79.10 mm, Γ = 0.9069 kNm, with stability judged for fixed p.
With these amounts successively introduced as three cases of Γ , results over a range of ambient pressures p 0 are given in Fig 7. Stability conclusions All equilibrium states in Fig. 7 were found stable based on the constrained Jacobian, even if the tangent stiffness matrix K had one negative eigenvalue in parts of the sequences. This illustrates an important distinction in stability conclusions, dependent on setting.
The constrained setting considers Eqs. (77) and (78) and studies stability with p 0 and Γ fixed, which demands p to vary with displacements in order to fulfil the Γ definition. The investigation of the tangent stiffness matrix is based on an analysis of the same problem, but based on Eq. (74) and (a re-formulation of) Eq. (76) with x = [p, p 0 ] T . This would yield the same equilibrium sequence but give stability conclusions for fixed p, p 0 and Γ .
As stability is in the latter case only related to expansion and over-pressure, Fig. 7a also clearly shows similar regions of instability in all three cases considered, cf. Fig. 4a. Figure 7c also shows how the three starting states are changing due to variable ambient pressure, with the amount of gas kept constant; the marks connect in each case the solutions for p 0 of 50, 100 and 200 kPa, but it is obvious that the curve of p versus u R is followed by continuous changes to p 0 . The three cases are shown to react differently to changes in ambient pressure from the initial 100 kPa, as Case 2 will be always unstable in this ambient pressure range, Case 3 always stable and Case 1 changing stability within the range.

Parameterized deviations from sphericity
Deviations from a spherical structure was modelled, by introducing a non-uniform scaling of the unstressed I h unit sphere model, with radii (R x , R y , R z ). Simulations were performed for variable R z , while R x = R y = 50 mm, with k = − 0.05, p 0 = 100 kPa, and other data as above. With successive hard-coded values for R z , parameters were introduced as c = [μ R , p], x = [Γ ]. With a potential energy from Eq. (77), constraint functions according to Eq. (78) were introduced, and a sequence selector function g * (y) used.
For ellipsoid geometry, the I h symmetry was reduced to a D 5d symmetry, when the sphere was stretched along an axis through two vertices of the basic icosahedron. Some results for different R z are shown in Fig. 8. It Some critical states shown in Fig. 8 are of multiplicity 2, or even higher for R z = 50 mm. Critical equilibrium sequences were traced for the same problem, now with parameters x = [Γ , R z ] T , and selector functions cf. Eq. (63)-and noting the different multiplicities of the sequences. Simulations used κ ε = 10 −11 , compared to common low eigenvalues of ≈ 10 −4 . The number of negative eigenvalues was checked at all obtained states as verification.
Results are shown in Fig. 9, indicating how similar critical modes are transformed by the non-sphericity.
Visibly identical critical sequences were obtained when solving for exact criticality κ ε = 0, after the introduction of the -for each critical state-relevant symmetry-preserving conditions of the form in Eq. (62). (c) Fig. 9 Critical equilibrium sequences for gas-pressurized ellipsoid membrane with constitutive constant k = − 0.05, and ambient pressure p 0 = 100 kPa. a Over-pressure p related to radial expansion at top point u R . b Gas amounts Γ related to u R . c Critical pressure related to radius R z . Thin background curves and "+" are from Fig. 8. Solid lines in c are onefold critical, dashed are twofold can be described, evaluated and interpreted. The work was partially motivated by the increasing interest in multi-stable equilibrium systems. A main starting point was that stability is a property of one particular equilibrium state, but can be evaluated as dependent on different parameter variations of not only the common load parameter. The stability conclusions were based on systematic usage of a total potential energy formulation, including constraints of different forms, and studied small, linearized vibrations around the equilibrium state, which gave the relevant eigenvalue problem. Included passive physical constraints on the solution were shown to affect the stability conclusions. In particular, the treatment shows how conclusions are dependent on which load parameter is fixed, and, thereby, which must be variable in the fulfilment of constraints.
The formulation of the constrained problems was in the form of a generalized equilibrium problem, where constraints were clearly distinguished from more general conditions, which can describe convenience transformations between different load parameters, or be selectors in the placement of evaluated equilibrium states. A generalization of the concept of stable equilibrium sequences (for a range of force values) allowed identification of stable parametric regions. The viewpoint thereby differed from common engineering practice in that it clearly separated the parameterization of an equilibrium sequence from the parameter seen as causing instability when passing a critical value.
The paper discussed settings for the evaluation of the needed eigensolutions for different forms of constraint equations, allowing efficient use of basic algorithms. One aspect was related to the inclusion of massless constraint equations. An important conclusion is that the evaluation of eigenvalues is normally straight-forward, while the corresponding eigenvectors, if needed, demand more elaborate procedures.
The present work shows how complex instability problems can be handled by a generalized equilibrium sequencefollowing algorithm, able to handle equilibrium states with multiple instabilities, and in particular cases with multiple simultaneously vanishing eigenvalues for the relevant Jacobian matrix. High accuracy is needed in all parts of the algorithm, with low tolerances for residual error, and low thresholds for what is considered as vanishing eigenvalues and tangent vector components. An improved algorithm for the treatment of the current problem classes would need more problem-oriented functions for the identification and handling of the multiple critical states. A basic idea shown in the paper is that equations enforcing the correct symmetry properties in a solution can stabilize the equilibrium iterations, but further work is needed to develop this as a general method.
Examples concerned with a simple truss structure and a gas-pressurized spherical membrane showed several aspects of parameterized stability investigations. Highly unstable equilibrium states could be found, isolated and identified. Situations where stability is regained after significant deformations of the structure were thereby amenable to analysis. Examples also showed how aspects of symmetry in a structure or a structural model can be of major importance in simulations, and how these aspects can introduce, hide or significantly modify instabilities.