A nodal variable ESO (BESO) method for structural topology optimization

Finite Elements in Analysis and Design 86 (2014) 34–40

Contents lists available at ScienceDirect

Finite Elements in Analysis and Design journal homepage: www.elsevier.com/locate/finel

A nodal variable ESO (BESO) method for structural topology optimization Fei Zhao Key Laboratory of Electronic Equipment Structure Design Ministry of Education, Xidian University, Xi'an 710071, China

art ic l e i nf o

a b s t r a c t

Article history: Received 7 November 2013 Received in revised form 23 March 2014 Accepted 29 March 2014 Available online 24 April 2014

A nodal variable ESO (BESO) method is proposed for topology optimization of continuum structures in this paper. The initial discrete-valued topology optimization problem is established as an optimization problem based on continuous design variables by employing a material density field into the design domain. The density field, with the Shepard family of interpolation, is mapped on the design space defined by a finite number of nodal density variables. The employed interpolation scheme has an explicit form and satisfies non-negative and range-restricted properties required by a physically significant density interpolation. It has the ability to deal with more complex spatial distribution of the material density within an individual element, as compared with the conventional elementwise design variable methods. Numerical examples demonstrate the validity and effectiveness of the improved method. & 2014 Elsevier B.V. All rights reserved.

Keywords: Topology optimization ESO (BESO) method Nodal density variables Density interpolation Shepard interpolation

1. Introduction It is of great importance for the development of new methods to create the best possible topology or layout for given design objectives and constraints at a very early stage of the design process (the conceptual and project definition phase). Thus, over the past two decades, substantial efforts of fundamental research have been devoted to the field of development of efficient and reliable procedures for solution of such problems [1]. During this period, the researchers have been mainly occupied with two different kinds of topology design processes. The first type can be regarded as material distribution models [2,3]. And another type is boundary description models, including the topology description function-based method, the level set-based method and its variants [4]. Meshfree methods have also been applied in a few topology optimization studies [5,6]. The material distribution models have become the most popular formulations in topology optimization. Two well-established approaches based on the material distribution concept are the homogenization method (and its variant, the Solid Isotropic Material with Penalization (SIMP) method) and the ESO (the evolutionary structural optimization) method [7,8]. In the homogenization method, the layout of materials is composed of porous materials with parametrically defined microstructures inside the design domain. An isotropic material model is employed in the SIMP method with power law penalization of the elastic constants. And the element densities are to be optimized during the optimization period. The ESO methods

E-mail address: [email protected] http://dx.doi.org/10.1016/j.finel.2014.03.012 0168-874X/& 2014 Elsevier B.V. All rights reserved.

introduce finite changes in a design on the basis of certain heuristic criteria, which may not be based on sensitivities. The basic idea of ESO method is based on a simple and empirical concept that a structure evolves towards an optimum by slowly removing elements with lowest stresses [9]. In order to maximize the stiffness of the structure, stress criterion was employed with elemental strain energy criterion in the ESO method. An extension of ESO method is a Bi-directional evolutionary structural optimization (BESO) method. New elements can be added in the locations next to those elements with highest stresses in the BESO method. It is worth noting that both of the two approaches are conventionally formulated by using elementwise design variables. It means that the FEM method is used to describe material distribution as well as the displacement fields inside the design domain simultaneously. The material properties are assumed constant within each element. The ESO (BESO) method has got fairly general success due to its efficient implementation and conceptual simplicity. These methods have been applied to a wide variety of industrial applications and have been also extended to design problems involving multimaterials and meshless method [10,11]. Although using element-based design variables seems natural in the model of optimal material distribution problems, it still presents some difficulties. Such kind of topology optimization methods often suffer from some well-known numerical instabilities, in particular the checkerboard pattern and mesh dependency, even when higher-order finite elements are used. It is noted that these numerical instabilities are partly attributed to the use of elementwise constant material density distribution [12]. The reason of checkerboards phenomenon is the inappropriate computational modeling of strains/stresses, which leads to over-stiffness evaluation

