Acknowledgements

I would like to thank to all the people that in one way or other have helped me to complete my PhD.

First of all, I would like to express my gratitude to my supervisor Prof. Xavier Oliver for his supportfi and availability from the beginning and for his inexhaustible energy devoted to achieving the results of this dissertation. Today, I can say that if I feel the profession of scientist as my own is undoubtedly because of him.

I really appreciate the confidence and patience of my co-advisor Juan Carlos Cante. I would like to thank all the fruitful scientific and personal discussions. I admire his ability to understand.

I would like also to express my gratitude to Alfredo Huespe. Not only for his availability and patience but also for his way on dealing with research. The enthusiasm with which you face new topics, at your age, is my role model.

My sincere thanks also go to Sebastián Giusti, who contributed to achieve the results of this dissertation. On top of a colleague, I get from this experience a friend.

I would like to sincerely thank Prof. Samuel Amstutz for his kind treatment in my stay in Avignon. All our discussions undoubtedly contributed to increasing my confidence and expertise in the topic. I personally believe that our collaboration will bring significant contribution in the topic in the near future.

I would like to thank the Spanish Ministry of Economy and Competitiveness for the FPI fellowship I has granted during my PhD. The research leading to these results has received also funding from the European Research Council under the European Unions Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement n. 320815 (ERC Advanced Grant Project Advanced tools for computational design of engineering materials COMP-DES-MAT).

I would like to express my gratitude to the Universitat Politècnica de Catalunya (UPC-BarcelonaTech), in which I could develop all my PhD studies. I would also thank CIMNE, where I feel at home. In addition, I would like to thanks ESEIAAT (UPC-BarcelonaTech) for giving me the opportunity of exploring my teaching vocation. It has been one of the most rewarding experiences.

My heartfelt gratitude to Chiara for all the conversations and for his readiness to help. My sincerely thank to Anna for how well she has behaved with me from the first moment. To my officemates throughout these years: Stefano, Lucia, Manuel, Vicent, Emmanuel and Marcelo for creating such an enjoyable environment inside and outside the office. To Joan, Ernesto and Arnau for all the conversations we daily enjoyed at lunch time. To Fermín and Ester for all our coffee conversations, for sharing all our PhD lamentations and hardships and for making this period a wonderful experience. You guys are awesome!

I would like to sincerely thank Celia, for her support throughout these years. It has been an incredible period, I will always bring you inside me.

I would like to give an special thank to my father, for his small gestures that show an enormous love. To my mother and Luis and the way they raised me, stood by me, and opened me to the world. To my sisters, Claudia, Julia and Marta for being there on my side and for the emotional support you provided to me.


Abstract

The present dissertation aims at addressing multiscale topology optimization problems. For this purpose, the concept of topology derivative in conjunction with the computational homogenization method is considered.

In this study, the topological derivative algorithm, which is non standard in topology optimization, and the optimality conditions are first introduced in order to a provide a better insight. Then, a precise treatment of the interface elements is proposed to reduce the numerical instabilities and the time-consuming computations that appear when using the topological derivative algorithm. The resulting strategy is examined and compared with current methodologies collected in the literature by means of some numerical tests of different nature.

Then, a closed formula of the anisotropic topological derivative is obtained by solving analytically the exterior elastic problem. To this aim, complex variable theory and symbolic computation are considered. The resulting expression is validated through some numerical tests. In addition, different anisotropic topology optimization problems are solved to show the macroscopic topological implications of considering anisotropic materials.

Finally, the two-scale topology optimization problem is tackled. As a first approach, an structural stiffness increase is achieved by considering the microscopic topologies as design variables of the problem. An alternate direction algorithm is proposed to address the high non-linearity of the problem. In addition, to mitigate the unaffordable time-consuming computations, a reduction technique is presented by means of pre-computing the optimal microscopic topologies in a computational material catalogue. As an extension of the first approach, besides designing the microscopic topologies, the macroscopic topology is also considered as design variables, leading to even more optimal solutions. In addition, the proposed algorithms are modified in order to obtain manufacturable optimal designs. Two-scale topology optimization examples display the potential of the proposed methodology

Resum

Aquest treball té com a objectiu abordar els problemes d'optimització de topologia de múltiples escales. Amb aquesta finalitat, es consideren els conceptes de derivada topologica juntament amb el mètode d'homogeneïtzació computacional.

En aquest estudi, es presenten primer les condicions d'optimalitat i l'algorisme d'optimització utilitzant la derivada topològica. A continuació, es proposa un tractament més precís dels elements de la interfície per reduir les inestabilitats numèriques i els elevats costos computacionals que apareixen quan s'utilitza l'algorisme de la derivada topològica. L'estratègia resultant s'examina i es compara amb les metodologies actuals, que es poden trobar sovint recollides a la literatura, mitjançant alguns assaigs numèriques.

A més, s'obté una fórmula tancada de la derivada topològica anisotròpica quan es resol analíticament el problema exterior d'elasticitat. Per aconseguir-ho, es considera la teoria de variable complexa i la computació simbòlica. L'expressió resultant es valida a través d'algunes proves numèriques. A més, es resolen diferents problemes d'optimització topològica anisotròpica per mostrar les implicacions de la topològia macroscòpica a considerar materials anisòtrops.

Finalment, s'aborda els problemes d'optimització topològica en dues escales. Com a primera estratègia, es consideren les topologies microestructurals com a variables de disseny del problema per obtenir un augment de la rigidesa de l'estructura. Es proposa un algoritme de direcció alternada per fer front a l'alta no linealitat del problema. A més, per mitigar els elevats càlculs computacionals, es presenta una tècnica de reducció per mitjà d'un precalcul de les topologies microestructural òptimes, que posteriorment són recollides en un catàleg de materials. Com a una extensió de la primera estratègia, a més del disseny de les topologies microestructurals, també es considera la topologia macroscòpica com una variable de disseny, obtenint així solucions encara més òptimes. A més, els algoritmes proposats es modifiquen per tal d'obtenir dissenys que puguin ser posteriorment fabricats. Alguns exemples numèrics d'optimització topològica en dues escales mostren el potencial de la metodologia proposada.

1 Introduction

1.1 Motivation

Topology optimization is an emerging field that aims to automate the design process of any structural domain. It seeks to propose optimal topological designs by means of the most leading computational tools. Certainly, topology optimization attempts to complement the design engineer work rather than replace it. On the one hand, due to the knowledge intensively developed in the last years, it can provide designs that offer equal, or even greater, performances. On the other hand, it presents optimal topologies that are often far from being intuitive. It contributes to expanding our creative design capabilities taking us to places often inaccessible to our mind. Admittedly, the possibility of designing complex topologies may seem to be unrealistic to manufacture. However, owing to the recent additive manufacturing techniques, they can be nowadays afforded in a reasonable time and cost.

As an exercise to exhibit the topological optimization scope, one could pose the following question: from a full-domain object, under certain loads and boundary conditions, which would be the best removing material strategy without undermining the structural response capacity? Topology optimization deals with giving the response to this question. In this sense, topology optimization problem can be properly addressed following the motivation quote of work [1]:

“The art of structure is where to put the holes”

Robert Le Ricolais, 1894-1977

Undoubtedly, material reduction is of particular interest in automotive and aeronautical industry. In the former, a decrease on the final structural mass results in a significant economic saving. In the latter, even more evidently, a decrease on the structural mass entails a considerable reduction on the fuel consumption. In order to present an order of magnitude of the economic saving impact, reference [Ipad] asserts that

American Airlines expected to save million a year replacing flight bags by 0.7 kg iPads.

In other words, a of Airbus-A320 structural weight reduction leads to a million saving per year. In addition, from the environment point of view, the fuel consumption reduction entails a considerable decrease on the emission impact.

Certainly, apart from designing the topology, weight reduction can be achieved by means of other techniques. The use of composite materials stands for a first example. These kind of materials, usually modeled by multiscale techniques, provide a reduced weight by appropriately arrange the microstructure with a stiff but heavy material (fibers) and a lighter but softer material (matrix). At this point, in terms of arrangement of fiber-matrix, is it possible to propose novel optimal configurations? Or even more stimulating, is it possible to devise optimal micro-structures by properly putting holes on the stiff material? It naturally leads us to wonder if, on top of designing the macroscopic topology, when the microscopic topology in each macroscopic point is also designed (by means of computational multiscale and topology optimization techniques), an outstanding impact on the mass reduction is achieved.

1.2 Goals

The aim of this study is to address multi-scale topology optimization problems. Under this perspective, through the computational homogenization and topological derivative techniques, the main goal consists in developing the necessary numerical tools to achieve a reduction of the cost function by designing the microscopic and/or macroscopic topologies.

Regarding the optimization problem, at the very onset, a robust and efficient strategy must be established for solving the topology optimization problem when considering the topological derivative. The strategy must intend, on the one side, to mitigate the spurious local minima that appear in the line-search method, and on the other side, to reduce the time-consuming re-meshing techniques.

Apart from the algorithm, computing the closed form of the anisotropic topological derivative yields crucial in this study to achieve the desired results. Since the homogenization of the constitutive tensor of a microscopic topology confers, in general, macroscopic anisotropic response, the current isotropic topological derivative must be extended to anisotropic materials. It stands as one of the main objective of the study. In addition, it is intended to examine the difference between the optimal macroscopic topologies in terms of using either isotropic or anisotropic materials.

