Evolutionary topology optimization using the extended finite element method and isolines

This study presents a new algorithm for structural topological optimization of two-dimensional continuum structures by combining the extended finite element method (X-FEM) with an evolutionary optimization algorithm. Taking advantage of an isoline design approach for boundary representation in a fixed grid domain, X-FEM can be implemented to improve the accuracy of finite element solutions on the boundary during the optimization process. Although this approach does not use any remeshing or moving mesh algorithms, final topologies have smooth and clearly defined boundaries which need no further interpretation. Numerical comparisons of the converged solutions with standard bi-directional evolutionary structural optimization solutions show the efficiency of the proposed method, and comparison with the converged solutions using MSC NASTRAN confirms the high accuracy of this method.


Introduction
In recent years, structural optimization has become a rapidly growing field of research, with application in many areas such as mechanical, civil and automotive engineering. Topology optimization is one of the most challenging aspects of structural optimization, in which one needs to find the best topology as well as shape of a design domain. The approaches that have been proposed for the topology optimization of continuous structures fall into two categories: first, mathematically based methods such as homogenization (Bendsøe and Kikuchi 1988), solid isotropic material with penalization (SIMP) (Bendsøe 1989;Zhou and Rozvany 1991) and level set method (Wang, Wang, and Guo 2003;Allaire, Jouve, and Toader 2004); and second, heuristic methods which are more intuitive and less mathematical, such as evolutionary structural optimization (ESO/BESO) methods (Xie and Steven 1993;Querin, Steven, and Xie 1998;Yang et al. 1999).
ESO is based on the assumption that the optimal layout of the design domain can be obtained by gradually removing inefficient material from the design domain (Huang and Xie 2009). In the original ESO method, the elements of the design space were ranked in terms of their sensitivity, and those with lower sensitivity were removed from the design domain until a converged solution was obtained. Bi-directional evolutionary structural optimization (BESO) is an extension of ESO in which the elements are allowed to be added and removed simultaneously. These heuristic methods are easy to program and provide a clear topology (no grey regions of intermediate *Corresponding author. Email: eaxma5@nottingham.ac.uk densities as in SIMP) in the resulting optimal designs. Conventional ESO/BESO algorithms have been successful since they can be easily combined with the finite element model of a structure. However, they suffer from a weak capability of boundary representation as they are defined by the finite element mesh, which is non-optimal with respect to the final converged solution. This limitation causes difficulties in combining these methods with computer-aided design (CAD), and the obtained solutions require postprocessing to manufacture a design with a smooth and manufacturable surface.
The fixed grid finite element method (FG-FEM) allows the boundaries of the design to cross over finite elements. This capability has been used in boundary-based optimization methods such as the level set method, and element-based optimization methods, such as fixed grid evolutionary structural optimization (FG-ESO). FG-ESO or isoline/isosurface approach (Victoria, Martí, and Querin 2009;Victoria, Querin, and Martí 2010) is an alternative to ESO in which the inefficient material is allowed to be removed or added within the elements of the design domain during an evolutionary process. The boundaries are defined by the intersection of the isoline plane with the criteria distribution of the design domain. Since in this approach the boundary of the design is no longer consistent with the fixed finite elements as in ESO, a classical finite element analysis may result in a poor finite element approximation on the boundary. Conventionally in the fixed grid finite element approach, the element stiffness is assumed to be proportional to the area fraction of the solid material within the element (also called the density scheme). Although this approach is widely accepted and implemented in many works (Allaire, Jouve, and Toader 2004;Victoria, Martí, and Querin 2009), studies have shown that it cannot provide accurate results for the boundary elements (Dunning, Kim, and Mullineux 2008;Wei, Wang, and Xing 2010).The extended finite element method (X-FEM) is another fixed grid approach which can be used to model void-solid interfaces. X-FEM extends the classical finite element approach by adding special shape functions which can represent a discontinuity inside finite elements. In this approach, the geometry of the discontinuity is often described by a level set method. A combination of the level set description of the geometry and the fixed mesh framework of X-FEM has been used in recent level set-based topology optimization work (Miegroet and Duysinx 2007;Wei, Wang, and Xing 2010).
This study presents a simple and effective evolutionary optimization approach in a fixed grid domain. The novelty of this work is to apply X-FEM to the evolutionary optimization algorithm. The proposed method does not require a level set framework for geometry description in the X-FEM and the boundaries of the design can be simply represented by isolines of a desired structural performance. The algorithm is implemented in the topology optimization of two test cases and the final solutions are compared to standard BESO solutions. To evaluate the accuracy of the proposed method, the solutions are imported to NASTRAN and reanalysed using the classical FEM.

