//Logo Image
Author: Ming-Hsiu Hsu, Yeh-Liang Hsu(2004-03-19); recommended: Yeh-Liang Hsu (2004-10-27).
Note: This paper is published in Engineering Optimization, Vol. 37, No. 1, January 2005, p. 83–102.

Generalization of two- and three-dimensional structural topology optimization

Abstract

This paper explores the fundamental issues on the quality of structural topology optimization results, and presents a generalized topology optimization process. In the generalized process, non-rectangular design domains with geometrical constraints can be accepted, and the use of an automatic mesh generator to mesh the design domain is allowed. The proper number of elements in the design domains to avoid the mesh-dependence problem is also suggested. Higher order elements are used to deal with the checkerboard problem, and a two-stage penalty function method is proposed for the topology optimization. Finally a continuity analysis is used to deal with the porous topology, and two filters are implemented to filter out the trivial solids and voids. The whole process is generalized to two- and three- dimensional structures and can be fully automated. Fourteen two-dimensional and six three-dimensional examples are used to demonstrate the effectiveness of this process.

Keywords: topology optimization, structure optimization

1.     Introduction

Topology optimization is a well-developed field. The purpose is to obtain the optimal layout of structural components to achieve a predetermined performance goal. The objective function and constraints commonly used in topology optimization research are to minimize the compliance of the structure under a given loading condition, subject to the amount of usable material. The same objective function and constraints are also used here.

Bendsfe and Kikuchi developed the homogenization method for topology optimization in 1988 [1], which is a milestone for topology optimization. The homogenization method was based on the assumption of a microstructure in which the properties are homogenized. Two other approaches were proposed for topology optimization in the early 90’s: the density distribution method [2, 3], in which the material density of each element was selected as the design variable, and the evolutionary structural optimization technique (ESO) [4], which gradually removed the low-stressed elements to achieve the optimal design. Rozvany provided a complete survey on topology optimization [5]. All three methods have their own advantages and have many followers. The density distribution method, the so-called Simple Isotropic Material with Penalization (SIMP) method, seems to be the most convenient method and is used here to obtain topology optimization result.

Figure 1 shows a general topology optimization process. The process begins with a user-defined design domain . Figure 2 shows a cantilever beam example that has a rectangular design domain. The design domain may also have geometrical constraints, namely, void regions and fixed regions as shown in the suspension arm example in Figure 3. A general design processor, for example, Pro/Engineering or AutoCAD, can be used to construct the design domain and to create a corresponding CAD model for topology optimization.

Figure 1. A general topology optimization process

Figure 2. The design domain and boundary conditions of the cantilever beam example

Figure 3. The suspension arm example with hole and support regions in design domain

An analysis processor, for example, ANSYS or NASTRAN, generates the finite element model from the CAD model. The boundary conditions, for example, displacement constraints on boundaries , loads and body force , and surface force , are applied to the finite element model.

The user defines the amount of usable material , and selects true material properties including Young’s modulus , material density , and Poisson’s ratio . The design variables in the material density distribution method are the “normalized material densities” ri,j (the normalized material density of the i-th element at the j-th iteration). Note that 0ri,j1, and ri,j = 0 denotes that the area of this element contains no material, while ri,j = 1 denotes that the area of this element contains material. The initial value of the normalized material density ri,0 is the ratio of the amount of usable material M0 to the total material if the whole design domain is occupied.

The density that is actually used in the finite element model can be calculated by

        = ri,j r0                                                                                             (1)

The Young’s modulus of the i-th element at the j-th iteration actually used in the finite element model is assumed to be:

                                                                                                    (2)

where  is an exponent usually between 2 and 4 to penalize the small  [3, 6, 7]. In this research, the value of  is taken as 2.

The finite element analysis solver, for example, ANSYS or NASTRAN, is used to analyze the finite element model under given boundary conditions and loads. The displacement of nodes and the strain energy of elements  are evaluated by the solver. At the j-th iteration, the compliance of the structure  and its sensitivities with respect to the normalized density of each element can then be evaluated by the following equations [7]:

                                                                               (3)

                                                                                           (4)