F. Zhao / Finite Elements in Analysis and Design 86 (2014) 34–40

of local regions with low-order finite elements. The mesh-dependence refers to the phenomenon that a finer finite element discretization will result in a more complex topological description of the structure. In this way, one and the same design with different mesh refinements will easily cause series local minima [13]. Nowadays, many researchers have made great efforts to overcome the above numerical instabilities for topology optimization problems. For example, a smoothing algorithm using the surrounding elements’ reference factors has been proposed to prevent checkerboarding. However, this smoothing procedure depends on the assigned mesh and fails to overcome the mesh-dependency problem. In order to overcome mesh-dependency, an additional restriction by perimeter control or filter has been introduced into the BESO method [10]. A sequential element rejection and admission (SERA) method is suggested in which the void element is replaced by a soft element with a very low density [14]. For the aforementioned numerical difficulties to be solved, many attempts have been made to improve the topology optimization problem by using nodal density variables of finite elements. Matsui and Terada were the first who studied the nodal density approximant for topology optimization problems by using finite element shape functions [15]. They proposed a so-called CAMD (continuous approximation of material distribution) model based on the homogenization method. In order to achieve checkerboard-free designs, they adopted a continuous approximation of material distribution model, which obtained the continuous design variables by interpolating the discrete nodal variables with the bilinear shape function of four-node quadrilateral (Q4) finite elements. A Heaviside projection method is employed to prescribe a minimum allowable length scale of solid/void phase to prevent sharp corners, onenode hinge, and over-thin members [16]. A Q4/Q4 implementation of the SIMP method is used to suppress checkerboard phenomenon with nodal density variables [17]. It is reported that the checkerboard patterns can thus be effectively avoided. However, a side effect, namely the so-called islanding phenomenon was introduced, especially when coarse meshes are used. Here, ‘islanding’ is referred to as the material distribution pattern with a layering-up manner in certain local areas. Paulino and Le [18] proposed a kind of hybrid low-order finite elements, in which the nodes of design variables are inconsistent with the nodes of the displacement vector. An internal averaging technique was further applied to suppress the islanding effect. Kang and Wang [19,20] discussed the merits of standard Shepard interpolation in preserving the physical meaning of density variables for topology optimization of structures. They proposed a Shepard function interpolation combined with higher-order elements for the elimination of numerical checkerboards. This study focuses on the proposal of a new topology optimization methodology using an alternative nodal variable model in the framework of ESO/BESO methods. Obviously, Shepard interpolation is suited for interpolating the material density field by using nodal density values. Therefore, a nodal design variable ESO/BESO methods for topology optimization is founded based on Shepard interpolation. 2. Nodal density interpolation using Shepard interpolation An appropriate interpolation formulation based on Shepard interpolation will be employed to create the nodal-value-based material density field in higher-order finite elements in this section. 2.1. The Shepard interpolation method First of all, the Shepard interpolation method is introduced in this section. Let ϕi , i ¼ 1; ⋯; n, denote a set of n non-negative data

35

values at the associated sampling points xi ¼ ðX i ; Y i ; Z i Þ, where X i ; Y i ; Z i define the ith point location in a Cartesian coordinate system. The interpolation scheme given by Shepard [21] and Brodlie et al. [22] is explicitly expressed as follows: n

ΦðxÞ ¼ ∑ ψ i ðxÞϕi

ð1Þ

i¼1

where the normalized interpolation function ψðxÞ can be formulated as ψ i ðxÞ ¼

ϖ i ðxÞ

ð2Þ

n

∑ ϖ j ðxÞ

j¼1

ϖ i ðxÞ ¼

1

ð3Þ

2

ddis ðxÞ

ddis ðxÞ ¼ jx  xi j2

ð4Þ