As a final goal, this study aims at proposing algorithms and appropiate numerical strategies to significantly decrease the cost function when designing microscopic topologies. Similarly, it is intended to obtain reductions in the cost function not only by designing microscopic topologies but also by designing both simultaneous macroscopic and microscopic topologies. This objective will also result in efficient strategies to tackle the unaffordable multiscale topology optimization problem. In addition, since this study has a practical purpose, developing optimal manufacturable multiscale topologies is the last goal.

1.3 Outline

This dissertation is organized as follows:

Chapter 2

The state of the art of the multiscale method is first described. In addition, the background of the computational homogenization method is presented by the variational multiscale framework including the Hill-Mandel principle, the boundary conditions and the homogenization of the constitutive tensor. Then, in a similar fashion, a brief revision of the topological optimization state of the art is described. The non-existence of solutions and numerical instabilities is addressed. In addition, a concise summary of the different methodologies to tackle topology optimization problems is presented, including the SIMP, shape optimization and topology optimization methods. The latter is devoted in more detail since it represents a core element of this study.

Chapter 3

This chapter is devoted to present the necessary numerical tools to address topological optimization problems when using topological derivative. First, an intuitive and mathematical description of the topological derivative concept is introduced. Then, the topological derivative for the two most relevant shape functions are examined. The optimality conditions in general and tailored to the use of a level-set function, are further explained. In addition, the Slerp algorithm, in the case of equality and inequality constraints, is pointed out. On top of that, a novel numerical technique to treat with the interface in these problems is then proposed and compared with the ones collected in the literature. Some numerical examples account for the proposed interface numerical technique.

Chapter 4

This chapter embraces a full analytical computation of the closed anisotropic topological derivative expression. First, the formulation of the problem is posed and the topological derivative is stated. A summary of the necessary steps to solve analytically the crucial isotropic and anisotropic exterior problems is presented. Full details are collected in Appendices 7 and 8. In addition, the remainders are estimated and the topological derivative is numerically validated. Finally, a wide number of numerical tests are computed for homogeneous and heterogeneous anisotropic materials.

Chapter 5

This chapter concerns multi-scale topology optimization problems. First, the stiffness of the structure is intended to be increased by means of designing the microstructure in each macroscopic point. For this purpose, adequate algorithms and reduction techniques are proposed and validated by some numerical examples. Likewise, a similar approach is proposed to fulfill manufacturability requirements and additional numerical examples are computed. On top of that, a two-scale topology optimization problems is then addressed. An extension of the material design strategies are proposed to achieve the desired results. Some numerical results exhibits the capability of the presented strategies.

Chapter 6

The conclusions of the study are collected in this last chapter. The achievements are first pointed out. Then, the chapter is focused in drawing the final conclusions and outlining the future work lines.

Note that, although a motivation section is devoted in this chapter, an specific motivation section is included at the beginning of each chapter so that each chapter lends itself to become self-contained. Likewise, a conclusion section, related to the specific content of each chapter, is presented.

1.4 Research dissemination

The work developed in this research gives rise to the following scientific publications:

A Ferrer, J Oliver JC Cante, and JA Hernández. Two-scale topology optimization: an efficient integrated structural optimization and material design approach. Draft, 2016.

'SM Giusti, A Ferrer and J Oliver. Topological sensitivity analysis in heterogeneous anisotropic elasticity problem. Theoretical and computational aspects. Computer Methods in Applied Mechanics and Engineering, 2016. http://dx.doi.org/10.1016/j.cma.2016.08.004

A Ferrer, J Oliver JC Cante, and O Lloberas. Vademecum-based approach to multi-scale topological material design. Advanced Modeling and Simulation in Engineering Sciences, 2016. doi:10.1186/s40323-016-0078-4

Additionally, the work has been presented in the following conferences and workshops:

JC Cante, A Ferrer, J Oliver. Numerical tools for Multi-scale material design and structural topology optimization. In ECCOMAS VII European Congress on Computational Methods in Applied Sciences and Engineering, Creta, 2016.

A Ferrer, J Oliver, JC Cante. Multi-scale topological design: a Vademecum-based approach. In New Challenges in Computational Mechanics (A Conference Celebration the 70th Birthday of Pierre Ladevèze), Paris, 2016.

J Oliver, JC Cante, A Ferrer. Computational design of engineering materials: An integrated multi-scale approach to structural topological optimization. In XIII International Conference on Computational Plasticity. Fundamentals and Applications. COMPLAS 2015, Barcelona 2015.

A Ferrer, JC Cante, J Oliver. An efficient tool for multi-scale material design and structural topology optimization. In XIII International Conference on Computational Plasticity. Fundamentals and Applications. COMPLAS 2015, Barcelona, 2015.

A Ferrer, JC Cante, J Oliver. On Multi-scale structural topology optimization and material design. In CMN-2016: Congress on numerical methods in Engineering, Lisbon, 2015.

A Ferrer, JC Cante, J Oliver. Towards real time analysis in multi-scale computational design of engineering materials. In CSMA-SEMNI Numerical techniques for computation speedup, Biarritz, 2015.

A Ferrer, J Oliver, A Huespe, JA Hernandez, JC Cante. On macrostructure and microstructure optimization techniques in multiscale computational material design. In 11th World Congress on Computational Mechanics, Barcelona 2014.

A Huespe, J Oliver, A Ferrer, A Huespe, JA Hernandez, JC Cante. Hierarchical multiscale optimization of microstructure arrangement and macroscopic topology in computational material design. In XII International Conference on Computational Plasticity. Fundamentals and Applications. COMPLAS 2013, Barcelona, 2013.

Finally, the work of this dissemination has been complemented by the developments achieved in the following research stay:

Université dAvignon, 4-month doctoral research stay. Worked under the supervision of Prof. Samuel Amstutz in the Laboratoire de Mathématiques dAvignon, Avignon, France. May-September 2016.


2 Background and review of the state of the art

Since the aim of this work is to solve multiscale topology optimization problems, this chapter is focused on describing the bases of the multi-scale and topology optimization techniques. On the one hand, the theoretical background is briefly summarized. On the other hand, as a state-of-the-art review, different theories are presented and their major advantages and disadvantages are highlighted.

First, the computational homogenization theory is introduced by detailing the foundations of the multi-scale variational framework, the Hill-Mandel Principle, the equilibrium and boundary conditions and finally the necessary steps for computing the constitutive tensor.

Second, the mathematical foundations of the the topology optimization problem are presented. The non-existence results and the numerical instabilities of that usually appear in topology optimization problems are . We also describe the most used topology optimization techniques, including the SIMP methodology and the shape derivative approach, and finally as a core element of this thesis, we feature the topological derivative for topology optimization problems applied to the macroscopic and microscopic domain.

2.1 Computational homogenization (FE)

In continuum mechanics, it is necessary to define the dependency between the stresses and strains in order to solve the standard problem of solid mechanics.

Usually, the main difficulty of the constitutive law ( relation) lies on how to model the non-linear behavior of the material. However, this is not the case of this work. Since the primary objective is to design materials and structures, we consider, throughout all this work, linear elasticity. That is, the stress tensor depends linearly on the deformation tensor at each point. Namely,

(2.1)

where is the fourth order constitutive tensor.

Apart from the non-linearity of the material, it is worth wondering if the constitutive law provides enough information on the material behavior. The answer would depend on the accuracy that we require. Since precisely we focus on designing materials, a high level accuracy on the constitutive law (through the relation) is required. Thus, the technique for setting the constitutive law must be powerful enough to provide an accurate relation. The existing approaches for modeling the constitutive law can be basically grouped in two currents, phenomenological techniques and multi-scale homogenization techniques.

In some applications, phenomenological constitutive laws are powerful enough to model the material and to define the relation. However, in highly demanding applications, it is necessary to make use of more sophisticated techniques such as the multi-scale homogenization method, in order to capture the complexity of the materials [2]. Phenomenological constitutive approaches are only able to capture these small variations, in the Finite Element (FE) context, with unaffordable fine meshes. By contrast, the computational homogenization method, usually called , through the definition of two different scales, is able to capture such small variations with a reasonable computational effort.

Regarding the computational homogenization method, in the last decades, it has gained considerable popularity in the computational mechanics community. Admittedly, providing an accurate relation with a reasonable computational effort represents a significant advantage. In addition, and more significantly, since the approach basically requires standard Finite Element (FE) techniques, the computational homogenization method suits naturally in the computational mechanics context from the formulation point of view and the implementation point of view.

In order to set up the corresponding mathematical framework of the computational homogenization method, different theories have been developed in the last years [3,4]. Apparently, asymptotic expansion and variational multi-scale approach are, nowadays, the most successful approaches in the context of computational mechanics. Although the asymptotic expansion approach is a rigorous mathematical theory and it has been used for a long time, the variational multi-scale theory seems more appropiate to extend to non-linear problems. Furthermore, variational approaches usually fit more naturally in the context of the Finite Element method.

The main concepts of the computational homogenization method, which holds for both the asymptotic expansion approach and variational multi-scale method, are first described. Further ahead, we describe the essential concepts of variational multi-scale method that are needed for achieving the results of all this work.

