Structural stability and artificial buckling modes in topology optimization

This paper demonstrates how a strain energy transition approach can be used to remove artificial buckling modes that often occur in stability constrained topology optimization problems. To simulate the structural response, a nonlinear large deformation hyperelastic simulation is performed, wherein the fundamental load path is traversed using Newton’s method and the critical buckling load levels are estimated by an eigenvalue analysis. The goal of the optimization is to minimize displacement, subject to constraints on the lowest critical buckling loads and maximum volume. The topology optimization problem is regularized via the Helmholtz PDE-filter and the method of moving asymptotes is used to update the design. The stability and sensitivity analyses are outlined in detail. The effectiveness of the energy transition scheme is demonstrated in numerical examples.


Introduction
Designs optimized for weight or compliance often consist of slender members which are inherently prone to buckle. This failure mode has been considered in the optimization of truss and beam structures, dating back to the work of Khot et al. (1976) where the weight of space truss structures was minimized subject to stability constraints and Olhoff and Rasmussen (1977) who presented closed form expressions for the sufficient optimality condition of statically indeterminate columns. More recently, Amir (2017, 2019) improved the accuracy of the buckling predictions by using geometrically nonlinear beam theory to optimize stable truss structures.
Although stability is crucial for a robust design, it is often neglected in topology optimization of continuum structures due to its complicated analysis. Indeed, an unstable equilibrium configuration is characterized by a singular tangent stiffness matrix, and therefore the most accurate buckling load estimate is obtained by performing an incremental loading until at least one eigenvalue of the tangent stiffness matrix becomes zero (see, e.g. Suleman and Sedaghati (2005) and Pedersen and Pedersen (2018)). In topology optimization, this is a costly computation as it must be performed for each design iteration, and its accuracy is adversely affected by the artificial buckling modes that arise due to the ersatz material model.
To lessen the computational burden associated with the buckling analysis, it is common to employ a so-called "linear prebuckling analysis", in which the buckling load estimate is obtained from a generalized eigenvalue problem about the structure's undeformed configuration. Research on topology optimization that incorporates linear prebuckling analysis include the seminal work of Neves et al. (1995) where the lowest critical load is maximized. Zhou (2004) and Bruyneel et al. (2008) discuss convergence issues encountered in density based topology optimization. The necessity of including a large number of buckling constraints to obtain monotonic convergence is scrutinized by Dunning et al. (2016), as well as the need to use an efficient Published online: 21 August 2021 / eigenvalue solver. Recently, Sigmund (2019, 2020) investigated the use of nonconforming finite elements and an aggregated stability constraint to solve large scale buckling constrained optimization problems.
Unfortunately, designs that are based on linear prebuckling analysis may not behave as expected due to the crude nature of the approximation. For more accurate buckling load estimates, a fully nonlinear analysis is required. In the nonlinear setting, there are two ways to identify the buckling load levels. The first, known as the extended system approach, directly computes the critical configurations by finding the unstable equilibrium solutions, cf. Wriggers (2008). The second approach, known as the one point approach, utilizes the information at a stable equilibrium configuration and approximates the buckling load level by solving a generalized eigenvalue problem based on a decomposition of the tangent stiffness matrix, cf. Brendel and Ramm (1980). The extended system approach provides an exact measure of the critical loads up to a numerical tolerance; however, it is only able to identify a single critical load. The one point approach does not suffer from this restriction and is therefore adopted in this work, despite the fact that it is an approximation. Reitinger and Ramm (1995) and Kemmler et al. (2005) employ the extended system approach to obtain critical buckling loads in their topology optimization work. The one point approach is used in continuum topology optimization by, e.g. Lindgaard and Dahl (2013), where different buckling load and compliance objectives are investigated, and by Lund (2010, 2011) to optimize stable composite structures. For reasons discussed below, a large number of buckling constraints must be considered to obtain a smooth convergence of the optimization problem and hence we employ the one point approach.
Unfortunately, the solution to eigenvalue problems in topology optimization with an ersatz material model is complicated due to the existence of the low volume fraction regions with minimal stiffness. Nonphysical, artificial instabilities appearing in the low volume fraction regions must be distinguished from the actual physical modes. Several methods for accommodating the artificial buckling modes exist. A simple remedy used in the linear buckling analysis is presented by Neves et al. (1995), wherein the geometrical stiffness of low stiffness elements is ignored. This method might, however, cause convergence issues in the optimization due to discontinuities in the gradient information, and therefore Bendsøe and Sigmund (2003) propose an interpolation scheme to gradually remove the low stiffness elements. Both Zhou (2004) and Gao and Ma (2015) discuss the difficulty of choosing the interpolation parameters such that only the physical buckling loads are obtained.
The fictitious domain-based topology optimization considering natural frequencies in Du and Olhoff (2007) and Pedersen (2000) use separate interpolations of the stiffness and mass matrices. This idea could potentially be applied in the one point buckling analysis. However, Lindgaard and Dahl (2013) found that different interpolations for different stiffness matrix contributions result in loss of convergence in the equilibrium equations. Instead, Lindgaard and Dahl (2013) used the same interpolation for all components of the stiffness matrix, and eliminated the artificial modes from the analysis by increasing the stiffness of the void regions. Of course, this method results in loss of precision in the sense that void regions will be artificially stiff. As an alternative, Gao and Ma (2015) distinguished between artificial and physical buckling modes by using a strain energy density criterion. We note that this approach is computationally intensive.
The low volume fraction regions pose additional problems when subjected to finite deformations, as excessively distorted finite elements might result in loss of convergence in the Newton equilibrium iterations. Therefore, in a fully nonlinear analysis, it is attractive to find an approach which deals with both the artificial modes and the equilibrium divergence simultaneously. One such approach is the element removal strategy, as done by Dalklint et al. (2020) in their topology optimization of preloaded elastic continuum structures subject to natural frequency constraints. Herein, we use the energy transition scheme presented by Wang et al. (2014), whereby the physical strain energy is replaced by a linearized strain energy in void regions. As such, we eliminate the artificial modes and avoid equilibrium convergence issues.
Summarizing, in our work, we obtain accurate estimates of the buckling loads of compressible neo-Hookean hyperelastic structures by using the one point approach wherein a generalized eigenvalue problem about the deformed equilibrium configuration is solved to approximate the buckling loads. The deformed configuration is found by incrementally loading the structure to approach, but not reach, an unstable configuration. Newton's method facilitates this analysis. The optimization minimizes the end-compliance at the operating load, subject to lower bound constraints on the n-lowest critical loads and an upper bound constraint on the mass. The method of moving asymptotes (MMA) nonlinear programming algorithm is used to update the design, and the sensitivity analysis is performed using the adjoint approach. To address the nondifferentiability of degenerate eigenvalues, we follow the work of Chin and Kennedy (2016) and Ferrari and Sigmund (2019) and replace the, possibly nondifferentiable, n eigenvalue constraints with a differentiable aggregated constraint.

