A strategy for avoiding spurious localized buckling modes in topology optimization

We introduce a strategy preventing the occurrence of spurious modes in the spectrum computed by linearized buckling analysis in the context of topology optimization. Spurious buckling modes may appear in low density regions, a well‐known and largely discussed phenomenon. However, localized modes may also appear in solid areas, where stress concentrations occur. This second phenomenon, which is due to the inherent limitations of the linearized buckling analysis when used for complex stress states, is hardly addressed in the topology optimization literature. The remedy we propose makes use of elementary operations in the topology optimization framework: filtering and erosion, but now applied to the stress field. We show how this simple strategy helps mitigating the occurrence of spurious modes, in turn regularizing the optimization process towards high performance designs, which are then verified by nonlinear analysis.

and stability demands, which complicates the optimization task especially in bending-dominated configurations, 11,12 the most acknowledged difficulty is perhaps the onset of spurious buckling modes during the optimization.
Spurious modes in low-density regions are characterized by having a high amount of the overall strain energy located over design regions with low relative density. 13 This is a well-known issue, both for vibrations and buckling TO problems, and much research has been devoted to its cure. The most common remedy is to properly tailor the material interpolations for stiffness and inertia (for vibrations problems), or stress (for buckling problems). [14][15][16][17][18] Other proposed remedies are based on ad-hoc shifting strategies for the eigensolver, 11,13 or the removal of low-density elements from the computational grid. 19,20 Non-physical deformations happen also in geometrically nonlinear analysis, where can be solved using linear modeling for the low-density elements, 21,22 or when modelling contact by a "third-medium approach", 23 where can be solved by penalizing bending-like deformations. 24 A more subtle phenomenon is the appearance of very localized buckling modes in solid regions. These, which are by no means linked to the presence of intermediate densities, further localize as the mesh is refined, and are not observed in a nonlinear analysis performed with an appropriate hyperelastic material model. 25 Also, compared to the low-density ones, these localized modes are more difficult to characterize. Methods currently used for such identification generally look at the ratio of the average and maximum values of the modal displacement 26,27 (or of its spatial gradient 28 ), as this will suddenly drop when a mode becomes localized. In reality, these regions will reach a limit stress before the onset of geometrically-related instabilities, calling for the modeling of the material failure. Therefore, in the linearized modeling of stability these localized modes should also be treated as non-physical artifacts, originating from the lack of a length scale in a continuum model, and the resulting stress localization. Until a decade ago, this second phenomenon was not mentioned by most of the works on buckling TO, likely because it was overshadowed by the earlier-acknowledged, low-density modes. In recent years, several authors have encountered these localized modes in solid regions while solving buckling TO problems on very fine discretizations, starting from the 3D cantilever design obtained by Dunning et al., 8 then the intricate design subjected to shear analyzed by Ferrari and Sigmund, 12 and more recently also by Russ and Waisman. 26 We aim at retaining the LBA as a cheap analysis, useful for evaluating the stability of a design with well-defined structural features; thus we propose a strategy that overcomes its limitations when applied to continuum-like material distributions appearing in the early-stages of TO. The method consists of regularizing the stress field entering the buckling eigenvalue problem, by applying filtering and projection operations. Doing so, we simultaneously alleviate stress concentrations, and achieve a slight erosion of the stress field within the regions where the stiffness has developed a high value.
The mitigation of high stress gradients is due to the filtering operation, and this suffices to remove the highly localized buckling modes happening in the solids, or at least to shift them to higher eigenvalues. However, using stress filtering alone has some shortcomings when considering a design with intermediate densities, as spreading the stress field over regions with low stiffness may potentially worsen the phenomenon of low-density modes. Therefore, combining filtering and projection results in a much more stable and flexible strategy, simultaneously achieving: (1) the removal of highly localized buckling modes in the solid, due to the stress filtering, (2) the removal of spurious modes in the low-density regions, by eroding the stress field within the solid boundaries, and (3) mitigates the competition between stiffness and buckling demands in the early stages of the optimization, thus stabilizing convergence towards a well-defined structural configuration.
The remainder of the manuscript is organized as follows. Section 2 reviews the density-based TO formulation used, and in Section 3 we show the occurrence of spurious buckling modes in the solid, introducing a rational measure for their characterization. The removal of these modes from the computed set, by stress filtering, is demonstrated in Subsection 4.1. Then, Subsection 4.2 discusses the overall stress regularization idea, combining stress filtering and erosion. The proposed strategy is tested on the buckling strength maximization problem, considering two popular TO benchmark examples, and compared to the "standard" procedure. The optimized designs are then post-evaluated by considering geometrical and material nonlinearities, verifying their good performance, in Subsection 5.3. Final remarks are given in Section 6, and details about stress filtering, and its shortcomings when used alone in the TO process, are discussed in Appendix A.