Since the computational homogenization method aims at considering small heterogeneities and small variations of the variables, it proposes to define, firstly, a macroscopic scale (characterized by the length scale ) which corresponds usually with the length of the domain and, secondly, a microscopic scale (characterized by the length scale ), which typically is of a smaller order of magnitude. It is assumed that the microscopic scale fulfills

(2.2)

Accordingly, one can define a macroscopic coordinate (macroscopic point) related to the macroscopic scale and a microscopic coordinate (microscopic point) related with its counterpart scale . The parameter

(2.3)

usually measures the jump between the scales. Note that if the coordinate is neglected, the standard one scale problems is recovered.

Thus, with the definition of these two different scales in mind, and with the idea of considering heterogeneities on the small scale, the constitutive tensor is modeled as a function of both macroscopic coordinate and the microscopic coordinate as

(2.4)

Since the variables (stresses , strains ) of a standard elasticity problem depend implicitly on the constitutive tensor through the equilibrium equation, in principle, they are also function of both macroscopic coordinate and the microscopic coordinate , that is

(2.5)

Conceptually, the main idea of the computational homogenization method is to collect the variation of the variables with respect the microscopic coordinate through an homogenization process. On the one hand, that allows capturing heterogeneities of the microscopic scale. On the other hand, after applying the homogenization process, the standard variables (stresses , strains ) of a macroscopic continuum mechanical problem may be retrieved. In mathematical terms, the explicit -dependence of the variables disappears after applying the homogenization process, this is

(2.6)

where corresponds to the homogenization constitutive law. In order to describe this homogenization process, the RVE (Representative Volume Element) concept is first introduced. It is usually defined as the microscopic domain (of order of magnitude ) in which the variations of the material properties are sufficiently representative. That allows associating the macroscopic variable to the coordinates of the macroscopic domain and the microscopic variable to the coordinates of the microscopic domain . When the jump between the scales is large enough, each integration/sampling point of the continuum macroscopic domain is associated to an RVE. A sketch of this concept is presented in Figure 1.

Representation of the macroscopic domain Ω and the microscopic domain Ωμ. In each macroscopic coordinate x, the associated RVE (Representative Volume Element) allows considering heterogeneities on the micro-scale through the microscopic coordinate y.
Figure 1: Representation of the macroscopic domain and the microscopic domain . In each macroscopic coordinate , the associated RVE (Representative Volume Element) allows considering heterogeneities on the micro-scale through the microscopic coordinate .

Due to the definition of both coordinates and , the heterogeneities of the material on the macroscopic domain and on the microscopic domain can be considered. At this point, the computational homogenization method proposes to replace the heterogeneous microscopic domain by an equivalent homogeneous microscopic domain, hence its name. See, in Figure 2, an sketch of this concept.

Computational homogenization representation. The variational multi-scale technique considers the heterogeneities on the microscopic scale through the coordinate y. However, after applying the homogenization process, the macroscopic variables σ(x), ɛ(x) and \mathbbCh(x) come to depend only to the macroscopic coordinate x.
Figure 2: Computational homogenization representation. The variational multi-scale technique considers the heterogeneities on the microscopic scale through the coordinate . However, after applying the homogenization process, the macroscopic variables , and come to depend only to the macroscopic coordinate .

The methodology of replacing the heterogeneous RVE by its homogeneous counterpart is what the variational multiscale method proposes. Note that, although the heterogeneous micro-scale becomes an homogeneous material, the macroscopic material can still be considered an heterogeneous material, i.e., it no longer depend on variable but it may still depend on variable .

As a first attempt of this homogenization process, one could think naturally on using the rule of mixtures. However, more sophisticated approaches may be employed [5], for instance, the multi-scale variational framework.

2.1.1 Multi-scale variational framework

The multi-scale variational approach has been extensively used in the last years [4] as a mathematical framework for the computational homogenization. On the one hand, it takes use of the powerful tools of calculus of variations and, on the other hand, it develops and presents the concepts nimbly. The methodology is based on: first, defining the kinematics, second, taking variations on the strain space of functions and, finally, postulating the Hill-Mandel principle. Henceforth, the variables related with the micro-scale are endowed by the sub-index .

Regarding the kinematics, the multi-scale variational approach asserts that the microscopic strain can be written as the sum of two terms, a macroscopic strain and a fluctuation strain , i.e,

(2.7)

The macroscopic strain is defined commonly as , where represents the displacements, and it depends only on the macroscopic coordinate , while the fluctuation strain depends on both coordinates and and must fulfill the constraint of zero mean value over the microscopic domain , that is,

(2.8)

Thus, both the splitting on two terms of the strain and the zero mean value constraint on the fluctuation strain are considered axioms of the multi-scale variational approach. Note that, this choice allows us to guarantee that the average of the microscopic strain will be the macroscopic strain, more specifically,

(2.9)

This relation corresponds to the homogenization process shown in Figure 2 applied to the strain field. Thus, the microscopic strain , which takes values over all the domain and depends on both and coordinates, is defined as the sum of its mean value over the microscopic domain and a fluctuation part, which is in charge of measuring the deviation over such mean value. As a remark, in the case where the fluctuations are canceled, the standard macroscopic problem is recovered.

Once the kinematic is defined, we move to the definition of the space of function of the strains. First, the space of function of the macroscopic and fluctuation strains are defined as

(2.10)
(2.11)

where stands for the symmetric second order tensor spaces. The macroscopic strain are only asked to belong to the symmetric second order spaces while fluctuation strains are additionally asked to satisfy the zero mean value constraint over the microscopic domain. Then, the space for the microscopic strain is defined as

(2.12)

with and .

2.1.2 Hill- Mandel principle

Taking variations in such spaces, the variational multi-scale approach makes use of the Hill-Mandel principle. Based on energy concepts, it postulates that the internal energy of a macroscopic point is equivalent to the average of the microscopic internal energy of the microscopic domain. In physical terms, it states that two different entities (the infinitesimal point of coordinate and the microscopic domain) are equivalent and replaceable if they are endowed with the same value of the internal energy. In mathematical terms, this statement is expressible as

(2.13)

where corresponds to the macroscopic stress and to the microscopic stress. By virtue of (2.12), a variation of the microscopic strain is given by the variation of the macroscopic and fluctuation strain as

(2.14)

Inserting such variation on the Hill-Mandel principle equation (2.13), we obtain the following equation

(2.15)

Since equation (2.15) holds for all , it also holds for . In consequence, the homogenization of the stress appears naturally as,

(2.16)

Clearly, this relation is, in stresses, the analogous equation of the strain equation (2.9). As in the strain case, equation (2.16) corresponds to the homogenization process shown in Figure 2, but in this case, applied to the stress field. Note that the strain and stress homogenization differ basically on how they have been stated, the former as an axiom, the latter as a consequence of the Hill-Mandel principle.

Similarly, since equation (2.15) holds for all , it also holds for . Inserting this condition into equation (2.15), the weak form of the micro-structure equilibrium equation is obtained as

(2.17)

In view of equation (2.17), the fluctuation strain does not produce internal energy on the RVE.

2.1.3 Micro-scale equilibrium equation and boundary conditions

The Hill-Mandel principle leads to solve the equilibrium equation (2.17) for each macroscopic point . This means that, through a discretization in , the macroscopic and a microscopic scale problem (in each macroscopic point) must be solved.

Regarding the micro-scale equilibrium equation, we first express the relation on the micro-scale, i.e, relation. Since the aim of the work is based on material design and structural optimization, we restrict the constitutive law to the elastic regime of materials. Thus, the macroscopic stress depends linearly on the strain through the micro-scale constitutive tensor as

(2.18)

Therefore, the micro-scale equilibrium equation (2.17) is rewritten as

(2.19)

and considering the split of the microscopic strain in equation (2.7), the equilibrium equation results

(2.20)

This last equation stands for the equilibrium equation in strain terms. To write it in terms of displacements, we must first apply the Gauss theorem to the fluctuation strain condition (2.8), that is,

(2.21)

where the fluctuation displacement has been introduced through the standard relation . Similarly, the micro-scale displacement can be defined from . Consequently, integrating the splitting equation (2.7), the micro-scale displacement can be written in terms of the macroscopic strain and the fluctuation displacement as,

(2.22)

Then, we can define the space function of the fluctuation displacement as

(2.23)

Unlike the macroscopic displacement field, the fluctuation field is requested not only to enjoy standard regularity of elliptic problems but also to fulfill condition (2.23).

Thus, the equilibrium equation (2.20) can be re-expressed in terms of the fluctuation displacement as,

(2.24)

This last equation corresponds to the equilibrium equation that must be solved in each micro-structure. Note that, the equilibrium equation (2.24) suggests that the micro-scale equilibrium could be interpreted as a standard macro-scale equilibrium problem where the fluctuation plays the role of the unknown and the macroscopic strain the role, after integration, of the right hand side.

Next step is to describe how the microscopic boundary conditions may be fulfilled. There are different approaches to satisfy these boundary conditions of the RVE. In literature, see [4], the most frequently used can be classified in Taylor, Linear, Periodic and Minimal condition.

