TOPOLOGICAL DERIVATIVE-BASED TOPOLOGY OPTIMIZATION OF STRUCTURES SUBJECT TO MULTIPLE LOAD-CASES

The topological derivative measures the sensitivity of a shape functional with respect to an infinitesimal singular domain perturbation, such as the insertion of holes, inclusions or source-terms. The topological derivative has been successfully applied in obtaining the optimal topology for a large class of physics and engineering problems. In this paper the topological derivative is applied in the context of topology optimization of structures subject to multiple load-cases. In particular, the structural compliance under plane stress or plane strain assumptions is minimized under volume constraint. For the sake of completeness, the topological asymptotic analysis of the total potential energy with respect to the nucleation of a small circular inclusion is developed in all details. Since we are dealing with multiple load-cases, a multi-objective optimization problem is proposed and the topological sensitivity is obtained as a sum of the topological derivatives associated with each load-case. The volume constraint is imposed through the Augmented Lagrangian Method. The obtained result is used to devise a topology optimization algorithm based on the topological derivative together with a levelset domain representation method. Finally, several finite element-based examples of structural optimization are presented.

In this paper the topological derivative is applied in the context of topology optimization of structures into two spatial dimensions subject to multiple load-cases. In particular, the structural compliance under plane stress or plane strain assumptions is minimized subject to volume constraint, which is imposed through the Augmented Lagrangian Method. For the sake of completeness, the topological asymptotic analysis of the total potential energy with respect to the nucleation of a small circular inclusion is developed in all details. These derivations can be found in the literature (see for instance Novotny and Sokołowski (2013) and references therein). However, our idea is to present them in a simplified and pedagogical manner by using simple arguments from the Analysis. In addition, we claim that the topological derivative obeys the basic rules of the Differential Calculus (see the examples on the introduction of the book by Novotny and Sokołowski (2013) and also the paper by  for application of these rules). Thus, since we are dealing with multiple load-cases, a multi-objective optimization problem is proposed and the topological sensitivity is obtained as a sum of the topological derivatives associated with each load-case. It is also worth to mention that the topological derivative is defined through a limit passage when the small parameter governing the size of the topological perturbation goes to zero. However, it can be used as a steepest-descent direction in an optimization process like in any method based on the gradient of the cost functional. Therefore, the obtained result is used to devise a topology optimization algorithm based on the topological derivative together with a level-set domain representation method, as proposed by Amstutz and Andrä (2006). The algorithm is presented in a pseudo-code format easy to implement. Finally, several finite element-based examples of structural optimization are presented. In summary, a comprehensive account on the application of the topological derivative in the context of compliance structural optimization is given. The theoretical development, interpretation of the results and computational aspects are discussed, and some misunderstandings currently found in the literature are elucidated. This paper is organized as follows. In Section 2 we firstly introduce the topological derivative concept and state the mechanical problem which we are dealing with, then two main results are derived: the existence and the closed formula for the topological derivative associated to the strain energy shape functional. In Section 3 the compliance topology optimization problem under volume constraint is stated together with its associated topological derivative. The topology optimization algorithm based on the topological derivative and a level-set domain representation method is described in Section 4. The numerical results are presented in Section 5. The paper ends with some concluding remarks in Section 6. Finally, the closed formula for the topological derivative presented in Section 2 is derived in Appendix A by using standard arguments from the Analysis.