TOPOLOGY OPTIMIZATION SETUP
We aim at best distributing a given amount of material within a design domain Ω, by using density-based TO. 18 The solid domain geometry is parametrized by the relative density field̂(x) ∈ [0, 1], linked to the design field (x) ∈ [0, 1] by the two operationŝ( Equation (2) is the so-called density filtering: 29,30 a linear convolution of (x) with the operator  r min , which is defined by the linearly decaying kernel k(x, ) = max{r min − ||x − || 2 , 0}, where r min > 0 is the filter radius. Then, (1) applies to the filtered field̃(x) the relaxed Heaviside function with inflection point at ∈ [0, 1] and curvature governed by ∈ [1, ∞). Hereafter, this will be referred to as the -projection, 31 and its effect is to push values̃(x) < closer to zero, and values̃(x) > closer to one, thus obtaining a sharper relative density field̂(x).
The material properties of the system depend on the relative densitŷ(x) through interpolation functions. 32 For the application at hand, we need to consider the interpolation of the stiffness, say i [̂(x)], and the one for the stresses, say i [̂(x)], whose specific form will be given in Section 5.

Discretization and definition of the response quantities
We consider a uniform discretization of Ω ≈ Ω h = ∪ m e=1 Ω e , consisting of m equi-sized elements Ω e , and a centroid-based discretization for ,̃and̂. Equations (1) and (2) read where  e is the set of points for which ||x i − x e || 2 ≤ r min and k ei = k(x e , x i ) = r min − ||x e − x i || 2 is the linearly decaying weighting. The volume fraction and its sensitivity are given by whereas the compliance due to a given load vector f, and its sensitivity read where u = K −1 (̂)f, is the pre-buckling equilibrium displacement. In the rightmost expression of (6) we used that̂e is locally defined on each element Ω e , and that the stiffness matrix K(̂) is assembled from the elemental contributions K e = i [̂e]K 0 , where K 0 is the constant matrix corresponding to the Young modulus E 0 and Poissons' ratio , of the solid material.
The non-dimensional stresses (i.e., strain combinations 33 ) are computed as where B is the linearized strain-displacement matrix. Equation (7) gives the physical stress only for̂e = 1 (i.e., on solid regions); otherwise, this must be obtained through interpolation We remark that also stresses are given a centroid-based discretization, and therefore (7) and (8) imply 0 e = 0 (x e ) and e = (x e ).
The element stress stiffness matrix is computed as where we assume that the same interpolation i is used for each stress component. The form of T[ 0 e ] is given in Reference 28, and we point out that each element matrix depends on local stresses.
The buckling load factors and the associated buckling modes ( i , i ) are computed by solving the generalized eigenvalue problem where i = 1∕ i can be seen as the "wavelength" associated with each buckling mode. The sensitivity of each i with respect to the relative density reads 12,34 is the adjoint vector. Finally, the gradients with respect to the design variables e are computed by using the inverse relationship to (2) and (1), namely where f may either be c, i or v f , and the expression for̃êe can be found in References 9,35.

CHARACTERIZATION OF THE LOCALIZED BUCKLING MODES IN SOLIDS
In this section, we demonstrate how spurious localized buckling modes appear in a simple solid column, and we introduce a measure of non-locality for their characterization.
We refer to the column shown in Figure 1A, clamped at the left edge and compressed from the right edge, with |f| = 0.01, E 0 = 1 and = 0.3. The domain is fully solid, thuŝ(x) = 1 everywhere. In the pre-buckling state, stress peaks are found at the clamped corners, as we can see from the von Mises stress distribution and its spatial gradient norm in Figure 1A. In order to enhance the visualization of spatial variations, these and other following plots are scaled as sign(f ) log 10 ( valid for both positive and negative f , and in each case we will list the value of c . In Figure 1B we plot the strain energy density (SED) for some of the buckling modes computed on a discretization of Ω h = 48 × 192. The lowest three modes are bending-dominated; then, interspersed with other bending modes, some very localized ones appear close to the loaded and clamped regions. Found on a fully solid design, this is clearly a different phenomenon than the classical one of spurious modes in low-density regions. 13,34 To characterize these spurious solid modes, we here get inspired by the mode-volume concept, commonly used in electro-magnetics and optics. 36,37 We consider the ratio between the average and maximum strain energy density (SED) where (x) stands for the system rigidity, ( (x)) is the strain associated with a given buckling mode, and dΩ is the smallest integration volume in the material domain. It is intuitive that (13) goes to zero as the field (x) = (x) ( (x)) 2 becomes more and more localized. In our discretized setup, (13) becomes as we have assumed the normalization ∫ Ω (x) ( (x)) 2 dΩ ≈ ∑ m e=1 T e K e e |Ω e | = 1. Equation (14) can be bounded as #e 1 ≥ [ ] ≥ #e 2 , where #e 1 is the number of elements where the SED field exceeds its mean value (see white contour lines in Figure 1B), and #e 2 is that corresponding to the so-called full width at half magnitude (see black contour line in Figure 1B). Most importantly, comparing Figure 1C with the SED plots shown in (B), we clearly see that highly localized modes are reflected by small values of [ ], meaning that high SED values occur over a few elements only.
For a few modes we have [ ] < 1 (see the below purple line in Figure 1C). This clearly proves their non-physical nature, as = 1∕ formally represents the mode "wavelength", which is then not captured on the current mesh, making the mode a numerical artifact. Moreover, whenever [ ] < 4, the corresponding mode would become an artifact if referred to the one-level coarser mesh, having elemental volume 4|Ω e |. Therefore, the condition [ ] < 4 can be used for identifying the mesh-dependency, and corresponding non-physical character, of these highly localized modes. This is shown in Figure 1D, for the modes computed on three finer grids. We acknowledge that, as the mesh is refined, more and more mesh-dependent modes appear, starting from lower eigenvalues. On Ω h = 12 × 48 we have two artificial modes, the first being 5 , whereas for finer discretizations, already 4 is artificial, and many more appear in the computed spectrum.

STRESS REGULARIZATION STRATEGY
We now describe the strategy for avoiding the spurious buckling modes characterized in the previous section. The idea is to replace the stress field used in the stress stiffness matrix G(̂, ) from the local one, computed by (7) and (8), with a regularized one, obtained by filtering and projection operations. We start by discussing the case when only linear filtering is applied to the stress field, in a similar fashion as density filtering (2), showing how this effectively removes the spurious modes in solid regions. Then, in Subsection 4.2 we discuss the whole stress erosion idea, which makes the regularization strategy more flexible and robust when applied to designs with intermediate densities.
We refer to Appendix A, for an in-depth discussion of some of the variants of the stress filtering, highlighting their merits and shortcomings, and their performance in the TO process.

Filtering of the stress field
We first consider filtering the original (non-physical) stress field 0 (x) according tõ where k(x, ) = max{r min − ||x − || 2 , 0} is the same as the kernel used for density filtering, now depending on the radius r min > 0, possibly different from r min . Equation (15)  The term (x) allows for a non-uniform spatial weighting of the filter, and choosing (x) as a function of the relative density will modify the kernel, unless we have a uniform design. On a fully solid design we can apply (15) to either 0 or , and setting (x) = 1 results in a smoothing of the physical stress field. However, if solid/void interfaces or intermediate densities are present, we must distinguish which quantity is filtered. Since non-physical stresses are large on low-density regions, applying (15) with (x) = 1 would increase the stresses on the solids, which must be avoided. Thus, we will use (x) =̂(x) to restore the correct, low weighting of the low-density regions in the filtering operation. We remark that this still allows the application of the material interpolation to recover the physical stresses (i.e.,̃= i [̂]̃0), since the result of (15) is still a non-dimensional stress, due to the kernel normalization.
Further details about the filter behaviour when using different, non-uniform weightings (x), both considering discrete 0/1 designs and designs with grayscales are given in Appendix A.
On the discretization grid, (15) is implemented as where i depends on the weighting function adopted, and  e is the set of points affecting element e, based on the radius r min . A PDE-based equivalent implementation is also possible, 38 perhaps allowing a more flexible choice of the boundary conditions. 39 However, our testing has not shown any relevant difference between the two forms, thus we do not discuss the PDE option any further. We now apply stress filtering to the uniform column discussed in the previous section (thus, = 0 ). Figure 2 shows the field distribution of the relative difference in the von Mises measure for the filtered and original stresses considering the two choices of boundary conditions for the filter operator: Neumann (N-BCs) and zero Dirichlect (D-BCs).
In the latter case, volume-preserving N-BCs must be maintained at the loaded and clamped boundaries, otherwise the modification of the stress induced by the zero D-BCs would violate the equilibrium. As a result of the filtering operation, high stress gradients are removed, whereas regions where the stress is nearly uniform are left untouched. The peak stress is reduced as r min is enlarged, and the D-BCs give a much higher stress reduction, as for this example the original stress peak was found near the boundaries. Stresses may slightly increase over the regions where they were originally low; however, just in the order of a few percent. Figure 3 shows how the filtering affects the results of the LBA. Already for the radius r min = 1 24 , many spurious modes are removed, and those with [ ] ∕4 < 1 are shifted to higher eigenvalues (i.e., 13 and 24 ). Then, for r min = 1 12 all spurious localized modes are purged from the spectral range of interest, and the surviving modes are bending waves with shorter and shorter wavelength (see Figure 3B,D). BLFs associated with localized modes are changed the most, whereas the low order ones increase a few % only (see Figure 3A). However, we should not choose a too large smoothing radius, as this would increase the magnitude of the BLFs too much.
Plot (c) displays only the BLFs associated with modes deemed as physical, that is, for which [ ] ∕4 > 1, based on Figure 1D. Taking the BLFs computed on Ω h = 12 × 48 as reference, we see that only a few of those computed on finer meshes survive as physical. However, these match closely together, which further proves their physical character, and the discrepancies are due to the larger approximation error on higher parts of the spectrum. 40 For example, the last two physical modes on the mesh Ω h = 192 × 768 correspond to the 18th and 21st originally computed (see Figure 1C), and therefore are affected by a higher approximation error. With the stress filtering, all the computed BLFs match almost exactly the physical ones computed on the coarsest grid without loosing accuracy, as they remain in lower regions of the spectrum.

Erosion of the stress field
We now combine linear and nonlinear filtering, the latter achieved by the -projection (1). Using the notation  ( , ) , to distinguish from the density projection, the whole stress regularization strategy is summarized in Figure 4. Starting from the same intermediate field,̃, we project two different relative density fields. The first, Then, setting > , we project another relative density field, , which is used for the interpolation of the physical stress,̂= i [̂( ) ]̃0, starting from the filtered, non-dimensional one.
To show the effect of this strategy, we consider the two intermediate designs of a cantilever optimized for minimum compliance, shown in Figure 5. We obtain̂( ) by using r min = 1 30 , = 0.5 and = 4, whereas for the stress filtering and projection we use r min = 1 30 , = 1 and = 4. The black continuous lines in the plots of Figure 5 bound the regions wherẽ e ≥ 0.75; using the SIMP interpolation with penalization p = 3, 32 this corresponds to E(̂e) ≥ 0.5E 0 . By looking at the maximum principal stress I , and its gradients' magnitude ||∇ x I (x)|| 2 , we can discuss the effect of the stress projection alone, and of the combination of filtering and projection.
We first look at the top row, referring to an early optimization stage, where we still have much grayscale in the design. The original stress field, computed from (8) usinĝ( ) , attains high values also on regions with low stiffness, potentially causing artificial buckling modes (cf. column (a)). If we apply the projection alone, thus computing = i [̂( ) ] 0 , all stresses outside the structural boundaries are projected to zero. However, the stress and stress gradient peaks are only mildly reduced, as they occur in a region with almost solid material (cf. column (b)); thus, spurious buckling modes are still likely to occur in the solids. The overall regularization strategy we propose combines filtering and projection. Filtering the non-dimensional stresses̃0 =  [ , 1 30 ] * 0 is responsible for reducing both the stress and stress gradient F I G U R E 4 (A) Stress erosion strategy adopted in the TO procedure. The field of design variables is filtered by the density filter; then subjected to two projections with different -values. The relative densitŷ( ) is used for interpolating the stiffness, whereaŝ( ) is used for interpolating physical stresses. (B) Influence of the parameters governing the erosion (r min ∕r min , ), setting = 4. We consider a 0/1 design variables field , where the interface is marked by the magenta continuous line. The dashed vertical lines mark the points where the filtered field̃goes to 0 and 1. The stress profile shown may be obtained with a uniformly compressed column (see Appendix A.1).

F I G U R E 5
Effect of the nonlinear stress filtering on an evolving design. In all the above plotŝ( ) is computed from the nonlinear projection (1), for = 1 and = 4. The colormap shows the distribution of the maximum principal stress I (top), and of its spatial gradient norm ||∇ x I || 2 (bottom). (A) Shows the original stress field, computed by (8), (B) the field obtained using the nonlinear projection alone, (C) the field obtained using filtering and then projection. The two rows, separated by a dotted line, refer to an early-stage and almost discrete design, respectively. Contour plots are scaled according to (12), for c = 4. peaks, as discussed in the previous section. Then, by using the projected density to interpolate the physical stress, = i [̂( ) ]̃0, all the stresses that have propagated outside the structural boundaries are projected to very low values (cf. column (c)).
The bottom row of Figure 5 shows the same comparison for a much later design stage, wherẽis close to a discrete 0/1 distribution. On such a design, the projection operation has little influence, and the noticeable effect is that of the stress filtering, reducing the stress concentrations over the solid domains and thus preventing spurious buckling modes in the solids.
We remark that the limit value = 1 was used here for illustration purposes. When applying the strategy in the TO process, we always choose ∈ (0, 1), which also makes the approach compatible with the PDE filter, as this never allows̃e = 1. The extent of the stress erosion clearly depends on the parameters r min , and , in combination with the corresponding parameters used for̂( ) . To summarize, we may refer to Figure 4B, showing the stress distribution across a 0/1 interface (marked by the vertical magenta line), where the original stress field is represented by the black continuous line. The number between parentheses in the legend correspond to the pair (r min ∕r min , ). We notice how, for a fixed , increasing the filter radius beyond r min does not help pushing the stress curve into the region wherẽ≥ 0.5. Also, using large values of r min leads to a too big reduction in the stresses, which is unwanted. A more practical choice is therefore to set r min ≤ r min , and to set a higher value for . A quantitative relationship between and the extent of the erosion, for a given radius r min , can be obtained from the analysis in References 41,42. Setting = 0.75 is sufficient to erode the stress field from roughly 10% of the elements within the solid boundaries, whereas for = 0.9, the eroded stress occupies only half of the material width.

NUMERICAL EXAMPLES
We focus on maximizing the fundamental BLF of a design, computed from the generalized eigenproblem (10), by minimizing the following function 12,43 (recall that 1 = 1∕ 1 ) which provides a smooth upper bound to 1 . In (18), KS is the aggregation parameter, governing the tightness of the bound, and  is the set of buckling modes included in the optimization. In order to check the occurrence of artificial modes within this set, we use the following criteria: • Low-density modes are characterized according to Gao and Ma, 13  • Highly localized solid modes are characterized with the help of the ratio in Equation (13). In this case, the modes for which [ ] ∕4 < 1 are considered artificial.
Given the volume fraction v f , the optimization problem reads where the link between and̂is given by (1) and (2), c opt is the compliance of the stiffest design with the same volume fraction, and c % ≥ 1 the allowed compliance increase. The compliance constraint may be inactive for the optimized design. However, we recall the inherent difficulty in reaching a well-defined structural configuration only with buckling and volume as driving criteria, especially when bending and shear effects are dominant. Therefore, we will always specify a compliance constraint, and show how the stress erosion described in Subsection 4.2 will promote its enforcement, helping the optimizer to proceed towards a well-defined design.
In the following subsections we consider two popular benchmark problems, comparing the results obtained by using the "original" approach, that is, using the original stress field, or the regularization strategy described in Section 4 for the LBA. In all cases, the interpolation functions for stiffness and stress used are where p is the penalization parameter and̂( ) e =̂( ) e for the original approach, whereas these two are distinct when using the strategy sketched in Figure 4. We use a continuation approach and start with p = 3, increasing it by Δp = 0.25 every 25 re-design steps, up to p = 6. Similarly, we start with = 2 and, after step 400, we increase it by Δ = 1 every 25 steps, up to = 16.
The design update is performed by the standard MMA routine, 44 and the parameters governing the asymptotes evolution are set as suggested in Reference 28.

Compressed column
We consider the compressed column sketched in Figure 6A, discretized by Ω h = 480 × 240, and we solve problem (19) for v f = 0.3 and c % = 2.5, starting from the minimum compliance design. Other dimensions, loads, boundary conditions, and optimization parameters are chosen as in Reference 28. In particular, the lowest 24 buckling loads and modes are included in the optimization, and we set KS = 150. Figure 6A displays the optimized design obtained for r min = 1 60 , and by using the original stress field in the LBA (i.e., without stress regularization). This design attains the BLF value 1 = 12.49 and compliance c = 1.37c opt , well below the imposed upper bound.
Starting from the minimum compliance design the optimization history is very smooth, and even if some of the higher modes are switching in the first iterations, this never affects the smooth increase of the lowest one (see Figure 6B). Low-density modes appear in four iterations only, and in the high part of the spectrum. Localized modes in solid regions happen more often, as we do not apply any stress filtering here; however, these very rarely affect one amongst the lowest five modes (see Figure 6C).
The optimization task becomes more challenging when starting from a uniform gray design, especially as we consider a loose compliance constraint. If we start from the uniform design with v f = 0.3 and set c % = 3, such that g C is violated in the beginning, the fulfillment of the compliance constraint initially drives the material towards a uniform column shape. Then, when g C is met, the design starts to increase its buckling resistance, reaching the BLF of 1 = 7.89, which is far lower than that of the design in Figure 6 (see Figure 7A). The compliance value is again well below the imposed bound (c = 1.162c opt ). Increasing the compliance upper bound to c % = 10 the optimization is never really driven by compliance, and this makes it very hard for the optimizer to navigate from the gray distribution to a well-defined design (see Figure 7B). Even if we reach a final design with high buckling resistance ( 1 = 12.05), the optimization history becomes very unstable, as grayscales are very persistent and cause the onset of low-density modes in most of the optimization iterations. From the rightmost plots of Figure 7, we clearly see that more and more spurious modes are developing in the spectrum as we start from a weak initial guess and we consider a loose compliance constraint. This both increases the computational cost, as we need to compute more eigenpairs to keep the same set of meaningful ones, and slows down the convergence due to the contribution of the spurious modes in the objective sensitivity. For instance, we point out that we needed to increase the number of modes included in (18) from 24 to 96, in order to converge to the design of Figure 7B.
On the other hand, even if using a tight compliance constraint is a way for regularizing the problem, this clearly shrinks the design space, and potentially drives the optimization towards a sub-optimal design from the beginning. This motivates the application of the regularization strategy discussed in Section 4, in the optimization process.

5.1.1
Use of the stress regularization strategy in the TO process Starting again from the same setting of Figure 7B  and =̂( ) . The relative densitŷ( ) e , used in the physical stress interpolation, is obtained with the projection parameters = 0.75 and = . All the other parameters are as before. The optimized design is shown in Figure 8A, and attains the BLF value, 1 = 12.94. Projecting the design obtained to a 0/1 distribution, and post-evaluating the buckling strength according to the original, local stresses, we obtain 1 = 13.51, which is almost 10% higher than the value attained by the design in Figure 6A.
This design is clearly efficient from the buckling perspective, consisting of two main struts braced by inner bars, and with a hierarchy of short linking members. Figure 8B,C, highlight the effect of the stress erosion strategy on the design and response evolution. The BLF starts from the very high value # 1 ≈ 18, as the design domain is initially filled with grayscales, and stresses are projected to zero almost everywhere. This quickly drives the design to activate the compliance constraint, and to develop a sparse truss structure. As the truss-like bars develop, and̃> , stresses are smoothly "activated" on them by the projection operation, and the buckling resistance of the bars enter the optimization. In this way, buckling influences the optimization in a more gradual way, simultaneously avoiding the competition of compliance and buckling and the appearance of the artificial modes, which would get the optimization process stuck, or to converge slowly. Figure 8E shows the computed spectrum, when using both the original and regularized̂stresses. For the design obtained from the optimization (i.e., having a non-discreteness measure m ND ≈ 0.22%), there is a gap in the high order eigenvalues corresponding to the two different stress distributions. However, for the design projected to a 0/1 distribution, this gap is significantly reduced, and it is entirely due to the stress filtering effect. In both cases, the lowest three BLFs are essentially not modified by the stress regularization. The performance of the design in Figure 8A has also been evaluated in COMSOL, 45 by using a body-fitted mesh. The critical load factor computed in COMSOL ( 1 = 14.19) is about 5% higher than that predicted in our fixed-grid framework (see Figure 8D), whereas the larger overestimation of the higher load factors is due to the better resolution of the boundaries given by the body-fitted mesh. Importantly, the buckling modes computed in the two frameworks match very well.
From the plot in (e), showing the trend of the parameter used for measuring the mode localization for the 0∕1 design, we see that all the modes are involving large portions of the domain. The two modes 16 and 17 , which attain the lowest value of [ ] ∕4, are also clearly physical, involving the buckling of one bar (see Figure 8A).

MBB beam
We consider the MBB beam discretized by Ω h = 480 × 160 elements, and a load of magnitude | f | = 2 ⋅ 10 −3 . Compared to the compressed column, this example is more challenging, as the large shear and bending stresses in the pre-buckling state introduce a tough competition between buckling and compliance. Also, this configuration is very sensitive to the The BLF for these two design is 1 = 0.438 and 0.615, respectively. Contour plots here are scaled according to (12) and c = 2. (B) plots the localization measure, based on the ratio given in (14). Here we also plot data corresponding to the design r min = 1 20 , for comparison. specific modeling of the boundary conditions. 46 To avoid stress concentrations, both the load and the support at the lower right corner are distributed over 16 nodes, and the support condition is enforced as a multi-point constraint through penalization, that is, prescribing the mean vertical displacement to be zero. 47 Figure 9 shows the compliance designs corresponding to v f = 0.35, obtained for the filter radii r min = 1 40 and 1 8 , and starting with p = 3 and = 2, then raising these up to 6 and 12, respectively, with a continuation scheme. The compliance and BLF of these two designs are (c, 1 ) = (8.955 ⋅ 10 −4 , 0.438) and (9.224 ⋅ 10 −4 , 0.615), respectively. To larger r min values correspond both higher compliance and BLFs; this is expected, since a design with thicker bars moves away from the stiffness optimal, and increases the buckling resistance of the individual members.
In both cases, there is a stress concentration in the upper left corner, near to the applied load. For the design with r min = 1 40 all the lowest 24 modes involve the deformation of one or more thin compressed bars. For the design with r min = 1 8 , the first six modes involve a physical buckling of the three compressed bars. Then, many artificial modes, localized in regions with stress concentrations, appear in the computed spectrum, as captured by the quantity [ ] ∕4 (see Figure 9C).
Using the compliance design with r min = 1 8 as initial guess, we now solve (19) allowing a compliance increase of 10% (thus, c % = 1.1), and with the density filter radii r min = { 1 20 , 1 40 , 1 80 } . By reducing r min , the optimizer should arrange the material to attain lower compliance, and have more freedom to improve buckling strength. The designs shown in Figure 10 are obtained using the original stress field in the LBA. In all three cases the compliance constraint is active and the buckling strength is significantly improved ( 1 = {3.805, 3.907, 3.925}). However, these design are very similar, and preserve the main features of the compliance one: (1) there are three main bars subjected to compression and three subjected to tension, and (2) the upper compressed bar runs almost horizontal for about 2/3 of the design domain. Compressed parts are just made thicker, and reinforcing bars appear in the upward triangular hole. However, the region near to the lower right support is weakly reinforced, and the compressed struts hardly develop any holes.
The middle row of Figure 10 shows the fundamental buckling mode corresponding to each design ( 1 ). As soon as the buckling mode shown in Figure 9B is prevented by the reinforcing bars in the central hole, the stress concentration close to the upper left corner becomes dominant, and the physical modes are shifted to higher eigenvalues. This causes the optimization process to stall. The bottom row shows the buckling mode ( # 1 ), computed in a post-processing phase using the filtered stresses̃. Removing the stress concentration by filtering, the fundamental buckling mode involves bending of the rightmost bar, with a slightly higher BLF. This further proves that we are in a local minimum, and the optimizer progresses are inhibited by a very local phenomenon. Figure 11 shows the designs obtained by using the stress erosion strategy, with r min = r min , =̂( ) , = 0.75 and = 4. These designs show an overall reinforcement of the structure, and their fundamental BLF ( 1 =  {3.936, 4.157, 4.428}), outperform those in Figure 10 by 3.5%, 6.4%, and 16.7%, respectively. In all three cases, the F I G U R E 10 Design obtained for maximum BLF, considering c % = 1.1, v f = 0.35 and three different density filter radii, starting from the compliance design of Figure 9B and using the original stresses in the LBA. The first row shows the final design, with contour plot of the maximum principal stress I . The second row shows the fundamental buckling mode corresponding to the original stress distribution ( ), which is the one used in the optimization. The third row shows the fundamental buckling mode computed, in a post-processing phase, using the filtered stresses (̃). compressed upper edge quickly tilts towards the supports, in a shallow arch-like configuration. The buckling of this upper arch is restricted by a truss-like set of bars, connected to the lower tension edge. This is a very efficient solution, compared with those in Figure 10, since making the upper edge more shallows both reduces the compression close to the support regions, and shortens the vertical bracing bars, thus avoiding the need to make them thicker. The bottom row of Figure 11 shows the designs obtained using the same set of parameters, but starting from the uniform gray design. The shear and bending-dominated pre-buckling stress state in the uniform material distribution complicates the task of the optimizer, as the trade-off between compliance and buckling strength makes the optimization task non trivial. Indeed, pushing too much on buckling strength from the beginning may drive towards a configuration with very poor pre-buckling stiffness, potentially leading to material disconnection. However, the designs in Figure 11 prove the ability of the stress regularization strategy to achieve a good design also when starting from a challenging situation. The BLFs reached here is even higher than before ( 1 = {4.157, 4.538, 4.825}), and we end up with designs having a slightly higher hierarchy, splitting some bars and introducing a few holes in the compressed struts.

Post-evaluation by geometric nonlinear analysis
The column designs are now post-evaluated considering geometric and material nonlinear modeling, in the commercial software COMSOL. 45 To avoid non-physical compressive stresses at large rotations, 25 we adopt the hyperelastic, neo-Hookean material law, and we trace the equilibrium path for the load parameter ∕ (LBA) 1 , normalized by the value predicted by the LBA. The critical load parameter is reached when the nonlinear iteration fails to converge, due to singularity of the tangent matrix.
At each equilibrium point u * , we also compute the nonlinear buckling load factor (NL)

1
, from the eigenvalue problem where G( * ) linearly depends on the current stresses, and K m (u * ) represents the material contribution to the current tangent stiffness matrix. 48 The results of the analysis, for the column design of Figure 8A and for the reference designs of Figure 6A, are summarized in Figure 12. The nonlinear force-displacement paths are plotted in log-log scale to enhance their variations, and the mean axial and transversal displacements at the loaded tip ( ⟨ u Tip ⟩ and ⟨ v Tip ⟩ ) are used as control deformations. For the reference design, the nonlinear critical point is attained for = 14.937, a value 16.6% higher than predicted by the F I G U R E 11 Design obtained for maximum BLF, considering c % = 1.1, v f = 0.35, and applying the stress erosion discussed in Section 4.
The two rows, above the dashed black line, refer to the compliance initial guess of Figure 9, and we reach 1 = 3.936, 4.157, and 4.428, respectively. The designs in the bottom row are obtained starting from a uniform gray initial guess and we reach 1 = 4.157, 4.538, and 4.825, respectively.
LBA. For the design obtained using the stress erosion, the critical point is = 16.801, 14.5% higher than predicted by the LBA. The black curves plotted in Figure 12 show the evolution of the ratio (NL)  Figure 12 also shows the deformed configuration in the last converged equilibrium point, and the fundamental buckling modes 1 predicted by (21), which are matching those predicted by the LBA. Therefore, the performance of the two columns, and the superior buckling strength of the design obtained by the stress erosion strategy, is confirmed also in the context of a full geometric and material nonlinear modeling.

FINAL REMARKS
We have presented a strategy for avoiding non-physical buckling modes when using linearized buckling analysis (LBA) for topology optimization. First, we introduced a measure, based on the mode-volume concept, 36,37 to characterize highly localized modes in solid regions. Then, we showed how these non-physical modes can be purged from the computed spectrum by filtering the stress field used in the LBA. Finally, by combining filtering and projection operations, we achieved the erosion of the stress field within the design regions that have already developed high stiffness. This prevents the onset of spurious buckling modes both in the solid and in the low density regions, and mitigates the competition between stiffness and buckling demands in the early stages of the optimization.
The numerical tests presented show the effectiveness of the method, and its ability to drive the optimization towards designs with higher buckling strength, even starting from a challenging situation, such as a uniform gray initial guess and a loose compliance constraint. In conclusion, the method makes the LBA more effective, for continuum TO, preserving its low computational cost, compared to the nonlinear buckling analysis. This latter, also accounting for material yielding, is a topic of ongoing research, and still poses several challenges.
The strategy can be easily extended to the 3D setting. In this context, it has been already observed that multigrid preconditioned iterative solvers, if not run to a fully accurate solution, provide an implicit smoothing of displacements F I G U R E 12 Evaluation of the nonlinear response for the reference column design (A) and for the design optimized considering the stress erosion strategy (B). The vertical load (F u ) is scaled by the factor ∕ (LBA) 1 , normalized by the linearly predicted BLF, and a perturbation F v = 10 −6 F u is applied at the tip. We consider to reach the critical point when the Newton solver fails to converge. Next to the load-displacement curves, we plot the deformation at the last converged step, and the fundamental buckling mode computed by the extrapolated eigenvalue problem (21). and stresses. 9 However, the method proposed here has more general validity, and the stress regularization is explicitly controlled by combining the filtering and projection parameters.

ACKNOWLEDGMENTS
This project is supported by the Villum Foundation through the Villum Investigator project "InnoTop".

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available on request from the corresponding author. The data are not publicly available due to privacy or ethical restrictions.

APPENDIX A. FURTHER DETAILS ON SOME VARIANTS OF THE STRESS FILTERING
We elaborate on the effect of some different choices for the kernel weighting in the stress filtering. An approach entirely relying on stress filtering has shortcomings when applied to more general designs than the uniform one: (1) a uniform weight increases the stress on low-density regions, potentially making the low-density modes phenomenon worse; (2) to avoid this, we may modify the filter operator. However, for designs with intermediate densities such modifications make the stress field non-smooth across the design boundaries. This motivates the projection discussed in Subsection 4.2 to achieve stress erosion. We also show a compressed column design, obtained with stress filtering alone, and achieving a sub-optimal solution, compared with the design in Figure 8.

A.1 Testing of different weightings on the stress filtering operation
The factor (x) in Equation (15) can be used to apply non-uniform stress filtering. For instance, by setting =̂the stresses on high density regions have more influence than those on low density ones, whereas = 1 −̂has the opposite effect. On a discrete design (i.e.,̂= {0, 1}), these choices correspond to dilation and erosion filters, excluding voids and solid regions from the kernel, respectively (see Figure A1).
On a design with intermediate densities, the erosion effect is achieved by ( ) = | , where and is a user-defined threshold. In this way, only points wherê> contribute to the filtering. Equation (A1) can also be used to modify the outcome of the filtering operation as such that the stresses are not modified over regions wherê(x) < . By combining ( ) = | and (A2), we apply a stress redistribution within the high density regions: stresses get smoothed on high density regions only, by a filter kernel involving high density regions only (see Figure A1D).

F I G U R E A1
Influence of some choices of the weight on the filter kernel, for a discrete design. The evaluation point x is marked by a circle, whereas the points in the neighborhood ∈  (x) are marked by squares. An empty circle means that the evaluation point itself is not considered in the filter kernel. We have three situations: both x and  (x) lie in the solid region (black); x is in the solid region and  (x) is partly in the solid and partly in the void region (blue); x is in the void region and  (x) is partly in the solid and partly in the void region (red).

F I G U R E A2
Effect of the stress filtering on the 0∕1 compressed column design. The 0∕1 interface is marked by the continuous magenta line in plot (A), and by the grey dotted line in plots (B,C). (A) Shows the von Mises stress for the original field and for the filtered fields̃, for different choices of the kernel weight . The contour plots are scaled according to (12) and c = 2. For each case, (B,C) show the cross section stress distributions at the mid-length and near the tip, respectively. Only half section is shown, due to symmetry.
Let us consider a compressed column, discretized by Ω h = 240 × 120 elements. First, we consider a discrete design, wherêe = 1 only over the central strip of 40 elements. The load is applied on the 24 elements symmetric to the midpoint, at the tip. We use N-BCs when applying (15); however, the solid-void interface is an internal boundary, not modeled in the filter operator, and therefore the boundary conditions influence only the loaded and clamped ends.
The effect of different choices of the weighting , considering r min = 1 15 , is summarized in Figure A2. The maximum stress, attained away from the solid-void interface, is reduced of about 4%, in all the cases. Figure A2B shows the cross section stress distribution at the mid span (i.e., at e x = 120). The uniform weight = 1 reduces the stress inside the solid boundary, but propagates it on the void region; this effect is exacerbated by the weight =̂, giving a stress propagation within a reach of r min in the void region. The weight = 1 −̂has the opposite effect, halving the stress within a reach of r min inside the solid boundary. The last two contour plots in (a) and the two bottom plots in (b,c) show the effects of the filter modification (A2), either keeping the kernel weight = 1, or setting = |0.5 . In the first case (cyan lines) we obtain the same result as with uniform filtering, inside the solid boundary, whereas the stress on the voids is zero. The second case (purple lines) corresponds to stress re-distribution in the solids, and has no effect on regions of uniform stress, whereas stress gradients are reduced near the tip. We point out that these two latter approaches introduce a non-smooth point in the stress field across the boundary.

F I G U R E A3
Effect of stress filtering on the column design with intermediate densities. The information shown in the three plots are the same as for Figure A2, but now (A) displays the relative difference between the von Mises measures of the filtered and original stress fields (see Equation 17). The contour plots are scaled according to (12) and c = 2. For these plots, we have set = 0.5 when using the filtering strategies (A1) and (A2).

F I G U R E A4
Post-evaluation of the design in Figure 6A, reduced to the 0/1̂distribution. (A) Shows the BLFs computed according to the original stress field, and to the filtered one, for r min = 1 60 and comparing some different kernel weights. The bottom plot shows the corresponding values of the localization measure, defined based on (13). Filtering on the solid purges all artificial modes, whereas the stress redistribution (purple line) still retains two of them, on the high part of the spectrum. (B) displays some representative buckling modes corresponding to the original stress field and to the one filtered using = 1 and the filter modification (A2). The colormap is scaled according to (12) and c = 4.

F I G U R E A5 (A)
Optimized design for the compressed column, obtained starting from a uniform gray material distribution and the stress filtering. The corresponding buckling modes are computed on the discrete 0/1 design, and using the local stress field. (B) Shows the evolution of the objective and constraints in the optimization progresses, (C) and (D) summarize the post-processing results, performed either in our framework and by using COMSOL.
For the same configuration, we now introduce a transition between solid and void by applying density filtering (2) with r min = 1 15 . Now all the choices of the filter weight amount to a stress increase somewhere in the intermediate density region (see Figure A3). Moreover, setting =̂or = 1 −̂does not correspond to dilation and erosion operators anymore, and the filter modification of (A2) may even introduce a stress discontinuity in the transition region, at points where >̂.
We now test the effect of these stress filtering alternatives on the modes and BLF computed by the LBA. We consider the design of Figure 6A, projected to a discrete 0/1̂distribution, and evaluate its buckling response for to the original and filtered stress fields, for r min = r min = 1 60 . The results are summarized in Figure A4. For = 1 (and for =̂as well), the stresses propagate outside the solid boundaries, leading to the onset of hundreds of buckling spurious modes in void regions, associated with zero BLFs. The weight = 1 −̂avoids this issue, and also removes the spurious solid modes. However, all the BLFs are substantially increased and one may argue against the overestimation of the lowest one (see green curve in Figure A4A). Filtering with = 1 and using the stress modification (A2) seems the best compromise, as the localized modes are effectively removed, and the lowest BLFs are not significantly changed (cf. cyan lines in Figure A4A and modes in Figure A4B). The stress re-distribution, achieved by combining (A2) with = |0.5 gives a similar result; however, some local modes survive at high frequencies, because excluding the voids from the kernel weight gives a smaller smoothing effect along the boundaries (cf. purple lines in Figure A4A). Figure A5 summarizes the outcomes of the TO process, when using stress filtering alone. The setup and all the parameters are the same as given in Subsection 5.1, and for the stress filtering we set r min = r min = 1 60 . The kernel weight is =f or the first 400 iterations. Then, when is increased, we switch to the filter modification (A2), setting = 0.75. This is motivated by the previous discussion, since the choice of the threshold may be arbitrary in the early optimization stages, whereas, once the main structural features have emerged, we need to avoid the back propagation of stresses on low-density regions.

A.2 Performance of stress filtering in the TO procedure
The buckling loads start from a high value due to the stress reduction caused by the filtering, and all the eigenvalues evolve smoothly, apart from small jumps in correspondence of continuation steps (see Figure A5B). In the end, the BLF attains # 1 = 8.80 according to the filtered stresses, whereas using the local stresses we have 1 = 8.36, about 5.2% lower. When post-evaluating the pure 0∕1 design, we obtain 1 = 9.76 and the BLF corresponding to the original and filtered stresses essentially coincide (see Figure A5C). The buckling modes computed by COMSOL, using a body-fitted mesh, match those computed in our framework, and the critical load 1 = 10.06 is about 7% higher. Moreover, from Figure 8D all the modes qualify as physical, and the ones that look localized still involve the failure of single, well-formed structural members.
The optimized design consist of two main struts, and a hierarchy of inner bars providing reinforcement against lateral bending. A narrowing of the struts at about 1/3 of the height resembles a hinged support with the clear effect of shortening the struts' inflection length, as we can see from the fundamental buckling mode (see Figure A5A). Despite its sound structural configuration, this design shows poor stability compared to those of Figures 6 and 7B, due to the lack of inner bracing elements, which would substantially increase the lateral stability of the struts. The inhibition in the development of inner bars is a downside of the stress filtering; in the early optimization stages, the smearing of the stress field caused by the filtering prevents the formation of bars close to the centerline, thus forcing the optimizer to take another path.