Taylor boundary conditions. Frequently, this model is commonly termed, in other contexts, rule of mixtures [5]. Intuitively, it homogenizes the properties by its volumetric contribution. In our terms, it turns into imposing zero fluctuation over all the domain (including the boundary), that is

(2.25)

Thus, the equilibrium equation is not necessary to be solved since the unknown is already known. Obviously, by definition, condition (2.21) fulfills, i.e.,

(2.26)

Linear boundary conditions. In comparison to the fluctuation on the Taylor conditions, the linear boundary conditions are imposed to be zero only on the boundary, i.e,

(2.27)

According to (2.22), the total displacement has, in this case, only the contribution of the macroscopic strain:

(2.28)

Thus, in the Linear boundary conditions, the micro-scale displacement depends linearly on the boundary with respect to coordinate , hence its name. In this case, if the fluctuation is zero over all the boundary, the integral of the symmetric open product between the fluctuation and the normal is also zero, that is,

(2.29)

Since less conditions are imposed, Linear boundary conditions are less stiffer than the Taylor boundary conditions. However, there is still room to impose softer ones.

Periodic boundary conditions. Alternatively, the periodic boundary conditions are the ones with better reputation in the multi-scale field. In the literature, there are numerical studies that suggest its use in periodic material [6]. The main advantage lies on the fact that, the size of the micro-structure in which the material is statistically representative, is the smaller size in comparison to the rest of boundary conditions. Thus, the condition on the jump of scales is easier to satisfy.

The periodic boundary conditions satisfy the fluctuation condition as follows. For some specific micro-scale geometries like square cells (hexagonal cells and others can be easily extended), the boundary is divided in , , and with outward unit normal such that

(2.30)
RVE periodic domain, square cell.
Figure 3: RVE periodic domain, square cell.

The different parts , , and of the boundary are represented in Figure 3. Considering the division on the boundary, we can re-express the fluctuation condition as

(2.31)

where and represents the fluctuations in and . Similarly, and stands for the fluctuations in and .

Finally, considering the opposite direction on the normals defined in Figure 3, we end up with the periodic boundary conditions

(2.32)
(2.33)

More specifically, the fluctuation on the left part of the square cell must be equal to the fluctuation on the right part and, similarly, on the up and bottom part. Physically, this feature permits considering other micro-cell surrounding the RVE, and thus, the fluctuation will vary periodically along the micro-cell, hence its name.

Minimal boundary conditions. The minimum boundary conditions appear as the last alternative. They are considered the weaker boundary conditions since, in contrast to other boundary conditions, they assume no extra condition, hence its name. For this purpose, the fluctuation conditions (2.21) are imposed directly.

Note that fluctuation condition (2.21) leads to impose that the integral over the boundary of the open product between the fluctuation and the outward unit normal is zero. That implies to impose six conditions in 3D problems and three conditions in 2D problems.

Selected boundary condition and implementation strategy. The choice of the boundary condition is a priori arbitrary, and it would depend on the addressed problem. In our study, and throughout all this work, we select the periodic boundary conditions since they can ensure a representative volume with the smaller length scale .

Regarding the way to impose the boundary conditions, there are two options. On the one hand, it can be imposed directly in the equilibrium problem and consequently the Lagrange multipliers appear as extra unknowns. On the other hand, it is possible to condensate some unknowns on the system through the boundary conditions.

The first option seems reasonable for small number of conditions like minimum conditions, however, for other kind of conditions, it enlarges the system of equations considerably. The option of condensing the unknowns in periodic and linear seems reasonable. As a difficulty, if iterative solvers are used, specially for large 3D problems, the condensation process confers worse conditioning to the matrix, lengthening the convergence process. If direct solvers are used, with the condensation technique, the matrix becomes less sparse bringing problems with memory. This features that the appropiate approach will depend on the specific problem to be solved. In our case, and through all this work, since not computationally high demanding meshes are required, the condensation process is considered.

Strong form of the microscopic equilibrium equation. Once the boundary conditions have been introduced, the strong form of the microscopic equilibrium equation can be stated. From the equilibrium equation in weak form (2.17), and undoing the steps of integration by parts (and assuming some regularity), the equilibrium equation requires divergence free of the microscopic stresses. This is,

(2.34)

To extend the formulation of the variational multi-scale method to microscopic equilibrium with body forces, the reader is refereed to work [7]. Thus, the strong form of the equilibrium equation jointly with the microscopic constitutive law (2.18) and the boundary conditions (periodic in our case) form the set of necessary equation of the micro-structural problem. Mathematically, it reads

(2.35)

2.1.4 Homogenized constitutive tensor

Once the microscopic equilibrium equation is presented, it is worth mentioning that the main aim of the multi-scale problem consists in obtaining the homogenized constitutive tensor rather than the microscopic fluctuations, microscopic strains and microscopic stresses. In fact, computing the homogenized constitutive tensor is the essential part of the computational homogenization method. Henceforth, it is denoted by As in other fields of mechanics, it is commonly defined as the variation of the stresses with respect to the strain , this is

(2.36)

Taking into account that the macroscopic stresses are related with its microscopic counterpart through the homogenization equation (2.16), and making use of the linear elasticity assumption (2.18), the constitutive tensor reads

(2.37)

The microscopic strain is related with the macroscopic strain through the assumption of the kinematics described in equation (2.7). By taking derivatives with respect to the macroscopic strain in equation (2.7), we obtain the following relation

(2.38)

where and stand for the identity fourth-order tensor and the so-called fourth-order localization tensor [8]. Unlike the homogenized fourth order constitutive tensor , both fourth order tensors and are dimensionless. Regarding this fourth-order localization tensor , note that, in the weak form of the equilibrium equation (2.24), the symmetric gradient of the fluctuations depends linearly (due to the assumption of linear material behavior) on macroscopic strain . Thus, we can relate with by writing the definition of the localization tensor as follows,

(2.39)

Since the localization tensor stands for a linear operator, its components can be obtained by solving, for each canonical base of [2], the symmetric gradient of the fluctuation from the weak form of the equilibrium equation (2.24). Once the is known, replacing equation (2.38) on the definition of the constitutive tensor, the homogenized constitutive tensor is reduced to

(2.40)

being and the volume average of the microscopic constitutive tensor and the fluctuation constitutive tensor. By definition, they take the following form

(2.41)

Note that when the fluctuations are null (for example with Taylor boundary conditions), the localization tensor is canceled. Consequently, the homogenized constitutive tensor depends only to the volume average of the microscopic constitutive tensor . That explains why is commonly called the Taylor counterpart of the homogenized constitutive tensor. At this point, it is worth stressing that the whole homogenization process presented in this section is basically reduced to the homogenization of the constitutive tensor described in equation (2.40).

Note that all the formulation of the variational multi-scale method has been introduced under elastic regime assumptions. To extend this formulation to non-linear problem the reader is referred to [4].

2.2 Topology optimization

In the last decades, topology optimization has been a wide active research topic. Nowadays, it is widely applied to Aeronautical [9], automotive and civil engineering industry. In addition, topology optimization tools are nowadays included in more than thirty commercial software packages [1], e.g. Abaqus [10], Altair HyperWorks [11]. As an example, the design of the A380 ribs shown in Figure RealRib through topological optimization techniques represents a prominent application of the method in the Aeronautical industry.

Draft Samper 118254298-real rib top opt.png Draft Samper 118254298-real rib top opt2.png
Topology optimization applications on aeronautical industry. The topology optimization fields is useful and developed enough for giving answers to real industry problems.
Figure 4: Topology optimization applications on aeronautical industry. The topology optimization fields is useful and developed enough for giving answers to real industry problems.

During these years of research and industrial development, different theories have been proposed. SIMP, shape optimization and topological derivative method are considered ones of the most convincing approaches.

The arguably most popular method, SIMP [1], is based on an heuristic regularization which leads to an appropriate (in terms of practical results) penalization. The root of this approach can be traced back to the seminal work [12] developed by Kikuchi and Bendsoe. In addition, due to this regularization, gradient-based methods can be used. However, although it can be sometimes interpreted in physical terms (micro-structures), intermediate values still appear and the selection of the penalization parameter is still an open issue [13]. This method has been successfully applied to material design in work [14] and to multi-scale topology optimization problems in the seminal work [15]. A more recent example can be found in work [16].

Boundary variation methods, based on classical shape sensitivity analysis, have appeared as a powerful alternative. Although shape differentiation has been long studied, the bases were established in the reference book [17]. The idea of using a level set method [18] for representing the topology is one of the strengths of the approach. The shape derivative is then used as a descend direction in a Hamilton-Jacobi equation which makes the level set evolves. Moreover, although strong mathematical theories have been developed, the inability of nucleating holes makes this approach limited.

Complementary, in the last years, due to the ever increase of the available computational power, discrete and evolutionary algorithms, like BESO [19], appeared as a clear way of defining interfaces in topology optimization. These methods, based in heuristics, offer the advantage that they are simple to code. However, they are limited to small cases and they are not computationally very efficient. For these reasons, they lack a great reputation in the topological optimization community [13].

In the last years, techniques based on topological derivative has become a significant tool for solving topology optimization problems. Its main merit is to provide sensitivity of the cost function when a small hole is created. The basis of topological derivative theory was first established in [20] by using the shape sensitivity results and then consolidated in the reference book [21]. Works [22] and [23] deserve also special attention. One of its main drawbacks, however, lies on the difficulty of obtaining such topological derivative, which can become in some cases burdensome. In fact, up to now, for some specific cost functions, topological derivatives are still missing. However, large advances have been achieved in the last years.