TOPOLOGICAL ASYMPTOTIC ANALYSIS
Let us consider an open and bounded domain denoted by Ω ⊂ R 2 . Associated with this domain, we introduce a characteristic function χ : where |Ω| is the Lebesgue measure of Ω. Suppose that Ω is subject to a singular perturbation confined in a small region ω ε ( x) = x + εω with size ε, where x is an arbitrary point in Ω and ω is a fixed domain in R 2 . We define a characteristic function χ ε ( x; x), x ∈ R 2 , associated with the topologically perturbed domain. In the case of a perforation, for instance, χ ε : R 2 → {0, 1}, χ ε ( x; x) = 1 Ω − 1 ωε and the perforated domain is obtained as Ω ε = Ω \ ω ε . Then we assume that a given shape functional ψ(χ ε ( x)) associated with the topologically perturbed domain, admits the following topological asymptotic expansion where ψ(χ) is the shape functional associated to the unperturbed domain Ω and f (ε) is a positive function such that, f (ε) → 0 when ε → 0. The function x → T ( x) is called the topological derivative of ψ at x. Therefore this derivative can be seen as a first order correction on ψ(χ) to approximate ψ(χ ε ). In addition, after rewriting (2.2) we obtain the classical definition for the topological derivative (Sokołowski and Żochowski (1999)): Note that the topological derivative is defined through the limit passage ε → 0. However, according to (2.2), it can be used as a steepest-descent direction in an optimization process similar to any gradient-based method, as shall be presented later.
In this paper the domain is topologically perturbed by the nucleation of a small inclusion, as shown in Fig. 1. More precisely, the perturbed domain is obtained when a circular hole ω ε ( x) := B ε ( x), the ball of radius ε > 0 and center at x ∈ Ω, is introduced in Ω. Next, this region is filled by an inclusion with different material property from the background. In particular, χ ε ( x) = 1 Ω − (1 − γ)1 Bε( x) and a piecewise constant function γ ε is introduced: where γ ∈ R + is the contrast in the material property.

Unperturbed problem
As mentioned before, the topological asymptotic expansion of the total potential energy associated with the elasticity system into two spatial dimensions is obtained. Thus, the unperturbed shape functional is defined as: where the vector function u is the solution to the variational problem: with b a constant body force distributed in the domain, the Cauchy stress tensor, the linearized strain tensor and C the constitutive tensor given by where I and I are the second and fourth identity tensors, respectively, µ and λ are the Lamé's coefficients, both considered constants everywhere. In the plane stress assumption we have (2.10) while in plane strain assumption they are where E is the Young modulus and ν the Poisson ratio. The set of admissible displacements U and the space of admissible displacements variations V are respectively defined as Here, Γ D and Γ N respectively are Dirichlet and Neumann boundaries such that ∂Ω = Γ D ∪ Γ N with Γ D ∩ Γ N = ∅, u is the displacement prescribed on Γ D and q is the load prescribed on Γ N , both assumed to be smooth enough. See details in Fig. 2. The strong system associated with the variational problem (2.6) is given by: (2.13)

Perturbed problem
Now let us state the associate topologically perturbed problem. In this case, the total potential energy is given by where the vector function u ε is the solution to the variational problem: where γ ε is defined by (2.4). The set U ε and the space V ε are defined as  to the variational problem (2.15) reads: (2.18)

Existence of the topological derivative
The following result ensures the existence of the topological derivative associated with the problem under analysis.
Lemma 1. Let u and u ε be solutions to (2.6) and (2.15), respectively. Then we have that the estimate Proof. We start by subtracting the variational problem (2.6) from (2.15) to obtain: From the above equation, we have: (2.20) Recalling that: (2.21) By taking η = u ε − u as test function in the above equation we obtain the following equality From the Cauchy-Schwartz and Poincaré inequalities it follows that where we have used the elliptic regularity of function u and the continuity of the function b at the point x ∈ Ω. Finally, from the coercivity of the bilinear form on the left-hand side of (2.22), namely we obtain the result with the constant C = C 6 /c independent of the small parameter ε.

2.4
The topological derivative formula According to Novotny and Sokołowski (2013) the topological asymptotic expansion of the energy shape functional takes the form (see also Appendix 6): where the polarization tensor P γ is given by the following fourth order isotropic tensor Finally, in order to extract the main term of the above expansion, we choose f (ε) = πε 2 , which leads to the final formula for the topological derivative, namely: (2.28) Remark 2. Formally we can take the limit cases γ → 0 and γ → ∞. For γ → 0, the inclusion leads to a void and the transmission condition on the interface of the inclusion degenerates to homogeneous Neumann boundary condition. In this case the polarization tensor is given by In addition, for γ → ∞, the elastic inclusion leads to a rigid one and the polarization tensor is stated as In the case of plane strain linear elasticity, the above formulas are written as while in the case of plane stress linear elasticity the formulas are explicitly given by Remark 3. The polarization tensor P γ and theirs associated particular representations lead to isotropic fourth order tensors because we are dealing with circular inclusions as topological perturbations. For arbitrary shaped inclusions the reader may refer to the book by Ammari and Kang (2007), for instance. On the other hand, there are two main advantages in using circular inclusions in the context of topology optimization, which are: • The associated topological derivative is given by a closed formula depending on the solution to the original unperturbed problem.
• There are optimality conditions rigorously derived in Amstutz (2011), allowing for use the topological derivative together with a level-set domain representation method as a steepestdescent direction in a topology optimization algorithm (Amstutz and Andrä (2006)).

