A Novel Topology Optimization Methodology Based on Energy Minimization Via α-β Swap Move

To facilitate the creation of new holes and to prevent simultaneously a checkerboard pattern in a topology optimization (TO) procedure, a novel TO methodology based on energy minimization via $\alpha-\beta $ swap is proposed. An energy function, measuring simultaneously the accuracy and piecewise smoothness of a new topology, is proposed and employed to determine the optimal way to update the current topology. Also, the proposed alpha-beta swap operations enable different materials to be exchanged with one another, and thus can create new holes naturally without including any additional mechanism. Moreover, the embedded smooth term in the energy function is capable of constraining the irregularity of the topology and thus hindering the checkboard pattern. Two case studies, an inductive power transfer system and a printed circuit board (PCB) are investigated to demonstrate the feasibility and advantage of the proposed methodology. The tools used for the two case studies include Matlab and Ansys. As demonstrated by the numerical results, the proposed methodology outperforms conventional methods in terms of both enhancing the performance parameter and constraining the checkerboard pattern.


I. INTRODUCTION
Topology Optimization (TO) is a conceptual design of a product and also the highest level in the design phase. It can create a novel topology, which is beyond one's imagination, for a product in its early conceptual design phase. Nowadays, TO has become a new paradigm to provide a quantitative design method for modern devices [1], [2].
The methodologies for TO can be categorized into two groups: the element-based method and the boundary-based method. Element-based methods, such as Solid Isotropic Microstructures with Penalization (SIMP) and ON/OFF method, split the design domain into discrete elements [3], [4]. The design variables then become the attribute parameters of materials within each element, enabling dramatic changes in shape and topology. One critical bottleneck problem of existing element-based methods is the checkerboard pattern in the final solution, which is regarded as numerical artifacts and unacceptable in engineering applications [5], [6]. Consequently, different methods have been proposed to eliminate the checkerboard patterns.
The associate editor coordinating the review of this manuscript and approving it for publication was Wai Keung Fung .
It is suggested that the checkerboards are typically prevented when 8 or 9-node quadrilaterals in combination with an element wise constant discretization of density is used [7]. A patch technique, controlling the formation of checkerboards in meshes of 4-node quadrilateral displacement elements coupled with constant material properties within each element, is proposed [8]. Also, in [7], a perimeter control technique is investigated. Constraining this perimeter clearly limits the number of holes that can appear in the domain, and the existence of solutions to the perimeter controlled TO is actually assured for both the discrete 0-1setting and the interpolated version using SIMP. However, the aforementioned mechanisms always involve additional constrains and excessive computational cost [9], [10].
Boundary-based method is another standard approach for TO. This paradigm inherently produces designs with clearly defined boundaries, and the solutions do not possess checkerboard patterns. However, one crucial issue in the boundary-based methods is the inefficient ability to insert or create new holes, where the advantage of TO over other optimization tools lies. In the existing TO procedures of this category, such as Level Set (LS) method, the merging and cancellation of holes can be operated in a natural way, while generating new holes is rather awkward and difficult. One remedy to resolve this shortcoming is to initiate the optimization from structures with many holes, as often used in LS approaches. However, this could lead to a strong dependence on the initial design [11]. To address this problem, different approaches have been proposed [12]- [14]. The topological derivative (TD) approach is a popular mechanism to enable hole creations in the direct LS method by combining a fixed point method with the natural extension of the classic shape gradient [15]. The main idea behind the TD is the evaluation of the sensitivity of the cost function to the creation of a hole inside the domain. However, many deficiencies exist in the TD methods. To start with, the number of iterations between the hole creations is unpredictable, and this can lead to a slow convergence or even to a suboptimal solution. Moreover, it is difficult to decide the instant either to start a new phase to create a hole or to continue the current boundary updating procedure [16]. In addition to TD methods, some other methodologies for hole insertions have been investigated. In [17], a combined shape and TO method is proposed, which uses a heuristic hole insertion rule that utilizes averaging over small element patches. However, this hybrid method is computationally inefficient, because only one hole can be inserted per Finite Element Analysis (FEA). A new method based on the boundary element and TD is proposed in [18]. This method selects the internal points with the lowest TD and creates holes by removing the selected points. However, the radius of the holes to be created is dependent on the spaced grid of internal points. A coarse grid of internal points leads to larger holes. Also, a degree of high complexities is always involved.
In a word, intractable issues remain in both the elementbased and the boundary-based methods for TOs. Table 1 tabulates the main cons and pros of the aforementioned two category methods. To solve these issues, a new TO methodology based on the Energy Minimization via α-β Swap Move (EMSM) theorem and related algorithm is firstly proposed. By converting the element-based sensitivity information of the finite element mesh into a weighted network, the new synthetic algorithm performs alpha-beta swap operations to enable the elements of different types of material to be exchanged to create new holes in the topology. Also, an inherent smooth term in the energy function is introduced and optimized to prevent the checkerboard patterns in the proposed methodology. As compared to the existing TO techniques, the effects of the aforementioned two approaches are twofold. The first is the increase of the algorithm diversity as coming from the new hole generation mechanism, and the second is the high computational efficiency as arising from the simple checkboard pattern prevention approach because there is no additional numerical operation for the proposed checkboard pattern prevention approach while complex and awkward numerical procedures are essential for existing checkboard pattern prevention approaches. In other words, a good balance between the global searching ability and the computational efficiency is obtained in the proposed TO methodology. Also, the numerical result confirms that this methodology is well-performed in terms of both generating new holes and constraining checkerboard patterns.
This paper is organized as follows: Section 2 introduces a new energy function to evaluate different directions to update the current topology, and explains the calculation of two main parts, the data term and the smooth term. Section 3 describes the implementation of α-β swap move, and establishes the relationship between minimizing energy function and finding the minimum cut in the network for α-β swap move. Section 4 shows case studies and the analysis of the corresponding results. Section 5 concludes this article.