This objective function value and its sensitivities are required when solving the optimization model.

User-defined termination criterion is used to check if the topology optimization process should be stopped. If the result is false, a topology optimization model that minimizes compliance of the structure and subjects it to the usable material constraint will be generated:

                                                              (5)

An optimization solver is used to find a new set of normalized material densities for the next iteration. The normalized material density will then be used to compute material properties  and  used in the finite element model in the next iteration using Equations (3) and (4).

If the termination criterion result is true, the topology optimization process will stop and output the set of final normalized material densities of each element,  in Figure 1. This topology optimization result can be plotted using different gray levels in the elements to represent the corresponding , as shown in Figure 4.

Figure 4. A cantilever beam example with poor topology optimization result

After the important advance by Bendsfe and Kikuchi [1] in the field of topology optimization, many researchers paid attention to integrating structural topology optimization and shape optimization since early 1990s. Some researchers treat structural design optimization as a three-phase design process, and aim at developing an automated design process [8, 9, 10]. Phase I is the topology generation process, in which the optimal topology of the structure is generated. Phase II is the topology interpretation process. Various approaches are used to interpret the topology optimization result. Phase III is the detailed design phase, in which the shape and size optimization are implemented.

The approaches used in the Phase II interpretation process can be roughly divided into three categories: image interpretation approach [8-13], density contour approach [14-16], and geometric reconstruction approach [17]. The image interpretation approach uses graphic facilities or computer vision technologies to represent the boundary of the black-and-white finite element topology optimization result. The density contour approach generates the boundaries of the structure by redistributing densities from the topology optimization result. In the geometric reconstruction approach, the boundaries are represented by the mathematical geometric reconstruction technique.

The quality of the topology optimization result generated in Phase I directly affects the implementation of the interpretation process in Phase II. For example, it will be very difficult to interpret the poor topology optimization result shown in Figure 4 using any interpretation technique.

Thus, several fundamental issues with respect to the quality of topology optimization result must be considered in the topology optimization phase:

(1)   Could the design domains be shapes other than rectangular ones?

(2)   Could the automatic mesh generation be used to create the finite element model?

(3)   What is the proper element size?

(4)   How can the checkerboard problem be dealt with?

(5)   How can a strictly black-and-white topology optimization result be obtained?

(6)   How can the result of topology optimization be interpreted to ensure that a continuous structure is obtained?

In the following sections, these issues are explored and suggestions are made for the generalization of the topology optimization process to general two- or three-dimensional structural design problems.

2.         The mesh generation and mesh-dependence problem

The structural design examples in the topology optimization literature often have rectangular design domains and are manually meshed using square elements, such as the cantilever beam example in Figure 4. Realistic structural design examples may have complicated design domains with geometry constraints, such as the suspension arm example shown in Figure 3. It is not easy to mesh the complicated design domain manually using square elements. Bendsfe [18] suggested that the use of an automatic mesh generator would, of course, simplify the treatment of problems with complicated design domains.

As the CAD model is developed or generated into finite element model, the design domain  is divided into  nodes and  finite elements. The topology optimization model in Equation (5) can be rewritten in a discretized form as:

                                                                    (6)

where  and are the body and surface forces applied to node k;  is the volume or area of element i.

The quality of topology optimization result is dependent on the discretization of the finite element model, the so-called “mesh-dependence problem.” Figure 5(a) shows a simply supported beam example to illustrate the mesh-dependence problem using the density distribution approach [19]. Figure 5(b) shows the topology optimization result for the discretization using 600 elements, and Figure 5(c) shows the result using 5,400 elements. The result in Figure 5(c) is much more detailed than that in 5(b), and more importantly, the two topologies are different in nature.

Figure 5. The simply supported beam example for mesh-dependence problem [19]

