Strojniški vestnik - Journal of Mechanical Engineering 60(2014)1, 43-50 © 2014 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2013.945
Original Scientific Paper
Received for review: 2013-01-02 Received revised form: 2013-05-23 Accepted for publication: 2013-11-12
Topology Optimization for Continua Considering Global Displacement Constraint
Jijun YiU.* - Tao Zeng1 - Jianhua Rong2 1 Central South University, School of Mechanical and Electrical Engineering, China 1 Changsha University of Science & Technology, Key Laboratory of Lightweight and Reliability Technology for Engineering
Vehicle, China
This paper presents a topology optimization method for continua that minimizes the volume subject to global displacement constraint. The method uses the p-norm displacement to represent the equivalent maximum displacement so as to avoid non-differentiability of the maximum function. Using nodal densities as design variables, a new topology optimization technique for controlling the global maximum displacement precisely is developed. Several examples are presented to demonstrate the effectiveness of the proposed method for achieving convergent optimal solutions of structures with global displacement constraint.
Keywords: topology optimization, global displacement constraint, nodal density variable, p-norm displacement
0 INTRODUCTION
Finding the best distribution of available material in the predetermined design domain satisfying various conditions is the target of topology optimization for continuum structures. In most topology optimization methods, the optimized goal is to find the structures with maximum stiffness. Stiffness is often closely related to the global displacement and, especially, to the maximum displacement of the structure. So a new topology method for minimizing the volume of the structure subject to global displacement is developed and a new approach to controlling the maximum displacement of the structure is proposed.
The topology optimization method based on elements is treated by dividing the design domain into finite elements and each element is taken as a design variable. The solid isotropic material and penalization (SIMP) method [1] and [2] transfers design variables from the discrete variables to contiguous ones with penalization. The evolutionary structural optimization (ESO) method [3] and [4], and its later version, the bi-directional ESO (BESO) method [5], remove inefficient material from the structure based on certain predefined criteria. The level set method represents the structure using a level set model which is embedded in a scalar function. Rong and Liang [6] and Wang et al. [7] investigated the level set method for topology optimization. However, the element-wise topology optimizations exhibit various numerical problems, such as grey-scale, checkerboard pattern and mesh dependency. Therefore, topology optimization methods based on nodal design variables were developed to avoid these problems. Matsui and Terada [8] proposed the concept of continuous distribution. A
Q4/Q4 continuum structural topology optimization is investigated by Rahmatalla and Swan [9].
A 3-D structural topology optimization and novel surface-smoothing scheme based on SIMP and subelement bilinear interpolation was developed using node densities by Song and Kim [10]. An adaptive density point refinement approach for continuum topology optimization on the basis of an analysis-mesh separated material density field description based on nodal design variables was presented by Wang et al. [11]. Wang et al. [12] proposed topological optimization of structures using a multilevel nodal density-based approximant. The nodal density field using the non-local Shepard function method is transformed to a practical elemental density field via a local interpolation with the elemental shape function.
The presented topology optimization method is based on nodal densities and utilizes the rational approximation for material properties (RAMP) interpolation scheme proposed by Stolpe and Svanberg [13]. The discrete nodal topological variables pt that only take value 0 or 1 are replaced by continuous topological variables between 0 and 1. Consequently, the difficulty of the discrete optimization is avoided by penalization. It is assumed:
f (P) =
P
1 + v
( ^ )
(1)
where px is the density of point x, i.e. p(x), v is the penalization factor (v is a parameter which in some sense corresponds to p in the SIMP approach), which makes the intermediate densities approach either 0 (void) or 1 (solid). Here, v = 5 is set. The relation
*Corr. Author's Address: School of Automobile and Mechanical Engineering, Changsha University of Science and Technology, Changsha, China, jijunyi@gmail.com
43
Strojniski vestnik - Journal of Mechanical Engineering 60(2014)1, 43-50
between the Young's modulus and the material density at point x is expressed by:
E (x ) = f (P) E0
(2)
where E0 is Young's modulus of the full solid material. The functionf(px) has the following properties:
f ( P) = 0 as P ^ 0+
df (P)/ d P =
1
1 + v
^ 0 as P ^ 0+ •
The volume of an element is given by:
V =fV0 P dV,	(3)
where V is the volume of the ith element, V° is the original volume of the ith element.
In this study, minimum volume with a reference domain Q in R3 is considered while satisfying the global displacement constraint for the structure.
Find: p e RN Minimize: V
. —f
, (4)
Subjectto:uf < U (j = 1,2,...,Ndof;f = 1,2,...,L) p <pt < 1 (/ = 1,2,...,N„)
where V is the structural volume being optimized, uj is the displacement of the jth degree of freedom of the structure under the fth load case, Ndof is the total degrees of freedom of the structure, Uf is its constraint limit, L is the number of the load cases acting on the structure, and p¡ is the density of the ¡th node, p¡ is its lower limit, Nn is the total number of nodes. Here, the small positive lower bound p¡ = 0.0001 is set so that the structure optimized is always kept non-singular in the optimization process.
1 EQUIVALENT MAXIMUM DISPLACEMENT
When the maximum displacement of the structure does not exceed a specified value, the global displacement constraint is satisfied. So, the maximum displacement is naturally the ideal design criterion of optimization models. However, the location of the maximum displacement usually varies following the change of material distribution in the optimization process, so the maximum function is not differentiable. To solve this problem,the maximum displacement needs to be smoothed, and the /»-norm or the Kreisselmeier-Steinhauser (KS) function could be used. Similar to
the /»-norm, geometric average displacement (GAD) is introduced by Kreisselmeier [14] and Qiao and Liu [15], which is expressed as:
( 1 N,	\llp
U =
—Y í spdn
i L e=i e
(5)
Here, Ne is the total number of elements, Q denotes the volume of the design region, p is the norm parameter and Qe is the element e solid volume and the displacement of any point inside the element e is 8.
On one extreme, as the displacement norm parameter p approaches infinity, the Ua approaches the maximum displacement, and there is no added smoothness. On the other extreme, whenp approaches 1, there is excessive smoothness but Ua approaches the average displacement. A good choice for p should, therefore, provides adequate smoothness so that the optimization algorithm performs well and an adequate approximation of the maximum displacement value so that the optimized design satisfies the imposed displacement constrains.
Only when p approaches infinity, can the Ua arrive at the constraint of the maximum displacement. To remedy this deficiency, a normalized global displacement measure is proposed to better approximate the maximum displacement. The normalized p-norm displacement uses information from the previous optimization iteration to scale, and the p-norm displacement as a| Ua|, so that it better approximates the maximum displacement. The maximum displacement max(|Mk-1) and the p-norm displacement values from the previous optimization iteration k-1 are used to define our evolving normalized p-norm displacement constraint at each iteration k as:
max(|uj) < U	< U,
(6)
where a is calculated at each optimization iteration k > 1 as follows:
a =-
x(| ^
\k-K
max( |w| )
Pa
(7)
As the design converges, | Ua|k—1 ~ |Ua|k , max(|M|k-1) ~ max(U|k), so that ak-1 ~ ak and hence the desired effect is achieved, i.e. a| Ua| ~ U.
Note that the constraint a|Ua| < U is non-differentiable because the value of a is changing in a discontinuous manner and results in a slightly different optimization problem at every iteration. However, as the optimization converges the changes between successive design iterations diminish and hence a
k-i
44
Yi, J.J. - Zeng, T. - Rong, J.H.
Strojniski vestnik - Journal of Mechanical Engineering 60(2014)1, 43-50
converges, thereby reducing the effects of the non-differentiability and inconsistency.
2 NODAL DESIGN VARIABLES AND THE INTERPOLATION SCHEME
Element-independent nodal densities are the design variables in this method. Thus, the relative density at any point is interpolated with an interpolation scheme, which determines the topology, stiffness and volume of material.
2.1	Identifying the Nodes
As proposed by Guest et al. [16] and Kang and Wang [17], the scale parameter rmin is set to identify the nodes that influence the density of point x. Nodes are included in the influence domain if they are located within a distance rmin of the point x. This can be visualized by drawing a sphere of radius rmin centered at the point x, thus generating the spherical subdomain Qx. Nodes located inside Qx contribute to the computation of density p(x) of point x. As the mesh is refined, rmin and consequently Qx do not change. The only difference between the two meshes is the number of nodes located inside Qx, and included in the interpolation function. This is essential to generate mesh-independent solutions.
2.2	Shepard Interpolation Scheme
Interpolation provides a continuum of density field and mesh-independence, which might alleviate numerical instability and checkerboard effects [18]. In implementing continuum structural topology optimization formulations, many functions are available to interpolate nodal density onto the points inside the element space; -for example, the standard C0 shape functions used in the finite element method. However, each node's shape function influences only the elements connected to that node. Mesh independency cannot be obtained when interpolation functions are influenced by mesh size. They should be based on a physical length scale that does not change following mesh refinement.
Shepard interpolation is proposed by Shepard [19], and used by Kang and Wang [17] to achieve mesh independency. Let pi (i=l,2,...,n) denote a set of density of nodes inside Qx at the associated point x= (X, Y, Z), where (X, Y, Z) define the point x location in the Cartesian coordinate system. Thus, the relative density at point x is interpolated by the nodal densities
inside the influence domain Qx with Shepard's interpolation method.
px =
£ NDi (x)pi (x * x )
ieSx
P
(8)
(x = x )
where Sx is the sub-domain of design variable located within the influence domain Qx of point x, and pi is the density value of the ith node. xi is the position of the point associated with the ith node. The corresponding interpolation function NDi(x) is defined as:
ND, ( x) =
R ( x) ÊR ( x)
(i = 1,2,..., n).
(9)
where Ri(x) = 1 / [r^x)] and ri(x) = |x - xi|2 being the Euclidean distance from the points x to x,-. n is the number of nodes inside the influence domain Qx.
In the method of element independent nodal variable density, the density in the element space is not constant, and the global density field of the structure has C0 continuity. It is easy simple to know from the bounded property of the Shepard interpolation that 0 <px< l holds if 0 <pi< l (i e Sx). Moreover, the property NDi(x) > 0 also guarantees that the derivative of the density with respect to the design variable will be always non-negative. This property is essential to guaranteeing a correct searching direction in seeking the optimal material distribution by a gradient-based algorithm.
3 SENSITIVITY ANALYSIS
The solution of the gradient-based optimization problem requires the computation of sensitivities of the objective function and the constraints. In a finite element analysis, the static behavior of a structure for any load case can be expressed by the following equilibrium equation:
K U = F
(10)
where K is the global stiffness matrix of a structure being optimized and, U and F are the global nodal displacement and nodal load vector, respectively.
In the finite element method, the displacement of any point S can be expressed by:
ô = N u
(11)
Topology Optimization for Continua Considering Global Displacement Constraint
45
Strojniski vestnik - Journal of Mechanical Engineering 60(2014)1, 43-50
Here, N, u denote the general shape function matrix and displacement vector of the /'th element, respectively.
The adjoint method can be used to determine the sensitivity of displacement by introducing a vector of Lagrange multiplier X. The modified /»-norm displacement can be expressed as:
( 1 N,
U„ =
y p
iÉUNuYdn, +XT (KU - F), (12)
^ e=1 e	y
where the term kT (KU - F) is equal to zero and, therefore, the modified displacement is identical to the original one. Taking derivatives of Eg. (12) with respect to the design variable p gives:
dUa 1 N.1-P 1 -r-s r	ou ou „
= ~(U°) ÖP(Nu) N——dQ +
dp p	Q '	dU dp
du dU
t'ôkU+K dU,.
^ dp dp The Eq. (13) can be rewritten as:
(13)
dp
1U ) '	(N. N dU ,O + K +
5U
dKu.
dp
(14)
To eliminate the unknown dU / dp from the sensitivity expression, let:
(U¿§J>ur N^dne + XtK = 0, (15)
p,=-u r i ±(dn
t!
dU
dU
Jn(Nu NTd
So the adjoint vector is defined as: K! = Px.
. (16)
(17)
Here, PA denotes the adjoint load of adjoint Eq. for yielding the adjoint vector k.
Thus, the sensitivity of the /»-norm displacement function is:
dUa _dUa dp _2 NC^t dKe
d(i/p) dp d(i/p)
= -p2 EM
dp
(18)
The sensitivity of the /-norm displacement requires the computation of the sensitivity of the stiffness matrix with respect to the design variable. The derivative of the elemental stiffness matrix with respect to the design variable is expressed by:
dK.
df (P )„ r
dp dp
df (px ) _ df(P ) p __ dp
B' D0BJ Qe
(1 + v)
dp dp ((1 -p )v +1)
(19)
-NDt, (20)
where B is the usual displacement-strain matrix and D0 corresponds to the constitutive matrix of the solid material.
For example, the formulation of the constitutive matrix for 3D isotropic solid structures is:
= _Em-jiL x (1 + v)(1 - 2^)
ß	0	0
1 -ß		
ß	0	0
1 -ß		
1	0	0
	1 - 2ß	0
	2(1-ß)	
		1 - 2ß
2(1-ß)
0 0 0
0
1 - 2ß 2(1 - ß)
. (21)
Since the stiffness matrix integrand is evaluated at the Gauss points, the densities at these Gauss points are directly computed from the design variables using the interpolation function. Numerical quadrature, such as Gaussian quadrature, is commonly reduced to the evaluation and summation of the stiffness integrand at specific Gauss points.
The sensitivity analysis of the objective function in Eq. (4) is calculated similar to that of Eq. (19). The derivative of the total material volume with respect to the design variables can be computed by the Gaussian quadrature method over the influence-domain.
=jnPdne = Jq X NDi (x)Pi dQe. (22)
1
X
e
where Nc is the number of element influenced by pt. K0 is the element stiffness matrix of the eth element of the solid material.	V = XV = XJ X NDi (x)pidD.e. (23)
e=1	e=1 e i^Sx
46	Yi, J.J. - Zeng, T. - Rong, J.H.
So the total volume is:
Nel	Nei
Strojniski vestnik - Journal of Mechanical Engineering 60(2014)1, 43-50
The derivative of V with respect to the design variable:
=	= -(p )2 Z(f	).(24)
0(1/pi) dp,, d(1 p,) ÙS'
4 TOPOLOGICAL OPTIMIZATION WITH VARYING DISPLACEMENT LIMITS
In order to make the approximation functions of the displacement constraint in Eq. (4) hold true at each iteration, an equivalent optimization model (Eq. (25)) with varying displacement constraint limits is built. At each iteration, a quadratic approximation model of the true objective function that satisfies the Taylor expansion is built around the current point. The model is assumed to be a good representative of the objective function in a so-called trust region [20]. Trust regions are used to ensure the robustness of the iteration and make progress toward feasibility and optimality [21].
min V = TjJn PdQe
s.t. aU{ <U( (f = 1,2,-,L; I = 1,2,-), (25)
Pi <Pi < 1
(i = 1,2,-, Nn )
U f =
where U( is expressed as:
Uf | + min(ß |uf |,(Uf - \ f |)), |uf | < Uf Ufa\ - min (ß |uf|, \Uf - \Ufa\ [J, |uf| > Uf ' f = 1,2,..., L; I = 1,2,...	(26)
In the above, ß is a displacement limit changing factor. Typical values of ß between 0.01 and 0.20 have been used for displacement constraints in the example problems in this paper. Ufak is the /-norm displacement of the structure under the fh load case. U( (j = 1,2,..., J) are varied by Eq. (26) at every iteration.
Assuming that only t, = lp is changed and is treated as a variable, the first-order series expansion for the /»-norm displacement Ua at t, (i = 1, 2, ..., Nnod) can be expressed as:
Ua = Ua' - £ (Ua tf + £ (U a t, , (27)
i=1 ' i=1
where Nn is the number of the nodes in the design domain.
Thus, the /»-norm displacement in the next iteration, Ufk+1, can be estimated by the current iteration k.
Nn
Let C =YUa' ¿sign (Ufak) (i= 1,-, Nn) and
i=1 i
Af = £(Uaresign(ufk) (i = 1,-,Nn), then Eq.
i=1 '
(25) can be transferred into Eq. (28):
min: K = £ ¡D 1
e=1 e
Nn	Nn
s.t. Yj Cfti < Ufk - Ufksign(Ufk ) + Y Af. (28)
i=1 i=1
(f = 1,2,-L; I = 1,2,-) 1 < ti < ^ (i = 1,2,-,Nn)
If the constant items in the objective function are omitted, solving Eq. (28) can be transferred to solving Eq. (29):
Find t e RN
Nn
mm
in: Y(b (ti )2 + at )
i=1
Nn
s.t. Y Cfti < ul - Ufksign(Ufk ) + Y A
i=1 i=1
1 < t, < ~t,
Nn
(29)
where
-3
' (t )
Y (i ND,dne ), b = ^ Y (f
ND.dQ. I.
The programming solving the problem of Eq. (29) can be transferred into solving dual programming problem by using the dual theory [22] and [23].
5 NUMERICAL EXAMPLES
This section illustrates the proposed approach with 3D applications. For simplicity, all the quantities are dimensionless. In addition, Young's modulus is chosen as 2.1*10n and Poisson's ratio as 0.3 for all examples. Let d denote the length of cuboid element diagonal. In the following examples, rmm is set to 1.5d while keeping the same displacement mesh size. As the computing platform, we have used a personal computer with the commercial software package ABAQUS has been used for FEA in this study.
Fig. 1a shows a 3D cantilever beam with length of 8, height of 5, and width of 2. The beam is fixed at
Topology Optimization for Continua Considering Global Displacement Constraint
47
Strojniski vestnik - Journal of Mechanical Engineering 60(2014)1, 43-50
the left end. A concentrated load P = 4800 is applied downward at the center of the right end. The initial displacement at the constraint point, the center of the right end, is 2.851x10-6. A displacement constraint is taken as 7.0x10-6. The cantilever domain is discretized into a mesh of 40x25x10 B8 elements which results in a total of 10,000 elements. There are 41x26x11 design variable points distributed within the initial design domain. The optimized topology is shown in Fig. 1b.
a) it
b)
Fig. 1. The 3D cantilever beam; a) design domain and boundary conditions; b) topology optimization result obtained from the proposed method
Then, the same problems are solved using the proposed method with varying p parameter in which p is set to even numbers (4, 10, 20, 40) as shown in Fig. 2. With the increase of p parameter, the displacement constraint converges to the same condition (Fig. 3), the material distribution (Fig. 2) and the volume (Fig. 4) move close to the convergence.
a)
d)
Fig. 2. Topology optimization results with: a) p = 4; b) p = 10; c) p = 20; d) p = 40
Fig. 5 shows a 3D, simply supported beam with a length of 4, height of 1, and width of 0.4. An external force P = 4000 is applied to the center of top area. The initial displacement at the constraint point, the center of the top area, is 1.1x10-7. A displacement constraint is taken as 2.0x10-7. The design domain is discretized into mesh size of 80x20x8 B8 elements. p is set to 20.
Fig. 3. The maximum displacement-varying history
Iteration
Fig. 4. The volume-varying history

w
Fig. 5. The design domain and boundary conditions of the 3D, simply supported beam
a)	~ b)
Fig. 6. Topology optimization results: a) optimal topology obtained from EIND; b) optimal topology obtained from the element-wise method
a)	b)
Fig. 7. Topology optimization results with the coarser mesh: a) optimal topology obtained from EIND; b) optimal topology obtained from the element-wise method



