We present a new method for the simultaneous topology optimization and material selection of structures made by the union of discrete geometric components, where each component is made of one of multiple available materials. Our approach is based on the geometry projection method, whereby an analytical description of the geometric components is smoothly mapped onto a density field on a fixed analysis grid. In addition to the parameters that dictate the dimensions, position, and orientation of the component, a size variable per available material is ascribed to each component. A size variable value of unity indicates that the component is made of the corresponding material. Moreover, all size variables can be zero, signifying the component is entirely removed from the design. We penalize intermediate values of the size variables via an aggregate constraint in the optimization. We also introduce a mutual material exclusion constraint that ensures that at most one material has a unity size variable in each geometric component. In addition to these constraints, we propose a novel aggregation scheme to perform the union of geometric components with dissimilar materials. These ingredients facilitate treatment of the multi-material case. Our formulation can be readily extended to any number of materials. We demonstrate our method with several numerical examples.

Introduction

Topology optimization has been used extensively in order to generate novel designs that improve structural performance while decreasing the cost. One important option to improve the structure is to employ multiple materials. By having an appropriate distribution of materials, multimaterial designs can outperform designs made of each of the materials separately. Moreover, besides decreasing cost, multimaterial designs can take advantage of substantial differences in properties to perform multiple functions.

Methods for multimaterial topology optimization were first introduced in density-based approaches. In [1], a method is developed to determine the optimal distribution of multiple phases to obtain composite materials with extreme thermal expansion behavior. This work employs an extension of the power-law interpolation used in solid isotropic material penalization (SIMP) [2,3] to three-phase material designs (two solid materials and void). The interpolation uses two design variables, one that indicates where to put material or void and another that determines which material to choose. This and similar formulations are applied to various problems, such as design of multiphase composites with extremal bulk modulus [4], design of multiphase piezoelectric actuators [5,6], combined optimization of material and voltage distribution [7], and optimal reinforcement of concrete structures [8]. As indicated in Ref. [9], this formulation violates the Hashin–Shtrikman bounds and renders different designs if the phases are interchanged. In Ref. [9], a formulation that combines the power-law penalization with an interpolation of every material property within its Hashin–Shtrikman bounds is proposed to circumvent the foregoing limitations to design multimaterial actuators. In all of these applications, the maximum number of phases involved is three (including two solid phases and a void phase) as the generalization of their interpolation schemes to more than three phases becomes quite involved.

A method for topology optimization of multi-material compliant mechanisms using an alternative material interpolation scheme that employs only one density variable is proposed in Ref. [10]. As opposed to the single-material topology optimization, where the optimal density indicates either void (zero) or solid (unity), in this method, the density variable in the optimal design can take any real value. The proximity of this value to one of the means of a sum of Gaussian distributions with given mean locations and with modes corresponding to the property values for different materials indicates the choice of material. To avoid premature convergence to undesired local minima, this method starts with a large standard deviation for each material's property distribution and decreases it gradually in order to sharpen the peaks.

The discrete material optimization method (DMO) [11] simultaneously optimizes the stacking sequence, reinforcement orientation, and choice of material of composite shell structures. In this method, there is a set of materials from which each finite element can be made in order to minimize the objective function. A material interpolation scheme is formulated, whereby an increase in the weighting factor for one material decreases the weighting factors for all other materials. The weighting factor for a given material indicates the relative influence of that material on the properties of the corresponding finite element. If the weighting factor for a given material is unity, it indicates that finite element is made solely of that material (e.g., pure phases). One challenge with this formulation is that the weighting factors do not add up to unity except for pure phases. The proposed solution is to normalize each weighting factor by the summation of all weighting factors. However, this alters the penalization effect and the monotonic convergence. The DMO method is used in Ref. [12] to optimize the buckling behavior of multimaterial composite shell structures, and in Ref. [13], to compare topology optimization of multimaterial structures with a mass constraint using both this method and the one presented in Ref. [1]. The work in Ref. [14] formulates simpler multimaterial interpolations by imposing one linear constraint per finite element that ensures that it is made of only one material (or no material).

Other density-based methods employ a heuristic rule for the multimaterial interpolation. The work in Ref. [15] uses a homogenization approach with a unit cell made of two materials and a square hole. The optimization allows for continuous mixtures of the three phases. Then, a postprocessing heuristic rule is applied to every finite element in the optimal design to determine if it is made of one of the two materials or void. In Ref. [16], material is iteratively removed from a fully solid initial design by applying a rule to each finite element that changes its material, or that makes it solid or void based on the elemental change in compliance incurred by these changes.

There also exist several level set methods for topology optimization of multimaterial structures. In the color level-set method [17], the boundaries of multiple regions made of different materials are given by the zero level-sets of multiple functions. Instead of using a level set function per material, this method employs m level-set functions whose combinations can represent up to $2m$ materials. This strategy circumvents the need for an equality constraint that would be required to avoid overlapping materials if one level set per material was used. Since the combination of level sets renders one and only one material at any point, this method eliminates the need for material interpolation. This method is applied to topology optimization of multimaterial compliant mechanisms in Ref. [18] and stress-based topology optimization of continuum structures involving multiple materials in Ref. [19]. In Ref. [20], the color level-set method is combined with a variational approach for optimization of structures made of functionally graded materials.

In the level set method of Ref. [21], only one level-set, piece-wise constant function is used as an “index” that indicates the choice of material. A constraint is applied to ensure that the level set function converges to the index values. This method is applied to design of piezoelectric actuators [22]. In Ref. [23], n level-set functions are used to represent n solid materials and void (i.e., n + 1 phases). This method employs a material interpolation similar to the one proposed in Ref. [1], except it uses the Heaviside of the level-set functions instead of density variables.

