Materials/structures with well-designed microstructures showing excellent properties are ubiquitous in nature, e.g., both high stiffness and toughness of the nacre [1], remarkable bending stiffness of bamboo [2], fascinating colors of the butterfly wings [3]. Furthermore, an increasing attention is seen in design of the metamaterials possessing extraordinary properties that are not commonly found in nature [4][5]. By taking an advantage of the modern manufacturing technology, complex designs of microstructures can be fabricated conveniently [6][7][8].
Inspired by the inverse homogenization approach [9], many works devoted to design microstructures to tailor properties of architected materials [6][7][8][9][10]. Taking this to the next step, multiscale topology optimization has been developed to simultaneously design a macroscopic structure and the associated material microstructure(s). Rodrigues et al. [11] obtained a hierarchical design by optimizing the periodic microstructures in every macroscale finite element. This method, however, would lead to a very high computational cost particularly for three dimensional cases [12]. As a more computational tractable strategy, Liu et al. [13] proposed the Porous Anisotropic Material with Penalization (PAMP) model to obtain a multiscale structure with a uniform microstructure. This approach has been extended to consider dynamic and thermomechanical effects [14,15]. Sivapuram et al. [16] recently proposed a more generalized simultaneous structure and material optimization formulation where any number of microstructures can be obtained. They proposed a linearization formulation to decompose the macroscale and microscale optimizations, thereby parallel and distributed computing can be adopted easily.
In majority of the multiscale optimization studies, the asymptotic homogenization theory has been used to obtain the effective material property of a periodic microstructure. The assumption of the scale separation and the periodicity in the homogenization theory [17] ignores the connectivity of the adjacent microstructures [11,12,16]. It should be noted that, besides poor manufacturability, more importantly, poor microstructural connectivity would lead to load transition issue as well as the deviation of effective property estimated by homogenization method. In order to obtain the optimal sizes of microstructures, consider boundary effect as well as guarantee the connectivity between different microstructures, Alexandersen and Lazarov [18] abandoned the homogenization theory and directly optimized the micro-structures with an extremely fine mesh. Even though the authors tried to reduce the computational resource requirements, the computation cost is substantially higher than the homogenization-based approaches. It has been shown that the optimal solution of a unit cell converges rapidly to that obtained by inverse homogenization as the number of repetitive cells increases (beyond five or six in the case of mechanical properties) [19,20]. In addition, for boundary effect, as shown in [21], the thickness of boundary layer has the same scale of the unit cell. Results in [19][20][21] suggested that the homogenization method can efficiently offer a reasonable approximation for a large number of repeated unit cells, when some global measures, e.g., mean compliance, are taken into consideration. Liu et al. [22] divided the structural domain into several subdomains and boundary layers. The optimum subdomain periodic microstructures were obtained by inverse homogenization and direct optimization was applied to the boundary layers with an extremely fine mesh to obtain the smooth transitions between different microstructures.
For microstructures with graded properties, the connectivity has been enforced implicitly via fixing some connective elements or applying a pseudo load or adding nonlinear diffusion term to the objective function, [23][24][25]. Whilst such implicit controls have been shown to be effective in many cases, they cannot guarantee the connectivity (Fig. 1(a)). It may also have an effect of over-constraining the design space.
Another approach for connecting microstructures is to apply a post-processing based on the metamorphosis technique. Wang et al. [26] generated a series of selfsimilar and connected microstructures by interpolating between a prototype cell and a solid cell. Such an interpolation technique is a well-established practice for the level set method in the field of image processing. However, it can create discontinuous member sizes such as shown in Fig. 1(b) which can lead to stress concentrations. Furthermore, such an interpolation method may lead to a floating topology that is no longer physical when applied to two cells with distinct topologies, e.g. Fig. 1(c). [16] is employed for numerical demonstrations and the level set topology optimization method of Dunning and Kim [27] is used for topology optimization at both macro and microscale.
The remainder of the paper is organized as follows. In Sections 2 and 3, the level set topology optimization method and the multiscale optimization formulations are outlined for completeness. To connect microstructures in multiscale design problems illustrated in Section 3, an explicit CI is proposed in Section 4 and applied to control the microstructural connectivity in Section 5. The CI-constrained optimization approach is applied to obtain optimal multiscale designs with well-connected microstructures in Section 6, followed by some concluding remarks.
The level set topology optimization method used in this paper follows Dunning and Kim [27] and is briefly summarized in this section.
In the level set method, a structure is implicitly represented by a level function 𝜙𝜙(𝒙𝒙), (1) and the Hamilton-Jacobi equation expressed by (2) advances the structural interface or boundary,
where D is the design domain, Ω is the structural domain
where 𝑘𝑘 is the iteration number, Δ𝜕𝜕 is the time step, 𝑉𝑉 𝑛𝑛,𝑖𝑖 is the normal velocity of the boundary point 𝑖𝑖 and 𝑍𝑍 𝑛𝑛,𝑖𝑖 is the distance of the associated normal movement.
A general topology optimization problem can be written as min
where 𝑓𝑓 is the objective functional and 𝑔𝑔 𝑗𝑗 denotes a constraint functional. With the help of shape derivative, linearizing (4) gives,
In the multiscale topology optimization framework [16], the design domain D is first divided to a certain number of subregions (i.e., D 1 , … , D 𝑁𝑁 with ⋃ D 𝑒𝑒 𝑁𝑁 𝑒𝑒=1
and ∀ D 𝑒𝑒 ⋂D 𝑙𝑙 = ∅, 𝑒𝑒 ≠ 𝑙𝑙; 𝑒𝑒, 𝑙𝑙 = 1, … , 𝑁𝑁). An illustrative example is shown in Fig. 2.
It is also assumed that uniform microstructures are distributed in each subregion and can be analyzed by the asymptotic homogenization theory. By simultaneously optimizing the macroscopic structure and microscopic unit cells, the design space is greatly extended to improve the functional performance.
For the case of the compliance minimization, the multiscale topology optimization problem is formulated as follows
The macroscale structure
where 𝜙𝜙 𝑀𝑀 is the macroscopic level set function defined on domain D; 𝜙𝜙 𝑚𝑚𝑒𝑒 is the 𝑒𝑒- For the multiscale topology optimization problems, Sivapuram et al. [16] suggested to linearize first both the objective function and constraints with the use of their
with 𝜀𝜀 𝑝𝑝𝑝𝑝 0 denoting the unit strain tensors and 𝒚𝒚 𝐵𝐵 denoting associated boundary point.
With use of the optimal boundary movements, the macroscopic and microscopic structures can be obtained via (3). It should be noted that, on one hand, the optimal microstructures are affected by the macroscopic strain through (11); on the other hand, the macroscopic strain is determined by the macroscopic structure and the effective properties of microstructures. The optimal macrostructure and microstructures are inherently coupled [29]. With such coupling being ignored, the effective properties of the optimal microstructures may not be consistent in optimality with respect to the current macroscopic structure at intermediate iterations. However, as discussed in [29],
it was observed that the inconsistency would vanish as a solution converges and has little effects on the final solution.
The L-beam example in Fig. 3 The obtained macroscale structure is very similar as the single scale design result [30],
which is reasonable. The three optimized microstructures, however, are not wellconnected. In addition, they cannot transfer loads well in reality as expected by the homogenization theory. The multiscale structures are analyzed in the structural scale using the effective properties of microscopic unit cells and the scale separation assumption makes the structure unware of microstructural topologies. As a result, the effects of the connectivity between spatially varying microstructures cannot be considered in analysis hence, cannot be considered in design optimization.
In order to address this challenge in the multiscale design optimization framework, we introduce a quantity to measure the microstructural connectivity. Connectivity Index (CI) is defined using a local interface region of the mutual boundary, as shown Fig. 5. Considering the two strips from each unit cell, the following normalized parameter, CI can be defined:
where Y 𝑖𝑖 s is the strip region colored in light red near the adjacent cell of 𝑖𝑖 -th microstructure. 𝐑𝐑 denotes the mirror reflection operator, which maps a point 𝒚𝒚 to its symmetric counterpart 𝐑𝐑𝒚𝒚 in the adjacent cell. CI is actually a symmetry measure of the material distribution in the interface region. Two unit cells are perfectly connected when 𝐶𝐶𝐶𝐶 = 0, which implies that their connection regions are symmetric about the interface, while 𝐶𝐶𝐶𝐶 = 1 is completely disconnected.
In numerical implementation, CI can be calculated as
, (𝑖𝑖, 𝑗𝑗 = 1,2; 𝑖𝑖 ≠ 𝑗𝑗) (13) with
and
. Furthermore, the corresponding shape derivative with respect to the 𝑖𝑖-th microstructure can be written as:
where Γ 𝑖𝑖 s denotes the microstructural boundary in the strip region Y 𝑖𝑖 s .
For the proposed CI, one is required to select a suitable width of strip region.
Generally, for a smaller width of the strip region, CI is more sensitive to the boundary movement; once the structural boundary moves out of the strip region Y 𝑖𝑖 s , the sensitivity ( 14) would have no effect. A larger width would be numerically more stable and effective; however, it can potentially reduce the objective function to a greater extent. This will be discussed further in the numerical investigation.
We apply this with the level set topology optimization method for structural mechanics in this paper for illustration. However, we note that the proposed CI formulation is only a function of the material distribution and is independent of the physics of the problem. This means this constraint is applicable potentially to any topology optimization problems. Since we are introducing an additional constraint function into the problem formulation, this approach would be applicable with any topology optimization methods such as SIMP [30], ESO [31,32] and MMC/MMV methods [33][34][35][36].
In this section, we investigate the use of CI in connecting two or more microstructures using homogenization-based topology optimization. Two formulations are investigated: (i) adding penalized CI to the objective function (ii) adding CI as a constraint. The following section first details the two formulations, followed by numerical examples.
In this method, CI is added to the objective function of (10), i.e.,
where 𝑀𝑀 � is a weighting factor, 𝑆𝑆 ̅ denotes the average of the sensitivity of 𝑓𝑓 and 𝑘𝑘 is the iteration number. In this way, the penalization term of CI is applied to the objective function gradually in order to not over-restrict the microstructures at the initial steps.
The second formulation is to add CI as a constraint, i.e., 𝐶𝐶𝐶𝐶 ≤ 𝐶𝐶𝐶𝐶 ��� in (10) to enforce the microstructural connectivity. When determining the optimal boundary movements of the e-th microstructure, the CI-related constraint reads
where 𝐶𝐶𝐶𝐶 ���
The purpose of this section is to numerically investigate the connectivity index function within microstructure optimization and compare the two formulations of Section 5.1. We therefore, first construct particularly challenging two cell optimization problems with a single interface, where the optimized topologies of the two adjacent cells are completely disconnected. They are shown in Fig. 6. The first problem, Fig. 6(a) optimizes the left hand side cell for the maximum shear modulus and the right hand side cell for the maximum bulk modulus, i.e., 𝑓𝑓 1 (𝔼𝔼
Following the procedure outlined in Section 2, optimal microstructures can be obtained.
Due to the periodic boundary condition, each cell optimization is unaware of its adjacent cell and the resulting topologies are completely disconnected (the interface regions are highlighted with dotted lines). In practice, these solutions cannot transfer loads between two adjacent cells and are only fictitious designs. For all solutions, 𝐸𝐸 = 1 and 𝜈𝜈 = 0.3 with 𝑤𝑤 1 = 𝑤𝑤 2 = 40% were used.
Tables 1 and2 summarize the optimized solutions obtained by adding CI as a penalized objective with a range of weighting factors, 𝑀𝑀 � . The design domains were discretized by 50 × 50 bilinear elements and the widths of strip regions Y 𝑖𝑖 s are set as 1 element-width. It is observed that as the weighting factor increased, the resulting topological solutions connected at the interface. It is anticipated that the effects of CI would reduce the moduli of the optimal solutions. As seen in Table 1, the reduction of the optimal objective function relative to those of the optimized solutions without CI is less than 2%. However, for inappropriate weight factors, step changing similar as Fig.
1(b) still exists. This is due to the fact that the penalization approach comes from the multi-objective optimization, and selecting appropriate weighting factor for converting it to a single-objective optimization is a trivial task and would be very difficult when multiple interfaces are involved.
We now apply the CI constraint to multiple interfaces. Three microstructural optimized solutions with a 40% volume constraint are used in this study: (1) maximum shear modulus, (2) maximum bulk modulus and (3) maximum 𝐷𝐷 22 H . For these results, we use the initial design with five circular holes, Fig. 7(a). As a result, the resulting optimal topologies of Fig. 7(b)-(d) are slightly different from those using the one-hole initial solution (Fig. 6) reflecting the non-convex nature of microstructural optimization [7,29]; however, the objective functions of these solutions are different only by around 1% (0.110 (maximum shear modulus), 0.138 (maximum bulk modulus) and 0.400
) for the five-hole initial solution in comparison with 0.109, 0.136 and 0.400 for the one-hole initial solution). These provide additional optimized topologies to challenge the CI function. 5 demonstrates that the connectivities are improved substantially. For all cases, the maximum moduli were reduced by no more than 3%.
The efforts above focused on improving the connectivity between two adjacent microstructures. By applying the connectivity constraint, the connection between different cells will be greatly improved while the symmetry of unit cells can be broken.
Subsequently, the connectivity (i.e., symmetry of the interface region) between the same unit cells (e.g., the connection of the right-side unit cells in Fig. 8) may not as perfect as the connectivity of different unit cells.
If one wants to further improve the self-connectivity of the microstructures, the CI constraint can also be applied to each cell as shown in Fig. 9(a), i.e., It should be noted that, due to numerical errors and limitation of optimizer, the obtained designs can hardly achieve a perfect connection (i.e. 𝐶𝐶𝐶𝐶 = 0). However, it is expected that, only a 'minor' postprocessing is required to smooth the part in the interface region to finally improve the connectivity, and such treatment would have very small influence on the optimality of the microstructures.
We return to multiscale optimization in Section 3 and apply additional CI constraints to enforce connectivity between the microstructural regions. This is achieved by adding the CI constraints only at the microstructural scale in the decomposed optimization problem. For the 𝑒𝑒-th microstructure with a number of 𝑇𝑇 boundaries interfacing different microstructures, ( 20) is added to obtain the optimal boundary movement at 𝑘𝑘-th step:
It is noted that the CI constraint of the 𝑒𝑒-th microstructure is dependent on its adjacent cells, in other words, the connectivity constraint couples the adjacent microstructures together. In the current investigation, such a coupling is ignored for simplicity and CI is updated iteratively, i.e., 𝐶𝐶𝐶𝐶 ��� 𝑙𝑙
The first example considers two cantilevered beams shown in The optimized solutions are summarized in Tables 6 and7 for a range of the strip width, 𝑛𝑛 𝑤𝑤 . As can be observed, all microstructures are connected and in that sense they are manufacturable and realistic. However, the discrete member size changes in lower 𝑛𝑛 𝑤𝑤 would lead to poor load transfer between different type of cells and reduce the mechanical performance of optimal designs. The member size changes become more continuous as 𝑛𝑛 𝑤𝑤 increases. As would be expected, this has an effect of increasing the objective function however, the increases are observed to be small (less than 2%). We note that the constraints are not satisfied in some solutions shown in Tables 6 and7. In order to investigate this, we examine the optimization history of the cantilevered beam of 2:3, Fig. 12 which is representative of both of the cantilevered beams. It can be seen that the structural mean compliances converge smoothly for all cases and the strip width 𝑛𝑛 𝑤𝑤 does not have a significant influence. The 𝑛𝑛 𝑤𝑤 constraint oscillates particularly when 𝑛𝑛 𝑤𝑤 is small, e.g., 𝑛𝑛 𝑤𝑤 = 1. This is because when 𝑛𝑛 𝑤𝑤 is small, even a small perturbation of the boundary can have a significant influence on the CI value. Increasing the strip width, e.g., 𝑛𝑛 𝑤𝑤 = 5 , has an effect of relaxing the constraint and the oscillation reduces leading to a more stable convergence and meeting the specified constraint. It can be deduced from this that an adaptive CI constraint may lead to a stable convergence and this will be investigated further in the following section.
We now consider the L-beam shown in Section 3 to minimize the overall structural compliance with three regions are specified to have different material architectures. In this case, the topologies are optimized both at the macroscopic and the microscopic scales. Since there are three material microstructures, two CI constraints are added to the microscale optimization (to get a better connection as well as preserving the original optimality, the CI bounds are set as 0.08 for this example). Three cases are considered using the following parameters with the last case adopting the adaptive CI constraints:
1. 𝑛𝑛 𝑤𝑤 = 2, 𝐶𝐶𝐶𝐶 ��� 1 = 𝐶𝐶𝐶𝐶 ��� 2 = 0.08; 8. Fig. 13 with a small strip region, 𝑛𝑛 𝑤𝑤 = 2, shows that the microstructures at interface 1 is not well connected (𝐶𝐶𝐶𝐶 1 = 0.24 and 𝐶𝐶𝐶𝐶 2 = 0.10). As seen in the previous example in Section 6.1, for Case 2, increasing the strip width to 𝑛𝑛 𝑤𝑤 = 10 leads to a solution that satisfy both of the CI constraints, i.e. 𝐶𝐶𝐶𝐶 1 = 0.08 and 𝐶𝐶𝐶𝐶 2 = 0.07, with relative increase of objective function value by 2.97%.
However, an additional horizontal bar (marked by dash-dotted circle) is generated to satisfy the CI constraint. Furthermore, the iteration history of Case 2 illustrated in Fig. 15 reveals that the optimizer tried to satisfy the CI constraints every step and this may be overly restricting the search space for the microstructure.
The adaptive constraint in Case 3 is introduced to avoid such issues as the CI constraints are not strictly enforced from the beginning. This offers a greater level of design freedom for the microstructures during the early stages of optimization. Fig. 16 shows the microstructural solutions that are well-connected at the interfaces with 𝐶𝐶𝐶𝐶 1 = 0.08 and 𝐶𝐶𝐶𝐶 1 = 0.04. Moreover, the overall compliance increase is only 0.77%. Fig. 17 shows the optimization history of the L-beam Case 3 with the adaptive CI constraints. The CI constraints are inactive at the beginning and the CI values grow quickly. After about 30 iterations, both of the CI values start to decrease and stabilize around 120 iterations. It should be point out that, due to the connectivity constraint which requires the material distribute symmetrically in the interface region, the material distribution in the interface region (e.g., variable thickness of unit cell 1) is not optimal from a point of view on pure-stiffness. This will be considered in the future work.
Connectivity Index (CI) is formulated as a quantified measure of the connectivity between two adjacent microstructures. The CI function was introduced to the optimization problem, firstly into the objective function and secondly as an explicit constraint. The numerical investigations were conducted to explore the effects of the two different formulations as well as the interface region widths for a range of problems in which the microstructures are optimized via the homogenization approach. The introduction of the CI functions was demonstrated to improve the connectivity at the interfaces substantially. As an additional constraint can reduce the design space and restrict the search, the objective function values of the final solutions are expected to be compromised. The numerical investigations show, however, that the increase is minimal, mostly in the order of 1%. It was observed there were cases that a solution that satisfied the strict CI constraint could not be found or the optimizer found a trivial solution where the connectivity is achieved by making the interface region solid. In such cases, an adaptive strategy where the CI constraint was relaxed in the early stages of optimization and this gave the optimizer the freedom required to find good microstructure topologies. As optimization progresses and the number of iteration increases, the CI constraint is enforced more strictly yielding a satisfactory solution with the minimal increase in the objective function. It is noted that the CI function is completely geometry dependent and independent of the physics of the environment, therefore, the CI function approach is applicable to multiphysics topology optimization and this will be explored in our follow up study.