Governing equations
To consider finite deformations, we need to distinguish between the undeformed, reference configuration, Ω, and the deformed, current configuration, Ω c . Assuming Lagrangian elasto-statics, each material point identified by the position vector X ∈ Ω in our n dim -dimensional space, is displaced to the position x ∈ Ω c , by the deformation mapping ϕ : Ω → Ω c , i.e.
where u : Ω → R n dim represents the displacement field.
To describe the strain in this nonlinear theory, we define the deformation gradient where ∇ is the material gradient operator on Ω and 1 is the second order identity tensor. The local deformation is described by the right Green-Cauchy deformation tensor C = F T F . For later purposes we also introduce the linearized strain = 1 2 (∇u) T + ∇u . Under the assumption of vanishing body forces, the potential energy of our hyperelastic body is where W = W (∇u) is the specific strain energy and the boundary ∂Ω with unit normal n is partitioned into complementary subsets ∂Ω t and ∂Ω u , over which the traction t =t and the displacement u = ϕ − X =ū, are prescribed, respectively. The equilibrium configurations of Ω are described by the admissible displacement fields u ∈ {u ∈ H 1 : u(X) =ū for X ∈ ∂Ω u } that minimize the energy, i.e. satisfy the stationary condition for all smooth virtual displacements δu ∈ {δu ∈ H 1 : δu(X) = 0 for X ∈ ∂Ω u } and where H 1 is a Hilbert space.