It has been shown that the mesh-dependence problem often relates to the problem of nonexistence of solutions. Several restriction methods were proposed in the literature to avoid this type of problems. In the “perimeter control method,” an upper-bound constraint on the perimeter is used to ensure a well-posed design problem [20, 21]. In “mesh independent filtering,” the filter modifies the design sensitivity of the elements within a specific zone to guarantee the radius of the members [22, 23]. In “global gradient constraint” method, the density function and its gradient variation are bounded to ensure the existence of a solution [18]. In “local gradient constraint” method, a local gradient constraint on the density variation is introduced to guarantee the existence of solution [24, 25]. Further details of these methods can be found in the survey by Sigmund and Petersson [19].

The model with a finer finite element mesh not only results in a better description of boundaries but also increases the complexity of the resulting topology. Zhou et al. [26] consider the manufacturability of the design from the engineering viewpoint, instead of worrying about the fact that different mesh densities may result in different final solutions. Therefore, in their work, the algorithm of minimum member size control is developed to improve the manufacturability of the design. Haber et al. [21] reported that the reduction in compliance obtained by increasing the complexity of the design topology is the modest. Effective designs can be obtained with relatively simpler topologies. Following this finding, a “convergence test” is used to find a proper topology optimization result. Table 1 shows the topology optimization results using different number of elements, and Figure 6 plots the number of elements vs. final compliance. It is observed that as the number of elements increases the resulting compliance decreases but converges after the number of elements exceeds 2,000. Table 1 also shows that, as the number of elements increases the topology of the structure becomes more complicated, which is impractical.

Table 1. The convergence test of the cantilever beam example

No. of Elements

Final compliance

Topology optimization results

No. of Elements

Final compliance

Topology optimization results

112

23286

2432

12702

135

20760

2680

12435

198

20289

3053

12524

240

18928

3375

12686

360

17174

3840

12690

540

16688

4472

12173

960

13973

5208

12434

2160

12461

6000

12604

MATLAB Handle Graphics

Figure 6. Final compliance vs. number of elements

The same convergence test is also applied in the example of automatic-meshed suspension arm. The convergence test results given in Figure 7 are similar to the result of cantilever beam example. The compliance tends to converge when the number of elements exceeds 2,000. Although this element number should be problem-dependent, for the 14 two-dimensional topology optimization examples and the 6 three- dimensional topology optimization examples presented later in Table 4, the element number around 2,000 has been a good reference number.

MATLAB Handle Graphics

Figure 7. Final compliance vs. number of elements

3.         The penalty function

Ideally, the normalized material densities of the elements should be either 0 or 1 after topology optimization. But there can still be intermediate normalized material densities at the end of the topology optimization, which results in an unreasonable structure. Using a penalty function to penalize the intermediate normalized material densities is the most common strategy to obtain a “black-and-white” topology.

Bensfe [18] embedded a penalization in the material property that is actually used in the finite element model by the following power law formulation:

                                                                                           (7)

where  is the Young’s modulus of the i-th element at the j-th iteration and  is the Young’s modulus of a given isotropic material. The penalization factor that is larger than 1 is represented by p. A similar penalization method was also used by Sigmund [23] and Zhou et al. [26]

Kumar and Gossard [14] added a penalty function to their objective function to penalize the intermediate densities:

                                             (8)

The last term of Equation (8) is the penalty function, in which cp is a penalty constant and is increased in steps. It can be verified easily that the penalty function has the maximum value when  and has the minimum value when  and .

Chen and Wu [7] presented a similar penalty function to minimize the number of elements with normalized densities not equal to either 0 or 1. The penalty function is included in the objective function, and is defined as

                                                                                                 (9)

where M is the total number of finite elements and R is a given penalty parameter which can vary iteration by iteration. This penalty function has several advantages. The function is continuous and differentiable and is symmetric about the most un-wanted density value of 0.5. Secondly, as a result of the nonlinear nature of this penalty function, heavier penalties are imposed on those densities that are close to 0.5. The sum of the penalties is an indicator of the clearness of the topology and is small if most densities are close to 0 or 1.

Borrvall and Petersson [27] treat the penalty function as an explicit constraint:

        g(r)=(- r)(r -)                                                                  (10)

