# Generalized Shape optimization based on the Level Set method

Laurent Van Miegroet<fnref n="1" />, Nicolas Moës<fnref n="2" />, Claude Fleury<fnref n="1" />, Pierre Duysinx<fnref n="1" />

## 1. Abstract

This paper describes a first step work devoted to applying XFEM and Level Sets methods in optimization of structures. This first step work is based on integrating an existing XFEM code within a general open optimization tool, SAMCEF BOSS QUATTRO. Unlike most of the existing works, this approach is more shape optimization oriented. A library of pre-formatted basic geometric entities (such as ellipses, squares, triangles, etc.) described by Level Sets functions are used. These basic Level Set features can be combined to represent many kinds of interfaces and holes. The construction parameters of the basic Level Sets are considered as the design variables. In order to evaluate the sensitivities, a finite difference scheme over the design variables is used in this first work. Different mechanical responses (energy, weight, displacement, ...) can be considered as objective functions or constraints in the problem formulation. Several academic 2D test cases of shape and topology optimization are presented within the XFEM and Level Set approach. In addition, a work by Missoum *et al.* [11], in which the shape and topology optimization of the structure is carried out by an optimal selection of holes characteristics with a genetic algorithm is presented.

**2. Keywords :** Topology and shape optimization, Level Set method, Extended finite element method.

## 3. Introduction

In the classical domain of structural optimization, two main techniques has been mainly studied and are now well known : the topology and the shape optimization. These two techniques have reached a certain degree of robustness and sophistication but still present some major drawbacks. Hence, in the method of shape sensitivity, the computational cost is very high due to the remeshing process. Also, the optimization procedure tends sometimes to fall into local minima (solution close to the initial design, no deep change allowed). However, this method presents the interesting aspect of using a classical finite element model and consequently can handle various objective functions and material laws. In opposition, the topology optimization method uses a fixed mesh thus avoiding the expensive mesh perturbation procedure. The method often presents good performance as the creation of new geometric entities is possible by varying the densities of the finite elements. However, it is not possible to use some specific material laws and, moreover, the aspect of the solution is always presented as a grayscale bitmap which needs to be interpreted in order to create a feasible finite element model for further analysis. Also, the only way to get a rather smooth solution is to increase the number of element leading to a growing computation time.

In order to avoid the main problems encountered with the classical optimization techniques, we propose in this paper a new method which may present all the advantages of this two ones and avoid the main difficulties associated with those. This method is based on two axes : the analysis method and the representation of the geometry. The structural analysis is realised by using the extended finite element method (XFEM) [1]. This enables the design to include features such as interfaces between material and void that are not coincident with the mesh. Hence, the mesh perturbation present in the shape optimization is suppressed as we keep the same mesh during all the optimization steps. The description of the geometry is represented by the zero iso-contour of an implicit function called the Level Set function [7]. This method is very convenient for the topology optimization as it allows deep and complex changes of the global shape of the model.

This work differs from other papers on the subject as we don't use the equation of motion of the Level Set method like [3], [4], [5] nor the nodal Level Set value as design variables like Belytschko *et al.* in [2]. In our work, a library of Level Set functions representing basic geometric entities (such as ellipses,

<footnote n="1">¹: Department of Mechanics and Aerospace, Université de Liège (ULg), Belgium. L.Vanmiegroet@ulg.ac.be</footnote>

<footnote n="2">²: Department of structures and simulations, École Centrale de Nantes (ECN), France</footnote>
squares, triangles, ...) can be combined together whereas the sensibilities are computed by finite difference over the parameters of these functions. This limitation comes mainly from the nature of this first work which objective was to rapidly integrate the existing XFEM code from Nicolas Moës with an general optimization tool BOSS QUATTRO from SAMCEF in order to evaluate the potential of the method.

This new optimization method based on XFEM and Level Sets exhibits really promising characteristics as it allows deep topological changes and a very flexible modelling of the geometry. However numerical applications pointed out some specific difficulties and problems to be handled. The origin of the problems is described in the paper and some remedies are presented to circumvent the problems.

The outline of this paper is as follows. In the next section, we give a brief introduction to the XFEM and the special case of modelling holes. In section 5, we introduce the basis of the Level Set method. Section 6 details two continuum variables test cases whereas section 7 presents a genetic algorithm approach. Finally, the conclusions and the perspectives at section 8.

## 4. The extended finite element method

The extended finite element method [1] is a recent method whose been firstly developed for the simulation and the analysis of structures presenting crack growth. The principal strength of this method is its capacity of including discontinuities inside the finite elements. Hence, this method enables the model to include geometric boundaries, material or phase changes that are not coincident with the mesh.