Material interpolation and design variables
The goal of topology optimization is to optimally distribute material in a design domain Ω. To quantify the material distribution, we use the material indicator function Γ : Ω → {0, 1}. However, as is standard in topology optimization, we replace the binary valued material indicator function with a continuous nondimensional volume fraction field z : Ω → [0, 1] and to obtain a well-posed optimization problem, we restrict the design space. Here, the restriction is performed by penalizing fine-scale oscillations of z via the Helmholtz PDE-filter proposed in Lazarov and Sigmund (2011). As such, z is replaced by the smoothed field ν : Ω → [0, 1], obtained from the solution of the boundary-value problem where Δ is the Laplacian with respect to Ω and the oscillations of z are further smoothed by increasing the value of l ∈ R + . The homogeneous Neumann boundary condition in (5) is well-known to cause an artificial "attraction" of material to the domain boundaries. This anomaly can be mitigated by applying a Robin condition, cf. Wallin et al. (2020). In our discretization, we parameterize z via the piece-wise uniform design field z : Ω → [0, 1] over the finite elements Ω e ∈ Ω, such that z(X) = z e for X ∈ Ω e . We use standard Lagrange finite elements to solve (5) for the continuous smoothed field ν, cf. Lazarov and Sigmund (2011) for details. Unfortunately, the filtering produces gray regions wherein ν(X) ∈ (0, 1). To limit the extent of these regions, we implement the thresholding technique proposed by Guest et al. (2004) and Wang et al. (2014), wherein a smoothed Heaviside function H β,η : R → [0, 1] is employed to define the thresholded filtered volume fraction where β and η are numerical parameters defined such that lim β→∞ H β,η (ν) = u s (ν − η) where u s is the unit step function. The filter (5) and the Heaviside approximation (6) are used to obtain the differentiable fieldν. Indeed, although z is a piecewise uniform function, ν andν are continuous functions that vary spatially throughout the domain Ω.
To further reduce the size of the gray regions for which ν(X) ∈ (0, 1) we employ the SIMP penalization scheme, cf. Bendsøe (1989), where a nonlinear scaling of the elastic material parameters, K and G, which describe the isotropic hyperelastic material model is introduced as In (7), G o and K o denote the shear and bulk moduli of the optimized structure's isotropic material, in the limit of infinitesimal strain. The SIMP scaling function χ : whereby increasing values of q > 1 enforce increasing levels of penalization and the lower threshold δ o = 10 −9 in this ersatz material model is used to avoid numerical problems otherwise caused by a vanishing material stiffness. The aforementioned material interpolation allows for nearly void regions to be part of the physical domain. As such, even minimal stress leads to excessive distortion in these regions resulting in loss of convergence in the finite element equilibrium simulations. To mitigate this issue, we use the energy transition strategy proposed by Wang et al. (2014). In this scheme, the deformation gradient is redefined as F = 1 + γ ∇u, where γ : [0, 1] → [0, 1] is another penalty function defined as whereβ = 500 andη = 0.01 replaces β and η in (6).
To formulate stress and stiffness relations consistent with the strain energy transition, we weight the physical strain energy, W P H , against the linearized energy, W L , such that the body's internal energy is By assigning points devoid of material γ = 0, we therefore obtain W = W L . The physical strain energy is modelled using an isotropic compressible neo-Hookean hyperelastic constitutive model such that Its linearization yields the linearized energy where D L is the linear stiffness tensor corresponding to an isotropic material, with bulk and shear modulii K o and G o .