Fixed grid approach
The fixed grid method is a technique to model the boundaries with a non-conforming mesh. It has been used for solving problems with moving boundaries such as those requiring structural optimization problems. Unlike the remeshing methods in which the design domain, or a narrow area around the boundary, is remeshed in every iteration, the fixed grid FEM does not require a time-consuming remeshing process and can be easily implemented. In this method the real structure is superimposed on the fixed finite elements of the design space and makes three types of element: the solid elements that lie inside the structure (S elements), the void elements that lie outside the structure (V elements) and the elements that lie on the boundary of the structure (B elements), as shown in Figure 1.

Density scheme
In a conventional fixed grid approach, the stiffness matrix of boundary elements is approximated by a density scheme in which the stiffness of the element is proportional to the area ratio of the solid part of the element: (1) where A (e) S , A (e) v and A (e) t represent the solid area, void area and the total area of the element, respectively, and K B and K S are the stiffness matrices of the boundary element and solid element, respectively.
In the density scheme, the material is considered to be uniformly distributed through the whole element and the variations in material distribution in an element are not taken into account in calculating the element stiffness matrix. For example, Figure 2 shows three different shapes for a boundary element where the area fraction of solid material within the element is 0.50. Using the density method (Equation 2) the same stiffness is calculated for all three elements. This method may cause errors near the boundary of the design during the optimization process.

Extended finite element method
X-FEM (Belytschko and Black 1999; Moës, Dolbow, and Belytschko 1999) is an alternative fixed grid approach in which the classical finite element approximation is enriched by special functions through the concept of partition of unity (Melenk and Babuška 1996). X-FEM was originally developed to represent crack growth in a fixed grid domain. Using traditional FEM for the simulation of crack propagation is very challenging because of the continuous changes in the topology of the domain. The application of X-FEM for this case has been very successful because the finite element mesh can be generated independent of the geometry of the crack, and remeshing is not required during crack propagation. For crack modelling, the idea of X-FEM is to use additional degrees of freedom in the usual finite element spaces by adding a discontinuous function (Heaviside step function) and the asymptotic crack tip displacement fields to the conventional finite element displacement field. Therefore, the standard finite element approximation is extended by adding enrichment functions to the continuous displacement field and defining additional degrees of freedom for the nodes which support the discontinuity. The X-FEM equation for the regions away from crack tip (ignoring tip enrichment) has the form: In the above equation, the first function on the right-hand side shows the conventional finite element approximation of the displacement field in an element where N i are the classical shape functions associated to the nodal degrees of freedom, u i . N j (x) H(x) supported by enriched degrees of freedom, a j , are discontinuous shape functions constructed by multiplying a classical N j (x) shape function with a Heaviside function H(x) presenting a switch value where the discontinuity lies. X-FEM has also been implemented for other kinds of discontinuities such as fluid-structure interaction (Gerstenberger and Wall 2008) and modelling holes and inclusions (Sukumar et al. 2001). In the present case the X-FEM scheme for modelling holes and inclusions can be implemented for modelling material-void interfaces during the optimization process. In this approach, the displacement field is approximated by the following equation where the Heaviside function H(x) has the following properties which implies that the value of Heaviside function is 1 for the nodes outside the void and 0 for the nodes in the interior of the void. Since there is no enrichment in the displacement approximation equation of X-FEM in modelling holes and inclusions, there will be no augmented degrees of freedom during optimization. The proposed X-FEM integration scheme will be discussed in the next sections.

