Stability constraints for geometrically nonlinear topology optimization

This paper compares four methods for formulating stability constraints in topology optimization with geometric nonlinearity. The methods are: a direct approach to compute the critical load factor, an approximation using an eigenvalue analysis at a load factor of 1, a new method based on an eigenvalue analysis at the constraint limit load factor, and an implicit method based on stiffness reduction, which has not previously been investigated for stability constraint formulation. These four methods are described in detail and then compared qualitatively and quantitatively (including optimization examples) in terms of accuracy, robustness, and computational efficiency. The results show that formulating the constraint using an eigenvalue analysis at a load factor of 1 is the most robust approach, as it is least likely to experience mode switching or mode skipping during optimization, which leads to poor convergence for the other three methods. It is also the most efficient, as it only requires a single eigenvalue solve, whereas other methods require additional linear solves to compute the constraint value. However, an eigenvalue analysis at a load factor of 1 only approximates the critical load factor, which may be over, or under-estimated. Therefore, none of the methods fully satisfy the criteria of accuracy, robustness, and efficiency, highlighting the need for further research, e.g., by improving the accuracy of the method based on an eigenvalue analysis at a load factor of 1, or by improving the robustness and efficiency of the direct approach.


Introduction
Stability constraints play an important role in structural optimization problems to ensure that optimized designs can safely support the intended loads without undergoing large displacements, or even collapse, due to instabilities such as buckling.Topology optimization problems that minimize linear compliance (maximize linear stiffness) often result in structures with thin members under compression that are susceptible to failure by buckling instability.This is a wellknown issue and has led to research on including stability constraints in topology optimization problems.This is usually achieved by estimating the critical load factor by a linear buckling analysis (Neves et al. 1995;Dunning et al. 2016;Ferrari and Sigmund 2019).However, buckling is inherently a geometrically nonlinear problem and the linear approximation may over or under-estimate the true load capacity of the structure.
Therefore, a geometrically nonlinear analysis can be used to compute the critical load factor more accurately (de Borst et al. 2012;Pedersen and Pedersen 2018;Dalklint et al. 2021).However, this approach has received less attention in topology optimization, possibly due to the high computational cost of the nonlinear analysis.Another issue is that accurately identifying the limit point for the critical load level is difficult because the tangent stiffness matrix becomes singular at a limit point.Thus, in practice, we can only obtain a solution near a limit point (de Borst et al. 2012;Kemmler et al. 2005).This leads to inaccuracies in the computed gradients, which can cause convergence issues when using a gradient-based optimizer, especially if the geometrically nonlinear stability constraint is difficult to satisfy (i.e., the optimal design with a stability constraint is significantly different to one without).
Several approaches have been proposed for formulating geometrically nonlinear stability constraints in the context of topology optimization.Kemmler et al. (2005) used an extended system to directly compute the critical load factor, and proposed a method for including geometric imperfection during optimization.Lindgaard and Dahl (2013) used a path tracing method, estimating upcoming unstable points using an eigenvalue analysis.However, this approach was only used to maximize the critical load factor as an objective, and was not used as a constraint.More recently, several studies approximate the critical load factor using an eigenvalue analysis at a nonlinear equilibrium point (Pedersen and Pedersen 2018;Dalklint et al. 2021;Zhang et al. 2023).Furthermore, Wang and Sigmund (2020) used an implicit method based on the reduction in stiffness to estimate the critical load factor of optimized infill structures, although this was only applied during post-optimization analysis, and not considered during optimization.Each of these approaches has benefits and drawbacks, but there is no study in the literature that directly compares these methods in the context of topology optimization.
Therefore, the aim of this paper is to compare different methods for formulating stability constraints in geometrically nonlinear topology optimization, focusing on their computational efficiency, numerical robustness, and accuracy in estimating the critical load factor.Four approaches are considered, including direct, eigenvalue, and implicit methods.These are detailed and qualitatively compared in section 2. A quantitative comparison is then performed in section 5 using the co-rotational and arc-length numerical methods for geometric nonlinear analysis (see section 3).Finally, optimization examples are presented in section 6, using a density-based topology optimization method, detailed in section 4.