Owing to the landmark work [24], topological derivative can be used as a descent direction in a level-set algorithm. Slerp (spherical linear interpolation) algorithm has been used as an efficient strategy for solving topology optimization. It provides clear boundary of the topology and, in contrast with shape optimization techniques, the nucleation of the holes appears naturally.

2.2.1 The topology optimization problem

In the following, the topology optimization problem is stated under the assumptions of linear elasticity and small strains. Generally speaking, the aim consists of obtaining an optimal topology such that it minimizes a desired functional and satisfies some particular constraints.

The description of the topology is determined by the characteristic function as follows

(2.42)

where the domain has been split into two parts . The sub-domains and are made of different materials, thus, the characteristic function is in charge of determining in the whole domain what part corresponds to either material. Such kind of problems are normally termed bi-material topology optimization problems. However, in most of the application, instead of dealing with two materials, the material corresponding to the domain is made of a very small stiffness characterizing the behavior of a void. The name of topological optimization problems usually refers to this case.

Regarding the cost functional, hereafter termed as , different options are possible. The compliance is widely used in the context of topology since it measures the stiffness of the structure. Other possibilities found in the literature are, for instance, the least square objective function [25], which attempt to achieve a certain displacement of the structure in norm.

Regarding the constraints, it is common to fix a desired volume of one of the materials, typically the strong one. In addition, it can be found, in works [26] and [27] of the literature, constraints imposing a threshold on the stresses. Perimeter constraints can also be found in order to alleviate numerical instabilities [28].

Concerning our work, we are interested in using the compliance as the cost function and the volume as the constraint. Accordingly, we state the topological optimization problem as follows,

(2.43)

where is a general cost function (normally the compliance), the volume value to achieve, and the characteristic function. The cost function, in the case of the compliance, has the form

(2.44)

where and represents the left hand side and the displacements solution of the following equilibrium equation

(2.45)

where represents the bilinear form obtained by taking the weak form of the standard elastic problem

(2.46)

The field stands for the stresses and stands for the constitutive fourth order tensor, usually defined in the bi-material problem as

(2.47)

or, more explicitly, as

(2.48)

where and are the constitutive tensor of the strong material and the weak material respectively.

2.2.2 Non-existence and numerical instabilities

In the topology optimization problems, some difficulties appear not only on the theoretical aspects but also on the numerical implementation. Regarding the theoretical aspects, we briefly describe the non-existence of optimal solutions by giving a representative example and we comment the standard remedies found in the literature. Regarding the numerical difficulties, we discuss the lack of unicity of the solution and the checkerboard instability.

Lack of existence of solutions

The existence of solution in topology optimization has been largely studied [29]. The question can be stated as: is there an optimal topology that minimizes the cost function and satisfies the constraints? This question is fully addressed [30] in the literature by providing some examples that show the non-existence of optimal solution. In the following, we are going to outline the counter-example described in book [29].

The aim lies on finding a topology that minimizes the compliance of a square domain under unit-axial loads.

Draft Samper 118254298-NonExistence.png Counter-example of non existence of solution for the topology optimization problem. Although the applied loads and the volume of the stiff material are the same, the example on the right is stiffer, i.e., it exhibits smaller compliance. Since there is no limit on increasing the number of inclusions, no optimal solution exists.
Figure 5: Counter-example of non existence of solution for the topology optimization problem. Although the applied loads and the volume of the stiff material are the same, the example on the right is stiffer, i.e., it exhibits smaller compliance. Since there is no limit on increasing the number of inclusions, no optimal solution exists.

Intuitively, the main idea consists in proposing a sequence of topologies with an increasing number of horizontal elliptic inclusions, smaller every time to preserve the volume. It can be seen that such sequence entails smaller compliance values. Since the sequence does not converge, the compliance value may always decrease. Thus, the optimal solution does not exist.

One remedy is to modify the problem by adding constraints. In several cases, it is very common to add a Perimeter constraint. In [28], it is shown how the Perimeter constraint may be used to limit the size of the inclusion. In addition, it can be proved with this Perimeter constraint, the existence of solutions is recovered. See works [31] and [28] for further information. Alternatively, it is also possible to introduce restrictions based on manufacturing issues to limit the size of the inclusions.

Lack of unicity

Besides the problems of existence, there are lack of unicity of solutions. In work [1], it is shown that one can get unicity of solutions if some specific regularizations of the problem are considered. However, most of the times, the obtained topology is full of gray regions and no physical meaning may be retrieved [13]. If this specific regularization is not considered, depending of the initial guess, several local minima may appear. This feature is due to the non-convexity of the problem. Thus, the initial value of the optimization problem influences on the optimal solution.

Checkerboard

In the context of topology optimization, it is well-known that, depending on the methodology, checkerboard solutions may appear. An intuitive way of understanding such phenomenon is writing the topology optimization problem (2.43) as a saddle point problem of two fields: the displacements and the characteristic function , this is

The checkerboard phenomenon usually appears, in other fields (e.g. Stokes flow), on problems. In that case, the Babushka-Brezzi condition must be satisfied. Similarly, in the topology optimization problem, the use of an appropiate interpolation of the displacement field and the characteristic function may circumvent the checkerboard problem. Alternatively, an other remedy for avoiding checkerboards may be the use of filters. A deep study on this topic is collected in work [32]. For further information, the reader is referred to works [32] and [33]. A representative example of the checkerboard phenomena is depicted in the standard Cantilever example of Figure 6.

Checkerboard phenomena in topology optimization problems (figure extracted from [1])
Figure 6: Checkerboard phenomena in topology optimization problems (figure extracted from [1])

2.2.3 Regularized topology optimization (SIMP)

Regarding the different approaches to address the topology optimization problem, we start by describing the SIMP method. As it was mentioned before, the SIMP method is the most common approach in topology optimization. The origin of the methodology comes from the seminal paper of Kikuchi and Bendsoe[34]. It has been applied to many fields with success and it has produced hundreds of publications. Certainly, all the publications of Sigmund and its group has strongly contributed on giving useful solutions to real industry problems. The reference book [1] of Bendsoe and Sigmund, certifies that, nowadays, the SIMP method is considered by the topology optimization community, a consolidated theory.

The approach is based on regularizing the discontinuous characteristic function . Instead of taking values zero or one, the characteristic function is allowed taking intermediate values. Schematically,

where denotes the regularized characteristic function and it is usually termed fictitious density.

In addition, in order to get black and white topologies, the definition on the constitutive tensor in equation (2.48) is modified as

(2.49)

where the parameter is arbitrary and leads to penalize the intermediate fictitious densities. Typically, in linear elasticity, the penalization parameter is taken as However, no penalization is considered in the volume constraint. It reads

(2.50)

The fact that the constitutive tensor is penalized with the power of law and the volume is not penalized leads to stimulate the problem on choosing zero or one values of the fictitious density . Intermediate values are expected to increase the value of the volume with a non strong increase on the stiffness. However, depending on the problem intermediate values could appear with the inconvenient that not clear interpretation can be done. There are attempts to interpret it as an homogenized microstructure made of both materials but not always is possible (since it falls outside the H-S bounds).

The SIMP method, as an advantage, entails a straightforward implementation and provides satisfactory results. In addition, due to the penalization parameter , well-established and powerful continuous optimization algorithms can be used.

However the SIMP method also entails some inconveniences. The value of the heuristic parameter is unclear. In addition, in many cases, checkerboards may appear and filters must be applied. Finally, in the optimal solution, the fictitious density takes frequently intermediate values. As a remedy, there are two different options: either thresholding techniques are applied to recover “black-white” solutions or the intermediate values are interpreted, when possible, in terms of micro-structures.

2.2.4 Shape derivative for topology optimization

Another very popular approach for topology optimization consists in using the shape derivative concept. Many works have been developed in the last decades, specially by the french school [35] , [36] and [25]. The idea lies on collecting in the shape derivative the measure of the change of a functional when applying a deformation on the domain. In the case that, the deformation is applied on a boundary in the normal direction and with unitary modulus, the shape derivative can be straightforwardly computed, see [29]. Shape derivative has been obtained for several functionals like the compliance, volume and least squares cost functions, among others.

After the computation of the shape derivative, the optimality conditions have also been defined in the literature [29]. Since the change of the cost is collected by the shape derivative, it may be used as a descent direction in an optimization algorithm. This idea was pioneered by Allaire in the seminal work [25]. The topology is defined by a level-set function and it evolves following a Hamilton-Jacobi scheme. With this methodology, encouraging results have been achieved.

Since this methodology is based on a level set function, it is a priori free of grays (intermediate values) and no checkerboards are expected to appear. In addition, it relies on a strong mathematical theory.

As a main drawback, the shape optimization problem is restricted, by construction, to boundary changes, and consequently, no new nucleations of holes (or inclusions) are accomplished. However, some remedies has been proposed in [37] by combining the shape derivative and topological derivative.

2.2.5 Topological derivative for topology optimization