The design of multimaterial structures has also been studied with phase-field topology optimization methods. The generalized SIMP interpolation scheme proposed in Ref. [11] in combination with the mutual material exclusion constraint of Ref. [14] is employed in the phase field methods of Refs. [2426]. In Ref. [27], the volume fraction of each phase in phase field model directly represents the contribution of the corresponding material to the material properties. A constraint on these volume fractions ensures that they add up to unity, while a penalization term added to the objective function ensures that the volume fractions are zero or unity upon convergence. Phase field methods for topology optimization of multimaterial structures require a very large number of iterations, typically in the order of thousands.

All the aforementioned methods produce organic designs that cannot be readily manufactured with stock material. The geometry projection method presented in Refs. [2830] generates designs that can be made of stock material, such as bars and plates, by smoothly mapping an analytical description of the geometric elements onto a density field over a fixed finite element grid. A size variable is ascribed to each geometric component that is penalized in the spirit of SIMP, allowing the optimizer to entirely remove that component from the design. In the moving morphable components method of Refs. [31] and [32], a parametric description of bars is mapped onto a level set representation (termed the “topological description function”) for the analysis. The foregoing methods for design with discrete geometric elements use a single material. Recently, the work in Ref. [33] proposed a geometry projection method for the topology optimization of multimaterial, three-dimensional (3D) lattice structures. This method adapts the material interpolation scheme of Ref. [1] to the geometry projection framework, by assigning additional variables to each geometric component that determine the choice of material. Although possible, the extension of this interpolation scheme to more than two materials is not straightforward and it is more prone to getting locked into undesired local minima, as reported in Ref. [11]. Also recently, the work in Ref. [34] extends the moving morphable components method to accommodate multiple materials by seeding the initial design with geometric components made of different materials (i.e., the material of a component does not change during the optimization). The material interpolation scheme in regions where two or more geometric components intersect chooses the stiffer material of any intersecting component.

In this paper, we extend the geometry projection method so that it can be readily applied to structures whose geometric components can be made of one of any number of available materials. We achieve this by means of several key ingredients. The first ingredient corresponds to a new aggregation function used to perform the union of geometric components, since the function used in the previous geometry projection schemes cannot accommodate components made of dissimilar materials. A second important ingredient pertains to the interpolation of properties from the multiple material candidates. Here, we adapt the DMO formulation to accommodate discrete geometric components made of different materials. Unlike DMO, however, we impose two separate aggregate constraints in the optimization to (a) penalize the size variables of geometric components (discreteness constraint), and (b) to ensure that each component is made of at most one of the available materials (mutual material exclusion constraint). We demonstrate the effectiveness of our method via numerical examples.

The rest of the paper is organized as follows. Section 2 describes the projection of geometric components onto the analysis mesh for analysis, including the aggregation function for the union of geometric components and the interpolation of material properties. Sections 3 and 4 detail the constraints to penalize intermediate values of size variables and to ensure that geometric components are made of at most one material, respectively. The modification of the geometry projection to enforce a symmetric design is presented in Sec. 5. Section 6 describes the optimization problem. We present numerical examples to demonstrate our method in Sec. 7 and draw conclusions of our work in Sec. 8.

Geometry Projection

To perform the analysis for any given design made of the union of geometric components, we use the geometry projection method [28,29], wherein a parametric description of the components is smoothly mapped onto a density field over a fixed grid. In this work, we model bars as the solid enclosed by the surface defined by the set of points equidistant to a straight line segment (the medial axis of the bar). In two-dimensional (2D), these bars correspond to semicircular-ended rectangles of width w and out-of-plane thickness t [29], and in three-dimensional, they correspond to cylinders with semi-spherical ends, both with diameter w. Here, we follow closely [29] to define the projected density for a single component; the reader is referred to that work for more details.