Stability constraint formulations
Equilibrium points for a geometrically nonlinear analysis with displacement independent loading are found iteratively by solving the following force residual equation: where r is the force residual vector, u is the displacement vector, f int (u) are the internal forces, f ext are the nominal external forces, and is the load factor (  > 0).
In a geometrically nonlinear analysis, an unstable equilibrium point is identified when the tangent stiffness matrix, K has at least one zero eigenvalue (de Borst et al. 2012;Kemmler et al. 2005): (1) where u * is the displacement vector, corresponding to a load factor * , and is an eigenvector.The critical load factor is then the lowest value of * : γ = min( * ).
In general, there are two types of instability: limit points and bifurcations (de Borst et al. 2012;Lindgaard and Dahl 2013).A limit point is a local maximum on the load-displacement equilibrium path, and a bifurcation point is when there are two (or more) equilibrium paths for an increment in the load (also known as loss of uniqueness), characterized by: T f ext = 0 .When a limit point is passed, one eigenvalue of the tangent stiffness matrix becomes negative, representing a decrease in load carrying capacity, for an increase in displacement (or negative gradient on the load-displacement path).However, if load capacity starts to increase, then the minimum eigenvalue becomes positive again.One, or more eigenvalues also become negative when a bifurcation point is passed, which are related to the other (non-trivial) equilibrium path, or paths.However, if the iterative solution finds an equilibrium point on the pre-buckling path, the minimum eigenvalue switches back to positive.
The aim of a stability constraint in geometrically nonlinear topology optimization is to ensure the critical load factor, γ , is greater than some limit value: γ ≥  lim > 1 , where we assume the nominal external forces, f ext , are scaled to the magnitude of the design load.The challenge is to formulate this constraint (or an equivalent) such that it is computationally robust, efficient, and accurate.Four possible constraint formulations are presented below and analyzed qualitatively, comparing their robustness and efficiency during optimization, and their accuracy in computing the critical load factor.Three of the formulations are taken from the literature, and one is proposed in this paper.

Direct method
The first approach is to directly compute the critical load factor by satisfying equations ( 1) and (2), as used by Kemmler et al. (2005).The constraint is then: Iterative schemes are used to directly compute the critical load factor.Although, in practice it is only possible to find a point close to the critical one, due to the tangent stiffness becoming singular at an unstable point.Kemmler et al. (2005) used an extended system to directly compute the critical point, where equations (1) and ( 2) are solved simultaneously, with an additional equation constraining the length of the eigenvector.A penalty function was added to avoid issues with a singular tangent stiffness matrix, and thus the approach converges on a point in the neighborhood of the critical point.
(3)  lim − γ ≤ 0 Alternatively, a path tracing method, such as the arclength method, can be used to incrementally compute several points on the equilibrium path.At each converged point, the eigenvalues of K can be checked to see if they are near zero, or negative (Lindgaard and Dahl 2013).Further iterations can be used to refine the estimate of the critical point.Whichever iterative scheme is used, there is significant computational cost in directly computing the critical load factor.Another issue is mode skipping, where the iterative scheme skips past the critical load factor due to large steps, as illustrated in Fig. 1.Thus, small steps may be required to make this approach robust, further increasing computational cost.Even if the critical load factor is computed accurately, the well-known problem of mode switching may occur, where the critical mode switches between optimization iterations (Lindgaard and Dahl , 2013) [as also observed in linear buckling problems (Bruyneel et al. 2008;Dunning et al. 2016)], which slows down, or even prevents convergence.Kemmler et al. (2005) suggested using an eigenvalue analysis at the equilibrium point to approximate higher modes, although this was not implemented in their work.
Derivatives with respect to a structural design parameter x are (Kemmler et al. 2005): where is an eigenvector of the tangent stiffness matrix, such that: K(u * ) = 0 .However, this condition is only approximately satisfied by the iterative procedure to directly locate the critical load factor.Furthermore, if γ is a bifurca- tion point then: T f ext = 0 , and an alternative method for computing the derivative is required: where p is the solution of the adjoint equation: In Eq. 6 the derivative of the tangent stiffness matrix with respect to displacement in the direction of the eigenvector can be efficiently computed using a directional derivative finite difference approach (Pajot and Maute 2006;Kemmler et al. 2005): where is a small number.Note that in Eq. ( 7) is scaled such that: ‖ ‖ 2 = ‖u * ‖ 2 , and: sign( T f ext ) = sign(u * T f ext ) , to ensure numerical robustness.Also, note that if there are multiple zero eigenvalues of the tangent stiffness matrix, then the eigenvector used in the derivatives, Eqs. ( 4) to (7), will be non unique.This is the well known problem of multiple, or coinciding modes (Seyranian et al. 1994).However, based on the authors experience, these are not often seen, and does not occur for any of the examples shown in this paper.Thus, no special treatment of coinciding modes is used.