II. THE OPTIMAL CHANGE DIRECTION IN TOS THROUGH ENERGY MINIMIZATION A. ENERGY MINIMIZATION IN TOS
The key to TO problem is to figure out the best way to change the current material distribution to a new one. Two indicators, accuracy and piecewise smoothness, are introduced in this paper to evaluate the change direction of the current topology. Accuracy implies whether the elements carrying the same level sensitivity value are assigned to the same material, and the piecewise smoothness is an evaluation of the checkerboard pattern. In this paper, the two indicators are measured simultaneously by the following proposed energy function where X 0 and X represent the current topology and the updated one, respectively, X 0 → X represents the way of changing from X 0 to X , N is the set of the interacting pairs of elements, D p measures how well label f p (material property assigned to element p) fits element p given the sensitivity data, and V p,q measures the extent to which a pair of adjacent elements is piecewise unsmooth.

B. SENSITIVITY BASED DATA TERM
In a TO, the sensitivity, SE, of an element is used to update the material attribute of the element. At each iteration, based on the classical discrete optimality criteria for a TO problem [20], the material region tends to evolve to maximize the sum of the sensitivities of all elements in this region. Since the objective of the data term, D p (f p ) in (1) is similar to that of the sensitivity, (1) is modified to where SE fp (p) is the calibrated sensitivity of the original one using where se r−k is the sensitivity changing from material r to k, SE r is the probability that the element is assigned to material r. Since the energy function is minimized, the element that should be turned from r to k has positive SE r and negative SE k .

C. SENSITIVITY BASED SMOOTH TERM
As is well known, the checkerboard pattern refers to the phenomena of alternating presences of solid and void elements ordered in a checkerboard like fashion. In (2), the smooth term V {p,q} (f p , f q ) measures the extent to which the material property of a pair of elements is piecewise unsmooth, in other words, the degree of checkerboard patterns. The topology with more severe checkerboard pattern results in larger smooth term. Therefore, it is reasonable to prevent the checkerboard pattern by minimizing the smooth term. Based on this analysis, the smooth term of (1) is modified as where K is a constant. The larger the K is, the greater the checkerboard pattern is constrained. Finally, finding the best way to update the current topology of TO problem is equivalent to minimizing the following function: For a two-phase materials design problem, the goal of TO is to find the best material distribution f that assigns each element x ∈ X a label f x ∈ L (finite set) (f is both piecewise smooth and consistent with the original sensitivity data) to divide the elements into two parts. A distribution f can be uniquely represented by a partition of elements P = {P l |l ∈ L} where P l = {p ∈ P|f p = l} is a subset of the elements assigned label l. Consequently, a distribution f is corresponding to a specific partition P [19]. For a pair of labels α, β, a move from a partition P (labeling f ) to a new partition P' (labeling f ) is called an α-β swap if Pl = P'l for any label l = α, β. This means that the only difference between P and P' is that some pixels that were labeled α in P are now labeled β in P', and some pixels that were labeled β in P' are now labeled α in P' as showed in Fig. 1. Clearly, the red material and the blue air are redistributed in the design domain by this operation, yielding a new holed (area 1 and 2) topology.