where  for pure material-void problem, in which  is either  or  at almost each element. The value ep can be small upon allowing some intermediate density. In order to guarantee the existence of solutions, a compact linear operator  is used to regularize the constraint. Thus, the constraint  is replaced by , and is called “regularized intermediate density control”. The  is called “regularized penalty function,” and is defined by a composite mapping .

In most approaches discussed above, there are penalty parameters that control the amount of penalty, for example, p in Equation (7),  in Equation (8), and R in Equation (9). The value of the penalty parameter has been shown to regularly increase with iteration times. If the value of the penalty parameter is small, the effect of penalizing the intermediate densities is not apparent. On the other hand, if the value is large, it will affect the objective function value, which could lead to non-optimal topology results. The initial value and the amount of increment of the penalty parameter are often decided by experience of the numerical process and are dependent on individual cases.

In order to avoid numerical problems, such as premature convergence and the difficulty of the densities crossing 0.5 arising from the addition of the penalty function in the early stage, Chen and Wu [7] presented a two-timing approach to add the penalty function. The first timing is at the fifth iteration, when the topologies have been formed roughly. The second timing is at the iteration when the number of densities whose values are greater than 0.95 and less than 0.05 is over 50% of the total number of the density values.

In the research presented here, the topology optimization process is divided into two stages. The purpose of the first stage is to reduce the objective function value, and the purpose of the second stage is to drive the normalized material densities to either 0 or 1. The amount of penalty added in each stage is planned ahead according to the different purposes of each stage.

The user defines the total number of iterations  as one of the termination criteria. The user also defines the ratio of the number of iterations of the first stage to the total number of the iterations . In order to have an adequate time to reduce the objective function, the ratio  is suggested to be 75%. Therefore, the first stage starts at the initial iteration  to 75% of the total number of iterations (, rounded to an integer).

The penalty function at the j-th iteration is defined as follows:

                                                                                      (11)

where  is a penalty parameter and is increased iteratively. The function  is the base of penalty function defined as

                                                                           (12)

where M is the total number of finite elements. As shown in Equation (12),  is a reward function to reward the normalized material densities close to either 0 or 1. Equation (12) is also a continuous and differentiable function, and is symmetric at about 0.5. The function receives maximum reward when the normalized material density is 0 or 1, and receives no reward when it is 0.5.

In Equation (11), the penalty parameter  is defined by the amount of “reward” intended for each stage, instead of numerical experience. At the end of the first stage, the magnitude of  is set to be equal to initial compliance. At the end of the second stage, the magnitude of  is set to be equal to 10 times of initial compliance. Thus,  can be derived from the following equations:

         for the first stage,

         for the second stage                                              (13)

The multipliers in Equation (13), 1 in the first stage and 10 in the second stage, can be adjusted if a heavier penalty is desired.

Figures 8 and 9 compare the topology optimization results with and without the penalty function. The effect of using a penalty function is more apparent in the suspension arm example shown in Figure 9. Figure 9(a) is the topology optimization result without using a penalty function. It is observed that the gray element covers a significant part of the design domain, and the members are not distinct in this region. In Figure 9(b), the penalty function is used. Although there are still gray elements, the members are much clearer than the result in Figure 9(a).

(a) The result without penalty function

(b) The result with penalty function

Figure 8. The comparison of the topology optimization results with and without penalty function in the cantilever beam example

(a) The result without penalty function

(b) The result with penalty function

Figure 9. The comparison of the topology optimization results with and without penalty function in the suspension arm example

4.         The checkerboard problem

The “optimal topology” often contains a checkerboard pattern, which is the alternately solid and void elements as shown in Figure 4. This is a typical result in topology optimization using finite elements. It was believed that some sort of microstructure exists in these regions. Bendsfe [18] demonstrated that the checkerboard problem is related to the features of finite element approximation as a numerical phenomenon. This interpretation is now widely accepted.