with i ¼ 1; ⋯; n and ddis is the Euclidean distance between points x and xi. It is easy to see that the interpolation function ψ i ðxÞ possesses the non-negative and range-bounded properties (0 rψ i ðxÞ r 1), and a normality (unity) property (∑ni¼ 1 ψ i ðxÞ ¼ 1). There is no conditional requirement that the Shepard method must satisfy the interpolation δij condition. Such kind of interpolation scheme is easy to implement in both two dimensional and three dimensional spaces straightforward. There is another important property of this interpolation method. It is worth while to note that the interpolation values at any points are bounded between minimum and maximum values of the sampling dataset, that is minfϕi g r ΦðxÞ r maxfϕi g i

i

ð5Þ

This characteristic is crucial for the material density interpolation in topology optimization applications. The reason is that the relative density is a ½0; 1 bounded quantity and any value not falling within this range is physically meaningless.

2.2. Density distribution based on the Shepard interpolation An elementwise nodal density based on the Shepard interpolation strategy is employed in this study. It should be noted that independent sets of points can be used for displacement and material density field interpolation simultaneously. Three kinds of possible arrangements of displacement points and material density points are displayed in Fig. 1. In this paper, for practical reasons, we focus on the case that the material density points coincide with the FEM nodes in a higher-order finite element. In particular, we employ the eight-node quadratic quadrilateral (Q8) element for discretization of the design domain, with the density point arrangement as shown in Fig. 1(c) [19]. It should be noted that the Shepard interpolation based on an eight-node quadratic quadrilateral (Q8) element is incorporated into the ESO (BESO) method. And it can be used to suppress checkerboard patterns and islanding phenomenon. Moreover, with the same mesh resolution, this method provides relatively finer topological details than other approaches and yields structural boundaries smoother than conventional nodal variable-based approach. Sometimes, the proposed method also results in a lower mean compliance. For more details about the merits of this method, please refer to [19,20]. The regular quadratic shape functions for higher-order elements are used for the displacement field. Therefore, the relative density at an arbitrary point ‘x’ within an element

36

F. Zhao / Finite Elements in Analysis and Design 86 (2014) 34–40

Displacement

Material Density

Both Fig. 1. Three examples of arrangement of displacement points (FEM nodes) and material density points: (a) D4/M4; (b) D4/M9; (c) D8/M8.

Shepard interpolation scheme as presented in Eq. (6). Therefore, the optimization problem can be clearly defined as follows [10]:

can be expressed as m

ρðxÞ ¼ ∑ ψ i ðxÞρi

ð6Þ

i¼1

where m ¼8 is the number of elemental nodes in this study. ρi is the density value of the ith elemental node. This interpolation scheme is able to reproduce a constant density field within the element. Another prominent characteristic of this interpolation method is satisfying the Kronecker delta criteria. It means that the essential boundary conditions can be imposed directly without extra techniques. It is also easy to know from the bounded property of Shepard interpolation that 0 r ρðxÞ r1 holds if 0 r ρi r 1 (i ¼ 1; ⋯; m). In addition, the property ψ i ðxÞ Z 0 also guarantees that the derivative of the density with respect to the nodal density variable will be always nonnegative. The desirable continuity of the approximation can be easy to inherit from its interpolation function of high continuity. The derivatives of ρðxÞ with respect to ρi is non-negative; this satisfies the inequality condition ∂ρðxÞ=∂ρi Z0 for any point x, which is important in ensuring the monotonic behavior and application of gradientbased algorithms. The above density approximant is included in the problem formulation rather than in the post-processing procedure. It means that the density approximant as a heuristic filtering scheme is used to enhance the smoothness of nodal densities. Obviously, the problem formulation may be similar to the filtering or projection methods in topology optimization algorithm. Thus, the elemental stiffness matrix can be expressed as Z Ke ¼ BT DðxÞBdΩ ð7Þ

1 Minimize : C ¼ uT Κu 2 Ne

Subject to : V 0  ∑ Ve ρe ¼ 0 e¼1

ρe ¼ ρmin or 1

ð9Þ