### 4.1. The basis of the method

In order to allow any types of discontinuities inside the elements and therefore to be able of representing discontinuities in the physics field of the problem's solution, it is necessary to add special properties to the shape functions. For example, in the case of cracked structures, the physical discontinuous field is the displacement field, hence, if one want to be able of modelling this discontinuity, we have to add discontinuous shape functions. The classical finite element approximation used is then extended to embed discontinuous shape function as in the following equation:

$$ \mathbf{u}(\mathbf{x}) = \sum_{i} u_i N_i(\mathbf{x}) + \sum_{j} a_j N_j(\mathbf{x}) H(\mathbf{x}) \quad (1) $$

where $N_i$ are the classical shape functions associated to degree of freedom $u_i$. The $N_j(\mathbf{x})H(\mathbf{x})$ are the discontinuous shape functions constructed by multiplying a classical $N_j(\mathbf{x})$ shape function with a Heaviside function (presenting a switch value where the discontinuity lies). These extended shape functions are supported only by the enriched (extended) degree of freedom $a_j$. Note that, generally, only the elements near the discontinuity support extended shape functions whereas the other elements remain unchanged. The modification of the displacement field approximation doesn't introduce a new form of the discretised finite element equilibrium equation but leads to an enlarged problem to solve:

$$ \mathbf{K} \cdot \mathbf{q} = \mathbf{g} \Leftrightarrow \begin{bmatrix} K_{uu} & K_{ua} \\ K_{au} & K_{aa} \end{bmatrix} \begin{bmatrix} u \\ a \end{bmatrix} = \begin{bmatrix} f_u^{ext} \\ f_a^{ext} \end{bmatrix} \quad (2) $$

With

$$ K_{uu} = \int_V B^T C B \quad (3) $$

$$ K_{ua} = \int_V \tilde{B}^T C B = K_{au} = \int_V B^T C \tilde{B} \quad (4) $$

$$ K_{aa} = \int_V \tilde{B}^T C \tilde{B} \quad (5) $$

As the elements can now present discontinuous shape functions, the numerical integration scheme has to be modified in order to take care about the discontinuity. In practice, the elements embedding a singularity are divided into subelements aligned with this discontinuity over which an under-integration is processed.
The modelling of material-void interfaces with XFEM (see [12]) is slightly different from the cracked structure's case. Hence, the displacement field is approximated by:

$$ \mathbf{u}(\mathbf{x}) = \sum_{i} u_i N_i(\mathbf{x}) V(\mathbf{x}) \qquad (6) $$

where $V(\mathbf{x})$ takes value 1 if the node lies inside the material and 0 unless. The elements lying outside the material are removed from the system of equations, whereas the partially filled elements are integrated using the XFEM under-integration procedure. Modelling holes with the XFEM is a very appealing method for the shape optimization but also for the topology optimization as no remeshing is needed and no approximation is done on the nature of the voids in opposition to the SIMP method.

## 5. The Level Set method

The explicit representation of the structure's shape that is used in the classical finite element method forbids deep boundary or topological changes such as creation of holes. This limitation is the main reason of the low performance generally associated to the shape optimization. In opposition, the Level Set method developed by Osher and Sethian [7] which consist of representing the boundary of the structure with an implicit method allow this kind of deep changes.

The Level Set method is a numerical technique first developed for tracking moving interfaces. It is based upon the idea of representing implicitly the interfaces as a Level Set curve of a higher dimension function $\psi(\mathbf{x}, t)$. The boundaries of the structure is then conventionally represented by the zero level ($\psi(\mathbf{x}, t)=0$) of this function $\psi$, whereas the filled region is then attached to the negative or the positive part of the $\psi$ function. In practice, this function is approximated on a fixed mesh by a discrete function which is usually the signed distance function:

$$ \psi(\mathbf{x}, t) = \pm \min_{\mathbf{x}_{\Gamma} \in \Gamma(t)} \|\mathbf{x} - \mathbf{x}_{\Gamma}\| \qquad (7) $$

The sign is positive (negative) if $\mathbf{x}$ is inside (outside) the boundary defined by $\Gamma(t)$. The evolution of the interfaces is then embedded in the evolution equation for $\psi$, which is given in [7] by:

$$ \frac{\partial \psi(\mathbf{x}, t)}{\partial t} + F \|\nabla \psi\| = 0 \qquad (8) $$

$$ \psi(\mathbf{x}, t) = 0 \quad \text{given} \qquad (9) $$

where $F$ is the speed function defined on the interface $\Gamma(t)$ in the outward normal direction to the interface. Applied to the XFEM framework, the Level Set is defined on the structural mesh and at each finite element node is associated a geometrical degree of freedom representing its Level Set function value. The Level Set is then interpolated on the whole design domain with the classical shape function:

$$ \psi(\mathbf{x}, t) = \psi_i N_i(\mathbf{x}) \qquad (10) $$

As example, the following figure (Fig. 1) illustrates the Level Set representation of a square and the XFEM model associated. The combination of different Level Sets is also one of the appealing characteristic of this method. Thus, for example, the union of two circles defined by their Level Set is easily given by the minimum of the two Level Sets.
<figure><img src="image_1.png" /><figcaption>Figure 1: Geometric representation and XFEM model</figcaption></figure>

# 6. Description of the method

Essentially two methods based on the Level Set applied to structural optimization has been proposed. The main difference between them is the method used for the sensitivity analysis :

1. the nodal Level Set values are the design variables and the sensitivity analysis is made upon these variables. This is the approach proposed by Belytschko *et al.* in [2];

2. the design variables are also the nodal Level Set values, but the sensitivity analysis is realised by a Frechet derivative : approach followed by Allaire *et al.* [6], Wang *et al.* [4].

In this first work, the evolution equation of the Level Set (eq. 8) is not used as we do not consider the nodal Level Set values as design variables but the parameters of the functions defining the Level Set. Thus, for example, the computation of a Level Set representing a circular or an elliptic void is realised by computing the signed distance between each mesh nodes to a cylinder. Hence, all the nodes have non-zero Level Set value except the nodes located on the cylinder. The design variables in this case can be the radius, the slope or the orientation of the cylinder. This method is convenient in order to construct and combine very basic Level Sets and allow us to modify the shape of the voids. However, it presents a certain limitation as we do not let the Level Set evolve freely but we limit it to a family of shape.

## 6.1. Algorithm

For realising this first study, we have used the algorithms available in the open optimization tool BOSS QUATTRO from SAMCEF. This one has several kind of different general algorithms and also some algorithms specially developed for the mechanical structures problem (Conlin [9], MDQA [10], GCM [8]). In our test cases, the GCM method has shown to be the most adapted and robust in the application treated in this paper. This algorithm has been chosen for the following applications, and the sensibilities are computed using a simple finite difference scheme. The following test cases (Section 6.3) will shown that this method does not introduce any difficulties. However, it limits the number of optimization variables as a full XFEM analysis is required for each variables per optimization step.

## 6.2. Objective functions and constraints

At this moment, all the constraints and the objective functions generally used in the classical shape optimization are available :

1. the potential energy;

2. the displacement;

3. the volume;

4. all stress components.
However it is only possible to take into account a scalar objective functions or a scalar constraints. This restriction is directly linked to the XFEM method. In fact, in the sentivity computation step, some nodes can be removed (appear) and so introduce a smaller (bigger) problem size. This increase leads to a variation in the dimension of the different results which cause some problems in the finite difference procedure.

### 6.3. Test cases

The first problem considered is a rather academic 2D case. A square plate with a traction free elliptic hole under a biaxial tension is considered : $\sigma=2\sigma_0$ on the x direction and $\sigma=\sigma_0$ on the y direction. The dimensions of the plate are 2 x 2 x 1, and the material properties associated are : Young modulus $E=1$, Poisson's ratio $\nu=0.3$. The plain stress state is assumed and the material is elastic. On the initial position, the ellipse is defined by a angle $\theta=45^\circ$, a great axis $a=0.5$ and a small axis $b=0.4$. The optimization objective is to minimize the potential energy of the structure with a constraint on the total surface of the plate. The variables are the angle $\theta$ and the long axis $a$. The discretization used is a Delaunay mesh with 21 nodes on each side.

<figure><img src="image_2.png" /><figcaption>Figure 2: Mesh used for the plate test case</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 3: Evolution of the objective function</figcaption></figure>

We can first notice that the approximation of the Level Set on the mesh gives a quite smooth boundary even with a normal mesh. We can see on the following figure (Fig. 3) that the optimal solution is reached after 11 iterations. These results are in accordance with the solution obtained in the classical shape optimization framework as compared on the table (Table 1.). However, a small difference appears between the solutions obtained with the two methods.
<table><caption>Table 1: optimization report XFEM - FEM</caption><thead><tr><th></th><th></th><th>Iteration 1</th><th>Iteration 11</th><th>Iteration 1</th><th>Iteration 9</th></tr></thead><tbody><tr><td>Objective function</td><td>Minimise <i>U</i></td><td>27.9</td><td>20.2</td><td>26</td><td>18.3</td></tr><tr><td>Constraint</td><td>Surface &lt; 3.45</td><td>3.59</td><td>3.45</td><td>3.50</td><td>3.45</td></tr><tr><td>Variable</td><td>10<sup>-4</sup> &lt; &theta; &lt; 90</td><td>45</td><td>10<sup>-4</sup></td><td>45</td><td>0</td></tr><tr><td>Variable</td><td>10<sup>-4</sup> &lt; a &lt; 1</td><td>0.5</td><td>1.06</td><td>0.5</td><td>0.88</td></tr></tbody></table>