At least four types of methods are proposed to prevent the checkerboard problem. In the “smoothing method,” image processing is used to smooth the output picture of topology optimization within the checkerboards. The use of higher-order finite elements is also suggested to avoid the checkerboard problem [28]. In the “patches method,” a “super-element” is introduced to the finite element formulation to damp the appearance of checkerboards [18]. Finally in the “filter method,” the modified design sensitivities of each iteration are used to prevent checkerboard patterns [22]. A survey of the checkerboard pattern problem was given by Sigmund and Petersson [19]. Of the four methods, the higher-order finite elements is probably the most convenient one. No external techniques are needed other than altering the element that is used to discretize the design domain to a higher-order finite element. Rodrigues and Fernandes [29] showed that the use of the nine-node element could prevent the occurrence of checkerboard pattern that occurred when using the four-node element. Díaz and Sigmund [28] and Jog and Haber [30] also showed that checkerboard pattern could be prevented efficiently when using eight- or nine-node elements.

Figure 4 shows the topology optimization result of the cantilever beam example using a four-node element. The finite element analysis solver used here is ANSYS, and the two-dimensional four-node element “Plane 42” provided by ANSYS [31] was used. It is observed that the checkerboard pattern appears in many parts of the “optimal structure” of this example. The same example is then solved again using the two-dimensional eight-node element “Plane 82” provided by ANSYS. The topology optimization result is shown in Figure 10. Comparing Figure 4 with Figure 10, the checkerboard patterns disappear resulting in a clear topology of the structure.

Figure 10. The topology optimization result of the cantilever beam example using eight-node element

Figure 11 shows the topology optimization result using four-node elements for the suspension arm example. The element distribution shown in Figure 11 are the elements in which the normalized material density is greater than 0.5. The checkerboard patterns are also distinct in this example. The topology optimization result when using eight-node elements is shown in Figure 9(b). The distributed elements shown in Figure 9(b) are also the elements in which the normalized material density is greater than 0.5. The checkerboard patterns disappear.

Figure 11. The topology optimization result of the suspension arm example using four-node element

These examples demonstrate that the use of higher-order finite elements can effectively prevent checkerboard patterns, even for the example with complicated design domain. Thus, the higher-order finite element is suggested and is actually used in the topology optimization process presented in this research.

5.         The continuity analysis of the topology optimization result

Topology optimization with penalty function attempts to generate the result with elements whose normalized material densities are either 0 or 1. But the “porous topology” [32], which is formed by elements with intermediate normalized material densities between 0 and 1, often occurs. In this paper, it is assumed that the material is isotropic, and materials with intermediate densities are not allowed.

Practically the most convenient way is to force the intermediate normalized material densities which are equal to or greater than a given threshold value of 1 directly. On the other hand, the densities that are less than the given threshold value are forced to be 0. The continuity of the structure should be of the most important consideration in deciding this threshold value. As shown in Figure 12, if the threshold value is high and a portion of porous topology A is discarded, a discontinuous structure (a structure with dangling ribs) will be generated, which may greatly affect the compliance of the structure. On the other hand, the continuity of the structure is not affected if the porous topology B is filtered out.

Figure 12. The illustration of the porous topology

A structure with continuous topology is stiffer than the discontinuous one, in other words, the compliance of the structure with continuous topology is lower. Thus, a proper threshold density value  is to be decided by comparing the compliances of the structure using different threshold density values.

Table 2 shows the compliance increase and the respective topology results of the cantilever beam example using different threshold density values. As shown in Table 2, when a discontinuity of the topology occurs (indicated by a circle) as the threshold value increases, there will be a corresponding “jump” in compliance. In this research, if the increase in compliance is more than 10% at certain threshold density value, the threshold density value is set to be . In the cantilever beam example, it was decided to set  at 0.2. In this way a proper threshold density value  can be selected to filter the porous topology without generating a discontinuous structure.

Table 2. The topology of different threshold density

Topology

Topology

0.1

0.6

Compliance increase: 0%

Compliance increase: 98.6%

0.2

0.8

Compliance increase: 34.9%

Compliance increase: 1937.4%

