A level set topology optimization method for the buckling of shell structures

Shell structures are some of the most widely used in engineering applications. Flat plates, stiffened panels, and wing ribs are each examples of components for which the design features may be dictated by the critical buckling load. Despite this practical significance, there exists only a handful of studies in the literature documenting applications of topology optimization which consider buckling performance. This is due to several issues innate to this domain, including mode switching, spurious behavior in void regions, and the presence of repeated eigenvalues. Herein, we propose a level set method capable of effectively optimizing structures despite these challenges in the context of linear buckling. We demonstrate the usefulness of such in the design of several common shell structures and explore the trade-off between stiffness and buckling load performance.

= change in φ values at the boundary points p j = vector to integrate element j nodal values ρ i = mass density of domain i T = linear map between φ and ψ u = displacement vector u a = adjoint vector v = buckling eigenvector w j = volume fraction of finite element j w j = vector to integrate element j nodal values z = vector used to identify deflection along z-axis

Introduction
There is recent and renewed interest in topology optimization methods that include consideration of buckling behavior. This results from the combination of two factors: (1) Topology optimization is beginning to find application in more complicated structures such as aircraft components (Aage et al, 2017;Townsend et al, 2018) which are often buckling-driven, and (2) There is a documented tendency for structural optimization methods considering only weight and stiffness to produce designs with poor stability (Rozvany, 1996;Rahmatalla and Swan, 2003). Numerous efforts have centered on buckling-driven optimization of frame and truss structures, where the design variables typically comprise the cross sectional features of members in a fixed ground structure (Zhou, 1996;Rozvany, 1996;Guo et al, 2001). Relatively few works appear in the literature concerning buckling-driven topology optimization of continua, though such structures stand to benefit substantially from these algorithms. In continuous structures, tension and compression paths become more difficult to discern, increasing the difficulty of intuitive design. A recent note by Ferrari and Sigmund (2018) exposes the issues hampering routine consideration of buckling behavior in this context, and those relevant to the present study will be detailed here, along with the works on buckling-driven topology optimization to date. It is noted that we herein deal only with the linear buckling response of structures. Whilst consideration of non-linear effects can be important in estimating realistic failure loads (Pedersen and Pedersen, 2018), the authors consider the simplicity and computational efficiency afforded by linear analysis retains its relevance in, at least, preliminary design stages. Furthermore, development of topology optimization algorithms which circumvent issues related to linear buckling may make non-linear buckling design more tractable in future works.
The first issue encountered in buckling-driven optimization pertains to the multi-mode nature of the linear buckling equations (Cook et al, 2007): where K, f and u denote the structural stiffness matrix, applied load and stationary deflection respectively. K s denotes the stress stiffness (also called geometric stiffness) matrix. λ and v represent the eigenvalue/eigenvector pair for a given buckling mode; there are as many pairs as there are degrees of freedom in the system. Eigen-vectors are orthonormalized such that The first positive buckling load factor λ 1 is typically the only of practical interest, since it signifies the onset of instability (λ 1 > 1 implies stability for a given f ); the simplest optimization formulations thus aim to either minimize weight subject to a lower bound on λ 1 (the constrained λ 1 formulations), or maximize λ 1 subject to an upper bound on weight (the max. λ 1 formulations). Both forms are problematic, since numerous modes may be grouped closely, and the critical mode may switch between design iterations. Thus in order to avoid erratic convergence in constrained formulations, a large number of modes must be computed during the analysis stage, along with the associated derivatives. Such provides the optimizer with sufficient information to constrain all modes which are currently active, and would become so on the next iteration. For example, Bruyneel et al (2008) minimized the weight of a stiffened panel subject to buckling load constraints using 42 discrete design variables; for such, 100 buckling modes were required to be computed in order to stabilize convergence. The work by Dunning et al (2016) showed a similar trend in the topology optimization of continua. In that work, a level set method (details of such methods are the subject of Section 2) was applied to the weight minimization of a 3D bracket structure, wherein 5 mode constraints failed to achieve convergence due to mode switching, and a significant reduction in the number of iterations was observed by increasing the number of modes computed from 10 to 25. It is worth noting that for such 3D continuum structures, computing eigenmodes is computationally-expensive, and eigensolvers other than the commonly-used ARPACK (Lehoucq et al, 1998) may become preferable (Dunning et al, 2016). To avoid mode switching in max. λ 1 formulations, one may convert the problem to a bound formulation, whereby the objective becomes to maximize an auxiliary variable, with constraints added such that said variable remains less than all of the computed eigenvalues (Bendøse and Sigmund, 2003). Others have advocated an eigenvalue separation approach, whereby additional constraints are added to the problem such that adjacent modes remain sufficiently separated (Stanford et al, 2014;Ferrari and Sigmund, 2018). It should be noted that this technique prevents the identity of the critical buckling mode changing throughout the optimization, as modes cannot switch from their baseline ordering, which in theory could lead to sub-optimal results. When applied to the topology optimization of panel stiffeners (Stanford et al, 2014), the eigenvalue separation technique was found to result in faster convergence without loss in performance when compared to the aforementioned auxiliary approach. The use of an aggregation function, such as the popular KS function (Kreisselmeier and Steinhauser, 1980), also alleviates the issue, since when aggregated, the order of eigenvalues is unimportant. Such aggregation techniques are common in stress-based topology optimization (Verbart et al, 2017;Picelli et al, 2018), and have begun to find application in buckling load control (Ferrari and Sigmund, 2018;Chin and Kennedy, 2016;Stanford, 2017). Judicious choice of the aggregation function parameters allows one to obtain a smooth lower bound to λ 1 , provided an adequate number of buckling modes are computed.
Compounding the many modes issue, standard practice in topology optimization is to approximate voids via ersatz material on a fixed finite element mesh. The presence of such weak material can cause additional, spurious buckling modes to appear, whereby the weak material dominates the deflection response. An arbitrary number of such modes can be grouped with a near-zero buckling load, and since a finite number of modes is computed, this can cause the real modes to be lost from the analysis. For the SIMP method, whereby each finite element is assigned a density variable which dictates Young's modulus, a modified interpolation scheme has been proposed (Bendøse and Sigmund, 2003;Thomsen et al, 2018;Gao and Ma, 2015;Stanford and Beran, 2013), which advocates using different Young's modulus values in computing K and K s . The aforementioned level set study by Dunning et al (2016) utilized a similar approach, whereby a penalization function was applied to reduce the stress measured in void elements, thus increasing the associated eigenvalues, making them noncritical.
The presence of repeated eigenvalues (λ i = λ i+1 ) may also hamper optimization, since at such a point, the standard eigenvalue derivative formula breaks down (such will be shown by example in section 2). Additionally, the eigenvalues are no longer Fréchet-differentiable; this is due to the re-ordering of buckling modes that may occur from one iteration to the next. This typically presents as oscillations in convergence when two eigenvalues become close (Manickarajah et al, 1998). Fortunately, these issues are well-documented (Seyranian et al, 1994): as will be detailed in section 2, an alternative derivative formula can be used when repeated eigenvalues are detected, and should the order of eigenvalues be important to the chosen optimization method, the concept of generalized gradients may be adopted (Rong et al, 2001;Seyranian et al, 1994;Thomsen et al, 2018).
Summarizing the above, the importance of including buckling performance in topology optimization methods is well-recognized, and the literature documents a number of approaches for circumventing mode-switching, spurious buckling and eigenvalue multiplicity, thus making the endeavor possible. The authors are aware of only one study to date which implements buckling constraints in level set topology optimization, namely that by Dunning et al (2016), which demonstrated use on a bracket-like component. The aim of the present work is to apply the level set method to the buckling behavior of shell structures, for which design tends to be strongly buckling-driven when under compression. We utilize a sequential level set topology optimization method which is shown to naturally avoid the issues related to mode ordering, and is capable of simultaneouslyincreasing several buckling loads at once if required. We apply the method to two-thickness shell structures (where each point on the structure is assigned one of two thicknesses), as well as those with voids/cut-outs. Inclusion of such topological features is shown to introduce a trade-off between buckling load and structural stiffness; such a relationship has been noted recently by others using SIMP methods for 2D solid structures (Gao and Ma, 2015;Gao et al, 2017). The method is shown to be effective in several design applications, including flat plates, stiffened panels and ribs in an aircraft wingbox.