Eigenvalue analysis at = 1
Directly computing the critical load factor can be computationally expensive if a large number of iterations are required.An alternative is to estimate the critical load factor using an eigenvalue analysis at = 1 (Pedersen and Pedersen  2018; Dalklint et al. 2021; Zhang et al. 2023).This relies on a partition of the tangent stiffness matrix: where K 0 is the initial tangent stiffness matrix (for u = 0 ), K u is the part dependent on displacement, and K s is the part dependent on stress.We then assume that part of the tangent stiffness matrix varies linearly with an increase in load factor, while the other parts remain constant, leading to an eigenvalue problem.For example, if the tangent stiffness matrix is computed at = 1 , and only K s is assumed to vary with an increase in load factor, we get the following eigenvalue problem (Lindgaard and Dahl 2013): where u 1 is the displacement vector at = 1 .Thus, the esti- mated critical load factor is: γ ≈  min and the constraint is then: (5) Note that some authors assume both K u and K s vary linearly with the load factor (Dalklint et al. 2021;Kemmler et al. 2005).
This method is similar to using a linear buckling analysis, except additional information is included via the tangent stiffness matrix computed at = 1 .However, the well-known issues of mode-switching (Bruyneel et al. 2008;Dunning et al. 2016) and coinciding eigenvalues (Seyranian et al. 1994) can still occur leading to convergence problems.Therefore, it is common to include several of the lowest modes in the problem.Here, the Kreisselmeier-Steinhauser constraint aggregation function is employed to estimate min in Eq. ( 10) using the lowest N eigenvalues, and an auxiliary variable: = 1∕ where is a parameter that controls the sharpness of the approximation, and λmin provides a lower bound on min .
The derivative of a single eigenvalue with respect to a structural design variable x is (Kemmler et al. 2005): which requires the solution of the following adjoint equation: where, again, a directional derivative finite difference approach is used to compute the right-hand side of the adjoint, see Eq. (7).

Eigenvalue analysis at = lim
The eigenvalue approximation of the critical load factor presented in the previous section may still over, or underpredict the true critical load factor.However, at any equilibrium point, the value of the lowest eigenvalue from Eq. ( 9) can be used to determine if an unstable point has been passed (  < 1 ), or not (  > 1 ) (Lindgaard and Dahl 2013).This feature is used to propose a new stability constraint formulation for geometrically nonlinear topology optimization, based on an eigenvalue analysis at = lim .If the lowest eigenvalue is greater than 1, then the structure is stable at the required load factor: Note that in Eq. ( 14) min is computed using Eq. ( 9), except the tangent stiffness matrix is evaluated at equilibrium for = lim , and the Kreisselmeier-Steinhauser function, Eq. ( 11), is again used to approximate min .
This formulation may provide a more accurate indicator of stability, compared with an eigenvalue analysis at = 1 , but mode skipping may occur.For example, if a limit point has been passed at  <  lim , but the gradient of the load-displacement curve is positive at = lim , then all eigenvalues will be positive, predicting the next critical point, and not the one that has already been passed.A similar issue can occur with bifurcation points, depending on which branch is followed after the bifurcation.Furthermore, additional computational cost is required to reach equilibrium at = lim , compared with = 1.

Stiffness reduction
Another approach is to implicitly indicate instability by monitoring the reduction in tangent stiffness along the equilibrium path.This approach was used by Wang and Sigmund (2020) when evaluating the geometrically nonlinear behavior of infill structures, although the author is unaware of it being used to formulate stability constraints for topology optimization.
The reduction in stiffness factor, , for a converged equilibrium point u , is defined relative to the initial stiff- ness under the applied loading: where u 0 and u are the solutions to: Therefore,  < 1 indicates a reduction in stiffness, compared to the initial stiffness at u = 0.
The constraint is formulated by defining an allowable reduction in stiffness ( 0 <  < 1 ).An iterative scheme is then required to solve Eqs. ( 1), ( 15) and ( 16) simultaneously for the load factor ( ) and displacement ( u ) that give the required stiffness reduction.This load factor is then used to define the constraint: The iterative scheme used to find in discussed in section 3.2.The derivative of with respect to a design variable x is: where p is an adjoint vector: Note that if = 0 , then this formulation is similar to the direct approach   with Eqs. (3, 5 and  6)), which is expected, as = 0 indicates zero stiffness and a limit point.However, they are not the same, as the stiffness reduction method cannot detect bifurcation points because eigenvalues, Eq. ( 2), are not evaluated.
This approach avoids the issue that the direct approach can only approximately locate the limit point.However, limit point mode skipping can still occur (see Fig. 1) and bifurcation points are not identified, as no eigenvalue analysis is used.If mode skipping or bifurcations do not occur, then the accuracy of the computed critical load factor depends on the value of used, with lower values providing a more accurate estimate.

Qualitative comparison summary
The discussion in this section on the accuracy, efficiency, and robustness of the four proposed geometrically nonlinear stability constraint formulations is summarized in Table 1.This shows that the direct approach and the eigenvalue analysis at = lim are the most accurate methods, if mode skip- ping is avoided.However, to avoid mode skipping, smaller steps in the iterative scheme are required, thus increasing computational cost.The method using an eigenvalue analysis at = 1 appears the most computationally efficient, as it only requires one additional eigenvalue solve, whereas existing implementations of the direct method require several (18) additional linear and eigenvalue solves, as does the formulation using an eigenvalue analysis at = lim .However, solv- ing large scale eigenvalue problems can still be computationally expensive.Using an eigenvalue analysis at = 1 is also the most robust, but may not be accurate, as the critical load is still approximated.The stiffness reduction formulation also gives an approximate value for the true critical load factor, but with potentially significantly more computational cost in simultaneously solving Eqs.(1, 15 and 16).It also cannot identify bifurcation points, and still may suffer from mode skipping and switching.These points are explored quantitatively using numerical examples in sections 5 and 6, using the geometrically nonlinear modeling and topology optimization method detailed in the next two sections.
3 Geometric nonlinear modeling