THE TOPOLOGY OPTIMIZATION PROBLEM
We consider the compliance topology optimization of structures subject to multiple load-cases under volume constraint. Therefore, the topology optimization problem can be stated as P1 : where NLC is the number of load-cases, M > 0 is the required volume at the end of the optimization process and u i is solution to: Here, D is used to denote a hold-all domain. In order to simplify the numerical implementation we consider that the elastic body D is decomposed into two sub-domains Ω and ω. The domain Ω = D \ ω represents the elastic part while ω ⊂ D is filled with a very complacent material, used to mimic voids. See the sketch in Fig. 4. This procedure allows us to work in a fixed computational domain. Therefore, we define a characteristic function of the form In addition, we introduce a piecewise constant function ρ, such that with 0 < ρ 0 1 used to mimic voids. That is, the original optimization problem, where the structure itself consists of the domain Ω of given elastic properties and the remaining empty part ω, is approximated by means of the two-phase material distribution given by (3.5) over D where the empty region ω is filled by a material (the soft phase) with Young's modulus, ρ 0 E, much lower than the given Young's modulus, E, of the structure material (the hard phase).

Augmented lagrangian
The volume constraint is imposed through the Augmented Lagrangian Method. It consists in transform the inequality constraint in Problem P1 given by (3.1) into an equality by introducing a slack function s, namely P2 : where g Ω = (|Ω| − M )/M . Note that for an optimal value s = s * the equivalence P2 ≡ P1 holds true. Let us introduce two Lagrange multipliers α and β. Then, the Problem P2 in (3.6) can be rewritten as P3 : subject to h Ω (s) = g Ω + s 2 . (3.7) Note that if h Ω (s * ) = 0, then P3 ≡ P2 ≡ P1. Thus, let us minimize (3.7) with respect to the slack function s, namely αh Ω (s) + βh Ω (s)h Ω (s) = 0, h Ω (s) = 2s. (3.8) Then, s = 0 or s 2 = −(α/β + g Ω ). Therefore It follows that: h Ω (s * ) = g Ω + max{0, −(α/β + g Ω )} = max{g Ω , −α/β} := g + Ω . (3.10) Finally, by setting s = s * we have P4 ≡ P3 ≡ P2 ≡ P1, that allows us to rewrite the constrained Problem P1, given by (3.1), in its equivalent unconstrained form, namely (3.11) In addition we have that β, the parameter associated to the quadratic term, controls the updating of the parameter α, which is associated with the linear term. Therefore, we can specify the required volume fraction at the final of the optimization process by solving the following recursive formula for the parameter α: where β > 0 is a fixed number. This updating process shall be repeated until the volume constraint is reached, as can be seen in Algorithm 1 at the end of Section 4.

Topological derivative
Let us consider two extremal situations. When ρ = 1, then γ = ρ 0 . On the other hand, if ρ = ρ 0 , then γ = 1/ρ 0 . Therefore, in order to simplify the numerical implementation, we can eliminate the contrast parameter γ by considering the limit cases discussed in Remark 2. In addition, since the topological derivative satisfies the basic rules of the Differential Calculus (Novotny and Sokołowski (2013)) and taking into account the linearity of the elasticity problem, we have (3.13) Therefore, the topological derivative associated with Problem P4 is obtained as a sum of the topological derivatives for each load-case together with the topological derivative of the augmented Lagrangian terms, which is trivially obtained by considering the limit cases in Remark 2.