Level Set and Analysis Methods
Level set methods (LSM) provide a way to define the position and track the movement of geometric boundaries. Though first applied to the motion of fluid interfaces (Sethian, 1999;Osher and Fedkiw, 2003), this approach is more recently gaining popularity in structural design, specifically for topology optimization (Wang et al, 2003;Allaire et al, 2004;Dunning and Kim, 2015). In contrast to the SIMP methods, LSM can avoid having to analyze designs with significant regions of intermediate material properties (gray elements), generally at the expense of increasing the number of design iterations required, since change is limited to the existing boundaries. We utilize a scalar level set field φ, which is defined everywhere in a given design space, and through which the location of domain boundaries is implicitly defined: In the above, domains 1 and 2 can refer to solid and void respectively, which represents a typical topology optimization scenario. However as will be shown below, domain 2 can just as easily represent another material, the same material at a different thickness, etc. The relationship between φ and the resulting structure is depicted in Figs. 1a and 1b. In this work, we discretize the level set function on a two dimensional mesh comprising regular square elements, with bilinear interpolation within a given element.  Herein we assess the behavior of structures using the finite element (FE) method. We make use of 4 node flat shell elements, comprised of a Mindlin-Reissner plate element (Bathe and Dvorkin, 1985) combined with a plane stress element, using code developed in-house. In order to avoid re-meshing, we employ a fixed grid analysis, whereby the entire design domain is meshed once prior to optimization, and the physical (plate thickness) and material (Young's modulus, Poisson's ratio and mass density) properties of the elements are assigned in order to approximate the behavior of an arbitrarilyshaped structure. We achieve this by first assigning a variable w to each finite element, in the same fashion as density variables are assigned in the SIMP method (note that we avoid using ρ in order to avoid confusion with material density). This idea is shown in Fig. 1c.
In our method, w j is computed in such a way as to represent the volume fraction of finite element j which lies in domain 1, namely: where the above represent volume integrals over element j, and H (·) denotes the Heaviside operator. The discrete integration is done as for the finite element analysis, such that the ψ j vector contains the value of φ at the finite element nodes, and w j is a vector formed by Gaussian quadrature and used to integrate nodal values over the element volume. It should be noted that finite element nodes need not coincide with the nodes of the level set grid. In this case, a linear map is obtained such that where φ and ψ denote φ values on the level set grid and finite element mesh nodes, respectively. We node that multiple level sets can be utilized on a given structure, as will be shown in section 3. In order to maintain smoothness in regard to design derivatives, the analytical Heaviside operator is substituted for a polynomial approximation: In the above, ∆ represents the width of the smooth transition zone, and a non-zero β acts to shift the curve on the φ axis, which we will discuss uses for below. Once w j has been computed, the stiffness matrix for element j is found as: In the above, K i j represents the stiffness matrix of element j computed using the physical and material properties of domain i; E, ρ, ν and h represent the Young's modulus, mass density, Poisson's ratio and thickness, respectively. By observation, elements within the Heaviside transition zone (closest to the structural boundary) will have 0 < w j < 1, and will thus retain intermediate properties between domains 1 and 2. As such, this method produces a region of gray elements and is not truly binary; such is a consequence of using a fixed grid analysis method in lieu of a fitted mesh. The same calculation is performed with the element mass matrices, which are needed for structural mass calculations. The stress stiffness matrices, however, require special treatment. In the case where domains 1 and 2 both represent non-void structures, such as in the design of two-thickness plates, the stress stiffness matrix is computed in the same manner as Eq. (6). If domain 2 represents a void, such that E 2 → 0, h 2 → 0 or both, computing K s in the manner of Eq. (6) is known to produce spurious buckling modes. In these cases, we use the concept of stress relaxation (Bendøse and Sigmund, 2003), whereby If domain 2 is void: Otherwise: It is noted that the above strategy is conceptuallyequivalent to using the modified interpolation scheme proposed for the SIMP buckling methods (Bendøse and Sigmund, 2003;Thomsen et al, 2018;Gao and Ma, 2015). The effect of using the relaxed formulation is depicted in Fig. 2. For the structure shown, E 1 = 200 GPa, E 2 = 10 −6 GPa, ν 1 = ν 2 = 0.3, h 1 = 30 mm, h 2 = 0.1 mm. That is, domain 2 was given properties to approximate a void region. Simply supported conditions were enforced on all boundaries and a uniformlydistributed compressive load was applied to the vertical edges. As for all subsequent analyses, the static problem Eq. (1a) was solved using the HSL MA57 solver (HSL, 2002), and Eq. (1b) with ARPACK (Lehoucq et al, 1998). Without stress relaxation, the first several buckling modes are spurious, dominated by deformation in the void region. Including relaxation confines the behavior to that of the solid structure. As noted above, it is possible to shift the soft Heaviside function on the φ axis using the parameter β in Eq. (5), and Fig. 3 demonstrates the motivation for doing so. Namely, as the width of a structural member (d 1 in the inset of 3) reduces towards zero, the nonshifted Heaviside function with β = 0.0 produces a structural discontinuity, since H (φ → 0 − ) = 0.5. Shifting the curve with β = ∆ means H (φ → 0 − ) = 1.0, and a resulting smooth reduction in buckling load as structural members vanish. As the width of structural members increases, the buckling load computed using either formulations converges. The effect of this on optimization convergence will be shown in section 3.  with β = ∆ removes the discontinuity in buckling load as structural members vanish. The structure shown comprised a 2m × 2m plate discretized using 100 × 100 shell elements.