After the continuity analysis, it is observed that trivial solids (such as islands or small salient) and trivial voids (the voids containing few elements whose normalized material densities are 0) appear in the topology result. Trivial solids and trivial voids are useless in mechanics and will complicate the implementation of the shape optimization. Thus, the trivial solids should be filtered out before implementing shape optimization. In this research two “filters” are implemented to filter out the trivial solids and trivial voids in the topology result after the continuity analysis.

The trivial solid is defined as the solid containing a group of elements that do not connect with other groups, or connect with other groups at only one node. Figure 13 is the result of cantilever beam example after applying the trivial solid filter. There can be many voids in the topology optimization result. In this research, if the number of elements contained in a void is less than 1% of the total number of elements multiplied by , the ratio of the amount of usable material, it is defined as a trivial void. Figure 14 shows the topology result of the suspension arm example after the trivial void filter. The checkerboard like trivial voids shown in Figure 14 are filtered out.

Figure 13. The result of the cantilever beam example after the trivial solid filter

Figure 14. The result of suspension arm example after the trivial void filter

6.         Two- and three-dimensional topology optimization examples

After a high-quality topology optimization result is obtained in the Phase I (topology generation process), it becomes much easier to interpret the result in the Phase II (topology interpretation process), then complete the shape and size optimization in the Phase III (detailed design phase). Table 3 shows 14 two-dimensional structural topology optimization examples commonly seen in literature. The topology optimization procedure developed in this research worked well for all 14 examples. The results are also presented in Table 3.

Table 3. Two-dimensional topology optimization examples

 

Examples

Topology optimization results in literatures

Topology optimization results

1

Bendsfe & Rodrigues [12]

2

Bendsfe & Rodrigues [12]

3

Bendsfe & Rodrigues [12]

4

Olhoff et al. [13]

5

Bremicker et al. [8]


6

Bremicker et al. [8]

7

Bremicker et al. [8]

8

Chirehdast et al. [9]

9

Chirehdast et al. [9]

10

Chirehdast et al. [9]

11

Tang & Chang, [17]

12

Youn & Park, [15]

13

Youn & Park, [15]

14

Youn & Park, [15]

Topology optimization has been extended to three-dimensional continuum structures in recent years. Cherkaev and Palais [33] considered the topology optimization of three-dimensional axisymmetric elastic structures. In their work, the method of finding optimal bounds upon the effective properties of a composite is used. Díaz and Lipton [34] presented a full relaxation strategy to find the optimal topology of minimum compliance for three-dimensional elastic structures subjected to the amount of available material. Olhoff et al. [35] and Jacobsen et al. [36] used three-dimensional microstructures developed by Gibianski and Cherkaev [37] to get the optimum topology for three-dimensional structures.

Beckers [38] developed a dual method in which the design variables can only be 1 (present) or 0 (absent) to deal with the topology optimization of continuous structure. In his work, two- and three-dimensional problems are solved. The examples of three-dimensional structural topology optimization using the density distribution method can be found in the papers by Yang and Chen [39] and Hsu et al. [16].

For three-dimensional structures, it is very difficult to interpret the topology optimization result manually. Therefore a generalized and automated process is even more important. The generalized topology optimization process shown in Figure 1 and the discretized optimization model Equation (6) can be applied to the topology optimization of three-dimensional structures without any modifications. Table 4 shows 6 three-dimensional structural topology optimization examples that are commonly seen in the literature. The same two-stage penalty function presented in Section 3 is used to penalize the intermediate normalized material density, and the higher order 20-nodes cubic element is used in order to avoid the checkerboard problem. The same continuity analysis and trivial solids and voids filters discussed in Section 5 are used to fine-tune the final topology. Both the topology optimization results reported in the literature and in this research are listed in Table 4 for comparison.

Table 4. Three-dimensional topology optimization examples

 

Examples

Topology optimization results in literatures

Topology optimization results

1

Olhoff et al. [35], Jacobsen et al. [36], Beckers, [38]

2

Olhoff et al. [35], Jacobsen et al. [36]

3

Olhoff et al. [35], Jacobsen et al. [36]

4

Olhoff et al. [35], Jacobsen et al. [36], Beckers, [38]

5

Díaz, & Lipton, [34]

6