The field of topological derivative has gained an increasing popularity in the last decades. Although much of its success lies on the advances developed in the context of the shape derivative, it was not until the work of Solowski & Zochowski in 1999 [20], where the basis of the topological derivative theory were rigorously established. In parallel, Masmoudi presented the basis of the shape and topological derivative in [38]. Special attention deserves the previous works of Schumacher in [39] and [40] and Cea in [41]. Some years later, Novotny and co-workers at [42] presented and established a clear relation between the shape derivative and the topological derivative through the Eshelby tensor and connecting it with the configurational forces fully described in the reference book of Gurtin [43]. As a consequence of that studies, the reference book Topological Derivatives in Shape Optimization [21] emerged.

The topological derivative is based on the asymptotic analysis of shape functionals and the analytical solution of classical elastic problems of an infinite domain with an elliptic inclusion. The studies of Nazarov in [44] helped to consolidate the theory of the asymptotic analysis. Regarding the analytical classical solution of elastic problems, the books of Little [45], Muskhelishvili [46], Lekhnitskii [47] and Saad [48] are classical references.

All the progress on topological derivative theory leads to an extended number of applications which could be summarized in the following list

From the author's point of view, topological derivative is an adequate, and probably, the most natural tool for solving topology optimization problems since it studies the sensibility of a shape functional when making a hole (or inserting an inclusion).

It is worth mentioning that there are two different and complementary procedures to compute the topological derivative. The first one is based on studying the shape functional on two different configurations, with and without and inclusion. Then, after few manipulations, the differences of such functional are related with the topological derivative. The second approach, according to the reference book [66], takes advantage of the shape derivative by making use of standard continuum mechanics tools like the material derivative concept. In this approach, it appears naturally the Eshelby tensor and it has more the flavor of the configurational forces school. This approach helps to understand intuitively the physical meaning of the topological derivative.

A brief introduction to the shape derivative concept is the following (see Figure 7). Let's assume that a circular hole (or inclusion) on the unperturbed domain with a radius of value is nucleated (or inserted). Then, the shape derivative of the desired objective function must by computed by taking the material derivative of the functional with zero velocity field on the boundary of the domain and with unitary and normal velocity field on the circular hole. This shape derivative could be seen as a limit of a small perturbation of the radius of the inclusion (limit in in Figure 7).

Then, once the shape derivative is obtained (which itself involves a limit), we expressed it in terms of , or more specifically, in terms of powers of , mimicking a standard Taylor expansion. It is worth mentioning that the most relevant part of computing the topological derivative remains on obtaining this shape derivative in terms of power. This procedure is well-explained and detailed in chapter 4 and in Appendices 7 and 8, for both isotropic and anisotropic materials. Finally, by taking the limit, we recover the topological derivative expression.

In Figure 7, we sketch these ideas of taking limits respect to the radius variation (interpretation of the shape derivative) and respect to the radius itself (interpretation of the topological derivative).

Illustration of the topological derivative interpretation by passing the shape (variation of the radius of the inclusion τ) and the topology (radius of the inclusion ϵ) to the limit.
Figure 7: Illustration of the topological derivative interpretation by passing the shape (variation of the radius of the inclusion ) and the topology (radius of the inclusion ) to the limit.

In mathematical terms, the topological derivative is defined, more precisely, as the linear operator that fulfills the following expansion

(2.51)

where and represent the domain with an without an inclusion. The term represents a circular inclusion in the point with a radius of value . Concerning the shape function , most of the times, it is taken similarly to the SIMP and shape optimization methods, i.e., the volume of the domain, the compliance of the structure or the norm difference between the displacement of a specific point of the domain with a displacement target value. However, the compliance functional plays usually a crucial role in the topological optimization problems. For this reasons, obtaining topological derivative of the compliance was an important stride in the topic of topological derivative. Specifically, the compliance shape functional is usually defined as

(2.52)

where and represent the stresses and displacements solution of the standard elastic problem (2.46). Following Novotny in work [21], the topological derivative of the compliance is typically [21] presented in the following terms

(2.53)

where stands for the fourth order polarization tensor. The expression of such polarization tensor, for the isotropic and anisotropic case, and the procedure to relate it with the topological derivative is fully explained in chapter 4.

3 Topological derivative and topology optimization

3.1 Motivation

As mentioned in Chapter 2, different methodologies are available to address the topology optimization problem. Certainly, the different formulations impact meaningfully on the numerical strategy proposed to solve the problem. For instance, the regularization of the characteristic function introduced in the SIMP method allows computing a continuous gradient. However, when using topological derivative, no continuous gradient is available. This difference results according to two significant considerations: on the one hand, standard KKT conditions are no longer imposed as optimality conditions; on the other hand, the usual continuous optimization algorithms (steepest descent, Newton Raphson, ...) must be replaced by alternative algorithms.

Regarding the optimality conditions, it is convenient to provide first an intuitive description of the topological derivative concept and, then, formalize it in mathematical terms. From that descriptions, the optimality conditions arise naturally.

The (non standard) topological derivative algorithm is also convenient to be described. The level-set updating, built to satisfy the optimality conditions, is not straightforward. In addition, this topological derivative algorithm becomes more complex when imposing constraints in the minimization problem, specially with inequality constraints.

Furthermore, the topological derivative algorithm presents two drawbacks of different nature. On the one hand, determining the line search parameter is not an easy task; significant oscillations appear leading to spurious local minima. On the other hand, for a threshold of the stopping criteria, the algorithm may not converge leading to time-consuming re-meshing processes. To alleviate both inconveniences, novel numerical strategies must be proposed.

3.2 Optimality conditions when using topological derivative

Once the topology optimization problem (2.43) has been stated in Chapter 2, the following question arises: how the optimality condition in the topology optimization problem must be imposed when using the topological derivative? Note that, strictly speaking, the topological derivative does not correspond to a continuous gradient, and consequently, non standard KKT conditions can be imposed.

3.2.1 Qualitative description of inserting an inclusion

For this purpose, an qualitative description of inserting an inclusion is first provided.

Draft Samper 118254298-InclusionA.png Representation of the initial topology and the update topology for a strong-to-weak material modification (Case A) and for a weak-to-strong material modification (Case B).
Figure 8: Representation of the initial topology and the update topology for a strong-to-weak material modification (Case A) and for a weak-to-strong material modification (Case B).

Let's assume the domain is divided into two sub-domains: , endowed with a strong material (or material ) and , endowed with a weak material (or material –). In Figure 8, the strong material is represented in gray and the weak material in white. The main idea is based on studying how a shape functional (compliance, volume, ...) is modified when a circle of one material (strong or weak) is replaced by a circle of the other material (weak or strong). In an attempt to give a qualitative idea of this relation, the following two cases are introduced:

  • Case A: when a small circle of the strong material is replaced by a small circle of the weak material, see Figure StrongToWeak.
  • Case B: when a small circle of the weak material is replaced by a small circle of the strong material, see Figure 8a.

In this work, the shape functional is restricted to the compliance, denoted by , and to the volume (or mass in general), denoted by . In addition, we denote the shape functionals with and when the center of the inclusion is inserted in or respectively. Thus, the change of the shape functionals in both cases is as follows:

  • Case A: Since the small circle made of the strong material is replaced by the weak material, i.e., to , the structure must behave less stiffer and, consequently, the compliance should increase as follows

(3.1)

On the contrary, when is replaced by , the mass of the domain should decrease. This is,

(3.2)
  • Case B: Since the small circle made of the weak material is replaced by the strong material, i.e., to , the structure must behave stiffer and the compliance should decrease as follows

(3.3)

On the contrary, when is replaced by , the mass of the domain should increase. This is,

(3.4)

Thus, it can be inferred from this analysis the opposite response of the compliance and volume functionals when inserting an inclusion. In practice, this is translated to obtain solution different from the trivial one (full domain of strong material or full domain of weak material). In addition, this trend of the inequalities will incur on understanding the sign of the topological derivative.

3.2.2 Mathematical description of inserting an inclusion

Let's analyze, more formally, the mathematical formulation of inserting an inclusion. The circular inclusion is represented by means of the ball defined as

(3.5)

where and stands for the center and the radius of the ball. The characteristic function on the ball centered in with radius as

(3.6)

and the sign function as

(3.7)

Bearing this in mind, a general description of the constitutive tensor from the initial topology to the modified topology (by only inserting one hole) can be mathematically expressed as

(3.8)

Note that if (case A), the constitutive tensor changes from to in the circular inclusion and remains the same in the rest of the domain. An opposite behavior occurs in case B. Similarly, we express the change of the density from the initial topology to the modified one as

(3.9)

3.2.3 Topological derivative of the volume

Let's express first the mass of the domain as the integral of the density over all the domain, i.e

(3.10)

Note that, in fact, it represents the mass of the domain before inserting the inclusion. Equivalently, the mass of the domain after inserting the inclusion corresponds to

(3.11)

By definition, see [67], the topological derivative holds in the following expansion,

(3.12)

Hence, by identifying terms, we obtain straightforwardly the expression of the topological derivative for the mass (or volume) shape functional as

(3.13)

Note that, for cases A and B, the topological derivative takes the following values,

(3.14)

(3.15)

3.2.4 Topological derivative of the compliance

In contrast to the volume functional, the compliance clearly depends on the constitutive tensor and is quite commonly expressed as

