- Email: [email protected]
The topology optimization design for continuum structures based on the element free Galerkin method
ARTICLE IN PRESS Engineering Analysis with Boundary Elements 34 (2010) 666–672
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
The topology optimization design for continuum structures based on the element free Galerkin method Zheng Juan a,b,, Long Shuyao a,b, Li Guangyao a a b
State Key Laboratory of Advanced Design and Manufacture for Vehicle Body, Hunan University, Changsha, China College of Mechanics and Aerospace Engineering, Hunan University, Changsha, China
a r t i c l e in f o
a b s t r a c t
Article history: Received 17 June 2009 Accepted 13 February 2010 Available online 21 March 2010
In this paper, the element free Galerkin method (EFG), combined with evolutionary structural optimization method (ESO), is applied to carry out the topology optimization of the continuum structures. Considering the deletion criterion based on the stresses, the mathematical formulation of the topology optimization is developed. The objective function of this model is the minimized weight. Several numerical examples are used to prove the feasibility of the approach adopted in this paper. And the examples show the simplicity and fast convergence of the proposed method. & 2010 Elsevier Ltd. All rights reserved.
Keywords: Element free Galerkin method (EFG) Evolutionary structural optimization method (ESO) Topology optimization design of continuum structures Deletion criterion
1. Introduction The topology optimization design of the continuum structures is one of the most challenging research topics in the field of the structural optimization [1]. The purpose of the topology optimization design is to find the optimal lay-out of a structure within a specified region. In this problem the only known quantities are the applied loads, the possible support conditions, the volume of the structure to be constructed and possibly some additional design restrictions, and the physical size and the shape of the structure are unknown. The topology optimization design of the continuum structures is essentially a discretized 0–1 variables problem. For the topology optimization design of the continuum structures, homogenization approach [2], variable density approach [3] and evolutionary structural optimization (ESO) approach [4] are often employed. ESO method is a simple and robust technique for optimization problems applicable to various types of structures [5–7]. The basis idea of ESO method is that by gradually removing the inefficient material from the design domain, the residual structure evolves towards an optimum. The significance of ESO method is its simplicity and generality. To date, the prevailing numerical method in topology optimization is the finite element method (FEM). For FEM, frequent remeshing is unavoidable when dealing with large deformation or
Corresponding author at: College of Mechanics and Aerospace Engineering, Hunan University, Changsha, China. Tel.: + 86 0731 8824724. E-mail address: [email protected] (Z. Juan).
0955-7997/$ - see front matter & 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2010.03.001
moving boundary problems. To ameliorate these difficulties, a group of meshfree methods have been developed and achieved remarkable progress, such as smooth particle hydrodynamics method (SPH) [8], element free Galerkin method (EFG) [9], reproducing kernel particle method (RKPM) [10] and meshless local Petrov–Galerkin method (MLPG) [11]. These are methods in which the approximate solution is constructed entirely in terms of a set of nodes, and no element or characterization of the interrelationship of the nodes is needed. EFG method is one of the popularly used meshless methods. The EFG method has been successfully applied to large variety of problems including twodimension and three-dimension linear and nonlinear elastic problems [12], fracture and crack growth problems [13], plate and shell structures [14], and so on. In EFG, the shape function is constructed by moving least square (MLS), and control equation is produced from the weak form of variational equation. Because EFG shape function constructed using the MLS approximation lacks the Kronecker delta function property, special techniques are needed in the implementation of essential boundary conditions, for example, penalty method, Lagrange multiplier method, and so on. EFG method possesses an excellent rate of convergence and reasonably accurate results. Currently, there are few research works exploring structural topology optimization using meshless methods. Cho and Kwak [15] formulated a structural topology optimization method using the meshless method for geometric non-linearly modeling. Zhou and Zou [16] introduced the implicit topology description function method into meshless method for optimization of continua. In this paper, the element free Galerkin method (EFG), combined with evolutionary structural optimization method
ARTICLE IN PRESS Z. Juan et al. / Engineering Analysis with Boundary Elements 34 (2010) 666–672
(ESO), is applied to carry out the topology optimization of the continuum structures. Considering the deletion criterion based on the stresses, the mathematical formulation of the topology optimization is developed. The objective function of this model is the minimized weight. Finally, the feasibility and efficiency of the proposed method are illustrated with several 2D examples that are widely used in the topology optimization design. And the examples show the simplicity and fast convergence of the proposed method.
667
the matrix B in Eq. (3) is defined as BðxÞ ¼ PT W ¼ ½wðx1 ,xÞPðx1 Þ,wðx2 ,xÞPðx2 Þ, . . . ,wðxn ,xÞPðxn Þ Solving Eq. (3) for a(x), we have aðxÞ ¼ A1 ðxÞBðxÞUs Substituting the above equation into Eq. (1), we have uh ðxÞ ¼
n X
fi ðxÞu^ i ¼ UT ðxÞUs
ð4Þ
i¼1
2. Element free Galerkin method for plane elasticity [9]
where U(x) is the vector of MLS shape functions corresponding to n nodes in the support domain of the point x, and can be written as
2.1. Moving least square approximation
UT ðxÞ ¼ ½f1 ðxÞ,f2 ðxÞ, . . . ,fn ðxÞ ¼ PT ðxÞA1 ðxÞBðxÞ
The moving least square (MLS) method is generally considered to be one of the best schemes to interpolate data with reasonable accuracy. If field variable is a function such as u(x), the interpolation function uh(x) in a sub-domain Os can be defined over a number of scattered local points xi (i¼1, 2, y, n) by uh ðxÞ ¼ PT ðxÞaðxÞ,
8x A Os
ð1Þ
where P(x) is the basis function of the spatial coordinates. In some special problems, enhancement functions can, however, be added to the basis to improve the performance of the MLS approximation. We use only pure polynomial bases in this paper. For twodimensional problems, the complete polynomial basis functions are chosen as linear basis function
PT ðxÞ ¼ ½1,x1 ,x2 ,
quadratic basis function
m¼3
PT ðxÞ ¼ ½1,x1 ,x2 ,x21 ,x1 x2 ,x22 ,
m¼6
In Eq. (1), a(x) is a vector containing coefficients which are functions of the global Cartersian coordinates [x1, x2]T, depending on the polynomial basis. These coefficients can be obtained by minimizing a weighted discrete L2 norm defined as JðxÞ ¼
n X
wðxi ,xÞ½PT ðxi ÞaðxÞu^ i 2
ð2Þ
i¼1
where n is the number of nodes in the support domain of x for which the weight function w(xi, x) a0; ui is the nodal parameter of u at x¼xi; and w(xi, x) is the weight function associated with the node i. In this paper, weight function is cubic spline weight function as 8 2 3 ri r0:5 > < 2=34ri þ4ri , 2 wðxi ,xÞ ¼ 4=34ri þ 4ri 4=3ri3 , 0:5 ori r 1 > : 0, ri 41 where ri ¼
di jxxi j ¼ rw rw
in which di ¼jx xij is the distance from node xi to the interest point x, and rw is the size of the support domain for the weight function. The stationary of J with respect to a(x) leads to the following set of linear relation: AðxÞaðxÞ ¼ BðxÞUs
ð3Þ
where Us is the vector that collects the nodal parameters of the field function for all the nodes in the support domain, Us ¼ ½u^ 1 ,u^ 2 , . . . ,u^ n T ; and A(x) is called the weighted moment matrix defined by AðxÞ ¼ PT WP ¼
n X i¼1
wðxi ,xÞPðxi ÞPT ðxi Þ
ð5Þ
2.2. Discrete equations in 2D problem Considering the following standard two-dimension problem of linear elasticity defined in the domain O bounded by the boundary G 8 T > < L r þ b ¼ 0 in O ð6Þ rn ¼ t on Gt > : u¼u on Gu where L is the differential operator; r is the stress vector; u is the displacement vector; b is the body force vector; t is the prescribed traction on the natural boundaries; u is the prescribed displacement on the essential boundaries; n is the vector of unit outward at a point on the natural boundary. The element free Galerkin method uses the moving least squares (MLS) shape functions. Because the MLS approximation lacks the Kronecker delta function property, the constrained Galerkin weak form should be posed as follows: Z Z Z dðLuÞT DðLuÞ dO duT b dO duT t dG O O Gt Z 1 d ð7Þ ðuuÞT aðuuÞ dG ¼ 0 Gu 2 where a ¼ ja1 , a2 , . . . , ak j is a diagonal matrix of penalty factors. Using the MLS shape functions constructed using n nodes in the local support domain, we have 8 9 u^ 1 > > > > > > > > #> ^ > > " < v1 > = f 0 . . . f 0 u 1 n ^ ¼ ð8Þ uh ¼ ¼ Uu 0 f1 . . . 0 fn > > v > > > > ^ u > > n > > > > :^ ; vn where, U is the matrix of the shape functions, u is the vector of the displacements at the field nodes in the support domain, and n is the number of nodes in the support domain of a sampling point at x. Using strain–displacement equations, we can have 3 2 8 9 @ u1 > > > 6 @x 0 7 > > > > > 7" > 6 #> > = 7 f1 0 fn 0 < v1 > 6 @ 7 60 h ^ e ¼ Lu ¼ LUu ¼ 6 7 @y 7 0 f1 0 fn > > 6 > > un > > 6 @ > > @ 7 > > 5 4 > : > ; vn @y @x 38 9 2 @f1 @fn u1 > 0 ... 0 7> > 6 @x > > > @x > > 7> > 6 > n < v1 > = 6 X @f1 @fn 7 7 6 0 ... 0 ^ BI uI ð9Þ ¼ Bu ¼ ¼6 7 @y @y > 7> 6 > I > > 7> 6 u > > n 4 @f1 @f1 @fn @fn 5> > > > : ; ... vn @y @x @y @x
ARTICLE IN PRESS 668
Z. Juan et al. / Engineering Analysis with Boundary Elements 34 (2010) 666–672
Substituting Eq. (8) into the second term of Eq. (7), and using the same arguments in deriving the stiffness matrix, we have
Similarly, Lduh ¼ LUdu ¼ Bdu ¼
n X
BI duI
Z
I
The stress vector can be obtained using the constitutive equations:
r ¼ De ¼ BDu ¼
n X
DBI duI
I
Substituting Eq. (9) into the first term of Eq. (7), we have Z X n n X ðLduÞT ðDLuÞ dO ¼ ð BI duI ÞT ð DBJ uJ ÞdO O
O
Z
I
O
I
duTI ½BTI DBJ uJ dO
ð10Þ
J
O
O
I
J
Move the integration inside the summations to arrive at Z N X N X ðL duÞT ðDLuÞ dO ¼ duTI ð ½BTI DBJ dOÞuJ
O
I
¼
J
N X N X I
I
N X
duTI FbI ¼ dUT Fb
ð13Þ
I
b
F is the global body force vector assembled using the nodal body force vectors for all the nodes in the entire problem domain, and is defined as 8 9 b > > > < F1 > = ^ Fb ¼ > > > : Fb > ; The treatment for the third term in Eq. (7) is exactly the same as that for the second term, except that the body force vector is replaced by the traction vector and the integrations are replaced by the boundary integrations. Hence, we can obtain
duT t dG ¼ Gt
ðLduÞ ðDLuÞ dO ¼ dU KU
ð14Þ
Gt
t
F is the global traction force vector assembled using the nodal traction force vectors which is defined as 8 t 9 > < F1 > = t ^ F ¼ > : Ft > ;
Gu
ð11Þ
N X N N X X 1 duTI KaIJ uJ duTI FaI ðuuÞT aðuuÞ dG ¼ 2 I J I
¼ dUT Ka UdUT Fa
J
T
duTI FtI ¼ dUT Ft
I
where FtI is the nodal traction force vector that is defined as Z UTI t dG FtI ¼
d duTI KIJ uJ
N X
Substituting Eq. (8) into the last term of Eq. (7), we have
O
ð12Þ
the global stiffness matrix in the form of 3 K12 . . . K1N 7 K22 . . . K2N 7 7 ^ & ^ 7 5 KN2 . . . KNN
In Eq. (12), the vector U is the global displacement vector that collects the nodal displacements of all the nodes in the entire problem domain, which has the form of 8 9 8 9 > u1 > > u1 > > > > > > > > > > > > < u2 > = > < v1 > = ^ U¼ ¼ > ^ > > > > > > > > > uN > > : ; > > > > > > uN : v ;
ð15Þ
where KaIJ is the nodal penalty stiffness matrix and FaI is the nodal penalty force vector, defined as Z Z KaIJ ¼ UTI aUJ dG, FaI ¼ UTI au dG Gu
O
N
UI uI ÞT b dO ¼
where FbI is the nodal body force vector that is defined as Z UTI b dO FbI ¼
Z
O
Finally, Eq. (11) becomes
where K is 2 K11 6 6 K21 K¼6 6 ^ 4 KN1
O
n X
N
where KIJ is called the nodal stiffness matrix and is defined as Z KIJ ¼ BTI DBJ dO
T
O
Z
Note that until this stage, I and J are based on the local numbering system for the nodes in the local support domain. We can now change the numbering system from the local one to the global one that records all the field nodes in the entire domain in a unique manner from 1 to N, the total numbering of nodes in the problem domain. Therefore, both I and J in Eq. (10) can now vary from 1 to N. When node I and node J are not in the same local support domain, the integrand vanishes. With this operation, Eq. (10) can be expressed as Z X Z N X N ðLduÞT ðDLuÞ dO ¼ duTI ½BTI DBJ uJ dO
Z
ðd
N
J
n X n X
¼
Z
Z
O
where BI is the strain matrix about the jth node and D is the material matrix for the plane stress problem. They are given by 2 3 2 3 fI,x 0 1 m 0 E 6 7 6 0 7 fI,y 5, D ¼ 0 BI ¼ 4 4m 1 5 1m2 fI,y fI,x 0 0 ð1mÞ=2 Z
duT b dO ¼
Gu
Ka is the global penalty stiffness matrix assembled using the nodal penalty stiffness matrices; Fa is the global penalty external force vector assembled using the nodal penalty force vectors. Substituting Eqs. (12)–(15) into Eq. (7), we can obtain
dUT KUdUT Fb dUT Ft þ dUT Ka UdUT Fa ¼ 0 or
dUT ðKU þKa UFb Ft Fa Þ ¼ 0 Because dU is arbitrary, the above equation can be satisfied only if ½K þ Ka U ¼ F þ Fa where F is the global force vector given by F ¼ Fb þ Ft Eq. (16) is the final discretized system equation.
ð16Þ
ARTICLE IN PRESS Z. Juan et al. / Engineering Analysis with Boundary Elements 34 (2010) 666–672
669
2.3. Enforcement of essential boundary conditions
3.3. Topology optimization procedure
In the present EFG formulation, the penalty method is used to enforce essential boundary conditions. The advantage of using the penalty method is that the dimension, symmetry and positive definite properties of the stiffness matrix are achieved, as long as the penalty factors chosen are positive. In addition, the symmetry and the bandness of the system matrix are preserved. Several other strategies have also been developed to enforce essential boundary conditions, such as, the Lagrange multiplier method [9], the method using the modified variational principle [17], the method coupling with the finite method [18], the orthogonal transform technique [19], and so on.
According to the theory above, an outline of the ESO algorithm based on EFG can be summarized as follows:
3. Formulation of the topology optimization design 3.1. ESO based on stress criteria 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. Stress is an important indication of material usage of a structure. In the plane stress problem, the Von Mises stresses is defined by qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð17Þ sVM ¼ s2x þ s2y sx sy þ 3t2xy The maximum Von Mises stress is used as the node removal criteria, which is expressed by VM sVM e =smax o RRi
Step 1: Define the design domain, the boundary conditions. Discretize the structure by using EFG method. Step 2: Analyse the structure, compute each node’s Von Mises stresses. Step 3: The nodes which satisfy the inequality (18) are removed. Step 4: The cycle of the meshless analysis and node removals is repeated using the same value of RRi until no more nodes can satisfy the inequality (18). At the steady state, no more nodes can be deleted from the structure. Step 5: Introduce evolutionary ratio ER, increase the rejection ratio by using Eq. (19). Step 6: Repeat steps (2)–(5), until a prescribed percentage of volume is reached.
4. Numerical examples 4.1. Example 1 A cantilever beam with size 5 12 1 m3, as shown in Fig. 1, is now discussed. The beam is fixed on the left side and loaded with a concentrated force F ¼15.6 kN at the middle of the right side. The design domain is discretized by 592 nodes uniformly
ð18Þ
where sVM is the Von Mises stress of the eth node; sVM e max is the maximum Von Mises stress in the current structure; RRi is the rejection ratio at the ith steady state used to control the node removal process. All nodes that satisfy inequality (18) are deleted from the structure. The cycle of node removals and meshless analysis is repeated until no more nodes can be deleted from the structure at the current steady state. At this stage, an evolutionary ratio ER is introduced and added to the rejection ratio, which becomes RRi þ 1 ¼ RRi þER
12 F
ð19Þ
With this increased rejection ratio, the evolutionary optimization process continues until a desired optimum is obtained.
3.2. Topology optimization model
5 Fig. 1. Cantilever beam model.
Considering the minimized weight as the objective function, the mathematical formulation of the topology optimization based on the stress criterion can be formulated as follows: 8 n X > > > wi xi > min W ¼ > > > i¼1 > > < sVM i ð20Þ s:t: 4 RRi VM > s > max > > > > V ¼ fV0 > > > : xi A f0,1g i ¼ 1, . . . ,n where W is the total weight of a structure; wi is the weight of the ith node; xi is the design variable which takes the value 0 or 1, xi ¼0 represents remove the node, and xi ¼1 represents retain the node; n is the total number of nodes; V is the material volume of the design domain; V0 is the given volume of the design domain; f is the prescribed volume fraction.
Fig. 2. Optimization result obtained by the present method.
ARTICLE IN PRESS 670
Z. Juan et al. / Engineering Analysis with Boundary Elements 34 (2010) 666–672
Fig. 3. Optimization results obtained at different iteration.
Weight of the beam
60 50 40 30 20 10 0
2
4
6
8
10
Iteration number Fig. 4. Convergence history of the cantilever beam using the present method.
Fig. 5. (a) Optimization result by FEM without sensitivity filtering and (b) optimization result by FEM with sensitivity filtering.
30
15
F Fig. 6. Simply supported beam model.
distributed, 540 integration cells would be used according to the location of the nodes, and 2 2 Gauss quadrature is used. The elastic material properties are chosen as Young’s modulus
Fig. 7. Optimization result obtained by the present method.
E ¼ 68:89 108 N=m2 , Poisson’s ratio m ¼0.3. The initial total weight of the structure is 60t. The initial rejection ratio is RR0 ¼0.02, the evolutionary ratio is ER ¼0.01, and the volume constraint is limited to 20%. The optimization result of the beam obtained by the present method is shown in Fig. 2, the number of iterations is 9. The optimization results obtained at different iteration are presented in Fig. 3. Fig. 4 gives the curve of convergence of the cantilever beam’s objective function. The almost monotonic and uniform convergence can be observed from this figure. In Ref. [20], considering the minimization of compliance as objective function, the same model is solved by the finite element
ARTICLE IN PRESS Z. Juan et al. / Engineering Analysis with Boundary Elements 34 (2010) 666–672
671
method (FEM). The optimization result by FEM without sensitivity filtering is shown in Fig. 5(a), the number of iterations is 14, and the optimization result by FEM with sensitivity filtering is shown in Fig. 5(b), the number of iterations is 23. From these results, it can be seen that similar optimization results were obtained by two different algorithms, and a smoother optimization result can be obtained by the present method, and the number of iterations of the present approach is less than the FEM. It also can be seen that the present approach can effectively avoid the checkerboard pattern.
4.2. Example 2 A simply supported beam with size 30 15 1 m3, as shown in Fig. 6, is now discussed. The two corner points on the bottom side are fixed, and a concentrated force F¼1 kN is loaded at the center of the bottom surface. The elastic material properties are chosen as Young’s modulus E¼30 Gpa, Poisson’s ratio m ¼0.3. The initial total weight of the structure is 450t. The initial rejection ratio is RR0 ¼0.01, the evolutionary ratio is ER¼0.01, and the volume constraint is limited to 30%. To determine an optimum structural layout, the design domain is discretized by EFG method using 31 16, 37 19, 43 22, 49 25 uniformly distributed nodes. And 30 15, 36 18, 42 21, 48 24 integration cells are used according to the
Fig. 9. Optimization result obtained by FEM with sensitivity filtering.
500
31 × 37 × 43 × 49 ×
450
Weight of the beam
400
16 nodes 19 nodes 22 nodes 25 nodes
350 300 250 200 150 100 0
5
10
15
20
25
30
35
Iteration number Fig. 8. Optimization result obtained by FEM without sensitivity filtering.
Fig. 10. Convergence history of the simply supported beam based on different nodes using the present method.
ARTICLE IN PRESS 672
Z. Juan et al. / Engineering Analysis with Boundary Elements 34 (2010) 666–672
location of the nodes and 2 2 Gauss quadrature is used. In Ref. [20], considering the minimization of compliance as objective function, and the volume of the structure as constraint, the same model is solved by the finite element method (FEM). The optimization results obtained by the present method is shown in Fig. 7, the optimization results obtained by FEM without sensitivity filtering is shown in Fig. 8, and the optimization results obtained by FEM with sensitivity filtering is shown in Fig. 9. From these results, it can be seen that similar optimization results were obtained by two different algorithms, and a smoother optimization result can be obtained by the present method, and the present approach can effectively eliminate the checkerboard pattern. Fig. 10 gives the curves of convergence of the simply supported beam’s objective function based on different nodes using the present method. The almost monotonic and uniform convergence can be observed from this figure. And for different nodes, their convergence characteristics are very similar.
5. Conclusions In this paper, the element free Galerkin method (EFG), combined with evolutionary structural optimization method (ESO), is applied to solve the topology optimization problem of the continuum structures. Considering the minimized weight as the objective function, the mathematical formulation of the topology optimization is developed based on the stress criterion. Several numerical examples are shown to prove the validity and feasibility of the present method. And the examples also show the simplicity and fast convergence of the proposed method.
Acknowledgements This work was supported by National 973 Scientific and Technological Innovation Project (2004CB719402), Natural Science Foundation of China (No. 10672055), The Key Project of NSFC of China (No. 60635020) and The National Science Foundation for Outstanding Youth of China (No. 50625519).
References [1] Bendsoe MP, Sigmund O. Topology optimization: theory, methods, and applications. Berlin, Heidelberg, New York: Springer; 2003. [2] Bendsoe MP, Kikuchi N. Generating optimal topology in structural design using a homogenization method. Computer Methods in Applied Mechanics and Engineering 1989;71:197–224. [3] Bendsoe MP, Sigmund O. Material interpolation schemes in topology optimization. Archive of Applied Mechanics 1999;69:635–54. [4] Xie YM, Steven GP. A simple evolutionary procedure for structural optimization. Computer and Structures 1993;49:885–96. [5] Xie YM, Stevel GP. Evolutionary structural optimization for dynamic problems. Computer and Structures 1996;58:1067–73. [6] Manickarajah D, Xie YM, Steven GP. An evolutionary method for optimization of plate buckling resistance. Finite Elements in Analysis and Design 1998;29:203–30. [7] Rong JH, Xie YM, Yang XY. Topology optimization of structures under dynamic response constraints. Journal Sound and Vibration 2000;234(2):l77–189. [8] Monaghan JJ. Smoothed particle hydrodynamics. Annual Review of Astronomy and Astrophysics 1992;30:543–74. [9] Belytschko T, Lu YY, Gu L. Element-free Galerkin method. International Journal for Numerical Methods in Engineering 1994;37:229–56. [10] Liu WK, Jun S, Li S. Reproducing kernel particle methods for structural dynamics. International Journal for Numerical Methods in Engineering 1995;38(10):1655–79. [11] Atluri SN, Zhu T. A new meshless local Petrov–Galerkin (MLPG) approach in computational mechanics. Computational Mechanics 1998;22:117–27. [12] Belytschko T, Krysl P, Krongauz Y. A three-dimensional explicit element-free Galerkin method. International Journal for Numerical Methods in Fluids 1997;24(12):1253–70. [13] Belytschko T, Gu L, Lu YY. Fracture and crack growth by element free Galerkin methods. Modelling and Simulation in Materials Science and Engineering 1994;2:519–34. [14] Liu GR, Chen XL. A mesh-free method for static and free vibration analyses of thin plates of complicated shape. Journal of sound and vibration 2001;241(5):839–55. [15] Cho S, Kwak J. Topology design optimization of geometrically non-linear structures using meshfree method. Computer Methods in Applied Mechanics and Engineering 2006;195:5909–25. [16] Zhou JX, Zou W. Meshless approximation combined with implicit topology description for optimization of continua. Structural and Multidisciplinary Optimization 2008;36(4):347–53. [17] Lu YY, Belytschko T, Gu L. A new implementation of the element free Galerkin method. Computer Methods in Applied Mechanics and Engineering 1994;113:397–414. [18] Krongauz Y, Belytschko T. Enforcement of essential boundary conditions in meshless approximations using finite elements. Computer Methods in Applied Mechanics and Engineering 1996;131:133–45. [19] Atluri SN, Kim HG, Cho JY. A critical assessment of the truly meshless local Petrov–Galerkin (MLPG), and local boundary integral equation (LBIE) methods. Computational Mechanics 1999;24(5):348–72. [20] Sigmund O. A 99 line topology optimization code written in Matlab. Structural and Multidiscipline Optimization 2001;21(2):120–7.
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
The topology optimization design for continuum structures based on the element free Galerkin method Zheng Juan a,b,, Long Shuyao a,b, Li Guangyao a a b
State Key Laboratory of Advanced Design and Manufacture for Vehicle Body, Hunan University, Changsha, China College of Mechanics and Aerospace Engineering, Hunan University, Changsha, China
a r t i c l e in f o
a b s t r a c t
Article history: Received 17 June 2009 Accepted 13 February 2010 Available online 21 March 2010
In this paper, the element free Galerkin method (EFG), combined with evolutionary structural optimization method (ESO), is applied to carry out the topology optimization of the continuum structures. Considering the deletion criterion based on the stresses, the mathematical formulation of the topology optimization is developed. The objective function of this model is the minimized weight. Several numerical examples are used to prove the feasibility of the approach adopted in this paper. And the examples show the simplicity and fast convergence of the proposed method. & 2010 Elsevier Ltd. All rights reserved.
Keywords: Element free Galerkin method (EFG) Evolutionary structural optimization method (ESO) Topology optimization design of continuum structures Deletion criterion
1. Introduction The topology optimization design of the continuum structures is one of the most challenging research topics in the field of the structural optimization [1]. The purpose of the topology optimization design is to find the optimal lay-out of a structure within a specified region. In this problem the only known quantities are the applied loads, the possible support conditions, the volume of the structure to be constructed and possibly some additional design restrictions, and the physical size and the shape of the structure are unknown. The topology optimization design of the continuum structures is essentially a discretized 0–1 variables problem. For the topology optimization design of the continuum structures, homogenization approach [2], variable density approach [3] and evolutionary structural optimization (ESO) approach [4] are often employed. ESO method is a simple and robust technique for optimization problems applicable to various types of structures [5–7]. The basis idea of ESO method is that by gradually removing the inefficient material from the design domain, the residual structure evolves towards an optimum. The significance of ESO method is its simplicity and generality. To date, the prevailing numerical method in topology optimization is the finite element method (FEM). For FEM, frequent remeshing is unavoidable when dealing with large deformation or
Corresponding author at: College of Mechanics and Aerospace Engineering, Hunan University, Changsha, China. Tel.: + 86 0731 8824724. E-mail address: [email protected] (Z. Juan).
0955-7997/$ - see front matter & 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2010.03.001
moving boundary problems. To ameliorate these difficulties, a group of meshfree methods have been developed and achieved remarkable progress, such as smooth particle hydrodynamics method (SPH) [8], element free Galerkin method (EFG) [9], reproducing kernel particle method (RKPM) [10] and meshless local Petrov–Galerkin method (MLPG) [11]. These are methods in which the approximate solution is constructed entirely in terms of a set of nodes, and no element or characterization of the interrelationship of the nodes is needed. EFG method is one of the popularly used meshless methods. The EFG method has been successfully applied to large variety of problems including twodimension and three-dimension linear and nonlinear elastic problems [12], fracture and crack growth problems [13], plate and shell structures [14], and so on. In EFG, the shape function is constructed by moving least square (MLS), and control equation is produced from the weak form of variational equation. Because EFG shape function constructed using the MLS approximation lacks the Kronecker delta function property, special techniques are needed in the implementation of essential boundary conditions, for example, penalty method, Lagrange multiplier method, and so on. EFG method possesses an excellent rate of convergence and reasonably accurate results. Currently, there are few research works exploring structural topology optimization using meshless methods. Cho and Kwak [15] formulated a structural topology optimization method using the meshless method for geometric non-linearly modeling. Zhou and Zou [16] introduced the implicit topology description function method into meshless method for optimization of continua. In this paper, the element free Galerkin method (EFG), combined with evolutionary structural optimization method
ARTICLE IN PRESS Z. Juan et al. / Engineering Analysis with Boundary Elements 34 (2010) 666–672
(ESO), is applied to carry out the topology optimization of the continuum structures. Considering the deletion criterion based on the stresses, the mathematical formulation of the topology optimization is developed. The objective function of this model is the minimized weight. Finally, the feasibility and efficiency of the proposed method are illustrated with several 2D examples that are widely used in the topology optimization design. And the examples show the simplicity and fast convergence of the proposed method.
667
the matrix B in Eq. (3) is defined as BðxÞ ¼ PT W ¼ ½wðx1 ,xÞPðx1 Þ,wðx2 ,xÞPðx2 Þ, . . . ,wðxn ,xÞPðxn Þ Solving Eq. (3) for a(x), we have aðxÞ ¼ A1 ðxÞBðxÞUs Substituting the above equation into Eq. (1), we have uh ðxÞ ¼
n X
fi ðxÞu^ i ¼ UT ðxÞUs
ð4Þ
i¼1
2. Element free Galerkin method for plane elasticity [9]
where U(x) is the vector of MLS shape functions corresponding to n nodes in the support domain of the point x, and can be written as
2.1. Moving least square approximation
UT ðxÞ ¼ ½f1 ðxÞ,f2 ðxÞ, . . . ,fn ðxÞ ¼ PT ðxÞA1 ðxÞBðxÞ
The moving least square (MLS) method is generally considered to be one of the best schemes to interpolate data with reasonable accuracy. If field variable is a function such as u(x), the interpolation function uh(x) in a sub-domain Os can be defined over a number of scattered local points xi (i¼1, 2, y, n) by uh ðxÞ ¼ PT ðxÞaðxÞ,
8x A Os
ð1Þ
where P(x) is the basis function of the spatial coordinates. In some special problems, enhancement functions can, however, be added to the basis to improve the performance of the MLS approximation. We use only pure polynomial bases in this paper. For twodimensional problems, the complete polynomial basis functions are chosen as linear basis function
PT ðxÞ ¼ ½1,x1 ,x2 ,
quadratic basis function
m¼3
PT ðxÞ ¼ ½1,x1 ,x2 ,x21 ,x1 x2 ,x22 ,
m¼6
In Eq. (1), a(x) is a vector containing coefficients which are functions of the global Cartersian coordinates [x1, x2]T, depending on the polynomial basis. These coefficients can be obtained by minimizing a weighted discrete L2 norm defined as JðxÞ ¼
n X
wðxi ,xÞ½PT ðxi ÞaðxÞu^ i 2
ð2Þ
i¼1
where n is the number of nodes in the support domain of x for which the weight function w(xi, x) a0; ui is the nodal parameter of u at x¼xi; and w(xi, x) is the weight function associated with the node i. In this paper, weight function is cubic spline weight function as 8 2 3 ri r0:5 > < 2=34ri þ4ri , 2 wðxi ,xÞ ¼ 4=34ri þ 4ri 4=3ri3 , 0:5 ori r 1 > : 0, ri 41 where ri ¼
di jxxi j ¼ rw rw
in which di ¼jx xij is the distance from node xi to the interest point x, and rw is the size of the support domain for the weight function. The stationary of J with respect to a(x) leads to the following set of linear relation: AðxÞaðxÞ ¼ BðxÞUs
ð3Þ
where Us is the vector that collects the nodal parameters of the field function for all the nodes in the support domain, Us ¼ ½u^ 1 ,u^ 2 , . . . ,u^ n T ; and A(x) is called the weighted moment matrix defined by AðxÞ ¼ PT WP ¼
n X i¼1
wðxi ,xÞPðxi ÞPT ðxi Þ
ð5Þ
2.2. Discrete equations in 2D problem Considering the following standard two-dimension problem of linear elasticity defined in the domain O bounded by the boundary G 8 T > < L r þ b ¼ 0 in O ð6Þ rn ¼ t on Gt > : u¼u on Gu where L is the differential operator; r is the stress vector; u is the displacement vector; b is the body force vector; t is the prescribed traction on the natural boundaries; u is the prescribed displacement on the essential boundaries; n is the vector of unit outward at a point on the natural boundary. The element free Galerkin method uses the moving least squares (MLS) shape functions. Because the MLS approximation lacks the Kronecker delta function property, the constrained Galerkin weak form should be posed as follows: Z Z Z dðLuÞT DðLuÞ dO duT b dO duT t dG O O Gt Z 1 d ð7Þ ðuuÞT aðuuÞ dG ¼ 0 Gu 2 where a ¼ ja1 , a2 , . . . , ak j is a diagonal matrix of penalty factors. Using the MLS shape functions constructed using n nodes in the local support domain, we have 8 9 u^ 1 > > > > > > > > #> ^ > > " < v1 > = f 0 . . . f 0 u 1 n ^ ¼ ð8Þ uh ¼ ¼ Uu 0 f1 . . . 0 fn > > v > > > > ^ u > > n > > > > :^ ; vn where, U is the matrix of the shape functions, u is the vector of the displacements at the field nodes in the support domain, and n is the number of nodes in the support domain of a sampling point at x. Using strain–displacement equations, we can have 3 2 8 9 @ u1 > > > 6 @x 0 7 > > > > > 7" > 6 #> > = 7 f1 0 fn 0 < v1 > 6 @ 7 60 h ^ e ¼ Lu ¼ LUu ¼ 6 7 @y 7 0 f1 0 fn > > 6 > > un > > 6 @ > > @ 7 > > 5 4 > : > ; vn @y @x 38 9 2 @f1 @fn u1 > 0 ... 0 7> > 6 @x > > > @x > > 7> > 6 > n < v1 > = 6 X @f1 @fn 7 7 6 0 ... 0 ^ BI uI ð9Þ ¼ Bu ¼ ¼6 7 @y @y > 7> 6 > I > > 7> 6 u > > n 4 @f1 @f1 @fn @fn 5> > > > : ; ... vn @y @x @y @x
ARTICLE IN PRESS 668
Z. Juan et al. / Engineering Analysis with Boundary Elements 34 (2010) 666–672
Substituting Eq. (8) into the second term of Eq. (7), and using the same arguments in deriving the stiffness matrix, we have
Similarly, Lduh ¼ LUdu ¼ Bdu ¼
n X
BI duI
Z
I
The stress vector can be obtained using the constitutive equations:
r ¼ De ¼ BDu ¼
n X
DBI duI
I
Substituting Eq. (9) into the first term of Eq. (7), we have Z X n n X ðLduÞT ðDLuÞ dO ¼ ð BI duI ÞT ð DBJ uJ ÞdO O
O
Z
I
O
I
duTI ½BTI DBJ uJ dO
ð10Þ
J
O
O
I
J
Move the integration inside the summations to arrive at Z N X N X ðL duÞT ðDLuÞ dO ¼ duTI ð ½BTI DBJ dOÞuJ
O
I
¼
J
N X N X I
I
N X
duTI FbI ¼ dUT Fb
ð13Þ
I
b
F is the global body force vector assembled using the nodal body force vectors for all the nodes in the entire problem domain, and is defined as 8 9 b > > > < F1 > = ^ Fb ¼ > > > : Fb > ; The treatment for the third term in Eq. (7) is exactly the same as that for the second term, except that the body force vector is replaced by the traction vector and the integrations are replaced by the boundary integrations. Hence, we can obtain
duT t dG ¼ Gt
ðLduÞ ðDLuÞ dO ¼ dU KU
ð14Þ
Gt
t
F is the global traction force vector assembled using the nodal traction force vectors which is defined as 8 t 9 > < F1 > = t ^ F ¼ > : Ft > ;
Gu
ð11Þ
N X N N X X 1 duTI KaIJ uJ duTI FaI ðuuÞT aðuuÞ dG ¼ 2 I J I
¼ dUT Ka UdUT Fa
J
T
duTI FtI ¼ dUT Ft
I
where FtI is the nodal traction force vector that is defined as Z UTI t dG FtI ¼
d duTI KIJ uJ
N X
Substituting Eq. (8) into the last term of Eq. (7), we have
O
ð12Þ
the global stiffness matrix in the form of 3 K12 . . . K1N 7 K22 . . . K2N 7 7 ^ & ^ 7 5 KN2 . . . KNN
In Eq. (12), the vector U is the global displacement vector that collects the nodal displacements of all the nodes in the entire problem domain, which has the form of 8 9 8 9 > u1 > > u1 > > > > > > > > > > > > < u2 > = > < v1 > = ^ U¼ ¼ > ^ > > > > > > > > > uN > > : ; > > > > > > uN : v ;
ð15Þ
where KaIJ is the nodal penalty stiffness matrix and FaI is the nodal penalty force vector, defined as Z Z KaIJ ¼ UTI aUJ dG, FaI ¼ UTI au dG Gu
O
N
UI uI ÞT b dO ¼
where FbI is the nodal body force vector that is defined as Z UTI b dO FbI ¼
Z
O
Finally, Eq. (11) becomes
where K is 2 K11 6 6 K21 K¼6 6 ^ 4 KN1
O
n X
N
where KIJ is called the nodal stiffness matrix and is defined as Z KIJ ¼ BTI DBJ dO
T
O
Z
Note that until this stage, I and J are based on the local numbering system for the nodes in the local support domain. We can now change the numbering system from the local one to the global one that records all the field nodes in the entire domain in a unique manner from 1 to N, the total numbering of nodes in the problem domain. Therefore, both I and J in Eq. (10) can now vary from 1 to N. When node I and node J are not in the same local support domain, the integrand vanishes. With this operation, Eq. (10) can be expressed as Z X Z N X N ðLduÞT ðDLuÞ dO ¼ duTI ½BTI DBJ uJ dO
Z
ðd
N
J
n X n X
¼
Z
Z
O
where BI is the strain matrix about the jth node and D is the material matrix for the plane stress problem. They are given by 2 3 2 3 fI,x 0 1 m 0 E 6 7 6 0 7 fI,y 5, D ¼ 0 BI ¼ 4 4m 1 5 1m2 fI,y fI,x 0 0 ð1mÞ=2 Z
duT b dO ¼
Gu
Ka is the global penalty stiffness matrix assembled using the nodal penalty stiffness matrices; Fa is the global penalty external force vector assembled using the nodal penalty force vectors. Substituting Eqs. (12)–(15) into Eq. (7), we can obtain
dUT KUdUT Fb dUT Ft þ dUT Ka UdUT Fa ¼ 0 or
dUT ðKU þKa UFb Ft Fa Þ ¼ 0 Because dU is arbitrary, the above equation can be satisfied only if ½K þ Ka U ¼ F þ Fa where F is the global force vector given by F ¼ Fb þ Ft Eq. (16) is the final discretized system equation.
ð16Þ
ARTICLE IN PRESS Z. Juan et al. / Engineering Analysis with Boundary Elements 34 (2010) 666–672
669
2.3. Enforcement of essential boundary conditions
3.3. Topology optimization procedure
In the present EFG formulation, the penalty method is used to enforce essential boundary conditions. The advantage of using the penalty method is that the dimension, symmetry and positive definite properties of the stiffness matrix are achieved, as long as the penalty factors chosen are positive. In addition, the symmetry and the bandness of the system matrix are preserved. Several other strategies have also been developed to enforce essential boundary conditions, such as, the Lagrange multiplier method [9], the method using the modified variational principle [17], the method coupling with the finite method [18], the orthogonal transform technique [19], and so on.
According to the theory above, an outline of the ESO algorithm based on EFG can be summarized as follows:
3. Formulation of the topology optimization design 3.1. ESO based on stress criteria 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. Stress is an important indication of material usage of a structure. In the plane stress problem, the Von Mises stresses is defined by qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð17Þ sVM ¼ s2x þ s2y sx sy þ 3t2xy The maximum Von Mises stress is used as the node removal criteria, which is expressed by VM sVM e =smax o RRi
Step 1: Define the design domain, the boundary conditions. Discretize the structure by using EFG method. Step 2: Analyse the structure, compute each node’s Von Mises stresses. Step 3: The nodes which satisfy the inequality (18) are removed. Step 4: The cycle of the meshless analysis and node removals is repeated using the same value of RRi until no more nodes can satisfy the inequality (18). At the steady state, no more nodes can be deleted from the structure. Step 5: Introduce evolutionary ratio ER, increase the rejection ratio by using Eq. (19). Step 6: Repeat steps (2)–(5), until a prescribed percentage of volume is reached.
4. Numerical examples 4.1. Example 1 A cantilever beam with size 5 12 1 m3, as shown in Fig. 1, is now discussed. The beam is fixed on the left side and loaded with a concentrated force F ¼15.6 kN at the middle of the right side. The design domain is discretized by 592 nodes uniformly
ð18Þ
where sVM is the Von Mises stress of the eth node; sVM e max is the maximum Von Mises stress in the current structure; RRi is the rejection ratio at the ith steady state used to control the node removal process. All nodes that satisfy inequality (18) are deleted from the structure. The cycle of node removals and meshless analysis is repeated until no more nodes can be deleted from the structure at the current steady state. At this stage, an evolutionary ratio ER is introduced and added to the rejection ratio, which becomes RRi þ 1 ¼ RRi þER
12 F
ð19Þ
With this increased rejection ratio, the evolutionary optimization process continues until a desired optimum is obtained.
3.2. Topology optimization model
5 Fig. 1. Cantilever beam model.
Considering the minimized weight as the objective function, the mathematical formulation of the topology optimization based on the stress criterion can be formulated as follows: 8 n X > > > wi xi > min W ¼ > > > i¼1 > > < sVM i ð20Þ s:t: 4 RRi VM > s > max > > > > V ¼ fV0 > > > : xi A f0,1g i ¼ 1, . . . ,n where W is the total weight of a structure; wi is the weight of the ith node; xi is the design variable which takes the value 0 or 1, xi ¼0 represents remove the node, and xi ¼1 represents retain the node; n is the total number of nodes; V is the material volume of the design domain; V0 is the given volume of the design domain; f is the prescribed volume fraction.
Fig. 2. Optimization result obtained by the present method.
ARTICLE IN PRESS 670
Z. Juan et al. / Engineering Analysis with Boundary Elements 34 (2010) 666–672
Fig. 3. Optimization results obtained at different iteration.
Weight of the beam
60 50 40 30 20 10 0
2
4
6
8
10
Iteration number Fig. 4. Convergence history of the cantilever beam using the present method.
Fig. 5. (a) Optimization result by FEM without sensitivity filtering and (b) optimization result by FEM with sensitivity filtering.
30
15
F Fig. 6. Simply supported beam model.
distributed, 540 integration cells would be used according to the location of the nodes, and 2 2 Gauss quadrature is used. The elastic material properties are chosen as Young’s modulus
Fig. 7. Optimization result obtained by the present method.
E ¼ 68:89 108 N=m2 , Poisson’s ratio m ¼0.3. The initial total weight of the structure is 60t. The initial rejection ratio is RR0 ¼0.02, the evolutionary ratio is ER ¼0.01, and the volume constraint is limited to 20%. The optimization result of the beam obtained by the present method is shown in Fig. 2, the number of iterations is 9. The optimization results obtained at different iteration are presented in Fig. 3. Fig. 4 gives the curve of convergence of the cantilever beam’s objective function. The almost monotonic and uniform convergence can be observed from this figure. In Ref. [20], considering the minimization of compliance as objective function, the same model is solved by the finite element
ARTICLE IN PRESS Z. Juan et al. / Engineering Analysis with Boundary Elements 34 (2010) 666–672
671
method (FEM). The optimization result by FEM without sensitivity filtering is shown in Fig. 5(a), the number of iterations is 14, and the optimization result by FEM with sensitivity filtering is shown in Fig. 5(b), the number of iterations is 23. From these results, it can be seen that similar optimization results were obtained by two different algorithms, and a smoother optimization result can be obtained by the present method, and the number of iterations of the present approach is less than the FEM. It also can be seen that the present approach can effectively avoid the checkerboard pattern.
4.2. Example 2 A simply supported beam with size 30 15 1 m3, as shown in Fig. 6, is now discussed. The two corner points on the bottom side are fixed, and a concentrated force F¼1 kN is loaded at the center of the bottom surface. The elastic material properties are chosen as Young’s modulus E¼30 Gpa, Poisson’s ratio m ¼0.3. The initial total weight of the structure is 450t. The initial rejection ratio is RR0 ¼0.01, the evolutionary ratio is ER¼0.01, and the volume constraint is limited to 30%. To determine an optimum structural layout, the design domain is discretized by EFG method using 31 16, 37 19, 43 22, 49 25 uniformly distributed nodes. And 30 15, 36 18, 42 21, 48 24 integration cells are used according to the
Fig. 9. Optimization result obtained by FEM with sensitivity filtering.
500
31 × 37 × 43 × 49 ×
450
Weight of the beam
400
16 nodes 19 nodes 22 nodes 25 nodes
350 300 250 200 150 100 0
5
10
15
20
25
30
35
Iteration number Fig. 8. Optimization result obtained by FEM without sensitivity filtering.
Fig. 10. Convergence history of the simply supported beam based on different nodes using the present method.
ARTICLE IN PRESS 672
Z. Juan et al. / Engineering Analysis with Boundary Elements 34 (2010) 666–672
location of the nodes and 2 2 Gauss quadrature is used. In Ref. [20], considering the minimization of compliance as objective function, and the volume of the structure as constraint, the same model is solved by the finite element method (FEM). The optimization results obtained by the present method is shown in Fig. 7, the optimization results obtained by FEM without sensitivity filtering is shown in Fig. 8, and the optimization results obtained by FEM with sensitivity filtering is shown in Fig. 9. From these results, it can be seen that similar optimization results were obtained by two different algorithms, and a smoother optimization result can be obtained by the present method, and the present approach can effectively eliminate the checkerboard pattern. Fig. 10 gives the curves of convergence of the simply supported beam’s objective function based on different nodes using the present method. The almost monotonic and uniform convergence can be observed from this figure. And for different nodes, their convergence characteristics are very similar.
5. Conclusions In this paper, the element free Galerkin method (EFG), combined with evolutionary structural optimization method (ESO), is applied to solve the topology optimization problem of the continuum structures. Considering the minimized weight as the objective function, the mathematical formulation of the topology optimization is developed based on the stress criterion. Several numerical examples are shown to prove the validity and feasibility of the present method. And the examples also show the simplicity and fast convergence of the proposed method.
Acknowledgements This work was supported by National 973 Scientific and Technological Innovation Project (2004CB719402), Natural Science Foundation of China (No. 10672055), The Key Project of NSFC of China (No. 60635020) and The National Science Foundation for Outstanding Youth of China (No. 50625519).
References [1] Bendsoe MP, Sigmund O. Topology optimization: theory, methods, and applications. Berlin, Heidelberg, New York: Springer; 2003. [2] Bendsoe MP, Kikuchi N. Generating optimal topology in structural design using a homogenization method. Computer Methods in Applied Mechanics and Engineering 1989;71:197–224. [3] Bendsoe MP, Sigmund O. Material interpolation schemes in topology optimization. Archive of Applied Mechanics 1999;69:635–54. [4] Xie YM, Steven GP. A simple evolutionary procedure for structural optimization. Computer and Structures 1993;49:885–96. [5] Xie YM, Stevel GP. Evolutionary structural optimization for dynamic problems. Computer and Structures 1996;58:1067–73. [6] Manickarajah D, Xie YM, Steven GP. An evolutionary method for optimization of plate buckling resistance. Finite Elements in Analysis and Design 1998;29:203–30. [7] Rong JH, Xie YM, Yang XY. Topology optimization of structures under dynamic response constraints. Journal Sound and Vibration 2000;234(2):l77–189. [8] Monaghan JJ. Smoothed particle hydrodynamics. Annual Review of Astronomy and Astrophysics 1992;30:543–74. [9] Belytschko T, Lu YY, Gu L. Element-free Galerkin method. International Journal for Numerical Methods in Engineering 1994;37:229–56. [10] Liu WK, Jun S, Li S. Reproducing kernel particle methods for structural dynamics. International Journal for Numerical Methods in Engineering 1995;38(10):1655–79. [11] Atluri SN, Zhu T. A new meshless local Petrov–Galerkin (MLPG) approach in computational mechanics. Computational Mechanics 1998;22:117–27. [12] Belytschko T, Krysl P, Krongauz Y. A three-dimensional explicit element-free Galerkin method. International Journal for Numerical Methods in Fluids 1997;24(12):1253–70. [13] Belytschko T, Gu L, Lu YY. Fracture and crack growth by element free Galerkin methods. Modelling and Simulation in Materials Science and Engineering 1994;2:519–34. [14] Liu GR, Chen XL. A mesh-free method for static and free vibration analyses of thin plates of complicated shape. Journal of sound and vibration 2001;241(5):839–55. [15] Cho S, Kwak J. Topology design optimization of geometrically non-linear structures using meshfree method. Computer Methods in Applied Mechanics and Engineering 2006;195:5909–25. [16] Zhou JX, Zou W. Meshless approximation combined with implicit topology description for optimization of continua. Structural and Multidisciplinary Optimization 2008;36(4):347–53. [17] Lu YY, Belytschko T, Gu L. A new implementation of the element free Galerkin method. Computer Methods in Applied Mechanics and Engineering 1994;113:397–414. [18] Krongauz Y, Belytschko T. Enforcement of essential boundary conditions in meshless approximations using finite elements. Computer Methods in Applied Mechanics and Engineering 1996;131:133–45. [19] Atluri SN, Kim HG, Cho JY. A critical assessment of the truly meshless local Petrov–Galerkin (MLPG), and local boundary integral equation (LBIE) methods. Computational Mechanics 1999;24(5):348–72. [20] Sigmund O. A 99 line topology optimization code written in Matlab. Structural and Multidiscipline Optimization 2001;21(2):120–7.