(2018

Recent years have seen a research revival in structural stability analysis. This renewed interest stems from a paradigm shift regarding the role of buckling instabilities in engineering design—from detrimental sources of catastrophic failure to novel opportunities for functionality. Novel nonlinear structures take the form of optimised thin-walled structures that operate safely in the post-buckling regime; shape-morphing structures that exploit multi-stability to snap and pop between different configurations; and meta-materials that derive novel material properties from a cascade of choreographed instabilities. Hence, elastic instabilities are no longer considered as structural failures but rather exploited for repeatable well-behaved adaptations. In this article we focus on shape-morphing—a bio-inspired design strategy that intends to conform structures to different operating conditions. Computational tools that integrate easily with established methods used in industry, and that are capable of capturing the full phase diagram of compound instabilities and entangled post-buckling paths typical of these structures, are limited. Such a capability is crucial, however, as confidence in predictive tools can be key in enabling non-conventional designs. One potential candidate in this regard is generalised path-following, which combines the computational robustness of numerical continuation algorithms with the geometric versatility of the finite element method. In this paper we collate an array of successful computational tools introduced by other researchers, and introduce our own developments, to present a modelling framework fit for analysing and designing with well-behaved nonlinear structures in industry and academia. Particularly, we show that the full complexity of multi-snap events of morphing composite laminates is robustly captured by generalised path-following algorithms, and that the ability to determine loci of singular points with respect to a set of parameters is especially useful for tracing the boundaries of bistability in parameter space. Furthermore, we shed new insight into the mechanics of multi-stable laminates, showing that the multi-stability and snapping behaviour of these structures is much richer than previously assumed, featuring many unstable post-buckling branches and localised regions of stability.

Indeed, there is an ongoing trend of shifting our intuitions about structural instabilities as sources of catastrophic failure to opportunities for functionality [21].It is well known that efficient and lightweight structures are prone to structural instabilities and collapse [22,23].In this case, the predominant design philosophy is to prevent buckling or at least make its effects benign.On the other hand, when reconfigurations in shape or large elastic displacements are required, buckling is often encouraged [24].Reis [21] recently reviewed the burgeoning research effort focusing on exploiting instabilities to enable novel designs, and therefore provided a new perspective on buckling, namely from buckliphobia to buckliphilia.
In the aerospace industry, shape-morphing structures are viewed as a promising technology to enable more structurally efficient designs [3][4][5][6][7][8][9][10][11].The premise behind this idea is simple-if structures can be designed to adapt their shape to more optimally conform to different loading conditions, then structural efficiency is improved as a result.In fact, this multi-functionality has a very strong empirical proponent: nature.Birds, for example, can adapt the camber and angle of attack of their wings to different flight scenarios.Even though some of the concepts found in nature are already being exploited in aircraft structures, such as slats and flaps, they often rely on rigid load-bearing components connected to heavy hydraulic or electric actuators.This is where multi-stable structures are particularly attractive.By applying a suitable force, a multi-stable structure can be snapped from one stable state to another, thereby considerably reconfiguring its shape.Because each stable state is self-equilibrated, it does not require external energy to hold its shape, and additionally, the sensing, actuation and control functions are embedded within the nonlinear mechanics of the structure (passive control) without adding additional mass via ancillary devices.
The uptake of these novel designs in industry is partly hampered by a lack of robust computational tools tailored to the design of structures whose characteristic feature is a form of spatial chaos [25], i.e. equilibrium manifolds featuring an entire series of bifurcations that give rise to many equilibrium branches and possible loading histories, especially in cases where dynamic snap-buckling is exploited for shape adaptation.Predicting these features reliably is pushing well-established finite element techniques to their limit, and is creating an acute requirement for new computational approaches for analysis and design [21,26].In recent years the focus has been on analytical and computational techniques constrained to the analysis of very specific morphing problems with particular load cases and geometries [27][28][29][30][31][32][33][34][35][36].A drawback of these tailored approaches is that simplifying assumptions are often made, which prevent applicability to a wider range of problems.Examples include restrictions on the geometric nonlinearity to von Kármán strains (small strains, small displacements and moderate rotations); posing the problem on a domain that is cumbersome to extend beyond simple geometries; and using nonlinear stability analyses without robust branchswitching.
Due to its geometric versatility and developmental maturity, the finite element method is the preferred technique for modelling complex structural problems.Commercial finite element packages may be used to analyse multi-stable structures [37][38][39], but most of the time, these analyses are rather ad hoc, because the full taxonomy of stable and unstable equilibria cannot be revealed robustly using the quasi-static implicit solvers implemented in these codes.Rather, the engineer needs to be aware of possible bifurcation points a priori, and then "coax" the algorithm to land on a specific post-buckling mode shape using initial imperfections.Such an approach is cumbersome and requires user intervention, such that it becomes difficult and inefficient to explore the entire design nonlinear space.
To illustrate this point, an example is shown in Fig. 1, which describes a typical snap-through load-displacement path.The behaviour is complicated by the presence of a secondary equilibrium path branching from the unstable region of the equilibrium curve.In a physical experiment under rigid loading, this hypothetical structure would snap to point (c) upon reaching limit load (a).A path-following solver without the ability to pinpoint singular points and then branch-switch when needed, would generally pass the bifurcation (b) and continue to traverse along the unstable fundamental path.Not being aware of the existence of point (c), an engineer interpreting the results would then suggest that the structure snaps from point (a) to point (d), rather than from point (a) to point (c).To path-follow along the secondary equilibrium path by means of the often-used imperfection method, the analyst needs to be aware Fig. 1.Hypothetical example structure with typical snap-through equilibrium path illustrating the necessity for robust branch-switching.Without the ability to detect bifurcations and branch-switch, an engineer would not be aware of point (c) and presume that the structure snaps to point (d).Source: Adapted from [33].
of the existence of bifurcation point (b), stop the path-following algorithm before reaching it, run a linear eigenvalue analysis and then apply the lowest eigenmode as an imperfection.Such a procedure is cumbersome to implement and restricted to simple bifurcation points (for compound bifurcations there is little control over which branch the solver will converge to).Finally, when conducting design sensitivity studies, the full equilibrium paths of different design iterations need to be traced, when the designer is, in fact, only interested in how certain points, e.g.critical points (a) and (b), vary with a particular design parameter.
Thus, the ideal computational framework for structures exploiting instabilities includes, but is not limited to, the following characteristics: 1. Applicable to arbitrary (thin-walled) geometries.2. Applicable to large displacements and rotations, e.g. total Lagrangian framework, such that a large variety of structural instabilities are accounted for.3. Able to detect singular points and branch-switch onto secondary paths without recourse to initial imperfections.4. Allows for rapid parametric studies of critical points with respect to any geometric, constitutive or secondary loading parameter.5. Is readily integrated with accepted computational methods used in industry, predominantly, the finite element method.
One possible framework that addresses all of these points is the so-called generalised path-following method, which provides the means to systematically explore the design space of multi-stable structures with respect to a given parameter set.The term generalised path was probably coined by Eriksson [40], where it referred to a multi-parametric setting of static equilibrium problems that extends the notion of applied load as the only active control parameter, and it is this definition that we refer to herein.Historically, generalised path-following has been used extensively in the fields of applied mathematics and physics [41][42][43][44][45][46], where the term numerical continuation is a more common designation.In engineering applications, however, path-following in load-displacement space, typically a variant of the Riks method [47], is a familiar term, such that generalised path-following is an intuitive extension to a multiparametric setting.
In the 1960's, Sewell introduced the notion of the "equilibrium surface" [48,49] -a surface whose shape could be used to identify the stability of a structure with respect to changes in the governing parameters (see Fig. 2)to structural mechanics.With the advent of catastrophe theory [51], interest in applying the topological concept of an equilibrium surface to elastic stability intensified, mostly in an analytical setting [52][53][54], but a generalised computational framework was not introduced until the seminal papers by Rheinboldt [55][56][57] in the 1980's.These Fig. 2.An equilibrium surface for the snap-through behaviour of an elastic arch.In a generalised path-following algorithm different paths on this equilibrium surface can be traced by defining particular auxiliary equations.For example, the figure shows a fundamental path in load-displacement space (height held constant), a parametric path in height-displacement space (load held constant), and a singular path in load-height-displacement space connecting the limit points for different arch heights.Source: Adapted from [50].concepts allowed pinpointing of singular points (bifurcation and limit points); branch-switching at bifurcation points; path-following with respect to any parameter, e.g.load, thickness, Young's modulus, dimensions, etc.; and tracing loci of bifurcation and/or limit points in multi-parametric space (see Fig. 2).In this sense, generalised path-following extends concepts from computational bifurcation theory [58][59][60], such as path-tracing, pinpointing of singular points and branch-switching, from a single-to a multi-parametric setting.Starting from the mid-1990's, Eriksson and co-workers [40,[61][62][63][64][65] established themselves as the main developers of generalised path-following, presenting numerous examples in structural engineering where the approach proved to be of great benefit, while also providing details on how the technique could be incorporated into commercial nonlinear finite element codes.Similarly, within the broader applied mathematics community, the Library of Continuation Algorithms (LOCA) was developed by Sandia National Laboratories to deal with large degree-of-freedom problems, typically arising from the discretisation of partial differential equations but not necessarily restricted to finite elements, to be run on parallel computing platforms [44].In engineering, numerical continuation has also been successfully applied to study nonlinear phenomena in the fields of aerodynamics [66][67][68], aeroelasticity [69][70][71] and the analysis of cylindrical shells [33,[72][73][74].
The purpose of this paper is to adapt the ideas and concepts presented particularly by Eriksson [40] and Salinger et al. [44] to the thermo-mechanical setting typical of many multi-stable structures; especially, for laminated orthotropic structures used for shape-morphing in the aerospace industry [75].As such, the current work is an amalgamation of computational techniques that have either been adapted or developed further to suitably apply them to the design of thermo-mechanical multi-stable structures.Section 2 provides the theoretical background in terms of the general setting (Section 2.1), together with implementation details for following paths with one changing parameter (Section 2.2), and its extension to two simultaneously varying parameters (Section 2.3).As a canonical example, Section 3 demonstrates the general capabilities of the algorithm using a simple problem.In Section 4 we show that generalised path-following is particularly well-suited to modelling morphing structures.This is demonstrated by means of two classical example problems from the literature-a flat nonsymmetric cross-ply laminate that is cooled-down from cure to room temperature and then snapped between two stable states (Section 4.1), and the effect of initial curvature on the bistable behaviour (Section 4.2).We provide new physical insight into the mechanics of the snap-through phenomenon, showing that the snapping behaviour is much richer than previously assumed, featuring many entangled secondary post-buckling paths with localised regions of stability, and using the bifurcation-tracking capability to delineate regions of bistability in parameter space.

Theory
A generalised path-following algorithm combines the mathematical domains of finite element analysis and numerical continuation.The mathematical methods used in numerical continuation are well established [41][42][43]45,56], but are not classically used for structural mechanics applications, where specialised arc-length techniques are predominant [47,76].The present formulation considers a discretised model of a slowly evolving, conservative and elastic structure, where the internal forces and tangential stiffness are uniquely defined from the current displacements by means of the first and second variations of an energy potential.Thus, non-conservative loading and historydependent problems such as plasticity are not included in the current formulation.
In Section 2.1 we proceed with the general framework and then discuss implementation-specific details in Sections 2.2-2.3.The interested reader is also especially encouraged to consult the excellent expositions by Eriksson [40,61,62] which form the basis of the theory.Finally, an accurate analysis using these numerical methods depends on suitable choices of nonlinear beam and shell elements that have sufficient fidelity to capture the full complexity of nonlinear instability phenomena (see [77,78]).

The general setting
In classical structural mechanics applications, equilibrium is expressed as a balance between internal and external forces, where, in a displacement-based finite element setting, this balance is written in terms of n discrete displacement degrees-of-freedom, u, and a scalar loading parameter, λ, The vectors p(λ) and f (u) are the external (non-follower) load and internal force, respectively.In the case of linear and proportional loading we have p(λ) ≡ λp ,λ 1 = λp, where p is a constant reference loading vector (dead loading).This system of n equations in (n + 1) unknowns -n displacement degrees-of-freedom and one loading parameter -is then solved for a solution point, x = (u, λ), by defining an additional scalar arc-length constraint, where n u and n λ take different forms depending on the nature of the arc-length constraint.By linearising about the current equilibrium state, x, and applying Newton's method for the iterative correction, δx, we can find a set of solution points that describes a continuous equilibrium curve.Note that the partial derivative of the residual with respect to the displacement vector, F ,u = f ,u (u), is equal to the tangential stiffness matrix K T (u).For generalised path-following, Eq. ( 1) is adapted to incorporate any number of additional parameters, where ⊤ is a vector containing p control variables.Λ 1 corresponds to parameters that influence the internal forces (e.g.material properties, geometric dimensions, temperature and moisture fields) and Λ 2 relates to externally applied mechanical loads (e.g.forces, moments, tractions).
The n number of equilibrium equations in Eq. ( 4), correspond directly to the n number of displacement degrees-offreedom in the system.Because the structural response is parametrised by p additional parameters, a p-dimensional solution manifold in R (n+ p) exists-the so-called equilibrium hypersurface [48].By defining additional auxiliary equations, g, specific solution subsets on this p-dimensional solution manifold are defined.Hence, we wish to evaluate solutions to the augmented system When r auxiliary equations are defined, the solution to Eq. ( 5) is ( p − r )-dimensional.Hence, r = p − 1 auxiliary equations are required to define a one-dimensional curve, or so-called subset curve [62], in R (n+ p) .
Posing the problem in this manner allows the structural response to be viewed not only as a function of a varying load but also as a function of other parameters that define the structure.By treating these additional parameters as "forcing" variables in an arc-length solver, the effect of these parameters on the structural response can be obtained.Hence, the computationally expensive approach of studying variations in geometry and material properties, by evaluating full load-displacement equilibrium curves for each additional model, is avoided.Instead, different load-displacement or parameter-displacement paths are described as cutsets of a higher-dimensional solution surface (see Fig. 2).
This treatment naturally leads to the notion of tracing loci of singular points in parameter space, such as the maxima and minima shown in Fig. 2. As the designer of multi-stable structures is predominantly interested in critical instability points, e.g.limit points that initiate snap-through or symmetry-breaking bifurcations, these singular curves are invaluable for rapidly exploring the design space.To constrain the system of n equilibrium equations to a locus of singular points, we simultaneously enforce the fulfilment of a criticality condition, for example K T φ = 0, i.e. at least one eigenvector φ of the tangential stiffness matrix K T spans the nullspace.In the most general form, not limited to but including the previous criticality condition, a vector of q auxiliary variables, v, may be added to the auxiliary equations g, Hence, Eq. ( 6) describes n equilibrium equations and r auxiliary equations in (n + p + q) unknowns leading to a ( p + q − r )-dimensional solution.To determine a one-dimensional subset curve of singular points, we thus require r = p +q −1 auxiliary equations to constrain the system.Following the example from above, when the n-dimensional null vector at the critical state is introduced as the auxiliary variable, v, a singular subset curve in two parameters, p = 2, is appropriately constrained by the associated r = n + 1 auxiliary equations K T v = 0 and ∥v∥ 2 = 1, where the scalar equation restricts the magnitude of the eigenvector.When evaluating one-dimensional subset curves (r = p + q − 1), one additional constraining equation is needed to uniquely solve the system of equations for a solution point y = (u, Λ, v) on the curve described by G(y).Hence, where N is a scalar equation that plays the role of a multi-dimensional arc-length constraint along a specific direction of the subset curve.Note that the system of equations for classical load-displacement equilibrium paths can be recovered by setting p = 1 and q = r = 0.A solution to Eq. ( 7) is determined by a consistent linearisation coupled with Newton's method, where the superscript denotes the jth equilibrium iteration and the subscript the kth load increment.The iterative correction cycle is typically started by a predictive forward Euler step.For most problems, the inversion of the iteration matrix, is significantly simplified by partitioning the system into blocks such that only the symmetric tangential stiffness matrix K T ≡ F ,u needs to be factorised in the solution process (see Appendix and e.g.[79][80][81]).
Within this framework, the meaning of the term generalised path-following becomes clear.It refers to the fact that any arbitrary curve can be traced on the equilibrium surface, as long as a pertinent auxiliary equation is defined that constrains the equilibrium equation to the locus of points required.These auxiliary equations can define, but are not necessarily limited to, the following interesting paths: • Classic equilibrium paths in load-displacement space (a loading parameter is varied).
• Parametric paths in parameter-displacement space (a geometric, constitutive or secondary loading parameter is varied).• Pinpointing singular points (bifurcation and limit points) on either of the two paths mentioned above.
• Bifurcated branches emanating from a bifurcation point.
• Singular paths that describe a locus of bifurcation and/or limit points in load-parameter-displacement space.
• Branch-connecting paths that connect points on distinct equilibrium curves.
Furthermore, auxiliary equations have also been devised to define optimality criteria [82] and to directly evaluate the imperfection mode shape that leads to the greatest knock-down in the buckling load of compressed cylindrical shells [83].

Path-following in one parameter
The parameters Λ = (λ 1 , . . ., λ p ) can be any set of p parameters that affect the equilibrium of the structure.In a structural engineering context, it is pertinent to treat parameter λ 1 as the fundamental loading parameter, here mechanical or thermal in nature, such that path-following in λ 1 traces the fundamental load-displacement equilibrium path of an idealised baseline model.Parameters Λ s = (λ 2 , . . ., λ p ) are then added as secondary parameters that perturb this idealised baseline problem, e.g.variations in geometry, changes to constitutive properties or the addition of secondary loadings.Any saved equilibrium solution on the fundamental path then serves as a starting point for path-following in a secondary parameter, λ s ∈ Λ s , thereby tracing a parametric path of perturbed equilibrium points with λ 1 constant.
Due to the practical limitation of visualising results in three dimensions, the present implementation is posed in terms of the displacement degrees-of-freedom, u, the fundamental loading parameter, λ 1 , and a secondary parameter, λ s ∈ Λ s , such that a pertinent norm of the displacement vector, e.g.u i for i = 1 . . .n, can be plotted against (λ 1 , λ s ).The results for different combinations of (λ 1 , λ s ) are therefore plotted in succession.For clarity, fundamental and parametric paths are first discussed individually in Sections 2.2.1 and 2.2.2, respectively, with all equations expressed in terms of (λ 1 , λ s ) only, and all other parameters removed from consideration (kept constant at baseline values).Section 2.2.3 then generalises the notation to any non-critical subset path with only one varying parameter.The topic of pinpointing singular points on these paths is then discussed in Section 2.2.4.

Fundamental paths
The system of equations for classical load-displacement equilibrium paths can be recovered by setting p = 1 and q = r = 0 such that λ 1 becomes the fundamental loading parameter.In the case of linear and proportional mechanical loading, In the case of thermal loading, λ 1 affects the internal force vector and p1 = 0, such that Introducing the scalar arc-length constraint N (u, λ 1 ), as defined in Eq. ( 2), and applying Newton's method where δu and δλ 1 are iterative corrections of the nodal displacements and loading parameter, respectively.F ,u ≡ K T is the tangential stiffness matrix and F ,λ 1 = −p 1 is the reference loading vector for linear and proportional mechanical loading.When mechanical loading is applied in terms of nodal forces and moments or element-wise tractions, the reference force vector p1 is readily assembled from the input data.In the case of non-zero prescribed displacements, the reference force vector is derived from p1 = −K C u c , i.e. the constrained portion of the tangential stiffness matrix multiplying the non-zero prescribed displacements.As the tangential stiffness matrix is generally a function of current displacements, the reference force vector p1 derived from prescribed displacements needs to be updated for every iteration of the solution scheme.
In the case of thermal loading, the effective strain method described by Parente Jr et al. [84] is applied.The effective strain is given by the difference between the classic geometric strain, here the Green-Lagrange strain, ϵ GL , and the free thermal strain, ϵ th .The thermal strain, represents the strain induced when thermal expansion of the material is not constrained, with α representing the thermal expansion coefficient vector and ∆T the change in temperature.For linear elastic materials, the stress-straintemperature relation can be written as where σ is the energetically conjugate second Piola-Kirchhoff stress tensor and ∆ T is a temperature change from a given strain-free reference temperature.In general, the tangent constitutive tensor, C T , is temperature dependent.
In the context of the principle of virtual displacements the internal force vector, used in Eqs. ( 10)- (11), is given by By assuming a general relation between the virtual geometric strain and the virtual displacements, δϵ GL (u) = B(u)δu, and substituting Eq. ( 14) for σ , the internal force vector for a thermal loading problem is The precise definition of the kinematic matrix, B, depends on the chosen finite element implementation and shape function interpolation, and is therefore not discussed in detail herein.The tangential stiffness matrix, K T , and the thermal load vector, pth , are derived from the variation of Eq. (11): In particular, the variation of Eq. ( 16) with respect to u gives: Eq. (18) shows that the elastic stiffness matrix, K e , is unchanged from the classical linear stiffness matrix and not affected by temperature changes.Even though the geometric stiffness matrix, K g , always depends on the particular element formulation chosen, Eq. (18) shows that K g is a function of temperature changes via the current stresses, σ .This means that the tangential stiffness matrix is only affected by temperature changes through the geometric stiffness component.
Similarly, the thermal load vector is given by the variation of Eq. ( 16) with respect to λ 1 , where we have assumed that the tangential constitutive tensor C T does not vary with temperature.If the constitutive tensor does vary with temperature then the term C T,λ 1 needs to be computed.In this case it is most convenient to evaluate an approximation to pth using a forward difference scheme, where ε is a small perturbation parameter typically in the range of 10 −5 to 10 −8 .The default choice of ε in Eq. ( 20), and many later forward difference approximations herein, is 10 −8 .In the authors' experience, the particular choice of ε within the range is robust, meaning that the convergence of Newton's method is not sensitive to changes in ε.Even though the accuracy of this numerical approach depends on the choice of ε, and is computationally less efficient than a direct approach, it has the advantage of being valid for all element types, and is therefore readily implemented in existing codes.In either case, it is important to remember that the thermal load vector needs to be iteratively updated during the entire thermal loading process, because B is a function of the current displacements u.
Finally, the form of the arc-length constraint depends on the particular method chosen, although in our experience, the planar constraint by Riks [47] and the cylindrical constraint by Crisfield [76] work most robustly.The details of different arc-length constraints are covered extensively in textbooks and many publications (see e.g.[84][85][86][87]), and are therefore not elucidated in further detail here.Some implementation-specific comments are nevertheless warranted.
• Given the quadratic nature of Crisfield's cylindrical arc-length constraint [76], two roots for the iterative correction δλ j 1k (of the jth equilibrium iteration and the kth load increment) arise, and if both of these are imaginary, or real and negative, it is easiest to repeat the current load increment with a reduced incremental arc-length.
• When both roots are real and positive, the root that is closest to the root of an analogue linear system (simply ignoring the quadratic term) is chosen.When one of the real roots is positive and the other negative, the positive real root is chosen.• The sign of δλ 1 1k , i.e. the load predictor in the first iteration of each increment k, is readily determined by taking the dot product between the total incremental displacement of the previous step, ∆u k−1 , and δu 1 k , i.e. the displacement predictor in the first iteration of the current step calculated from the unit load vector 1 1k is defined to take the sign of this dot product, then limit points are easily traversed.Noting these caveats, Crisfield's cylindrical arc-length constraint [76] performs very robustly for both fundamental and parametric paths, including the traversal of turning and snap-back points, and is therefore explicitly used for the problems studied in Section 3 and in Section 4.

Parametric paths
Any equilibrium solution on a fundamental path can be used as a starting point for path-following in one of the secondary parameters, λ s ∈ Λ s .In this manner, the sensitivity of the baseline design to variations in λ s is assessed.For example, if the basic topology of the problem is to be kept constant we may parametrise a specific geometric dimension as L(λ s ) = L 0 (1 + cλ s ), which can be used to analyse the effects of thickness, length or height on the structural behaviour.Equally, we may add a mode shape, β, to the nodal coordinates of the basic geometry, x 0 , such that x 0 (λ s ) = x 0 + cλ s β, which can be useful for imperfection sensitivity studies or parametrically varying the topology of the problem.Similarly, we may study the effects of constitutive properties, such as Young's modulus E(λ s ) = E 0 (1 + cλ s ).Finally, a secondary loading p s = cλ s ps may be added to study the combined load case p = λ 1 p1 + cλ s ps , which may be interpreted as studying the effect of either loading under the perturbing influence of the other.In all cases, it may be necessary to scale the secondary parameter, λ s , by a factor, c, such that the order of magnitude of the two parameters is similar, i.e.O(λ 1 ) ≈ O(λ s ).
The methods introduced in the previous section can readily be adapted to parametric paths.In the current notation there are three different combinations of λ 1 and λ s , with λ 1 held constant in each case.For combined mechanical loading Alternatively, for fundamental mechanical loading with secondary thermal loading, geometrical changes or constitutive variations Finally, for fundamental thermal loading with secondary thermal loading, geometrical changes or constitutive variations Introducing the scalar arc-length constraint N (u, λ s ) and applying Newton's method where δu and δλ s are iterative corrections of the nodal displacements and secondary parameter, respectively.For secondary linear and proportional mechanical loading F ,λ s = −p s .Otherwise, F ,λ s is calculated from a forward difference scheme, where ε is a small perturbation parameter typically in the range of 10 −5 to 10 −8 .Finally, the caveats regarding Crisfield's cylindrical arc-length constraints raised in the previous section are equally valid here.

Single-parameter non-critical subset paths
To illustrate the general concepts of continuing in load-displacement and parameter-displacement space, the equations in the previous two sections are purposely restricted to fundamental parameter λ 1 and secondary parameter λ s .These concepts are now combined to generalise the notation to any path with only one varying parameter (fundamental or parametric), also known as a non-critical subset paths in one parameter [40].Hence, for fundamental mechanical loading, and for fundamental thermal loading.Ps is a matrix of column-wise secondary load vectors, and Λ i s and Λ e s parametrise internal and external force vectors, respectively.Each parameter takes the value that describes a baseline design and then one target parameter, , is varied to trace a particular equilibrium path, while all other parameters where Σ is the set of prescribed constant parameters and N (u, where 1 is the identity matrix.The 2nd row and 3rd column of the iteration matrix can generally be omitted as δΛ c = 0 by definition.Similar to classical arc-length methods [76], the iteration matrix in Eq. ( 29) is never inverted in its entirety but split into the blocks shown, and then solved via a partitioning procedure and back-substitution such that only F ,u ≡ K T needs to be inverted [76] (see Appendix).

Pinpointing singular points
Pinpointing singular points is useful for evaluating the exact value of snap-through loads (to within a certain tolerance) and for ascertaining the existence of bifurcations onto other branches.Furthermore, unfolding of these singular points with respect to other parameters can provide invaluable insights into imperfection and design parameter sensitivity.
Remark (Singular Points).Using a Taylor series expansion, a small change in the energy potential, Π (u, Λ), of a conservative elastic system is with all parameters Λ held constant.By definition, F(u, Λ) = 0 for equilibrium, such that the sign of δu ⊤ K T δu is a sufficient condition for ascertaining the stability of the equilibrium state for small perturbations, δu.Assuming a symmetric tangential stiffness matrix, K T , a singular point on an equilibrium path with target parameter λ t ∈ Λ and unvaried parameters Λ c , may thus be determined by δu ⊤ K T δu = 0 for all arbitrary δu.This coincides with the condition that Here, µ = 0 is the eigenvalue that corresponds to the critical eigenvector φ at the singular point (u * , λ * t , Λ * c ). Parametrising the equilibrium equations further in terms of an arc-length parameter, s, and differentiating with respect to this curve parameter gives where a superimposed dot denotes differentiation with respect to s.Given the symmetry of K T and the singularity condition Eq. ( 31), pre-multiplication of Eq. ( 32) by φ ⊤ at a singular point yields Hence, on an equilibrium path in one parameter (fundamental or parametric) there can only be two types of singular points-a limit point, i.e. a local extremum with λt = 0 and φ ⊤ F ,λ t ̸ = 0, or a bifurcation point, i.e. an intersection between two or more distinct equilibrium curves with φ ⊤ F ,λ t = 0.For numerical consistency, the latter condition is typically implemented as The choice of ε is not calculated automatically but from experience based on a "good guess" of 10 −3 .From the equilibrium curve, it is typically straightforward to ascertain if a critical point should be a limit point or not, and in the rare cases where the choice of ε = 10 −3 is not sufficient, the tolerance can be reduced ad hoc.
The second derivative F = 0 can be used to determine the two tangent vectors at a simple bifurcation point -one tangent to the primary path and the other tangent to the secondary bifurcated path -and also to classify the type of bifurcation-pitchfork, transcritical, isola formation point or cusp point.The disadvantage is that the derivative K T,u needs to be computed approximately via finite differences, with an associated penalty in computational cost.For computational efficiency, these calculations are not performed herein and the interested reader is directed to Refs.[84,88].
To pinpoint singular points, an augmented system of the form described in general by Eq. ( 7) is formulated.A number of different auxiliary equations are defined in the literature for this purpose, e.g.g = µ = 0, g = det(K T ) = 0, g = K T φ = 0, etc., as outlined in the comparative review by Melhem & Rheinboldt [89].The advantage of this method is that the singularity condition forces Newton's method to converge to the singular point directly in a single increment.Although bisection techniques have also been developed [59,88], they require the calculation of multiple intermediate equilibrium points to hone in on the singularity and are thus computationally more expensive.
In the present finite element setting, the nullvector approach described by Wriggers et al. [80,81] -generally first presented by Seydel [90] and Moore & Spence [91] -and the minimally augmented method introduced by Griewank & Reddien [92], and further developed by Eriksson [61] and Battini et al. [93], are used.The minimally augmented technique is chosen by default in the developed computer program, and the nullvector method used as an alternative option when the former has trouble converging.
The nullvector method is based on the fact that the tangential stiffness matrix, F ,u ≡ K T , has at least one zero eigenvalue, µ = φ ⊤ K T φ, at a singular point.Therefore, the associated eigenvector, φ, is in the nullspace of K T .Thus, the augmented system is written as where the norm of the nullvector is required to eliminate the trivial solution φ = 0. Eq. ( 34) features (2n +1) equations in (2n + p) variables, and the ( p − 1) extra equations required to solve the system are implicit in the definition that the added control parameters in Λ c are held constant, i.e.Λ c j = Σ j for j = 2 . . .p (see Eq. ( 29)).The resulting system of equations can therefore be solved in the usual manner via Newton's method for a singular point, (u * , λ * t , Λ * c ), as well as the associated nullvector, φ.The nullvector can then be used to distinguish between limit and bifurcation points via Eq.(33).The solution process follows the typical predictor-corrector scheme with the iteration matrix derived from the linearisation of Eq. ( 34), Following the reasoning by Wriggers & Simo [81], approximate directional derivatives of the tangential stiffness matrix can be computed as follows, with ε in the range of 10 −5 to 10 −8 .As noted before, a seemingly robust default choice is ε = 10 −8 , although it is possible that at some bifurcation points Eqs. ( 36)-( 37) may become sensitive to changes in ε, and in these cases a central difference scheme may be more appropriate for good convergence of Newton's method.When λ t is a primary or secondary displacement-independent mechanical loading parameter, (K T φ) ,λ t vanishes and F ,λ t equals −p 1 or −p s .Otherwise, F ,λ t is approximated by When solving the system of equations in Eq. ( 35), the iteration matrix is not inverted in its entirety, but split into individual blocks (as shown), and solved by a partitioning procedure in such a manner that only the symmetric tangential stiffness matrix needs to be factorised (see Appendix and Wriggers & Simo [81] for the general algorithm).This means that the computational cost of pinpointing is similar to that of a standard arc-length continuation step, especially for large degree-of-freedom systems, where the factorisation of K T dominates.The present computer implementation proceeds as follows.While continuing along a fundamental or parametric equilibrium path, the 20 smallest magnitude eigenvalues of the tangential stiffness matrix are monitored.This can be done efficiently for large yet sparse tangential stiffness matrices by using the FORTRAN Arnoldi package ARPACK [94], embedded in many popular numerical computing environments such as MATLAB and SCIPY.When the number of negative eigenvalues between two consecutive converged equilibrium solutions changes, a singular point must exist between these two converged equilibria and the pinpointing procedure is started.The number of singular points present depends on the change in the number of negative eigenvalues, N * .The set of eigenvectors, Φ, associated with the smallest N * eigenvalues at the last converged equilibrium state, (u l , λ l t , Λ l c ), are then extracted, and each φ j ∈ Φ for j = 1 . . .N * consecutively seeded alongside (u l , λ l t , Λ l c ) as the starting point for the iterative pinpointing procedure.In our experience, this eigenmode seeding procedure works reliably for one singular point, as well as multiple distinct or coincident (compound) singular points.If the solver does not converge, then an additional equilibrium point between the two previously determined equilibria is determined and the process is repeated.This is typically the case when the nonlinearity between the last converged equilibrium state, i.e. the starting point, and the singular point is too high.
In the minimally augmented method [92], the vector equation K T φ = 0 is replaced by the scalar equation where µ is an eigenvalue of the symmetric tangent stiffness matrix, K T , and φ the associated eigenvector.By definition, µ must vanish at a singular point such that φ becomes a nullvector of K T .The scalar equation ( 39) is derived by pre-multiplying the eigenvalue problem, K T φ = µφ, with φ ⊤ and enforcing the normalisation constraint φ ⊤ φ = 1.Rather than computing the singular point (u * , λ * t , Λ * c ) and the associated nullvector φ simultaneously (as is done in the nullvector method), µ and φ are determined separately.Thus, as a first step the eigenvalue problem is written as an iterative problem [93] where the latter equation (40b) now constrains the updated eigenvector to be parallel to the previous one with Note, by pre-multiplying Eq. (40a) with φ ⊤ k+1 and heeding Eq. (40b), we recover the scalar condition µ k+1 = φ ⊤ k+1 K T φ k+1 .Eqs. (40a)-(40b) are rearranged into matrix form as follows: Eq. ( 41) is solved efficiently using a partitioning procedure by first computing an immediate vector ϕ = K −1 T φ k , and then back-substituting to find µ k+1 = 1/(φ ⊤ k ϕ) and subsequently φ k+1 = µ k+1 ϕ.With an arbitrary random starting vector, φ 0 , the iteration scheme converges to the smallest magnitude eigenvalue and associated eigenvector typically within 2-3 iterations.When starting the algorithm with an informed, non-arbitrary approximation of a particular eigenvector, the solver refines the eigenvector and returns the associated eigenvalue often within a single iteration.
As a second step, the equilibrium equations are augmented with the scalar singularity condition of Eq. ( 39), Eq. ( 42) features (n + 1) equations in (n + p) variables, and the ( p − 1) extra equations required to solve the system are again implicit in that parameters Λ c take prescribed values.A singular point (u * , λ * t , Λ * c ) is calculated via Newton's method using the iteration matrix The differentials of the scalar singularity condition are computed via finite difference approximations, with ε in the range of 10 −5 to 10 −8 .As mentioned previously, if the choice of ε = 10 −8 does not provide good convergence for some bifurcation points, then a central difference scheme may be used for better results.When λ t is a primary or secondary displacement-independent mechanical loading parameter, µ ,λ t vanishes and F ,λ t equals −p 1 or −p s .Otherwise, F ,λ t is approximated by Eq. (38).Eq. ( 43) is solved in a typical manner using a partitioning procedure by first determining δu 1 = −K −1 T F and δu 2 = −K −1 T F ,λ t , and then back-substituting to find δλ t = − µ+µ ⊤ ,u δu 1 µ ,λ t +µ ⊤ ,u δu 2 and subsequently δu = δu 1 + δλ t δu 2 .
The minimally augmented method therefore requires two iteration steps to be performed in succession-the main iteration loop performed on Eq. ( 43) to converge the solution (u, λ t ) closer to a singular point, followed by an embedded iteration loop via Eq.( 41) to update the eigenvalue and eigenvector for each iteration of (u, λ t ).The pinpointing algorithm is again initiated by consecutively seeding a particular eigenvector φ 0 = φ j ∈ Φ, j = 1 . . .N * in the set of eigenvectors, Φ, associated with the smallest N * eigenvalues of K T that changed sign between the last two converged equilibrium solutions (see description for the nullvector method on page 13).Hence, the N * eigenvalues and eigenvectors at the last converged equilibrium state, obtained in a computationally efficient manner from ARPACK [94], serve as starting points for pinpointing a set of N * singular points.

Branch-switching at bifurcations
There are many ways of branch-switching between two equilibrium paths that intersect at a bifurcation point (see [59]).In the case of a simple degeneracy, i.e. exactly one zero eigenvalue with two paths intersecting at a bifurcation point, switching from one path to the other is relatively straightforward.The simplest method, based on the notion of the nullvector, 2 is to inject the critical eigenvector into the displacement field at the bifurcation point [58].Thus, the critical eigenvector at the bifurcation point, φ, is used as a perturbation to the solution, u * , at the bifurcation point as follows: such that the perturbed configuration u p is now used as the predictor (with δΛ = 0) for the first step on a new path starting from the bifurcation point.The magnitude of the scaling factor, ζ , is determined from where the sign of ζ controls the direction of path-following along the bifurcated path, and ρ is a problem-specific constant in the range of 1-100.If ρ is too small, then the algorithm may continue on the primary path, and if ρ is too large, then the solution may not converge.Thus, it is beneficial to implement a restarting function into the algorithm that varies the magnitude of ρ for certain scenarios, although in our experience, ρ = 100 works reliably for most cases.At a compound singularity, i.e. exactly two or more zero eigenvalues, the situation is more complicated.Within the framework of the eigenvector insertion method, Wagner & Wriggers [58] propose a linear combination of the critical eigenvectors to perturb the known equilibrium state, where φ j is the jth eigenvector of a total of N * c compound/coincident critical eigenvectors.The drawback of this approach is that there is no general method for determining different combinations of the scaling factors, ζ j , such that the solver lands on all paths branching from the bifurcation point.Numerous methods for dealing with compound bifurcations have been proposed [95][96][97], although the perturbation method by Huitfeldt [98] seems to be the most promising, as it can determine the full set of equilibrium states in the vicinity of a bifurcation point by means of a simple algorithm that integrates nicely into the generalised path-following framework (see Section 2.3.1).

Path-following in two parameters
Section 2.2 restricted the notation to a single varying parameter, λ t , such that path-following was only possible along fundamental or parametric paths.It is often useful and/or computationally efficient to vary two parameters simultaneously.The example of a branch-connecting path, alluded to in the previous Section 2.2.4, is particularly useful for determining all bifurcation branches at compound bifurcation points, but can also be used to uncover adjacent equilibrium paths such as broken pitchforks.Evaluation of a branch-connecting path, described in Section 2.3.1, requires continuation in two parameters simultaneously-a fundamental and a perturbing load.Furthermore, computing the sensitivity of singular points with respect to a second parameter, be it a secondary loading, geometry or constitutive property, is achieved more efficiently by tracing the locus of singular points with respect to the secondary parameter, than by evaluating a set of full fundamental paths for different values of the secondary parameter.Such curves are known as foldlines or critical subset paths [61], and the exposition here distinguishes between tracking general singular points in Section 2.3.2,valid for both limit and bifurcation points, and bifurcation points only in Section 2.3.3.Finally, the theory concludes in Section 2.3.4 with the computation of tangent spaces needed for the predictor step of Newton's method, which is trivially done for non-singular problems, such as branch-connecting paths, but more intricate for foldlines.The following sections on foldlines follow closely from Eriksson [61], although more recently Moghaddasie & Stanciulescu [99] and Rezaiee-Pajand & Moghaddasie [100] have proposed similar algorithms using different auxiliary singularity conditions.

Branch-connecting paths
Huitfeldt [98] introduced the notion of using a perturbing load vector to search for adjacent equilibrium states of a baseline system.This approach is very useful for determining all secondary bifurcation paths branching from a compound bifurcation, but can equally be used to determine non-trivial solutions that are not connected to a known Fig. 3. Two equilibrium surfaces parametrised by a fundamental loading parameter, λ 1 , and a perturbing loading parameter, λ s , intersecting to produce a simple bifurcation.All equilibrium states of the unperturbed system within the vicinity of the bifurcation point, y * , are determined by tracing the intersection of the solution manifold and a hypersphere describing a region around the bifurcation point.primary path.For the simple case of a fundamental mechanical loading parameter, λ 1 , and a perturbation parameter, λ s , a perturbed equilibrium state is defined as follows: where ps is an arbitrary perturbing vector, not in the same direction as the fundamental loading vector, i.e. p1 •p s It is obvious from Eq. ( 49) that the unperturbed equilibrium equation is recovered for λ s = 0.The intersection of the two-dimensional solution manifold, described by Eq. ( 49), and a hypersphere of radius γ centred around a bifurcation point, y * = (u * , λ * 1 , 0), is a closed one-dimensional curve (see Fig. 3).This curve can be traced by defining the hypersphere as an auxiliary equation, This system can be solved via Newton's method once a pertinent arc-length constraint, As this curve is followed, additional equilibrium solutions of the unperturbed system, either on a fundamental or bifurcated equilibrium path, are determined every time λ s changes sign (an approximate solution to λ s = 0 can be determined by linear interpolation).Due to the fact that this curve connects solutions on adjacent, yet distinct equilibrium curves, it is often called a branch-connecting path.
In an analogous way to Crisfield's spherical and cylindrical arc-length constraints [76], Eriksson [40] reports that the simplified cylindrical auxiliary equation g = (u − u * ) ⊤ (u − u * ) − γ 2 performs more robustly than the spherical equation.Eriksson [40] further recommends choosing one of the critical eigenvectors at the bifurcation point as the perturbing vector, ps , and taking a value γ ≈ 0.1∥u * ∥ 2 as the radius.
The advantage of the branch-connecting approach is that the computational expense is basically the same no matter how many additional equilibrium paths exist.Furthermore, the same algorithm can also be used to determine the broken-away curve of a broken pitchfork bifurcation.In fact, the method can generally be used to uncover adjacent equilibrium paths that do not intersect with a primary loading path.

General foldlines
Path-following of foldlines could be posed as a general two-parameter problem with target parameters Λ t = (λ t1 , λ t2 ), Λ t ⊂ Λ.In a structural mechanics setting, however, our main interest is to compute the sensitivity of critical points on a fundamental load-displacement path with respect to an added parameter.The equilibrium equations, F(u, λ 1 , λ s , Λ c ) = 0, are therefore written in terms of the active parameters λ 1 and λ s -the former being a fundamental loading parameter (thermal or mechanical), and the latter parametrising a secondary loading, a change in geometry or a variation in constitutive properties.The baseline problem corresponds to the baseline secondary parameter, λ s = λ b s , and all other parameters, Λ c ⊂ Λ, that define the problem are held constant at their baseline values.The choice of secondary parameter, λ s ∈ Λ s , can be changed one at a time, such that a system with parameters p > 2 is framed as a collection of individual problems, each with p = 2.
Depending on the nature of λ 1 (mechanical or thermal loading) and λ s (secondary loading or changes to the structure), the equilibrium equations can be written in one of three forms as defined by Eqs. ( 21)- (23).The foldline algorithm then constrains the equilibrium equations to a curve that unfolds a baseline singular point with respect to another parameter, i.e. a curve describing λ * 1 = λ * 1 (λ s ).The unfolding of limit and bifurcation points often results in a sequence of singular points of the same nature, limit and bifurcation, respectively.For the classic pitchfork bifurcation of the elastica, however, small geometric imperfections can break the pitchfork, thereby transforming the bifurcation point into a limit point.This means that for specific values of the secondary parameter, the unfolding of limit and bifurcation points can lead to sequences of singular points of the opposite nature.The foldline algorithm presented in this section can handle both types of singular points, such that sequences of limit points only, bifurcation points only and combinations of limit and bifurcation points can be traced.
In the present computational implementation both the nullvector and minimally augmented methods (see Section 2.2.4) are used to formulate the auxiliary singularity condition required to trace along a foldline.The difference to the pinpointing procedure for individual singular points is that λ s is now introduced as a second parameter to be varied.This means that an additional equation needs to be specified to uniquely solve the system.As with all pathfollowing techniques, this equation takes the form of a path-following constraint, and either Crisfield's cylindrical [76] or Riks' hyperplane constraint [47] may be used for this purpose.In general, the arc-length constraint may be written as where σ is a constant that constrains the arc-length.For the nullvector method, the general augmented system of Eq. ( 7) reduces to Eq. ( 52) features (2n + 2) equations in (2n + p) variables, and the ( p − 2) extra equations required to solve the system are implicit in the definition that the parameters Λ c are held constant, i.e.Λ c j = Σ j for j = 3 . . .p. Starting from a singular state, (u * , λ * t , Λ * c ), with associated critical eigenvector, φ, on a fundamental equilibrium path of the baseline problem, λ s = λ b s , a locus of singular points is continued using Newton's method, Approximate directional derivatives of the tangential stiffness matrix, K T , are computed using Eqs.( 36) and (37).
If the directional derivatives of the equilibrium equations, F ,λ 1 and F ,λ s , are not equal to displacement-independent mechanical load vectors, −p 1 and −p s , respectively, then approximate directional derivatives are computed using Eq. ( 38).
In the same manner, the secondary parameter, λ s , is introduced into the minimally augmented system of Eq. ( 42), and the system further augmented using an arc-length constraint.Hence, which defines (n + 2) equations in (n + p) variables with the ( p − 2) extra equations being implicit in the definition that control parameters Λ c are constants.Linearisation of Eq. ( 54) gives where approximate directional derivatives of the scalar singularity equation are computed using Eqs.( 44) and (45).
To solve the nullvector and minimally augmented systems of Eqs. ( 53) and ( 55) efficiently, a partitioning procedure is utilised such that only the tangential stiffness matrix needs to be factorised.The partitioning procedure is entirely algebraic (as shown in Appendix), and for a more detailed incremental-iterative bordering algorithm the interested reader is directed to Moghaddasie & Stanciulescu [99] and Rezaiee-Pajand & Moghaddasie [100].

Foldlines of bifurcation points
In some cases, it is useful to constrain the foldline solver purely to bifurcation points, e.g. when foldlines of limit and bifurcation points intersect and we want to prevent the solver from jumping from one curve to the other.To explicitly constrain the solver to symmetry-breaking bifurcation points, the equilibrium equations are perturbed by a vector ϕ that is antisymmetric with respect to the symmetry in the displacement vector u [44].The perturbing vector ϕ remains constant throughout the iterations of a foldline increment, and is chosen as the nullvector of the previously converged increment along the foldline.The magnitude of the asymmetry is controlled by a scalar variable, τ , and the additional equation h = u ⊤ φ = 0 enforces orthogonality of the displacement vector and nullvector throughout the equilibrium iterations.The unperturbed equilibrium equations are recovered for τ = 0, and this is input as an initial predictor to start the equilibrium iterations.In our experience, τ rarely exceeds 10 −7 throughout the iterative corrector procedure.
Note that u does not change in Eq. ( 58) because the state, (u, λ 1 , λ s ), is updated in a separate iterative procedure of the augmented equilibrium equations.These augmented equilibrium equations are which are solved in the following manner using Newton's method, The eigenvalue, µ, and eigenvector, φ, are updated at the end of every iteration of Eq. ( 60) using the iterative eigenvalue problem in Eq. ( 58).All directional derivatives in the augmented systems Eqs. ( 57) and ( 60) are calculated using the expressions referenced in the previous section, and are again solved efficiently using an algebraic partitioning procedure (see Appendix).For this purpose, the bordering algorithm presented in Refs.[44,99,100] may also be used as a template.

Tangent vectors to curves
The evaluation of the tangent space is important for predicting the direction of new solutions along a solution path.In a multi-parametric setting, the tangent space can also be used to determine the direction of greatest descent, i.e. the parameter or combination of parameters that has the greatest effect on the load-carrying capacity.The tangent space, T, of the general augmented system, G(y) = 0, as defined in Eq. ( 6), is given by the nullspace of the differential matrix, The dimension of T depends on the number of equations and the number of variables, y, in the augmented system, G(y).As outlined in Section 2.1, G generally describes a system of n equilibrium equations and r auxiliary equations in (n + q + p) variables (n degrees-of-freedom, q auxiliary variables and p parameters).In this general setting, the nullspace is (q + p − r )-dimensional.For the one-dimensional subset curves studied here, however, we define r = q + p − 1 such that the nullspace is typically one-dimensional and therefore describes the tangent direction to the subset curve.At a singular point the nullspace may be of higher dimension due to the singularity of K T within G ,y , such that the direction of the curve is not unique, e.g. two or more tangent directions at a bifurcation point.First, we address the case of non-critical subset curves (q = 0) with possibly more than one parameter ( p ≥ 1).The tangent space is now a function of y = (u, Λ) and defined by the differential matrix of Eq. ( 5): with a non-singular tangential stiffness matrix, K T , such that a single tangent vector, τ = [τ ⊤ u , τ ⊤ Λ ] ⊤ , spans the nullspace, T. Solving the first row of Eq. ( 62), where each of the p columns in T u corresponds to a displacement vector associated with a specific parameter, Λ i ∈ Λ, set to unity and all other parameters equal to zero.Hence, τ Λ = 1 i and each column i in T u describes the displacement response of the structure to an applied load vector −F ,Λ i .The relationship between the different parameters in τ Λ is established by solving the second row of Eq. ( 62) and substituting Eq. ( 63), For one-dimensional subset curves, Eq. ( 64) defines r = p − 1 equations in p parameters and thus we require one additional equation to uniquely compute the curve tangent.The simplest solution, and the one implemented here, is to define a unit value for one of the tangent components in τ Λ .This constrains the tangent vector to a unique direction but leaves its magnitude undefined up to a scalar factor.Many readers will be familiar with the partitioning procedure of classical arc-length solvers in load-displacement space (r = 0 and p = 1), whereby a solution to K T δu− p1 δλ 1 = −F is written as δu = δu 1 + δλ 1 δu 2 , with δu 2 = K −1 T p1 the tangent displacement vector (the only column of T u ) and δλ 1 assigned a value of unity.
For critical subset curves (foldlines), Eq. ( 63) has no unique solution because the tangent stiffness matrix, K T , is singular at a limit or bifurcation point.As discussed by Eriksson [61,62], the important consideration then becomes which columns i of F ,Λ are in the range of K T .Returning to Eq. ( 63), we know that the equation This condition is formally enforced by subtracting from F ,Λ i its projection in the direction of the nullvector φ.Hence, where By choosing the particular solution, τ ′ u i , to be orthogonal to the nullvector, we can combine the condition φ ⊤ τ ′ u i = 0 with Eq. ( 65) to write a system of simultaneous equations where we have defined the projection of F ,Λ i on the critical eigenvector as Π = −φ ⊤ F ,Λ i , and made use of the fact that αK T φ = 0.The system in Eq. ( 67) is invertible because the column-and rowspace of K T are expanded by the nullvector.By solving this system for each control parameter, i = 1 . . .p, we can assemble each displacement response, τ ′ u i , into the ith column of a tangent displacement matrix T ′ u such that the tangent displacement vector is, Naturally, Eq. ( 68) needs to satisfy the auxiliary equations (second row of Eq. ( 62)), Eq. ( 69) defines r = p − 1 equations in p + 1 unknowns ( p parameters and α), and therefore two more equations are needed.The first equation is the trivial case of prescribing a unit value for one of the tangent components in τ Λ .The second equation is derived from pre-multiplying the first row of Eq. ( 62) with φ ⊤ and heeding φ ⊤ K T = 0. Hence, By introducing a unit value for one of the tangent components in τ Λ , Eq. ( 70) defines a unique proportional relation between all parameter values.For most analyses, Eqs. ( 67), ( 69) and ( 70) can be solved for a unique tangent vector.In cases where g ,u φ = 0 or φ ⊤ F ,Λ = 0, to within a predefined numerical tolerance, the system of equations is singular (see Eriksson [62]), such that the tangent vector is not uniquely defined.In particular, if g ,u φ = 0, then α cannot be determined in Eq. ( 69).Similarly, when φ ⊤ F ,Λ = 0 we lose Eq. ( 70) such that there is one free variable in the system.In both these cases, we set α = 0 and solve for τ Λ using Eq. ( 69) with one parameter component assigned to unity.A non-unique tangent vector can, for example, occur at a hilltop-branching point, where a limit point and a bifurcation point coincide (two tangent directions).This is typically the case when a structure has been optimised to fail in two modes simultaneously.Due to the non-uniqueness of the tangent vector, it is more difficult to trace a locus of hilltop-branching points accurately, but an approximate method that bounds the foldline within a certain tolerance is possible (see Eriksson [61]).
When using the augmented system for bifurcation points (see Section 2.3.3),some of the non-uniqueness problems are eliminated because limit points are excluded explicitly from the tangent space.Considering the relevant augmented system in Eq. ( 57), the scalar condition of Eq. ( 70) is now transformed into which is not trivially satisfied when φ ⊤ F ,Λ = 0 due to the presence of the slack variable tangent component τ τ .
Similarly, the auxiliary equation ( 69) is augmented by the condition that where Eq. ( 68) has been used for τ u .Thus, even if g ,u φ vanishes in Eq. ( 69), α still features in auxiliary equation ( 72) and the system is non-singular.

The toggle frame: a canonical example
With the algorithms described above, a comprehensive investigation of structural stability and design parameter sensitivity can be performed.We start from an idealised model with predefined geometry, material properties and loading, and evaluate the fundamental load-displacement behaviour, including pinpointing of all relevant singular points.Additional non-critical and critical subset curves are then traced by starting from a chosen saved solution on the fundamental path.Hence, in addition to the idealised situation, the sensitivity to variations in certain modelling assumptions can be investigated, either as small variations of an imperfection sensitivity analysis or as larger variations to explore the surrounding design space.
As an example, consider the snap-through equilibrium path of a centrally loaded toggle frame, originally analysed by Eriksson [62], with clamped ends as shown in Fig. 4a, modelled using finite strain Reissner beam elements [101].Note that no units are provided for this particular example, because they are arbitrary for the purposes of this illustrative example.The toggle frame initially deforms symmetrically on the fundamental equilibrium path but this deformation mode becomes unstable at a symmetry-breaking bifurcation just before the maximum limit point on the curve.Because the connected nonsymmetric bifurcation path, branching from the bifurcation point, is unstable, the toggle frame snaps dynamically into the inverted stable shape by means of this nonsymmetric mode.As shown in Fig. 4a, the generalised path-following algorithm detects and then computes the location of all singular points exactly (to within a predefined numerical precision), and then automatically branches onto the bifurcated path to provide a complete picture of the nonlinear behaviour.
Fig. 4a restricts path-following to the classical displacement-load space.To illustrate generalised path-following capabilities, Fig. 4b extends the analysis to changes in the height, H , of the toggle frame.Fig. 4b shows an isometric view in displacement-load-height space of the fundamental and bifurcation paths discussed above, and two additional parametric paths.For these parametric paths, the applied load is held constant at P = 37.4 and P = 64.8,respectively, and the relationship between height (H ) and central displacement (w) traced using an arc-length solver.These parametric paths can be useful when the design load on the toggle frame is fixed and changes in the displacement want to be explored as a function of the toggle frame height (or any other parameter).
By imposing a singularity condition in the generalised path-following algorithm, the locus of limit and bifurcation points can be traced illustrating how changes in the height of the toggle frame affect the load-displacement solution of these singular points.The utility of these foldlines is threefold.First, they can be used to identify interesting points such as (i) the coincidence of limit and bifurcation points -the hilltop-branching points at H = 0.567 and H = 0.581 -or (ii) points where bifurcation and limit points cease to exist-the cusp catastrophes at H = 0.506 and H = 0.346.These points are clearly marked in the orthographic projections of Fig. 4c (w vs. H ) and Fig.(b) Isometric view of the fundamental and bifurcation paths in displacement-load-height space, two additional parametric paths that show the relationship between height (H ) and central displacement (w) at applied loads of P = 37.4 and P = 64.8, and the locus of limit and bifurcation points with changing height; (c) and (d) Orthographic projections of (b) in displacement-height and load-height space, respectively, that clearly denote cusp catastrophes and hilltop-branching points.(Note that no units are defined for this particular example, because they are arbitrary for the purposes of this illustrative example.Deformation mode shapes are not to scale.)(P vs. H ). Second, foldlines can be used in design studies to determine the sensitivity of singular points with respect to design parameters without having to perform computationally expensive parametric studies.Finally, foldlines can be used for optimisation purposes.For example, the displacement at the first instability load can be maximised by reducing the height of the toggle frame to coincide with the hilltop-branching point at H = 0.567 (see Fig. 4c).

Bistable composite laminates for morphing applications
Fibre-reinforced composite materials play a fundamental role in enabling shape-morphing because their orthotropy can be exploited to design structures with high stiffness in one direction, typically the loading direction, and low stiffness in another direction, a potential actuating direction.Modern technology even allows fibres to be steered in smooth curvilinear trajectories such that high-and low-stiffness regions can be connected seamlessly [6,102].In this manner, composites create the opportunity for servo-elastic tailoring.The laminated nature of advanced composites also facilitates a union of materials with dissimilar thermal expansion coefficients which can be used to build devices that exhibit large shape changes in response to variations in the surrounding temperature [103].When layers of fibre-reinforced plastic are laminated non-symmetrically, e.g.into a laminate of two layers with perpendicular fibredirections, and then cooled to room temperature after curing in an oven, differences in thermal expansion between layers cause warping into curved configurations [104].At room temperature such non-symmetric laminates are often cylindrical and can be snapped between two inverted shapes.As such, the nature of laminated composites allows residual stresses to be coupled conveniently with the geometric nonlinearity of curved shells to design bistable [105][106][107][108][109] or even tristable shell structures [29,38].
The main aim herein is to show that generalised path-following is particularly well-suited to modelling and designing morphing structures of arbitrary geometry, loading and constitutive properties.This is illustrated here by means of a classical example problem from the literature (see Refs. [104,33,110]).A flat nonsymmetric cross-ply laminate that is cooled-down from cure to room temperature initially deforms into a saddle-like shape but very soon bifurcates into one of two cylindrical solutions.At room temperature, the ensuing cylindrical panel can be snapped from one cylindrical shape into its inverted counterpart, thereby creating a bistable plate that is useful for morphing applications.The associated snap-through behaviour is intricate with multiple instability points, and its full complexity is exposed here for the first time (Section 4.1).Finally, the phase diagram of bistability can be altered significantly by introducing initial curvature during curing (Section 4.2).

The initially flat laminate
Consider a four-layer [90 2 /0 2 ] square carbon/epoxy plate with material properties shown in Table 1, generic inplane dimensions L = L x = L y and fully clamped at its planar midpoint.As a baseline model, we assume panel dimensions L x = L y = 0.25 m, although the effect of varying these dimensions is investigated later.The composite laminate is discretised into a 43 × 43 node mesh of 196 fully integrated 16-node, total Lagrangian shell elements using the shell director parametrisation of Ramm [111] (two rotations per node, no drilling around director).The effects of shear and membrane locking is minimised by using cubic isoparametric interpolation and by refining the mesh sufficiently until convergence with respect to the results by Pirrera et al. [33] is obtained, who used a 51 × 51 node mesh of 2500 reduced integration 4-node S4R elements in the commercial finite element software ABAQUS.
As a first load step, the post-cure cool-down of the originally flat [90 2 /0 2 ] laminate from 200 • C to room temperature of 20 • C is simulated.Fig. 5a shows the behaviour of the laminate in terms of the absolute value of the change in temperature |dT | = [0, 180] K and the out-of-plane displacement w of one of the corners of the laminate.On the fundamental equilibrium path the laminate cools into a saddle shape with zero transverse deflection of the corners.However, this fundamental saddle shape soon becomes structurally unstable at a temperature change of dT = −4.15K.The eigenmode associated with this bifurcation point is cylindrical, meaning that, as cooling proceeds, the laminate transitions in a structurally stable manner from the initial saddle shape into a semi-cylinder at room temperature.By nature of being a symmetric pitchfork bifurcation, this cylindrical shape can take one of two shapes (shown in Fig. 5a), curving up in one direction or curving down in the orthogonal direction.Fig. 5a shows that there are further bifurcation points on the fundamental saddle-shape path, but because the branches emanating from these points are all structurally unstable, they are not explored in more detail herein.However, as is shown later, these additional bifurcation points are important for the snap-through behaviour.
Upon cooling to room temperature (dT = −180 K), snap-through from one cylindrical shape (up or down) to the other (down or up, respectively, with the curvature rotated through 90 • ) is initiated by applying a transverse point load at each of the four corners of the laminate.In the generalised path-following algorithm this loading is readily incorporated as an additional parameter such that the load-displacement plot of snap-through can be superimposed on the previously defined pitchfork diagram.Fig. 5b shows this plot in three dimensions, with the transverse force F at one corner plotted versus the associated transverse displacement w at room temperature dT = −180 K.The load-displacement diagram shows the typical characteristics of a bistable structure with three displacement-axis intercepts, two stable and one unstable, and the necessary limiting maxima where the structure loses stability and snaps from one stable shape to the other.Inset A of Fig. 5b confirms the computational findings of Pirrera et al. [33], and the experimental findings by Potter et al. [107], that snap-through is in this case characterised by a multi-snap event-the equilibrium path reaches a first limit point (L1) causing it to snap to the adjacent stable region, at which point the load increases slightly before the structure loses stability entirely and snaps into the inverted and rotated shape.Contrary to the computational findings by Pirrera et al. [33], this loss of stability does not occur at the second maximum (L3) but at a preceding bifurcation point (B1).The semi-analytical Ritz model by Pirrera et al. used the continuation solver EPCONT but bifurcation point (B1) was not recognised by the solver due to symmetry conditions imposed on computational efficiency grounds.The ABAQUS model used by the same authors simply did not possess the capabilities to detect this bifurcation point.In fact, inset A of Fig. 5b shows two further bifurcation points (B2) and (B3) that are here determined for the first time.Finally, inset B shows a fourth limit point (L4) in the region of the unstable self-equilibrated state of zero applied corner force and zero corner displacement.Due to the symmetry between the two cylindrical deformation modes, it is no surprise that the snap-through equilibrium path is symmetric about the origin.Therefore, the four limit points and three bifurcation points have corresponding inverted points, which in the following are identified by the letter "i", e.g.(L1i), (L2i), (B1i), etc., and referred to as "inverted analogues".The branches emanating from bifurcation points (B1), (B2), (B3) and their inverted analogues are now explored in detail.Fig. 5c shows the bifurcation branch emanating from (B1), which connects to the corresponding bifurcation point (B1i) of the inverted and rotated cylindrical shape.This bifurcation path is unstable throughout, apart from a very localised region close to, but disconnected from, (B1) and (B1i).Furthermore, the regions of the curve close to (B1) and (B1i) are tangled and punctured by multiple other limit and bifurcation points.It would be possible to explore the secondary branches emanating from these bifurcation points, but as most of these are likely to be unstable and only add further complexity to the equilibrium diagrams, they are not continued further in this work.Fig. 5d presents the path branching from bifurcation point (B2).This branch does not connect to its inverted analogue (B2i) but to (B3i) instead.An identical path mirrored about the origin (symmetric about both axes) also exists which connects (B2i) with (B3), thereby maintaining the expected symmetry of the bifurcation points.These bifurcation branches are again quite tangled close to points (B2) and (B2i), and multiple secondary bifurcation points exist that are not explored further herein.However, there exist six localised regions that are structurally stable.The mode shapes for the two centrally located regions on either side of w = 0 are shown and these are rotationally out-of-phase by 180 • .The six localised regions of stability are surrounded by unstable regions and can, in an experimental setting, not be reached by simply increasing the applied loading.Instead, a force of 1 N could be applied to all four edges and the laminate then forced manually into the depicted mode shape by hand or ideally using a mould.Forcing the structure into this shape requires traversal of an energy hump, but once the required energy threshold has been passed, the laminate should naturally snap into the envisioned shapes.
Given the unstable nature of the depicted bifurcation paths, the reader might be led to believe that information regarding these paths is of no use to the practising engineer.However, it is well known that such subcritical bifurcations can lead to localisation phenomena [72] and/or extreme sensitivity to initial imperfections [52] that can detrimentally effect the load carrying capability of the structure.Furthermore, the unstable regions can play a crucial role in dynamic behaviour.Hence, an awareness of the existence of unstable post-buckling branches is necessary if such bistable laminates are to be used safely and reliably in engineering design.
A pertinent question to ask at this point is how the multi-stability of the laminate changes with temperature.Hence, if the laminate only cools to 50 • C rather than to room temperature, then how is the snap-through behaviour affected?As the snap-through behaviour is governed by limit and bifurcation points on the snap-through equilibrium path, we can use the limit point and bifurcation point continuation capability to determine how each of the previously identified points (L1)-(L4), (B1)-(B3) and their inverted analogues change with temperature.Indeed, Fig. 6a reveals that many of these points are directly related to bifurcation points on the fundamental cool-down path of Fig. 5a.Individual limit point and bifurcation point loci are shown in more detail in Fig. 6b-6d.Fig. 6b shows that limit point (L1) and its inverse analogue (L1i) arise as a direct continuation of the first bifurcation point on the fundamental cool-down path.This is not surprising as without cooling beyond this first critical temperature, the two cylindrical shapes do not exist, and hence there is no snap-through between them.It is also evident that the transverse corner forces at limit loads (L1) and (L1i) increase monotonically in magnitude with increased post-cure cooling, meaning that the greater the degree of post-cure cooling, the greater the required force to snap between the two cylindrical shapes.Furthermore, limit points (L1) and (L1i) remain the first singular points throughout the entire temperature range and are therefore expected to initiate snap-through at all times.
Fig. 6b-6c show that all three bifurcation points (B1)-(B3) and their inverted analogues (B1i)-(B3i) arise as a result of instabilities on the fundamental cool-down path.What is interesting is that bifurcation points (B2) and (B3) arise concurrently at a double bifurcation point as shown in the inset of Fig. 6c.On a snap-through diagram (F vs. w) bifurcation points (B2)/(B3i) and (B3)/(B2i) first appear at the two cusp catastrophes dT = −17.0K as four unique points.With further cooling to dT = −17.9K two of these bifurcation points, (B3) and (B3i), then coincide at a single point.Note that the connection of (B2) to (B3i) and (B3) to (B2i), depicted in the loci of bifurcation points of Fig. 6c, mirrors their connection by the bifurcation branches on the snap-through diagram of Fig. 5d.Finally, bifurcation points (B1) and (B1i) arise on the fundamental cool-down path at dT = −19.8K such that beyond this temperature, the snap-through equilibrium path is governed by one pair of limit points, (L1) and (L1i), and three pairs of bifurcation points, (B1)/(B1i)-(B3)/(B3i).On the other hand, limit points (L2)/(L2i) and (L3)/(L3i) do not arise because of an instability on the fundamental cool-down path.Fig. 6d shows that limit points (L2) and (L3) are connected, and by symmetry, so must be their inverted analogues (L2i) and (L3i).As the extent of post-cure cool-down is decreased, these two pairs of limit points move closer to each other, finally colliding and vanishing at the cusp catastrophe at dT = −155.1 K.This means that the multi-snap event described above does not exist for cool-down temperatures > −155.1 K. Fig. 6d shows another cusp catastrophe, this time involving limit points (L4) and (L4i) which collide at the final bifurcation point on the cool-down path dT = −167.7 K. Consequently, the localised feature of the snap-through equilibrium path defined by limit points (L4) and (L4i) is straightened out for dT > −167.7 K. Finally, we would like to investigate the snap-through behaviour evolves with the side length, L, of the laminate.The baseline design considered above a length of L = L x = L y = 0.25 m and this is now systematically varied in the generalised path-following algorithm using two possible approaches.From the cooldown pitchfork diagram (Fig. 5a) we know that snap-through can only occur beyond the first instability point on the fundamental cool-down path, i.e. when the panel transitions from the saddle-like shape to one of two cylindrical shapes.Hence, we use the generalised path-following algorithm to trace the locus of this first pitchfork bifurcation point with decreasing side length.Fig. 7a shows this locus as a three-dimensional view in L vs. dT vs. w space, whereas Fig. 7b shows the same plot as an orthographic projection in L vs. dT space (w = 0 along the entire locus).It is clear that the shorter the side length, the greater the degree of cool-down to initiate bistability.In fact, the relationship between side length and the first critical temperature is exponential.Initially, decreasing the side length only requires a marginal increase in cooling before the saddle-like shape becomes unstable, but for side lengths smaller than 0.1 m, the degree of cooling required increases rapidly and approaches infinity as the side length tends to zero.In a similar manner, we could pose the question of how the two cylindrical shapes at room temperature (dT = −180 K) evolve with decreasing side length.The V-shaped response in Fig. 7a shows that with decreasing side length the two cylindrical shapes become shallower until these shapes vanish completely for a side length of 0.038 m, which incidentally, is the point where this path intersects the locus of bifurcation points.Hence, a representative panel with side length of 0.035 m is not bistable at room temperature.

The initially curved laminate
Generally speaking, morphing devices require different snap-through responses and geometries for different applications.First, curing the laminate on a curved tool plate allows the geometry to be adapted to different scenarios.Second, by forcing the stress-free curing configuration to be curved the post-cure cooling phase does not lead to two equal and opposite cylindrical shapes as is the case for a flat laminate.As a result, the snap-through and snap-back response at ambient temperature can be adjusted to serve specific design requirements.Although the layup of the composite may also be changed to influence these characteristics, this approach changes the thickness and stiffness of the laminate in a manner that may not be beneficial or intended in terms of overall component mass and structural stiffness.
Previous papers studying the effect of initial curvature on the bistability of laminates, e.g.[34,112,113], have studied the effect of curvature in a classic parametric manner.Here we show that the sensitivity of important features -such as initiation of bistability -to initial curvature can be established more efficiently using generalised pathfollowing.
Consider a six-layer [0 3 /90 3 ] square carbon/epoxy laminate with in-plane dimensions L x = L y = 0.23 m and fully clamped at its planar midpoint.The laminate is cured on a curved tool plate with curvature of κ y = −2.5 m −1 in the 90 • fibre direction and κ x = 0 m −1 in the 0 • fibre direction.This problem is representative of the analytical and experimental study presented by Eckstein et al. [34] only that in the present study, following Lamacchia et al. [28], the temperature-independent material properties shown in Table 2 are assumed.The same mesh and elements as described in the previous section are used to model the structure.
The initially curved plate is cooled through a temperature differential of dT = −180 • C from a stress-free curing temperature of T = 180 • C. Because the original unstressed shape of the laminate is not flat the laminate does not initially deform into a saddle shape upon cooling (see Section 4.1).Rather, as shown in Fig. 8a, the laminate begins to flatten out its original cylindrical shape until this shape becomes unstable at a symmetry-breaking pitchfork bifurcation.At this point (dT = −57 K) laminate bifurcates into one of two twisted mode shapes which is a transition mode that transforms the negative y-direction curvature at curing temperature into negative x-direction curvature at temperature.Hence, the second bifurcation point (dT = −126 K) connects the twisted mode with the second cylindrical mode which is rotated by 90 • with respect to the original curing curvature.The two temperatures T = 123 • C and T = 54 • C match with the findings of Lamacchia et al. [28] By tracing a foldline of these bifurcation points we can establish the range of initial curvatures for which the twisted transition mode exists.Fig. 8b shows that the two bifurcation points that mark the beginning and end of the twisting transition mode move closer to each other and then collide at a cusp.For initial curvatures κ y > −0.67 m −1 , the twisted transition mode ceases to exist, and for increasing negative curvature κ y < −0.67 m −1 the temperature differential between the two bifurcation points increases.
Fig. 8c shows the cooling behaviour of the flat laminate which, as discussed in the previous section, initially deforms in a saddle shape and then transitions into one of two cylindrical modes.For the particular material properties and geometric dimensions assumed here, this occurs for dT = −50 K and no higher-order bifurcations are present in the temperature range analysed here.Fig. 8c also shows how the bifurcation point that initiates the saddle-to-cylinder transition unfolds with variations in the initial curvature of the laminate.It is evident that this foldline does not intersect the configuration κ y = −2.5 m −1 in the temperature range considered here, and this is why no saddle-to-cylinder transition exists for the curved laminate in Fig. 8a.Finally, this foldline is not symmetric with respect to the line κ y = 0 m −1 , i.e. the flat laminate, as the non-symmetric stacking sequence of the laminate means that negative and positive curvatures |κ y | lead to different post-cure internal stresses.
The introduction of initial curvature breaks a geometric plane of symmetry such that the perfect pitchfork bifurcation of the flat laminate is transformed into a broken pitchfork.Fig. 8d shows one of these broken pitchforks for an initially curved laminate of κ y = −0.25 m −1 , alongside the perfect pitchfork of the flat laminate.As expected, unfolding the saddle-to-cylinder bifurcation of the flat laminate with respect to κ y (also shown in Fig. 8c) intersects the broken away equilibrium path at a limit point.In general, Fig. 8 can be used intuitively to determine the combination of post-cure cooling and initial curvature required to render a bistable structure.

Conclusions
Generalised path-following combines a numerical continuation algorithm with the geometrical versatility of the finite element method.The advantage of this technique over most path-following algorithms in commercial finite element codes that trace curves in load-displacement space is that any parameter influencing the structure, be it load, geometrical dimensions or constitutive properties, etc., can be followed on a multi-dimensional solution manifold.This means that the sensitivity of critical points with respect to design parameters is efficiently established by bifurcation tracking in parameter space.Hence, computationally expensive parametric studies, that trace the full nonlinear equilibrium diagrams for many iterations of a baseline design, are replaced by the more efficient process of tracing the locus of specific design points with respect to potentially any number of varying parameters.
Such a capability is especially useful for a new class of well-behaved nonlinear structures that exploit buckling instabilities for novel functionality and shape adaptation.This is because tracking of the design-driving snap-through instabilities and bifurcations enables rapid exploration of the design space.The present paper combines, adapts and extends concepts introduced by Eriksson [40] and Salinger et al. [44] to apply them to the design of thermo-mechanical multi-stable structures that are especially promising for the aerospace industry.The techniques outlined herein are efficiently implemented in existing research codes or programming languages that have access to ARPACK [94], and can even be programmed to function as a wrapper around existing codes [66].
We have shown that a generalised path-following algorithm is ideally suited to the analysis and design of bistable composite laminates for shape-morphing applications.The results herein show that the full complexity of multisnap events is captured robustly by this algorithm, and that the ability to follow loci of singular points with respect to additional parameters can yield invaluable insights into the underlying physics of multi-stable phenomena.The bifurcation-tracking capability has been especially valuable to trace the boundary of bistable regions in parameter space.Finally, the ability to branch-switch to additional bifurcation branches has revealed a number of localised regions of stability, which could be exploited as intermediate shapes in a multi-event shape-adaptation process.Future work will focus on exploiting these localised islands of stability for potential multi-shape adaptation.

Data Statement
All data required to reproduce the figures in this paper can be found on the data repository of the University of Bristol via the URL: https://data.bris.ac.uk/data/.

Appendix. Path-following
Equilibrium curves are evaluated using typical predictor-corrector schemes of the constrained Newton's method [47,76].The equilibrium curve is "path-followed" by incrementally evaluating different solution points along the curve.Fig. A.1 illustrates this concept schematically for a curve parametrised by two parameters, λ 1 and λ s .For such a curve, the predictor is evaluated using the tangent space formulated in Section 2.3.4.The corrector is computed in an efficient manner by algebraic partitioning-also known as a block elimination bordering algorithm [44,79].
For the general system defined in Eq. ( 7) with n displacement degrees-of-freedom, u, p parameters, Λ, and q auxiliary variables, v, the corrector, δy = (δu, δΛ, δv), is computed from The r = p + q − 1 auxiliary equations, g, can generally be split into q equations, g v , that pertain to the auxiliary variables, v, and p−1 constraining equations, g c , that either constrain the magnitude of v and/or some of the parameters in Λ.For example, when following foldlines in two parameters (see Section 2.3.2), the auxiliary variable, v, is the eigenvector of the tangential stiffness matrix, F ,u ≡ K T , with g v = K T v = 0 and g c = ∥v∥ 2 = 1.To reflect this partitioning of g, we rewrite Eq. (A.1) as g c,u g c,Λ g c,v N ⊤ ,u N ⊤ ,Λ 0 1×q where Eq. (A.3) has been used to replace δu, and both δv 1 and δv 2 are computable quantities using the inverse of g v,v , which has dimensions q × q.In many cases, such as the foldline example given above, g v,v ≡ K T such that solving Eqs.(A.3) and (A.4) relies on the inversion of the same matrix K T .Finally, the parameter corrections, δΛ, are determined from the third and fourth rows of Eq. (A.where Eqs.(A.3) and (A.4) have been used for δu and δv, respectively, and the system in Eq. (A.5) is solvable as long as the preceding p × p matrix is invertible.In most cases, the system in Eq. (A.5) does not need to be inverted numerically but can be solved algebraically a priori. 4d

Fig. 4 .
Fig. 4. (a) Fundamental and bifurcation equilibrium paths of load (P) versus central displacement (w) for a toggle frame of height H = 0.65;(b) Isometric view of the fundamental and bifurcation paths in displacement-load-height space, two additional parametric paths that show the relationship between height (H ) and central displacement (w) at applied loads of P = 37.4 and P = 64.8, and the locus of limit and bifurcation points with changing height; (c) and (d) Orthographic projections of (b) in displacement-height and load-height space, respectively, that clearly denote cusp catastrophes and hilltop-branching points.(Note that no units are defined for this particular example, because they are arbitrary for the purposes of this illustrative example.Deformation mode shapes are not to scale.)

Fig. 5 .
Fig. 5. (a) Pitchfork bifurcation diagram of an unsymmetric [90 2 /0 2 ] laminate.Upon cool-down from cure temperature at 200• C the laminate initially deforms into a saddle shape, but this deformation mode is structurally unstable beyond a change in temperature of −4.15 K (all mode shapes not drawn to scale).Beyond this temperature, a symmetry-breaking cylindrical shape is added as a perturbing eigenmode such that at room temperature of 20 • C, the laminate takes one of two cylindrical mode shapes; (b) Snap-through diagram from one cylindrical shape to the inverted and rotated cylindrical shape at room temperature.As the inset shows, the snap-through event is not governed by a simple limit point.Rather, the panel loses stability at (L1) and then finally snaps to the inverted and rotated shape at bifurcation point (B1); (c) Snap-through equilibrium path with additional bifurcation branch connecting bifurcation point (B1) and its symmetrically inverted analogue (B1i); (d) Snap-through equilibrium path with additional bifurcation branch from (B2) to (B3i).The bifurcation branch includes six localised stable regions with two of the associated deformation modes shown.

Fig. 7 .
Fig. 7. (a) Isometric view of the locus of first instability points on the fundamental cool-down path with changing laminate side length.The plot also shows that the two cylindrical shapes become shallower with decreasing side length; (b) An orthographic view of figure (a) in side length versus temperature change space.As the side length decreases, there is an exponential increase in cooling required to induce a bistable laminate.

Fig. 8 .
Fig. 8. (a) Bifurcation diagram of an unsymmetric [0 3 /90 3 ] laminate with initial cylindrical curvature κ y = −2.5 m −1 .Upon cooling from cure temperature the laminate flattens out until the initial cylindrical shape becomes structurally unstable and transitions into one of two twisting modes.These twisting modes connect to a second rotated cylindrical shape at ambient temperature (all mode shapes not drawn to scale); (b) Foldline describing the variation of the two bifurcation points in (a) with respect to κ y .Note that the twisting mode disappears for an initial curvature of κ y > −0.67 m −1 ; (c) Orthographic projection of (b) in temperature-curvature space.The cool-down equilibrium path and unfolding of the saddleto-cylinder bifurcation of the flat laminate are also shown; (d) Broken and unbroken pitchforks of an initially curved laminate (κ y = −0.25 m −1 ) and an initially flat laminate cooled from curing to ambient temperature.

F 2 )Fig. A. 1 .
Fig. A.1.Schematic illustrating path-following along a solution curve using a typical predictor-corrector scheme.The solver starts from a known equilibrium solution A and iteratively converges to a new solution B along the curve.This procedure is then repeated incrementally to progress along the curve.The figure shows a spherical arc-length constraint but planar, cylindrical and elliptical constraints are equally possible.

Table 1
Typical material properties for a carbon/epoxy composite.