Inside the cutted elements, the Level Set is interpolated linearly when we use first order finite element. As a consequence, the representation of the boundaries can overestimate or underestimate the surface (volume) of the structure. In our test case, this overestimation of the surface leads to an XFEM optimization problem slightly different from the FEM problem (more restrictive on the XFEM). Now, if the two objective functions are corrected in order to take into account this difference, the results are perfectly similar.

<figure><img src="image_4.png" /><figcaption>Figure 4: Approximation of the interfaces inside an element</figcaption></figure>

However, the approximation of the Level Set boundary as a piecewise linear segment introduce another problem. The estimation of the surface depends on the number of cutted elements and on the position of the interface inside the cutted element (Fig. 4). To illustrate this phenomenon, we have computed the variation of the surface with respect to the position of the hole on a square plate (Fig. 5).

<figure><img src="image_5.png" /><figcaption>Figure 5: Evolution of the surface</figcaption></figure>

We can easily conclude that optimization with an objective function on the volume can cause some serious problems.

The second test case is more topology oriented. The structure is a square plate with two overlapping
circular holes under uniform biaxial sollicitation. The material is the same as the previous test case, and the plane stress state is also assumed. The optimization objective is here to minimize the potential energy of the structure with a lower bound on the surface. The design variables are the position of the two circles under the *x* direction. These two variables are bounded to the mesh domain and an lower bound is imposed on the total surface of the structure. The goal is to clearly show that the method is able to merge two circles without any problems. We can already remark on the following figures (Fig. 6) that a fine mesh is necessary in order to approximate with a rather good precision the point where the two holes are joining. This limitation comes from the fact that, inside an element, the interfaces can only be modelled by one segment. Thus, no singular point is accepted inside an element.

<table><caption>Table 2: optimization report</caption><thead><tr><th></th><th></th><th>Iteration 1</th><th>Iteration 12</th></tr></thead><tbody><tr><td>Objective function</td><td>Minimise <i>U</i></td><td>26,6</td><td>14,9</td></tr><tr><td>Constraint</td><td>Surface> 7.8</td><td>6.9</td><td>7.95</td></tr><tr><td>Variable</td><td>-0.5&lt;<i>x</i><sub>1</sub> position&lt;0.5</td><td>0.5</td><td>-0.0662076</td></tr><tr><td>Variable</td><td>-0.5&lt;<i>x</i><sub>2</sub> position&lt;0.5</td><td>-0.5</td><td>0.0457915</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>Figure 6: Evolution of the shape</figcaption></figure>

Here again, the optimal solution is reached in a few number of iterations. As expected, the two holes are placed in a configuration near to perfectly overlap. We can observe that the positioning is not perfect, this error is due to the precision of the XFEM analysis and the Level Set representation. However, this example shows the potential of this method and its capabilities of holding several merging interfaces.

<figure><img src="image_7.png" /><figcaption>Figure 7: Objective function</figcaption></figure>
## 7. Application with genetic algorithms

In this section, we revisit rapidly a paper by Missoum et al [11] in which topology optimization were carried out by an optimal selection of holes characteristics such as size or position with an genetic algo-rithm.

Most of simple optimal design (mitchell-type structures, mbb beam, ...) present a truss structure. Following this observation, the optimization scheme is to place a set holes allowing a variety of truss configurations and let the genetic algorithm choose the best configuration. In our first implementations, the holes characteristics are limited to the presence or not of an entity. This restriction could be extended to the position and the size of an entity. However, the combination of different natures of variables such as geometric characteristics and presence of holes needs a better tuning of the genetic algorithm.

### 7.1 Test case

We have chosen to realise the well kown problem of the mbb beam. The dimension of this one is 18 × 3 × 1, and the material is again an elastic material with Young modulus of 1 and Poisson's ratio of 0.3. The structure is loaded with a unit force on the y direction situated on the middle of the beam. The objective function used is to minimize the compliance with a constraint on the total volume of the beam. The initial configuration chosen is a full mbb beam where no holes are present in the structure and the constraint is included into a multiobjective function.

Regarding to the solution obtained by topology optimization, 8 equilateral triangles are placed on the desgin domain (9 × 3 × 1) at fixed positions : x=0,3,6,9 and y=1.5. Four of these triangles lie on their edge whereas the four others lie on a vertice. The goal is to obtain a solution as close as possible to the topolgy solution.

<figure><img src="image_8.png" /><figcaption>Figure 8: Evolution of the objective function</figcaption></figure>

The convergence is obtained after 36 iterations (Fig. 8) which is very acceptable and related to the number of variables. As expected, the structure obtained is very similar to the topology optimization solution. However, three differences appear between the solutions :

* in our solution, all the truss sections have the same thickness whereas it is not the case in topology;

* the left and right top corner of the structure present some material, not in the topology case;

* the orientation of the triangle's edges are always the same in our case.

These three differences could have been naturally expected as no size variables allowing the optimization of the truss section were present. This simple example shows also that many attentions have to be done when some parts of the structure become very thin as the XFEM analysis may become very inaccurate
<figure><img src="image_9.png" /><figcaption>Figure 9: Comparison between topology optimization and XFEM</figcaption></figure>

if the number of elements per section is very low. The mesh refinement, which would increase the computation time, seems to be the only remedy to this problem.

## 8. Conclusions and perspectives

In this paper, we have presented a new method based on the Level Set method and the XFEM for shape optimization of structures. No important problems have been encountered in the optimization procedure and quite good results were obtained with a simple finite difference scheme. We have shown that the method can exhibit really promising characteristics as it allows topological changes. The XFEM method has proven to be very useful : no remeshing process is needed in our applications and void is not approximated as a smooth material in opposition to the SIMP method.

Now, the main limitations present in our method are the sensitivity analysis and the use of rather simple and basic Level Sets for the representation of the holes. Our ongoing work will be to extend the results already in hand. First, the sensitivity analysis have to be replaced by a more powerful semi-analytic method, this will allow us to study problems involving much more variables. The next step will certainly focus on the development of a "real" Level Set framework in order to allow much more performance in the sense of geometrical modification. Moreover, this development could also improve the calculation of the geometric variables with a much more accuracy and could even remove the oscillation observed in the XFEM computation. The implementation of higher order extended finite element may also reduce this phenomenon as the approximation of the interfaces will hold on more points inside an element including thus a curvature. We also have to study the case of structures becoming very thin as in the mbb beam example.

Clearly, the coupling of the XFEM with the Level Set has proven to be a really promising method for the shape optimization but also for the topology optimization.

## 9. Acknowledgements

The authors are grateful for the research support of the French community through project ARC (Action de recherche concertée) 03/08-298. We also would like to thank Nicolas Moës who allowed us to use the XFEM code from the ECN.

## 10. References

[1] N. Moës, J. Dolbow, N. Sukumar, T. Belytschko, A finite element method for crack growth without remeshing, *International journal for numerical methods in engineering* 1999, no. 46, 131-150.

[2] T. Belytschko, S. P. Xiao, C. Parimi, Topology optimization with implicit functions and regularization. *International Journal for numerical methods in engineering* 57, 2003, 1177-1196.

[3] J. A. Sethian, A. Wiegmann, Structural boundary design via level set and immersed interface methods.
ods, *J. Compt. Phys.* (1999)

[4] M. Yu Wang, X. Wang, D. Guo, Structural Shape and Topology Optimization in a Level Set Based Implicit Moving Boundary Framework

[5] F. Jouve, G. Allaire, A. M. Toader, optimisation de forme par la mthode des lignes de niveaux.

[6] G. Allaire, F. Jouve, A. M. Toader, Structural optimization using sensitivity analysis and a level-set method.

[7] J. Sethian, Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision and Materials Science, *Cambridge University Press*, 1999.

[8] M. Bruyneel, P. Duysinx, C. Fleury, A family of MMA approximations for structural optimization, *Structural & Multidisciplinary Optimization* (2002), no. 24, 263-276

[9] C. Fleury, V. Braibant, Structural optimization : a new dual method using mixed variables, *Internal journal for numerical methods in engineering* (1986), no. 23, 409-428

[10] W. H. Zhang, C. Fleury, A modification of convex approximation methods for structural optimization, *Computers & Structures* (1997), no. 64, 89-95

[11] S. Missoum, Z. Gürdal, P. Hernandez, J. Guillot, A Genetic Algorithm Based Topology Tool for Continuum Structures. *AIAA-2000-4943*.

[12] N. Sukumar, D. L. Chopp, N. Noës, T. Belytschko, Modelling holes and inclusions by level set in the extended finite element method. *Comput. Methods Appl. Mech. Engrg. 190 (2001) 6183-6200*