- Email: [email protected]
Evolutionary topology optimization of continuum structures with a global displacement control
Computer-Aided Design 56 (2014) 58–67
Contents lists available at ScienceDirect
Computer-Aided Design journal homepage: www.elsevier.com/locate/cad
Evolutionary topology optimization of continuum structures with a global displacement control✩ Zhi Hao Zuo, Yi Min Xie ∗ Centre for Innovative Structures and Materials, School of Civil, Environmental and Chemical Engineering, RMIT University, GPO Box 2476, Melbourne 3001, Australia
highlights • • • •
The proposed method achieves optimal designs with globally controllable deflections. Solutions with the proposed method can maintain aerodynamic shape when deformed. Economic lightweight designs are achieved with a volume minimization scheme. Robust and effective numerical algorithms are proposed for clear 0–1 designs.
article
info
Article history: Received 17 August 2013 Accepted 14 June 2014 Keywords: Topology optimization Displacement control Bi-directional evolutionary structural optimization (BESO) Aircraft design Aerodynamics
abstract The conventional compliance minimization of load-carrying structures does not directly deal with displacements that are of practical importance. In this paper, a global displacement control is realized through topology optimization with a global constraint that sets a displacement limit on the whole structure or certain sub-domains. A volume minimization problem is solved by an extended evolutionary topology optimization approach. The local displacement sensitivities are derived following a powerlaw penalization material model. The global control of displacement is realized through multiple local displacement constraints on dynamically located critical nodes. Algorithms are proposed to secure the stability and convergence of the optimization process. Through numerical examples and by comparing with conventional stiffness designs, it is demonstrated that the proposed approach is capable of effectively finding optimal solutions which satisfy the global displacement control. Such solutions are of particular importance for structural designs whose deformed shapes must comply with functioning requirements such as aerodynamic performances. © 2014 Elsevier Ltd. All rights reserved.
1. Introduction Topology optimization is a powerful tool resulting in high structural performance with great material saving. It is capable of speeding up the structural design process and producing reliable solutions for structural design problems. This area has been extensively investigated in the past three decades since the modern formulation of optimal layout theories by Prager and Rozvany [1]. Several Finite Element (FE) based methods have been developed for topology optimization of continuum structures. One most popular technique is the Solid Isotropic Material with Penalization (SIMP) method [2,3], where the isotropic material property in each
✩ This paper has been recommended for acceptance by Michael Yu Wang.
∗
Corresponding author. Tel.: +61 3 99253655; fax: +61 3 96390138. E-mail address: [email protected] (Y.M. Xie).
http://dx.doi.org/10.1016/j.cad.2014.06.007 0010-4485/© 2014 Elsevier Ltd. All rights reserved.
element is determined by the element relative density. By penalizing the element relative density using a power-law interpolation scheme, a nearly void–solid design is expected after an iterative procedure such as the Method of Moving Asymptotes [4]. A popular group of methods for topology optimization is the Evolutionary Structural Optimization (ESO) method and its descendent versions. The ESO method was first proposed by Xie and Steven [5,6] in the early 1990s. This original version follows a straightforward algorithm of iteratively removing inefficient material and thus drives the structure to evolve towards an optimum. The bi-directional evolutionary structural optimization (BESO) method was proposed as an extended version of ESO in the late 1990s [7,8]. The BESO method allows material to be added to the most demanding places while inefficient material is removed simultaneously. The ESO/BESO methods use binary design variables and directly deliver void–solid solutions that are very desirable. The ESO/BESO methods can be easily implemented in a wide range of engineering applications such as [9–13].
Z.H. Zuo, Y.M. Xie / Computer-Aided Design 56 (2014) 58–67
Deflection constrained optimization problems are related with stiffness optimization problems [14] to a large extent. These two sets of problems are in fact equivalent under single point loaded systems and with deflection control at the load point in the load direction. However generally, the stiffness optimal design is not the same with the optimal deflection design. The stiffness optimization increases the structural overall stiffness and thus the deflection is roughly reduced as an indirect effect. Very often in practice, it is not the stiffness but the deflection over parts of or whole of the structure that is the crucial factor. Despite vast research on stiffness optimization, the local displacement constraint on continuum structural optimization is found in relatively less literature such as [15–17]. For some structures such as an aircraft wing, the exterior surface should undergo minimal shape change under deformation in order to maintain the aerodynamic performance [18]; this requires that the displacements of the surface shall remain within a certain limit. In such cases, the displacements of a group of local nodes are of concern; or more generally, the displacement limit is addressed as a global constraint. The deflection constrained optimization problem addressed in the present paper is a global control problem. The common way of obtaining the displacement sensitivity is to apply a unit virtual load on the original model and get the displacement vector from the virtual system [16]. For global displacement control, the maximum deflections occur unpredictably at different nodes through the optimization history. Therefore the virtual loads must be built dynamically after the maximum deflections are identified in the real system at each stage; and the real and virtual systems must be analysed in two subsequent FEAs. This increases the computational cost significantly. More importantly, numerical instabilities caused by dynamical relocation of control nodes are hard to handle; as a consequence, stable solution convergence is highly difficult to achieve. Multiple displacement constraints optimization can be found in literature such as [19–23]. The common strategy of this art is to control the local displacement at pre-defined locations, usually at the load point and in the load direction, such as in [19–21]; in other words, a modified compliance minimization problem is actually dealt with. Pre-defined displacement constraints are also applied to specific structural systems such as truss structures [22,23]. The complexity of solving multiple constraints is reduced due to the fact that less design variables are present in the design domain. Since the displacement control locations and directions are limited (by pre-definition or by the applied load), a genuine ‘‘global control’’ for displacements is not realized in these methods. In some cases, the critical displacement may occur at a location completely different from the load points; such fixed multi-constraint approaches will not be able to tackle the maximum displacement. Multi-constraint optimization for continuum structures is usually a highly non-linear problem that is theoretically very difficult to solve, especially in cases as the global displacement constraint that controls displacements on dynamic locations. This paper presents a global control method for displacements of continuum structures. A topology optimization problem of volume minimization is formulated and solved by a new BESO approach. The global displacement constraint directly imposes a maximum allowable limit for nodal displacements within the whole domain or some user-specified sub-domains of the structure. The locations and directions of maximum nodal displacements are dynamically detected, the numerical instabilities of which are adaptively dealt with by robust stabilization algorithms. The subsequent sections are organized as follows: Section 2 formulates the topology optimization problem for global displacement control; Section 3 presents the sensitivity calculation for displacements; Section 4 is dedicated to the global displacement constraint algorithms for stabilizing the evolution history and solution convergence; the numerical implementation of the proposed method is outlined in this section; Section 5 shows numerical examples with discussions; the conclusions are drawn in Section 6.
59
2. Problem statement The global displacement control can be realized by confining the maximum displacement of the structure. With the volume being the objective function, the problem is formulated as follows. Minimize V , subject to V =
N
Vi x i ,
xi ∈ {xmin , 1}
(1)
i =1
dj ≤ dmax ≤ d∗ ,
j = 1, . . . , L
P = KU
(2) (3)
where V is the volume of the structure, Vi is the volume of the ith element, N is the total number of elements; xi is the binary design variable, i.e. either 1 or 0 denoting the element status of ‘‘solid’’ or ‘‘void’’; the design variable can also be treated as the element relative density [16] with a very small value such as 10−6 representing the void status; d∗ is the maximum allowable displacement value, dmax is the maximum displacement in the structure under loading, and dj is the displacement of the jth element within a prescribed domain Ωdisp where the displacement constraint is active. Eq. (3) is the static equilibrium with P as the applied load vector, K as the global stiffness matrix and U as the global displacement vector. The displacements above can be addressed in any certain directions, or can be treated as relative displacements simply when some nodes in Ωdisp are fixed. Constraining the relative displacements is of particular importance when the deformed shape is concerned. 3. Displacement sensitivity calculation 3.1. Material interpolation scheme Material interpolation schemes usually express the material properties as a function of the design variables in order to facilitate the sensitivity analysis. In the power-law material model SIMP [24], the Young’s modulus of an element is determined by the element relative density through the following penalization formulation. p
E (xi ) = xi E 0
(4)
0
where E is the Young’s modulus of the solid material and p denotes the penalty exponent which is usually set as p ≥ 3. Note that the Poisson’s ratio is usually assumed to be irrelevant in SIMP. Therefore, the stiffness matrix can be expressed in a similar way as follows. K=
p
xi K0i
(5)
i
where K0i denotes the stiffness matrix for the solid element, i.e. when xi = 1. 3.2. Displacement sensitivity in axial directions The displacement kth component can be obtained by multiplying the displacement vector with a unit virtual load vector Fk , of which the kth component is unity while all other components are zero. uk = Fk · u
(6)
F = {0, 0, . . . , 1, 0, . . .} = {f1 , f2 , . . . , fk , fk+1 , . . .}. k
k
(7)
With the virtual load F and the applied load P being constant, differentiating the kth displacement component with respect to
60
Z.H. Zuo, Y.M. Xie / Computer-Aided Design 56 (2014) 58–67
the ith element and substituting Eq. (3) gives the following k,T
The sensitivity of the absolute displacement is then immediately derived as
k,T
∂u ∂u ∂F ∂F = u + Fk,T = u ∂ xi ∂ xi ∂ xi ∂ xi −1 ∂P ∂K P + K−1 + Fk,T ∂ xi ∂ xi k
p−1
αi = −
u21 + u22 + u23 p−1
∂ K −1 ∂K ∂K = Fk,T P = −Fk,T K−1 u = −uk,T u ∂ xi ∂ xi ∂ xi
(8)
where uk = K−1 Fk is system response (displacement field) under the unit virtual load. Substituting the material interpolation in Eq. (5) into the above derivative gives finally the element sensitivity for the kth displacement vector with respect to the ith element.
αi =
pxi
∂u = −pxpi −1 uik,T K0i ui ∂ xi k
(9)
k,T
where ui and ui are the element displacement vector under the unit virtual load and the real load conditions respectively. By setting the penalty exponent p as infinity, assigning an element with the minimum relative density is equivalent to removing that element (E = 0), and the sensitivity is simply zero. In such a ‘‘hard-kill’’ manner, the model size is reduced and the computation cost can be saved significantly.
= −
pxi u21
+
u22
1 ,T
u1 ui
+ u2 u2i ,T + u3 u3i ,T K0i ui
k,T
+
u23
ui K0i ui .
(14)
3.4. Filtering sensitivities for mesh-independence Topology optimization usually encounters chequerboard and mesh-dependence problems [25]. To circumvent these problems, a mesh-independency filter [26] is used by averaging the element sensitivities with its neighbourhood elements using a low-pass filter of radius rmin . Firstly, nodal sensitivity numbers αjn which do not carry any physical meaning on their own are calculated by evenly extrapolating the sensitivities of connected elements; then these nodal sensitivity numbers are converted back into elements using a weight function w based on rij —the distance between the centre of the element i and node j, and on the filter radius rmin , in the following.
w(rij )αjn , αˆ i = w(rij )
3.3. Displacement sensitivity in arbitrary directions The displacement sensitivity in non-axial directions can be derived based on the displacement sensitivity in axial directions. The real displacement at the kth node is expressed by a vector u = {u1 , u2 , u3 } with u1 , u2 and u3 being the components in the x, y and z axial directions. Then the projective displacement on a specified direction e = {e1 , e2 , e3 } can be treated as the function of the three components. u˜ = u · e = u˜ (u1 , u2 , u3 ) = u1 e1 + u2 e2 + u3 e3 .
(10)
The derivative of u˜ is obtained as the dot product of the gradient of u˜ and the derivative of u.
∂ u˜ ∂ u1 ∂ u˜
∂ u1 ∂ xi ∂u
∂ u˜ ∂u 2 = ∇ u˜ · = · . ∂ xi ∂ xi ∂ x ∂ u i 2 ∂ u3 ∂ u˜ ∂ xi ∂ u3
(11)
Substituting Eq. (9) into the above equation, we finally obtain the displacement sensitivity in the specified direction e defined by {e1 , e2 , e3 }.
u1i ,T K0i ui e 1 ∂ u˜ = −pxpi −1 e2 · u2i ,T K0i ui αi = ∂ xi 3,T 0 e3 ui Ki ui = −pxpi −1 e1 ui1,T + e2 u2i ,T + e3 u3i ,T K0i ui (12)
e=
u21 + u22 + u23
,
u2 u21 + u22 + u23
,
The local displacement sensitivity in Section 3 forms the basis for global displacement control. To realize such a global control, straightforwardly, all the displacements within the concerned domain Ωdisp shall be addressed, i.e. the global displacement constraint in Eq. (2) shall be satisfied. 4.1. From local control to global control Ideally, the global displacement constraint is satisfied if the globally maximum displacement dmax is within the allowable limit d∗ . In other words, the global displacement constraint might be simplified to constraining the maximum displacement d∗ over the set Ωdisp . In such a ‘‘worst case’’ algorithm, the location of the maximum displacement might vary from iteration to iteration. As a consequence, the displacement of other nodes might increase significantly to become the maximum displacement. The jumping control node will lead to numerical instabilities and make the solution difficult to converge. Therefore, compared to controlling only one node of the maximum displacement, it is more reasonable to control a set of nodes whose displacements are close to dmax . when
dmax − dn dmax
≤t (16)
Note that according to the linear superposition of displacements, uki is the system response under the virtual load Fk = {e1 , e2 , e3 }. If the unit vector e is in the direction of the displacement u, namely u1
(15)
4. Global displacement control
dn ≤ dn+1 ≤ · · · ≤ dn+Q −2 ≤ dmax ≤ d∗ ,
= −pxpi −1 uik,T K0i ui .
where w = max(rmin − rij , 0).
u3
u21 + u22 + u23
.
(13)
where Q is the number of nodes to be controlled; t is a pre-defined tolerance that is set to a relatively small value such as 5% used in the current work. Increasing the tolerance t means putting more nodes into control which increases the smoothness of the evolution history. However the final optimum might degrade since the maximum displacement is not of enough focus. To solve this issue, the individual virtual loads are calculated using a weighted function that reflects their priorities based on the relevant displacements.
Fj = qj |1| = dj . dmax
(17)
Z.H. Zuo, Y.M. Xie / Computer-Aided Design 56 (2014) 58–67
Note that this mechanism is equivalent to formulating a global objective function for displacement [27] which is composed by local displacement functions with weighted factors qj =
dj dmax
.
61
threshold can be easily determined using a simple and efficient bisection algorithm. Starting from the full design of the structure, the target volume of each iteration is progressively reduced as described below. V k = V k−1 (1 − ER · f )
4.2. Numerical stabilization
(20)
where V is the target volume of the current iteration and V k−1 is the volume of the previous iteration; ER is the user-prescribed evolutionary volume ratio which determines the volume variation between two adjacent iterations, e.g. ER = 2% means the remaining volume will be changed by 2% (if f = 1.0) for each iteration; note that with a greater ER the volume will be reduced faster, but at the same time numerical instabilities may also arise that may require more iterations for the solution to converge; usually ER is selected between 1% and 5%; f is a function indicating the relative gap between the current maximum displacement and the imposed constraint value d∗ . f serves as a factor to adaptively fine-tune the volume change step size in order to stabilize the optimization objective and constraints. k
Still, the nodes under control may still be subject to change through the optimization history. The sensitivity history-averaging [26] can be applied to further stabilize the evolution. However, a refined modification shall be made here to achieve more reasonable smoothing. Instead of simple averaging [26], a weighted combination of sensitivities of two adjacent iterations is introduced.
α = w1 α k + w2 α k−1 where w1 + w2 = 1.
(18)
Here α is the completed element sensitivity after smoothing, α k and α k−1 are the sensitivities of the current and the immediately previous iteration, w1 and w2 are the weighted factors for both items. During the optimization procedure, the locations of the maximum displacements, i.e. the control nodes of two adjacent iterations might vary significantly. Thus the evolutionary history may suffer from numerical instabilities caused by the sudden change of sensitivities. The usage of the sensitivity weighted factors is able to stabilize the optimization procedure by setting w1 and w2 adaptively. As simple sensitivity averaging [26], the initial values can be set as w1 = w2 = 0.5; when the maximum displacement locations changes significantly, it is helpful to set w2 greater than w1 in order to avoid excessive disturbance of the optimization trend. The weighted factors w2 can take the following general form:
w2 = 0.5 +
q Q
b
(19)
where q (q < Q ) is the number of new nodes in the total Q control nodes within Ωdisp , and b is a user-prescribed gap factor chosen between 0 and 0.5. Note that the greater b is chosen, the less the solution will suffer from the numerical instabilities, but it might take more iterations for the solution to converge to the final optimum; 0.3 is an appropriate candidate for b. 4.3. Element update and numerical stabilization The element sensitivity indicates the changing trend of the maximum global displacement due to changing the design variable of one element. Changing the design variable xi from 0 to 1 for one element with negative sensitivity decreases the maximum global displacement, while changing xi from 1 to 0 of one element with positive sensitivity decreases the maximum global displacement. For optimization methods using continuous design variables with no restriction imposed, the ideal optimality condition is certainly that the sensitivities of all elements remain the same as implied by the Kuhn–Tucker necessary condition. However when the design variables are restricted to be either 0 or 1, the optimality criterion for a general minimization problem can be described as: the void element sensitivities are always higher than the solid element sensitivities which are less than or equal to zero. Therefore in order to reduce the maximum displacement, xi shall be changed from 1 to 0 for elements with highest displacement sensitivities and from 0 to 1 for elements with lowest sensitivities. Not that in a ‘‘hard-kill’’ approach, switching xi from 1 to 0 is equivalent to switching off the relevant solid elements, while switching xi from 0 to 1 is equivalent to re-admitting the relevant element. In the numerical implementation, a target volume is set for each iteration. A sensitivity threshold is determined according to the target volume such that elements with sensitivities higher than the threshold are removed and vice versa. In practice, this sensitivity
f = 1.0
when (d∗ − dmax )/d∗ ≥ tol,
(d − dmax )/d ∗
f =
or
(21)
∗
when (d∗ − dmax )/d∗ < tol
(22) tol where tol is a user-defined tolerance and chosen in the current work as 5%. Note that the gap function can be negative indicating that the displacement constraint is violated. When the displacement constraint is reached, the volume change stops while the element switching (solid or void) still continues which intends to reduce the maximum displacement according to the element sensitivities. Thus in the next iteration the gap function f may become positive and the total volume fraction continues to be reduced. This way, the proposed approach asymptotically pushes the maximum global displacement towards the allowable value while reducing the volume simultaneously. 4.4. Solution convergence Solution convergence is usually required for iterative methods to confirm the final solution has been achieved to terminate the iteration. For a volume minimization problem with a global displacement constraint, convergence shall be confirmed on both the volume and the maximum displacement. Solution convergence is checked by evaluating the relevant variation of the objective function and the constraints in the last few iterations. The convergence criteria for the volume and maximum displacement are presented in the following form.
error =
N k−i+1 k−N −i+1 (V − V ) i=1
N
≤ t∗
(23)
V k−i+1
i =1
error =
N k−i+1 k−N −i+1 (d − dmax ) max i=1
N
≤ t∗
(24)
k−i+1 dmax
i =1 k−i+1 where V k−i+1 and dmax denote the volume and the maximum displacement of the (k − i + 1)th iteration respectively, t ∗ is the convergence tolerance which can be set as small as 1 × 10−3 . N is an integral number chosen as 5 in the current work, i.e. the solution convergence is checked via the performance variation in the last 10 iterations. If the relative variation error is within the convergence tolerance t, the optimization procedure terminates and an optimum solution is obtained.
62
Z.H. Zuo, Y.M. Xie / Computer-Aided Design 56 (2014) 58–67
5. Numerical examples and discussions As the present method can be easily implemented as a postprocessor to commercial FEA software packages, we solve the structural static equilibrium with ABAQUS for all the examples hereafter. 5.1. A distributed-loaded cantilever
Fig. 1. Flowchart of the proposed method.
Fig. 2. A cantilever with non-designable deck: design domain.
4.5. Numerical procedure The proposed approach for volume minimization with a global displacement control can be summarized into the flowchart in Fig. 1.
The first example addresses the topology optimization of a cantilever shown in Fig. 2. A uniform distributed load is acting over the top non-designable deck illustrated in dark grey. This example aims to minimize the volume while the maximum absolute displacement is not allowed to exceed a limit of 1.48 × 10−6 m. The designable domain is discretized into 60 × 40 plane stress quadrilateral elements with a unity thickness. The material used has the property: E = 2.0 × 1011 Pa, v = 0.3. The BESO optimization performs with the evolutionary ratio ER = 2.0% and filter radius rmin = 1.0 m. The evolutionary history for the maximum absolute displacement and the volume fraction is plotted in Fig. 3. During the evolution, small jumps of the maximum displacement are identified which are caused by switching of the local controlled nodes. However due to the present stabilization strategies, the maximum displacement still converges immediately below the imposed value. The optimization run takes totally 45 iterations before the final solution convergence. Note that in each iteration the proposed approach requires two FEAs (one on the real system and the other on the virtual system) that consume a large portion of the computational effort. Compared with FEA, the computational cost spent on the optimization algorithms, e.g. calculating element sensitivities and removing elements, is practically trivial. Therefore, the actual computational time will highly depend on the efficiency of the FEA solver used. The final topology is shown in Fig. 4. A minimized volume fraction of 49.8% is obtained while the maximum displacement of 1.478 × 10−6 m satisfies the absolute displacement constraint of 1.48 × 10−6 m. Note that the contour in the deformed view represents the displacement magnitude. A comparison is made with the final optimal stiffness design shown in Fig. 5. This rival is produced by a conventional stiffness optimization with the volume constraint of v = 49.8% (to compete with the displacement design above), for the same model. The difference between these two final topologies is obvious by surveying the deformed shapes. The displacement design has a strong bar supporting the end of the non-designable deck while no bar is left at the middle, in order to reduce the deflection of the deck end point which is the critical position for the maximum displacement
Fig. 3. Evolutionary history of the cantilever with displacement constraint.
Z.H. Zuo, Y.M. Xie / Computer-Aided Design 56 (2014) 58–67
(a) Undeformed view.
63
(b) Deformed view (scaled).
Fig. 4. Final topology of the cantilever with displacement constraint (V = 49.8%, dmax = 1.478 × 10−6 m).
(a) Undeformed view.
(b) Deformed view (scaled).
Fig. 5. Comparison for the cantilever from conventional stiffness optimal design (V = 49.8%, dmax = 1.834 × 10−6 m).
Fig. 6. A simply-supported bridge: design domain.
Fig. 7. Evolutionary histories for the simply-supported bridge.
to occur. On the other hand, as the stiffness design considers the structural overall stiffness, the supporting bars are positioned relatively evenly under the non-designable deck, complying with the distributed loading condition. The two final solutions have a significant difference in the maximum displacement: with the same volume fraction, the maximum displacement of the stiffness design is 24% higher than that of the displacement design.
5.2. A simply-supported bridge The second example deals with a simply-supported bridge shown in Fig. 6. A uniform distributed load is acting on a nondesignable deck layer at the bottom. The deck is simulated by beam elements with the height of 0.5 m and has the same thickness and material properties as the bridge designable parts. For a bridge, it is the deflection of the deck that is the crucial functional factor.
64
Z.H. Zuo, Y.M. Xie / Computer-Aided Design 56 (2014) 58–67
(a) Undeformed view.
(b) Deformed view (scaled). Fig. 8. Final topology of the bridge with displacement constraint (V = 57.8%, dmax = 6.495 × 10−5 m).
(a) Undeformed view.
(b) Undeformed view (scaled). Fig. 9. Comparison for the bridge from conventional stiffness optimal design (V = 57.8%, dmax = 7.291 × 10−5 m).
a 0
b
c
d
Fig. 10. Optimal designs with different global displacement constraint: (a) d∗ = 5.8 × 10−5 m; (b) d∗ = 6.3 × 10−5 m; (c) d∗ = 6.5 × 10−5 m; (d) d∗ = 6.7 × 10−5 m.
Therefore for the optimal design, the volume is to be minimized with the deck deflection in the vertical direction (Y -direction) within a limit of 6.50 × 10−5 m. The design domain is discretized into 200 × 40 plane stress quadrilateral elements with the elastic properties of: E = 2.0 × 1011 Pa and v = 0.3. An evolutionary ratio ER = 1.0% and filter radius rmin = 1.0 m are used in the BESO optimization. The evolutionary histories for the maximum vertical deflection and the volume fraction are plotted in Fig. 7. An obvious jump of the maximum vertical deflection is found during the evolution. Accordingly, the volume fraction is increased automatically in order to fix the violation of the displacement constraint. After the maximum deflection gets again within the allowable limit, the volume fraction continues to decrease, until finally it converges to the minimum of 57.8%. At the same time, the maximum deflection satisfies the convergence criterion and converges immediately below the imposed displacement limit. The final topology in the undeformed and deformed views is shown in Fig. 8. Note that the contour in the deformed view depicts
the displacement in the Y -direction. The final maximum vertical deflection is 6.495 × 10−5 m which satisfies the constraint of 6.50 × 10−5 m. Stiffness optimization for this simply-supported bridge is carried out as a comparison. The final optimal topology is shown in Fig. 9. The stiffness design is produced by conventional stiffness optimization under a volume constraint. With the same volume fraction (57.8%), the maximum vertical deflection of the stiffness design (7.291 × 10−5 m) is 12.3% higher than that of the deflection design (6.495 × 10−5 m). Comparing the two deformed optimal designs, it is found that the deformation of the middle part of the deck in the deflection design is flat which reduces the maximum vertical deflection and outperforms the stiffness design in this regard. It is seen that the volume variation is controlled adaptively with the constrained displacement and the element update scheme works well to achieve an optimum with minimized volume fraction which outperforms conventional stiffness optimal designs.
Z.H. Zuo, Y.M. Xie / Computer-Aided Design 56 (2014) 58–67
65
a
b
c
Fig. 11. A 3D spoiler design domain from different views: (a) top; (b) side; (c) perspective.
Perspective-1
Perspective-2
Fig. 12. Smoothed final topology of the spoiler core (V = 20.9%, dmax = 1.495 × 10−3 cm).
In order to examine the effect of the prescribed limit of the global displacement on the final design, additional imposed values are tried on the displacement constraint: 5.80 × 10−5 m, 6.30 × 10−5 m and 6.70 × 10−5 m. All the four solutions are summarized in Fig. 10. It is seen that generally a smaller volume can be achieved with a more relaxed control on the global displacement. 5.3. A three dimensional spoiler As the proposed BESO acts as a post-processor to commercial FEA packages, it is easy to extend to three dimensions. This example deals with topology optimization of a generic 3D spoiler for a typical aircraft wing. Fig. 11 shows the spoiler which is covered by a non-designable skin of 0.1 cm thickness, a uniform-distributed load is acting on the bottom surface of the skin. In this example, the volume of the spoiler core is to be minimized. A limitation on the vertical deflection (in Y -direction) d∗ = 1.50 × 10−3 cm for the skin is applied in order to maintain the aerodynamic shape of the spoiler. Since the spoiler model is symmetric due to the middle X –Y plane, only a half of the model is optimized to save
the computational cost. The skin is simulated with 11 620 fournode shell elements and the core is discretized into 92 664 trilinear cube elements. The material properties for the structure are: E = 2.0 × 1011 Pa and v = 0.3. The optimization performs with an evolutionary ratio ER = 2.0% and filter radius rmin = 2.5 cm. The final topology is obtained with a resultant volume fraction 20.9% and the maximum vertical deflection 1.495 × 10−3 cm satisfying the constraint. The final design of the spoiler core is shown in Fig. 12 to clearly exhibit the optimized spoiler after being smoothed. Note that a technique is used here to smooth the final model surface by slightly moving the surface nodes in order to produce a design-friendly topology. The final full spoiler is shown in Fig. 13 with the skin attached. 6. Conclusions This paper addresses a topology optimization problem of volume minimization with a global displacement control. A global displacement constraint is a difficult but important constraint in practical structural design where the displacements of some parts
66
Z.H. Zuo, Y.M. Xie / Computer-Aided Design 56 (2014) 58–67
Perspective-1
Perspective-2
Fig. 13. Final full spoiler model with skin.
or of the whole structure are not allowed to exceed a certain limit. In the proposed approach, the global displacement control is realized by confining the maximum local displacements that may occur at different locations in the structure during the optimization history. In order to overcome the numerical instabilities caused by the multiple local maximum displacements at unpredictable locations, stabilization algorithms are proposed: priority-based control of multiple nodal displacements, adaptive history averaging of element sensitivities, and adaptive volume variation based on the constrained displacement. The hard-kill manner that gradually reduces the finite element model size saves the computational cost for the displacement constraint problem that needs to analyse the additional virtual system. Numerical tests through various 2D and 3D examples are carried out to examine the effectiveness and efficiency of the proposed BESO approach. The evolutionary history for the maximum global deflection is found not always smooth, local jumps of the performance can occur within adjacent iterations. This is caused by switching of critical nodes with the maximum deflection from iteration to iteration. With the adaptive volume variation algorithm, volume increase can be observed as the proposed approach tries to amend automatically. On the other hand, multiple control of critical nodes and history-averaging sensitivities with weighted factors further reduce the displacement fluctuation. As a result, the evolutionary history becomes smooth and stable to reach a final convergent optimum. The global displacement control has very important implementations in practical structural design. With the aim for economic light-weight design, the global displacement control technique proposed in this paper is especially valuable for the design of structures whose deformed exterior shapes are a crucial factor, such as an aircraft wing or a submarine body. References [1] Prager W, Rozvany GIN. Optimization of structural geometry. In: Dynamical systems. New York: Academic Press; 1977. p. 265–93. [2] Bendsøe MP. Optimal shape design as a material distribution problem. Struct Optim 1989;1:193–202. [3] Rozvany GIN, Zhou M, Birker T. Generalized shape optimization without homogenization. Struct Optim 1992;4:250–4. [4] Svanberg K. The method of moving asymptotes—a new method for structural optimization. Internat J Numer Methods Engrg 1987;24:359–73. [5] Xie YM, Steven GP. A simple evolutionary procedure for structural optimization. Comput Struct 1993;49:885–96.
[6] Xie YM, Steven GP. Evolutionary structural optimization. London: Springer; 1997. [7] Querin OM, Steven GP, Xie YM. Evolutionary structural optimization (ESO) using a bi-directional algorithm. Eng Comput 1998;15(8):1031–48. [8] Yang XY, Xie YM, Steven GP, Querin OM. Bidirectional evolutionary method for stiffness optimization. AIAA J 1999;37(11):1483–8. [9] Manicharajah D, Xie YM. An evolutionary method for optimization of plate buckling resistance. Finite Elem Anal Des 1998;29:205–30. [10] Li Q, Steven GP, Xie YM. Shape and topology design for heat conduction by evolutionary structural optimization. Int J Heat Mass Transfer 1999;42: 3361–71. [11] Querin OM, Steven GP, Xie YM. Evolutionary structural optimization using an additive algorithm. Finite Elem Anal Des 2000;34:291–308. [12] Huang X, Zuo ZH, Xie YM. Evolutionary topological optimization of vibrating continuum structures for natural frequencies. Comput Struct 2010;88:357–64. [13] Huang X, Xie YM. Evolutionary topology optimization of continuum structures: methods and applications. Chichester: John Wiley & Sons; 2010. [14] Bendsøe MP, Sigmund O. Topology optimization: theory, methods and applications. Berlin (Heidelberg): Springer; 2003. [15] Yin L, Yang W. Optimality criteria method for topology optimization under multiple constraints. Comput Struct 2001;79:1839–50. [16] Huang X, Xie YM. Evolutionary topology optimization of continuum structures with an additional displacement constraint. Struct Multidiscip Optim 2010;40: 409–16. [17] Zuo ZH, Xie YM, Huang X. Evolutionary topology optimization of structures with multiple displacement and frequency constraints. Adv Struct Eng 2012; 15(2):385–98. [18] Sobieszczanski-Sobieski J, Haftka RT. Multidisciplinary aerospace design optimization: survey of recent developments. Struct Optim 1997;14:1–23. [19] Liang QQ, Xie YM, Steven GP. A performance index for topology and shape optimization of plate bending problems with displacement constraints. Struct Multidiscip Optim 2001;21:393–9. [20] Rong JH, Yi JH. A structural topological optimization method for multidisplacement constraints and any initial topology configuration. Acta Mech Sin 2010;26:735–44. [21] Suresh K, Ramani A, Kaushik A. An adaptive weighting strategy for multiload topology optimization. In: Proceedings of the ASME 2012 international design engineering technical conferences & computers and information in engineering conference. [22] Farahpour I. Constrained optimization of structures with displacement constraints under various loading conditions. Adv Eng Softw 2010;41:580–9. [23] Sonmez M. Discrete optimum design of truss structures using artificial bee colony algorithm. Struct Multidiscip Optim 2011;43:85–97. [24] Bendsøe MP, Sigmund O. Material interpolation schemes in topology optimization. Arch Appl Mech 1999;69:635–54. [25] Sigmund O, Peterson J. Numerical instabilities in topology optimization: a survey on procedures dealing with checkerboards, mesh-dependencies and local minima. Struct Optim 1998;16:68–75. [26] Huang X, Xie YM. Convergent and mesh-independent solutions for the bidirectional evolutionary structural optimization method. Finite Elem Anal Des 2007;43:1039–49. [27] Taylor JE, Bensøe MP. A mutual energy formulation for optimal structural design. Struct Multidiscip Optim 2001;22:95–101.
Z.H. Zuo, Y.M. Xie / Computer-Aided Design 56 (2014) 58–67 Zhi Hao Zuo obtained his bachelor degree in Civil Engineering from Tongji University in China. After he completed his Ph.D. in topology optimization of periodic structures in 2009 in RMIT University in Australia, he has been working as a research fellow in the Centre for Innovative Structures and Materials, at RMIT University.
67
Yi Min Xie obtained his bachelor degree in Engineering Mechanics from Shanghai Jiao Tong University in China and Ph.D. in Civil Engineering from Swansea University in Britain. He is one of the pioneers of the evolutionary structural optimization (ESO) method. He is an elected fellow of the Australian Academy of Technological Sciences and Engineering.
Contents lists available at ScienceDirect
Computer-Aided Design journal homepage: www.elsevier.com/locate/cad
Evolutionary topology optimization of continuum structures with a global displacement control✩ Zhi Hao Zuo, Yi Min Xie ∗ Centre for Innovative Structures and Materials, School of Civil, Environmental and Chemical Engineering, RMIT University, GPO Box 2476, Melbourne 3001, Australia
highlights • • • •
The proposed method achieves optimal designs with globally controllable deflections. Solutions with the proposed method can maintain aerodynamic shape when deformed. Economic lightweight designs are achieved with a volume minimization scheme. Robust and effective numerical algorithms are proposed for clear 0–1 designs.
article
info
Article history: Received 17 August 2013 Accepted 14 June 2014 Keywords: Topology optimization Displacement control Bi-directional evolutionary structural optimization (BESO) Aircraft design Aerodynamics
abstract The conventional compliance minimization of load-carrying structures does not directly deal with displacements that are of practical importance. In this paper, a global displacement control is realized through topology optimization with a global constraint that sets a displacement limit on the whole structure or certain sub-domains. A volume minimization problem is solved by an extended evolutionary topology optimization approach. The local displacement sensitivities are derived following a powerlaw penalization material model. The global control of displacement is realized through multiple local displacement constraints on dynamically located critical nodes. Algorithms are proposed to secure the stability and convergence of the optimization process. Through numerical examples and by comparing with conventional stiffness designs, it is demonstrated that the proposed approach is capable of effectively finding optimal solutions which satisfy the global displacement control. Such solutions are of particular importance for structural designs whose deformed shapes must comply with functioning requirements such as aerodynamic performances. © 2014 Elsevier Ltd. All rights reserved.
1. Introduction Topology optimization is a powerful tool resulting in high structural performance with great material saving. It is capable of speeding up the structural design process and producing reliable solutions for structural design problems. This area has been extensively investigated in the past three decades since the modern formulation of optimal layout theories by Prager and Rozvany [1]. Several Finite Element (FE) based methods have been developed for topology optimization of continuum structures. One most popular technique is the Solid Isotropic Material with Penalization (SIMP) method [2,3], where the isotropic material property in each
✩ This paper has been recommended for acceptance by Michael Yu Wang.
∗
Corresponding author. Tel.: +61 3 99253655; fax: +61 3 96390138. E-mail address: [email protected] (Y.M. Xie).
http://dx.doi.org/10.1016/j.cad.2014.06.007 0010-4485/© 2014 Elsevier Ltd. All rights reserved.
element is determined by the element relative density. By penalizing the element relative density using a power-law interpolation scheme, a nearly void–solid design is expected after an iterative procedure such as the Method of Moving Asymptotes [4]. A popular group of methods for topology optimization is the Evolutionary Structural Optimization (ESO) method and its descendent versions. The ESO method was first proposed by Xie and Steven [5,6] in the early 1990s. This original version follows a straightforward algorithm of iteratively removing inefficient material and thus drives the structure to evolve towards an optimum. The bi-directional evolutionary structural optimization (BESO) method was proposed as an extended version of ESO in the late 1990s [7,8]. The BESO method allows material to be added to the most demanding places while inefficient material is removed simultaneously. The ESO/BESO methods use binary design variables and directly deliver void–solid solutions that are very desirable. The ESO/BESO methods can be easily implemented in a wide range of engineering applications such as [9–13].
Z.H. Zuo, Y.M. Xie / Computer-Aided Design 56 (2014) 58–67
Deflection constrained optimization problems are related with stiffness optimization problems [14] to a large extent. These two sets of problems are in fact equivalent under single point loaded systems and with deflection control at the load point in the load direction. However generally, the stiffness optimal design is not the same with the optimal deflection design. The stiffness optimization increases the structural overall stiffness and thus the deflection is roughly reduced as an indirect effect. Very often in practice, it is not the stiffness but the deflection over parts of or whole of the structure that is the crucial factor. Despite vast research on stiffness optimization, the local displacement constraint on continuum structural optimization is found in relatively less literature such as [15–17]. For some structures such as an aircraft wing, the exterior surface should undergo minimal shape change under deformation in order to maintain the aerodynamic performance [18]; this requires that the displacements of the surface shall remain within a certain limit. In such cases, the displacements of a group of local nodes are of concern; or more generally, the displacement limit is addressed as a global constraint. The deflection constrained optimization problem addressed in the present paper is a global control problem. The common way of obtaining the displacement sensitivity is to apply a unit virtual load on the original model and get the displacement vector from the virtual system [16]. For global displacement control, the maximum deflections occur unpredictably at different nodes through the optimization history. Therefore the virtual loads must be built dynamically after the maximum deflections are identified in the real system at each stage; and the real and virtual systems must be analysed in two subsequent FEAs. This increases the computational cost significantly. More importantly, numerical instabilities caused by dynamical relocation of control nodes are hard to handle; as a consequence, stable solution convergence is highly difficult to achieve. Multiple displacement constraints optimization can be found in literature such as [19–23]. The common strategy of this art is to control the local displacement at pre-defined locations, usually at the load point and in the load direction, such as in [19–21]; in other words, a modified compliance minimization problem is actually dealt with. Pre-defined displacement constraints are also applied to specific structural systems such as truss structures [22,23]. The complexity of solving multiple constraints is reduced due to the fact that less design variables are present in the design domain. Since the displacement control locations and directions are limited (by pre-definition or by the applied load), a genuine ‘‘global control’’ for displacements is not realized in these methods. In some cases, the critical displacement may occur at a location completely different from the load points; such fixed multi-constraint approaches will not be able to tackle the maximum displacement. Multi-constraint optimization for continuum structures is usually a highly non-linear problem that is theoretically very difficult to solve, especially in cases as the global displacement constraint that controls displacements on dynamic locations. This paper presents a global control method for displacements of continuum structures. A topology optimization problem of volume minimization is formulated and solved by a new BESO approach. The global displacement constraint directly imposes a maximum allowable limit for nodal displacements within the whole domain or some user-specified sub-domains of the structure. The locations and directions of maximum nodal displacements are dynamically detected, the numerical instabilities of which are adaptively dealt with by robust stabilization algorithms. The subsequent sections are organized as follows: Section 2 formulates the topology optimization problem for global displacement control; Section 3 presents the sensitivity calculation for displacements; Section 4 is dedicated to the global displacement constraint algorithms for stabilizing the evolution history and solution convergence; the numerical implementation of the proposed method is outlined in this section; Section 5 shows numerical examples with discussions; the conclusions are drawn in Section 6.
59
2. Problem statement The global displacement control can be realized by confining the maximum displacement of the structure. With the volume being the objective function, the problem is formulated as follows. Minimize V , subject to V =
N
Vi x i ,
xi ∈ {xmin , 1}
(1)
i =1
dj ≤ dmax ≤ d∗ ,
j = 1, . . . , L
P = KU
(2) (3)
where V is the volume of the structure, Vi is the volume of the ith element, N is the total number of elements; xi is the binary design variable, i.e. either 1 or 0 denoting the element status of ‘‘solid’’ or ‘‘void’’; the design variable can also be treated as the element relative density [16] with a very small value such as 10−6 representing the void status; d∗ is the maximum allowable displacement value, dmax is the maximum displacement in the structure under loading, and dj is the displacement of the jth element within a prescribed domain Ωdisp where the displacement constraint is active. Eq. (3) is the static equilibrium with P as the applied load vector, K as the global stiffness matrix and U as the global displacement vector. The displacements above can be addressed in any certain directions, or can be treated as relative displacements simply when some nodes in Ωdisp are fixed. Constraining the relative displacements is of particular importance when the deformed shape is concerned. 3. Displacement sensitivity calculation 3.1. Material interpolation scheme Material interpolation schemes usually express the material properties as a function of the design variables in order to facilitate the sensitivity analysis. In the power-law material model SIMP [24], the Young’s modulus of an element is determined by the element relative density through the following penalization formulation. p
E (xi ) = xi E 0
(4)
0
where E is the Young’s modulus of the solid material and p denotes the penalty exponent which is usually set as p ≥ 3. Note that the Poisson’s ratio is usually assumed to be irrelevant in SIMP. Therefore, the stiffness matrix can be expressed in a similar way as follows. K=
p
xi K0i
(5)
i
where K0i denotes the stiffness matrix for the solid element, i.e. when xi = 1. 3.2. Displacement sensitivity in axial directions The displacement kth component can be obtained by multiplying the displacement vector with a unit virtual load vector Fk , of which the kth component is unity while all other components are zero. uk = Fk · u
(6)
F = {0, 0, . . . , 1, 0, . . .} = {f1 , f2 , . . . , fk , fk+1 , . . .}. k
k
(7)
With the virtual load F and the applied load P being constant, differentiating the kth displacement component with respect to
60
Z.H. Zuo, Y.M. Xie / Computer-Aided Design 56 (2014) 58–67
the ith element and substituting Eq. (3) gives the following k,T
The sensitivity of the absolute displacement is then immediately derived as
k,T
∂u ∂u ∂F ∂F = u + Fk,T = u ∂ xi ∂ xi ∂ xi ∂ xi −1 ∂P ∂K P + K−1 + Fk,T ∂ xi ∂ xi k
p−1
αi = −
u21 + u22 + u23 p−1
∂ K −1 ∂K ∂K = Fk,T P = −Fk,T K−1 u = −uk,T u ∂ xi ∂ xi ∂ xi
(8)
where uk = K−1 Fk is system response (displacement field) under the unit virtual load. Substituting the material interpolation in Eq. (5) into the above derivative gives finally the element sensitivity for the kth displacement vector with respect to the ith element.
αi =
pxi
∂u = −pxpi −1 uik,T K0i ui ∂ xi k
(9)
k,T
where ui and ui are the element displacement vector under the unit virtual load and the real load conditions respectively. By setting the penalty exponent p as infinity, assigning an element with the minimum relative density is equivalent to removing that element (E = 0), and the sensitivity is simply zero. In such a ‘‘hard-kill’’ manner, the model size is reduced and the computation cost can be saved significantly.
= −
pxi u21
+
u22
1 ,T
u1 ui
+ u2 u2i ,T + u3 u3i ,T K0i ui
k,T
+
u23
ui K0i ui .
(14)
3.4. Filtering sensitivities for mesh-independence Topology optimization usually encounters chequerboard and mesh-dependence problems [25]. To circumvent these problems, a mesh-independency filter [26] is used by averaging the element sensitivities with its neighbourhood elements using a low-pass filter of radius rmin . Firstly, nodal sensitivity numbers αjn which do not carry any physical meaning on their own are calculated by evenly extrapolating the sensitivities of connected elements; then these nodal sensitivity numbers are converted back into elements using a weight function w based on rij —the distance between the centre of the element i and node j, and on the filter radius rmin , in the following.
w(rij )αjn , αˆ i = w(rij )
3.3. Displacement sensitivity in arbitrary directions The displacement sensitivity in non-axial directions can be derived based on the displacement sensitivity in axial directions. The real displacement at the kth node is expressed by a vector u = {u1 , u2 , u3 } with u1 , u2 and u3 being the components in the x, y and z axial directions. Then the projective displacement on a specified direction e = {e1 , e2 , e3 } can be treated as the function of the three components. u˜ = u · e = u˜ (u1 , u2 , u3 ) = u1 e1 + u2 e2 + u3 e3 .
(10)
The derivative of u˜ is obtained as the dot product of the gradient of u˜ and the derivative of u.
∂ u˜ ∂ u1 ∂ u˜
∂ u1 ∂ xi ∂u
∂ u˜ ∂u 2 = ∇ u˜ · = · . ∂ xi ∂ xi ∂ x ∂ u i 2 ∂ u3 ∂ u˜ ∂ xi ∂ u3
(11)
Substituting Eq. (9) into the above equation, we finally obtain the displacement sensitivity in the specified direction e defined by {e1 , e2 , e3 }.
u1i ,T K0i ui e 1 ∂ u˜ = −pxpi −1 e2 · u2i ,T K0i ui αi = ∂ xi 3,T 0 e3 ui Ki ui = −pxpi −1 e1 ui1,T + e2 u2i ,T + e3 u3i ,T K0i ui (12)
e=
u21 + u22 + u23
,
u2 u21 + u22 + u23
,
The local displacement sensitivity in Section 3 forms the basis for global displacement control. To realize such a global control, straightforwardly, all the displacements within the concerned domain Ωdisp shall be addressed, i.e. the global displacement constraint in Eq. (2) shall be satisfied. 4.1. From local control to global control Ideally, the global displacement constraint is satisfied if the globally maximum displacement dmax is within the allowable limit d∗ . In other words, the global displacement constraint might be simplified to constraining the maximum displacement d∗ over the set Ωdisp . In such a ‘‘worst case’’ algorithm, the location of the maximum displacement might vary from iteration to iteration. As a consequence, the displacement of other nodes might increase significantly to become the maximum displacement. The jumping control node will lead to numerical instabilities and make the solution difficult to converge. Therefore, compared to controlling only one node of the maximum displacement, it is more reasonable to control a set of nodes whose displacements are close to dmax . when
dmax − dn dmax
≤t (16)
Note that according to the linear superposition of displacements, uki is the system response under the virtual load Fk = {e1 , e2 , e3 }. If the unit vector e is in the direction of the displacement u, namely u1
(15)
4. Global displacement control
dn ≤ dn+1 ≤ · · · ≤ dn+Q −2 ≤ dmax ≤ d∗ ,
= −pxpi −1 uik,T K0i ui .
where w = max(rmin − rij , 0).
u3
u21 + u22 + u23
.
(13)
where Q is the number of nodes to be controlled; t is a pre-defined tolerance that is set to a relatively small value such as 5% used in the current work. Increasing the tolerance t means putting more nodes into control which increases the smoothness of the evolution history. However the final optimum might degrade since the maximum displacement is not of enough focus. To solve this issue, the individual virtual loads are calculated using a weighted function that reflects their priorities based on the relevant displacements.
Fj = qj |1| = dj . dmax
(17)
Z.H. Zuo, Y.M. Xie / Computer-Aided Design 56 (2014) 58–67
Note that this mechanism is equivalent to formulating a global objective function for displacement [27] which is composed by local displacement functions with weighted factors qj =
dj dmax
.
61
threshold can be easily determined using a simple and efficient bisection algorithm. Starting from the full design of the structure, the target volume of each iteration is progressively reduced as described below. V k = V k−1 (1 − ER · f )
4.2. Numerical stabilization
(20)
where V is the target volume of the current iteration and V k−1 is the volume of the previous iteration; ER is the user-prescribed evolutionary volume ratio which determines the volume variation between two adjacent iterations, e.g. ER = 2% means the remaining volume will be changed by 2% (if f = 1.0) for each iteration; note that with a greater ER the volume will be reduced faster, but at the same time numerical instabilities may also arise that may require more iterations for the solution to converge; usually ER is selected between 1% and 5%; f is a function indicating the relative gap between the current maximum displacement and the imposed constraint value d∗ . f serves as a factor to adaptively fine-tune the volume change step size in order to stabilize the optimization objective and constraints. k
Still, the nodes under control may still be subject to change through the optimization history. The sensitivity history-averaging [26] can be applied to further stabilize the evolution. However, a refined modification shall be made here to achieve more reasonable smoothing. Instead of simple averaging [26], a weighted combination of sensitivities of two adjacent iterations is introduced.
α = w1 α k + w2 α k−1 where w1 + w2 = 1.
(18)
Here α is the completed element sensitivity after smoothing, α k and α k−1 are the sensitivities of the current and the immediately previous iteration, w1 and w2 are the weighted factors for both items. During the optimization procedure, the locations of the maximum displacements, i.e. the control nodes of two adjacent iterations might vary significantly. Thus the evolutionary history may suffer from numerical instabilities caused by the sudden change of sensitivities. The usage of the sensitivity weighted factors is able to stabilize the optimization procedure by setting w1 and w2 adaptively. As simple sensitivity averaging [26], the initial values can be set as w1 = w2 = 0.5; when the maximum displacement locations changes significantly, it is helpful to set w2 greater than w1 in order to avoid excessive disturbance of the optimization trend. The weighted factors w2 can take the following general form:
w2 = 0.5 +
q Q
b
(19)
where q (q < Q ) is the number of new nodes in the total Q control nodes within Ωdisp , and b is a user-prescribed gap factor chosen between 0 and 0.5. Note that the greater b is chosen, the less the solution will suffer from the numerical instabilities, but it might take more iterations for the solution to converge to the final optimum; 0.3 is an appropriate candidate for b. 4.3. Element update and numerical stabilization The element sensitivity indicates the changing trend of the maximum global displacement due to changing the design variable of one element. Changing the design variable xi from 0 to 1 for one element with negative sensitivity decreases the maximum global displacement, while changing xi from 1 to 0 of one element with positive sensitivity decreases the maximum global displacement. For optimization methods using continuous design variables with no restriction imposed, the ideal optimality condition is certainly that the sensitivities of all elements remain the same as implied by the Kuhn–Tucker necessary condition. However when the design variables are restricted to be either 0 or 1, the optimality criterion for a general minimization problem can be described as: the void element sensitivities are always higher than the solid element sensitivities which are less than or equal to zero. Therefore in order to reduce the maximum displacement, xi shall be changed from 1 to 0 for elements with highest displacement sensitivities and from 0 to 1 for elements with lowest sensitivities. Not that in a ‘‘hard-kill’’ approach, switching xi from 1 to 0 is equivalent to switching off the relevant solid elements, while switching xi from 0 to 1 is equivalent to re-admitting the relevant element. In the numerical implementation, a target volume is set for each iteration. A sensitivity threshold is determined according to the target volume such that elements with sensitivities higher than the threshold are removed and vice versa. In practice, this sensitivity
f = 1.0
when (d∗ − dmax )/d∗ ≥ tol,
(d − dmax )/d ∗
f =
or
(21)
∗
when (d∗ − dmax )/d∗ < tol
(22) tol where tol is a user-defined tolerance and chosen in the current work as 5%. Note that the gap function can be negative indicating that the displacement constraint is violated. When the displacement constraint is reached, the volume change stops while the element switching (solid or void) still continues which intends to reduce the maximum displacement according to the element sensitivities. Thus in the next iteration the gap function f may become positive and the total volume fraction continues to be reduced. This way, the proposed approach asymptotically pushes the maximum global displacement towards the allowable value while reducing the volume simultaneously. 4.4. Solution convergence Solution convergence is usually required for iterative methods to confirm the final solution has been achieved to terminate the iteration. For a volume minimization problem with a global displacement constraint, convergence shall be confirmed on both the volume and the maximum displacement. Solution convergence is checked by evaluating the relevant variation of the objective function and the constraints in the last few iterations. The convergence criteria for the volume and maximum displacement are presented in the following form.
error =
N k−i+1 k−N −i+1 (V − V ) i=1
N
≤ t∗
(23)
V k−i+1
i =1
error =
N k−i+1 k−N −i+1 (d − dmax ) max i=1
N
≤ t∗
(24)
k−i+1 dmax
i =1 k−i+1 where V k−i+1 and dmax denote the volume and the maximum displacement of the (k − i + 1)th iteration respectively, t ∗ is the convergence tolerance which can be set as small as 1 × 10−3 . N is an integral number chosen as 5 in the current work, i.e. the solution convergence is checked via the performance variation in the last 10 iterations. If the relative variation error is within the convergence tolerance t, the optimization procedure terminates and an optimum solution is obtained.
62
Z.H. Zuo, Y.M. Xie / Computer-Aided Design 56 (2014) 58–67
5. Numerical examples and discussions As the present method can be easily implemented as a postprocessor to commercial FEA software packages, we solve the structural static equilibrium with ABAQUS for all the examples hereafter. 5.1. A distributed-loaded cantilever
Fig. 1. Flowchart of the proposed method.
Fig. 2. A cantilever with non-designable deck: design domain.
4.5. Numerical procedure The proposed approach for volume minimization with a global displacement control can be summarized into the flowchart in Fig. 1.
The first example addresses the topology optimization of a cantilever shown in Fig. 2. A uniform distributed load is acting over the top non-designable deck illustrated in dark grey. This example aims to minimize the volume while the maximum absolute displacement is not allowed to exceed a limit of 1.48 × 10−6 m. The designable domain is discretized into 60 × 40 plane stress quadrilateral elements with a unity thickness. The material used has the property: E = 2.0 × 1011 Pa, v = 0.3. The BESO optimization performs with the evolutionary ratio ER = 2.0% and filter radius rmin = 1.0 m. The evolutionary history for the maximum absolute displacement and the volume fraction is plotted in Fig. 3. During the evolution, small jumps of the maximum displacement are identified which are caused by switching of the local controlled nodes. However due to the present stabilization strategies, the maximum displacement still converges immediately below the imposed value. The optimization run takes totally 45 iterations before the final solution convergence. Note that in each iteration the proposed approach requires two FEAs (one on the real system and the other on the virtual system) that consume a large portion of the computational effort. Compared with FEA, the computational cost spent on the optimization algorithms, e.g. calculating element sensitivities and removing elements, is practically trivial. Therefore, the actual computational time will highly depend on the efficiency of the FEA solver used. The final topology is shown in Fig. 4. A minimized volume fraction of 49.8% is obtained while the maximum displacement of 1.478 × 10−6 m satisfies the absolute displacement constraint of 1.48 × 10−6 m. Note that the contour in the deformed view represents the displacement magnitude. A comparison is made with the final optimal stiffness design shown in Fig. 5. This rival is produced by a conventional stiffness optimization with the volume constraint of v = 49.8% (to compete with the displacement design above), for the same model. The difference between these two final topologies is obvious by surveying the deformed shapes. The displacement design has a strong bar supporting the end of the non-designable deck while no bar is left at the middle, in order to reduce the deflection of the deck end point which is the critical position for the maximum displacement
Fig. 3. Evolutionary history of the cantilever with displacement constraint.
Z.H. Zuo, Y.M. Xie / Computer-Aided Design 56 (2014) 58–67
(a) Undeformed view.
63
(b) Deformed view (scaled).
Fig. 4. Final topology of the cantilever with displacement constraint (V = 49.8%, dmax = 1.478 × 10−6 m).
(a) Undeformed view.
(b) Deformed view (scaled).
Fig. 5. Comparison for the cantilever from conventional stiffness optimal design (V = 49.8%, dmax = 1.834 × 10−6 m).
Fig. 6. A simply-supported bridge: design domain.
Fig. 7. Evolutionary histories for the simply-supported bridge.
to occur. On the other hand, as the stiffness design considers the structural overall stiffness, the supporting bars are positioned relatively evenly under the non-designable deck, complying with the distributed loading condition. The two final solutions have a significant difference in the maximum displacement: with the same volume fraction, the maximum displacement of the stiffness design is 24% higher than that of the displacement design.
5.2. A simply-supported bridge The second example deals with a simply-supported bridge shown in Fig. 6. A uniform distributed load is acting on a nondesignable deck layer at the bottom. The deck is simulated by beam elements with the height of 0.5 m and has the same thickness and material properties as the bridge designable parts. For a bridge, it is the deflection of the deck that is the crucial functional factor.
64
Z.H. Zuo, Y.M. Xie / Computer-Aided Design 56 (2014) 58–67
(a) Undeformed view.
(b) Deformed view (scaled). Fig. 8. Final topology of the bridge with displacement constraint (V = 57.8%, dmax = 6.495 × 10−5 m).
(a) Undeformed view.
(b) Undeformed view (scaled). Fig. 9. Comparison for the bridge from conventional stiffness optimal design (V = 57.8%, dmax = 7.291 × 10−5 m).
a 0
b
c
d
Fig. 10. Optimal designs with different global displacement constraint: (a) d∗ = 5.8 × 10−5 m; (b) d∗ = 6.3 × 10−5 m; (c) d∗ = 6.5 × 10−5 m; (d) d∗ = 6.7 × 10−5 m.
Therefore for the optimal design, the volume is to be minimized with the deck deflection in the vertical direction (Y -direction) within a limit of 6.50 × 10−5 m. The design domain is discretized into 200 × 40 plane stress quadrilateral elements with the elastic properties of: E = 2.0 × 1011 Pa and v = 0.3. An evolutionary ratio ER = 1.0% and filter radius rmin = 1.0 m are used in the BESO optimization. The evolutionary histories for the maximum vertical deflection and the volume fraction are plotted in Fig. 7. An obvious jump of the maximum vertical deflection is found during the evolution. Accordingly, the volume fraction is increased automatically in order to fix the violation of the displacement constraint. After the maximum deflection gets again within the allowable limit, the volume fraction continues to decrease, until finally it converges to the minimum of 57.8%. At the same time, the maximum deflection satisfies the convergence criterion and converges immediately below the imposed displacement limit. The final topology in the undeformed and deformed views is shown in Fig. 8. Note that the contour in the deformed view depicts
the displacement in the Y -direction. The final maximum vertical deflection is 6.495 × 10−5 m which satisfies the constraint of 6.50 × 10−5 m. Stiffness optimization for this simply-supported bridge is carried out as a comparison. The final optimal topology is shown in Fig. 9. The stiffness design is produced by conventional stiffness optimization under a volume constraint. With the same volume fraction (57.8%), the maximum vertical deflection of the stiffness design (7.291 × 10−5 m) is 12.3% higher than that of the deflection design (6.495 × 10−5 m). Comparing the two deformed optimal designs, it is found that the deformation of the middle part of the deck in the deflection design is flat which reduces the maximum vertical deflection and outperforms the stiffness design in this regard. It is seen that the volume variation is controlled adaptively with the constrained displacement and the element update scheme works well to achieve an optimum with minimized volume fraction which outperforms conventional stiffness optimal designs.
Z.H. Zuo, Y.M. Xie / Computer-Aided Design 56 (2014) 58–67
65
a
b
c
Fig. 11. A 3D spoiler design domain from different views: (a) top; (b) side; (c) perspective.
Perspective-1
Perspective-2
Fig. 12. Smoothed final topology of the spoiler core (V = 20.9%, dmax = 1.495 × 10−3 cm).
In order to examine the effect of the prescribed limit of the global displacement on the final design, additional imposed values are tried on the displacement constraint: 5.80 × 10−5 m, 6.30 × 10−5 m and 6.70 × 10−5 m. All the four solutions are summarized in Fig. 10. It is seen that generally a smaller volume can be achieved with a more relaxed control on the global displacement. 5.3. A three dimensional spoiler As the proposed BESO acts as a post-processor to commercial FEA packages, it is easy to extend to three dimensions. This example deals with topology optimization of a generic 3D spoiler for a typical aircraft wing. Fig. 11 shows the spoiler which is covered by a non-designable skin of 0.1 cm thickness, a uniform-distributed load is acting on the bottom surface of the skin. In this example, the volume of the spoiler core is to be minimized. A limitation on the vertical deflection (in Y -direction) d∗ = 1.50 × 10−3 cm for the skin is applied in order to maintain the aerodynamic shape of the spoiler. Since the spoiler model is symmetric due to the middle X –Y plane, only a half of the model is optimized to save
the computational cost. The skin is simulated with 11 620 fournode shell elements and the core is discretized into 92 664 trilinear cube elements. The material properties for the structure are: E = 2.0 × 1011 Pa and v = 0.3. The optimization performs with an evolutionary ratio ER = 2.0% and filter radius rmin = 2.5 cm. The final topology is obtained with a resultant volume fraction 20.9% and the maximum vertical deflection 1.495 × 10−3 cm satisfying the constraint. The final design of the spoiler core is shown in Fig. 12 to clearly exhibit the optimized spoiler after being smoothed. Note that a technique is used here to smooth the final model surface by slightly moving the surface nodes in order to produce a design-friendly topology. The final full spoiler is shown in Fig. 13 with the skin attached. 6. Conclusions This paper addresses a topology optimization problem of volume minimization with a global displacement control. A global displacement constraint is a difficult but important constraint in practical structural design where the displacements of some parts
66
Z.H. Zuo, Y.M. Xie / Computer-Aided Design 56 (2014) 58–67
Perspective-1
Perspective-2
Fig. 13. Final full spoiler model with skin.
or of the whole structure are not allowed to exceed a certain limit. In the proposed approach, the global displacement control is realized by confining the maximum local displacements that may occur at different locations in the structure during the optimization history. In order to overcome the numerical instabilities caused by the multiple local maximum displacements at unpredictable locations, stabilization algorithms are proposed: priority-based control of multiple nodal displacements, adaptive history averaging of element sensitivities, and adaptive volume variation based on the constrained displacement. The hard-kill manner that gradually reduces the finite element model size saves the computational cost for the displacement constraint problem that needs to analyse the additional virtual system. Numerical tests through various 2D and 3D examples are carried out to examine the effectiveness and efficiency of the proposed BESO approach. The evolutionary history for the maximum global deflection is found not always smooth, local jumps of the performance can occur within adjacent iterations. This is caused by switching of critical nodes with the maximum deflection from iteration to iteration. With the adaptive volume variation algorithm, volume increase can be observed as the proposed approach tries to amend automatically. On the other hand, multiple control of critical nodes and history-averaging sensitivities with weighted factors further reduce the displacement fluctuation. As a result, the evolutionary history becomes smooth and stable to reach a final convergent optimum. The global displacement control has very important implementations in practical structural design. With the aim for economic light-weight design, the global displacement control technique proposed in this paper is especially valuable for the design of structures whose deformed exterior shapes are a crucial factor, such as an aircraft wing or a submarine body. References [1] Prager W, Rozvany GIN. Optimization of structural geometry. In: Dynamical systems. New York: Academic Press; 1977. p. 265–93. [2] Bendsøe MP. Optimal shape design as a material distribution problem. Struct Optim 1989;1:193–202. [3] Rozvany GIN, Zhou M, Birker T. Generalized shape optimization without homogenization. Struct Optim 1992;4:250–4. [4] Svanberg K. The method of moving asymptotes—a new method for structural optimization. Internat J Numer Methods Engrg 1987;24:359–73. [5] Xie YM, Steven GP. A simple evolutionary procedure for structural optimization. Comput Struct 1993;49:885–96.
[6] Xie YM, Steven GP. Evolutionary structural optimization. London: Springer; 1997. [7] Querin OM, Steven GP, Xie YM. Evolutionary structural optimization (ESO) using a bi-directional algorithm. Eng Comput 1998;15(8):1031–48. [8] Yang XY, Xie YM, Steven GP, Querin OM. Bidirectional evolutionary method for stiffness optimization. AIAA J 1999;37(11):1483–8. [9] Manicharajah D, Xie YM. An evolutionary method for optimization of plate buckling resistance. Finite Elem Anal Des 1998;29:205–30. [10] Li Q, Steven GP, Xie YM. Shape and topology design for heat conduction by evolutionary structural optimization. Int J Heat Mass Transfer 1999;42: 3361–71. [11] Querin OM, Steven GP, Xie YM. Evolutionary structural optimization using an additive algorithm. Finite Elem Anal Des 2000;34:291–308. [12] Huang X, Zuo ZH, Xie YM. Evolutionary topological optimization of vibrating continuum structures for natural frequencies. Comput Struct 2010;88:357–64. [13] Huang X, Xie YM. Evolutionary topology optimization of continuum structures: methods and applications. Chichester: John Wiley & Sons; 2010. [14] Bendsøe MP, Sigmund O. Topology optimization: theory, methods and applications. Berlin (Heidelberg): Springer; 2003. [15] Yin L, Yang W. Optimality criteria method for topology optimization under multiple constraints. Comput Struct 2001;79:1839–50. [16] Huang X, Xie YM. Evolutionary topology optimization of continuum structures with an additional displacement constraint. Struct Multidiscip Optim 2010;40: 409–16. [17] Zuo ZH, Xie YM, Huang X. Evolutionary topology optimization of structures with multiple displacement and frequency constraints. Adv Struct Eng 2012; 15(2):385–98. [18] Sobieszczanski-Sobieski J, Haftka RT. Multidisciplinary aerospace design optimization: survey of recent developments. Struct Optim 1997;14:1–23. [19] Liang QQ, Xie YM, Steven GP. A performance index for topology and shape optimization of plate bending problems with displacement constraints. Struct Multidiscip Optim 2001;21:393–9. [20] Rong JH, Yi JH. A structural topological optimization method for multidisplacement constraints and any initial topology configuration. Acta Mech Sin 2010;26:735–44. [21] Suresh K, Ramani A, Kaushik A. An adaptive weighting strategy for multiload topology optimization. In: Proceedings of the ASME 2012 international design engineering technical conferences & computers and information in engineering conference. [22] Farahpour I. Constrained optimization of structures with displacement constraints under various loading conditions. Adv Eng Softw 2010;41:580–9. [23] Sonmez M. Discrete optimum design of truss structures using artificial bee colony algorithm. Struct Multidiscip Optim 2011;43:85–97. [24] Bendsøe MP, Sigmund O. Material interpolation schemes in topology optimization. Arch Appl Mech 1999;69:635–54. [25] Sigmund O, Peterson J. Numerical instabilities in topology optimization: a survey on procedures dealing with checkerboards, mesh-dependencies and local minima. Struct Optim 1998;16:68–75. [26] Huang X, Xie YM. Convergent and mesh-independent solutions for the bidirectional evolutionary structural optimization method. Finite Elem Anal Des 2007;43:1039–49. [27] Taylor JE, Bensøe MP. A mutual energy formulation for optimal structural design. Struct Multidiscip Optim 2001;22:95–101.
Z.H. Zuo, Y.M. Xie / Computer-Aided Design 56 (2014) 58–67 Zhi Hao Zuo obtained his bachelor degree in Civil Engineering from Tongji University in China. After he completed his Ph.D. in topology optimization of periodic structures in 2009 in RMIT University in Australia, he has been working as a research fellow in the Centre for Innovative Structures and Materials, at RMIT University.
67
Yi Min Xie obtained his bachelor degree in Engineering Mechanics from Shanghai Jiao Tong University in China and Ph.D. in Civil Engineering from Swansea University in Britain. He is one of the pioneers of the evolutionary structural optimization (ESO) method. He is an elected fellow of the Australian Academy of Technological Sciences and Engineering.