THE TOPOLOGY OPTIMIZATION ALGORITHM
In this section a topology optimization algorithm based on the topological derivative together with a level-set domain representation method is presented. It has been proposed by Amstutz and Andrä (2006) and consists basically in looking for a local optimality condition for the minimization problem (3.11), written in terms of the topological derivative and a level-set function. Therefore, the elastic part Ω as well as the complacent material ω are characterized by a level-set function Ψ ∈ L 2 (D) of the form: where Ψ vanishes on the interface ∂ω. A local sufficient optimality condition for Problem (3.11), under the considered class of domain perturbation given by circular inclusions, can be stated as (Amstutz (2011)) T (x) > 0 ∀x ∈ D.
(4.2) Therefore, let us define the quantity allowing for rewrite the condition (4.2) in the following equivalent form We observe that (4.4) is satisfied wether the quantity g coincides with the level-set function Ψ up to a strictly positive number, namely ∃ τ > 0 : g = τ Ψ, or equivalently θ := arccos g, Ψ L 2 (D) which shall be used as optimality condition in the topology design algorithm, where θ is the angle between the functions g and Ψ in L 2 (D).
Algorithm 1: The topology design algorithm input : NLC, D, Ψ 0 , M , α 0 , β, κ , θ , M output: The optimal topology Ω 1 n ← 0; 2 Ω n ← Ψ n ; E pt ← 0; 10 end for 11 Compute F Ωn according to (3.11); 12 Compute T ( x) using (3.13); 13 Compute g n according to (4.3); 14 θ n ← arccos gn,Ψn gn L 2 (D) Ψn L 2 (D) Let us now explain the algorithm. We start by choosing an initial level-set function Ψ 0 ∈ L 2 (D). In a generic iteration n, we compute function g n associated with the level-set function Ψ n ∈ L 2 (D). Thus, the new level-set function Ψ n+1 is updated according to the following linear combination between the functions g n and Ψ n Ψ 0 ∈ L 2 (D), where θ n is the angle between g n and Ψ n , and κ is a step size determined by a linear-search performed in order to decrease the value of the objective function F Ωn , with Ω n used to denote the elastic part associated to Ψ n . The process ends when the condition θ n ≤ θ and at the same time the required volume |1 − |Ω n |/M | ≤ M are satisfied in some iteration, where θ and M are given small numerical tolerances. In particular, we can choose and by construction Ψ n+1 ∈ S, ∀n ∈ N. If at some iteration n the linear-search step size κ is found to be smaller then a given numerical tolerance κ > 0 and the optimality condition is not satisfied, namely θ n > θ , then a uniform mesh refinement of the hold all domain D is carried out and the iterative process is continued. The resulting topology design algorithm is summarized in a pseudo-code format, see Algorithm 1.

NUMERICAL EXAMPLES
Since we are dealing with multiple load-cases, two situations denoted by C1 and C2 are considered. In the first case (C1), the loads are applied simultaneously (single load-case) and the associated topological derivative is evaluated. On the other hand, in the second case (C2), the loads are applied separately (multiple load-cases) and the resulting topological derivative is obtained as a sum of the topological derivatives associated with each load-case. In all numerical examples the stopping criterion, the optimality threshold and the tolerance of the volume requirement are given respectively by κ = 10 −3 , θ = 1 o and M = 1%. The angle θ has converged to a value smaller than 1 o , namely, the optimality condition has been satisfied in all cases up to a small numerical tolerance. Furthermore, the mechanical problem is discretized into linear triangular finite elements and three steps of uniform mesh refinement were performed during the iterative process.
We assume that in the first three examples the structures are under plane strain assumption while in the last two examples the structures are under plane stress assumption. Finally, the material property threshold is set as ρ 0 = 10 −4 , while the Young's modulus is given by E = 1.0.

Example 1
Let us consider a prismatic bar submitted to a pair of loads, as shown in Fig. 5. The hold-all domain is given by a square section of size 1 × 1 supported on the two opposites bottom corners. The loading consists of a pair of forces q 1 = (−1.0, 0.0) and q 2 = (1.0, 0.0) applied on the two opposites top corners. The hold-all domain is discretized into a uniform mesh with 1600 elements and 841 nodes.
The Poisson ratio is given by ν = 0.3. The required volume fraction is set as M = 50%, while the parameters of the augmented Lagrangian method are given by α 0 = 0.0 and β = 10.0. The final topologies are obtained after 41 and 31 iterations for the cases C1 and C2, respectively, as shown in Figs. 6(a) and 6(b). We observe that the obtained result for multiple load-cases is feasible while the result for single load-case doesn't make sense from physical point of view.