(3.16)

where stands for the external boundary forces and represents the displacements, solution of the weak form of the equilibrium equation (2.46), i.e.,

(3.17)

Note that the body forces are neglected for the sake of simplicity. Expression (3.16) corresponds in fact to the compliance of the initial topology of the domain. The difference between the modified topology and the initial topology is typically written, see work [21] and [68], in terms of the polarization tensor , as follows

where are the stresses, solution of (2.46). In the context of linear elasticity, they are described as . Thus, by identifying terms, the topological derivative is reduced to

(3.18)

At this point, it is worth stressing that the polarization tensor has the following properties

(3.19)

(3.20)

The full expressions of the polarization tensor for both isotropic and anisotropic material are obtained in Chapter 4

As a summary of the variables and properties described so far, in Table 1, we show the initial and modified topology properties, the shape functional increments and the topological derivative for the Case A and B in a compact form.


Table 1. Interpretation of the topological derivative interpretations
Case A (Strong to weak) Case B (Weak to strong)

3.2.5 Optimality conditions

Once the topological derivatives are described for the volume and compliance and the topological optimization problem is stated, we present the optimality conditions of the topological optimization problem when using topological derivative. Let's assume that a shape functional (the compliance or the volume in our case) depends on a material parameter (the constitutive tensor or density) that takes values on and on . Defining an arbitrary direction as

(3.21)

where determine an specific hole and the number of holes, we say, by definition, that is a local minimizer if, for any direction (and for any number of holes ), the shape function will always increase, this is

(3.22)

If can be expanded asymptotically (the topological derivative exists), the difference of the shape functional is expressed, by definition, as

Thus, for small enough values , the necessary optimality conditions are set by imposing positivity on the topological derivative, i.e.,

(3.23)

Accordingly, the topological optimization algorithm must fulfill condition (3.23). A similar description of the optimality conditions can be found in work [67].

3.3 The Slerp algorithm

The lack of continuous gradient, and consequently, of standard continuous algorithms (steepest descent, Newton...) gave rise to propose, as an alternative, the Slerp algorithm for solving topology optimization problems when using the topological derivative. We recall that it was resourcefully achieved by Amstutz and co-workers in the seminal work [49].

3.3.1 The level-set function

The methodology lies essentially on defining the topology via a continuous function, usually called level-set function. More specifically, the zero level of that function determines the characteristic function and, consequently, the topology. As a main advantage, this method allows obtaining complex and very different topologies by means of a slight continuous change of the level set function, as we can observe in Figure 9.

Level-set function representation and its relation with the characteristic function (topology) [LevelSet2016]. Slight variations of the level-set function (assuming a small “off-set”) entail large differences on the topology.
Figure 9: Level-set function representation and its relation with the characteristic function (topology) [LevelSet2016]. Slight variations of the level-set function (assuming a small “off-set”) entail large differences on the topology.

More formally, the characteristic function is defined on the domain usually by the level-set function as

(3.24)

where represents the Heaviside function. Note that, owing to this definition, the design variable of the topology optimization problem switches, in practice, from to .

3.3.2 Optimality condition when using a level-set function

Let's now relate the level-set function with the optimality condition (3.23). The function is defined as the scalar function satisfying

(3.25)

Hereafter, we will use or indistinctly when referring to the topological derivative. Thus, with this definition in mind, the necessary optimality condition (3.23) becomes as follows

(3.26)

Consequently, both the level-set function (by definition) and the topological derivative function (in order to fulfill optimality conditions) satisfies

(3.27)

In view of this result, one can set the optimality condition for the level-set function as

(3.28)

Note that, since the topological derivative depends on the topology, and consequently on the level set, the above equation is highly non-linear.

3.3.3 Slerp algorithm for unconstrained optimization problems

From the optimality conditions of the level-set function, one can naturally impose parallelism between the level-set function and the topological derivative at the minimum, that is

(3.29)

where . Note that this relation fulfills automatically the optimality condition (3.28) for the level-set function. One could think on using equation (3.28) in a fix-point algorithm to get the optimality conditions. However, since is an arbitrary parameter, the level set function can increase unlimitedly and, consequently, the algorithm may not converge. In order to mitigate such inconvenient, one could fix such that the level-set function is enforced to have unitary norm. Thus, by taking , equation (3.29) becomes

(3.30)

This relation satisfies equation (3.28) and can be understood as a fix-point algorithm: given a topology through a level-set function, compute its topological derivative and, by normalizing it, obtain the new level-set function and, consequently, the new topology. However, this fix point method could be highly aggressive leading to no descent direction during the iterations.

As a remedy, instead of updating the level-set in terms only of the new normalized topological derivative, one could combine it with the previous value of the level set, in such a way that the updated level-set function has unitary norm. This strategy is, in fact, what the slerp algorithm proposes.

The name of slerp arises from the computer graphics community and is shorthand for spherical linear interpolation. Shoemake [69] introduced this concept in the quaternion interpolation context for the propose of animating 3D rotation. In the topological optimization context, the slerp algorithm was proposed by Amstutz in work [49]. Note that in the context of quaternion interpolations, the objects to be interpolated are vectors of dimension four while in the context of topology optimization, the objects to be interpolated are fields defined over the domain.

In general, the slerp algorithm can be understood as the interpolation of two different functions on the unit sphere. In Figure (10), we show schematically the relation between the new level set function , the actual level-set function and the topological derivative .

Draft Samper 118254298-Slerp.png Representation of the slerp algorithm. The updated ψₙ₊₁ level-set function is computed by interpolating the actual level-set function ψₙ and the topological derivative gₙ on the unit sphere. The topological derivative gₙ plays the role of a descent direction on a steepest descent algorithm and the variable κₙ plays the role of a line-search parameter.
Figure 10: Representation of the slerp algorithm. The updated level-set function is computed by interpolating the actual level-set function and the topological derivative on the unit sphere. The topological derivative plays the role of a descent direction on a steepest descent algorithm and the variable plays the role of a line-search parameter.

According to the triangle relation described in Figure (10), the following equation must hold

(3.31)

where the scalar numbers and can be computed by imposing the following law of sinus

(3.32)

Note that this last relation is fulfilled due to the unitary norm of , and . In Figure (11), we depicted that trigonometric relation.

On the left, the vector relation between the updated ψₙ₊₁ level-set function, the actual level-set function ψₙ and topological derivative gₙ . On the right, the triangular relation that leads to find the scalar values αₙ and βₙ.
Figure 11: On the left, the vector relation between the updated level-set function, the actual level-set function and topological derivative . On the right, the triangular relation that leads to find the scalar values and .

Thus, the new level set function can be written as a combination of the actual level-set function and the topological derivative as follows

(3.33)

where is a line search-like parameter and the angle between and which is written as

(3.34)

Note that, an alternative way of imposing parallelism between the level-set function and the topological derivative s achieved by requires zero vale of the angle . In this respect, the stopping criteria of the algorithm can be set by imposing the tolerance as a threshold of the angle .

In addition, it is remarkable that, by using the slerp algorithm, the updated level-set function has automatically unit norm. Note that since , we have chosen the norm for both, the norm and the scalar product of equations (3.33) and (3.34). However, other norms (or ) can be used.

Moreover, the line-search parameter allows controlling the size of the step. Initially, it is usually set to and it is divided by 2 until the new topology provides a smaller cost function. Clearly, if takes unitary values, we recover the aggressive fix-point algorithm proposed in equation (3.30).

Thus, the slerp algorithm is a fix-point algorithm with a line-search parameter that interpolates in the unit sphere the actual level-set with the topological derivative. It is worth mentioning that the slerp algorithm can be seen in the optimization context as a standard steepest descent method with the particularity that the update variable must have unit norm. Consequently, the topological derivative plays the role of the gradient in the topology optimization problem, hence its importance. Although, in this work, an exhaustive numerical analysis has not been considered (see [70] for this purpose), one can expect linear convergence of the algorithm.

3.3.4 Slerp algorithm combined with an augmented Lagrangian scheme for constrained optimization problems

So far, we have described the slerp algorithm for unconstrained topology optimization problem. However, most of the applications require fulfilling some constraints. In this sub-section we present how to deal with the case of a volume constraint. To the author's knowledge, there is no a vast amount of algorithms to deal with constrained topology optimization problems when using topological derivative. The lack of continuous gradient undermine the possibility of using standard continuous optimization algorithms. However, the augmented Lagrangian algorithm is exempt of such limitation since it uncouples the update of the design variables (the topology in our case) and the update of the Lagrange multipliers [71]. Thus, in view of this property, the slerp algorithm can be combined with an augmented Lagrangian scheme in constrained topology optimization problems. The reader is referred to works [72], [73] and [74] for further information.

Equality constraints

Considering the level-set function as the design variable, the topology optimization problem with volume constraint (2.43) may be expressed as

(3.35)

According to the reference book [71], and following the works [75,72,74], the augmented Lagrangian scheme for equality constraints proposes to solve the minimization problem (3.35) by introducing the following saddle point problem

(3.36)

where is the Lagrangian functional, is the Lagrange multiplier and the penalty parameter. The augmented Lagrangian scheme is based upon solving the primal-dual problem sequentially with the particularity that an extra term is added in order to convexify the problem [71]. The algorithm considers first, a single (or multiple) iteration for minimizing the cost function and then a single iteration as