Co-rotational method
Most geometrically nonlinear topology optimization methods employ nonlinear strain measures to model the nonlinearity, such as Green's strain, e.g., Buhl et al. (2000); Kemmler et al. (2005); Lindgaard and Dahl (2013); Pedersen and Pedersen (2018), or hyperelastic material models, e.g., Wang et al. (2014); Luo et al. (2015); Dalklint et al. (2021).However, these usually require additional techniques to avoid convergence issues caused by highly distorted void elements.An alternative is the co-rotational formulation, where large displacements are modeled by attaching a local coordinate system to each element.This allows an element's rigid body rotation to be computed and removed from the local element displacement vector, leaving just the straining part, where strains are assumed small such that linear theory can be used.Thus, this approach is suitable for large displacements and rotations, but small strain.It also avoids convergence issues with void elements, as shown by Dunning (2020).In this work, 2D plane stress problems are used to quantitatively compare stability constraint formulations, using the co-rotational element developed by Crisfield and Moita (1996).The element is based on a 4-node bilinear plane stress element with incompatible modes, which improves accuracy in bending by avoiding shear locking (Cook et al. 2001).

Iterative schemes
The Newton-Raphson method with automatic load increment adjustment and under-relaxation is used as the primary method for iteratively finding equilibrium points for the two constraint formulations based on eigenvalue analysis (sections 2.2 and 2.3).
For the other two constraint formulations, an arc-length path tracing method is utilized to directly locate the critical load factor (section 2.1), or its approximation by stiffness reduction (section 2.4).The arc-length method utilized in this work follows the implementation proposed by Lam and Morley (1992).This includes a method for automatically adjusting the arc-length for each step, based on a desired number of iterations to find equilibrium within each step: where Δl n+1 and Δl n are the arc lengths for the next and cur- rent steps, respectively, T n is the number of iterations used to find equilibrium for the current step, and T d is the desired number of iterations.The parameter T d is set by the user to optimize for efficiency, where a lower value leads to smaller increments in the arc-length each step, but fewer iterations within each step to find equilibrium.It can also be set to a smaller value to achieve a higher resolution of the equilibrium path, which can help avoid mode skipping.
For the direct method, at each converged equilibrium point an eigenvalue analysis is performed using Eq. ( 9) and two checks are made.If the load factor has decreased, then we have passed a limit point.However, if the load factor has increased, but the lowest positive eigenvalue is less than 1, (20) then we have passed a bifurcation point.In both cases we return to a point just before the critical one and use the information at that point to closely approximate the actual critical load factor ( γ ≈  1  ), and compute gradients using either Eq. ( 4) for a limit point, or Eq. ( 5) for a bifurcation point.
The approach for finding a critical point using the direct method is shown in algorithm 1.Note that the approximation could be improved with further iterations, e.g., using the bisection method.However, in practice a small value of T d is used in Eq. ( 20) to help avoid mode skipping, which usually results in small increments of near the critical point.Thus, the approximation is usually close, indicated by 1 ≈ 1 .An alternative is to use the extended system (Kemmler et al. 2005).This approach aims to solve Eqs.(1 & 2) simultaneously by iteratively solving an extended system of equations that also includes a normalization of the eigenvector.To solve the extended system, a penalty formulation is often required to avoid issues with the tangent stiffness matrix becoming singular near the limit point.It is possible that this approach is more efficient and robust than the approach used in this work, but it depends on the initial guess and penalty formulation used, and investigating this is beyond the scope of this paper.The value of , Eq. ( 15), for the stiffness reduction method is computed at each converged equilibrium point.When is lower than the target value, the load factor at the current and previous converged equilibrium points provide an upper and lower bound on .The bi-section method is then used to refine the error on to be within 0.01 of the target value.
For all methods, the force residual is used to determine convergence to equilibrium: where is a small number, chosen conservatively as 5 × 10 −9 .