Total Lagrangian FE-formulation and solution procedure
Our equilibrium equations are discretized using standard finite elements, see, e.g. Crisfield (1993). We use to represent the finite element assembly operator, lower case indices e ∈ N e to indicate quantities associated with finite element Ω e , and lower case n ∈ N n to denote the finite element degrees-of-freedom. The sets of finite element degrees-of-freedoms and finite elements, are denoted N n and N e , respectively.
We seek the discretized displacement field u ∈ R n dim that approximates the displacement u via standard Lagrange shape functions N within the element Ω e such that u(X) = N(X)a e , where a e is the element nodal degrees-of-freedom vector. The finite element formulation is based on the Galerkin approach, implying δu(X) = N(X)δa e , where δa e is the virtual element nodal degrees-of-freedom vector.
To discretize (4), we first note that from (10) the virtual stress power is where are the physical and linearized stress tensors. In (13), we define the virtual Green-Lagrange strain δE = γ 2 (∇δu) T F + F T ∇δu and the virtual linearized strain δ = 1 2 (∇δu) T + ∇δu . The discretization of δu is used to express these tensors as δE = B(γ a e )γ δa e = B c + B l (γ a e ) γ δa e and δ = B c δa e , where B is expressed as the sum of a constant B c and a linear function of the displacement B l (γ a e ), cf. Crisfield (1993).
Under the assumption of dead loading, we express the external load as λP where P ∈ R n is a unit magnitude external load vector and λ ∈ R + is the scalar load intensity factor. Using the arbitrariness of the virtual nodal displacements, δa, and (13), the discretization of (4) is expressed as where r(a) ∈ R n is the residual vector.
We find the degree-of-freedom vector a in (15) for the given λ. To this end, we employ Newton's method to solve (15). As such, we iterate solving K(a)Δa = −r(a) for Δa and updating a → a + Δa until convergence, where the tangent stiffness matrix K(a) ∈ R n×n is and where D = 4 ∂ 2 W P H ∂C∂C , G contains the gradient of the element shape functions, Y is a symmetric block matrix of S P H (γ a) and D, D L are the Voigt representations of D, D L respectively. After convergence, we increase λ and again evaluate a. We repeat this process until we reach the operating load, i.e. until λ = λ op .

Structural stability
A conservative system is stable if the potential energy has a (local) isolated minimum. This condition is equivalent to the requirement that the second variation of the potential energy is positive definite. For our system, this implies that all eigenvalues of the tangent stiffness matrix (15) are strictly positive. On the other hand, the tangent stiffness matrix has negative eigenvalues if the system is unstable. A system that transitions from being stable to unstable must pass a critical, i.e. singular point, for which at least one eigenvalue of K vanishes, i.e. when has a nontrivial solution φ φ φ j . Unfortunately, the zero eigenvalue condition is difficult to enforce in topology optimization due to the appearance of artificial eigenmodes attributed to the ersatz material. As such, the search for the zero eigenvalue would require many eigenvalue solves and a strategy to distinguish between artificial and physical eigenmodes.
In our alternate strategy, we find the critical physical eigenvalues by approximating the eigenvalue problem (17), as a generalized eigenvalue problem that is based on a decomposition of the stiffness matrix. Motivated by Wriggers (2008), we assume that K is the sum of the stress dependent contribution K g that varies linearly with respect to the load level in the vicinity of the critical point, and the constant contribution K o . To leverage this assumption, a scalar π ∈ R + is introduced such that λ c = πλ op where λ op ∈ R denotes the operating load level at which we evaluate where and In this work, we limit the analysis to buckling load factors π > 1, i.e. we assume that the decomposition (18) is performed in a stable equilibrium position, prior to reaching any instabilities.
We use the assumption (18), to approximate the eigenvalue problem (17) via the generalized eigenvalue problem 1 where (π j , φ φ φ j ) ∈ (R × R n ) are the eigenpairs. The eigenvalues π j are ordered in ascending order such that 1 In the computations we rewrite (21) as K g φ φ φ j = − 1 π j K o φ φ φ j , j ∈ N n since K g is indefinite and K o positive definite. π 1 ≤ π 2 ≤ . . . ≤ π j ≤ π n and the eigenvectors are normalized such that It is emphasized that in the case of N-fold degenerate eigenvalues, the eigenvectors are not uniquely defined as they span an N-dimensional hyperplane, i.e. any linear combination of φ φ φ j , j ∈ N N , satisfy (21) where the set N N ⊆ N n contains the N indices j for which π j = π j +1 = . . . = π j +N−1 .