(3.37)

for maximizing the cost function, i.e., an Usawa-like scheme [76]. Certainly, the augmented Lagrangian scheme has an impact in the slerp algorithm. On the one hand, the cost function is replaced by the Lagrangian functional when determining the line-search parameter . On the other hand, the topological derivative of the Lagrangian functional must be computed by considering an extended topological derivative as

(3.38)

where stands for the topological derivative of the cost function ). Note that, the topological derivative of the volume constraint is equal to and has been consequently omitted. Hereafter, for simplicity, the extended topological derivative is also called topological derivative.

Inequality constraints The augmented Lagrangian scheme can be extended to minimization problems with inequality constraints. The main idea consists in retrieving the minimization problem with equality constraints by adding an extra variable , often termed slack variable, to the minimization problem with inequality constraints. In mathematical terms, it reads as

(3.39)

where the constraint function has been defined. Then, following the augmented Lagrangian scheme for equality constraints, the augmented Lagrangian is defined and the following saddle-point problem

(3.40)

must be solved. The procedure consists in imposing one of the KKT conditions to isolate the slack variable in terms of the design variable and the Lagrange multiplier , and then inserting its expression in the augmented Lagrangian . For this purpose, we enforce that the partial derivative of the Lagrangian functional with respect to the slack variable must be canceled. This is,

(3.41)

The two possible solution of the above equation are written as

(3.42)

In order to determine the solution that leads to a smaller augmented Lagrangian value, we first insert and into the definition of the constraint function as follows

(3.43)

and then to the Lagrangian functional as

(3.44)

(3.45)

Comparing both Lagrangian functionals

(3.46)

we obtain that solution always provide a smaller value of the Lagrangian functional. Thus, the optimal slack variable takes always value, provided of course that exists, this is when . Otherwise, the optimal slack variable takes the other solution . Therefore, it is useful to write the square of the optimal slack variable as

(3.47)

or, more compactly, as

(3.48)

By inserting this last result in the constraint function , we obtain

(3.49)

where we have conveniently defined the constraint. Note that the problem no longer depends on the slack variable appears. Thus, the saddle-point problem for solving minimization problem with inequality constraints (3.40) becomes

(3.50)

which represents a standard saddle-point problem for solving minimization problem with equality constraints (see equation (3.36)) with the particularity that the constraint depends explicitly on the Lagrange multiplier . At this point, it is convenient to examine if problem (3.50) can be treated as a standard minimization problem with equality constraints. For this purpose, the KKT

(3.51)

are first imposed with slight abuse of notation. The is not strictly imposed on the topology optimization problem when using topological derivative. Instead, the optimality conditions (3.23) are considered. In any case, from equation (3.49) the and terms are determined as follows

(3.52)

In addition, it can be shown that

(3.53)

and, consequently, the product between these last two equations leads to the following result

(3.54)

Thus, the KKT condition of the augmented Lagrangian with inequality constraint becomes

(3.55)

which coincide with the KKT condition of the augmented Lagrangian with equality constraint when the the constraint is replaced by the constraint . This is,

(3.56)

This last consideration extends the augmented Lagrangian scheme with equality constraint neatly to the inequality-constrained case. Thus, the Lagrange multiplier is updated in the same manner, i.e.,

(3.57)

or, more explicitly,

It is worth stressing that this last equation shows how the augmented Lagrangian cancels the Lagrange multiplier when the inequality is not active. Likewise, the extended topological derivative is defined in the same manner of equation (3.38), i.e.,

(3.58)

As an advantage, this scheme allows using the same implementation for both cases; equalities or inequalities; by means of exchanging the constraint with the constraint . For further information, the reader is referred to the reference book [71] for a rigorous description and work [72] for applying to topological optimization problem.

Penalty value and algorithm. It is worth mentioning that giving an adequate value to the penalty variable is not an easy task. In the reference book [71], it is suggested increasing the penalty during the iterations (when the constraint is not tightened enough). Similarly, in work [77], it is proposed to increase the penalty every five iterations. However, our experience shows us that the penalty cannot increase unlimitedly. A very small value of the penalty will make the problem converge very slowly but if, on the contrary, is very high the solution oscillates with the risk of non-convergence. So, as it was mentioned before, the appropiate value of this parameter will depend on the problem and on the sensibility of the cost when varying the constraint. The numerical experiments, further presented, suggest to normalize the cost function and take a small value of the penalty.

A detailed scheme of the Augmented Lagrangian Slerp algorithm is presented in Algorithm 1.

Init: choose initial values of , , tol, , and
Compute and from (2.46).
Compute from (3.56) with from (3.10).
Compute from (3.58) with computed from (3.18).
while do
Set , , , .
while (line search) do
Update from (3.33) with , and , and update from (3.24).
Compute and from (2.46).
Compute from (3.50), with from (3.16) and set and .
Set , , ,
Compute from (3.34)
Compute from (3.56) with from (3.10).
Update from (3.57).
Compute from (3.58) with computed from (3.18) and set .
Algorithm. 1 Augmented Lagrangian slerp algorithm. The Lagrange multiplier is updated in every topology iteration, i.e., Usawa-like scheme is used.

3.4 Treatment of the interface

The stresses and the topological derivative deserve special attention when dealing with the interface elements. In this section, we analyze an specific treatment of the interface elements in a bi-material elastic problem. Note that, although in topology optimization the aim is where to nucleate holes and where to leave material, in practice, the aim is how to distribute two materials, provided that the weak material takes significant lower stiff (usually times less).

3.4.1 Bi-material elastic problem

We proceed to examine how to deal with the bi-material problem, from the continuous and discrete point of view, and how the different treatments of the interface affect the cost and the topological derivative.

Formulation

In the bi-material problem, the domain is first divided in two parts, one sub-domain with a strong material and the other sub-domain with a weak material . Then, the classical form of the elastic problem (without considering, for simplicity, the body forces) is written as

(3.59)

and by multiplying times the test function and integrating by parts, we obtain the weak form of the bi-material elastic problem as

(3.60)

where the domain has been split into t and . Let's consider a a Finite Element discretization of equation (3.60). The domain is equipped with a conforming, triangular mesh composed of triangles , , and vertices , . Recalling the classical Finite element spaces, is the finite-dimensional space of Lagrange Finite Element functions, i.e. of affine functions in restriction to each triangle . A basis of is composed of the functions , , where is the unique element in such that if and otherwise. We take advantage of defining as the finite-dimensional space of Lagrange Finite Element functions on , i.e. of constant functions in restriction to each triangle . A basis of is composed of the functions , , where on and on , .

Thus, defining as the diameter of the triangles , the discretized displacement is the unique solution of the discretized weak form

(3.61)

Let's split the domain into three different sub-groups: the elements (endowed with constitutive tensor), the elements (endowed with constitutive tensor), and the elements intersecting the interface that share both materials and are, for the time, denoted by the interface constitutive tensor . The different possible definitions of that interface constitutive tensor are further described and represents the main ingredient of the bi-material elastic problems. Bearing this in mind, we can define the regularized constitutive tensor as

(3.62)

Consequently, the weak form equation (3.60), is re-written, after discretization, as

(3.63)

The right side of the above equation has not been split into the elements because no it does not depend on the constitutive tensor.

Undoubtedly, the main difficulty of equation (3.63) relies on how to deal with the interface elements . Since we use Finite Element, the strains are element-wise constant and, consequently, the integration of the last term of the left hand side of 3.63 requires the value of on the Gauss points of the interface elements , in terms of material and . From definition (2.48), the regularized constitutive tensor on the interface is expressed in terms of the characteristic function as

(3.64)

In view of equation (3.64), the description of the interface is based on the treatment of the characteristic function in the Gauss Points of the interface elements and, in turn, on the treatment of the level-set function . The characteristic function evaluated on the Gauss Points of the interface elements is, hereafter, termed as . For this purpose, two different approaches are commonly used in topology optimization when using topological derivative [49,78]: the In or Out approach and the P1-projection approach.

In or Out approach

The In or Out approach is based on taking the characteristic function on the interface in terms of the level-set function evaluated on the Gauss point . This is,

(3.65)

Note that, in this case, the characteristic function on the Gauss point takes binary values, i.e.,

(3.66)

Consequently, the constitutive tensor is restricted to and , i.e., it belongs to

(3.67)

and, in terms of the level-set function, is written as

(3.68)

P1-projection approach

Alternatively, a -projection approach can be considered. It consists on projecting the characteristic function from to the smaller finite-dimensional space , composed by the Finite Element functions . In mathematical terms, the -projection reads as

(3.69)

where the values on the nodes are determined from the nodal level-set values , i.e.,

(3.70)

Thus, since we are dealing with Finite Element functions, the value of the characteristic function on the Gauss point is reduced to

(3.71)

where , and are the values of the characteristic function on the nodes and of the element . In the following, depending on the values of the level-set function, we examine the four different cases that may appear. Without loose of generality, we order the nodal values of the level-set function as follows: where , and are the position of the the nodes and .

Case A: All the values of the level-set function are negative and consequently all the values of the characteristic function are equal to , including the value on the Gauss point. This is,

(3.72)

In fact, in this case, the element is not considered an interface element.

Case B: