PDF File - MDO Lab
An Isoparametric Approach to Level Set Topology Optimization Using a Body-Fitted Finite Element Mesh Kai A. Jamesa , Joaquim R.R.A. Martinsb a University b University
of Toronto Institute for Aerospace Studies, Toronto, Ontario, Canada of Michigan Department of Aerospace Engineering, Ann Arbor, Michigan, United States
Abstract We present a new method for performing level set topology optimization of structures confined to irregularly-shaped physical domains. The shape sensitivities are computed in physical space using a body-fitted, non-uniform finiteelement mesh. These values are then mapped to the Euclidean computational space using the Jacobian matrix of the isoparametric transformation taking place within each cell. The method is demonstrated on a three-dimensional problem in which we minimize the compliance of a contoured wingbox structure. Results show that the method converges consistently, and generates solutions comparable to those obtained using the SIMP method. Keywords: topology optimization; level set method; body-fitted mesh; finite element analysis; isoparametric finite elements 1. Introduction The level set method has been shown to be an effective tool for topology optimization of structures. This is evidenced by the large number of papers published on the topic within the past several years [1, 2, 3]. As a result, the level set method now rivals more established topology optimization methods, such as Solid Isotropic Material with Penalization (SIMP) [4], in terms of popularity. However due to the relative difficulty in implementing the level set method when compared with SIMP, the vast majority of the research devoted to the level set method has focused on a small class of problems, while seeking to improve the numerical performance of the method [5]. Consequently, there remain large classes of problems to which the level set method has yet to be applied. One such area involves the optimization of structures confined to irregularly-shaped physical domains. These problems require the use of body-fitted finiteelement meshes, which are typically non-uniform and non-rectilinear. Such problems have been solved successfully by SIMP researchers on a number of occasions [6, 7, 8, 9]; however, they present a unique challenge for users of the level set method. The Hamilton– Jacobi equation, which is at the heart of the level set method, is solved on a uniform Cartesian grid, which typically coincides with the finite-element mesh. By usPreprint submitted to Computers & Structures
ing a Cartesian computational grid, the sensitivities can be calculated at the Gauss points of the finite elements and passed directly into the Hamilton–Jacobi equation as advection velocities, using a straightforward one-toone mapping. Some authors have developed techniques for solving the Hamilton-Jacobi equation in cases where the finiteelement mesh is either unstructured or simply nonuniform [10, 2]. One approach is to construct a continuous representation of the sensitivity field by interpolating the sensitivity values calculated at the Gauss points of the finite elements [10]. This sensitivity field is then sampled at the nodes of the uniform Cartesian grid in order to obtain a set of advection velocities with which to solve the Hamilton–Jacobi equation. Interpolation of the sensitivity field can be handled in a variety of ways including the use of radial basis functions [10], as well as least squares fitting [1] or boundary elements [11]. While these methods allow for application of the level set method to any arbitrary structural mesh, they do not address situations where the working domain of the problem is non-rectangular. These situations arise frequently in engineering optimization, as it is often necessary to design a structure that must fit entirely within a non-rectangular region. In the following sections of this paper, we describe a new method for optimizing structures confined to such irregular domains. As in previous methods, the Hamilton–Jacobi equation is solved November 21, 2011
on a uniform rectilinear grid; however, the finite element modeling of the structure is performed using a fixed non-uniform structured mesh whose cells represent units in actual physical space. The structural sensitivities are calculated on the non-uniform mesh and then mapped into computational space using an isoparametric transformation within each cell. Because each element in the finite element mesh has its own unique cell in the rectilinear computational grid, we retain a oneto-one mapping between the structural sensitivities and advection velocities and avoid significant loss of robustness. This method is demonstrated in two ways: first, on a series of benchmark examples based on a variation of the classic two-dimensional L-bracket problem, and second, on a three-dimensional problem in which the working domain is given by the outer surface of an aircraft wing. We also introduce a novel approach to enforcing constraints based on the augmented Lagrangian formulation.
Figure 1: Sample structural design problem using level set parametrization. Here ΓD denotes the fixed boundary and ΓN denotes the boundary on which surface tractions are applied
min J Ω
subject to: c = 0 c˜ ≤ 0 Kd − F = 0 ( 1, ρ(xi ) = 10−3 ,
2. The Level Set Method
(1) xi ∈ Ω ¯ xi ∈ Ω
k(x) = ρ(x)k0 2.1. Problem Formulation Here, J is an arbitrary objective function that is dependent on the design variable, Ω. The structure may be subject to equality and inequality constraints, c and c˜ , and must satisfy the governing equation Kd − F = 0, where K is the global stiffness matrix, d is the vector of nodal displacements and F is the vector of applied forces. Elements located inside the solid region Ω are assigned a material density of unity, while those lying outside the boundary are given some small non-zero density, ρmin , so that they mimic the behaviour of a void space but can still be included into the global stiffness matrix without causing singularities. The densities of elements that are bisected by the structural boundary, ∂Ω, are interpolated based on the fraction of that element’s volume that lies inside Ω. The elasticity modulus, E, of a given element is the product of the element’s relative density and the elasticity of the element in the solid phase, E0 . The material properties are taken to be piecewise constant throughout each element.
The level set method was developed in the late 1980s by Osher and Sethian [12] as a method for tracking front propagation. It has been especially popular among researchers in fluid dynamics and computer vision, though an increasing number of authors have begun to apply the technique to structural topology optimization. Although similar to SIMP methods in some respects, the level set method differs in the way in which the structural shape and topology are parameterized. In both methods, the design domain is discretized into a series of finite elements whose relative material density, ρ, is allowed to vary continuously between 0 and 1 with the two extremes representing void and solid elements, respectively. In the SIMP approach, one optimizes these material densities directly, whereas in the level set method, one optimizes the location of the material boundary, ∂Ω, and from the location of this boundary, one determines which elements are solid and which are void, as shown in (1). Due to this boundary-based parameterization, the level set method avoids meshdependence and checkerboarding [3, 13] – two of the main numerical challenges associated with the SIMP method [14, 15].
Like the SIMP method, the level set approach begins by defining a bounded domain D ⊂ 0; x ∈ (D ∩ Ω). ψ
5 0
By implicitly representing the material boundary using a higher-order level set function, one allows for changes in topology, such as the merging of holes (in two dimensions) or cavities (in three dimensions) in the structure. The level set function is defined discretely at the nodes of a structured Cartesian mesh. These data are periodically interpolated in order to extract the precise location of the zero level set. While the level set method is mesh-independent, the converged solution of the optimization problem is dependent upon the initial topology [3, 13]. Therefore, one must carefully select the number of holes in the initial design according to the desired length scale for the final solution. Once the initial topology and material boundary are defined, the level set function is initialized as the signed distance of each point in the domain from the boundary, with negative distances used for points inside the solid region Ω and positive distances outside.
−5 −10 40 30 40 30
10
20 0
Y
10 0
X
Figure 2: Two-dimensional example of an optimized level set function, with contour lines displayed in the xy-plane
this way, the Hamilton–Jacobi equation replaces the optimization algorithms used in element-based topology optimization methods such as SIMP. The optimization reaches convergence once the advection velocities are within a small tolerance of zero along the material interface. As noted by Allaire et al. [3], the constrained optimization problem is generally not convex despite the absence of a penalization on intermediate densities. As a result, convergence to a global minimum cannot be guaranteed. Figure 2 shows an optimized level set function for a two-dimensional problem with the contour plots displayed in the xy-plane. If one takes the contour line corresponding to the surface ψ = 0 (the cyan-coloured contour) this plot will give the material boundary of the optimized structure. Because the Hamilton–Jacobi update assumes that the spatial derivative is valid over the distance vdt, it is necessary to periodically reinitialize the level set function to a form that validates this assumption. Therefore after every few updates, one must recover the signed distance form of the level set function, while maintaining the location of the current zero level set [3]. The results presented rely on an implementation of the fast sweeping method [17, 18, 19].
2.2. The Hamilton–Jacobi Equation In each optimizing iteration, the level set function is updated using the Hamilton–Jacobi equation (2). ∂ψ + v|∇ψ| = 0, (2) ∂t Rearrangement of this equation leads to the following scalar formula for updating the level set function at each point in the domain: ψt+dt = ψt − v|∇ψt |dt
60 50
20
(3)
Here the divergence ∇ψ is computed numerically at each point using finite differencing of the values of the level set function at adjacent nodes. The advection velocity v is given by the shape sensitivity of the objective function at each point. The variable t is a fictitious time parameter introduced to track the evolution of the level set function over the course of the optimization. The update described in Eq. (3) corresponds to moving the material boundary by a distance of vdt in the outward normal direction from the interface. The time step dt is typically chosen so that the Courant-Friedrichs-Lewy (CFL) condition is satisfied [16]. When the advection velocity is given by the shape sensitivities of the objective function, each Hamilton– Jacobi update is equivalent to a steepest descent step. In
2.3. Sensitivity Analysis For an arbitrary objective function, the shape sensitivity is defined according to the following equation, which uses the Fr´echet functional derivative [20], J 0 (Ω)(θ) = 3
Z vθ · nds, ∂Ω
(4)
where θ is an arbitrary small vector field, and the advection velocity v is equal to the local shape sensitivity of the objective J at a given point on the boundary. The shape sensitivity is usually some function of the displacement state u and can be derived using the following identity for generalized shape derivatives. For objectives of the form Z
J(Ω) =
Ω
j(x)dx +
L(Ω) = J(Ω) − λ · c(Ω),
where c is the vector of equality constraints satisfying c(Ω) = 0. This allows us to perform unconstrained optimization on the Lagrangian with respect to the design variable Ω as well as the vector of Lagrange multipliers λ. Convergence is reached once the first order KKT conditions are satisfied. The shape variable Ω is updated using the Hamilton– Jacobi equation, while the vector λ is updated using the following heuristic,
Z (5)
l(x)ds, ∂Ω
where Ω is a smooth, bounded open set, we have
J 0 (Ω)(θ) =
Z
+
Z
⇒ J 0 (Ω)(θ) = +
Ω
θ(x) · n(x)
Z∂Ω
λk+1 = λk + rc(Ω).
(6)
div(θ(x) j(x))dx
! ∂l(x) + H(x)l(x) ds ∂n(x)
θ(x) · n(x) j(x)ds
∂l(x) + H(x)l(x) ds, ∂n(x) !
θ(x) · n(x) ∂Ω
(9)
This update corresponds to a steepest descent step whose length is obtained by taking the first derivative of the Lagrangian with respect to the Lagrange multiplier. The step size r is a fixed value chosen to obtain the best trade-off between convergence time and stability. Consequently, this method is relatively inexpensive and works well for problems involving a single constraint.
(7)
∂Ω
Z
(8)
where H is the mean curvature of ∂Ω given by H = ∇n [3]. In the current study, all sensitivities are evaluated using the adjoint method derived by Allaire et al. [3].
3. The Isoparametric Formulation The shape sensitivities derived in Section 2.3 assume a smooth, continuous displacement state u, and strain field , which can be evaluated at any location throughout the working domain. In practice, the finite element method is used to generate a discrete representation of the displacement state, and the shape sensitivities are evaluated at the Gauss quadrature points for each element [21]. Alternatively the sensitivities can be averaged over the domain of the element, in which case the sensitivity is assumed to be piecewise constant across each element. Whereas most level set frameworks use square and cubic finite elements to coincide with the Cartesian computational grid, the isoparametric formulation allows for the use of arbitrary quadrilateral and hexahedral elements. Therefore, the isoparametric level set method relies heavily upon the isoparametric finite element formulation [24]. As with the standard approach, the Hamilton–Jacobi equation is solved on a uniform Cartesian grid with unit spacing between the nodes in all directions. This uniform, orthogonal spacing means that the spatial derivatives of the level set function can be evaluated simply by taking finite differences of the function evaluated at adjacent grid points. Each node in this computational grid corresponds to a four-node element (eight-node in three dimensions) in the finite
2.4. Constraint Handling Few papers address the topic of constrained optimization within the context of the level set method. Most authors account for constraints by adding a Lagrangianstyle constraint term to their objective. Some use a constant Lagrange multiplier [3], while others update the Lagrange multiplier using a bisection method [21]. The bisection approach is similar to the optimality criteria method, which has also been used in combination with SIMP frameworks [22]. However, like the optimality criteria method, the bisection approach is only applicable to problems involving a single constraint. Given this limitation, other researchers have proposed alternative methods for problems with multiple constraints [23], but no method has yet been widely adopted, as most lack flexibility, robustness or, both. Here, we introduce a new method based on an augmented Lagrangian approach, where the Lagrange multiplier is updated according to the Karush-Kuhn-Tucker (KKT) conditions. To apply equality constraints, we construct a Lagrangian function L by adding an additional term to the objective for each constraint to be handled, i.e., 4
Consider an infinitesimal patch located at some point p on the surface of the material interface. Since the orientation of the local coordinate axes is arbitrary, we can choose them so that the ζ axis points in the direction normal to the patch. Therefore, the area of the patch, measured in computational space, is given by dξdη. If the patch moves a distance dζ in the outward normal direction, it will trace a prism whose volume is given by dξdηdζ, when measured in computational space. This prism represents a total volume of dξdηdζ|J| in physical space. We now define a new quantity that is hereafter referred to as the relative impact. The relative impact of a shape perturbation at a point p on the material surface is defined as the product of the incremental volume Q caused by that perturbation, and the shape sensitivity of the objective function at the point p. Given a shape sensitivity, v p , computed in physical space at point p, we choose an equivalent computational sensitivity, vc , such that the relative impact of an infinitesimal shape perturbation at p, when measured in computational space, is maintained. Therefore we can deduce the following relationship:
element mesh (see Fig. 3). The shape sensitivities are evaluated in physical space using the non-uniform finite-element mesh (one sensitivity value for each element), and each of these values is transformed into its representative analog in computational space before being passed into the Hamilton–Jacobi equation. This transformation is based on the isoparametric transformation required to morph a unit square (or cube) into the quadrilateral (or hexahedral) shape of the element in question. The Jacobian for this transformation is defined in Eq. (12). Within each element one defines a local coordinate system ξ, η, ζ, where the orthogonal basis vectors ~ξ, ~η, ~ζ represent the directions in computational space. Figure 4 shows the two-dimensional local coordinate system associated with the quadrilateral element whose nodes are located at {xi , yi }, as expressed in physical space. Note that in the figure, the computational coordinates ξ, η range from −1 to 1. This choice of domain facilitates Gauss quadrature integration. This implies that each element occupies a volume of 4, as opposed to unity, in computational space. However, this quantity is inconsequential as long as all elements have the same volume in computational space. The physical location {x, y, z} of any point within the element can be expressed as a weighted sum of the node locations, x=
8 X
Ni (ξ, η, ζ)xi ,
where the weights, Ni , are the linear shape functions of the finite element: (1 + ξξi )(1 + ηηi )(1 + ζζi ) , ξi , ηi , ζi = ±1, 8 (11) which depend on the local coordinates {ξ, η, ζ}. Given this relationship between the physical coordinates and the local or computational coordinates, we can define the Jacobian matrix of the coordinate transformation within each element: Ni (ξ, η, ζ) =
∂x ∂ξ ∂y ∂ξ ∂z ∂ξ
∂x ∂η ∂y ∂η ∂z ∂η
∂x ∂ζ ∂y ∂ζ ∂z ∂ζ
.
(13)
⇒ dξdηdζvc = dξdηdζ|J|v p
(14)
⇒ vc = |J|v p .
(15)
Because the shape sensitivities are computed at discrete locations corresponding to the finite-element mesh, we approximate the Jacobian as being piecewise constant throughout each element. To do this we approximate |J| as being equal to the determinant of J averaged over the element volume, which turns out to be equal to the volume of the element itself, i.e., Z vole = det(J(ξ, η, ζ))dξdηdζ. (16)
(10)
i=1
J(ξ, η, ζ) =
Q c vc = Q p v p
Ωe
Therefore given some shape sensitivity, v p , evaluated in physical space, the resulting advection velocity in computational space is given by vc = vole v p . Thus, using the set of transformed sensitivities as the advection velocities in the Hamilton–Jacobi equation corresponds to taking a step in the direction of steepest descent.
(12)
3.1. The Constitutive Relation The constitutive matrix, which defines the relationship between the stress tensor, σ, and the strain tensor, ε, is given in Eq. (18) for the two-dimensional case
From here, we can derive the relationship between the sensitivities as calculated in physical space and the equivalent sensitivities to be used in computational space.
σ = Dε 5
(17)
Figure 3: Two-dimensional mapping of an arbitrary structured mesh to the uniform rectilinear computational mesh.
Figure 4: Two-dimensional illustration of the mapping from local to global coordinates for an arbitrary quadrilateral element
6
σ xx σyy σ xy
1 E = ν (1 − ν2 ) 0
ν 1 0
0 ε xx 0 εyy 1−ν ε xy 2
dividing the structure into equally-sized blocks of elements and assigning each block to a processor. (18)
3.2. Compliance Minimization Compliance minimization is commonly used among researchers of the level set method [3, 1, 16] to demonstrate new concepts, largely because of the ease with which it can be implemented. As shown by Allaire et al. [3], for the compliance objective function (i.e. the work done by external forces, Eq. (23)) the shape sensitivity reduces to the form shown in Eq. (24):
Here E and ν are the Young’s modulus and Poisson’s ratio of the isotropic material. Note that in all twodimensional examples, the plane stress condition is assumed. As is typical among previous studies, we use the ersatz material approach [3, 25], in which constitutive relation is tailored to account for those elements which lie, either partially or fully, outside the material domain. For these elements, the Young’s modulus is interpolated according to the portion of the element’s volume that falls inside the boundary, i.e., Ee = ρe E0 ⇒ De = ρe D0 ,
Jcomp (Ω) = =
1 vole
ρe ≈
1 vole
Z Ωe np X
Ω
g(x) · u(x)dx +
Z Ω
Z ΓN
f (x) · u(x)ds
εT (x)DεT (x) dx
(23)
(19) (20)
0 Jcomp (Ω)(θ) =
where, Ee is the effective Young’s modulus of the element, and E0 is the Young’s modulus of the solid material. The variable ρe represents the element volume fraction, which can also be interpreted as the relative material density. Under the isoparametric formulation, ρe is evaluated using the integral, ρe =
Z
h(−ψ(ξ, η))det(J(ξ, η))dξdη
(21)
h(−ψ p )|J p |.
(22)
Z Γ0
θ(x) · n(x) εT (x)DεT (x) ds
(24)
In the above equations, g and f denote body forces and surface tractions respectively, while u represents the continuous displacement field. The domain Γ0 refers to the traction-free surface, which is the only portion of the material boundary that varies during the optimization. Despite this restriction, the resulting advection velocity formula v = εT (x)DεT (x) is used for all points within the working domain when solving the Hamilton–Jacobi equation so that continuity is maintained. The compliance shape sensitivity, which is equal to twice the local strain energy density, is treated as being piecewise constant within each element. Therefore the advection velocity for each element is approximated as being proportional to the average strain energy density within that element. Under the isoparametric formulation, the calculation proceeds as follows for the twodimensional case: Z 1 Ue = εT (ξ, η)De ε(ξ, η)|J|dξdη, (25) 2 Ωe
p
The symbol h denotes the Heaviside function, which is equal to 1 wherever the level set function, ψ, is negative, and is equal to 0 everywhere else. In order to evaluate this integral numerically, the element is divided into a uniform, Cartesian grid of pixels, where the number of pixels, n p is equal to or greater than 1/ρmin . Each pixel is either solid or void depending on the value of the level set function at that pixel’s location, ψ p . In order to evaluate the level set function at these off-node locations, we simply interpolate using the linear shape functions shown in Eq. (11). The integral is then approximated as a weighted sum of solid pixels, with the weighting given by the determinant of the Jacobian matrix. Interpolating the level set function in this way can become computationally expensive, especially for threedimensional problems. This challenge can be mitigated significantly, however, through the use of parallel processing. We have vastly improved computation time by
where Ue is the total strain energy for element e. In order to obtain the average strain energy density, we divide this integral by the element volume given in Eq. (16). The expression can be simplified by introducing the strain-displacement matrix, B, which is comprised of the set of first derivatives of the element shape functions (Eq. (11)). We can therefore express the strain in terms of the strain displacement matrix and the vector of element nodal displacements de . 7
ε(ξ, η) = BT (ξ, η)de
This problem is solved using two different geometries. In the long L-bracket problem, the vertical and horizontal segments of the bracket have an aspect ratio of 2. In the short L-bracket problem, shown in Fig. 6, this aspect ratio is reduced to 1. This geometry causes the individual elements to have larger aspect ratios, particularly near the bottom-left corner of the domain. In both cases the applied force has a magnitude of 1.
(26)
Finally we introduce the element stiffness matrix, ke , which is defined as
k0 =
Z
ke =
Z
Ωe
Ωe
BT (ξ, η)D0 B(ξ, η)|J|dξdη
(27)
BT (ξ, η)De B(ξ, η)|J|dξdη
(28)
⇒ ke = ρe k0 ,
(29)
and the advection velocity becomes vp =
ρe dTe k0 de . vole
(30)
In order to obtain the computational velocity to be used in performing the Hamilton–Jacobi update, we multiply this expression by the element volume. Therefore the computational advection velocity reduces to
F vc =
ρe dTe k0 de .
(31)
In order to enforce a volume constraint we must find the shape sensitivity of the structural volume. From Eq. (6), we see that the shape sensitivity of the volume function is simply equal to unity. The equivalent sensitivity in computational space is simply the element volume. Combining the two functions according to the Lagrangian formulation defined above, the computational advection velocity for the compliance minimization problem subject to a constraint on structural volume is given by Eq. (32):
(a) Initial topology and loading conditions
90 80 70 60 50 40
vc = ρe dTe k0 de − λvole .
(32)
30 20
4. Numerical Examples 10
The isoparametric formulation derived above is demonstrated in a series of two-dimensional and threedimensional sample problems. The following examples are based on the classic L-bracket problem, which has been used frequently in topology optimization studies [5, 26]. In all L-bracket sample problems, the mesh is comprised of a lattice of trapezoidal elements as shown in Figure 5(b). This mesh is mapped to a single uniform rectangular mesh using a one-to-one mapping. The L-bracket problem is solved for minimum compliance subject to a constraint on the structural volume so that the total volume of the optimized structure does not exceed 40% of the volume of the working domain.
0
0
20
40
60
80
(b) Finite-element mesh
Figure 5: The long L-bracket problem
Figs. 7 and 8 show the optimized topologies for the long and short L-bracket problems. In each figure the zero contour of the level set function (i.e. ψ = 0) is plotted. In both cases, the isoparametric level set method produces converged, feasible solutions. The 8
F
Figure 7: Optimized topology for the short L-bracket problem (a) Initial topology and loading conditions 70
60
50
40
30
20
10
0
0
10
20
30
40
50
60
Figure 8: Optimized topology for the long L-bracket problem
70
(b) Finite-element mesh 2000
1
Figure 6: The short L-bracket problem
Compliance
zero contour of the level set function reveals straight, well-defined members and is unaffected by the lack of uniformity in the finite-element mesh. The convergence plots for both problems are shown in Figures 9 and 10. The plots indicate smooth convergence to a local minimum. Also, while the constraint function does oscillate during the optimization, the objective is reduced monotonically. Careful selection of the coefficient r in Eq. (9) can be used to reduce the amplitude and number of these oscillations. An effective choice is to use r ≈ 0.1λ0 , where λ0 is the initial value of the Lagrange multiplier. We now compare the solutions displayed in Figs. 7 and 8 with solutions obtained using the SIMP method. In both SIMP solutions, the finite-element mesh is iden-
1600
0.8
1200
0.6
800
0.4
400
0.2
0
0
20
40
60
80 100 120 Iteration Number
140
160
Volume Fraction
Compliance Volume Fraction
0 180
Figure 9: Convergence history for the long L-bracket compliance minimization problem
tical to that used above. We use a density filtering 9
500
1
400
0.8
300
0.6
200
0.4
100
0.2
Volume Fraction
Compliance
Compliance Volume Fraction
(a) Level set solution
(b) SIMP solution
Figure 12: Optimized density distribution for the short L-bracket 0
0
50
100 150 Iteration Number
200
0 250 800
Figure 10: Convergence history for the short L-bracket compliance minimization problem
level set method SIMP
700 600
Compliance
method to avoid checkerboarding [27, 28], in combination with the optimality criteria method [29] to carry out the constrained optimization problem. In the case of the short L-bracket, the isoparametric level set method generates a solution nearly identical to that of the SIMP method, as illustrated in Figure 12. In contrast, the two solutions to the long L-bracket problem differ greatly from one another, although both solutions exhibit similar compliance, as shown in Table 1.
500 400 300 200 100 0
0
50
100 150 Iteration Number
200
250
Figure 13: Comparison of the convergence histories of the SIMP and level set methods for the short L-bracket problem
(a) Level set solution
ture, where material is distributed throughout a circular domain, which contains a void region in the center. Fig. 14 shows the problem geometry, as well as the support and loading conditions. The optimized density distribution for the cantilevered ring problem is shown in Fig. 15. As is the case with the short L-bracket, the isoparametric level set method and the SIMP method yield very similar structures. Table 1 confirms that in this case, as well as in both L-bracket problems, the isoparametric level set solution exhibits very similar compliance to the SIMP solutions. Also, from the convergence plot given in Fig. 13, we see that the two methods have comparable computational cost in terms of the number of structural analyses required to reach convergence. Another problem that requires the use of a nonuniform mesh is the topology optimization of an aircraft wingbox. In the following example, the working domain of the problem is roughly determined by the shape of the wing’s outer surface, and material can be placed anywhere within this three-dimensional region.
(b) SIMP solution
Figure 11: Optimized density distribution for the long L-bracket
While the above examples make effective use of the isoparametric level set method, both L-bracket problems could have been solved using uniform Cartesian meshes. In addition to these problems, the isoparametric formulation enables us to apply the level set method to an entire class of problems that would otherwise be beyond the reach of most level set algorithms. We provide two examples of such problems. The first example is a novel problem designed specifically to test the effectiveness of the isoparametric formulation. This problem involves the optimization of a cantilevered ring struc10
level set method
SIMP
long L-bracket
364.44
366.18
short L-bracket
83.84
84.31
cantilevered ring
6.07
5.84
Table 1: Minimum compliance values for the SIMP and level set solutions. In all cases, the volume fraction is constrained at 0.4. (Note that the compliance values for the SIMP solutions have been calculated without any penalization.) (a) Initial topology 30
20
10
0
−10
Figure 16: Working domain and finite-element mesh for the wingbox optimization problem
−20
−30 −30
−20
−10
0
10
20
30
(b) Finite-element mesh
Figure 14: The two-dimensional cantilevered ring problem Figure 17: Loading conditions for the wingbox optimization problem
(a) Level set solution
The cross-section of the domain is based on a symmetric two-dimensional airfoil with the leading and trailing edges removed in order to replicate the shape of a wingbox. The resulting cross-section is then extruded to an aspect ratio of 3.02. The structure also has a taper ratio of 0.91, and a leading edge sweep angle of 9.4◦ . Fig. 16 shows the working domain and the finiteelement mesh for the problem, while Fig. 17 shows the loading conditions. The top surface of the structure is subject to a uniform distributed load acting in the downward direction. To the bottom surface we apply another distributed load, which is constant in the chordwise (x)
(b) SIMP solution
Figure 15: Optimized designs for the cantilevered ring problem
11
direction, and tapers off elliptically in the spanwise (z) direction. The structure is fixed along the face located at the root of the wing, where the wing connects with the fuselage of the aircraft. The structural mesh is comprised of 32 × 16 × 96 linear, eight-node hexahedral elements in the x, y, and z directions respectively. In this example, we use a fixed Lagrange multiplier and begin with an initial volume fraction of 60%. The Lagrange multiplier is set according to λ = 0.65C0 /V f0 , where C0 is the initial compliance and V f0 is the initial volume fraction. A minimum skin thickness is enforced on the top and bottom surfaces of the structure. At a given location along the span, the local skin thickness is equal to 1/2 the thickness of the local surface element measured in the y-direction. This thickness requirement is enforced by setting the minimum relative material densities of all top and bottom surface elements to ρmin = 0.5. The structure is optimized for minimum compliance, with a 40% volume constraint.
tom surfaces of the structure. The optimized structure can be interpreted as being comprised of two main types of component. A large portion of material is devoted to reinforcing the top and bottom skin of the structure, particularly near the root, where the fixed boundary condition is applied and where the internal bending moment of the structure is the largest. Also, as illustrated in Fig. 20, the remainder of the material is used to create a series of spar-like structures. These consist of truss-like members aligned at roughly 45◦ from the horizontal xzplane, in order to resist internal shear force in the z direction. Fig. 21 also reveals that these members are slanted slightly in the chordwise direction. This reflects the asymmetry of the structure due to the sweep angle. Because the wing is swept, any upward bending also results in a torsional moment. The chordwise slant of the spar members allows the structure to resist the twisting caused by this torsional moment.
Figure 20: Span-wise view of the optimized wing structure
Figure 18: Initial wing topology
Figure 21: Chord-wise view of the optimized wing structure
5. Conclusions In this paper we have presented a robust and computationally efficient method for performing structural topology optimization using the level set method combined with a non-uniform finite-element mesh. We have demonstrated how this isoparametric formulation can be used to effectively optimize structures that are confined to irregularly-shaped domains. In addition to the isoparametric formulation, we have also presented a
Figure 19: Optimized wing topology
Figs. 18 and 19 show the initial and optimized internal topologies for the wingbox problem. The figures do not include the wing skin, which covers the top and bot12
new method for enforcing equality constraints based on the augmented Lagrangian approach. These methods have been demonstrated on a series of two- and threedimensional structures optimized for minimum compliance subject to weight constraints. The results indicate that the method generates converged, feasible solutions that are comparable in performance to solutions obtained using the SIMP method. Optimization of structures modeled using body-fitted meshes represents a large and vitally important class of problems in structural optimization. The developments presented are aimed at extending the range of applicability of the level set method, thereby allowing researchers to use the method for tackling problems that more closely resemble those encountered in real-world applications. An important component of the push toward handling more realistic design problems is the implementation of additional objective and constraint functions. Although the current study focuses on compliance minimization subject to a volume constraint, the isoparametric method can also be applied to other objectives without any modification to the general algorithm. Going forward, some useful problems to explore include mass minimization subject to yield failure constraints — as is often done in the design of aerospace structures — as well as maximization of failure load with respect to elastoplastic collapse, as one might encounter in a structural engineering context.
[6] L. Krog, A. Tucker, Topology optimization of aircraft wing box ribs, in: 10th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference, Albany, New York, 2004, pp. AIAA 2004–4481. [7] T. H. Lee, S. Y. Han, J. K. Lim, Topology optimization of the inner reinforcement for an automobile hood using modal design sensitivity analysis, Key Eng. Mat. 183–187 (2000) 439–444. [8] S. Jin-fa, S. Jian-hui, Overview on innovation of topology optimization in vehicle cae, in: Electronic Computer Technology, 2009 International Conference on, 2009, pp. 457–460. [9] D.-F. Wang, B. Zhang, J. Chen, S.-P. Wang, S.-Q. Liu, Application of topology optimization to commercial vehicle cab design, in: International Conference on Measuring Technology and Mechatronics Automation, Vol. 2, Zhangjiajie, Hunan, China, 2009, pp. 763–766. [10] S.-H. Ha, S. Cho, Level set based topological shape optimization of geometrically nonlinear structures using unstructured mesh, Comput. Struct. 86 (13-14) (2008) 1447 – 1455. [11] K. Abe, S. Kazama, A boundary element approach for topology optimization problem using the level set method, Commum. Numer. Meth. Eng. 23 (5) (2007) 405–416. [12] S. Osher, J. A. Sethian, Fronts propagating with curvaturedependent speed: algorithms based on Hamilton–Jacobi formulations, J. Comput. Phys. 79 (1988) 12–49. [13] H. Jia, H. G. Beom, Y. Wang, S. Lin, B. Liu, Evolutionary level set method for structural topology optimization, Comput. Struct. 89 (2011) 445–454. [14] K. A. James, J. S. Hansen, J. R. R. A. Martins, Structural topology optimization for multiple load cases using a dynamic aggregation technique, Eng. Optimiz. 41 (2009) 1103–1118. [15] O. Sigmund, Numerical instabilities in topology optimization: A survey on procedures dealing with checkerboards, meshdependencies and local minima, Struct. Optimiz. 16 (1998) 68– 75. [16] S. Wang, M. Y. Wang, A moving superimposed finite element method for structural topology optimization, Int. J. Numer. Meth. Eng. 65 (2006) 1892–1922. [17] T. Oberhuber, Numerical recovery of the signed distance function, in: Czech-Japanese Seminar in Applied Mathematics, Prague, Czech Republic, 2004, pp. 148–164. [18] J. A. Sethian, A fast marching level set method for monotonically advancing fronts, in: Proc. Nat. Acad. Sci, 1995, pp. 1591– 1595. [19] Y. R. Tsai, L. Cheng, S. Osher, H. Zhao, Fast sweeping algorithms for a class of hamilton-jacobi equations, SIAM J. Numer. Anal. 41 (2003) 673–694. [20] J. Simon, Differentiation with respect to the domain in boundary value problems, Numer. Funct. Anal. and Optimiz. 2 (1980) 649–687. [21] S. Wang, K. Lim, B. Khoo, M. Wang, An extended level set method for shape and topology optimization, J. Comp. Phys. 221 (1) (2007) 395 – 421. [22] O. Sigmund, A 99 line topology optimization code written in matlab, Struct. Multidiscip. Optimiz. 21 (2001) 120–127. [23] M. Yulin, W. Xiaoming, A level set method for structural topology optimization with multi-constraints and multi-materials, Acta Mech. Sinica 20 (2004) 507–518. [24] B. M. Irons, Engineering applications of numerical integration in stiffness methods, AIAA J. 4 (11) (1966) 2035–2037. [25] Z. Luo, L. Tong, Z. Kang, A level set method for structural shape and topology optimization using radial basis functions, Comput. Struct. 87 (7-8) (2009) 425–434. [26] C. Le, J. Norato, T. Bruns, C. Ha, D. Tortorelli, Stress-based topology optimization for continua, Struct. Multidiscip. Optimiz. 41 (2010) 605–620.
Acknowledgements This research was supported by the Natural Sciences and Engineering Research Council of Canada and the Ontario Graduate Scholarship program.
References [1] P. Dunning, H. Kim, Investigation and improvement of sensitivity computation using the area-fraction weighted fixed grid fem and structural optimization, Finite Elem. Anal. Des. 47 (2011) 933–941. [2] Z. Liu, J. G. Korvink, Adaptive moving mesh level set method for structure topology optimization, Eng. Optimiz. 40 (6) (2008) 529–558. [3] G. Allaire, F. Jouve, A.-M. Toader, Structural optimization using sensitivity analysis and a level-set method, J. Comput. Phys. 194 (1) (2004) 363–393. [4] G. I. N. Rozvany, M. Zhou, T. Birker, Generalized shape optimization without homogenization, Struc. and Multidisc. Optim. 4 (1992) 250–252. [5] G. Allaire, F. Jouve, Minimum stress optimal design with the level set method, Eng Anal. Bound. Elem. 32 (11) (2008) 909– 918.
13
[27] T. E. Bruns, D. A. Tortorelli, Topology optimization of nonlinear elastic structures and compliant mechanisms, Comput. Method Appl. Mech. Eng. 190 (26-27) (2001) 3443–3459. [28] B. Bourdin, Filters in topology optimization, Int. J. Numer. Meth. Eng. 50 (2001) 2143–2158. [29] M. P. Bendsøe, Optimization of structural topology, shape and material, Springer, Berlin–Heidelberg–New York, 1995.
14
of Toronto Institute for Aerospace Studies, Toronto, Ontario, Canada of Michigan Department of Aerospace Engineering, Ann Arbor, Michigan, United States
Abstract We present a new method for performing level set topology optimization of structures confined to irregularly-shaped physical domains. The shape sensitivities are computed in physical space using a body-fitted, non-uniform finiteelement mesh. These values are then mapped to the Euclidean computational space using the Jacobian matrix of the isoparametric transformation taking place within each cell. The method is demonstrated on a three-dimensional problem in which we minimize the compliance of a contoured wingbox structure. Results show that the method converges consistently, and generates solutions comparable to those obtained using the SIMP method. Keywords: topology optimization; level set method; body-fitted mesh; finite element analysis; isoparametric finite elements 1. Introduction The level set method has been shown to be an effective tool for topology optimization of structures. This is evidenced by the large number of papers published on the topic within the past several years [1, 2, 3]. As a result, the level set method now rivals more established topology optimization methods, such as Solid Isotropic Material with Penalization (SIMP) [4], in terms of popularity. However due to the relative difficulty in implementing the level set method when compared with SIMP, the vast majority of the research devoted to the level set method has focused on a small class of problems, while seeking to improve the numerical performance of the method [5]. Consequently, there remain large classes of problems to which the level set method has yet to be applied. One such area involves the optimization of structures confined to irregularly-shaped physical domains. These problems require the use of body-fitted finiteelement meshes, which are typically non-uniform and non-rectilinear. Such problems have been solved successfully by SIMP researchers on a number of occasions [6, 7, 8, 9]; however, they present a unique challenge for users of the level set method. The Hamilton– Jacobi equation, which is at the heart of the level set method, is solved on a uniform Cartesian grid, which typically coincides with the finite-element mesh. By usPreprint submitted to Computers & Structures
ing a Cartesian computational grid, the sensitivities can be calculated at the Gauss points of the finite elements and passed directly into the Hamilton–Jacobi equation as advection velocities, using a straightforward one-toone mapping. Some authors have developed techniques for solving the Hamilton-Jacobi equation in cases where the finiteelement mesh is either unstructured or simply nonuniform [10, 2]. One approach is to construct a continuous representation of the sensitivity field by interpolating the sensitivity values calculated at the Gauss points of the finite elements [10]. This sensitivity field is then sampled at the nodes of the uniform Cartesian grid in order to obtain a set of advection velocities with which to solve the Hamilton–Jacobi equation. Interpolation of the sensitivity field can be handled in a variety of ways including the use of radial basis functions [10], as well as least squares fitting [1] or boundary elements [11]. While these methods allow for application of the level set method to any arbitrary structural mesh, they do not address situations where the working domain of the problem is non-rectangular. These situations arise frequently in engineering optimization, as it is often necessary to design a structure that must fit entirely within a non-rectangular region. In the following sections of this paper, we describe a new method for optimizing structures confined to such irregular domains. As in previous methods, the Hamilton–Jacobi equation is solved November 21, 2011
on a uniform rectilinear grid; however, the finite element modeling of the structure is performed using a fixed non-uniform structured mesh whose cells represent units in actual physical space. The structural sensitivities are calculated on the non-uniform mesh and then mapped into computational space using an isoparametric transformation within each cell. Because each element in the finite element mesh has its own unique cell in the rectilinear computational grid, we retain a oneto-one mapping between the structural sensitivities and advection velocities and avoid significant loss of robustness. This method is demonstrated in two ways: first, on a series of benchmark examples based on a variation of the classic two-dimensional L-bracket problem, and second, on a three-dimensional problem in which the working domain is given by the outer surface of an aircraft wing. We also introduce a novel approach to enforcing constraints based on the augmented Lagrangian formulation.
Figure 1: Sample structural design problem using level set parametrization. Here ΓD denotes the fixed boundary and ΓN denotes the boundary on which surface tractions are applied
min J Ω
subject to: c = 0 c˜ ≤ 0 Kd − F = 0 ( 1, ρ(xi ) = 10−3 ,
2. The Level Set Method
(1) xi ∈ Ω ¯ xi ∈ Ω
k(x) = ρ(x)k0 2.1. Problem Formulation Here, J is an arbitrary objective function that is dependent on the design variable, Ω. The structure may be subject to equality and inequality constraints, c and c˜ , and must satisfy the governing equation Kd − F = 0, where K is the global stiffness matrix, d is the vector of nodal displacements and F is the vector of applied forces. Elements located inside the solid region Ω are assigned a material density of unity, while those lying outside the boundary are given some small non-zero density, ρmin , so that they mimic the behaviour of a void space but can still be included into the global stiffness matrix without causing singularities. The densities of elements that are bisected by the structural boundary, ∂Ω, are interpolated based on the fraction of that element’s volume that lies inside Ω. The elasticity modulus, E, of a given element is the product of the element’s relative density and the elasticity of the element in the solid phase, E0 . The material properties are taken to be piecewise constant throughout each element.
The level set method was developed in the late 1980s by Osher and Sethian [12] as a method for tracking front propagation. It has been especially popular among researchers in fluid dynamics and computer vision, though an increasing number of authors have begun to apply the technique to structural topology optimization. Although similar to SIMP methods in some respects, the level set method differs in the way in which the structural shape and topology are parameterized. In both methods, the design domain is discretized into a series of finite elements whose relative material density, ρ, is allowed to vary continuously between 0 and 1 with the two extremes representing void and solid elements, respectively. In the SIMP approach, one optimizes these material densities directly, whereas in the level set method, one optimizes the location of the material boundary, ∂Ω, and from the location of this boundary, one determines which elements are solid and which are void, as shown in (1). Due to this boundary-based parameterization, the level set method avoids meshdependence and checkerboarding [3, 13] – two of the main numerical challenges associated with the SIMP method [14, 15].
Like the SIMP method, the level set approach begins by defining a bounded domain D ⊂ 0; x ∈ (D ∩ Ω). ψ
5 0
By implicitly representing the material boundary using a higher-order level set function, one allows for changes in topology, such as the merging of holes (in two dimensions) or cavities (in three dimensions) in the structure. The level set function is defined discretely at the nodes of a structured Cartesian mesh. These data are periodically interpolated in order to extract the precise location of the zero level set. While the level set method is mesh-independent, the converged solution of the optimization problem is dependent upon the initial topology [3, 13]. Therefore, one must carefully select the number of holes in the initial design according to the desired length scale for the final solution. Once the initial topology and material boundary are defined, the level set function is initialized as the signed distance of each point in the domain from the boundary, with negative distances used for points inside the solid region Ω and positive distances outside.
−5 −10 40 30 40 30
10
20 0
Y
10 0
X
Figure 2: Two-dimensional example of an optimized level set function, with contour lines displayed in the xy-plane
this way, the Hamilton–Jacobi equation replaces the optimization algorithms used in element-based topology optimization methods such as SIMP. The optimization reaches convergence once the advection velocities are within a small tolerance of zero along the material interface. As noted by Allaire et al. [3], the constrained optimization problem is generally not convex despite the absence of a penalization on intermediate densities. As a result, convergence to a global minimum cannot be guaranteed. Figure 2 shows an optimized level set function for a two-dimensional problem with the contour plots displayed in the xy-plane. If one takes the contour line corresponding to the surface ψ = 0 (the cyan-coloured contour) this plot will give the material boundary of the optimized structure. Because the Hamilton–Jacobi update assumes that the spatial derivative is valid over the distance vdt, it is necessary to periodically reinitialize the level set function to a form that validates this assumption. Therefore after every few updates, one must recover the signed distance form of the level set function, while maintaining the location of the current zero level set [3]. The results presented rely on an implementation of the fast sweeping method [17, 18, 19].
2.2. The Hamilton–Jacobi Equation In each optimizing iteration, the level set function is updated using the Hamilton–Jacobi equation (2). ∂ψ + v|∇ψ| = 0, (2) ∂t Rearrangement of this equation leads to the following scalar formula for updating the level set function at each point in the domain: ψt+dt = ψt − v|∇ψt |dt
60 50
20
(3)
Here the divergence ∇ψ is computed numerically at each point using finite differencing of the values of the level set function at adjacent nodes. The advection velocity v is given by the shape sensitivity of the objective function at each point. The variable t is a fictitious time parameter introduced to track the evolution of the level set function over the course of the optimization. The update described in Eq. (3) corresponds to moving the material boundary by a distance of vdt in the outward normal direction from the interface. The time step dt is typically chosen so that the Courant-Friedrichs-Lewy (CFL) condition is satisfied [16]. When the advection velocity is given by the shape sensitivities of the objective function, each Hamilton– Jacobi update is equivalent to a steepest descent step. In
2.3. Sensitivity Analysis For an arbitrary objective function, the shape sensitivity is defined according to the following equation, which uses the Fr´echet functional derivative [20], J 0 (Ω)(θ) = 3
Z vθ · nds, ∂Ω
(4)
where θ is an arbitrary small vector field, and the advection velocity v is equal to the local shape sensitivity of the objective J at a given point on the boundary. The shape sensitivity is usually some function of the displacement state u and can be derived using the following identity for generalized shape derivatives. For objectives of the form Z
J(Ω) =
Ω
j(x)dx +
L(Ω) = J(Ω) − λ · c(Ω),
where c is the vector of equality constraints satisfying c(Ω) = 0. This allows us to perform unconstrained optimization on the Lagrangian with respect to the design variable Ω as well as the vector of Lagrange multipliers λ. Convergence is reached once the first order KKT conditions are satisfied. The shape variable Ω is updated using the Hamilton– Jacobi equation, while the vector λ is updated using the following heuristic,
Z (5)
l(x)ds, ∂Ω
where Ω is a smooth, bounded open set, we have
J 0 (Ω)(θ) =
Z
+
Z
⇒ J 0 (Ω)(θ) = +
Ω
θ(x) · n(x)
Z∂Ω
λk+1 = λk + rc(Ω).
(6)
div(θ(x) j(x))dx
! ∂l(x) + H(x)l(x) ds ∂n(x)
θ(x) · n(x) j(x)ds
∂l(x) + H(x)l(x) ds, ∂n(x) !
θ(x) · n(x) ∂Ω
(9)
This update corresponds to a steepest descent step whose length is obtained by taking the first derivative of the Lagrangian with respect to the Lagrange multiplier. The step size r is a fixed value chosen to obtain the best trade-off between convergence time and stability. Consequently, this method is relatively inexpensive and works well for problems involving a single constraint.
(7)
∂Ω
Z
(8)
where H is the mean curvature of ∂Ω given by H = ∇n [3]. In the current study, all sensitivities are evaluated using the adjoint method derived by Allaire et al. [3].
3. The Isoparametric Formulation The shape sensitivities derived in Section 2.3 assume a smooth, continuous displacement state u, and strain field , which can be evaluated at any location throughout the working domain. In practice, the finite element method is used to generate a discrete representation of the displacement state, and the shape sensitivities are evaluated at the Gauss quadrature points for each element [21]. Alternatively the sensitivities can be averaged over the domain of the element, in which case the sensitivity is assumed to be piecewise constant across each element. Whereas most level set frameworks use square and cubic finite elements to coincide with the Cartesian computational grid, the isoparametric formulation allows for the use of arbitrary quadrilateral and hexahedral elements. Therefore, the isoparametric level set method relies heavily upon the isoparametric finite element formulation [24]. As with the standard approach, the Hamilton–Jacobi equation is solved on a uniform Cartesian grid with unit spacing between the nodes in all directions. This uniform, orthogonal spacing means that the spatial derivatives of the level set function can be evaluated simply by taking finite differences of the function evaluated at adjacent grid points. Each node in this computational grid corresponds to a four-node element (eight-node in three dimensions) in the finite
2.4. Constraint Handling Few papers address the topic of constrained optimization within the context of the level set method. Most authors account for constraints by adding a Lagrangianstyle constraint term to their objective. Some use a constant Lagrange multiplier [3], while others update the Lagrange multiplier using a bisection method [21]. The bisection approach is similar to the optimality criteria method, which has also been used in combination with SIMP frameworks [22]. However, like the optimality criteria method, the bisection approach is only applicable to problems involving a single constraint. Given this limitation, other researchers have proposed alternative methods for problems with multiple constraints [23], but no method has yet been widely adopted, as most lack flexibility, robustness or, both. Here, we introduce a new method based on an augmented Lagrangian approach, where the Lagrange multiplier is updated according to the Karush-Kuhn-Tucker (KKT) conditions. To apply equality constraints, we construct a Lagrangian function L by adding an additional term to the objective for each constraint to be handled, i.e., 4
Consider an infinitesimal patch located at some point p on the surface of the material interface. Since the orientation of the local coordinate axes is arbitrary, we can choose them so that the ζ axis points in the direction normal to the patch. Therefore, the area of the patch, measured in computational space, is given by dξdη. If the patch moves a distance dζ in the outward normal direction, it will trace a prism whose volume is given by dξdηdζ, when measured in computational space. This prism represents a total volume of dξdηdζ|J| in physical space. We now define a new quantity that is hereafter referred to as the relative impact. The relative impact of a shape perturbation at a point p on the material surface is defined as the product of the incremental volume Q caused by that perturbation, and the shape sensitivity of the objective function at the point p. Given a shape sensitivity, v p , computed in physical space at point p, we choose an equivalent computational sensitivity, vc , such that the relative impact of an infinitesimal shape perturbation at p, when measured in computational space, is maintained. Therefore we can deduce the following relationship:
element mesh (see Fig. 3). The shape sensitivities are evaluated in physical space using the non-uniform finite-element mesh (one sensitivity value for each element), and each of these values is transformed into its representative analog in computational space before being passed into the Hamilton–Jacobi equation. This transformation is based on the isoparametric transformation required to morph a unit square (or cube) into the quadrilateral (or hexahedral) shape of the element in question. The Jacobian for this transformation is defined in Eq. (12). Within each element one defines a local coordinate system ξ, η, ζ, where the orthogonal basis vectors ~ξ, ~η, ~ζ represent the directions in computational space. Figure 4 shows the two-dimensional local coordinate system associated with the quadrilateral element whose nodes are located at {xi , yi }, as expressed in physical space. Note that in the figure, the computational coordinates ξ, η range from −1 to 1. This choice of domain facilitates Gauss quadrature integration. This implies that each element occupies a volume of 4, as opposed to unity, in computational space. However, this quantity is inconsequential as long as all elements have the same volume in computational space. The physical location {x, y, z} of any point within the element can be expressed as a weighted sum of the node locations, x=
8 X
Ni (ξ, η, ζ)xi ,
where the weights, Ni , are the linear shape functions of the finite element: (1 + ξξi )(1 + ηηi )(1 + ζζi ) , ξi , ηi , ζi = ±1, 8 (11) which depend on the local coordinates {ξ, η, ζ}. Given this relationship between the physical coordinates and the local or computational coordinates, we can define the Jacobian matrix of the coordinate transformation within each element: Ni (ξ, η, ζ) =
∂x ∂ξ ∂y ∂ξ ∂z ∂ξ
∂x ∂η ∂y ∂η ∂z ∂η
∂x ∂ζ ∂y ∂ζ ∂z ∂ζ
.
(13)
⇒ dξdηdζvc = dξdηdζ|J|v p
(14)
⇒ vc = |J|v p .
(15)
Because the shape sensitivities are computed at discrete locations corresponding to the finite-element mesh, we approximate the Jacobian as being piecewise constant throughout each element. To do this we approximate |J| as being equal to the determinant of J averaged over the element volume, which turns out to be equal to the volume of the element itself, i.e., Z vole = det(J(ξ, η, ζ))dξdηdζ. (16)
(10)
i=1
J(ξ, η, ζ) =
Q c vc = Q p v p
Ωe
Therefore given some shape sensitivity, v p , evaluated in physical space, the resulting advection velocity in computational space is given by vc = vole v p . Thus, using the set of transformed sensitivities as the advection velocities in the Hamilton–Jacobi equation corresponds to taking a step in the direction of steepest descent.
(12)
3.1. The Constitutive Relation The constitutive matrix, which defines the relationship between the stress tensor, σ, and the strain tensor, ε, is given in Eq. (18) for the two-dimensional case
From here, we can derive the relationship between the sensitivities as calculated in physical space and the equivalent sensitivities to be used in computational space.
σ = Dε 5
(17)
Figure 3: Two-dimensional mapping of an arbitrary structured mesh to the uniform rectilinear computational mesh.
Figure 4: Two-dimensional illustration of the mapping from local to global coordinates for an arbitrary quadrilateral element
6
σ xx σyy σ xy
1 E = ν (1 − ν2 ) 0
ν 1 0
0 ε xx 0 εyy 1−ν ε xy 2
dividing the structure into equally-sized blocks of elements and assigning each block to a processor. (18)
3.2. Compliance Minimization Compliance minimization is commonly used among researchers of the level set method [3, 1, 16] to demonstrate new concepts, largely because of the ease with which it can be implemented. As shown by Allaire et al. [3], for the compliance objective function (i.e. the work done by external forces, Eq. (23)) the shape sensitivity reduces to the form shown in Eq. (24):
Here E and ν are the Young’s modulus and Poisson’s ratio of the isotropic material. Note that in all twodimensional examples, the plane stress condition is assumed. As is typical among previous studies, we use the ersatz material approach [3, 25], in which constitutive relation is tailored to account for those elements which lie, either partially or fully, outside the material domain. For these elements, the Young’s modulus is interpolated according to the portion of the element’s volume that falls inside the boundary, i.e., Ee = ρe E0 ⇒ De = ρe D0 ,
Jcomp (Ω) = =
1 vole
ρe ≈
1 vole
Z Ωe np X
Ω
g(x) · u(x)dx +
Z Ω
Z ΓN
f (x) · u(x)ds
εT (x)DεT (x) dx
(23)
(19) (20)
0 Jcomp (Ω)(θ) =
where, Ee is the effective Young’s modulus of the element, and E0 is the Young’s modulus of the solid material. The variable ρe represents the element volume fraction, which can also be interpreted as the relative material density. Under the isoparametric formulation, ρe is evaluated using the integral, ρe =
Z
h(−ψ(ξ, η))det(J(ξ, η))dξdη
(21)
h(−ψ p )|J p |.
(22)
Z Γ0
θ(x) · n(x) εT (x)DεT (x) ds
(24)
In the above equations, g and f denote body forces and surface tractions respectively, while u represents the continuous displacement field. The domain Γ0 refers to the traction-free surface, which is the only portion of the material boundary that varies during the optimization. Despite this restriction, the resulting advection velocity formula v = εT (x)DεT (x) is used for all points within the working domain when solving the Hamilton–Jacobi equation so that continuity is maintained. The compliance shape sensitivity, which is equal to twice the local strain energy density, is treated as being piecewise constant within each element. Therefore the advection velocity for each element is approximated as being proportional to the average strain energy density within that element. Under the isoparametric formulation, the calculation proceeds as follows for the twodimensional case: Z 1 Ue = εT (ξ, η)De ε(ξ, η)|J|dξdη, (25) 2 Ωe
p
The symbol h denotes the Heaviside function, which is equal to 1 wherever the level set function, ψ, is negative, and is equal to 0 everywhere else. In order to evaluate this integral numerically, the element is divided into a uniform, Cartesian grid of pixels, where the number of pixels, n p is equal to or greater than 1/ρmin . Each pixel is either solid or void depending on the value of the level set function at that pixel’s location, ψ p . In order to evaluate the level set function at these off-node locations, we simply interpolate using the linear shape functions shown in Eq. (11). The integral is then approximated as a weighted sum of solid pixels, with the weighting given by the determinant of the Jacobian matrix. Interpolating the level set function in this way can become computationally expensive, especially for threedimensional problems. This challenge can be mitigated significantly, however, through the use of parallel processing. We have vastly improved computation time by
where Ue is the total strain energy for element e. In order to obtain the average strain energy density, we divide this integral by the element volume given in Eq. (16). The expression can be simplified by introducing the strain-displacement matrix, B, which is comprised of the set of first derivatives of the element shape functions (Eq. (11)). We can therefore express the strain in terms of the strain displacement matrix and the vector of element nodal displacements de . 7
ε(ξ, η) = BT (ξ, η)de
This problem is solved using two different geometries. In the long L-bracket problem, the vertical and horizontal segments of the bracket have an aspect ratio of 2. In the short L-bracket problem, shown in Fig. 6, this aspect ratio is reduced to 1. This geometry causes the individual elements to have larger aspect ratios, particularly near the bottom-left corner of the domain. In both cases the applied force has a magnitude of 1.
(26)
Finally we introduce the element stiffness matrix, ke , which is defined as
k0 =
Z
ke =
Z
Ωe
Ωe
BT (ξ, η)D0 B(ξ, η)|J|dξdη
(27)
BT (ξ, η)De B(ξ, η)|J|dξdη
(28)
⇒ ke = ρe k0 ,
(29)
and the advection velocity becomes vp =
ρe dTe k0 de . vole
(30)
In order to obtain the computational velocity to be used in performing the Hamilton–Jacobi update, we multiply this expression by the element volume. Therefore the computational advection velocity reduces to
F vc =
ρe dTe k0 de .
(31)
In order to enforce a volume constraint we must find the shape sensitivity of the structural volume. From Eq. (6), we see that the shape sensitivity of the volume function is simply equal to unity. The equivalent sensitivity in computational space is simply the element volume. Combining the two functions according to the Lagrangian formulation defined above, the computational advection velocity for the compliance minimization problem subject to a constraint on structural volume is given by Eq. (32):
(a) Initial topology and loading conditions
90 80 70 60 50 40
vc = ρe dTe k0 de − λvole .
(32)
30 20
4. Numerical Examples 10
The isoparametric formulation derived above is demonstrated in a series of two-dimensional and threedimensional sample problems. The following examples are based on the classic L-bracket problem, which has been used frequently in topology optimization studies [5, 26]. In all L-bracket sample problems, the mesh is comprised of a lattice of trapezoidal elements as shown in Figure 5(b). This mesh is mapped to a single uniform rectangular mesh using a one-to-one mapping. The L-bracket problem is solved for minimum compliance subject to a constraint on the structural volume so that the total volume of the optimized structure does not exceed 40% of the volume of the working domain.
0
0
20
40
60
80
(b) Finite-element mesh
Figure 5: The long L-bracket problem
Figs. 7 and 8 show the optimized topologies for the long and short L-bracket problems. In each figure the zero contour of the level set function (i.e. ψ = 0) is plotted. In both cases, the isoparametric level set method produces converged, feasible solutions. The 8
F
Figure 7: Optimized topology for the short L-bracket problem (a) Initial topology and loading conditions 70
60
50
40
30
20
10
0
0
10
20
30
40
50
60
Figure 8: Optimized topology for the long L-bracket problem
70
(b) Finite-element mesh 2000
1
Figure 6: The short L-bracket problem
Compliance
zero contour of the level set function reveals straight, well-defined members and is unaffected by the lack of uniformity in the finite-element mesh. The convergence plots for both problems are shown in Figures 9 and 10. The plots indicate smooth convergence to a local minimum. Also, while the constraint function does oscillate during the optimization, the objective is reduced monotonically. Careful selection of the coefficient r in Eq. (9) can be used to reduce the amplitude and number of these oscillations. An effective choice is to use r ≈ 0.1λ0 , where λ0 is the initial value of the Lagrange multiplier. We now compare the solutions displayed in Figs. 7 and 8 with solutions obtained using the SIMP method. In both SIMP solutions, the finite-element mesh is iden-
1600
0.8
1200
0.6
800
0.4
400
0.2
0
0
20
40
60
80 100 120 Iteration Number
140
160
Volume Fraction
Compliance Volume Fraction
0 180
Figure 9: Convergence history for the long L-bracket compliance minimization problem
tical to that used above. We use a density filtering 9
500
1
400
0.8
300
0.6
200
0.4
100
0.2
Volume Fraction
Compliance
Compliance Volume Fraction
(a) Level set solution
(b) SIMP solution
Figure 12: Optimized density distribution for the short L-bracket 0
0
50
100 150 Iteration Number
200
0 250 800
Figure 10: Convergence history for the short L-bracket compliance minimization problem
level set method SIMP
700 600
Compliance
method to avoid checkerboarding [27, 28], in combination with the optimality criteria method [29] to carry out the constrained optimization problem. In the case of the short L-bracket, the isoparametric level set method generates a solution nearly identical to that of the SIMP method, as illustrated in Figure 12. In contrast, the two solutions to the long L-bracket problem differ greatly from one another, although both solutions exhibit similar compliance, as shown in Table 1.
500 400 300 200 100 0
0
50
100 150 Iteration Number
200
250
Figure 13: Comparison of the convergence histories of the SIMP and level set methods for the short L-bracket problem
(a) Level set solution
ture, where material is distributed throughout a circular domain, which contains a void region in the center. Fig. 14 shows the problem geometry, as well as the support and loading conditions. The optimized density distribution for the cantilevered ring problem is shown in Fig. 15. As is the case with the short L-bracket, the isoparametric level set method and the SIMP method yield very similar structures. Table 1 confirms that in this case, as well as in both L-bracket problems, the isoparametric level set solution exhibits very similar compliance to the SIMP solutions. Also, from the convergence plot given in Fig. 13, we see that the two methods have comparable computational cost in terms of the number of structural analyses required to reach convergence. Another problem that requires the use of a nonuniform mesh is the topology optimization of an aircraft wingbox. In the following example, the working domain of the problem is roughly determined by the shape of the wing’s outer surface, and material can be placed anywhere within this three-dimensional region.
(b) SIMP solution
Figure 11: Optimized density distribution for the long L-bracket
While the above examples make effective use of the isoparametric level set method, both L-bracket problems could have been solved using uniform Cartesian meshes. In addition to these problems, the isoparametric formulation enables us to apply the level set method to an entire class of problems that would otherwise be beyond the reach of most level set algorithms. We provide two examples of such problems. The first example is a novel problem designed specifically to test the effectiveness of the isoparametric formulation. This problem involves the optimization of a cantilevered ring struc10
level set method
SIMP
long L-bracket
364.44
366.18
short L-bracket
83.84
84.31
cantilevered ring
6.07
5.84
Table 1: Minimum compliance values for the SIMP and level set solutions. In all cases, the volume fraction is constrained at 0.4. (Note that the compliance values for the SIMP solutions have been calculated without any penalization.) (a) Initial topology 30
20
10
0
−10
Figure 16: Working domain and finite-element mesh for the wingbox optimization problem
−20
−30 −30
−20
−10
0
10
20
30
(b) Finite-element mesh
Figure 14: The two-dimensional cantilevered ring problem Figure 17: Loading conditions for the wingbox optimization problem
(a) Level set solution
The cross-section of the domain is based on a symmetric two-dimensional airfoil with the leading and trailing edges removed in order to replicate the shape of a wingbox. The resulting cross-section is then extruded to an aspect ratio of 3.02. The structure also has a taper ratio of 0.91, and a leading edge sweep angle of 9.4◦ . Fig. 16 shows the working domain and the finiteelement mesh for the problem, while Fig. 17 shows the loading conditions. The top surface of the structure is subject to a uniform distributed load acting in the downward direction. To the bottom surface we apply another distributed load, which is constant in the chordwise (x)
(b) SIMP solution
Figure 15: Optimized designs for the cantilevered ring problem
11
direction, and tapers off elliptically in the spanwise (z) direction. The structure is fixed along the face located at the root of the wing, where the wing connects with the fuselage of the aircraft. The structural mesh is comprised of 32 × 16 × 96 linear, eight-node hexahedral elements in the x, y, and z directions respectively. In this example, we use a fixed Lagrange multiplier and begin with an initial volume fraction of 60%. The Lagrange multiplier is set according to λ = 0.65C0 /V f0 , where C0 is the initial compliance and V f0 is the initial volume fraction. A minimum skin thickness is enforced on the top and bottom surfaces of the structure. At a given location along the span, the local skin thickness is equal to 1/2 the thickness of the local surface element measured in the y-direction. This thickness requirement is enforced by setting the minimum relative material densities of all top and bottom surface elements to ρmin = 0.5. The structure is optimized for minimum compliance, with a 40% volume constraint.
tom surfaces of the structure. The optimized structure can be interpreted as being comprised of two main types of component. A large portion of material is devoted to reinforcing the top and bottom skin of the structure, particularly near the root, where the fixed boundary condition is applied and where the internal bending moment of the structure is the largest. Also, as illustrated in Fig. 20, the remainder of the material is used to create a series of spar-like structures. These consist of truss-like members aligned at roughly 45◦ from the horizontal xzplane, in order to resist internal shear force in the z direction. Fig. 21 also reveals that these members are slanted slightly in the chordwise direction. This reflects the asymmetry of the structure due to the sweep angle. Because the wing is swept, any upward bending also results in a torsional moment. The chordwise slant of the spar members allows the structure to resist the twisting caused by this torsional moment.
Figure 20: Span-wise view of the optimized wing structure
Figure 18: Initial wing topology
Figure 21: Chord-wise view of the optimized wing structure
5. Conclusions In this paper we have presented a robust and computationally efficient method for performing structural topology optimization using the level set method combined with a non-uniform finite-element mesh. We have demonstrated how this isoparametric formulation can be used to effectively optimize structures that are confined to irregularly-shaped domains. In addition to the isoparametric formulation, we have also presented a
Figure 19: Optimized wing topology
Figs. 18 and 19 show the initial and optimized internal topologies for the wingbox problem. The figures do not include the wing skin, which covers the top and bot12
new method for enforcing equality constraints based on the augmented Lagrangian approach. These methods have been demonstrated on a series of two- and threedimensional structures optimized for minimum compliance subject to weight constraints. The results indicate that the method generates converged, feasible solutions that are comparable in performance to solutions obtained using the SIMP method. Optimization of structures modeled using body-fitted meshes represents a large and vitally important class of problems in structural optimization. The developments presented are aimed at extending the range of applicability of the level set method, thereby allowing researchers to use the method for tackling problems that more closely resemble those encountered in real-world applications. An important component of the push toward handling more realistic design problems is the implementation of additional objective and constraint functions. Although the current study focuses on compliance minimization subject to a volume constraint, the isoparametric method can also be applied to other objectives without any modification to the general algorithm. Going forward, some useful problems to explore include mass minimization subject to yield failure constraints — as is often done in the design of aerospace structures — as well as maximization of failure load with respect to elastoplastic collapse, as one might encounter in a structural engineering context.
[6] L. Krog, A. Tucker, Topology optimization of aircraft wing box ribs, in: 10th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference, Albany, New York, 2004, pp. AIAA 2004–4481. [7] T. H. Lee, S. Y. Han, J. K. Lim, Topology optimization of the inner reinforcement for an automobile hood using modal design sensitivity analysis, Key Eng. Mat. 183–187 (2000) 439–444. [8] S. Jin-fa, S. Jian-hui, Overview on innovation of topology optimization in vehicle cae, in: Electronic Computer Technology, 2009 International Conference on, 2009, pp. 457–460. [9] D.-F. Wang, B. Zhang, J. Chen, S.-P. Wang, S.-Q. Liu, Application of topology optimization to commercial vehicle cab design, in: International Conference on Measuring Technology and Mechatronics Automation, Vol. 2, Zhangjiajie, Hunan, China, 2009, pp. 763–766. [10] S.-H. Ha, S. Cho, Level set based topological shape optimization of geometrically nonlinear structures using unstructured mesh, Comput. Struct. 86 (13-14) (2008) 1447 – 1455. [11] K. Abe, S. Kazama, A boundary element approach for topology optimization problem using the level set method, Commum. Numer. Meth. Eng. 23 (5) (2007) 405–416. [12] S. Osher, J. A. Sethian, Fronts propagating with curvaturedependent speed: algorithms based on Hamilton–Jacobi formulations, J. Comput. Phys. 79 (1988) 12–49. [13] H. Jia, H. G. Beom, Y. Wang, S. Lin, B. Liu, Evolutionary level set method for structural topology optimization, Comput. Struct. 89 (2011) 445–454. [14] K. A. James, J. S. Hansen, J. R. R. A. Martins, Structural topology optimization for multiple load cases using a dynamic aggregation technique, Eng. Optimiz. 41 (2009) 1103–1118. [15] O. Sigmund, Numerical instabilities in topology optimization: A survey on procedures dealing with checkerboards, meshdependencies and local minima, Struct. Optimiz. 16 (1998) 68– 75. [16] S. Wang, M. Y. Wang, A moving superimposed finite element method for structural topology optimization, Int. J. Numer. Meth. Eng. 65 (2006) 1892–1922. [17] T. Oberhuber, Numerical recovery of the signed distance function, in: Czech-Japanese Seminar in Applied Mathematics, Prague, Czech Republic, 2004, pp. 148–164. [18] J. A. Sethian, A fast marching level set method for monotonically advancing fronts, in: Proc. Nat. Acad. Sci, 1995, pp. 1591– 1595. [19] Y. R. Tsai, L. Cheng, S. Osher, H. Zhao, Fast sweeping algorithms for a class of hamilton-jacobi equations, SIAM J. Numer. Anal. 41 (2003) 673–694. [20] J. Simon, Differentiation with respect to the domain in boundary value problems, Numer. Funct. Anal. and Optimiz. 2 (1980) 649–687. [21] S. Wang, K. Lim, B. Khoo, M. Wang, An extended level set method for shape and topology optimization, J. Comp. Phys. 221 (1) (2007) 395 – 421. [22] O. Sigmund, A 99 line topology optimization code written in matlab, Struct. Multidiscip. Optimiz. 21 (2001) 120–127. [23] M. Yulin, W. Xiaoming, A level set method for structural topology optimization with multi-constraints and multi-materials, Acta Mech. Sinica 20 (2004) 507–518. [24] B. M. Irons, Engineering applications of numerical integration in stiffness methods, AIAA J. 4 (11) (1966) 2035–2037. [25] Z. Luo, L. Tong, Z. Kang, A level set method for structural shape and topology optimization using radial basis functions, Comput. Struct. 87 (7-8) (2009) 425–434. [26] C. Le, J. Norato, T. Bruns, C. Ha, D. Tortorelli, Stress-based topology optimization for continua, Struct. Multidiscip. Optimiz. 41 (2010) 605–620.
Acknowledgements This research was supported by the Natural Sciences and Engineering Research Council of Canada and the Ontario Graduate Scholarship program.
References [1] P. Dunning, H. Kim, Investigation and improvement of sensitivity computation using the area-fraction weighted fixed grid fem and structural optimization, Finite Elem. Anal. Des. 47 (2011) 933–941. [2] Z. Liu, J. G. Korvink, Adaptive moving mesh level set method for structure topology optimization, Eng. Optimiz. 40 (6) (2008) 529–558. [3] G. Allaire, F. Jouve, A.-M. Toader, Structural optimization using sensitivity analysis and a level-set method, J. Comput. Phys. 194 (1) (2004) 363–393. [4] G. I. N. Rozvany, M. Zhou, T. Birker, Generalized shape optimization without homogenization, Struc. and Multidisc. Optim. 4 (1992) 250–252. [5] G. Allaire, F. Jouve, Minimum stress optimal design with the level set method, Eng Anal. Bound. Elem. 32 (11) (2008) 909– 918.
13
[27] T. E. Bruns, D. A. Tortorelli, Topology optimization of nonlinear elastic structures and compliant mechanisms, Comput. Method Appl. Mech. Eng. 190 (26-27) (2001) 3443–3459. [28] B. Bourdin, Filters in topology optimization, Int. J. Numer. Meth. Eng. 50 (2001) 2143–2158. [29] M. P. Bendsøe, Optimization of structural topology, shape and material, Springer, Berlin–Heidelberg–New York, 1995.
14