Structural optimization problem
The topology optimization problem where the objective is to minimize the strain energy can be written as Subject to : where c is the total strain energy, and U and K are the global displacement and global stiffness matrices, respectively. N denotes the number of finite elements in the design domain. v (e) S is the volume of the solid part of the element, V 0 is the design domain volume and V * is the prescribed volume fraction. While in ESO/BESO methods, the presence or absence of each element in the design domain is considered as a design variable, in the proposed method the material distribution inside each element is considered as a design variable. Strain energy density (SED) is a reliable criterion to indicate inefficient use of material in a design space (Huang and Xie 2009). Elemental SED can be obtained from with u e the element displacement vector and k e the element stiffness matrix, which is calculated using an X-FEM scheme. The topology optimization operates by gradually removing material from low SED regions and adding it to high SED regions during an evolutionary procedure. The effective removal/redistribution of material within the design domain can be achieved through an isoline topology optimization approach by assigning a weak material property to low SED regions and solid material property to high SED regions (soft-kill scheme).

Isoline topology optimization
The basic idea of isoline design is to represent the shape and topology of the structure using the contours of the desired structural behaviour. This idea has been suggested in several studies (Maute and Ramm 1995;Lee, Park, and Shin 2007). The isoline boundary representation differs from the level set description of the design boundary. In the level set method the boundary is described by a zero level set function (usually a signed distance function) ( Figure 3) and the evolution of boundary is represented by solving a Hamilton-Jacobi equation. In the isoline method, the design boundary is represented by a minimum level of a criterion (such as von Mises stress or SED) which is iteratively updated during the design process ( Figure 4). The isoline optimization algorithm that is used in this article originates from the isoline topology design (ITD) algorithm proposed by Victoria, Martí, and Querin (2009). Unlike ESO/BESO methods in which the optimization operates at an elemental level, in the ITD approach the optimization is performed at a global level, based on structural performance. In the present study, the distribution of the desired structural behaviour is obtained in the analysis phase through the use of X-FEM. The ITD approach can be summarized in the following steps: (1) An extended finite element analysis is performed to find the distribution of SED within the design domain. Steps 1-3 are repeated by gradually increasing the MSL until a desired optimum is obtained.

The evolutionary procedure in isoline topology optimization
The nodal level of SED in any part of the design domain can be obtained by performing an extended finite element analysis and comparing the SED of that node with the maximum SED of the entire domain. Using the nodal SED level, the elements are categorized into three groups: void elements, solid elements and boundary elements: boundary element if neither of the above conditions exists (9) where RF i is the current redistribution factor (RF), and j denotes the element's node number. The boundary of design can be obtained by the intersection of the SED distribution and the minimum SED level, which is calculated by With the current redistribution factor, the iterative process of the extended finite element analysis and material removal/redistribution takes place until the change in volume fraction is less than a minimum value V , which means that a steady state is almost reached. An evolutionary rate, ER, is added to the redistribution factor, such that With the new redistribution factor, the extended finite element analysis and material removal/redistribution are repeated until a new steady state is reached. The evolutionary process continues until a desired optimum, such as a final volume fraction or a maximum strain energy level, is reached. In this study, the volume fraction and the objective function (total strain energy) are used to check the convergence. The evolutionary process continues until the volume fraction condition is satisfied. From this time forwards, the optimization process runs with a constant redistribution factor (zero evolution rate) until the changes in the objective function in the last five iterations are within a 0.1% tolerance. The number of iterations in the evolutionary process in the proposed method is affected by the value of evolution rate as well as the maximum strain energy density (SED max ). Selecting a high evolution rate can reduce the computational time. However, very high values of evolution rate may result in local optima or non-convergent solutions. A typical value for the evolution rate can be obtained from where SED ave is the average SED for the fully solid design domain.