It is worth to note that C ¼ ð1=2ÞuT Κu is the mean compliance rather than the compliance, C ¼ uT Κu used in the current ESO (BESO) method. K is the global stiffness matrix of the structure, and u is the displacement vector. Ve and V0 are the volumes of an individual element and the prescribed structural volume, respectively. The binary design variable ρe denotes the density of the eth elemental material density. In the initial ESO (BESO) method, a solid element is completely removed from the design domain that could result in theoretical difficulties in topology optimization. Obviously, it could be rather irrational when the design variable (an element) is directly eliminated from the topology optimization problem. Hence, a small value of ρmin is used to denote the void materials. Under normal conditions, ρmin is equal to 0.0001 for avoiding numerical difficulties associated with zero density elements. The objective function can be further stated as 1 1 Ne C ¼ uT Κu ¼ ∑ uTe Κ e ue 2 2e¼1

ð10Þ

where ue is the elemental displacement vector and Ke is the elemental stiffness matrix.

Ωe

where B is the strain displacement matrix, D is the elasticity constant matrix, and Ωe is the volume occupied by the eth element. Here, a 4  4 Gauss numerical integration is employed for computing the elemental stiffness matrix.

The sensitivity of Young's modulus of the material with respect to eth design variables is expressed as ∂EðxÞ ¼ ψ e ðxÞE0 ∂ρe

3. Formulation of topology optimization and numerical implementation The ESO method is based on the simple idea that the optimal structure can be produced by gradually removing the ineffectively used material from the design domain. Bi-directional evolutionary optimization (BESO) is an extension of that idea which allows for efficient material to be added to the structure at the same time as the inefficient one to be removed. The relation between the Young modulus and the material density at point x is thus defined as [20] Eðxe Þ ¼ E0 ρðxe Þ

4. Sensitivity analysis and sensitivity numbers

ð8Þ

where E0 denotes Young's modulus for solid material. Suppose that Ne design variable points ρe (e ¼ 1; ⋯; N e ) are considered in the design domain, the material density field ρðxÞ in above equation can be related to the design variable points with

ðe ¼ 1; …; N e Þ

ð11Þ

The sensitivity of the objective function with respect to design variables can be easily derived by the adjoint method [23] and can be defined as ∂C 1 ∂K 1 Ne ∂K e ¼  uT u ¼  ∑ uTe ue ∂ρe 2 ∂ρe 2 e ¼ 1 ∂ρe Substituting Eq. (7) into above equation yields Z  ∂C 1 Ne ∂DðxÞ ¼  ∑ uTe BT BdΩ ue ∂ρe 2e¼1 ∂ρe Ωe

ð12Þ

ð13Þ

where the sensitivity of the elastic matrix ∂DðxÞ=∂ρe can be easily derived from Eq. (11). The discrete design variables are used to optimize a structure in the ESO (BESO) method. It should be noted that only two bound materials are allowed in the design. Therefore, the relative ranking of elemental sensitivities for both solid and soft elements can be

F. Zhao / Finite Elements in Analysis and Design 86 (2014) 34–40

