Minimal compliance topologies for maximal buckling load of columns

This paper presents a new approach to the topology optimization of columns exposed to a loss of stability. The idea is to replace a conventional maximization of a buckling load by a locally formulated topology optimization problem based on compliance minimization. In order to do this, the standard instability analysis of a compressed column is performed first and the buckling mode is determined. Then the compressive loading is replaced by a transverse one which is selected so as to generate a bending moment, the distribution of which coincides with the one representing the considered buckling mode. Minimization of compliance is performed for the bent structure and optimal topology is generated. Finally, the critical load for the optimal column is calculated. The selected numerical results obtained with the use of the above technique, under the assumption that the buckling load is unimodal, are presented. The above approach allows generating optimal topologies using Cellular Automata (CA) method, which by nature requires local formulation of the design problem.

paper by Bendsoe and Kikuchi (1988) one can find in the literature numerous approaches to generating optimal topologies based both on optimality criteria and evolutionary methods. A general overview as well as a broad discussion on topology optimization concepts are provided by many survey papers e.g. Rozvany (2001), Eschenauer and Olhoff (2001), Bendsoe and Sigmund (2003), Arora and Wang (2005) or Rozvany (2008). At the same time hundreds of papers present numerous solutions including classic Michell examples as well as complicated spatial engineering structures, implementing specific methods ranging from gradient based approaches to biologically inspired algorithms and level set method, e.g. Wang et al. (2003).
Topology optimization of structures exposed to a loss of stability or vibrating rarely appears in the literature, hence the number of publications dealing with this subject is rather modest. One can point out early papers by Neves et al. (1995) where an approach to introduce buckling load control into topology optimization model was presented, Ma et al. (1994) in which topology optimization was applied for obtaining desired eigenfrequencies, or Min and Kikuchi (1997) devoted to optimal reinforcement design of structures under buckling load. Later on Cheng et al. (2000) reported design of elastic structures for stability, Zhou (2004) presented optimum designs of shell topological optimization under linear buckling response, while Bian and Sui (2009) studied topology optimization of plane elastic structures with buckling and displacement constraints. Recently Lindgaard and Dahl (2013) discussed topology optimization of static geometrically nonlinear structures experiencing snap-through behaviour. As selected from among papers discussing generation of optimal topologies for vibrating elements Pedersen (2000) and Gong et al. (2012) optimized free vibration eigenvalues whereas Du and Olhoff (2007) also optimized frequency gaps. It is worth noting that buckling constraints more often appear in topology optimization of truss structures. The papers by Birker (1996) and Bojczuk and Mróz (1999) and Guo et al. (2005) may serve here only as exemplary ones.
This paper presents a new approach to topology optimization of columns exposed to a loss of stability. The idea is to replace the conventional maximization of buckling load by a locally formulated topology optimization problem based on compliance minimization.
The Cellular Automata method is used in this paper as optimal topology generator. The first application of CA to optimal structural design, and to topology optimization in particular, was proposed by Inou et al. (1994) and Inou et al. (1997). In these papers the design domain was divided into cells, the states of which were represented by the Young moduli of the material as design variables. By applying iteratively the local CA rules the values of the elastic moduli for all cells were updated based on the difference between current and target stress values. The cells with low values of elastic moduli were removed. The idea of implementation of CA to optimal design was described also by Kita and Toyoda (2000), Haleja and Kim (2001), and by Tatting and Gurdal (2000), where the authors proposed a new scheme for CA, in which analysis and design were performed simultaneously -simultaneous analysis and design. This technique has been modified and extended, for example by Cortes et al. (2005), Setoodeh et al. (2006), Canyurt and Haleja (2007). During the last two decades implementation of CA in structural design has been under permanent development, and numerous papers related to application of CA to topology optimization, see e.g. Missoum et al. (2005), Abdalla et al. (2006), Hassani and Tavakkoli (2007), Penninger et al. (2009), Sanaei and Babaei (2011), Bochenek and Tajs-Zielińska (2012 or Du et al. (2013), have been published. In addition, the series of papers by Tovar and co-workers should be mentioned: Tovar et al. (2004a, b), Tovar et al. (2006, Penninger et al. (2011), in which a new CA technique inspired by a process of functional adaptation taking place in bones has been implemented.
It is worth noting that CA approach to structural optimization requires a local formulation of the design problem by nature, and therefore applying it to the maximization of a structure buckling load, which is a global quantity, is not straightforward. Fortunately one can observe that for the optimal column, for which critical load has been maximized, the maximal bending stress is uniformly distributed along column axis. Taking that into account it is possible to replace conventional maximization of buckling load by a problem formulated as the fully stressed design. This concept was presented by Bochenek and Tajs-Zielińska (2012), where the locally formulated problem of optimal sizing of columns prone to instability was considered. In fact while conducting research on CA application to structural design under stability constraints the mentioned above new idea has come out. The detailed description of the proposed concept, namely how to generate topologies for the maximal buckling load of columns by adapting for this purpose a typical compliance minimization problem, is given in the following section.

Basic idea
As the basic problem, the standard maximization of the buckling load of a column subject to a given total volume, is considered. If a cross-section area A(x) of the column is chosen as the design function, and its relation to moment of inertia complies with I (x) = cA 2 (x), then the optimal solution for the problem, takes the following form e.g. Seyranian and Privalova (2003): (1) It has been assumed that the optimal solution is represented by only a single buckling mode, which is described by the bending moment function M g (x). The quantity c 1 is a constant, the value of which can be selected so as to fulfill the volume constraint. If as an alternative the optimization problem is posed as minimization of compliance of a column/beam of length L under transverse loading with respect to the constant volume constraint: where λ stands for the Lagrange multiplier, then the optimality condition applied to (2) results in: One can observe that the above relation is analogous to (1). This observation leads to the following conclusion. It is possible to replace a conventional maximization of buckling load by the optimization problem based on compliance minimization. The implementation of this idea is a main goal of the present paper. The concept of the implementation is as follows. First the standard instability analysis of a compressed column is performed and the buckling mode is determined. Then the compressive loading is replaced by a transverse one which is selected so as to generate the bending moment, the distribution of which coincides with the one representing particular buckling mode. In order to do that a number of forces distributed along column axis is specified and then their values are sought for with a view to get expected distribution of bending moment in several nodal points. For the bent structure minimization of compliance is performed, and optimal topology is generated. Finally the critical load for the optimal column topology is calculated.

Illustrative example
In order to illustrate this methodology the following example which makes use of the above idea is presented. A simply supported column of rectangular cross-section, presented in Fig. 1 has been chosen. The transverse loading is applied in one plane which stands also for the buckling one (Fig. 2). The well known solution of the Euler buckling problem leads to instability mode described by a sine function. This refers to both the deflection and the bending moment.
The optimal topology has been generated for reference column which consists of 14 × 14 × 200 cells (1 cm × 1 cm × 1 cm). The Young modulus of the column material is E = 200 GPa and the volume fraction of the final topology is set to 0.5. The minimal compliance column is given in Fig. 3. For the obtained topology the critical load 12592 kN has been calculated numerically using finite element system. The buckling mode is presented in Fig. 4. It is worth noting that the critical load calculated for prismatic column of the same volume is 3790 kN.
While considering the possibility of buckling in both planes the above solution found for only one plane taken into account may not be correct if the critical load referring to out-of-plane buckling has a lower value. Indeed for the considered example out-of-plane buckling load is 5521 kN.  In order to eliminate this drawback the formulation of topology optimization problem is extended so as to take into consideration the loadings applied in both planes, see Fig. 5. As the result the topology shown in Fig. 6 has been generated, for which critical load value is 11701 kN. To make details of the final topology better recognized, a section view is added in Fig. 7, together with a magnified view of a part of generated topology shown in Fig. 8.
It is worth noting that assumption of doubly symmetrical column cross-section implies that instability considerations are limited to buckling in principal planes.
The above presented introductory example shows that it is possible to achieve the effect of buckling load maximization via minimization of compliance performed for a specially matched loading scheme. It has been assumed for the above that the buckling load is unimodal.
At the same time it is worth pointing out that for clamped-clamped column the unimodal formulation of the problem may be not sufficient to obtain optimal solution. If such situation occurs the modified -bimodal formulation of the problem has to be implemented.
In order to generate optimal topologies the efficient algorithm is needed. The one used in this paper is based on Cellular Automata rules. The rules of Cellular Automata are briefly described in the next section.

A Cellular Automata based topology generator
As it has been mentioned earlier in this paper the problem of topology optimization is formulated in a manner allowing to solve it with Cellular Automata method. The performance of CA algorithms, reported in literature, is based mostly on heuristic local rules. Similarly, in the present paper the efficient heuristic algorithm, introduced by Bochenek and Tajs-Zielińska (2012) has been implemented.
The power law approach defining solid isotropic material with penalization (SIMP) with design variables being relative densities of a material has been utilized. The elastic modulus of each cell element is modelled as a function of Fig. 4 The in-plane buckling mode The local update rule applied to design variables d i associated with central cells is now constructed based on the information gathered from adjacent cells forming the Moore or von Neumann type neighborhood. It is set up as linear combination of design variables corrections with coefficients, the values of which are influenced by the states of the neighborhood surrounding each cell, as presented in (5): The compliance values calculated for central cell U i and N neighboring cells U ik are compared to a selected threshold value U * , e.g. a global average for whole structure in every iteration. Depending on whether U i and U ik are larger or smaller than the selected U * positive or negative coefficients C α 0 for central cell and C α for surrounding cells are chosen and transferred to the design variable update, as given by (6) and (7):  Practical implementation of proposed local rules requires specification of introduced parameters. The values of C α 0 and C α are selected so as to provide that their sum for all cells within the neighborhood equals one. Based on numerical tests performed for von Neumann type neighborhood all coefficients equal to 0.2 for plane problems seems to be a good choice. The move limit m (m = 0.2 is proposed) implemented in the above algorithm controls the allowable changes of the design variables values.
The numerical algorithm has been build in order to implement the above proposed design rule. As to the optimization procedure the sequential approach, has been adapted, meaning that for each iteration, the structural analysis performed for the optimized element is followed by the local updating process. Simultaneously a global volume constraint can be applied for specified volume fraction. If so the generated optimal topology preserves a specified volume fraction of a solid material.
Below, pseudo code of the algorithm build for generation of compliance based topologies is presented: The above rules have been already successfully applied to the generation of plane and spatial topologies of elastic structures, and numerous results are presented and discussed in Tajs-Zielińska (2012, 2013).

Generation of column topologies for maximal buckling load
Following the preliminary example of Section 2 selected column topologies for various boundary conditions are generated.

A simply supported column
The simply supported column of Section 2 is taken into account once again. This time calculations are performed for the value of the volume fraction set to 0.35. The result is given in Fig. 9. For the obtained topology the critical load 10083 kN has been calculated. Concurrently the out-of-plane buckling load of 2328 kN has been found.
Considering two buckling planes leads to topology presented in Fig. 10, for which the critical load value is 8760 kN. To make details of the final topology more visible, a section view is given in Fig. 11. It has to be pointed out that the critical load calculated for prismatic column of the same volume is 1910 kN, which shows a significant buckling load increase.
The above calculations are now repeated for a slenderer column 14 × 14 × 400 cells (1 cm × 1 cm × 1 cm). The resulting topology for volume fraction 0.5 is given in Fig. 12. For the obtained topology the critical load 3367 kN has been calculated. For the considered example out-of-plane buckling load is 2109 kN.
Taking into account two buckling planes proceeds to topology presented as a section view in Fig. 13, for which the critical load value is 3122 kN. Supplementary to that a magnified view of a part of generated topology is shown in Fig. 9 One buckling plane. Critical load 10083 kN

A clamped-simply supported column
The next example regards a clamped-simply supported column, see Fig. 15.
The optimal topology has been generated for a reference column which consists of 20 × 20 × 200 cells (1 cm × 1 cm × 1 cm). The volume fraction of the final topology is set to 0.5. The result is given in Fig. 16, for which the critical load 89595 kN has been calculated. For out-of-plane buckling critical load is 53687 kN.
For two buckling planes considered, generation of optimal topology conducts to the result presented in Fig. 17. The critical load is 85327 kN. A section view is shown in Fig. 18, whereas a magnified view of a part of generated topology is presented in Fig. 19.
The critical load calculated for a prismatic column of the same volume is 33650 kN.

A clamped-clamped column
The optimal topology has been generated for a reference column which consists of 20 × 20 × 200 cells (1 cm × 1 cm × 1 cm). The volume fraction of the final topology is set to 0.5. The result of topology generation for in-plane buckling is given in Fig. 20. For the final topology the critical load is 149814 kN. Simultaneously critical load for out-of-plane buckling is 98330 kN.
Considering two buckling planes results in the topology presented in Fig. 21. The critical load is 148230 kN.
To show the details of the final topology, a section view is given in Fig. 22, together with a magnified view of a part of a generated topology shown in Fig. 23.
The critical load calculated for a prismatic column of the same volume is 65691 kN.     The above topologies have been obtained for the fundamental, symmetric buckling mode of the column (see Fig. 24). It is known, however, that for this type of clamping the buckling mode which is fundamental for a prismatic column not necessarily remains the first one for optimized structure. This problem was pointed out already by Olhoff and Rasmussen (1977) for optimal sizing of clampedclamped column. Having this in mind critical loads for the first two buckling modes in each of the buckling planes of the column in Fig. 21 have been evaluated. It occurred indeed that critical load 127528 kN obtained for an antisymmetric second buckling mode (see Fig. 25) is smaller than the load 148230 kN calculated for the symmetric buckling mode shown in Fig. 24. This means that the column in Fig. 21 is not optimal for the problem under study, and indicates the necessity of bimodal formulation of an With a view to resolve the above case, the following problem can be formulated. Let us search for a volume fraction for which both critical loads would have equal values. Solving the above formulated problem leads to the volume fraction 0.67, for which bimodal critical load equals 190513 kN. The topology of this bimodal column is presented in Fig. 26.

Concluding remarks
The new approach to topology optimization of columns against instability has been introduced, and some results of its application to generation of column topologies for various boundary conditions have been presented.
The critical loads for generated topologies are significantly higher as compared to prismatic columns of the same volume. It should be also stressed that critical loads obtained for columns of optimal topology are also higher than the ones found for columns for which the standard approach of optimal sizing against buckling has been applied.
While building a computational model of the column prone to buckling, within the approach proposed in this paper, compressive loading is replaced by a distributed transverse one. Distribution of this transverse loading along column axis, requires a number of forces to be chosen. It is worth noting that for different numbers of forces the final topologies may slightly differ from each other. On the other hand, values of critical loads, evaluated for various intensities of transverse loading, practically do not vary. To confirm that a comparison has been made for clampedsimply supported column discussed in Section 4.2. The numbers of transverse forces are selected in sequence as 3, 5, 9, 19 and 39, for which values of buckling loads 85465 kN, 86500 kN, 85672 kN, 85327 kN and 84915 kN are obtained, respectively.
It was shown in Section 4.3 that for a clamped-clamped column unimodal formulation of optimization problem may be not sufficient. In such a case the bimodal approach which allows generating topologies for which both critical loads have the same value is necessary. The exemplary result of such an approach has been presented in the paper.