c)
48
Yi, J.J. - Zeng, T. - Rong, J.H.
Strojniski vestnik - Journal of Mechanical Engineering 60(2014)1, 43-50
When the nodal density values are used to determine a smooth iso-line that describes the boundary of the optimization layout, as a result, a smooth optimal topology can be obtained in Fig. 6a. Its final volume is 0.947. The results obtained from the element-based approach is shown in Fig. 6b. Its final volume is 1.062. When the mesh is coarser, and the mesh size of 50*12*5 is used, the results obtained from the element-based approach and the proposed approach based on element independent nodal density are shown in Figs. 7a and 3b, respectively. Their final volumes are 1.071 and 1.155.
These figures show that, for the different displacement mesh size, the topology obtained from element independent nodal density has a much better resolution and is smoother than that of the element-based approach. The solution attained by the proposed method exhibits no checkerboard patterns or mesh dependency.
When the same mesh is used, the computational cost for the topology optimization based on element independent nodal density is higher than the element-based approach. This is mainly attributable to the large number of density nodes in the influence domain. However, the topology resolution resulting from the proposed approach based on the proposed method is higher than that of the element-based approach. To improve the efficiency of the proposed approach, especially for a 3D large-scale optimization problem, the parallel programming technique could be used to carry out the finite element analysis and the optimization procedure.
The total CPU time and the CPU time spent on the sensitivity analysis in every optimization iteration are 144 and 123 seconds, respectively. When the codes are reprogrammed with the OpenMP and four threads are used, the total CPU time and the CPU time spent on the sensitivity analysis decrease to 50 and 42 seconds, respectively. By taking this approach, it is possible to obtain benefits from parallelization without the need for extensive modification to the code structure.
6 CONCLUSIONS
This paper has developed a topology optimization method for minimizing the volume of a structure subject to the global displacement constraint. In contrast to the element-based procedure, here we take the nodal density as the design variable, which is interpolated into any point by Shepard functions. This technique avoids checkerboard patterns and mesh-dependency for low order finite elements. With the
help of the global displacement constraint, an optimal structure with appointed deformation can be obtained, and it is unnecessary to know where the maximum displacement is. The proposed method is highly useful with regard to practical engineering applications. The numerical examples demonstrate the effectiveness of the proposed method with respect to the optimal solution and convergence.
7 ACKNOWLEDGEMENTS
This work has been supported by National Natural Science Foundation of China (11372055, 11302033), Hunan Provincial Natural Science Foundation of China (12JJ3044), the Key Laboratory of Lightweight and Reliability Technology for Engineering Vehicle, Education Department of Hunan Province (Changsha University of Science & Technology) (2012KFJJ0 2), the Huxiang Scholar Fund (Changsha University of Science & Technology), the Open Fund of State Key Laboratory of Automotive Simulation and Control (20121105).
8 REFERENCES
[1]	Bendsoe, M.P. (1989). Optimal shape design as a material distribution problem. Structural and Multidisciplinary Optimization, vol. 1, no. 4, p. 193202, D0I:10.1007/BF01650949.
[2]	Rozvany, G., Zhou, M., Birker, T. (1992). Generalized shape optimization without homogenization. Structural and Multidisciplinary Optimization, vol. 4, no. 3, p. 250-252, D0I:10.1007/BF01742754.
[3]	Xie, Y.M., Steven, G.P. (1993). A simple evolutionary procedure for structural optimization. Computers & Structures, vol. 49, no. 5, p. 885-896, D0I:10.1016/0045-7949(93)90035-C.
[4]	Xie, Y.M., Steven, G.P.(1997). Evolutionary Structural Optimization. Springer, Berlin, DOI:10.1007/978-1-4471-0985-3.
[5]	Huang, X., Xie, Y.M. (2010). Evolutionary Topology Optimization of Continuum Structures: Methods and Applications. John Wiley & Sons, Chichester, D0I:10.1002/9780470689486.
[6]	Rong, J.H., Liang, Q.Q. (2008). A level set method for topology optimization of continuum structures with bounded design domains. Computer Methods in Applied Mechanics and Engineering, vol. 197, no. 1718, p. 1447-1465, D0I:10.1016/j.cma.2007.11.026.
[7]	Wang, M.Y., Wang, X., Guo, D. (2003). A level set method for structural topology optimization. Computer Methods in Applied Mechanics and Engineering, vol. 192, no. 1, p. 227-246, D0I:10.1016/S0045-7825(02)00559-5.
[8]	Matsui, K., Terada, K. (2004). Continuous approximation of material distribution for topology
Topology Optimization for Continua Considering Global Displacement Constraint
49
Strojniski vestnik - Journal of Mechanical Engineering 60(2014)1, 43-50
optimization. International Journal for Numerical Methods in Engineering, vol. 59, no. 14, p. 1925-1944, DOI:10.1002/nme.945.
[9]	Rahmatalla, S.F., Swan, C.C. (2004). A Q4/ Q4 continuum structural topology optimization implementation. Structural and Multidisciplinary Optimization, vol. 27, no. 1, p. 130-135, DOI:10.1007/ s00158-003-0365-9.
[10]	Song, J.-H., Kim, C. (2012). 3-D topology optimization based on nodal density of divided sub-elements for enhanced surface representation. International Journal of Precision Engineering and Manufacturing, vol. 13, no. 4, p. 557-563, DOI:10.1007/s12541-012-0071-x.
[11]	Wang, Y., Kang, Z., He, Q. (2013). An adaptive refinement approach for topology optimization based on separated density field description. Computers & Structures, vol. 117, p. 10-22, DOI:10.1016/j. compstruc.2012.11.004.
[12]	Wang, Y., Luo, Z., Zhang, N. (2012). Topological optimization of structures using a multilevel nodal density-based approximant. Computer Modeling in Engineering and Sciences, vol. 84, no. 3, p. 229, DOI:10.3970/cmes.2012.084.229.
[13]	Stolpe, M., Svanberg, K. (2001). An alternative interpolation scheme for minimum compliance topology optimization. Structural and Multidisciplinary Optimization, vol. 22, no. 2, p. 116-124, DOI:10.1007/ s001580100129.
[14]	Kreisselmeier, G. (1979). Systematic control design by optimizing a vector performance index. International Federation of Active Control Symposium on Computer-Aided Design of Control Systems, p. 113-117.
[15]	Qiao, H., Liu, S. (2013). Topology optimization by minimizing the geometric average displacement. Engineering Optimization, vol. 45, no. 1, p. 1-18, DOI: 10.1080/0305215X.2012.654789.
[16]	Guest, J.K., Prevost, J.H., Belytschko, T. (2004). Achieving minimum length scale in topology optimization using nodal design variables and projection functions. International Journal for Numerical Methods in Engineering, vol. 61, no. 2, p. 238-254, DOI:10.1002/nme.1064.
[17]	Kang, Z., Wang, Y.Q. (2011). Structural topology optimization based on non-local Shepard interpolation of density field. Computer Methods in Applied Mechanics and Engineering, vol. 200, no. 49, p. 35153525, DOI:10.1016/j.cma.2011.09.001.
[18]	Diaz, A., Sigmund, O.(1995). Checkerboard patterns in layout optimization. Structural and Multidisciplinary Optimization, vol. 10, no. 1, p. 40-45, DOI:10.1007/ BF01743693.
[19]	Shepard, D. (1968). A two-dimensional interpolation function for irregularly-spaced data. Proceedings of the 1968 23rd ACM national conference, p. 517-524, DOI:10.1145/800186.810616.
[20]	Conn, A.R., Gould, N.I., Toint, P.L. (1987). Trust Region Methods. Vol. 1. Society for Industrial and Applied Mathematics, Philadelphia.
[21]	Byrd, R.H., Gilbert, J.C., Nocedal, J. (2000). A trust region method based on interior point techniques for nonlinear programming. Mathematical Programming, vol. 89, no. 1, p. 149-185, DOI:10.1007/PL00011391.
[22]	Beckers, M. (1999). Topology optimization using a dual method with discrete variables. Structural and Multidisciplinary Optimization, vol. 17, no. 1, p. 14-24, DOI:10.1007/BF01197709.
[23]	Rong, J.H., Li, W.X., Feng, B. (2010). A structural topological optimization method based on varying displacement limits and design space adjustments. Advanced Materials Research, vol. 97-101, p. 36093616, DOI: 10.4028/www.scientific.net/AMR.97-101.3609.
50
Yi, J.J. - Zeng, T. - Rong, J.H.