Void element treatment
The C-shaped benchmark problem, introduced by Yoon and Kim (2005), is used to demonstrate that distorted void elements do not prevent convergence when using the co-rotational method.The problem is solved with, and without the strain energy interpolation method introduced by Wang et al. (2014), which has become a popular technique for avoiding convergence issues caused by distorted void elements, e.g., Dalklint et al. (2021).Details of how the strain energy interpolation method is adapted for the co-rotational method are included in Appendix 1. Comparing the solutions with and without strain energy interpolation, Figure 2, we can see the converged position of the solid region is almost identical (and almost identical to the solution where the void region use arc-length method to solve for: u i+1 , γ i+1 solve eigenvalue problem, Eq. ( 9), at: Limit point, use Eq. ( 4) for sensitivities Bifurcation point, use Eq.(5) for sensitivities else Algorithm 1 Direct method for finding the critical load factor is not modeled).Furthermore, 20 iterations are required to reach convergence in each case.This demonstrates that the co-rotational method does not suffer from convergence issues due to distorted void elements.However, we can see that some void elements are distorted in the solution without strain energy interpolation, Fig. 2b, and the issue of fictitious unstable modes in void regions is not addressed by the corotational method.Therefore, for eigenvalue analysis only, the tangent stiffness matrix is formulated using the strain energy interpolation method, which is effective at increasing the eigenvalues of fictitious modes so they are not included in the set of lowest modes used for the constraint formulation, as observed by Dalklint et al. (2021).