Díaz, & Lipton, [34]

7.         Conclusion

A generalized topology optimization process is presented in this paper. Several fundamental issues on the quality of the topology optimization result are considered in order to achieve a clear topology optimization result. They are mesh-dependence, elements with intermediate densities, the checkerboard problem, and porous topology.

A convergence analysis is implemented to find the minimum number of elements that does not significantly affect the compliance. A two-stage penalty function strategy is developed in this research to determine a proper topology optimization result in the first stage using a light penalty, then heavily penalize the elements with intermediate densities in the second stage. Higher order finite element is used, which is convenient to the users, and the checkerboard pattern is avoided effectively. Finally a continuity analysis is used to deal with the porous topology, and two filters are implemented to filter out the trivial solids and voids.

The whole process is generalized to two- and three- dimensional structure and can be fully automated. Fourteen two-dimensional and six three-dimensional examples are used to demonstrate the effectiveness of this process.

References

[1]     Bendsfe, M. P., Kikuchi, N. (1988) Generating Optimal Topologies in Structural Design Using Homogenization Method. Computer Methods in Applied Mechanics and Engineering, 71, 197-244.

[2]     Mlejnek, H. P., 1992, Some Aspects of the Genesis of Structures. Structural Optimization, 5, 64-69.

[3]     Yang, R. J., Chuang, C. H. (1994) Optimal Topology Design Using Linear Programming. Computers & Structures, 52, 265-275

[4]     Xie. Y. M., Steven, G. P. (1993) A Simple Evolutionary Procedure for Structural Optimization. Computers & Structures, 49, 885-896.

[5]     Rozvany, G. I. N. (2001) Aim, Scope, Methods, History and Unified Terminology of Computer-Aided Topology Optimization in Structural Mechanics. Structural and Multidisciplinary Optimization, 21, 250-252.

[6]     Yang, R. J., Chahande, A. I. (1995) Automotive Applications of Topology Optimization. Structural Optimization, 9, 245-249.

[7]     Chen, T. -Y., Wu, S. –C. (1998) Multiobjective Optimal Topology Design of Structures. Computational Mechanics, 21, 483-492.

[8]     Bremicker, M., Chirehdast, M., Kikuchi, N., Papalambros, P. Y. (1991) Integrated Topology and Shape Optimization in Structural Design. Mechanics of Structures and Machines, 19, 551-587.

[9]     Chirehdast, M., Gea, H-C., Kikuchi, K., Papalambros, P. Y. (1994) Structural Configuration Examples of an Integrated Optimal Design Process. Journal of Mechanical Design, 116, 997-1004.

[10] Lin, C. Y. Chao, L. S. (2000) Automated Image Interpretation for Integrated Topology and Shape Optimization. Structural and Multidisciplinary Optimization, 20, 125-137.

[11] Papalambros, P., Chirehdast, M., (1990) An Integrated Environment for Structural Configuration Design. International Journal for Numerical Methods in Engineering, 41, 1417-1434.

[12] Bendsfe, M. P., Rodrigues, H. C. (1991) Integrated Topology and Boundary Shape Optimization of 2-D Solid. Computer Methods in Applied Mechanics and Engineering, 87, 15-34.

[13] Olhoff, N., Bendsfe, M. P., Rasmussen, J. (1991) On CAD-Integrated Structural Topology and Design Optimization. Computer Methods in Applied Mechanics and Engineering, 89, 259-279.

[14] Kumar, A. V., Gossard, D. C. (1996) Synthesis of Optimal Shape and Topology of Structures. Journal of Mechanical Design, 118, 68-74.

[15] Youn, S. K., Park, S-H. (1997) A Study on the Shape Extraction Process in the Structural Topology Optimization Using Homogenized Material. Computers & Structures, 62, 527-538.

[16] Hsu, Y. L., Hsu, M. S., Chen, C. T. (2001) Interpreting Results from Topology Optimization Using Density Contours. Computers & Structures, 79, 1049-1058.

[17] Tang, P. S., Chang, K. H. (2001) Integration of Topology and Shape Optimization for Design of Structural Components. Structural and Multidisciplinary Optimization, 22, pp.65-82.