B. α-β SWAP MOVE AND NEW TOPOLOGY
Define G(V , E) as a weighted network with two distinguished vertices 's' and 't' called the terminals. An s-t 'cut' C is a set of edges that requires 's' and 't' to be in different subsets in the induced graph G(V , E-C), and it consists of edges whose removal divides the graph into two disjoint parts (Fig. 2). To connect the α-β swap move with an s-t cut, the meshed elements in the TO will firstly be transformed into a network. Each element p ∈ P is linked with elements in its neighborhood that touch one of its edge by e{p, q} called e-link. Additionally, each element is connected with terminals VOLUME 8, 2020 'α' and 'β', which denote two materials to be assigned to each element, by t α p and t β p called t-links, respectively. All kinds of mesh-grid could be transformed to a network regardless of the shape of the element. Fig. 3 illustrates the transformation of a quadrilateral and triangular mesh-grid [20]. An s-t cut in the network for α-β swap move is defined as a set of edges whose removal separates the graph into two disjoint parts which, respectively, contains 'α' and 'β'. Furthermore, the correspondence between cut C and the assignment of materials for each element is [4] where f C p indicates the new material attribute of element p determined by cut C.
In the network for α-β swap move, each element p ∈ P is connected with terminals 'α' and 'β' by t α p and t β p , and either of the two edges has to be included in cut C based on the definition of s-t cut. In other words, each element will be assigned one material once a cut C is determined. As shown in Fig. 4, each element is updated to either 'α' or 'β' according to (1), thus a new topology with updated material distribution is established.

C. α-β SWAP MOVE AND ENERGY MINIMIZATION
The edge weights of the network for α-β swap move are defined in Table 2. Specifically, the edge weights of two types of edge: t-link and e-link are, respectively, derived from data term and smooth term [20].

TABLE 2. Edge weights of t-links and e-links for EMGS in TO.
It has been proved in [20] that the cost of a cut |C| in the weighted network for the α-β swap move, the sum of the edge weights of it, equals to the energy function in (6) [20]. Also, it has been explained previously that a cut C in the network for α-β swap move is corresponding to a way to change the current topology. Thus, finding the best way of the material variation based on current topology, evaluated by energy function in (6), can be achieved by finding the min cut in the network for α-β swap move, as illustrated by Fig. 5. Consequently, the overall procedure of the proposed methodology is presented in Fig. 6.

D. INTUITIVE PROCEDURE FOR FINDING THE BEST WAY TO UPDATE THE CURRENT TOPOLOGY IN TO
To intuitively explain the mechanism of the proposed methodology, a step by step instruction is given using a 2 × 2 swap network as an example.
As explained previously, finding the best α-β swap move can be completed by transferring it to a min-cut problem to find a cut C with the minimum cost among all possible cuts. Moreover, the solution to a min-cut problem is dual to that of a max-flow problem of the same graph in combinational optimizations. In this regard, the following paragraphs will explain the procedures to find the solution of the max-flow problem of Fig. 7(a) [21]- [23]. The algorithm repeats the following three stages intuitively: 1) Construct a directed graph G(u, v) based on the original undirected swap graph of Fig. 7(a). Each edge carries both flow value f (u, v) and capacity c(u, v), which comes from the values for t-links and e-links. Initialize the graph and set the flow on each edge to be zero.
2) Find the augmenting path p in the corresponding residual graph G f . The residual graph consists edges with a residual capacity c f representing the flow on the edge to be changed. The edge can only be put into residual graph G f from the original graph G carrying a positive residual capacity c f . The rules for updating c f of one-way edge and two-way edges are given in (7) and (8), respectively.
An augmenting path p is a simple path connecting α and β. From the definition of the residual graph, the flow on an edge (u, v) in an augmenting path can be increased by up to c f (u, v) without violating the capacity constraints. The maximum amount by which the flow on each edge in an augmenting path can be increased is called the residual capacity of p, given by It should be pointed out that the flow on edge (u, v) can also be decreased by up to f (u, v) representing flow sending back through the edge (v, u). Fig. 7(b)-(d) show the three successive augmenting paths and the corresponding residual graphs. Taking the 2 nd iteration in Fig. 7(c) as an example, an augmentation path can be established consisting edges {(α, a), (a, b), (b, β)}, the residual capacity c p of the path equals to min{30,6,24} = 6, and the value of the max-flow is added up to 1+6 = 7, the residual capacity of the edges (α, a), (a, b), (b, β), and (b, a) can be updated to 24, 0, 18, and 12 respectively. This phase will be terminated when the residual graph G f contains no more augmenting path.
3) Find the minimum cut C dual to maximum flow path. In a residual graph, the edges carrying the zero residual capacity is considered as saturated ones. When the residual graph contains no more augmenting path, the vertices that can only be connected to α or β with unsaturated edges naturally become two separated parts. In Figs. 6(e)-(f), the vertices {a, d} and {b, c} can only be linked to α and β with unsaturated edges respectively, and thus they become two separated subparts.

