Bi-directional evolutionary structural optimization with buckling constraints

Buckling is a critical phenomenon in structural members under compression, which could cause catastrophic failure of a structure. To increase the buckling resistance in structural design, a novel topology optimization approach based on the bi-directional evolutionary structural optimization (BESO) method is proposed in this study with the consideration of buckling constraints. The BESO method benefits from using only two discrete statuses (solid and void) for design variables, thereby alleviating numerical issues associated with pseudo buckling modes. The Kreisselmeier-Steinhauser aggregation function is introduced to aggregate multiple buckling constraints into a differentiable one. An augmented Lagrangian multiplier is developed to integrate buckling constraints into the objective function to ensure computational stability. Besides, a modified design variable update scheme is proposed to control the evolutionary rate after the target volume fraction is reached. Four topology optimization design examples are investigated to demonstrate the effectiveness of the buckling-constrained BESO method. The numerical results show that the developed optimization algorithm with buckling constraints can significantly improve structural stability with a slight increase in compliance.


Introduction
In the past decades, the development of topology optimization has been rapidly progressing. Structures with high performance and minimum material utilization can be achieved by using topology optimization. Most of the early works on topology optimization focused on maximizing structural stiffness (Bendsøe and Kikuchi 1988;Rozvany et al. 1992;Xie and Steven 1993). Later, other effects have also been considered by researchers, such as natural frequencies , multiple materials (Huang and Xie 2009;Li and Xie 2021), structural complexity (He et al. 2020(He et al. , 2022, additive manufacturability (Bi et al. 2020(Bi et al. , 2022, and principal stress (Amir 2017;Chen et al. 2021). In addition, topology optimization under buckling constraints is regarded as a challenging and important topic. The neglect of the buckling effect in topology optimization may lead to poor stability and unexpected failure of a structure with optimum stiffness.
At present, the solid isotropic material with penalization (SIMP) method has been commonly used in buckling-constrained topology optimization. A typical issue of the SIMP material model is the appearance of pseudo buckling modes in low-density regions, which would generate incorrect sensitivities and poor convergence in the optimization process (Bruyneel et al. 2008). Many researchers have attempted to alleviate the difficulties in the eigenvalue-related optimization caused by the pseudo-buckling modes. For example, Neves et al. (1995) proposed to ignore the stress stiffness of the low-density elements. However, Kemmler et al. (2005) pointed out that such a cut-off method could lead to nondifferentiable functions. Bendsøe and Sigmund (2004) suggested a smooth version of this approach by employing two different material interpolation schemes for stiffness and stress stiffness matrices. This approach has been proven to be effective in eliminating pseudo buckling modes in many studies (Munk et al. 2017;Yi et al. 2019;Dalklint et al. 2021;Ferrari et al. 2021). Besides, Gao and Ma (2015) proposed a method that combined the eigenvalue shift and pseudo mode identification to avoid the influence of pseudo buckling modes.
Another challenge in buckling-constrained topology optimization is the complexity in sensitivity analysis of buckling load factors (BLFs). Buckling analysis based on a non-linear equilibrium equation was adopted in a few studies for topology optimization (Rahmatalla and Swan 2003;Lindgaard and Lund 2010;Lindgaard and Dahl 2013). Despite the fact that non-linear analysis provides an accurate description of the buckling phenomenon, the high computational cost hinders its widespread utilization in topology optimization. To improve the computational efficiency of nonlinear analysis, Pedersen and Pedersen (2018) suggested the use of non-incremental analysis and proposed a simple sensitivity analysis along with recursive redesign. In contrast, linear buckling analysis has been commonly employed in topology optimization due to its ease of implementation and acceptable computational cost (Ferrari and Sigmund 2019). The stress stiffness matrix in the buckling sensitivity is a function of design variables and displacement through the stress field . To reduce computational costs, Bruyneel et al. (2008) proposed to carry out the sensitivity analysis in a simplified way by neglecting the stress state variations. However, Luo and Tong (2015) indicated that by doing so, the sensitivity analysis of buckling constraints was treated in the same way as that of natural frequencies, which could cause errors in the optimization process. To avoid this issue, sensitivity analysis including the stress variations was used in more recent studies (Gao and Ma 2015;Ferrari and Sigmund 2019;Gao et al. 2020;Zhang et al. 2022). To overcome the difficulties in setting up buckling analysis, Ferrari et al. (2021) presented a 250-line MATLAB code that significantly reduced computational costs. Besides, BLFs may be non-differentiable when mode switching and repeat eigenvalues exist in the optimization process (Seyranian et al. 1994). Thus, a large number of buckling modes need to be considered in the optimization problem (Bruyneel et al. 2008). To reduce the multiple constraints into a single and differentiable one,  took the average of the sensitivities of the multiple eigenvalues with similar values. In addition, the Kreisselmeier-Steinhauser (KS) aggregation function (Kreisselmeier and Steinhauser 1979) was employed in many studies to approximate the maximum buckling constraint (Gao et al. 2020;Ferrari et al. 2021;Zhang et al. 2022).
The method of moving asymptotes (MMA) algorithm has been commonly used to solve buckling-constrained problems (Lindgaard and Dahl 2013;Gao and Ma 2015). Besides, Ferrari et al. (2021) employed an optimizer called 'ocUpdate' in which the buckling constraints were incorporated into the objective function using a Lagrangian multiplier determined by the bisection method. Alternatively, the augmented Lagrange multiplier has an additional penalty term, that is, the augmentation, which can be adjusted to control the contribution of the constraints to the overall sensitivity (Otomori et al. 2015;Wei et al. 2018). To the best of the authors' knowledge, no augmented Lagrangian multipliers have been developed for buckling-constrained topology optimization.
This study proposes a novel methodology based on bidirectional evolutionary structural optimization (BESO) for buckling-constrained topology optimization to alleviate the impact of pseudo buckling modes. The BESO method benefits from using two statuses (void and solid) for the design variables, thus avoiding the numerical problems associated with intermediate densities in buckling-constrained optimization. The KS function is employed to aggregate multiple buckling constraints into a single one. Furthermore, an augmented Lagrangian multiplier is developed for the consideration of convergence. Besides, one obstacle to the BESO method is that there is no control over the evolutionary rate once the volume fraction is reached. Therefore, a modified design variable update scheme is proposed to solve this problem.
The remainder of the paper is structured as follows. Formulas and procedures of the buckling-constrained BESO method are introduced in Sect. 2, followed by four numerical examples in Sect. 3 to demonstrate the effectiveness of the proposed approach. Concluding remarks are given in Sect. 4.

Buckling-constrained BESO method
In this section, the proposed buckling-constrained BESO method is introduced. In this method, an augmented Lagrangian multiplier is developed to combine the objective function and buckling constraints. The KS aggregation function is used to approximate the maximum buckling constraint so as to reduce multiple buckling constraints into a single one. Besides, a modified design variable update scheme is proposed to control the evolutionary rate after the volume fraction is reached.

Problem statement
The mathematical model for minimizing compliance subjected to buckling constraints is stated as where C is the mean compliance, K and u are the global stiffness matrix and the displacement vector, respectively. F is the load applied to the structure. v i is the volume of the (1) corresponding element and V * is the predefined volume fraction. 1 is the lowest BLF and − is the prescribed buckling constraint. x i is the binary design variable of an individual element, which equals to either x min (void) or 1 (solid).

Sensitivity analysis
In this study, two distinct material interpolation schemes are employed for stiffness and stress stiffness matrices to restrain pseudo-buckling modes (Bendsøe and Sigmund 2007;, which are given by Eq. (2).
where p is the penalty exponent, and E is the Young's modulus of the solid element. The sensitivity of the compliance with respect to the change in the design variable can be expressed as Linear buckling analysis is defined as where K and G are the stiffness and stress stiffness matrices, respectively, and j and j are the jth BLF and the corresponding buckling mode vector, respectively. To maximize the BLF j , its inverse value j shall be minimized, thus Eq. (4) can be restated as j is the eigenvalue of the above equation, therefore, the minimum BLF 1 is equal to 1/max j = 1/ 1 , where The eigenvectors j fulfil the orthonormalization condition: where jk is Kronecker delta. Then the maximum buckling constraint is approximated by the KS aggregation function (Kreisselmeier and Steinhauser 1979) where q is the number of computed eigenvalues. The number of buckling modes to be considered is dependent on the tendency to eigenvalue coalescing of each case. In this study, the lowest 12 buckling modes are computed, which (2) � seem to be sufficient for the problems considered. is the aggregation factor. In this paper, is set to 160, which has been shown to be appropriate for buckling optimization problems considered by Ferrari et al. (2021). To be specific, the numerical stability of the KS function is not negatively influenced by such a high value, and a good approximation of constrained aggregations can be achieved.
The approximated buckling load factor sensitivity, i.e., the derivative form of Eq. (8), is written as (Raspanti et al. 2000) where the buckling load factor sensitivity is calculated as The second term on the right side is the so-called adjoint term, and j is the jth adjoint vector, which can be solved by The sensitivity of buckling constraints [Eq. (9)] is added to the sensitivity of compliance [Eq. (3)] via an augmented Lagrangian multiplier Λ . Thus, the overall sensitivity can be found as Therefore, the sensitivity of element i in the BESO technique can be written as The augmented Lagrangian multiplier Λ is determined as where k is the additional penalty term (the augmentation) of the Lagrangian multiplier in the k th iteration of the optimization. is the status of buckling constraints. Specifically, = 0 indicates that the buckling constraints have never been activated; otherwise, is switched to 1 once buckling constraints have been activated during the optimization process. In this way, the continuity of the multiplier can be ensured. k is updated using the following scheme.
where is a predefined update coefficient. 0 = 1 and = 1.03 are set as default values, which have been numerically tested as suitable for most cases. A smalle r or larger value of can be selected to decrease or increase the update step. is calculated as = − 1 . Δ 1 is defined as the normalized value of the change of 1 over the last two iterations, i.e., is an update threshold and typically = 10 −3 is adopted. error ≤ is the convergence criterion that is introduced in Sect. 2.4. To improve computational efficiency and achieve convergence, two criteria and schemes are devised in this study for updating k . When the difference between 1 and − is large , the update criterion is tightened ( Δ 1 ≤ and error ≤ ), and the increment of k becomes smaller.
The aforementioned sensitivity members shall be averaged by a filter scheme to enforce a mesh-independent solution (Sigmund and Petersson 1998). The filter scheme is determined by generating a circular sub-domain Ω i centred at the centroid of the ith element with a radius of r min . Nodes and elements inside the sub-domain Ω i are used to calculate the averaged sensitivity number of the jth element ̃ j , as given by Eq. (16). The same method was also used by (Huang and Xie 2007).
where N denotes the number of the neighbouring elements.
ij is the linear weight factor related to the distance ( r ij ) between the centroids of the centred element j and its neighbouring elements i where r min is the filter radius. The sensitivity number is averaged by the previous three historical information for stability and convergence concerns. This process is expressed by where k is the current iteration number. Compared to the conventional BESO method, where the sensitivity number is averaged by the previous two historical information , the use of previous three historical sensitivity values in this study could enhance the stability during the optimization process. A further study of the dependence of the optimization history on the number of averaged steps can be systematically investigated in future research.

Design variable update scheme
The target volume for the next iteration can be determined by adding or deleting elements as where ER is named the evolutionary volume ratio. Once the volume constraint is satisfied, the volume will be kept constant for the subsequent iterations as Then solid elements are switched from 1 to x min if the following criterion is met Void elements are switched from x min to 1 if the following criterion is met where th is the threshold of the sensitivity number that can be determined by the bisection algorithm . However, it should be noted that there is no control over ER once the volume fraction is reached. This issue can lead to problematic designs in constrained optimization problems when starting from a design domain that has already reached the volume fraction. Therefore, a modified scheme is proposed in this study for determining the update threshold, which includes three steps as follows.
1. Calculate the percentage of the solid elements ( s ) that have been changed to voids by using the bisection algorithm (Huang and Xie 2010). 2. If the volume does not reach the target volume, or s ≤ ER , skip step 3. Otherwise, the threshold sensitivity numbers for removing ( del th ) and adding ( add th ) elements are recalculated in step 3. 3. Sort the sensitivity numbers of solid elements from low to high. del th is the sensitivity number of the mth element ranked in the sorted array, where m is calculated by multiplying ER and the number of solid elements. Similarly, add th is the sensitivity number of the mth element ranked in the sorted sensitivity numbers of void elements (from high to low). For example, there are 1000 solid elements and 1000 void elements in the design domain, and it is assumed that ER = 1% . The sensitivity numbers of solid elements are sorted as s The sensitivity numbers of void elements are sorted as v 1 > v 2 ⋯ > v 1000 , then add th = v 10 . Finally, the solid elements with sensitivity numbers ≤ del th are switched from 1 to x min , and the void elements with sensitivity numbers ≥ add th are switched from x min to 1.

Convergence criterion
The optimization process stops when buckling constraints, volume constraint and convergence criterion are all satisfied. The convergence criterion is defined in terms of the change in the objective function as where k is the current iteration number, is an allowable convergence tolerance, and M is an integer number. In this study, M = 5 is employed, denoting that the change in the mean compliance over the last 10 iterations is sufficiently small, e.g., = 10 −3 to 10 −5 .

Overall procedure
The evolutionary iteration procedure of the presented buckling-constrained BESO method is given as follows.
Step 1. Discretise the design domain, define boundary conditions and assign initial parameters, including volume fraction V * , evolutionary ratio ER , material penalty exponent p , filter radius r min , buckling constraint , and update coefficient ; Step 2. Obtain displacement U and compliance C by static analysis; Step 3. Obtain the BLFs j and buckling mode vector j by linear buckling analysis; Step 4. Determine the augmented Lagrangian multiplier; Step 5. Calculate the sensitivity numbers of all elements; Step 6. Improve the sensitivity numbers using the filter scheme; Step 7. Stabilize the evolutionary process by averaging the sensitivity numbers of the past three iterations; Step 8. Update design variables by using the proposed update scheme; Step 9. Return to Step 2 if the volume constraint, buckling constraints, or the convergence criterion is not satisfied. (23)

Numerical examples
In this section, four numerical examples are investigated to illustrate the effectiveness of the proposed algorithm. In all examples, the parameters are set as E 1 = 1 , = 0.3 , p = 3 , and x min = 10 −6 .

Compressed column
This example aims to maximize the buckling resistance of a compressed column. As shown in Fig. 1a, the rectangular design domain is discretized by 240 × 480 elements. Plane stress elements are employed, and the thickness is set as unity. A downward force F = 1 × 10 −3 is evenly distributed over a length of l = b∕15 at the top. The bottom of the design domain is fixed.
The stiffness optimization is conducted with b = 1 , V * = 25% , ER = 1% and r min = 4 × the element side length. The optimized design is shown in Fig. 1b, which appears to be a simple column and has the smallest BLF 1 = 0.89 , meaning the stiffness design would buckle under the applied load.
Therefore, the buckling resistance is improved by maximizing the lowest BLF 1 . In this case, the sensitivity of the BLFs [Eq. (9)] is the only term in the overall sensitivity, and the convergence criterion [Eq. (23)] applies to the KS aggregation of the BLFs. Figure 1c shows the evolutionary histories of the first four BLFs and the corresponding modal strain energy ratios r v j (j = 1-4) proposed by Gao and Ma (2015) for identifying pseudo buckling modes. r v j is defined as the contribution of the modal strain energy (in the jth buckling mode) of the element nodes in the void region to the total modal strain energy. It can be seen that the modal strain energy ratios are very close to zero throughout the optimization process, meaning the nodes in the void region contribute little to the total modal strain energy, and no pseudo-buckling modes occur in the first four buckling modes (Gao and Ma 2015). 1 (the blue curve) increases from 0.89 to 9.21 with compliance increased by 44.4% from 8.52 × 10 −6 to 1.23 × 10 −5 . The structural topologies at different iterations are shown in Fig. 1d-i. The initial column is separated into two bars connected by thin cross-like bars, and the distance between the separated bars is gradually increased in later iterations. It can be noticed that a sudden drop of 1 occurs at iteration 200 (Fig. 1c), which is due to the elimination of thin bars in the upper structure (see the bars highlighted by red circles in Fig. 1f). Afterwards, 1 gradually recovers as the remaining upper bars become thicker (see the bars highlighted by the red circle in Fig. 1g). The BLF maximization design of the compressed column obtained by Ferrari et al. (2021) using the SIMP method is shown in Fig. 1. As can be seen, the geometry of the present structure ( Fig. 1) is similar to the design of the SIMP method (Fig. 1).

Two-bar frame
The second example shows the geometrical characteristics of the compression and tension members in optimal designs obtained with different buckling constraints. Figure 2a shows the boundary condition and the design domain, which is discretized by 90 × 210. Plane stress elements are employed, and the thickness is set as unity. A downward force F = 2 × 10 −2 is spread over a length of l = b∕10 at the middle of the right side. The left side of the domain is fixed near the upper and lower ends with a length of l . b = 1 , V * = 20% , ER = 1% and r min = 2 × the element side length are employed in this example. The stiffness design is shown in Fig. 2b, consisting of upper and lower bars. The first three buckling modes of the twobar frame are illustrated in Fig. 2c, coloured from blue to red, which corresponds to the logarithm of the normalized strain energy density from low to high. It can be seen that the lower bar has a higher level of strain energy as a compressed member. On the contrary, the upper bar is under tension and has lower strain energy. The lowest BLF 1 = 0.89 , which means that the buckling occurs on the lower bar at the critical load F cr = 1 × F = 1.78 × 10 −2 .
The buckling-constrained optimization is then carried out based on the stiffness design (Fig. 2b). ER is reduced to 0.4% for stabilization consideration. The optimized designs with increasing buckling constraints are shown in Fig. 3. It can be seen that the material is relocated from the upper bar to the lower bar, which is divided into two bars connected by thin members to resist buckling. Besides, a bar connecting the upper bar and the lower bar appears, and it becomes thicker and longer with the increase of (see the bars highlighted by red circles in Fig. 3b-f). This bar increases the moment of inertia of the global structure and prevents the structure from rotating at the junction of the upper and lower bars (near the force).  Ferrari et al. (2021) With the increment of BLFs, the compliance of the buckling-reinforced design is higher than that of the stiffness design. The increment of the compliance compared to the stiffness design is ΔC = 2.6% when = 1.50, while ΔC can be 43.1% when = 4.00. The stiffness reduction is mainly due to the thinner upper bar as well as the separated and bent lower bar in the buckling-reinforced design. Figure 4 shows the evolutionary history of the compliance and the first four BLFs when = 4.00 (Fig. 3f). Generally speaking, 1 smoothly increases and the other BLFs gradually get close to 1 . It is noticed that the curves of 1 and 2 almost overlap at later stage, which is referred to as a multimodal case. As mentioned in Sect. 1, the sensitivity of the multiple BLFs is not unique, so that the nondifferentiable sensitivity may be obtained due to repeat eigenvalues and mode switching. Thus, the KS aggregation function is employed to approximate the maximum buckling constraint.
The first twelve buckling modes of the buckling design with = 4.00 (Fig. 3f) are shown in Fig. 5. It can be seen that these buckling modes are all real buckling modes, and no pseudo buckling modes are observed.   Table 1 shows the effect of the predefined update coefficient on compliance and required optimization iterations. When = 2.0 , by defining = 1.03 , 29.7% fewer iterations with only 0.1% extra compliance increment can be achieved compared to the case with = 1.01 . When = 3.0 , = 1.04 results in 0.7% extra compliance increment but 11.8% fewer iterations than = 1.02 . Larger leads to a larger k when it is updated, and it means a larger update step of the multiplier Λ . As shown in Eq. (12), the overall sensitivity is equal to the compliance sensitivity plus the buckling sensitivity multiplied by Λ . Thus, a larger update step of the multiplier Λ leads to a larger contribution of the buckling sensitivity to the overall sensitivity. In the design variable update process, the criterion of adding/deleting elements is determined by the sensitivity value (as shown in Sect. 2.3), which is more dependent on the buckling sensitivity when the multiplier Λ is larger.

Hollow square tube
This case shows a hollow square tube under external pressure, as illustrated in Fig. 6a. The width of the cross-section b = 1.9 , and the thickness of the tube t = b∕10 . The tube length is much larger than the dimension of the cross-section, hence plane strain elements are employed. Each side of the tube is under a pressure load P = 1.05 × 10 −3 per unit area. The cross-section is optimized using the proposed algorithm. As the cross-section of the tube is symmetrical, to save computational costs, a quarter of the cross-section (shaded blue area in Fig. 6a) is used, and the corresponding boundary conditions are shown in Fig. 6b. The width of the quarter section is b∕2 , which is divided into 250 elements.  First twelve buckling modes of buckling design with − = 4.00 . Buckling modes are coloured from blue to red, corresponding to the logarithm of the normalized strain energy density from low to high As a quarter model, the vertical degrees of freedom of the left edge and the horizontal degrees of freedom of the bottom edge are fixed. The cross-section is designed as a sandwich structure, i.e., the inner and outer layers are non-design domains, each of which has a thickness of t∕10 (discretized using 5 elements), and the core area is the design domain.
Set V * = 50% , ER = 1% , and r min = 6 × the element side length for stiffness optimization. The stiffness design and the recovered full model are shown in the first row of Fig. 7a and b. It can be seen that the optimum designs are oblique symmetric due to the even-distributed pressure. The principal stress distribution (first row of Fig. 7c and d) also shows the symmetric properties of the optimum designs. The lowest BLF 1 of the stiffness design is 0.54, corresponding to the critical load P cr = 1 × P = 5.67 × 10 −4 per unit area.
Buckling-constrained optimization is then conducted based on the stiffness design with ER = 0.2% , r min = 3 × the element side length, and = 1.00 . The optimization results are shown in the second row of Fig. 7a and b. As can be seen, the material is relocated from the corners of the tube to the middle of each side (second row of Fig. 7b) and from tension members to the compression members ( Fig. 7c and d). The lowest BLF 1 and the compliance increase from 0.54 to 1.00 and 7.47 × 10 −5 to 8.12 × 10 −5 , respectively (Fig. 7e). Figure 8 presents the buckling analysis of stiffness and buckling designs. It can be seen that the middle part at each side of the stiffness design (first row of Fig. 8) has a higher level of strain energy density. However, in the buckling-reinforced design (second row of Fig. 8), as the material is relocated to the middle of each side, the strain energy density distribution is more uniform, and the buckling phenomenon can be well suppressed.
Furthermore, in the fundamental buckling mode of the 3D model, only in-plane buckling occurs, and no buckling is observed in the longitudinal direction. Therefore, an optimized design for preventing buckling in the crosssection is also of great significance. Figure 9 shows the buckling-constrained optimization when starting from the stiffness design obtained with r min = 3 × the element side length (other parameters remain the same). The stiffness design (first row of Fig. 9) possesses a better mechanical performance compared to that of Fig. 7, i.e., a higher fundamental BLF 1 of 0.85 and a lower compliance of 7.07 × 10 −5 . The buckling-constrained optimization design (second row of Fig. 9) has a compliance of 7.27 × 10 −5 and the geometry is slightly different from that obtained when starting from the stiffness design with r min = 6 × the element side length (the second row of Fig. 7a and b). The most distinct differences are the extra bars in the middle of each side, highlighted by red circles in Fig. 9b.
Oscillations can be observed in the history of 1 , indicated by red circles in Fig. 9c. While no oscillation occurs when starting from the stiffness design obtained with r min = 6 × the element side (Fig. 7e). The use of the stiffness design obtained with a larger r min as the initial design domain for buckling-constrained optimization can enhance structural redundancy. A higher level of structural redundancy may help prevent the BLF oscillations caused by the elimination of structural components during the optimization process. Figure 10 shows the buckling-constrained optimization design obtained by starting from a full design domain with V * = 50% , ER = 1% , r min = 3 × the element side length, and = 1.00 . The final design has a compliance of 7.43 × 10 −5 and the geometry ( Fig. 10a and b) is similar to that shown in Fig. 9. The evolutionary histories of compliance and the first two BLFs are shown in Fig. 10c. The volume fraction is reached at the 69th iteration and therefore the compliance Fig. 7 Results of stiffness design (first row) and buckling-constrained optimization (second row): a a quarter of cross-section; b recovered full model; c maximum principal stress ( 1 ); d minimum principal stress ( 2 ); e evolutionary histories of compliance and the first four BLFs changes little in the subsequent optimization process. The buckling constraints are activated at the 42nd iteration when 1 < and satisfied at the last (116th) iteration. Slight oscillations can be observed in the evolutionary history of 2 in Figs. 7e, 9c and 10c (and 3 in 7e). The amplitude between two consecutive iterations is within 0.2, and such small amplitudes have no harm to the KS aggregation of BLFs and the stability of the optimization process. Besides, it can be seen that the result of the buckling-constrained optimization is dependent on the structural property of the initial design. In this case, the buckling-constrained design has the lowest compliance when starting from the stiffness design obtained with r min = 3 × the element side length.

L-Shape beam
This case shows the buckling-constrained topology optimization of a long L-shape beam subjected to a linear load on the edge. As shown in Fig. 11a, the height of the section shape h = 1 , and the thickness t = 0.4 × h . The length of column is much larger than the height, and the plane strain elements are used. The top surface of the column is fixed, and a linear load q = 3 × 10 −3 per unit length is applied to the edge. As shown in Fig. 11b, the section shape is used as the design domain, and the height and the width are discretized using 300 elements, respectively. The force is spread on a length of h∕20 on the right edge to avoid stress concentration near the load point.
The first row of Fig. 12a shows the stiffness design, which is obtained by using V * = 40% , ER = 1% , and r min = 10 × the element side length. The maximum and the minimum principal stress distributions (Fig. 12b and c) show the tension and compression members of the structure. It can be seen that the members on the left and bottom are in compression, while those in the middle and on the top right side are in tension. The lowest BLF of the stiffness design is only 0.30, and the critical load q cr = 1 × q = 9 × 10 −4 per unit length.
To obtain a design with higher buckling resistance, the buckling-constrained optimization is carried out using ER = 0.2% , r min = 3 × the element side length, and = 1.00 . The obtained results are shown in the second row of Fig. 12. It can be seen that the material is transferred from the tension members to the compression members, and the moments of inertia of compression members are larger, leading to better buckling resistance. Figure 13 shows the first four buckling modes of the stiffness design (first row) and the buckling-reinforced design (second row). The buckling phenomenon is observed in the compression members from both designs. The stiffness design exhibits a more localized buckling effect, whereas the buckling design exhibits a more global buckling effect. For instance, buckling occurs on multiple bars simultaneously in the first buckling mode ( 1 ), and the strain energy density distributes more evenly in buckling design. Compared to the stiffness design, the compliance of the buckling design increases by 4.2%.

Fig. 8
First four buckling modes of stiffness design (first row) and buckling-constrained optimization (second row). Buckling modes are coloured from blue to red, corresponding to the logarithm of the normalized strain energy density from low to high The numerical examples show that stiffness designs may exhibit poor stability in practical applications. In comparison, the proposed topology optimization with the consideration of buckling constraints can greatly improve the buckling resistance of a structure with only a small increase in compliance.

Conclusions
In this study, a novel buckling-constrained topology optimization algorithm is developed based on the BESO method. An augmented Lagrangian multiplier is proposed to integrate buckling constraints into the objective function. A design variable update scheme is developed to control the evolutionary rate when the volume fraction is reached. Four topology optimization designs are presented to show the validity and effectiveness of the proposed algorithm. Based on the numerical results, the main conclusions are summarized as follows.
(1) By using the proposed augmented Lagrangian multiplier, the contribution of the buckling sensitivity to the overall sensitivity can be controlled by the augmentation , which is adjusted based on a user-defined coefficient . A larger can save the required number of iterations but slightly increase the compliance.
(2) With the implementation of the augmented Lagrangian multiplier, the modified design variable update scheme and the enhanced stabilization scheme, excellent stability and convergence can be obtained from the proposed optimization process.

Fig. 9
Results of the stiffness design obtained with r min = 3 × the element side length (first row) and buckling-constrained optimization (second row): a a quarter of cross-section; b recovered full model; and c evolutionary histories of compliance and first two BLFs   (3) The lowest BLF 1 of the hollow square tube is improved from 0.54 to 1.00 with a small increment of compliance (8.7%). 1 of the L-shape beam is increased from 0.30 to 1.00 while the compliance is increased by only 4.2%.
This study contributes to the knowledge of bucklingconstrained topology optimization by providing a novel algorithm and a series of numerical examples. As most current studies associated with buckling-constrained optimization focus on 2D examples, further research should extend the algorithm into 3D to solve more practical engineering problems. Experimental investigations on the buckling and post-buckling behaviour of structural components are needed to confirm the effectiveness of the algorithm in real engineering applications.
Funding Open Access funding enabled and organized by CAUL and its Member Institutions. The authors gratefully acknowledge the financial support from the Australian Research Council (FL190100014, DP200102190).

Conflict of interest
The authors declare that they have no conflict of interest.
Ethical approval This article does not contain any studies with human participants or animals performed by any of the authors.

Replication of results
The results of the optimized designs and the basic code of this work are available from the corresponding author upon reasonable request.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.