Topology optimization
A three-field density-based topology optimization method is used.First a design variable, x e ∈ [0, 1] is assigned to each element, e, in the mesh.A density-based filter is then applied to avoid well-known issues in density-based topology optimization (such as checkerboard patterns and mesh dependent solutions) (Bruns and Tortorelli 2001;Bourdin 2001).
where N i,e is the set of elements i where the distance between the center of element i and center of element e, d(e, i), is less than the filter radius r min .H i,e is a weight factor, defined as: The filtered field x creates a band of intermediate density values around the boundary of the structure.Thus, a Heaviside projection function is used to obtain a more crisp design description (Xu et al. 2010;Wang et al. 2011 where controls the sharpness of the projection.Initially, = 2 for the first 40 iterations, when it is then increased to 4. It is then increased by 2 every 20 iterations to a maximum of 12. Hence, the optimization runs for at least 120 iterations.This continuation approach helps avoid slow progress and getting stuck in poor local minima near the start of the optimization.The optimization process is stopped when the relative change in objective function between two iterations is less than 10 −4 , and all constraints are satisfied. The value of xe is used to compute the physical properties of the element, where the volume is directly proportional to xe and the modified SIMP (Simple Isotropic Material with Penalization) method is used to compute the Young's Modulus: where E 0 is the Young's modulus of the material, E min is a small value to prevent the stiffness matrix becoming singular ( E min = 10 −9 E 0 ) and p is the penalization factor, where the typical value of p = 3 is used in this paper.Note that the physical pseudo-density values xe are used to plot solutions.All optimization problems are solved using the Method of Moving Asymptotes (MMA) (Svanberg 2002).

Quantitative comparison
Two benchmark problems are used to quantitatively compare the four different methods for computing geometrically nonlinear stability constraints, as described in Sect. 2. These are a cantilever, and a problem where the linear solution exhibits snap-through behavior (Buhl et al. 2000), Fig. 3.The cantilever is discretized by 160 × 40 square elements, with a filter radius of 3 times the element edge length, and f ext = 100 kN.The snap-through problem is discretized by 240 × 80 square elements, with a filter radius of 2 times the element edge length, and f ext = 30 kN.In both cases, ele- ment thickness is 0.1 and material properties are a Young's modulus of 3 × 10 9 and Poisson's ratio of 0.4.In this section we compare stability constraint formulations by first solving the linear problem -minimization of linear compliance, subject to a volume constraint.The volume constraints are 40% and 10% for the cantilever and snap-through problem, respectively, and the linear compliance solutions are shown in Fig. 3.The solutions are then analyzed using each of the four geometrically nonlinear stability constraint formulations, comparing accuracy, efficiency, and robustness for computing the critical load factor only.Numerical experiments with geometrically nonlinear stability constraints applied during the optimization are presented in the next section.
Geometrically nonlinear analysis of the linear compliance solution of the cantilever is shown in Figure 4.The load-displacement paths are assembled from a number of nonlinear analyses and reveal complex behavior.The first unstable point is at ≈ 3.0 , where a bifurcation occurs, with another bifurcation at ≈ 4.4 .Thus, there are multiple pos- sible equilibrium states for some load factor values, demonstrating the difficultly in robustly computing the critical load factor, as may end up on different paths, depending on the solution method and its settings.However, the direct approach used here does accurately compute the critical load factor as 2.964.
Critical load factors estimated using eigenvalue analysis at various load factors are shown in Figure 5, where γ =  1 ×  .For load factors less than the critical value, the eigenvalue analysis provides a good estimate on the critical value, with errors less than 3%.However, for load factors larger than the critical value, the estimate depends on which branch of the equilibrium path is taken by the iterative method.If the lower branch is taken, then the estimate is still reasonable, with errors growing to about 9%, but if the higher branch is taken, then the eigenvalue analysis starts to identify the next bifurcation point as the critical load factor.At load factors larger than this, the eigenvalue analysis may also identify another unstable point as the critical one, at ≈ 6.0 .This demonstrates how the eigenvalue analy- sis approach may not be robust for load factors larger than the critical one, which may cause convergence problems, especially for the approach using an eigenvalue analysis at = lim (section 2.3).
Results for the stiffness reduction approach are shown in Fig. 6 for various values of .In this case, the iterative scheme locates the first bifurcation point by identifying the stiffness reduction associated with the lower branch.The accuracy of the estimated critical load factor increases, as decreases (from 2.2% for = 0.9 to 0.2% for = 0.1 ), which is expected.The efficiency of each method is evaluated by comparing the number of linear solves required to compute the critical load factor, in addition to the 4 solves required to reach equilibrium at = 1 .The direct method required 117 solves when using a very conservative step size in the arc-length method, by setting T d = 2 in Eq. ( 20).However, if the step size is increased by setting T d = 3 , then mode skipping occurs along the lower branch and the direct method computes a critical load factor of 5.54.This highlights the difficult balance between robustness and efficiency when using the direct approach, where a smaller step size leads to improved robustness, at the cost of more computational time.For the eigenvalue analysis methods, no additional linear solves are required for an evaluation at = 1 , and the additional cost is relatively small for load factors lower than the critical one (ranging from 3 at = 1.2 , to 10 solves at = 2.6 ).However, sig- nificantly more solves are required if the iterative scheme starts to move along the second branch, as shown in Figure 5.For the stiffness drop method, significant additional linear solves are required to satisfy Eqs.(1, 15 and 16) using T d =2, as shown in Figure 6.The number of solves for the stiffness drop method can be reduced by by increasing T d , but even if this parameter is tuned, it still typically requires 10's of additional solves.
Geometrically nonlinear analysis for the linear compliance solution to the snap-through problem is shown in Figure 7.A bifurcation occurs at = 2.93 , where the two members either bend down, or up.In the latter case, snapthrough occurs at = 5.09 .The direct approach accurately identifies the bifurcation point as the critical load factor, with 18 linear solves ( T d = 2).For this example, eigen- value analysis accurately identifies the bifurcation point for load factors up to the snap-through point, as shown in Fig. 8.This is achieved with only a few extra linear solves (3 to 6) until = 5 , where convergence becomes difficult due to being near the snap-through limit point.For load factors greater than the snap-through limit point, the accuracy of the method depends on which branch is followed after the bifurcation.If the branch with the members bending down is followed, then the eigenvalue analysis still identifies the bifurcation point, otherwise, a post snap-through instability is identified as the critical point.Therefore, both eigenvalue approaches for formulating the stability constraint will work in this example, assuming the constraint value, lim , is less than the snap-through limit load value.However, as expected, the stiffness drop method skips the bifurcation point and computes a critical load factor close to the snap-through limit point, with increasing accuracy (but greater computational cost) as gets smaller, Fig. 9.
The analysis of these two numerical examples confirm the qualitative observations summarized in Table 1.The cantilever example shows how mode skipping can occur when using the direct method, and eigenvalue analysis approaches -especially when eigenvalue analysis is performed past a limit point.Whereas the snap-through problem shows how the stiffness drop approach can skip over bifurcation points.Therefore, no method can guarantee robustness, but approximating the critical load factor using an eigenvalue analysis at = 1 is the most robust approach, as it is the least likely to suffer from mode skipping.It is also the most efficient, only requiring an additional eigenvalue solve (assuming equilibrium at = 1 has already been computed).However, approximation by eigenvalue analysis at = 1 may not be as accurate as the direct approach, or eigenvalue analysis at = lim .This is investigated further in the next section using optimization examples.

Optimization examples
Two benchmark problems are solved considering geometric nonlinearity during the optimization.The problems are solved using each of the four nonlinear stability constraint methods detailed in section 2, and also without any stability constraint.Note that T d =3 in Eq. ( 20) for the arc-length method, 6 modes are used in the K-S aggregation function for eigenvalue methods, Eq. ( 11), and the stiffness reduction method uses = 0.5 , Eq. (15).

Cantilever optimization results
The first example is the cantilever problem, Fig. 3a.The objective is to minimize end compliance: C = f T ext u 1 , with a 40% volume constraint and a stability constraint, with lim = 2 .It is solved for three values of the applied load: 100kN, 200kN, and 300kN, to test the constraint formulations under increasing amounts of nonlinear deformation.
The results for 100kN are shown in Figure 10.In this case the stability constraint is inactive and the solutions are almost identical, with objective function values within 0.2%, and the results are obtained in a similar number of optimization iterations (between 124 and 129).When the load is increased to 200kN the designs are again similar, with objective values within 5%, but there are some differences, see Fig. 11.Table 2 shows that the result using an eigenvalue analysis at = 1 took 127 iterations, whereas all other methods took more than 200 iterations.Slow convergence is caused by mode switching and mode skipping, demonstrated by the erratic optimization histories for the stability constraint shown in Fig. 12.
All final designs are analyzed using each of the four constraint formulations, with predicted critical load factors shown in Table 2. Firstly, the design obtained without any stability constraints does not meet the constraint ( lim = 2 ) using any of the methods, demonstrating that the constraint is active for this problem.The result using the direct method computes a critical load factor of 4.491, but further analysis shows that a bifurcation mode is skipped at ≈ 2 .The design is feasible, but mode skipping causes oscillations in the critical load factor during optimization, leading to slow convergence, see Fig. 12a.The design obtained using an eigenvalue analysis at = 1 predicts a critical load factor of 2.035 using this method, but the true critical load factor is approximately 5% lower, as computed by the direct method (with similar values computed using all other methods).However, using an eigenvalue analysis at = 1 is the most efficient approach (127 optimization iterations and no additional linear solves per iteration, Table 2), and avoids mode skipping and mode switching, resulting in a smooth convergence, Figure 12b.The design obtained using an eigenvalue analysis at = 2 is not feasible, as this method predicts a critical load factor of 1.986.For this problem, when the stability constraint becomes active during the optimization, it becomes more difficult for the iterative scheme to find an equilibrium point at = 2 due to the existence of a limit point.The method then suffers from the same drawback as the direct approach, as the tangent stiffness matrix close to = 2 becomes singular.The result using the stiffness reduc- tion stability constraint with = 0.5 skips a bifurcation point at ≈ 1.8 , which is identified by all other methods.
Designs obtained for the 300kN load are shown in Fig. 13 which have some similarities, and objective values are within 5%.The main difference is the design obtained using an eigenvalue analysis at = 1 , which has some addi- tional thin bracing members.The design obtained without stability constraints has a bifurcation point at ≈ 0.7 .The designs obtained using stability constraints are all feasible for their own method of computing the critical load factor, but none are feasible for all constraint formulation methods, Table 3.For example, the result obtained using the direct method skips a bifurcation at ≈ 1.2 .A similar issue occurs for the results obtains using an eigenvalue analysis at = 2 and stiffness reduction with = 0.5 .The accuracy of the direct method could be improved by taking smaller steps in the iterative scheme, but this would increase computational cost.Note that this method already uses an additional 226 linear solves per iteration (on average), Table 3.The result obtained using an eigenvalue analysis at = 1 (Figure 13c) overestimates the true critical load factor (as predicted by the direct method) by approximately 15%.However, this method does not suffer from mode switching, or skipping, as demonstrated by the reasonably smooth convergence (Figure 14b).
In summary, the observations in sections 2 and 5 are also seen during optimization.Mode skipping and switching affect the convergence of the direct method, mode skipping affects the convergence of the method using = lim , and these two methods also encounter difficulties with evaluating accurate derivatives near a limit point.The stiffness reduction method is also affected by mode skipping, and may not identify bifurcation points.The method using eigenvalue analysis at = 1 is the most robust, and most efficient, but may overestimate the true critical load factor.

Column optimization results
The second example is a compressed column, shown in Fig. 15a, and discretized by 120 × 240 square elements.The objective is to minimize volume, with an end compliance  The result without a stability constraint is a simple column (Fig. 15b), which has a bifurcation point at = 0.345 .Note that, in this case, including geometric nonlinearity in the stiffness (end compliance) constraint does not avoid the bifurcation instability, which highlights the need for a stability constraint.The result using the direct method, Fig. 15b, has a wider base and additional reinforcement members, which increases the predicted critical load factor to 2, just satisfying the constraint, although there is a bifurcation at ≈ 1.97 .This occurs because the implementation of the direct approach in this paper approximates the critical load factor from a pre-critical point i using γ ≈  i 1  i , which slightly overestimates the true value in this case.For this example, the direct method has slow convergence, taking 372 iterations to find a solution, Table 4, which is caused by mode switching, also seen in the oscillations of in the optimization history, Fig. 16a.When using the eigenvalue constraint formulation at = 1 , mode switching and skipping is avoided and a solu- tion is found in 148 iterations, Table 4.The design has similar features to the one obtained using the direct method, with a wider base and additional reinforcements, Figure 15d.However, when the direct method and eigenvalue analysis at = 2 are used to estimate the critical load factor for this design, they both show a bifurcation point at  < 1.9 , Table 4.This again shows that an eigenvalue analysis at = 1 provides an approximate value for the critical load factor.
For this example, no feasible solution is found when using the eigenvalue constraint formulation at = lim .This is caused by a form of mode switching where the eigenvalue analysis and sensitivity information depend on which path is followed after a bifurcation is passed at  < 2 , leading to oscillations in , see Figure 16c.In this case, the column either leans to the left, or right post bifurcation.The design shown in Figure 15e is obtained after 500 iterations, which has some similar features to the design obtained using the direct method and eigenvalue analysis at = 1 , but also some clearly unresolved regions of intermediate density.
The solution obtained using the stiffness reduction method is almost identical to the solution obtained without imposing a stability constraint (compare Figure 15b with 15f).This is because the stiffness reduction method does not identify the critical bifurcation point, and the stability constraint is effectively inactive, Figure 16d.

Conclusion
This paper compares four approaches for formulating stability constraints in geometrically nonlinear topology optimization, in terms of accuracy, robustness, and computational efficiency.The methods are the direct approach, eigenvalue  analysis at a load factor of 1, eigenvalue analysis at a load factor equal to the constraint limit value, and an implicit method based on a reduction in stiffness (compared with the initial stiffness).Note that using an eigenvalue analysis at the constraint limit load is newly proposed in this paper, and the stiffness reduction method has not previously been used to The qualitative and quantitative analysis shows that the direct approach and eigenvalue analysis at the constraint limit load can be the most accurate, but are not robust, as they may not identify the critical load value due to mode skipping.The direct approach and stiffness reduction method may also experience mode switching.In addition, the stiffness reduction method may miss bifurcation points.The method using an eigenvalue analysis at a load factor of 1 is the most robust, as it is least likely to experience mode skipping, or switching, and is the most efficient, as no additional linear solves are required to evaluate the constraint (only a single eigenvalue analysis).However, is may over, or underestimate the true critical load factor due to the eigenvalue approximation.
In conclusion, none of the methods studied in this paper completely satisfy all criteria (robustness, efficiency, and  accuracy).However, from the methods studied, it recommended that the approach using an eigenvalue analysis at a load factor of 1 is used to formulate stability constraints for geometrically nonlinear topology optimization, due to its robustness and efficiency.However, a challenge remains to improve its accuracy, which is recommended as a future research direction.Alternatively, techniques for improving the robustness and efficiency of the direct approach could also be investigated.For example, deflation methods (Xia et al. 2020) could be used to avoid issues with mode skipping, and several eigenvalues could be used to augment the constraint formulation to reduce issues with mode switching (Kemmler et al. 2005).To improve the efficiency of the direct approach, alternative methods for locating the critical load factor could be explored, such as using the extended system (Kemmler et al. 2005), and reanalysis or surrogate modeling approaches could also be developed.

Appendix 1: Strain energy interpolation
The strain energy interpolation method introduced by Wang et al. ( 2014) is based on a smooth interpolation of strain energy density within an element, where a fully solid element has strain energy associated with a geometrically nonlinear analysis, and a void element has strain energy associated with a linear analysis: where i (u e ) is the strain energy density for element e, (⋅) is the nonlinear strain energy, L (⋅) the linear strain energy, and e is a parameter dependent on the element physical pseudo-density, xe : For the co-rotational method, the element tangent stiffness matrix with strain energy interpolation becomes: where K(⋅) is the nonlinear tangent stiffness matrix, and K L is the linear stiffness matrix.Note that because K L is inde- pendent of the load factor, void elements have no contribution to K u or K s when formulating the eigenvalue problem, Eq. ( 10), which effectively eliminates fictitious modes in void regions.

Fig. 1
Fig.1Illustration of mode skipping when using an iterative method to directly compute the critical load factor Fig. 2 C-shaped benchmark test.a problem setup, b solution without strain energy interpolation, c solution with strain energy interpolation

Fig. 3
Fig. 3 Benchmark problems, with linear compliance solutions.a cantilever, b snap-through problem

Fig. 5
Fig. 5 Cantilever linear compliance solution -critical load factor estimation using eigenvalue analysis at different load factors

Fig. 7
Fig. 7 Geometrically nonlinear analysis of the snap-through solution for linear compliance

Fig. 10
Fig. 10 Results of cantilever with f ext = 100kN using: a no stability constraints, b direct method, c eigenvalue analysis at = 1 , d eigenvalue analysis at = lim , e stiffness reduction ( = 0.5)

Fig. 12
Fig. 12 Optimization history for cantilever with f ext = 200kN using: a direct method, b eigenvalue analysis at = 1 , c eigenvalue analysis at = lim , d stiffness reduction ( = 0.5)

Fig. 13
Fig. 13 Results of cantilever with f ext = 300kN using: a no stability constraints, b direct method, c eigenvalue analysis at = 1 , d eigenvalue analysis at = lim , e stiffness reduction ( = 0.5)

Fig. 16
Fig. 16 Optimization history for compressed column using: a direct method, b eigenvalue analysis at = 1 , c eigenvalue analysis at = lim , d stiffness reduction ( = 0.5)

Table 1
Qualitative comparison of geometrically nonlinear stability constraints -summary of key points

Table 2
Analysis of cantilever optimization results for 200kN load

Table 3
Analysis of cantilever optimization results for 300kN load

Table 4
Analysis of compressed column optimization results