IV. NUMERICAL RESULTS
In the following two case studies, the tools used for performance computations and figures productions include Matlab and Ansys software. Also, since the performance computation using high fidelity model (finite element model in the following case studies) takes the most computational sources and the CPU times for an inverse or optimization problem, the number of iterations that call for high fidelity model computations is used as a measure of the computational burden of a TO procedure, as commonly practiced in literature.

A. PRINTED CIRCUIT BOARD
A Printed Circuit Board (PCB) substrate is used to validate the feasibility and advantage of the proposed method. One function of a PCB substrate is to dissipate the maximum amount of thermal energy with limited amount of material while at the same time withstanding the induced thermal stresses [25]. As showed in Fig. 8, in a prototype PCB substrate discretized into 53 × 33 quadrilateral elements, four steady heat flux inputs, F1, F2, F3 and F4, are set to be 0.2KW/mm 2 . The temperature on the exterior boundary is maintained at 0 • C. The mechanical boundary condition is: a fully clamped on the four edges is considered. The objective is to maximize the average heat flux (HF, J/s) passing through the material elements: where V e = 1, element e is real material 0, element e is air (11) and where hf x , hf y and hf z are the heat flux of element e in x, y, and z direction.   9 illustrates the variation of material distribution of the PCB substrate in the process of TO. It is shown that the holes (represented by air) emerge in early stages of the optimization and then evolve iteratively. Fig. 10 shows the evolution of the objective, the average heat flux, and material percentage in the process of TO. Fig. 11 depicts the contour of the heat flux of initial design and optimized one, respectively. Fig. 12 demonstrates the comparison of flux histogram between the initial design and the optimized one, from which the efficiency of the thermal flux can be observed.
It should be noted that no significant difference can be ascertained between the optimized design obtained by the proposed methodology and the one in [21] using the Evolutionary Structural Optimization (ESO) method. However, the total iteration number by the proposed method is about 14, which is dramatically lower than that of ESO (about 50).

B. INDUCTIVE POWER TRANSFER
The magnetic core of an inductive power transfer system is solved by the proposed TO methodology to validate the effectiveness and efficiency. A wireless power transfer model is given in Fig. 13, where the transmitter and the receiver windings are placed in parallel in the vacant region. The current of 0.5A is applied to a 10 turns transmitter winding. The design objective is to maximize the mutual inductance between the two windings by distributing the magnetic material in the design domain. The design domain is divided into two regions, separated by the air gap [24]. The mutual inductance is maximized by maximizing the flux linkage of receiver winding. Thus, the objective function in this example is defined as the flux linkage of receiver winding; where and A w is the area of the receiver winding. Fig. 14 illustrates the process of TO using the proposed method. After the 1 st iteration, the material is added into the design domain which is wholly occupied by air initially, indicating that the holes are created (holes are represented by the iron in this example). Then, the topology evolves iteratively to the optimal one with no checkerboard-pattern. Fig. 15 depicts the flux line of the initial design and the optimized one, respectively. Fig. 16 shows the optimized topologies obtained by the hole sensitivity method [24], ON/OFF method, and the   conventional C-shape one. From Fig. 16, it is clear that the proposed method outperforms other conventional methods in VOLUME 8, 2020 terms of both enhancing the objective function and constraining checkerboard pattern.
It is deserving to mention that the coercivity and remanence of the magnetic core, if considered in the high fidelity model (finite element model), will affect the performance parameters and consequently the absolute solution time or efficiency in this case study. However, since all the TO procedures use the same high fidelity model in performance calculations, the relative performances comparison results as given in the manuscript will be fair.

V. CONCLUSION
In this paper, the EMSM algorithm and theorem is firstly employed to a TO problem to introduce a new TO procedure with the exclusive characteristics of new hole generating ability and checkerboard suppression mechanism. An energy function, measuring simultaneously the accuracy and piecewise smoothness of a new topology, is employed to determine the optimal way to update the current topology. Without adding random terms artificially throughout the optimization process, the holes emerge and update in a natural way by the α-β swap moves. Moreover, some modified smooth term is included to constrain the checkerboard pattern. Based on the numerical results reported, the proposed method is capable of both significantly enhancing the performance parameter and keeping a good balance between inserting new holes and constraining checkerboard pattern. Since the proposed graph model is focused on two-dimensional topology optimizations, the future direction of this research is to extend it to three-dimensional TO problems.