Optimization formulation
The optimization problem is stated as finding the design z that solves where g 0 is the cost function to be minimized and g i (z) ≤ 0 are the i ∈ N c = {1, 2, . . .} constraint inequalities. A common cost function to minimize in topology optimization, and the one used herein, is the end-compliance defined by Since P is independent of the load magnitude, the objective minimizes the displacement over the regions where the load is applied. Alternatively, one could maximize the tangent stiffness as done in Wallin et al. (2018).
For all examples considered, the amount of available material is limited by imposing the constraint where V max ≤ V = Ω dv is the maximum volume available for the design. Ideally, the stability constraints are defined as where λ c i = π i λ op represent the critical load level, and N λ ⊂ N n is a subset of eigenvalues, i.e. buckling load levels, considered. In (26), we introduce constraints on the N λ first critical loads to avoid convergence issues associated with mode-crossing, cf. Dunning et al. (2016). The operating load level λ op and the safety factor s ≥ 1 are problem specific parameters. We follow Ferrari and Sigmund (2019) and Maharaj and James (2020) and employ an alternative representation of (26), whereby the N λ -constraints (26) are replaced by the single constraint In (28) (27) constraint is conservative. The use of (27) instead of (26) will be justified when the sensitivity analysis is discussed.

Sensitivity analysis
To evaluate the gradients of the objective and constraint functions, we use the adjoint procedure. In the followin, we obtain the sensitivities ∂g i ∂ν ν ν . But remember that the optimization requires the sensitivities with respect to the design variables z, i.e. ∂g i ∂z . These are obtained from the chain rule where the product ∂g i ∂ν ν ν ∂ν ν ν ∂z is obtained from an adjoint sensitivity analysis using (5), cf. Lazarov and Sigmund (2011).

Sensitivity of objective
The sensitivity of the objective function, i.e ∂g o ∂ν ν ν is obtained by augmenting g o with the zero valued residual vector by means of the adjoint vector ϕ ∈ R n , such that Differentiating (30) with respect to ν ν ν results in To eliminate the implicitly defined derivative ∂a ∂ν ν ν , we require ϕ to solve the adjoint problem where we utilize the symmetry of K. The final expression for the sensitivity of g 0 is thus reduced to 2 Other feasible choices of the norm exist, cf. Ferrari and Sigmund (2019).
where the notation( ·) here and henceforth denotes quantities which are treated as constants in the differentiation.

Sensitivity of critical loads
To obtain the sensitivity of the stability constraint, g 2 , we start by differentiating (27) ∂g which follows from (28), and where Λ j = 1 λ c j . As the operating load is constant, the sensitivity cf. the discussion preceding (17). The eigenvalue derivative ∂π j ∂ν ν ν is obtained by manipulating (21) and using the equilibrium residual equation to obtain the equality where μ ∈ R n is another adjoint vector. Differentiation of (36) with respect to ν ν ν results in Similar to (32), we require the adjoint μ j vector to solve which annihilates the implicit sensitivity, ∂a ∂ν ν ν , in (36). Insertion of (38) into (37) and making use of the normalization (22) results in The explicit computations in the right hand sides of (38) and (39) appear in the Appendix.
It is well known that (39) does not hold for degenerate eigenvalues as they are not differentiable. Nonetheless, the simple eigenvalue result is adequate because both K o and K g are real, symmetric and continuously differentiable with respect to ν ν ν and g 2 is a symmetric polynomial of the eigenvalues π j for j ∈ N λ . As such, the gradient ∂λ c j ∂ν ν ν exists even for the case of degenerate eigenvalues, by evaluating the sensitivity ∂g 2 ∂ν ν ν via (34), (35) and (39), cf. Torii and De Faria (2017) and Gravesen et al. (2011).

Artificial buckling modes
If the material properties in void regions are allowed to strictly assume the zero lower bound, i.e δ o = 0 in (8), the equilibrium equation cannot be solved. To avoid this, it is common to rely on an ersatz material model such that material points outside the physical domain are assigned a "small" stiffness, i.e. 0 < δ o 1. For eigenvalue problems in topology optimization, the small stiffness in void regions renders nonphysical, i.e. artificial eigenmodes.
This artificial mode anomaly is readily shown in the following example, wherein a portal frame structure, similar to that studied in Gao and Ma (2015), is analyzed, cf. Fig. 1. The design domain is discretized using 40 × 20 × 1 3D 8-node fully integrated displacement based finite elements. Following Gao and Ma (2015), plane stress conditions are simulated by assigning the out-of-plane displacement on the front x-y plane to zero and enforcing traction free boundary conditions on the opposing face. The operating load level is λ op = 60 N and the elastic parameters are E = 3 GPa and ν = 0.4 which provide G o = E/(2(1 + ν)) and K o = E/3(1 − 2ν). In this example, we disable the filter (5), and assign the outer frame, i.e. the dark grey physical domain of thickness 5 mm, a fully dense, i.e. ρ = 1, material, whereas Fig. 1 Portal frame structure. The thickness is 10 mm ρ = 0.001 for the inner "void" region (light grey). The SIMP penalization parameter is q = 3 and the Heaviside projection parameters are β = 1 and η = 0.5, cf. (8) and (6). We search for eigenvalues π j in the interval I = [1, 100]. The physical buckling modes and load factors, are evaluated in three different ways. Firstly, the computations are performed on a conforming mesh over the physical domain, i.e. excluding the inner grey region. Secondly, the grey void domain is modelled via the SIMP ersatz model, cf. (8), without the energy transition scheme, i.e. by setting γ = 1. Lastly, the void domain is modelled via the combined SIMP and strain energy transition schemes of (8) and (10). Various matrix penalizations in the generalized eigenvalue problem (21) exist, cf., e.g. Lindgaard and Dahl (2013). However, as noted by Gao and Ma (2015) and Lindgaard and Dahl (2013), these methods deteriorate the convergence of the equilibrium analysis and we therefore exclude them.
The smallest buckling load factor, π 1 , together with the number of modes, j max , found in the interval I , are presented in Table 1. As expected, we find a large number of artificial modes in the inner "void" region when using the second method. This behaviour is readily predicted by inspecting the Rayleigh quotient of (21) When employing the strain energy transition scheme, the denominator in the Rayleigh quotient (40) tends to zero as γ → 0, thereby removing the artificial modes from the interval I . The idea of identifying artificial modes using the Rayleigh quotient also appears in Gao and Ma (2015). The first buckling modes for the three simulations are depicted in Fig. 2. We emphasize that Fig. 2b displays an out-of-plane artificial buckling mode. Table 1 The magnitude of the smallest load factor π 1 and the number of eigenvalues j max found in the interval I for the three modeling schemes π 1 j max Conforming 6.0077 46 SIMP ersatz 2.3705 3325 SIMP strain energy transition 6.0078 46 Fig. 2 The first buckling mode when using a conforming mesh (a), the SIMP ersatz material (b) and combined SIMP and strain energy transition models (c) To conclude the investigation, we compute the ratio between the greatest positive physical eigenvalue and smallest artificial eigenvalue obtained using the proposed method. To this end, we extend I = [1, 10 20 ] and compute all positive physical eigenvalues using the first method, whereby we find 391 positive physical eigenvalues. Subsequently, we compute the eigenvalues in I using the third method. As expected, the first 391 eigenvalues are physical, whereas the rest are artificial. The first artificial eigenvalue is approximately 11 magnitudes larger than the greatest physical eigenvalue, which further demonstrates the efficacy of the proposed method.

Numerical examples
We solve the two optimization problems using the strain energy transition model and the MMA-scheme with default parameters, cf. Svanberg (1987). The problems are discretized with 2D quadrilateral plane strain fully integrated elements. The constitutive parameters are the same as in the previous example. The 10 lowest eigenvalues, i.e. buckling loads are considered in the (26) constraint. The filter (5) is enabled and a continuation scheme for the penalty factor, q, is utilized wherein the initial value of q = 2 is increased to the final value q = 6 in increments of 0.05 every 5th design iteration. We chose a conservative continuation scheme of q as this can be used for all numerical examples. After the terminal value of q = 6 is reached, we wait 20 iteration before the Heaviside parameter β, cf. (6), is increased from β = 1 to β ≈ 12 following the exponential continuation scheme β → 1.3 · β every β incr design iteration. The threshold parameter η = 0.5 in (6) remains fixed throughout the iterations.

Stubby cantilever example
The first optimization problem is the stubby cantilever beam design illustrated in Fig. 3 subject to a transverse load. The filter length scale is l = 40/2 √ 3 mm, cf. (5), and the maximum allowed volume is V max = 0.2V . The operating load is initially λ op = 20 kN, and increased by 20 kN each design iteration to the terminal load of λ op = 140 kN. After the terminal load is reached, we wait three additional iterations before introducing the stability constraint (27). The reason for not initially activating the stability constraint is to limit the occurrence of artificial buckling modes associated with stress concentrations in the support and load regions, see Ferrari and Sigmund (2019) for further discussion. The β update increment parameter is β incr = 10. The geometry is discretized using 90 × 200 = 18, 000 elements. The load vector, P, is uniformly applied to the 11 vertical degrees-of-freedom symmetrically located about the right center point.
To investigate the influence of the safety factor, we provide 10 examples with s in the range s ∈ [1.1, 2.0], cf.
(27). The optimized designs appear in Fig. 4. The results show that as the safety factor increases, i.e. minimum buckling load increases, the lower beam maintains its volume while increasing its second moment of inertia with respect to its neutral axis by developing a web-like structure.   Scaled first ten buckling modes φ φ φ j , j ∈ N λ for the Fig. 8, s = 2.0, design

Arch example
In this example, we design the arch depicted in Fig. 7. The filter length scale is l = 32/2 √ 3 mm, cf. (5), the maximum allowed volume is V max = 0.4V and β incr = 30. The operating load is initially λ op = 1.5 kN, and increased by 0.2 kN each design iteration until λ op = 2.7 kN. The artificial modes associated with stress localizations in load and support regions are not problematic in the arch compared to the stubby cantilever problem, and therefore the stability constraint is enabled after the first design iteration. The The design histories of the ten first load scaling factors π j for the Fig. 8, s = 2.0, design geometry is discretized using 666 × 50 elements. The structure is subject to the load defined by the load vector, P, which is uniformly applied to the 19 vertical degrees-offreedom symmetrically located about the top center point. We again investigate the influence of the safety factor s on the optimized design, cf. Fig. 8, where it is again apparent that as s increases so too does the second moment of inertia with respect to the neutral axes of the beams. Figure 9 depicts the first 10 buckling modes of the Fig. 8, s = 2.0, design. Once again we emphasize the effectiveness of the strain energy transition scheme in removing artificial modes. Figure 10 displays the histories of the ten first load scaling factors π j , j ∈ N λ for the Fig. 8, s = 2.0, design. We again observe oscillatory behaviour due to the q and β continuation strategies, but otherwise a stable convergence.

Conclusions
In this work, we demonstrate how the energy transition scheme presented by Wang et al. (2014) can be applied in buckling constrained topology optimization to eliminate artificial eigenmodes. When considering finite strain loadings, this scheme also improves convergence in the finite element simulations. The nondifferentiability of degenerate critical points is handled by p-norm aggregation, as done by Torii and De Faria (2017).
The numerical examples validate the efficacy of the energy transition scheme to remove artifical eigenmodes and produce stiff designs that are stable.

Replication of results
The authors state that all the data necessary to replicate the results are presented in the manuscript. Relevant parts of the code can be shared by contacting the corresponding author.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.