The projection calculates a density at each point p in the design domain as the area fraction (for 2D components) or volume fraction (for 3D components) of the portion of the sample window $Bpr$ that intersects the solid structure ω (cf. Fig. 1), i.e.,
$ρ(x,r):=|Bpr∩ω|Bpr$
(1)
where
$Bpr:={x| ‖p−x‖≤r}$
(2)
Fig. 1
Fig. 1
Close modal
Assuming that the sampling window is small in relation to the geometric component, the projected density of Eq. (1) for component q can be approximated for 2D components as the area ratio of the circular segment of height $r−ϕq$ (cf. Fig. 1)
$ρq(ϕq,r):={0ifϕq>r1πr2[r2arccos(ϕqr)−ϕqr2−ϕq2 ]if−r≤ϕq≤r1ifϕq<−r$
(3)
and for 3D components as the volume ratio of the circular cap of height $r−ϕq$
$ρq(ϕq,r):={0ifϕq>r12+ϕq34r3−3ϕq4rif−r≤ϕq≤r1ifϕq<−r$
(4)

where $ϕq(p)$ is the signed distance from point $p$ to the boundary $∂ωq$ of geometric component q. That is, the projected density at a point $p$ with respect to a component q is entirely dictated by its signed distance to the component (for a fixed sampling window radius r).

For the bars considered in this work, a bar q is parameterized by the positions of the endpoints of its medial axis, $xqo$ and $xqf$. Thus, the signed distance is given by
$ϕq(dq,w):=dq(xqo,xqf,p)−w2$
(5)
where dq is the unsigned distance from point $p$ to the medial axis of bar q. As indicated by the arguments in the equation above, the minimum distance dq is a function of the bar parameters (cf. Fig. 1) and can be computed as
(6)
where
$a:=xqf−xqo$
(7)
$b:=p−xqo$
(8)
$e:=p−xqf$
(9)
$Pa⊥:=I−1||a||2a⊗a$
(10)
$g:=Pa⊥b$
(11)

In these expressions, $‖·‖$ denotes the Euclidean 2-norm, ⊗ denotes the dyadic (or tensor) product, and $Pa⊥$ is the perpendicular projector on $a$. We fix w as we desire to design structures with fixed-width members, but note that this could readily be a design parameter.

Since several components can intersect, in the previous work we aggregated the corresponding densities in the intersection by using a p-norm, which smoothly approximates the maximum projected density of any of the intersecting components. While this approach is effective for single-material structures, it does not accommodate components made of multiple materials. In this work, we consider a different aggregation strategy, which we explain in the following.

Let us first consider components made of the same material. We want the contribution of a bar q to the projected density at point $p$ to be 1.0 if $p$ is in the interior of the bar, and 0.0 if it is outside the bar. We can achieve this by using the Heaviside function of the signed distance to express the effective density as
$ρeff(z,p)=∑q=1NbH(−ϕq(z,p))ρq(z,p)∑q=1NbH(−ϕq(z,p))$
(12)
In the expression above, $z=[z1T z2T…zNbT]T$ is the vector of design variables, with $zq=[xq0T xqfT]T$ being the vector of design variables for bar q. Nb is the number of bars in the design. With the Heaviside weighting factor, a component contributes to the effective density when $ϕq≤0$, i.e., if $p$ is inside the component. The sum of all weighting factors in the denominator ensures $ρeff∈[0,1]$. Since the exact Heaviside function is not differentiable, we replace it with a smooth approximation $H̃ε$ [35] so that design sensitivities are well defined and we can use gradient-based optimizers
$H̃ε(x)={0 if x<−ε[12+x2ε+12πsin(πxε)]p if −ε≤x≤ε1 if x>ε$
(13)

We introduce the exponent p > 1 as a parameter to sharpen the Heaviside (i.e., to reduce the range over which it attains an intermediate value). In the case of single-material structures, the material properties are modified by some function of the foregoing effective density, as in our previous work. With the smooth Heaviside, it is clear that material properties vary continuously within a narrow band around the boundary of the bars.

To extend this interpolation to components made of one of Nm materials, we express the effective elasticity tensor at point $p$ as
$ℂ(z,p)=ℂmin+∑i=1Nm(ℂi−ℂmin) ρeffi(z,p)$
(14)
where $ℂi$ is the elasticity tensor for material i, $ℂmin$ is the elasticity tensor of a weak isotropic material that prevents the analysis from being ill-posed, and $ρeffi$ is an effective density per material given by
$ρeffi(z,p)=∑q=1NbH̃ε(−ϕq(z,p))ρq(z,p)wiq(z)∑q=1NbH̃ε(−ϕq(z,p))$
(15)
In this expression, $wiq$ is a weighting factor for bar q and material i. Here, we adapt the DMO interpolation scheme to our method, and define these weights as
$wiq(z)=(αiq)∏j=1Nm(1−(αj≠iq))$
(16)

The main difference between this expression and the original DMO approach for density-based topology optimization is in the variables α. In the density-based approach, these variables correspond to element-wise densities per material. In our approach, we ascribe a size variable $αiq∈{0,1}$ to geometric component q that indicates if it is made of material i when $αiq=1$. Correspondingly, the vector of design variables for bar q is now given by $zq=[xq0T xqfT αqT]T$ (note that $wiq$ only depends on $αq$ and not on $xq0T$ or $xqfT$, but we use the foregoing notation for simplicity).

These size variables are discrete, which precludes the use of nonlinear programming methods; therefore, we relax them so that they can take any value between 0.0 and 1.0. To ensure that the optimizer converges to a design with pure phases, we penalize intermediate values in the same spirit as the SIMP method. The second difference between our approach and the original DMO method lies in that the penalization is not imposed through the weights of Eq. (16), but through an optimization constraint, as described in Sec. 3. Furthermore, as we desire each component to be made of one and only one material, a mechanism to ensure that $αiq$ is 1.0 for at most one material is required. We detail a constraint to enforce this requirement in Sec. 4.

One problem with the expression for $ρeffi$ in Eq. (15) is that in the intersection between a solid bar made of material i (i.e., with $αiq=1$ and $αjq=0$ for $j≠i$) and a void bar (i.e., with $αiq=0 ∀i$), the effective density for material i is 0.5. Consequently, the effective elasticity tensor of Eq. (14) equals $0.5(ℂi+ℂmin)$, which is clearly incorrect. This occurs because the denominator is counting all the bars in the intersection and not only the solid bars. To remedy this, we modify the denominator so that it only counts the nonvoid bars (i.e., those for which $∑iNmαiq>0$) as follows:
$ρeffi(z,p)=∑q=1NbH̃ε(−ϕq(z,p))ρqwiq(z)∑q=1Nb(H̃ε(−ϕq(z,p))∑i=1Nmαiq)$
(17)
However, another problematic situation arises: if all bars in an intersection have a zero size variable for all materials, the denominator would be zero. In this case, we want the denominator to be replaced by 1, hence we redefine the effective density per material as
$ρeffi(z,p)=∑q=1NbH̃ε(−ϕq(z,p))ρqwiq(z)A+B$
(18)
wherein
$A=∑q=1Nb(H̃ε(−ϕq(z,p))∑i=1Nmαiq)$
(19)
$B=1−max(H̃ε(−ϕq(z,p))∑i=1Nmαiq)$
(20)
The term B equals 1 if all intersecting bars have zero size variables, and 0 if at least one size variable equals 1. We note that both situations occur only when all intersecting bars have pure phases. Finally, the maximum function is not differentiable, hence we replace it with a Kreisselmeier–Steinhauser (KS) approximation [36]
$ρeffi(z,p)=∑q=1NbH̃ε(−ϕq(z,p))ρqwiq(z)A+C$
(21)
where
$C=1−KSq(H̃ε(−ϕq(z,p))∑i=1Nmαiq)$
(22)
and
$KSi(x):=1kln(∑iekxi)$
(23)

The KS function approximates better the maximum as the parameter k increases.

Discreteness Constraint

Since the relaxed size variables $αiq$ can take any value between 0.0 and 1.0, a penalization scheme is needed to push them toward these values throughout the optimization in order to have a physical meaning. Unlike the DMO method, which performs the penalization through the weighting factors directly, we enforce this penalization via an equality constraint in the optimization
$gd(z):=4KSi,q(αiq(1−αiq))=0$
(24)
in which we use the lower-bound KS function
$KSi(x):=1kln(1N∑iekxi)$
(25)

where $αiq$ is the size variable for component q and material i. The lower-bound KS function approaches the maximum from below, therefore we circumvent the need to adaptively adjust the constraint limit (cf. [37]) and guarantee that $gd∈[0,1]$. As before, we denote gd being dependent on $z$ for notational simplicity, however it only depends on the size variables. If all size variables are 0.5, the above constraint will attain its largest value, 1.0. If all size variables are either zero or unity, the constraint value will be zero.

Since enforcing equality constraints is more difficult than enforcing inequality constraints in nonlinear programming methods, we replace the constraint of Eq. (24) with the inequality constraint
$gd(z)≤εd≪1$
(26)
For the penalization to be effective, the positive bound $εd$ should be small enough. However, if we start the optimization with a very small value, the size variables will quickly converge to 0.0 or 1.0 as the optimizer tries to reach the feasible region, which means that a geometric component can prematurely get “locked” into one of the available materials. Consequently, the optimization can get locked into undesirable local minima. To prevent this, we use a continuation strategy, whereby we start the optimization with a relatively large value of $εd(0)$. Then, we slowly start decreasing it at every iteration by an amount $Δεd$ once the relative change in the objective function in consecutive iterations falls below a specified value $Δf*$, i.e.,
$If Δf(I+1)≤Δf* then εd(I+1)←max(εd(I)−Δεd,εd*)$
(27)
where $εd*$ is the final constraint limit we want to attain, and the relative change in the objective function at iteration I + 1 is defined as
$Δf(I+1):=|f(I+1)−f(I)|f(I)$
(28)

Mutual Material Exclusion Constraint

By applying the discreteness constraint, we can ensure pure phases in the optimal design. However, this constraint does not prevent more than one material from having a size variable of 1.0 for a bar. To avoid this situation, we introduce a mutual material exclusion constraint in the optimization, defined as
$gm(z):=KSq(∑i=1Nmαiq)−1≤0$
(29)
in which we again use the lower-bound KS function of Eq. (25) with N = Nb.
When all the size variables $αiq$ for bar q satisfy the discreteness constraint, the term in parenthesis equals 1.0 if the bar is solid (in which case it is made of one and only material) or 0.0 if the bar is void. We employ a similar continuation strategy as in Sec. 3 to avoid premature convergence to an undesired local minimum. To this end, we write
$gm(z)≤εm≪1$
(30)
As before, we start with a relatively large value of $εm(0)$, and start decreasing it once the relative change in the objective function in consecutive iterations is smaller than a specified value, i.e.,
$If Δf(I+1)≤Δf* then εm(I+1)←max(εm(I)−Δεm,εm*)$
(31)

Note that $εm(0)$ must be smaller than 1.0, as otherwise it is possible for a bar to be made of more than one pure phase. In this situation, a decrease in one of the material size variables can increase the violation of the discreteness constraint, making it difficult for the optimizer to find feasible designs with bars made of at most one material.

Symmetry

If the boundary conditions and design envelope are such that a symmetric design is expected, density-based and level set topology optimization methods have significant design freedom that allows them to readily produce symmetric designs. However, as discussed in Ref. [29], the more restrictive design representation enforced by discrete geometric components can lead to asymmetric designs as the optimizer satisfies exactly the resource constraint. Specifically, if the design can only be made of a finite number of discrete geometric components, it is entirely possible to find an asymmetric design that exactly satisfies the weight fraction constraint and it has lower compliance than the best symmetric design that can be obtained with the available set of components.

In cases where a symmetric design is expected or desired, we employ a simple approach to enforce symmetry that consists of reflecting the geometry projection for elements on the opposite side of the symmetry plane. This approach was introduced in Ref. [33]. We define geometric components only on one side of the symmetry plane and bound their location to ensure they remain on that side. To compute the projected density for a point on the reflected side, we reflect the point with respect to the symmetry plane, and then compute the projected density for the reflected point as usual using Eq. (3). The sensitivities are modified accordingly to account for this reflection.

In the case of a symmetry line for two-dimensional problems, the reflected point is given by
$p̂:=Rp+2s$
(32)
$R:=[ cos(2ϕ) sin(2ϕ) sin(2ϕ)−cos(2ϕ)]$
(33)

where $s$ is a vector from the origin to the closest point on the line of symmetry, $ϕ$ is the angle between the symmetry line and the $e1$ axis, and $R$ is the reflection matrix. We use this approach to enforce symmetry in the example of Sec. 7.3.

Optimization Problem

In this work, we aim to minimize the structural compliance. Since we have multiple materials, a total volume constraint results in the exclusive use of the stiffest material. To avoid this, most of the works cited in Sec. 1 impose individual volume constraints for each one of the available materials. While this is an effective strategy, the determination of constraint limits for each material is somewhat arbitrary in practice. Instead, here we impose a weight fraction constraint. The space occupied by the structure and the design envelope are denoted by ω and Ω, respectively, with $ω⊂Ω$. We consider linearly elastic problems without a body load. The compliance minimization problem is then stated as
$minzf(u(z)):=∫Γtu(z)·tds$
(34)
$subject to$
$wf:=1γref|Ω|∑i=1Nmγi∫Ωρeffi(z,p)dv≤wf*$
(35)
$a(u(z),v)=l(v), ∀v∈UΩ,u∈UΩ$
(36)
$gd(z)≤εd(I)$
(37)
$gm(z)≤εm(I)$
(38)
$xq0,xqf∈Ω$
(39)
$0.0≤αiq≤1.0$
(40)
where γi is the physical density for material i, γref is a reference density that we choose to be 1 here, $|Ω|$ denotes the volume of the design region, $v$ denotes the virtual displacement, $t$ denotes the design-independent traction, and wf is the weight fraction. $UΩ:={u|u∈H1(Ω),u|Γu=0}$ is the set of admissible displacements, and $u$ is the displacement obtained from the solution to Eq. (36), with $a$ and $l$ being the energy bilinear and load linear forms, respectively, given by
$a(u,v):=∫Ω∇v·ℂ(z,p)∇udv$
(41)
$l(u,v):=∫Γtv·tds$
(42)
We impose move limits on the design variables at every iteration to improve convergence by normalizing the design variables by some reference value so that $0≤ẑ≤1$, where $ẑ$ denotes the scaled design variable. The coordinates of the end points $xq0$ and $xqf$ are normalized by the corresponding dimensions of the design region. The size variables $αiq$ do not need normalization, as they already lie within the desired magnitude range. A single move limit $0≤m≤1$ is imposed on all normalized variables as
$max(0,ẑ(I−1)−m)≤ẑ(I)≤min(1,ẑ(I−1)+m)$
(43)

We stop the optimization when the discreteness and mutual material exclusion constraints are feasible, and when the relative change of the objective function between consecutive iterations falls below the specified value $Δf*$.

Examples

To illustrate the effectiveness of our method, we present several examples. In all examples, we employ bilinear quadrilateral elements for the analysis. The entire code, including the finite element analysis, the sensitivities calculation, and the optimization, is implemented in matlab. To solve the optimization problem, we use the method of moving asymptotes [38,39] with the default parameters described in Ref. [39], i.e., $a0=1$ for the objective function, and al = 0, cl = 1000, and dl = 1 for every constraint l in the optimization (we refer the reader to [39] for a description of these parameters). Unless noted, the stopping criterion on the relative change in compliance between consecutive iterations is $Δf*=10−4$. Also, unless specified, we employ a move limit of m = 0.3. We use $εd0=1.0$ and $εd*=0.01$ for the discreteness constraint of Eq. (27), and $εm0=0.3$ and $εm0=0.01$ for the mutual material exclusion constraint of Eq. (31). We use a power of p = 2 in the smooth Heaviside approximation of Eq. (13). The constant k = 25 is used for the KS functions of Eqs. (23) and (25). All the materials considered are homogeneous, isotropic, and linearly elastic with Poisson's ratio $ν=0.3$, but with different Young's moduli and material densities.

Two-Bar Cantilever Beam.

The first example is a short cantilever beam made of two bar. The design envelope, boundary conditions, and initial design are shown in Fig. 2. The design envelope is meshed with a regular grid of 80 × 80 elements. There are two available materials with Young's moduli $E1=10$ and $E2=5$, and physical densities $γ1=0.9$ and $γ2=0.45$, respectively. The initial design, shown in the same figure, consists of two horizontal bars of width w = 0.25 and with $α1q=α2q=0.5$, q = 1, 2. For this problem, we use a looser stopping criterion on the relative change in compliance between consecutive iterations of $Δf*=10−3$.

Fig. 2
Fig. 2
Close modal

We perform the optimization for several weight-fraction limits $wf*$, and the results are presented in Table 1. The choice of $wf*$ for each run corresponds to expected configurations. For example, $wf*=0.0520$ corresponds to a design made of a single horizontal bar made of material 2 and completely inside the design envelope. Some of the weight-fraction limits account for the fact that half of the horizontal bar may be outside of the design envelope.

Table 1

Optimization results for two-bar cantilever beam problem. Red (×) indicates material 1 and green (○) indicates material 2. The last column indicates the number of iterations to convergence.

Design IDOptimal design$wf*$wfC$α11$$α21$$α12$$α22$Its.
10.05200.05061.00300.00001.00000.00100.0000191
20.10200.10190.03471.00000.00000.00090.9999164
30.11500.11500.01931.00000.00000.00080.999865
40.16000.15980.01631.00000.00000.00001.000035
50.18000.18000.01301.00000.00000.99960.000681
60.23000.22580.01031.00000.00001.00000.000036
Design IDOptimal design$wf*$wfC$α11$$α21$$α12$$α22$Its.
10.05200.05061.00300.00001.00000.00100.0000191
20.10200.10190.03471.00000.00000.00090.9999164
30.11500.11500.01931.00000.00000.00080.999865
40.16000.15980.01631.00000.00000.00001.000035
50.18000.18000.01301.00000.00000.99960.000681
60.23000.22580.01031.00000.00001.00000.000036

Two interesting cases are worth noting. For the second run, the weight-fraction limit $wf*=0.1020$ corresponds to the weight fraction of a single horizontal bar made of material 1 and completely inside the design envelope. However, the optimization produces a better design by using a short diagonal bar made of material 2. For the third run, the weight-fraction limit $wf*=0.115$ corresponds to a V-shape design made of the lighter material. However, the optimizer finds a better design whereby the horizontal bar is partially outside the design envelope.

In Fig. 3, we plot the compliance versus the weight fraction for the optimal designs of all six runs to illustrate that the compliance decreases as the weight fraction increases. As expected, for the lowest weight-fraction limit, the optimal design corresponds to a single bar made of the lighter (and weaker) material. As the weight-fraction constraint increases, the optimizer first introduces a second bar made of the heavier (and stiffer) material and eventually obtains a two-bar design made of the stiffest material. Finally, Fig. 4 shows the effective density for each of the two materials for the initial and optimal designs of run 4 in Table 1.

Fig. 3
Fig. 3
Close modal
Fig. 4
Fig. 4
Close modal

Messerschmitt–Bölkow–Blohm Beam.

The second example corresponds to the Messerschmitt–Bölkow–Blohm (MBB) beam widely studied in topology optimization. The design envelope, boundary conditions, and initial design are shown in Fig. 5. We note that, as we have discussed in the previous work (cf., [29]), geometry projection methods are more prone to converging to different local minima than free-form topology optimization methods due to the more restrictive design representation. The design envelope is meshed with a regular grid of 160 × 40 elements. We use the same two materials of the preceding example. The initial design consists of 21 bar of width w = 0.4 and with $α1q=α2q=0.5$, $q=1,…,21$. We consider two configurations of bars. In one configuration, as in the previous example, the endpoints $xq0$ and $xqf$ for each bar are independent from other bars, so that bars are “floating” inside the design envelope. In the second configuration, bars share common endpoints, so that the design remains connected at all times.

Fig. 5
Fig. 5
Close modal

We perform the optimization for several weight-fraction limits, $wf*=0.1,0.11,…,0.19$. The optimal designs for the floating and connected configurations are shown in Figs. 6 and 7, respectively. Several of the designs obtained resemble known solutions for the MBB beam. As expected, the runs with floating bars produce better designs than their connected bars counterparts since they have more design freedom. Also, the smallest “solid” size variable $αiq$ is 0.988249 for all floating bars runs and 0.983790 for all connected bars runs; and the largest “void” size variable $αiq$ is 0.036276 for all floating bars runs and 0.038084 for all connected bar runs. This is an indication that the discreteness constraints are very effective in penalizing intermediate values of the size variables.

Fig. 6
Fig. 6
Close modal
Fig. 7
Fig. 7
Close modal

Michell Cantilever.

We now present an example that shows that the proposed formulation can be readily extended to any number of materials. It corresponds to a cantilever frame considered in Michell's work [40], with design envelope and boundary conditions as shown in Fig. 8. The design envelope is meshed with 9700 elements. We enforce symmetry of the design with respect to the horizontal center line shown in Fig. 8. For this example, we use a tighter move limit of m = 0.2 to improve convergence. The initial design is made of 12 near-zero length bars of width w = 0.25 and with $αiq=0.5, q=1,…,12$.

Fig. 8
Fig. 8
Close modal

We perform the optimization using four materials with Young's moduli $E1=6.5, E2=5.0, E3=4.5$, and $E4=3.5$ and physical densities $γ1=0.55, γ2=0.4, γ1=0.35$, and $γ2=0.25$. We use a weight-fraction limit of $wf*=0.028$ that allows us to obtain a design with four materials, as otherwise we would get designs made only of the stiffer material(s) if $wf*$ is larger, or made only of the weaker material(s) if $wf*$ is smaller. Figure 9 shows the optimal designs using one, two, three, and four materials. As expected, the designs improve as we increase the number of materials available. Also, the smallest solid size variable $αiq$ is 0.999817 for all runs, and the largest void size variable $αiq$ is 0.003536 for all runs, once again indicating the effectiveness of the penalization scheme.

Fig. 9
Fig. 9
Close modal

Three-Dimensional Cantilever Beam.

In the last example, we perform the optimization for a 3D cantilever beam. For this example, we migrated our matlab code to C++ using the deal.II library [41,42]. The design envelope, boundary conditions, and initial design are shown in Fig. 10. The design envelope is meshed with a regular grid of $64×32×32$ elements. We use the same two materials of the first example. The initial design consists of 42 bar of width w = 0.2 and with $α1q=α2q=0.5, q=1,…,42$. We consider an initial design with connected bars as in the previous MBB beam example of Sec. 7.2. The optimal designs corresponding to weight-fraction limits of $wf*=$ 0.01, 0.02, 0.03, and 0.04 are shown in Table 2.

Fig. 10
Fig. 10
Close modal
Table 2

Optimization results for 3D cantilever beam problem. Red (dark) indicates material 1 and green (light) indicates material 2.

C$wf*$IsoSideTopFront
3.89141%
2.524622%
1.407823%
1.105114%
C$wf*$IsoSideTopFront
3.89141%
2.524622%
1.407823%
1.105114%

Conclusions

We presented a method for the design via topology optimization of structures constructed as the union of geometric components, where each component is made of one of several available materials or removed from the design. Several examples that minimize the structural compliance subject to a constraint on the weight fraction demonstrate the proposed method. The available materials have different moduli but also different physical densities, hence a combination of materials is most advantageous for some weight fraction limits.

The examples demonstrate our method's effectiveness in producing structurally efficient multi-material designs. By penalizing intermediate size variables and enforcing the mutual material exclusion requirement as constraints in the optimization and not through the interpolation scheme, our technique makes it easier to incorporate any number of materials. Unlike density- and level set-based topology optimization methods for design with multiple materials, which produce material phases with geometries that are difficult to manufacture and assemble, the use of geometric components that are readily fabricated makes it easier to physically realize the multimaterial structures. Also, instead of imposing arbitrary volume fraction constraints on each of the available materials, we directly constrain the weight, which is a more natural design requirement.

We have demonstrated our method for two- and three-dimensional bars modeled with offset surfaces. Since the computation of the effective density ρeff of Eq. (21) is uncoupled from the calculation of the projected density ρq of Eqs. (3) and (4), the scheme to interpolate the material properties from various discrete geometric components should work for any geometry representation, and thus, we believe that our method can be applied to other component geometries; this will be demonstrated in future work. Future work will also focus on incorporating other important structural criteria.

Acknowledgment

Support from the National Science Foundation, award CMMI-1634563 to conduct this work is gratefully acknowledged. We also thank Professor Krister Svanberg for kindly providing his MMA Matlab optimizer to perform the optimization.

Funding Data

• National Science Foundation (Grant No. CMMI-1634563).

References

1.
Sigmund
,
O.
, and
Torquato
,
S.
,
1997
, “
Design of Materials With Extreme Thermal Expansion Using a Three-Phase Topology Optimization Method
,”
J. Mech. Phys. Solids
,
45
(
6
), pp.
1037
1067
.
2.
Bendsøe
,
M. P.
,
1989
, “
Optimal Shape Design as a Material Distribution Problem
,”
Struct. Multidiscip. Optim.
,
1
(
4
), pp.
193
202
.
3.
Rozvany
,
G.
, and
Zhou
,
M.
,
1991
, “
The COC Algorithm, Part I: Cross-Section Optimization or Sizing
,”
Comput. Methods Appl. Mech. Eng.
,
89
(
1–3
), pp.
281
308
.
4.
Gibiansky
,
L. V.
, and
Sigmund
,
O.
,
2000
, “
Multiphase Composites With Extremal Bulk Modulus
,”
J. Mech. Phys. Solids
,
48
(
3
), pp.
461
498
.
5.
Carbonari
,
R. C.
,
Silva
,
E. C.
, and
Nishiwaki
,
S.
,
2007
, “
Optimum Placement of Piezoelectric Material in Piezoactuator Design
,”
Smart Mater. Struct.
,
16
(
1
), p.
207
.
6.
Luo
,
Z.
,
Gao
,
W.
, and
Song
,
C.
,
2010
, “
Design of Multi-Phase Piezoelectric Actuators
,”
J. Intell. Mater. Syst. Struct.
,
21
(
18
), pp.
1851
1865
.
7.
Kang
,
Z.
,
Wang
,
R.
, and
Tong
,
L.
,
2011
, “
Combined Optimization of Bi-Material Structural Layout and Voltage Distribution for in-Plane Piezoelectric Actuation
,”
Comput. Methods Appl. Mech. Eng.
,
200
(
13–16
), pp.
1467
1478
.
8.
Luo
,
Y.
, and
Kang
,
Z.
,
2013
, “
Layout Design of Reinforced Concrete Structures Using Two-Material Topology Optimization With Drucker–Prager Yield Constraints
,”
Struct. Multidiscip. Optim.
,
47
(
1
), pp.
95
110
.
9.
Sigmund
,
O.
,
2001
, “
Design of Multiphysics Actuators Using Topology Optimization–Part II: Two-Material Structures
,”
Comput. Methods Appl. Mech. Eng.
,
190
(
49–50
), pp.
6605
6627
.
10.
Yin
,
L.
, and
Ananthasuresh
,
G.
,
2001
, “
Topology Optimization of Compliant Mechanisms With Multiple Materials Using a Peak Function Material Interpolation Scheme
,”
Struct. Multidiscip. Optim.
,
23
(
1
), pp.
49
62
.
11.
Stegmann
,
J.
, and
Lund
,
E.
,
2005
, “
Discrete Material Optimization of General Composite Shell Structures
,”
Int. J. Numer. Methods Eng.
,
62
(
14
), pp.
2009
2027
.
12.
Lund
,
E.
,
2009
, “
Buckling Topology Optimization of Laminated Multi-Material Composite Shell Structures
,”
Compos. Struct.
,
91
(
2
), pp.
158
167
.
13.
Gao
,
T.
, and
Zhang
,
W.
,
2011
, “
A Mass Constraint Formulation for Structural Topology Optimization With Multiphase Materials
,”
Int. J. Numer. Methods Eng.
,
88
(
8
), pp.
774
796
.
14.
Hvejsel
,
C. F.
, and
Lund
,
E.
,
2011
, “
Material Interpolation Schemes for Unified Topology and Multi-Material Optimization
,”
Struct. Multidiscip. Optim.
,
43
(
6
), pp.
811
825
.
15.
Buehler
,
M. J.
,
Bettig
,
B.
, and
Parker
,
G. G.
,
2004
, “
Topology Optimization of Smart Structures Using a Homogenization Approach
,”
J. Intell. Mater. Syst. Struct.
,
15
(
8
), pp.
655
667
.
16.
Mirzendehdel
,
A. M.
, and
Suresh
,
K.
,
2015
, “
A Pareto-Optimal Approach to Multimaterial Topology Optimization
,”
ASME J. Mech. Des.
,
137
(
10
), p.
101701
.
17.
Wang
,
M. Y.
, and
Wang
,
X.
,
2004
, “
Color” Level Sets: A Multi-Phase Method for Structural Topology Optimization With Multiple Materials
,”
Comput. Methods Appl. Mech. Eng.
,
193
(
6–8
), pp.
469
496
.
18.
Wang
,
M. Y.
,
Chen
,
S.
,
Wang
,
X.
, and
Mei
,
Y.
,
2005
, “
Design of Multimaterial Compliant Mechanisms Using Level-Set Methods
,”
ASME J. Mech. Des.
,
127
(
5
), pp.
941
956
.
19.
Guo
,
X.
,
Zhang
,
W.
, and
Zhong
,
W.
,
2014
, “
Stress-Related Topology Optimization of Continuum Structures Involving Multi-Phase Materials
,”
Comput. Methods Appl. Mech. Eng.
,
268
, pp.
632
655
.
20.
Wang
,
M. Y.
, and
Wang
,
X.
,
2005
, “
A Level-Set Based Variational Method for Design and Optimization of Heterogeneous Objects
,”
Comput.-Aided Des.
,
37
(
3
), pp.
321
337
.
21.
Wei
,
P.
, and
Wang
,
M. Y.
,
2009
, “
Piecewise Constant Level Set Method for Structural Topology Optimization
,”
Int. J. Numer. Methods Eng.
,
78
(
4
), pp.
379
402
.
22.
Luo
,
Z.
,
Tong
,
L.
,
Luo
,
J.
,
Wei
,
P.
, and
Wang
,
M. Y.
,
2009
, “
Design of Piezoelectric Actuators Using a Multiphase Level Set Method of Piecewise Constants
,”
J. Comput. Phys.
,
228
(
7
), pp.
2643
2659
.
23.
Wang
,
Y.
,
Luo
,
Z.
,
Kang
,
Z.
, and
Zhang
,
N.
,
2015
, “
A Multi-Material Level Set-Based Topology and Shape Optimization Method
,”
Comput. Methods Appl. Mech. Eng.
,
283
, pp.
1570
1586
.
24.
Zhou
,
S.
, and
Wang
,
M. Y.
,
2007
, “
Multimaterial Structural Topology Optimization With a Generalized Cahn–Hilliard Model of Multiphase Transition
,”
Struct. Multidiscip. Optim.
,
33
(
2
), pp.
89
111
.
25.
Tavakoli
,
R.
, and
Mohseni
,
S. M.
,
2014
, “
Alternating Active-Phase Algorithm for Multimaterial Topology Optimization Problems: A 115-Line Matlab Implementation
,”
Struct. Multidiscip. Optim.
,
49
(
4
), pp.
621
642
.
26.
Tavakoli
,
R.
,
2014
, “
Multimaterial Topology Optimization by Volume Constrained Allen–Cahn System and Regularized Projected Steepest Descent Method
,”
Comput. Methods Appl. Mech. Eng.
,
276
, pp.
534
565
.
27.
Wallin
,
M.
,
Ivarsson
,
N.
, and
Ristinmaa
,
M.
,
2015
, “
Large Strain Phase-Field-Based Multi-Material Topology Optimization
,”
Int. J. Numer. Methods Eng.
,
104
(
9
), pp.
887
904
.
28.
Bell
,
B.
,
Norato
,
J.
, and
Tortorelli
,
D.
,
2012
, “
A Geometry Projection Method for Continuum-Based Topology Optimization of Structures
,”
AIAA
Paper No. AIAA 2012-5485..
29.
Norato
,
J.
,
Bell
,
B.
, and
Tortorelli
,
D.
,
2015
, “
A Geometry Projection Method for Continuum-Based Topology Optimization With Discrete Elements
,”
Comput. Methods Appl. Mech. Eng.
,
293
, pp.
306
327
.
30.
Zhang
,
S.
,
Norato
,
J. A.
,
Gain
,
A. L.
, and
Lyu
,
N.
,
2016
, “
A Geometry Projection Method for the Topology Optimization of Plate Structures
,”
Struct. Multidiscip. Optim.
,
54
(
5
), pp.
1173
1190
.
31.
Guo
,
X.
,
Zhang
,
W.
, and
Zhong
,
W.
,
2014
, “
Doing Topology Optimization Explicitly and Geometrically—a New Moving Morphable Components Based Framework
,”
ASME J. Appl. Mech.
,
81
(
8
), p.
081009
.
32.
Deng
,
J.
, and
Chen
,
W.
,
2016
, “
Design for Structural Flexibility Using Connected Morphable Components Based Topology Optimization
,”
Sci. China Technol. Sci.
,
59
(
6
), pp.
839
851
.
33.
Watts
,
S.
, and
Tortorelli
,
D. A.
,
2017
, “
A Geometric Projection Method for Designing Three-Dimensional Open Lattices With Inverse Homogenization
,”
Int. J. Numer. Methods Eng.
,
112
(
11
), pp.
1564
1588
.
34.
Zhang
,
W.
,
Song
,
J.
,
Zhou
,
J.
,
Du
,
Z.
,
Zhu
,
Y.
,
Sun
,
Z.
, and
Guo
,
X.
, “
Topology Optimization With Multiple Materials Via Moving Morphable Component (MMC) Method
,”
Int. J. Numer. Methods Eng.
,
113
(11), 1653–1675.
35.
Faure
,
A.
,
Michailidis
,
G.
,
Parry
,
G.
,
Vermaak
,
N.
, and
Estevez
,
R.
,
2017
, “
Design of Thermoelastic Multi-Material Structures With Graded Interfaces Using Topology Optimization
,”
Struct. Multidiscip. Optim.
,
56
(
4
), pp.
823
837
.
36.
Kreisselmeier
,
G.
,
1979
, “
Systematic Control Design by Optimizing a Vector Performance Index
,” IFAC Symposium, Zürich, Switzerland, Aug. 29–31, pp. 113–117.
37.
Le
,
C.
,
Norato
,
J.
,
Bruns
,
T.
,
Ha
,
C.
, and
Tortorelli
,
D.
,
2010
, “
Stress-Based Topology Optimization for Continua
,”
Struct. Multidiscip. Optim.
,
41
(
4
), pp.
605
620
.
38.
Svanberg
,
K.
,
2002
, “
A Class of Globally Convergent Optimization Methods Based on Conservative Convex Separable Approximations
,”
SIAM J. Optim.
,
12
(
2
), pp.
555
573
.
39.
Svanberg
,
K.
,
2007
, “
MMA and GCMMA, Versions September 2007
,” Optimization and Systems Theory, KTH, Stockholm, Sweden, accessed July 12, 2018, https://people.kth.se/~krille/mmagcmma.pdf
40.
Michell
,
A. G. M.
,
1904
, “
LVIII. the Limits of Economy of Material in Frame-Structures
,”
Philos. Mag. J. Sci.
,
8
(
47
), pp.
589
597
.
41.
Bangerth
,
W.
,
Hartmann
,
R.
, and
Kanschat
,
G.
,
2007
, “
Deal.II—A General Purpose Object Oriented Finite Element Library
,”
ACM Trans. Math. Software
,
33
(
4
), pp.
24/1
24/27
.
42.
Bangerth
,
W.
,
Davydov
,
D.
,
Heister
,
T.
,
Heltai
,
L.
,
Kanschat
,
G.
,
Kronbichler
,
M.
,
Maier
,
M.
,
Turcksin
,
B.
, and
Wells
,
D.
, “
The Deal.II Library, Version 8.4
,”
J. Numer. Math.
,
24
(3), 135–141.