Example 2
In this case the hold-all domain consists in a rectangular beam of size 2 × 1 supported on the two opposites bottom corners and subject to a pair of forces q 1 = (−0.5, −1.0) and q 2 = (0.5, −1.0), as shown in Fig. 7. The initial mesh is uniform with 1600 elements and 861 nodes.  The convergence of the θ n angle and volume fraction, for the cases C1 and C2, are respectively presented in Figs. 9 and 10. The jumps on the values of θ n in the graph of Fig. 9 are due to the mesh refinement. The multiple load-cases C2 allow for a more realistic combination of the applied forces resulting in topologies that best satisfy the conditions of the problem. Thus in the next three examples only the results for multiple load-cases are presented.

Example 3
This example simulates the design of a barrage. The hold-all domain is a square of size 1 × 1 clamped on its bottom edge, as illustrated by Fig. 11(a). The hydrostatic pressure is represented by a distributed load on the left-hand side of the barrage. There is also a constant body force b = (0.0, −0.1) acting in the whole domain. The initial mesh is uniform with 1600 elements and 841 nodes. The Poisson ratio is given by ν = 0.2. The required volume fraction is set as M = 50%, while the parameters of the augmented Lagrangian method are given by α 0 = 0.0 and β = 1.0. The final topology is obtained with 29 iterations, as shown in Fig. 11(b).

Example 4
Let us consider now a bridge design. The hold-all domain is given by a rectangle of size 3 × 1 supported on the two opposites bottom corners. The bridge is submitted to three uniformly distributed traffic loading q 1 = (0.0, −10.0), q 2 = (1.0, 0.0) and q 3 = (−1.0, 0.0) applied on the dark strip of height h = 0.05 positioned at the distance c = 0.45 from the top of the hold-all domain. This strip represents the road, which is simply supported on their opposites bottom corners, and therefore remains unchanged throughout the optimization process. We also consider a body force given by b = (0.0, −0.4). See the sketch in Fig. 12. The domain is discretized into a uniform mesh with 4800 elements and 2481 nodes.

Example 5
In this last example the design of an alloy wheel is considered. The hold-all domain is given by a ring of radii 0.2 and 1.0. The dark strip remains unchanged during the optimization process. The wheel is clamped on the smaller holes (little circles of radius 0.04). A uniformly distributed shear load q 9 of intensity 1.0 and eight normal loads q i , i = 1, ..., 8, of intensity 10.0 are applied on the contour of the wheel. All the details are presented in Fig. 14(a). The initial mesh is nonuniform with 7596 elements and 3937 nodes. The Poisson ratio is set as ν = 0.3, while the required volume fraction is M = 73%, where the parameters of the augmented Lagrangian method are given by α 0 = 0.0 and β = 12.0. The final topology is obtained with just 20 iterations and can be seen in Fig. 14(b).

CONCLUSION
In this paper the topological derivative was applied in the context of topology optimization of structures subject to multiple load-cases. The structural compliance into two spatial dimensions was minimized under volume constraint. Since the topological derivative obeys the basic rules of the Differential Calculus, the topological sensitivity of a multi-objective shape functional was obtained as a sum of the topological derivatives associated with each load-case. In addition, the obtained sensitivity has been used as a steepest-descent direction similar to any gradient-based method. In particular, a fixed-point algorithm based on the topological derivative together with a level-set domain representation method has been presented, which converges to a local optimum. The resulting algorithm has been summarized in a pseudo-code format easy to implement and several finite element-based examples of structural optimization were presented. Finally, for the reader convenience, the topological asymptotic analysis of the total potential energy with respect to the nucleation of a small circular inclusion was developed in a simplified and pedagogical manner by using standard arguments from the Analysis. Therefore, we believe that this paper would be useful for the readers interested on the mathematical aspects of topological asymptotic analysis as well as on applications of topological derivatives in structural optimization.

APPENDIX A ASYMPTOTIC ANALYSIS
In this appendix the proof of the result (2.25) together with the estimation for the remainder on the topological asymptotic expansion are presented. The results are derived by using simple arguments from the Analysis.