[18] Bendsfe, M. P. (1995) Optimization of Structural Topology, Shape, and Material, Springer-Verlag Berlin Heidelberg.

[19] Sigmund, O., Petersson, J. (1998) Numerical Instabilities in Topology Optimization: A Survey on Procedures Dealing with Checkerboards, Mesh-dependencies and Local Minima. Structural Optimization, 16, 68-75.

[20] Ambrosio, L., Buttazzo, G.. (1993) An Optimal Design Problem with Perimeter Penalization. Calculus of Variations and Partial Differential Equations, 1, 55-69.

[21] Haber, R. B., Jog, C. S., Bendsfe, M. P. (1996) A New Approach to Variable-Topology Shape Design Using a Constraint on Perimeter. Structural Optimization, 11, 1-12.

[22] Sigmund, O. (1994) Design of material structures using topology optimization, Ph.D thesis, Department of Solid Mechanics, Technical University of Denmark.

[23] Sigmund, O. (1997) On the Design of Compliant Mechanisms Using Topology Optimization. Mechanics of Structures and Machines, 4, 493-524.

[24] Niordson, F. I. (1983) Optimal Design of Plates with a Constraint on the Slope of the Thickness Function. International Journal of Solids and Structures, 19, 141-151.

[25] Petersson, J., Sigmund, O. (1998) Slope Constrained Topology Optimization. International Journal for Numerical Methods in Engineering, 41, 1417-1434.

[26] Zhou, M., Shyy, Y. K., Thomas, H. L. (2001) Checkerboard and Minimum Member Size Control in Topology Optimization. Structural and Multidisciplinary Optimization, 21, 152-158.

[27] Borrvall, T., Petersson, J. (2001) Topology Optimization Using Regularized Intermediate Density Control. Computer Methods in Applied Mechanics and Engineering, 190, 4911-4928.

[28] Díaz, A. R., Sigmund, O. (1995) Checkerboard patterns in layout optimization. Structural Optimization, 10, 40-45.

[29] Rodrigues, H., Fernandes, P. (1995) A Material Based Model for Topology optimization of Thermoelastic Structures. International Journal for Numerical Methods in Engineering, 38, 1951-1965.

[30] Jog, C. S., Haber, R. B. (1996) Stability of Finite Element Models for Distributed-parameter Optimization and Topology Design. Computer Methods in Applied Mechanics and Engineering, 130, 203-226.

[31] ANSYS Inc. (1999) PLANE42 2-D Structural Solid. ANSYS Element Reference, 11 edition, ANSYS Inc.

[32] Rozvany, G. I. N., Zhou, M., Birker, T., (1992) Generalized Shape Optimization without Homogenization, Structural Optimization, 4, 250-252.

[33] Cherkaev, A., Palais R. (1996) Optimal design of three-dimensional axisymmetric elastic structures. Structural Optimization, 12, 33-45.

[34] Díaz, A., Lipton, R. (1997) Optimal material layout for 3D elastic structures, Structural Optimization, 13, 60-64.

[35] Olhoff, N., Rønholt, E., Scheel, J. (1998) Topology optimization of three-dimensional structures using optimum microstructures. Structural Optimization, 16, 1-18.

[36] Jacobsen, J. B., Olhoff, N., Rønholt, E. (1998) Generalized shape optimization of three-dimensional structures using materials with optimum microstructures. Mechanics of Materials, 28, 207-225.

[37] Gibiansky, L. V., Cherkaev, A. V. (1987) Microstructures of composites of extremal rigidity and exact estimates of provided energy density. A. F. Ioffe Physico-Technical Institute, 1115. English translation to appear in: Kohn, R. (1994) Topics in the mathematical modeling of composite materials, Birkhauser.

[38] Beckers, M. (1999) Topology optimization using a dual method with discrete variables. Structural Optimization, 17, 14-24.

[39] Yang R. J., Chen, C. J., 1996 Stress-based Topology Optimization, Strucutral Optimization, 12, 98-105.