X-FEM integration scheme
Equation (4) defines a zero displacement field for the void part of the element, which means that only the solid part of the element contributes to the element stiffness matrix. Thus, it is possible to use the same displacement function as FEM (first term in Equation 3) and simply remove the integral in the void subdomain of the element: where S is the solid subdomain, B is the displacement differentiation matrix, D S is the elasticity matrix for the solid material, and t is the thickness of the element. When an element is cut by the boundary, the remaining solid subdomain is no longer the reference rectangular element. In order to numerically calculate the integral given by Equation (13), the solid part of the boundary element is partitioned into several subtriangles ( Figure 5) and the Gauss quadrature method is used: where n is the number of subtriangles inside the element and T denotes the triangle domain. Figure 5. The solid subdomain of the boundary elements is partitioned into several subtriangles.

Triangulation of boundary elements
The topology and shape of a boundary element can be found using the values of nodal relative SED, which is defined by: where j denotes the node number. Nodes with negative relative SED belong to the void part of the domain and nodes with positive relative SED are located in the solid domain. Therefore, an element which has at least one node with negative SED rel and one with positive SED rel is a boundary element. The intersection point of the boundary and element edge between two neighbouring nodes i and j can be found using bilinear interpolation of nodal SED rel and shape functions: where x im is the distance between node i and the intersection point m, and l ij is the element length between the nodes i and j. Depending on its topology, a boundary element may have two, three or four intersection points. The subtriangles can be defined by defining an extra point inside the solid subdomain of the element (typically in the centre of solid area) and connecting it to the solid nodes as well as the intersection points. Figure 6 shows a boundary element with two intersection points and typical values for SED rel .

Gauss quadrature method
In order to apply the Gauss quadrature method to triangles, the two-dimensional (2D) integrals in terms of the physical coordinates are transferred to the triangle's natural coordinates and represented as a series of weighted functions: Figure 6. A boundary element with typical values for nodal relative strain energy density.  where m is the number of gauss points, ξ are the coordinates of the gauss points, and A T is the area of the triangle. Substituting Equation (17) into Equation (14), the element stiffness matrix can be obtained by with n the number of subtriangles in the solid domain of the element. In this study, the second order gauss rule with three midline gauss points was implemented (Figure 7). To validate this X-FEM scheme, the stiffness matrix of a fully solid rectangular element with Young's modulus E = 1 and Poisson's ratio υ = 0.3 was calculated by two different approaches: first, using the classical finite element approximation; and secondly, using the X-FEM scheme described above in which the element is divided into subtriangles and integration is performed using gauss quadrature for triangles, as shown in Figure 8. Both methods resulted in exactly the same stiffness matrix for the element, thus validating the X-FEM scheme. Figure 9 illustrates the topology optimization procedure used, which in general consists of the following steps:

Combining X-FEM and the optimization algorithm
(1) Initialization: in this step, the dimensions of the design domain, fixed grid mesh and initial material distribution within the design domain are defined; boundary and loading conditions is checked by comparing the convergence criteria with the defined convergence threshold. If the convergence condition is satisfied, the algorithm jumps to step 9, otherwise it progresses to step 7. (7) X-FEM structural analysis: an X-FEM analysis is performed on the fixed grid design domain.
Using nodal SED numbers, the elements are categorized into three groups: solid, void and boundary elements. Solid and void elements are treated using classical finite element approximation. The stiffness matrices of the boundary elements are calculated by partitioning the solid subdomain into several subtriangles and applying the gauss quadrature integration scheme described in Section 3.2. The global stiffness matrix is calculated by assembling the element stiffness matrices. (8) Go to step 3. (9) Stop the optimization process.

Numerical examples
The proposed method of combining the X-FEM and evolutionary optimization algorithm was implemented in a MATLAB code to present the topology optimization of 2D rectangular domains as a first validation stage prior to full 3D implementation. Two test cases are used in this study ( Figure 10). A consistent dimensionless set of parameters is used for both test cases. Test case 1 was a short cantilever beam with length 60, height 30 and thickness 1, where a unit concentrated load was applied in the middle of the free end (symmetric case). Test case 2 was a cantilever beam with the same dimensions as test case 1 but with the load applied at the bottom of the free Figure 11. Evolution history of objective function, strain energy (SE) and volume fraction (VF) for (a) test case 1 and (b) test case 2.
end (non-symmetric case). The material properties of the solid material were Young's modulus E = 1 and Poisson's ratio υ = 0.3. The target volume was 50% of the initial design domain.
To avoid singularity issues with the concentrated loading, the strain energies inside the loading regions shown in Figure 8 were not used in the calculation of the total strain energy and the tip displacements were measured from outside the loading region along the line of loading.

Preliminary examination of convergence
The initial design domain was discretized using a 60 × 30 mesh. The optimization started with a fully solid design domain. The evolution histories of the objective function and volume fraction for test cases 1 and 2 are shown in Figure 11. It can be seen that the strain energy increases as material is gradually removed from the design domain, then reaches a constant value at convergence. The development of the topology in the iterative optimization processes for the two test cases is illustrated in Figures 12 and 13. It can be seen that initially a number of holes appears as the volume fraction decreases. After a certain number of iterations some of the holes merge to make larger holes, thus reducing the final complexity of the topology. Figure 12 shows that the topology remains symmetric throughout the optimization. It can be seen in the final topologies that despite using a coarse mesh for this optimization problem, the final designs have clearly defined, smooth boundaries and need no further interpretation (unlike standard SIMP and ESO/BESO methods).  In the next two sections, the accuracy of the results is studied and the method is compared with standard BESO.

Evaluating X-FEM solutions
In order to accurately evaluate the performance of the final solutions and the accuracy of the proposed method, the obtained solutions were discretized by a converged, fine structured finite element mesh and solved using the commercial finite element solver NASTRAN from MSC Software (Santa Ana, California, USA) ( Figure 14). Table 1 compares the X-FEM solutions and the converged NASTRAN solutions in terms of their strain energies and tip displacements. It can be seen that the X-FEM solutions are very close to the regenerated NASTRAN solutions with percentage error less than 0.7. The small difference in the X-FEM and NASTRAN results may be attributed to the different mesh size used in the two approaches and would be expected to decrease by reducing the mesh size of the design space. However, it could be argued that the accuracy obtained using the coarse mesh is sufficient for the topology optimization and the added accuracy of a finer mesh will unnecessarily increase computational time. This is an increasingly important consideration when the method is used for the optimization of real-life 3D structures (as will be shown in future works).

Comparison with BESO solutions
The solutions obtained from the proposed method are compared with BESO solutions for a range of mesh sizes. The BESO solutions were obtained using a soft-kill BESO MATLAB code (Huang and Xie 2009). In order to overcome the chequerboard problem (Jog and Harber 1996) in the BESO solutions and retain the complexity of the converged solutions, a small filter radius of 1.2 times the element length was used. The selected evolution rate for BESO was 0.004, which was chosen to give approximately the same number of iterations to converge as the X-FEM-isoline optimization approach (180-200 iterations).

Comparison of strain energies
Comparing the topologies obtained using the two approaches for test case 1 (the symmetric problem), one may notice that the converged solutions for the same mesh size have similar topologies (Table 2); however, the BESO solutions tend to have higher strain energies than the X-FEM solutions. The probable reason for this is the poorer edge representation in the BESO method, which has reduced the performance of the converged solutions. At high mesh density the two methods had very similar performance, as would be expected; however, this would obviously be more computationally expensive. To increase the performance of these BESO solutions, additional postprocessing is required to smooth the boundaries. In test case 2, which is a non-symmetric problem, the two approaches have generated different topologies (Table 3). The strain energies of the X-FEM solutions are again lower than the BESO ones, indicating better performance for the X-FEM solutions. It can be seen in Tables 2 and 3 that both methods result in final topologies that are mesh dependent, which is generally the case for element-based topology optimization methods. Pseudomesh-independent topologies can be obtained by the use of coarsening actions such as filtering the sensitivities or increasing the filter radius (Huang and Xie 2007). However, these methods still have their limits in terms of mesh independency and will tend to result in coarser solutions that can have lower performance than more refined solutions. Generally speaking, increasing the mesh density increases the complexity of the converged solutions and reduces the strain energies. Therefore, increasing the mesh density can improve the performance of the final optimized result. This issue has also been studied in earlier investigations for BESO (Aremu et al. 2013). Figure 15 illustrates the strain energy of the converged solution as a function of mesh density for test case 2. It shows that the X-FEM optimization approach is more robust than the BESO method and that the strain energies of the X-FEM solutions converged earlier than the BESO ones. However, as the mesh gets finer the strain energies of BESO solutions become closer to the X-FEM solutions. It can be seen that increasing the mesh density results in the strain energies produced using the X-FEM-isoline approach and NASTRAN finite element analysis become close, relative to the fluctuations in the data beyond 1/h = 80. Figure 17. Extended finite element method (X-FEM) solutions before and after smoothing. Table 4 compares the computational time for the two optimization methods for test case 1 after 100 iterations. It can be seen that at low mesh densities BESO is much faster than X-FEM; however, the computational time ratio decreases by increasing the mesh density. It is shown in the article that BESO solutions require higher mesh densities and more postprocessing to obtain a smooth topology; therefore, the total time for design could be decreased using the X-FEM-isoline method.

Comparison of surface roughness
As a postprocessing stage, a Laplacian smoothing algorithm can be used to create smoother boundaries for the optimized topologies, if this is required for manufacture, for instance. Laplacian smoothing is an iterative smoothing technique, commonly used in image processing and to improve the quality of finite element meshes (Cannan et al. 1993;Freitag 1997). In imageprocessing applications, Laplacian smoothing operates by replacing the grey value of a pixel with an average of the grey values of neighbouring pixels. In finite element meshing applications, the location of a vertex is modified using the average of the locations of neighbouring vertices (Vollmer et al. 2001). Figure 16 shows the boundaries before and after smoothing for the BESO solutions of test case 2 for a range of mesh densities. The average surface roughness of the optimized topologies before smoothing and the number of iterations in the Laplacian smoothing are also included in the figure. The root mean square roughness (R q ) was determined by comparison of the topology boundaries before and after the Laplacian smoothing. It can be seen in Figure 16 that the surface roughness of the topologies from the BESO optimization increases as the mesh density decreases. Figure 17 shows the boundaries of X-FEM solutions for the same test case before and after smoothing. Compared to the BESO solutions, the surface roughness of the X-FEM solutions is much lower and they need far fewer iteration steps in the Laplacian smoothing. It is also apparent that unlike the BESO solutions, the surface roughness of the X-FEM solutions have little dependency on the mesh density.

Conclusions
In this study, X-FEM and the isoline design method are combined with an evolutionary optimization algorithm to determine the optimal topology of 2D continuum structures. The results suggest that using X-FEM has significant advantages, as it not only avoids time-consuming remeshing techniques, but also generates structures that have smooth boundaries requiring little or no further interpretation or postprocessing. Using simple test geometries, it has been shown that X-FEMbased topology optimization has the potential for greater accuracy and more robust solutions with less dependence on mesh size than BESO, although this needs to be established for more realistic, 3D geometries. It is anticipated, however, that this method will be relatively simple to extend to three dimensions and a possible route to achieving this may involve using eight-node brick elements where the solid part of the boundary element could be represented by subtetrahedrons.