A.1 Polarization tensor in elasticity
In order to calculate the difference between the functionals ψ(χ) and ψ(χ ε ), respectively defined through (2.5) and (2.14), we start by taking η = u ε − u as test function in the variational problems (2.6) and (2.15), leading respectively to So that the shape functionals ψ(χ) and ψ(χ ε ) can be rewritten as: Now, after subtracting (.3) from (.4) we obtain: Note that the resulting terms of the difference between ψ(χ) and ψ(χ ε ) are given by integrals concentrated over the inclusion B ε . Now, in order to evaluate the limit in (2.3), we need to know the behavior of the function u ε with respect to ε → 0. Thus, let us introduce the following ansätz: where w ε (x) is the solution to an auxiliary exterior problem andũ ε (x) is the remainder. After applying the operator σ ε in (.6), we obtain where σ(u(x)) has been expanded in Taylor series around the point x and ξ is used to denote a point between x and x. From the transmission condition on the interface ∂B ε , we have Therefore, according to Fig. 3, n = (x − x)/ε, and Thus, we can choose σ ε (w ε ) such that: and by a changing variables, we write w ε (x) = εw(x/ε) and y = x/ε, ( which implies ∇ y w(y) = ε∇w(x/ε). In the new variable the following exterior problem is considered: with u = −(1 − γ)σ(u( x))n. The above boundary value problem admits an explicit solution. In fact, since the tress σ y (w) is uniform inside the inclusion, it can be written in a compact form making use of the Eshelby's Theorem (Eshelby (1957(Eshelby ( , 1959): where T is a fourth order isotropic tensor written as with the constants α 1 and α 2 given by (2.27). Now we can construct σ ε ( u ε ) in such a way that it compensates the discrepancies introduced by the higher-order terms in ε as well as by the boundary-layer σ y (w) on the exterior boundary ∂Ω. It means that the remainder u ε must be solution to the following boundary value problem: with h = −(1 − γ)(∇σ(u(ξ))n)n. Moreover, we can obtain an estimate for the remainder u ε of the form O(ε). In fact, before proceeding, let us state the following result, which can be found in the book Novotny and Sokołowski (2013) in its optimal version, namely O(ε 2 ): Lemma 4. Let u ε be solution to (.15). Then, the following estimate holds true: with the constant C independent of the small parameter ε.
Proof. From the expansion for u ε and making use of the triangular inequality, we can write where we have used the change of variables (.11), the equivalence between the semi-norm |·| H 1 (Ω) and the norm · H 1 (Ω) and the estimate in Lemma 1. Finally, the result comes out from the Poincaré inequality.

(.21)
Note that both tensors P γ and T are isotropic because we are dealing with circular inclusions. For arbitrary shaped inclusions the reader may refer to Ammari and Kang (2007).

A.2 Estimation of the remainders
The first remainder term E 1 (ε) is estimated as follow: with ϕ := σ(u) · ∇ s u, where we have used the Cauchy-Schwartz inequality and the elliptic regularity of u. From the fact that the stress σ y (w) is uniform inside the inclusion and using the same arguments, we have the following estimate for the second remainder term: Once again, from the Cauchy-Schwartz inequality and the elliptic regularity of u we have the estimate bellow for the third remainder E 3 (ε) = Bε σ( u ε ) · ∇ s u = Bε ∇ s u ε · σ(u) ≤ ∇ s u ε L 2 (Bε) σ(u) L 2 (Bε) ≤ C 4 ε ∇ s u ε L 2 (Bε) .

(.24)
In addition, note that the right-hand side of (.15) depends explicitly on the small parameter ε. Since this problem is linear and in view of Lemma 4, we can write u ε = εϕ 0 . Therefore Now, using again the elliptic regularity of u and the continuity of function b at the point x ∈ Ω, the remainder term E 4 is estimate as follows Finally, the last remainder E 5 (ε) is defined as: Therefore, from the Cauchy-Schwartz inequality and Lemma 1 we obtain where we have considered the fact that the function b is assumed to be constant in the neighborhood of the point x ∈ Ω. By making use of the Hölder inequality together with the Sobolev embbeding theorem, there is where 1/p + 1/q = 1, with q > 1, and δ = 1/q. Therefore Bε b · (u ε − u) ≤ C 10 ε 1+δ u ε − u H 1 (Ω) ≤ C 11 ε 2+δ = o(ε 2 ), ( where we have used Lemma 1 and the fact that 0 < δ < 1.