Problem Formulation and Derivatives
In each of the numerical examples given in Section 3, we will aim to maximize the first buckling mode λ 1 . As noted in Section 1, several modes are expected to change order or become simultaneously-active during the optimization; as such, objectives based solely on λ 1 are non-smooth. Also noted above, strategies to circumvent this issue include the use of aggregation functions and eigenvalue separation techniques. A third option, which we adopt here, is to employ a bound problem formulation, whereby maximization of λ 1 is replaced by that of an auxiliary variable λ aux , with constraints added such that said variable remains less than all of the computed eigenvalues (Bendøse and Sigmund, 2003): In the above, m and m * represent actual and target mass values, respectively; P represents a perimeter measure, added for regularization. The constant c p controls the strength of the regularization. For all the numerical examples presented in Section 3, we found that computing the first 10 buckling modes was sufficient to ensure optimization stability. That is, less than 10 modes were active during the optimization of all structures depicted in Section 3. Such a formulation retains the information related to all computed eigenvalues, though mode switching can still affect convergence if the order of eigenvalues is relevant to the chosen optimization algorithm (such as when performing line search, etc.). In such cases, mode tracking can be employed, which allows modes to cross without being re-ordered (Eldred et al, 1995). As will be explained in section 2.3, we herein utilize a sequential algorithm for which the order of eigenvalues is unimportant. In order to solve problem (8), we require the derivatives of m, P and λ i with respect to the topology, implicitly defined via φ. For m and λ i , we do this via the chain rule, first deriving with respect to w, then converting to φ via Eq. (3). The mass m is defined via the mass matrix: where the vector z contains 1 for deflection degrees of freedom along the z-axis direction and zero elsewhere (note that any translational direction could be used).
For element j we thus have: The derivative of λ i can be computed by pre-multiplying Eq. (1b) by the eigenvector v i , then differentiating: Note that [∇ u K s ] represents the gradient of the stress stiffness matrix with respect to the displacement vector.
In order to avoid computing the expensive ∂u ∂wj term in Eq. (12), we use the adjoint method. Pre-multiplying Eq. (1a) by an adjoint vector u a , then differentiating: Since Eq. (14) is zero, we can add it to Eq. (12) without changing the result. Doing this, and setting the terms in ∂u ∂wj to zero, we retrieve the standard eigenvalue sensitivity formula: where where we have made use of Eq.
(2). The term ∂f ∂wj is present when the structure is subjected to designdependent loading. Herein, some of our examples include both static and self-weight loading, such that f = f s + f w ; the latter can be computed as: where n g accounts for gravity (for example n g = −9.81 m/s 2 for 1g loading at sea level). As such, the force derivative is: When repeated eigenvalues are present (λ i = λ i+1 , etc.), Eq. (15) does not give a unique value for ∂λi ∂wj , since more than one eigenvector becomes associated with it. That is, one could choose to substitute either of the v i in (11) for v i+1 , and ultimately conclude a different value for ∂λi ∂wj . In such cases it has been shown that in place of Eq. (15), one should solve the following sub-eigenvalue problem (Seyranian et al, 1994): where A and B are n × n matrices, where n is the eigenvalue multiplicity, with element values as follows: where In practice, retrieved eigenvalues are unlikely to be exactly equal, and the switch from Eq. (15) to (19) can be made according to a pre-specified tolerance, for example when λ i and λ i+1 differ by less than 1%, and are considered close.
To demonstrate the behavior of Eqs. (15) and (19), consider the following simple test problem: System (23) has a repeated eigenvalue at w = 1, as shown in Fig. 4a. Note that mode tracking has been applied in order to color the lines. The solid and dashed lines in Fig. 4b depict the values of ∂λi ∂w computed using Eqs. (15) and (19) respectively; also shown are values obtained via central finite difference. As shown, Eq. (15) is appropriate for all values of w other than the exact point of multiplicity, where a discontinuity occurs. In fact, for this system w = 1 ± 10 −16 was far enough from the multiplicity that Eq. (15) agreed with the finite difference. The value obtained from Eq. (19) remains smooth through the multiplicity, though quickly diverges from Eq. (15), and the finite difference result, away from this point.
We note that while this simple system demonstrates a clear difference in the value of ∂λi ∂wj using Eqs. (15) and (19) in the presence of repeated eigenvalues, we were unable to produce such a difference using actual finite element models of structures. For reference, each of the numerical examples presented in Section 3 were run using both derivative formulations (15) and (19), and found to produce the same results to numerical precision. Based on the simple example given above, this is most likely due to eigenvalues differing by at least numerical precision at all times, which is enough to validate Eq. (15). It has been noted that many results and algorithms can be found in the literature which ignore the possibility of repeated eigenvalues altogether for this reason (Bendøse and Sigmund, 2003).  Derivatives with respect to w j are mapped to the level set field via Eq. (3), for example where φ k is the φ value at level set mesh node k. The perimeter P is estimated for element j as: where, as in the calculation of w j , the discrete integration is done as for the finite element analysis, such that the ψ j vector contains the value of φ at the finite element nodes, and p j is a vector formed by Gaussian quadrature and used to integrate nodal values over the element volume. We note that this is not the true perimeter -for an exact measure, the integrand should read (δ (φ) |∇φ|) -though the amendment greatly-simplifies the computations. As for the w j calculation, we use a smooth Delta function defined as the derivative of Eq. (5a). The perimeter derivative is then simply Lastly, for all the numerical examples given in section 3, the compliance of the structure will be measured, and as we will show, in some cases must also be incorporated into the optimization method. Compliance and its derivative is given as (Bendøse and Sigmund, 2003):

Optimization Algorithm
It is common-practice to maintain φ as a signed-distance function, such that |∇φ| = 1 everywhere in the computational domain; this ensures a well-behaved boundary, both in terms of the percentage of finite elements with intermediate w values and merging/splitting behavior of boundaries. In order to convert an arbitrary φ field to a signed distance function with the same boundary locations, we use a combination of the marching squares and fast marching algorithms (Osher and Fedkiw, 2003); the former locates a discrete set of boundary points x b , defined such that φ (x b ) = 0, while the latter propagates the |∇φ| = 1 property outwards to the remainder of the field. In order to preserve the signed distance property through design changes, it is apparent that the entire φ field is not free to change arbitrarily; φ should only change in ways that maintain the signed distance property. The fast velocity extension algorithm proposed by Adalsteinsson and Sethian (1999) details a method to propagate changes at the boundary ∆φ b to changes in the rest of the field ∆φ which achieves this preservation. The method detailed by Adalsteinsson and Sethian (1999) assumes that the values of change at the boundary nodes is known a priori, however since all the operations required by the method are linear on the level set grid, one simply needs to track which boundary nodes would contribute changes to the value of φ at each of the level set nodes; then a linear operator can be formed which allows the result of arbitrary changes to be computed. That is to say, the operator ∂φ ∂φ b can be computed in the following relation: such that |∇ (φ + ∆φ)| = 1. Computationally, both fast marching and fast velocity extension can be undertaken simultaneously. Eq. (29) allows the formulation of optimization problems in terms of changes at the boundary points only, which are subsequently propagated to the entire field. We note that since the fast velocity extension algorithm is only first order accurate, we also re-initialize the φ field using fast marching after each design update. It is noted that the above procedure differs from other reported level set topology optimization methods, which formulate equations in terms of a velocity normal to the boundary, which is propagated to the remainder of the field via the Hamilton-Jacobi equation. The difference arises due to the use of the smooth Heaviside to map structural properties, which makes it possible to work with changes in φ directly. In order to solve problem (8), our optimization method proceeds by solving sequential linear approximations to the objective and constraints, where the variables comprise the changes in φ value at the boundary points, ∆φ b : At each iteration, we utilize IPOPT (Wächter and Biegler, 2006) to solve the following sub-problem: where J 0 , m 0 and λ 0,i denote values at the current iteration. The value of γ denotes the move limit for a given boundary point φ value, and can be understood in the same way as a trust region in other SLP algorithms. In this work, we initially set γ to be half the width of a typical finite element in the given mesh, and update the value as follows: If the projected objective value -J 0 + ∂J ∂φ b T ∆φ b -differs by the real objective value by more than 1% for 10 sequential iterations, γ is halved; such facilitates convergence for SLP algorithms.
Although not done here, logic could be employed such that the move limit increases under certain conditions to facilitate faster convergence. By inspection of sub-problem (30), mode tracking need not be employed, since the optimizer is required to maintain all of the computed eigenvalues above a given value at the next iteration, regardless of their ordering. Instability could be encountered in cases where all of the computed modes are being actively-constrained; in such cases, modes which are not computed at a given iteration could switch places with one of the computed modes. In order to avoid this, a sufficient number of modes must be computed during the buckling analysis, such that not all are simultaneously-active. As noted above, we found 10 modes to be sufficient for all numerical examples presented in Section 3.
The optimization procedure can be summarized as follows: We begin by choosing an initial design and the accompanying φ field. A set of boundary points are then extracted using the marching squares algorithm, followed by the fast marching method to reinitialize φ as a signed distance function. At the same time, the relational matrix in Eq. (29) is computed. The w j are then assigned to all elements via Eq. (3), followed by finite element analysis Eq. (1) using the HSL MA57 and ARPACK solvers. Once complete, the objective function, constraint values and the associated derivatives are computed first with respect to w, then to φ as described above, and finally to the boundary points via Eq. (29). We then formulate the linearized sub-problem Eq. (30) and employ IPOPT (Wächter and Biegler, 2006) to select the optimal design change ∆φ b . Such is repeated until convergence, which herein is deemed to occur when all constraints are satisfied and the improvement in objective value is less than 0.1% for 10 sequential design iterations. This is shown diagrammatically in Fig. 5.

Numerical Examples
We begin this section by studying two-thickness flat plates, for which there are several benchmark designs available in the literature. We then proceed to study the topology optimization analogue; that is, plates with voids/cut-outs. Panel stiffeners are then explored, for which the SIMP method has been applied previously. Finally, we apply the method to the ribs in an aircraft wingbox.

Flat Plates
The first model comprises a 2 × 2 m flat steel plate (E = 200 GPa, ν = 0.3, ρ = 7850 kg/m 3 ), which is edge-loaded uniformly in the x−direction, with either simply supported or clamped conditions on all edges as shown in Fig. 6a. The finite element mesh comprised 100 × 100 square shell elements. Manickarajah et al (1998) applied the ESO topology optimization method to the buckling load maximization of such a model, whereby each point on the structure could take one of two pre-defined thicknesses h min or h max > 0. The authors also compared their results to theoretical optima identified previously in the literature, and this problem is thus considered to serve as an appropriate benchmark for the present work. Figs. 6b and 6c depict the first buckling modes for a uniform 15 mm thick steel plate with simply-supported and clamped boundary conditions, respectively. Note that subsequent optimized designs will be compared to these reference performance values. Fig. 7 depicts optimization results for the simply supported plate with h min = 10 mm, h max = 20 mm using problem formulation (8), with the target mass m * set equal to that of the uniform 15 mm thick reference design, and regularization parameter c p = 0.01. This small value for c p was found to have negligible effects on the optimization progression, though produced qualitatively smoother designs than were retrieved without its inclusion. As mentioned above, the first 10 buckling modes were computed during the finite element analysis Eq. (1b). As shown in Fig. 7d, all performance metrics were scaled by the constant reference values throughout the optimization in order to have magnitudes on the order of 1. The compliance C, as well as the second buckling load λ 2 are also included in the plot, though were not actively constrained in this case.
As per Fig. 7a, the initial design chosen comprised a regular 2 × 2 array of holes, with radii chosen in order to exactly satisfy the mass constraint. Figs. 7b and 7c depict the design progression through to the optimized design shown in Fig. 7e, which possesses a buckling load 39% higher than the reference, uniform thickness design in Fig. 6b. The optimized design bears a close resemblance to the analogous design obtained by Manickarajah et al (1998), which reported a 37% increase in buckling load. Those authors also compared the performance of such two-thickness plates to results reported in the literature for continuous thickness distributions; that is, where every point in the plate can take on h ∈ [h min , h max ], in contrast to the binary distributions depicted herein. The optima given by Levy (1996) and Pandey and Sherbourne (1992) perform only 29% and 24% better, respectively, than our reference design; it is thus apparent that binary thickness distributions may be advantageous compared to continuous distributions in regard to flat plate buckling loads. It is worth noting that this design problem was observed to be remarkably independent to the choice of initial design. Fig. 8 depicts the design progression for the above design problem when using different initial designs. Fig. 9 depicts the optimization for the same plate model with clamped boundary conditions on all sides. In contrast to the simply supported counterpart, this design becomes bimodal, with λ 1 ≈ λ 2 after iteration 10. As noted in section 2, the eigenvalues are not exactly equal, and in this case differed by approximately 0.01% from iteration 10 until convergence. As mentioned previously, this was enough separation to validate the use of Eq. (15) in computing ∂λ i /∂w, and indeed the results in Fig. 9 were obtained using Eq. (15)  accurate throughout. This is in contrast to the convergence history reported by Manickarajah et al (1998), which underwent marked oscillation once the buckling loads coalesced. The final designs obtained in that work, however, are again very similar when compared visually and performance-wise to that in Fig. 9e, which achieved a buckling load increase of 63% with respect to the reference, Fig. 6c.
As noted above, the convergence histories shown in Figs 7d and 9d included compliance, which was measured though not constrained as part of the optimization. By inspection, the value of compliance varied by  only 1-2% throughout the optimizations, even though the buckling load varied by as much as 63%. This trend does not continue, however, when topological change is introduced; that is, when the plates contain voids. Fig. 10 presents the static deflection and first buckling modes for two plates: The first is the uniform 15 mm thick plate from earlier; the second is a 30 mm thick plate with a single central hole, with radius chosen to give the same mass as the uniform plate. The void region was effected in this example, and all those subsequent, via ersatz material with a Young's modulus E void /E solid = 10 −6 . As shown, a vast (431%) increase in buckling load is easily obtained by introducing the circular void, though this accompanies a much higher (1561%) increase in compliance, and would most likely render the structure unusable. Such a trade-off relationship has been identified previously for 2D solid structures (Gao and Ma, 2015;Gao et al, 2017). For this reason, for the remainder of this section, we adopt a modified version of formulation (8), which includes the compliance measure in the optimization objective: where c b ∈ [0, 1]. Formulation (31) facilitates the study of Pareto trade-offs between compliance and buckling load. Prior to doing so, we note that numerical experiments indicated a much greater dependence of optimization on the initial design choice when topology is included. For example, Fig. 11 shows the performance and topology of optima obtained using formulation (31) with c b = 0.00, c p = 0.01; that is, minimum compliance designs, using different initial designs and simply supported boundary conditions. In general, the size of the features in the initial design dictate those in the final design, and smaller features facilitate stiffer optima. However, there are diminishing returns on the performance increase, and also likely a reduction in simulation accuracy, since for the designs with small features, only a few finite elements are present across the width of each structural member. In light of this, we undertook subsequent investigation of these plates using an initial 4 × 4 hole array, which strikes a compromise between performance and feature size.  As per Fig. 12a, the trend for both simply supported and clamped boundary conditions is similar: compliance rises modestly at lower buckling loads and exponentially at higher buckling loads. The minimum compliance designs for both boundary conditions (obtained with c b = 0.0) are unimodal, with the first two modes coalescing at higher buckling loads. Figs. 12b -12d and 12e -12g depict a selection of simply supported and clamped optima, respectively. Topologically, the minimum compliance designs (with c b = 0.0) comprise stiffening elements almost entirely in the load direction; cross-bracing members become prominent as c b > 0 and designs with higher buckling loads are favored by the optimizer. The maximum buckling load designs (with c b = 1.0) vastly reduce the volume of stiffening elements in the load direction, thus reducing the stress in the center of the plate, where buckling occurs; such increases the buckling load at the expense of stiffness. Note that the y-axis of Fig. 12a has been truncated for clarity; for reference, the maximum buckling load designs shown in Figs. 12d and 12g recorded C/C ref values of 12.3 and 14.0 respectively. Also shown in Figs. 12h -12j are the mode shapes for the clamped optima. The minimum compliance design is unimodal, and the mode shape takes the form of a half-sine wave. The design with c b = 0.3 is bimodal, whereby the half and full-sine waves coalesce. The maximum buckling load design has shed the half-sine mode entirely, and instead matches the full and one-and-a-half sine waves. Fig. 13 presents detailed optimization results for one of the simply supported optima, with c b = 0.2. Fig. 13e shows the convergence history, both with and without implementation of the soft Heaviside shift discussed in section 2. Without the Heaviside shift, the cross members are prone to snapping (see Figs. 13g and 13h), since as per the discussion on Fig. 3 this represents a discontinuity in the buckling load: The optimizer is unaware that buckling load will sharply decrease. Thus, two sharp drops in buckling load are observed in 13e, which ultimately leads to reduced buckling performance compared to the case where the shift is present.

Panel Stiffeners
We next present examples for the topology optimization of stiffening members on a blade-stiffened panel. Bedair (2009) reviews the importance of stiffened plates and shells in a vast array of engineering applications, including aircraft fuselages and wings, ship hulls, bridges and off-shore structures. The structural system comprises a face sheet/skin, reinforced by a series of stiffeners attached longitudinally (sometimes also orthogonally) with respect to the applied load, as shown in Fig. 14. Dimensions of the model studied herein were chosen to match those used by Stanford et al (2014), who undertook buckling load maximization using a SIMP method. As mentioned above, the first 10 buckling modes were computed during the finite element analysis Eq. (1b). Fig. 15 presents the results of one such optimization. The initial design was chosen such that the level set boundary traced the rectangular perimeter of the stiffeners. As such, the height and width of the stiffeners is free to vary, and topological features can form on the lower edge, where the stiffeners meet the face sheet. This is observed as the optimization progresses in Figs. 15a -15c. The motivation for this initial design choice, rather than introducing numerous holes in the body of the stiffeners, is that only 20 shell elements are present across the stiffener height, which is considered too few to support complicated topologies. The chosen initial design prevents the formation of small structural features, and is considered more likely to accurately model the stiffener behavior.
The convergence history Fig. 15d documents an initial reduction in buckling load which is caused by the optimizer satisfying the mass constraint. Once m = m * , the buckling load increases and becomes trimodal (with λ 1 = λ 2 = λ 3 ). These three active modes are depicted alongside the optimized design in Fig. 15e. Topologically, the final design comprises two solid outside stiffeners with reduced height compared to the reference, with the central stiffener becoming arched. Although 15% lighter than the uniform reference design, 93% of the buckling load is maintained. The topology reported by Stanford et al (2014) is markedly different to that shown here, with far more structural features present and resembling a truss. This notwithstanding, the same 93% buckling load performance was also reported by those authors. Such results reinforce our previous observations regarding the strong dependence of the optimization on the initial design choice, and demonstrate the non-convex nature of this problem.
To further explore this notion, we include optimization results of the stiffened panel using two alternative initial designs. The first is using the same FE mesh from above, with an initial level set boundary at the ends and top of the stiffeners, but not on the lower edge. Thus, the width and height of the stiffeners are free to change, though the arch-like features noted above are not able to form. The second is using an FE mesh with seven stiffeners instead of three. The results for both cases are presented in Fig. 16. As shown in Figs. 16a -16c and 16d -16f, both cases progress to the same final design, shown in Fig. 16h, and comprising three stiffeners, the central being slightly higher than those outer and all tapering noticeably at the load edges of the face sheet. Such demonstrates the ability of the proposed method to remove stiffening elements entirely, a feat not easily achieved using parametric formulations. In this case, the optimized design is bimodal, and performs significantly worse than the earlier design, this time retaining only 71% of the reference buckling load.

Wingbox Ribs
The last example we provide is for the topology optimization of ribs in an aircraft wingbox. The model, shown in Fig. 17, consists of an internal leading and trailing edge spar and ten ribs. Skin, in the shape of a NACA 0012 airfoil is fixed atop the ribs and spars, and four stringers span the length of the wing on both the upper and lower skins. The wing has a root chord of 1.0 m, aspect ratio of 3.0, taper ratio of 0.5 and leading edge sweep angle of 18.4 degrees. The stringer height is 10 mm. We note that this model was chosen quite arbitrarily; this was due to the lack of a true benchmark in the literature, as could be found for the flat plates and stiffened panels above. Whilst wingbox ribs have been subjected to topology optimization methods by others ( Stanford and Dunning, 2014;Krog et al, 2002), we are unaware of an example which was buckling driven. Recently, Stanford (2017) applied a SIMP method to wingbox stringers and included several constraints including buckling loads, though in that work ribs were not optimized. The finite element mesh for our chosen model comprised 69600 shell elements and for simplicity, the entire structure was assigned 2.0 mm thick aluminum (E = 70 GPa, ν = 0.3, ρ = 2800 kg/m 3 ). The root edge was fixed, and a static elliptic load was applied in the z−direction to every node in the skin elements as: where y is the y−coordinate of the node and f root z , is the load value at the root; since all optima are compared to reference values, this was set as 1.0 for sim-plicity. Whilst the above is a coarse approximation of aerodynamic wing loading, it is considered sufficient to demonstrate the method herein. It is noted that higherfidelity aerodynamic loads induce twisting behavior and may significantly alter the optimized designs reported herein. Self-weight was also present in this example, with n g = −9.81 in Eq. (17). We note that numerical experiments revealed the presence of self-weight to be largely inconsequential for this structure. In higherfidelity wing models, where engines, fuel tanks, etc. are present, as well as consideration given to vibration behavior and ground loads (e.g. a taxi bump), self-weight is expected to carry more significance. We retain it here for completeness. To undertake topology optimization on the ribs of this structure, we require ten level set functions; the grid aligned with rib 5 is present in Fig.  17a.   17c depicts the deflection of the reference design -where all ribs are present and uniform -under the applied elliptic and self-weight loading. Fig. 17e depicts the corresponding first buckling mode, which presents on the upper skin towards the root, where compressive stress is highest. Out of interest, we present in Figs. 17d and 17f the behavior of the structure when all ribs are removed. The static stiffness has increased by only 5%, however the absence of ribs results in a larger part of the upper skin buckling (a larger buckling wavelength), since it is no longer being broken into shorter panels; the buckling load consequently reduced by a significant 23%.
We note that in the detailed design stage, the individual panels, including the skin thickness and stringer geometry, would need to be optimized considering numerous loading conditions and many closely-spaced buckling modes are likely to be encountered; as will be shown below, our examples demonstrate well-spaced buckling modes. Whilst this may be considered non-representative of real structures, we include this example in order to demonstrate simultaneous optimization of numerous level sets (one per rib), as well as the potential to smoothly change the number of ribs while respecting the active buckling constraints. Fig. 18 presents the results of topology optimization using problem formulation (31) with both c b = 1.0 and c b = 0.0. As earlier, the regularization was constant at c p = 0.01. The mass constraint was set as m ≤ 0.96 m ref , where m ref = 30.8 kg was the weight of the reference design; such was chosen to force the optimizer to remove a significant portion of rib structure. As mentioned above, the first 10 buckling modes were computed during the finite element analysis Eq. (1b). The initial design was chosen such that level set boundaries were present at the leading and trailing edges of each rib. As for the stiffened panel results above, complicated topologies could not be supported on the reasonably coarse rib mesh, and the initial design was chosen to induce rib width changes instead of small-featured structures which were found, via numerical experiment, prone to form from initial designs with internal holes.
Topologically, the optimized design for c b = 1.0 in Fig. 18b (maximum buckling load) retains most of the rib structure near the root -enough to act as panel breakers, and increase the buckling load by 1% with respect to the reference design. Ribs 4 -7 have been completely removed, with structure returning in ribs 8 -10 at the wing tip. As for the stiffened panel example above, the convergence history Fig. 18a documents smooth performance changes, even when entire structural components disappear. The compliance of this design increased by 2% with respect to the reference.
In contrast, the optimized design for c b = 0.0 in Fig.  18c (minimize compliance) retains part of each rib (non were removed entirely). The buckling load dropped by only 5% with respect to the reference, due to ribs 1 -3 closest to the root remaining largely present, forcing a similar mode shape to the maximum buckling load design. Compliance increased by 1%, half the increase of the maximum buckling load design.

Conclusion
We herein documented and applied a level set topology optimization method to the buckling of shell structures. The method is equally-applicable to structures with binary thickness distributions, such as two-thickness plates, as well as those with void regions such as cutouts. In the latter case, consideration of static stiffness should be made, as a significant trade-off between buckling and compliance becomes evident. Stress relaxation was employed to effectively remove spurious buckling modes from void regions, allowing entire members to be removed entirely during optimization; such was demonstrated on panel stiffeners and wingbox ribs. A smooth, shifted Heaviside function was used in the definition of material properties, which was shown to avoid singularities associated with vanishing structural features, and facilities smooth optimization convergence. The sequential optimization method adopted avoids the need for mode-tracking, while allowing several buckling modes to be increased simultaneously.