expressed as ( 1 T ∂K 2u ∂ρe u αe ¼ 0

37

Start

ρe ¼ 1 ρe ¼ ρmin

ð14Þ Define design structure

where αe is named as the sensitivity number for the eth element. This equation indicates that the sensitivity values of solid elements and soft elements are different. The values of solid elements and soft elements are equal to the elemental strain energy and zero, respectively. According to the above material interpolation approach, Young's modulus of soft elements also becomes zero. Due to the above reasons, when elements are equivalent to void elements and can be completely removed from the design domain. However, the ESO (BESO) algorithm using above sensitivity numbers is hard to convergent because the sensitivity numbers are calculated based on the different status of elements (1 or xmin ). Computational experience has shown that averaging the sensitivity number with its historical information is an effective way to overcome the drawback. It can be simply defined as 1 αe ¼ ðαe;k þ αe;k  1 Þ 2

Define ESO/BESO parameters

Analyze the structure and calculate sensitivity numbers

Construct a new structure

Averaging the sensitivity number

ð15Þ

where k is the iteration number. Then let αe;k ¼ αe ; thus the modified sensitivity number considers the sensitivity information in the previous iterations.

No

Is volume constraint satisfied ? Yes

5. Optimality algorithm In the ESO (BESO) method, the strain energy densities of all elements should be equal [7]. It means that the elements with higher strain energy density should have increased ρe . On the other hand, the elements with lower strain energy density should have decreased ρe [8]. Hence, the optimality criteria can be defined as: strain energy densities of solid materials are always higher than those of soft materials. Therefore, the update scheme can be described as that of the design variables ρe by changing from 1 to ρmin for elements with lowest sensitivity numbers, and, on the contrary, from ρmin to 1 for elements with highest sensitivity numbers. During the optimization progress, the threshold of sensitivity number can be easily determined by the target volume for the next iteration and relative ranking of sensitivity numbers [9]. Fig. 2 shows the flowchart of the proposed ESO (BESO) method.

No

Converged ? yes End

Fig. 2. An ESO (BESO) flowchart.

L

Design Domain

F

H

6. Numerical examples In this section, two typical numerical examples in the area of topology optimization for structural mean compliance are used to demonstrate the effectiveness of the proposed ESO (BESO) method. In the numerical implementation of the proposed method, Q8 elements are employed for the displacement field analysis. The ‘Shepard function’ is included in the mathematical formulation of the topology optimization problem. No additional treatment (such as sensitivity filtering, penalty exponent material interpolation scheme and density averaging) for numerical instability control is applied in the optimization process. For the sake of numerical simplicity, an artificial material model is adopted in both examples. The material properties are defined as follows: Young's modulus for the full-solid region is E0 ¼ 1, for weak material to fill the void area is E ¼ 0:0001, and Poisson's ratio for all material inside the design domain is nu ¼ 0:3. During the design optimization, the units for the artificial material model can be defined as flexible, but all the units are required to be unchanged.

Fig. 3. Design domain of the cantilever beam case.

6.1. Examples for a cantilever beam case A long cantilever shown in Fig. 3 is selected as a test example because it involves a series of bars broken during the evolution process of the ESO (BESO) topology optimization. The geometrical dimension of the design domain with a length and width ratio of L:H is equal to 2:1, as shown in Fig. 3. The left side of the domain is fixed as the Dirichlet boundary, and a downward force (F¼1) applied in the middle of the right border. The right side is treated as a nonhomogenous Neumann boundary. The objective function is to minimize the mean compliance, and the constraint is to limit the maximum material usage less than 50%. A 180  90 mesh is used for discretization of the design domain. The evolutionary ratio is set to be ER ¼ 2%.

38

F. Zhao / Finite Elements in Analysis and Design 86 (2014) 34–40

After 48 iterations, the final topology optimization result is obtained. The structural mean compliance is C ¼ 31:2357. Fig. 4 shows the corresponding plots of the topologies at different design stages (e.g. v¼ initial, 0.9, 0.8, 0.7, 0.6 and final). The initial guess of the design variables is shown as ρe ¼ 1 (e ¼ 1; …; N e ) in Fig. 4(a). The initial guess does not need to be assumed with extra holes at suitable locations. Under normal conditions, the initial guess work needs rich experience and skillful capacity. It is easy to implement the derivative algorithm without extra complex auxiliary techniques. Fig. 4(b)–(e) shows the plots at different iteration stages. Fig. 4(f) displays the ultimate result. The final design outcome is similar to that of the well known method in relative literatures. There are distinct topology contours with smooth design boundaries of solid or void interfaces. The new method is more capable in describing the local geometric details of topology layout. Finer meshes are more capable of describing local substructures during the design. This is a heuristic filtering scheme to enhance the smoothness of material densities rather than in the post-processing procedure. In the meantime, the islanding phenomenon can be prevented with higher-order approximation; because using higher-order displacement shape functions will substantially improve the accuracy of the mechanics analysis in case of nonuniform material property. Without intermediate densities and checkerboard phenomena is another merit of the proposed method. Fig. 5 shows the evolutionary histories of the mean compliance and volume fraction for the proposed method. It indicates that the mean compliance increases initially as the volume fraction decreases, then it converges to an almost constant value after the objective volume is achieved. The final structure is convergent to a stable topology after 48 iterations. The mean compliance of the proposed method increases, with occasionally abrupt jumps (due to breaking up of some bars), during which the total volume gradually decreases. After about 33 iterations the volume fraction reaches its target of 50%. In subsequent optimization process, while the volume remains unchanged the mean compliance gradually converges to a constant value.

6.2. Examples for an MBB beam case To demonstrate the proposed method, a topology optimization problem for a miniature bending beam structure is shown in Fig. 6. The beam length and width size ratio is L : H ¼ 6 : 1. The bottom left corner is a hinged support, and the bottom right corner supports as a roller. A vertical force (F ¼ 1) is applied at the center

Fig. 5. Histories of mean compliance and volume fraction.

F Design Domain

L

Fig. 6. Design domain of the MBB beam case.

Fig. 4. Topology results of the proposed method.

H

F. Zhao / Finite Elements in Analysis and Design 86 (2014) 34–40

point of the lower side. The objective function is to minimize the mean compliance, and the allowed material usage is limited to 50%. A 360  60 mesh and evolutionary ratio of ER ¼ 2% are used to implement the example, respectively. The optimal solution is converged after 68 iterations. The corresponding structural mean compliance is C ¼ 47:1385. The different design stages (e.g. v¼ initial, 0.9, 0.8, 0.7, 0.6 and final) of the corresponding distributions are demonstrated and revealed in Fig. 7. The corresponding density fields contain only two types of values. The first is ρ ¼ 1 for solid regions with full solid materials. The second is ρ ¼ 0:0001 used to determine the void fields with weak materials. It is an inevitable result that the intermediate densities phenomenon is overcome straightforward. Fig. 7(a) shows the initial layout distribution of uniform density variables (ρ ¼ 1) topology for the whole design domain. In Fig. 7 (b)–(e), the corresponding optimization procedure is shown one after another. The eventual optimal result is displayed in Fig. 7(f). It is unobvious difference between the ultimate solution in the proposed method and the traditional methods in relevant literatures. The concise interface and smooth boundary can be obtained obviously. Using higher-order displacement shape functions will substantially enforce the accuracy of the mechanics analysis in case of nonuniform material property. The islanding phenomenon can be suppressed naturally by using higher order interpolation because of better predicting the stiffness with layering material distribution. The initial topology is one of most simple layouts with uniform distributions. It is easy to carry out more complicate topological changes.

39

Fig. 8. Histories of mean compliance and volume fraction.

The histories of the mean compliance and volume are displayed in Fig. 8. It can be seen that a number of 36 iterations are used to compress the volume back to satisfy the constraint at first. The mean compliance increases initially as the volume fraction decreases, with occasionally abrupt jumps (because of breaking up of some bars), then it is convergent to an almost constant value. So the rest iterations are consumed for further shape variations. It is obvious that this process occupies nearly half of the computational procedure, but this effort is essential as it is used to gradually modify the distribution of the material inside the design domain until the optimal criteria is satisfied. The overall iterative history discovers that the initial topology process is used to yield a coarse topology layout of the material distribution, and the shape optimization is applied to modify the specified structural performance around local regions.

7. Conclusions Using a nodal design variable method of topology optimization based on Shepard interpolation of material density, a new ESO (BESO) methodology for topology optimization of continuum structures is investigated. The Shepard interpolation is applicable for physically meaningful density interpolation. Therefore, the density field is mapped on nodal design variables based on the Shepard family of interpolations. The initial topology optimization problem is transformed into a finite-dimensional continuous optimization problem, in which nodal density values were considered as the design variables. Numerical examples show that the present approach has a good convergence and performs well in resolving fine details of the optimal structural layout. At the same time, without any additional instability compress technique skill, the present method can effectively suppress the checkerboard patterns as well as the islanding phenomenon.

Acknowledgment

Fig. 7. Topology results of the proposed method.

This work was supported by the National Natural Science Foundation of China under Grant nos. 51035006, 51305321 and 51105290.

40

F. Zhao / Finite Elements in Analysis and Design 86 (2014) 34–40

References [1] M.P. Bendsøe, N. Kikuchi, Generating optimal topologies in structural design using homogenization method, Comput. Methods Appl. Mech. Eng. 71 (1988) 197–224. [2] H.A. Eschenauer, Topology optimization of continuum structures: a review, Appl. Mech. Rev. 54 (4) (2001) 331–389. [3] G.I.N. Rozvany, A critical review of established methods of structural topology optimization, Struct. Multidiscip. Optim. 37 (3) (2009) 217–237. [4] K.A. James, E. Lee, R. Joaquim, R.A. Mattins, Stress-based topology optimization using an isoparametric level set method, Finite Elem. Anal. Des. 58 (2012) 20–30. [5] H.Y. Xu, L.W. Guan, X. Chen, L.P. Wang, Guide-weight method for topology optimization of continuum structures including body forces, Finite Elem. Anal. Des. 75 (2013) 38–49. [6] F. Zhao, A meshless Pareto optimal method for topology optimization, Eng. Anal. Bound. Elem. 37 (2013) 1625–1631. [7] Y.M. Xie, G.P. Steven, Evolutionary Structural Optimization, Springer, London, 1997. [8] X. Huang, Y.M. Xie, A new look at the ESO/BESO optimization methods, Struct. Multidiscip. Optim. 35 (1) (2008) 89–92. [9] X. Huang, Y.M. Xie, A further review of ESO type methods for topology optimization, Struct. Multidiscip. Optim. 41 (5) (2010) 671–683. [10] X. Huang, Y.M. Xie, Bi-directional evolutionary topology optimization of continuum structures with one or multiple materials, Comput. Mech. 43 (3) (2009) 393–401. [11] J. Zheng, S.Y. Long, G.Y. Li, The topology optimization design for continuum structures based on the element free Galerkin method, Eng. Anal. Bound. Elem. 34 (2010) 666–672.

[12] O. Sigmund, J. Petersson, Numerical instabilities in topology optimization: a survey on procedures dealing with checkerboards, mesh-dependencies and local minima, Struct. Multidiscip. Optim. 16 (1998) 68–75. [13] A.R. Diaz, O. Sigmund, Checkerboard patterns in layout optimization, Struct. Multidiscip. Optim. 10 (1995) 40–45. [14] G.I.N. Rozvany, O.M. Querin, Combining ESO with rigorous optimality criteria, Int. J. Veh. Des. 28 (2002) 294–299. [15] K. Matsui, K. Terada, Continuous approximation of material distribution for topology optimization, Int. J. Numer. Methods Eng. 59 (2004) 1925–1944. [16] J.K. Guest, J.H. Prévost, T. Belytschko, Achieving minimum length scale in topology optimization using nodal design variables and projection functions, Int. J. Numer. Methods Eng. 61 (2004) 238–254. [17] S.F. Rahmatalla, C.C. Swan, A Q4/Q4 continuum structural topology optimization implementation, Struct. Multidiscip. Optim. 27 (2004) 130–131. [18] G.H. Paulino, C.H. Le, A modified Q4/Q4 element for topology optimization, Struct. Multidiscip. Optim. 37 (2009) 255–264. [19] Z. Kang, Y.Q. Wang, A nodal variable method of structural topology optimization based on Shepard interpolant, Int. J. Numer. Methods Eng. 90 (2012) 329–342. [20] Z. Kang, Y.Q. Wang, Structural topology optimization based on non-local Shepard interpolation of density field, Comput. Methods Appl. Mech. Eng. 200 (2011) 3515–3525. [21] D. Shepard, A two-dimensional interpolation function for irregularly-spaced data, in: Proceedings of the 23rd National Conference, ACM, New York, 1968, pp. 517–523. [22] K.W. Brodlie, M.R. Asim, K. Unsworth, Constrained visualization using the Shepard interpolation family, Comput. Graph. Forum 24 (4) (2005) 809–820. [23] M.P. Bendsøe, O. Sigmund, Topology Optimization: Theory, Methods, and Applications, Springer, Berlin/Heidelberg/New York, 2003.