Graded-material design based on phase-field and topology optimization

(will be inserted by the editor)

Graded-material Design based on Phase-eld and Topology Optimization

Massimo Carraturo^1,2 · Elisabetta Rocca^3,4 · Elena Bonetti^4,5 · Dietmar Hömberg^6,7 · Alessandro Reali¹ · Ferdinando Auricchio¹

Received: 19 November 2019 / Revised version: June 3, 2019

Abstract In the present work we introduce a novel graded-material design based on phase-eld and topol- ogy optimization. The main novelty of this work comes from the introduction of an additional phase-eld vari- able in the classical single-material phase-eld topology optimization algorithm. This new variable is used to grade the material properties in a continuous fashion.

Two dierent numerical examples are discussed, in both of them, we perform sensitivity studies to asses the ef- fects of dierent model parameters onto the resulting structure. From the presented results we can observe that the proposed algorithm adds additional freedom in the design, exploiting the higher exibility coming from additive manufacturing technology.

Keywords Phase-eld·functionally graded material· multi-material design·topology optimization·additive manufacturing

Massimo Carraturo

[email protected]

1 Dipartimento di Ingegneria Civile ed Architettura (DICAr), Universitá degli Studi Pavia, via Ferrata 3, 27100 Pavia, Italy

2Chair for Computation in Engineering, Technical University of Munich, Arcisstr. 21, 80333 Munich, Germany

3Dipartimento di Matematica, Universitá degli Studi Pavia, via Ferrata 5, 27100 Pavia, Italy

4IMATI-CNR, via Ferrata 1, 27100 Pavia, Italy

5 Dipartimento di Matematica "F.Enriques", Universitá degli Studi di Milano, via Saldini 50, 20133 Milano, Italy

6 Weierstrass Institute for Applied Analysis and Stochastics, Mohrenstr. 39, 10117 Berlin, Germany

7Department of Mathematical Sciences, NTNU, Alfred Getz vei 1, 7491 Trondheim, Norway

1 Introduction

Structural topology optimization (TO) [1] is a numeri- cal method which aims, by means of a density function, at optimally distributing a limited amount of mate- rial within a volume, representing the initial geometry of a body undergoing specic loads and displacement boundary conditions.

Structural TO was originally introduced as a discrete formulation where areas of dense material and voids are alternated without any transition region [2]. This rst approach, also known as the 0-1 topology optimization problem, leads to many diculties from both an ana- lytical and a numerical point of view [3].

A possible alternative approach is based on homog- enization methods, where the macroscopic material properties are obtained from microscopic porous ma- terial characteristics [4, 5]. Such a strategy leads to op- timized structures with large grey-scale regions of per- forated, porous material, which are in general undesired due to their elevated manufacturing complexity, in par- ticular when classical manufacturing processes, such as milling or molding are adopted. To obtain a clear black and white design, Solid Isotropic Material Penalization (SIMP) method has been introduced in [6]. The SIMP method consists of penalizing the density region, dier- ent from the void or bulk material, by choosing a suit- able interpolation scheme for material properties at the macroscopic scale [7, 8]. This approach has been suc- cessfully employed in many engineering applications be- side of structural problems, e.g. uid analysis [9], uid- and acoustic-structure interaction [10, 11], heat con- duction [12], multi-physics [13], and composite struc- tures [14].

An alternative to the SIMP method is a TO based on the phase-eld method, for the rst time intro- duced by Bourdin and Chambolle [15]. Successively, this method has been employed by Burger and Stainko [16]

for stress constrained problems and by

Takezawa et al. [17] in the shape and topology opti- mization context for minimum compliance and eigen- frequency maximization problems. More recently, Pen- zler et al. [18] have solved nonlinear elastic problems by means of the phase-eld approach, while Dedè et al. [19]

have been the rst to apply this method in the context of isogeometric analysis. Similarly to the SIMP method, phase-eld based TO penalizes an approximation of the interface perimeter, such that, by choosing a very small positive penalty term, one can obtain a sharp interface region separating solid materials and voids [20].

Inspired by the aforementioned works on phase-eld and TO, we aim here at developing an optimization procedure to obtain structures which exploit the pos- sibility of additive manufacturing (AM) technology to distribute material through a body locally varying the material density, leading to the so called functionally graded materials (FGM). Numerical models to simulate manufacturing processes of FGM products have been proposed and validated in [2125] for dierent AM tech- nologies and materials. All these contributions show how the manufacturing process plays a crucial role in the resulting mechanical properties of FGM structures.

Since including the inuence of the process parameters within topology optimization schemes is an extremely challenging task (cf. [26] for a rst attempt in this direc- tion) we prefer to neglect the inuence of AM process on the optimized structure. The presented results have to be thus considered only as an initial starting point towards the nal design of optimized FGM structures.

Topology optimization routines to obtain FGM design have been recently investigated by many researchers (see, e.g. [2729]). Nevertheless, at the best knowledge of the authors a phase-eld approach has not been em- ployed to obtain optimized graded-material structures yet. The phase-eld method proposed in this contribu- tion allows deriving a rigorous mathematical analysis of the problem for FGM structures and to obtain a clear separation between areas of material and voids.

Moreover, this work aims at investigating by means of 2D examples the sensitivity of the proposed phase-eld method to various numerical and physical parameters.

On one hand, the choice of considering only plane stress, two-dimensional structures may limit the validity of the presented results while, on the other hand, it allows us to gain a clear insight into the sensitivity of the method w.r.t. the considered parameters with a negligible com- putational eort.

Even if our approach can potentially include a multi- material case, in this work we consider only a graded- material design, i.e., where a single material is gradually distributed through the body. The result of such an op- timization routine is a structure with graded stiness values, i.e., a material with stiness continuously vary- ing within the body, alternating regions of soft mate- rial with other regions of stier material. The approach proposed in the present paper reintroduces the typical grey-scale regions of early topology optimization meth- ods but within a controlled and numerically stable for- mulation. This choice is justied by the fact that mod- ern AM technologies allow grading the density of a body in an almost continuous fashion, varying the amount of distributed material point-by-point during the printing process.

The outline of the work at hand is organized as fol- lows. In Sec. 2 we recall the formulation for a single- material phase-eld based TO. Sec. 3 introduces the novel phase-eld approach for graded material struc- tures. Then, Sec. 4 discusses two-dimensional numeri- cal examples, carrying out sensitivity studies for dier- ent choices of numerical and physical parameters and presenting also a possible solution to obtain a manu- factured product from the numerical results. Finally, in Sec. 5 , we draw our conclusions on the present work.

2 Single-material phase-eld topology optimization

In this section we recall the classical formulation for a phase-eld TO of a single-material homogeneous struc- ture, closely following [30].

2.1 State equations

We consider a domain Ω ⊂R^d where material is dis- tributed by means of a scalar phase-eld variableφ, rep- resenting a material density fraction, hence φ ∈ [0,1]

with φ ≡ 0 corresponding to voids (i.e., no material) and φ ≡1 to bulk material. Adopting a linear elastic model, the state equations are as follows:

−div(σ) =0 in Ω (1)

u=0 on Γ_D (2)

σ·n=g on Γ_N (3)

with σ = σ(φ) = C(φ) : ε(u), C = C(φ) the fourth- order linear material tensor, u the displacement eld vector,ε(u)the symmetric strain dened asε=∇^Su=

∇u+∇u^T

/2, g the external load on the boundary Γ_N ⊂ ∂Ω, n the unit normal vector, and Γ_D ⊂ ∂Ω,

|Γ_D|6= 0, the portion of the boundary where homoge- neous Dirichlet boundary conditions are applied.

Assuming the material tensor C to depend on φ, the solution of problem Eq.s (1) to (3) depends on the dis- tribution of the scalar eldφ(i.e.,u=u(φ)). We treat the void as a very soft material, adopting the following expression forC:

C(φ) =Cbulkφ^p+Cvoid(1−φ)^p

where Cbulk is the positive denite material tensor of the bulk, dense material, Cvoid is the positive denite material tensor of an idealized very soft material (rep- resenting the voids), and pcan be any positive value;

for simplicity, we assumeCvoid=γ²Cbulk, withγ1, while, following [7], we set p= 3.

The weak form of the linear elastic problem Eq.s (1) to (3) can be written as:

Z

Ω

σ(φ) :ε(v)dΩ= Z

ΓN

g·vdΓ. (4) withv∈ H¹_D(Ω)a virtual displacement eld. Referring to [30] we can prove that for any given g ∈ L²(Γ_N) and φ ∈ L^∞(Ω), there exists a unique u ∈ H_D¹(Ω) fullling Eq. (4) , withH_D¹(Ω) :={v∈H¹(Ω) : v= 0 onΓ_D}.

2.2 Single-material topology optimization as a minimization problem

The goal of our TO process is to properly minimize the compliance of a given structure, by optimally distribut- ing a limited amount of material.

To properly minimize the compliance, we introduce an objective functionalJ(φ,u(φ))dened as:

J(φ,u(φ)) = Z

ΓN

g·u(φ)dΓ +κ Z

Ω

γ

2 k ∇φk²+1 γψ₀(φ)

dΩ (5) where the rst integral represents a measure of the global system compliance, dened as the inverse of the stiness, while, assuming κ >0 and a double-well po- tential functionψ₀(φ) = (φ−φ²)², the second integral is an approximation of the perimeter of the interfaces between regions withφ= 0 and φ= 1. In Eq. (5) γ corresponds to the thickness of the diuse interface, i.e., the region where 0< φ < 1, the termγ/2 | ∇φ|² pe- nalizes jumps betweenφ= 0andφ= 1, whileψ₀(φ)/γ represents the double-well potential function penalizing phases with φdierent from 0 and 1. We remark that

following [30] we choose the same scaling parameterγ to penalize the sharp interface region and to dene the void soft material; this choice is justied by the assump- tion that when one of the two values goes to zero also the other one has to vanish.

The minimization of the functional in Eq. (5) is im- posed under the assumption of distributing a limited constant quantity of material inside the domain, hence, we introduce the constraint:

Z

Ω

φdΩ=m|Ω|

with 0< m ≤1 representing a target domain volume fraction. Clearly, the displacement eldu(φ)solving the TO problem should also be the solution of the linear elastic problem of Eq. (4) .

In conclusion, the minimization problem we aim to solve is the following.

Problem(P): minφ J(φ,u(φ))

such that the following constraints are satised:

Z

Ω

σ(φ) :ε(v)dΩ= Z

ΓN

g·vdΓ. (6) M(φ) =

Z

Ω

φdΩ−m|Ω|= 0, (7) withφ∈H¹(Ω)satisfying the constraint:

0≤φ≤1 a.e. inΩ. (8)

Following the argument by [30], we can prove that the minimum constrained problem(P)has at least one so- lution (cf. [30, Thm. 4.1]). In particular, to solve prob- lem(P)we introduce the Lagrangian functionalL, de- ned as:

L(φ,u, λ,p) =J(φ,u) +λM(φ) +S(φ,u,p), (9) whereλis the Lagrange multiplier introduced to impose the volume constrain of Eq. (7) and the operatorS is dened as:

S(φ,u,p) = Z

Ω

σ(φ) :ε(p)dΩ− Z

ΓN

g·pdΓ,

which we introduce together with the adjoint variable p. The solution of problem(P)is equivalent to the mini- mization of Eq. (9) subjected to constraint in Eq. (8) ; this last problem can be seen as an optimal control

problem, with solutions( ¯φ,u,¯ ¯λ,p)¯ that have to satisfy the rst order optimality conditions dened by:

D_uL φ,¯ u,¯ λ,¯ p¯

= 0, D_pL φ,¯ u,¯ λ,¯ p¯

= 0, D_λL φ,¯ u,¯ λ,¯ p¯

= 0, DφL φ,¯ u,¯ ¯λ,p¯

φ−φ¯

≥0 ∀φ∈Φad,

where Φ_ad is the set of admissible controls dened as follows:

Φ_ad:={φ∈H_D¹(Ω) : 0≤φ≤1 a.e. inΩ}.

We also note that for the problem under investiga- tionD_pL=D_uL, hence the so-called adjoint equation (holding true for everyv∈H_D¹(Ω)):

Z

Ω

C(φ)ε(¯p) :ε(v)dΩ= Z

ΓN

g·vdΓ,

is identical to the weak form of the linear elastic prob- lem Eq. (4) , which implies thatp¯ = ¯u. We refer to [30]

for the complex analysis of optimality conditions.

To obtain a more compact formulation, we dene here the energy density of the system and its derivative w.r.t.

the scalar eldφas:

E(φ,u) =σ(φ) :ε(u),

and

∂E(φ,u)

∂φ =∂σ(φ)

∂φ :ε(u), where

∂σ(φ)

∂φ = ∂C(φ)

∂φ :ε(u).

To discretize our continuous problem we employ a gradient ow dynamics, namely Allen-Cahn gradient ow [31], a steepest descent pseudo-time stepping method with a time-step incrementτ. Thus the optimal control problem (P)can be now rewritten as follows:

D_uLv= 0, (10)

D_λLv_λ=Mvλ= 0, (11)

γ τ

Z

Ω

(φ_n+1−φ_n)v_φdx=−D_φLv_φ, (12) where

D_φL=∂J

∂φ +λ∂M

∂φ +∂S

∂φ, withv_λ∈Randv_φ∈Φ_ad.

The problem dened in Eq.s (10) to (12) can be written in the following weak extended formulation:

Z

Ω

σ(φ) :ε(v)dΩ= Z

ΓN

g·vdΓ, (13) Z

Ω

v_λ(φ_n+1−m)dΩ= 0, (14) γ

τ Z

Ω

(φ_n+1−φ_n)v_φdΩ+κγ Z

Ω

∇φ_n+1· ∇v_φdΩ +λ

Z

Ω

v_φdΩ− Z

Ω

v_φ∂E(φ_n,u_n)

∂φ dΩ +κ

γ Z

Ω

∂ψ0(φn)

∂φ v_φdΩ= 0. (15)

2.3 Single-material nite element formulation

We derive here a nite element approximation of the phase-eld TO problem dened in Eq.s (13) to (15).

To this end we discretize the physical domainΩ using two dierent meshesQ_u,Q_φcorresponding to the eld variables uand φ and their variationsv and vφ. The Lagrange multiplier λused to constrain the volume is applied using a constant scalar value on Ω. On each mesh, we interpolate the nodal values of the eld vari- ables and their variations by means of piecewise linear basis functions, such that:

u≈N_uu,˜ v≈N_uv,˜ φ≈Nφφ,˜ vφ≈Nφ˜vφ,

Introducing the proposed discretization in Eq.s (13) to (15) the discrete version of the optimal control prob- lem becomes:

1 τ









0 0 M^φλ 0 M^φφ 0 0 M^λφ 0

















˜ u φ˜ λ˜







 +









K^uu 0 0 0 K^φφ 0 0 0 0

















˜ u φ˜ λ˜









=



 f q^φ+q^s+q^ψ

q^λ





(16)

with the matrix and vector terms dened as follows:

K^uu= Z

Ω

∇NuT

C∇NudΩ, M^φφ =γ

Z

Ω

N^T_φN_φdΩ, K^φφ=κγ

Z

Ω

∇N^T_φ∇NφdΩ, M^λφ =τ

Z

Ω

N_φdΩ= M^φλT , f =

Z

ΓN

N_u^TgdΓ, q^φ= γ

τ Z

Ω

N^T_φNφφ˜_ndΩ=M^φφφ˜_n, q^λ=

Z

Ω

mdΩ, q^s=

Z

Ω

N^T_φ∂E(φ˜_n,˜u_n)

∂φ dΩ, q^ψ=−κ

γ Z

Ω

N^T_φ∂ψ₀(φ˜_n)

∂φ dΩ.

The discrete linear system in Eq. (16) can be solved using a staggered approach, i.e. solving rst the state equation system:

K^uuu˜=f, (17)

and then the discretized optimization problem:

1 τ

"

M^φφ+τK^φφM^φλ M^λφ 0

# "φ˜ λ˜

#

=

"

q^φ+q^s+q^ψ q^λ

# .

(18) Following this approach Alg. 1 solves the problem in two separate steps: rst it solves the state equa- tion system Eq. (17) to get the solution vectoru˜n+1

(line 3), secondly, the linear system of Eq. (18) is solved to obtain the phase-eld vector φ˜^∗_n+1 and the Lagrange multiplier vector λ˜_n+1 (line 4). Finally, the vector φ˜^∗_n+1 is projected within the interval [0,1] to obtain the phase-eld solution vectorφ˜_n+1fullling the constraints (line 5). We use the increment∆_φbased on theL²-norm and dened as:

∆φ=

φ˜_n+1−φ˜_{n L2}

φ˜_{n L2}

, (19)

as a criterion to assert the convergence of the algorithm, which is otherwise stopped when user dened maximum number of iterationsmax_iter is reached.

Algorithm 1: Single-material optimization al- gorithm

input :Qu,Q_φ,φ˜₀. output: Optimal topology

1 φ˜_n←φ˜₀

2 while∆φ≥tolandn≤maxiterdo

3 u˜n+1←solve(17)

4 (˜λn+1,φ˜^∗_n+1)←solve(18)

5 φ˜_n+1←rescale φ˜^∗_n+1

to[0,1]

6 update(∆φ)

7 φ˜_n←φ˜_n+1

8 end

3 Graded-material phase-eld topology optimization

In the following section, we extend the previously pre- sented formulation of TO to the case of a graded ma- terial denition. We refer to this approach as graded- material phase-eld topology optimization. The math- ematical analysis of the corresponding optimization problem will be the subject of a forthcoming paper [32].

3.1 State equation

We now consider the case of an inhomogeneous material distribution; in particular, we assume that the material elastic fourth-order tensorCcan vary linearly through a material grading scalar variableχ ∈[0, φ], with φ∈ [0,1], such that:

C(χ) =Cbulkχ+ 1

βCbulk(φ−χ),

with 0 < β ≤ 1 a so-called softening factor, used to dene the soft material tensor as a fraction of the bulk, fully dense material tensor Cbulk. In such a way, the stiness of the body can continuously vary from a full dense material (χ=φ) to a softer one (χ= 0). There- fore, the denition of the fourth-order material tensor C(φ, χ), previously dened in Sec. 2 , can now be mod- ied as:

C(φ, χ) =C(χ)φ^p+γ²_φC(χ)(1−φ)^p, (20) where 0< γ_φ 1 and again we choose a penalty pa- rameterp= 3.

We would like to remark here that FGM structures are intrinsically heterogeneous but, as recently demon- strated by Cheng et al. [33], an asymptotic homoge- nization method can be eectively employed. In this

work, we assume a simple linearly interpolated homog- enized material tensor but a more complex homogeniza- tion could be directly introduced within the proposed numerical scheme.

The denition of the material tensor in Eq. (20) leads to an optimized structure where, as in the single- material case, the perimeter of the body is dened by the sharp interface of the phase-eld variable φ, while the stiness of the material continuously varies within the structure, following the distribution of the material grading variable χ. Hence, the graded-material weak form of the linear elastic problem of Eq. (4) reads:

Z

Ω

σ(φ, χ) :ε(v)dΩ= Z

ΓN

g·vdΓ (21) with the virtual displacement v ∈ H_D¹(Ω) and where σ(φ, χ) =C(φ, χ) :ε(u).

3.2 Graded-material topology optimization as a minimization problem

We want now to dene an objective functional which optimizes a structure with an inhomogeneous material distribution. This new graded-material objective func- tionalJ^M(φ, χ,u(φ, χ))can be dened as:

J^M(φ, χ,u(φ, χ)) = Z

ΓN

g·u(φ, χ)dΓ+ κ_φ

Z

Ω

γ_φ

2 | ∇φ|²+1 γ_φψ₀(φ)

dΩ+κ_χ Z

Ω

γ_χ

2 | ∇χ|²dΩ, with κ_φ, κ_χ > 0 and γ_χ > 0, and where the rst two integrals are the same of the objective func- tional in Eq. (5) , while the additional integral term γ_χ/2 (| ∇χ|)² is introduced to penalize the gra- dient of the scalar eld χ.

Following the same approach described for the single- material case, the global graded-material minimization problem can be now written as follows:

minφ,χ J^M(φ, χ,u(φ, χ)), such that:

Z

Ω

σ(φ, χ) :ε(v)dΩ= Z

ΓN

g·vdΓ, M(φ) =

Z

Ω

φdΩ−m|Ω|= 0, whereφ, χ∈H¹(Ω), under the constraint 0≤φ≤1 a.e. inΩ,

and the additional constraint on χ: 0≤χ≤φ a.e. inΩ.

We can now dene the graded-material LagrangianL^M as:

L^M =J^M +λM+S^M, explicitly written as:

L^M(φ, χ,u, λ,p) =J^M(φ, χ,u) +λM(φ) +S^M(φ, χ,u,p), where, the operatorS^M for the graded-material formu- lation is calculated as:

S^M(φ, χ,u,p) = Z

Ω

σ(φ, χ) :ε(p)dΩ− Z

ΓN

g·pdΓ.

Analogously to the previously introduced set of admissi- ble controlsΦ_adfor the phase-eld variableφ, we dene now the set of admissible controlsΞ_ad for the grading variableχas:

Ξad:={χ∈H¹(Ω) : 0≤χ≤φ a.e. inΩ}.

Clearly, also in the graded-material case, we want that the optimal control solutions φ¯ and χ¯ have to satisfy the rst order necessary optimality conditions, which can be derived as:

D_φL^M( ¯φ,χ,¯ u,¯ λ,¯ p)¯ φ−φ¯

≥0 ∀φ∈Φ_ad and

D_χL^M φ,¯ χ,¯ u,¯ λ,¯ p¯

(χ−χ)¯ ≥0 ∀χ∈Ξ_ad,

where u¯ and p¯ are solutions of the graded-material state equation Eq. (21) and of the corresponding ad- joint problem, respectively. As in the previous case, the displacement eldu is self-adjoint and hence we have p¯ = ¯u. For a complete analysis of necessary rst or- der optimality conditions we refer to the forthcoming paper [32].

Analogously to the single-material case, we can dene the energy density of the system and its derivatives w.r.t. both the scalar eldφ and the material grading variableχas:

E^M(φ, χ,u) =σ(φ, χ) :ε(u),

∂E^M(φ, χ,u)

∂φ = ∂σ(φ, χ)

∂φ :ε(u) = 3C(χ)φ²+ 3γ_φ²C(χ)(1−φ)²

ε(u) :ε(u) and

∂E^M(φ, χ,u)

∂χ = ∂σ(φ, χ)

∂χ :ε(u) =

Cbulk−1 βCbulk

φ³+γ_φ²(1−φ)³

ε(u) :ε(u).

The optimal control problem can be solved as in the single-material case by means of the Allen-Cahn gradi- ent ow, leading to the following set of equations:

γ_φ τ

Z

Ω

(φn+1−φn)vφdΩ+κφγφ

Z

Ω

∇φ· ∇vφdΩ+

Z

Ω

v_φλdΩ− Z

Ω

v_φ∂E^M(φ_n, χ_n,u_n)

∂φ dΩ+

κφ

γ_φ Z

Ω

∂ψ0(φn)

∂φ v_φdΩ= 0, (22) γ_χ

τ Z

Ω

(χ_n+1−χ_n)v_χdΩ+κ_χγ_χ Z

Ω

∇χ· ∇v_χdΩ−

Z

Ω

v_χ∂E^M(φ_n, χ_n,u_n)

∂χ dΩ= 0, (23) to be solved under the volume constraint

Z

Ω

v_λ(φ−m)dΩ= 0. (24)

In order to estimate the total amount of material in the structure, we dene a material fraction indexmχ as:

m_χ = 1

|Ω| Z

Ω

χdΩ,

which can be considered as a measure of the global amount of material used to print the structure. The equivalent material fraction index for the single- material case m_φ is equal to the volume fraction m, such that:

m_φ=m= 1

|Ω| Z

Ω

φdΩ.

3.3 Graded-material nite element formulation

We aim now at obtaining a discrete formulation for the graded-material phase-eld topology optimization problem. To this end, the displacement eld u, the phase-eld variable φ, the Lagrange multiplier λ and their corresponding variations are approximated using the same discretization already dened in Subsec. 2.3 . Additionally, we need to discretize the material grad- ing variableχon the domainΩ; such a discretization is obtained introducing an additional meshQ_χ, such that the material grading variableχand its variationv_χ can be written as:

χ≈N_χχ˜ and v_χ≈N_vχv˜_χ, whereN_χ andN_vχ are the piecewise linear shape func- tions which interpolate the nodal degrees of freedoms

˜

χandv˜_χ, respectively.

The discrete form of Eq.s (22) to (24) can thus be writ- ten in a compact notation as:

1 τ









0 0 0 0 0 M^φφ 0 M^φλ 0 0 M^χχ 0 0 M^λφ 0 0















 u˜ φ˜

˜ χ λ˜







 +









K^uu 0 0 0 0 K^φφ 0 0 0 0 K^χχ0 0 0 0 0















 u˜ φ˜

˜ χ λ˜









=







 f q^φ+q^s0+q^ψ

q^χ+q^t q^λ







 , (25)

where the newly dened matrix and vector terms are:

M^χχ=γ_χ Z

Ω

N^T_χN_χdΩ, K^χχ=κ_χγ_χ

Z

Ω

∇N^T_χ∇N_χdΩ, q^χ= γχ

τ Z

Ω

N^T_χNχχ˜_ndΩ, q^s0 =

Z

Ω

N^T_φ∂E^M(φ˜_n,χ˜_n,˜u_n)

∂φ dΩ,

q^t= Z

Ω

N^T_φ∂E^M(φ˜_n,χ˜_n,˜u_n)

∂χ dΩ.

Alg. 2 describes the iterative procedure to obtain the graded-material optimized structure discussed so far.

The adopted solution scheme is very similar to Alg. 1 but in this case, we have to solve at each iteration the graded-material linear system dened in Eq. (25) to obtain the phase-eld solution vector φ˜_n+1 and the grading scalar variable vector χ˜_n+1. As in the single material case, the system can be solved following a stag- gered scheme, since Eq. (25) can be split into two separate systems as follows:

K^uuu˜=f, (26)

and

1 τ









M^φφ+τK^φφ 0 M^φλ 0 M^χχ+τK^χχ 0

M^λφ 0 0















 φ˜

χ˜

˜λ









=









q^φ+q^s0+q^ψ q^χ+q^t

q^λ







 . (27)

Finally, we use the relative increment of χ in the L²- norm, dened as:

∆_χ =

χ˜_n+1−χ˜_n L²

kχ˜_nk_L2

. (28)

as an additional criteria to check the convergence of Alg. 2 .

Algorithm 2: Graded-material optimization algorithm

input :Qu,Q_φ,Qχ,φ₀,χ₀ output: Optimal topology

1 φ_n←φ₀

2 χ_n←χ₀

3 while (∆φ≥tolor∆χ≥tol) andn≤maxiterdo

4 ˜un+1←solve(26)

5 (φ^˜∗_n+1,χ˜^∗_n+1,˜_λn+1)←solve(27)

6 φ˜_n+1←rescale φ˜^∗_n+1

to[0,1]

7 χ˜_n+1←rescale χ˜^∗_n+1 to[0, φ]

8 update(∆φand∆χ)

9 φ_n←φ_n+1

10 χ_n←χ_n+1

11 end

4 Numerical Examples

In this section two numerical examples are presented: in the rst one, we consider a cantilever beam structure while in the second one we study a simply-supported beam structure. For each example, we decided to run two sensitivity studies w.r.t. a numerical and a physi- cal problem parameter to asses the robustness of Alg. 2 to these changes. In the cantilever beam example, we discuss two sensitivity studies. Firstly, we vary the val- ues of the graded-material interface parameterγ_χ (i.e., the parameter which represents the thickness of the material grading variable interface) and, secondly, we change the slenderness of the structure. In the simply- supported beam example, again we perform two stud- ies. In the rst one, we use dierent values for the soft- ening factorβ while in the second one we increase the load acting on the structure. Finally, in Subsec. 4.3 we present a possible solution to obtain an AM prod- uct from our numerical results.

4.1 Cantilever beam

We consider the cantilever beam problem depicted in Fig. 1 , with dimensions a = 2mm and b = 1mm and a traction forceg= (0,−600)N/mm applied at the right-end of the lower edge of the structure, while the left edge is xed. We assume the initial material being a dense isotropic material, i.e.,Cbulk= (λ+2µ)1⊗1+2µI, where the Lame's parametersλandµcan be expressed in terms of the Young modulusE and the Poisson co- ecientν as follows:

λ= Eν

(1 +ν)(1−2ν), (29)

x y

Γ

D

Ω

a

b Γ

N

g

Fig. 1: Cantilever beam: Initial conguration and prob- lem domain.

and

µ= E

2(1 +ν). (30)

The softening factor β is chosen equal to 4, i.e., the soft material tensor is four times softer thanCbulk. We choose a dense material havingE = 12.5GPa andν = 0.25. We discretize the domain Ω using a mesh with 128×64quadrilateral elements and we set m = 0.45, κ_χ =κ_φ = 4, γ_φ = 0.02, a time step increment ∆τ = 1.0×10⁻⁶,φ0= 0.5as initial solution, and a tolerance equal to 0.01.

4.1.1 Sensitivity study of the graded-material interface parametersγ_χ

In this rst sensitivity study, we investigate the dier- ent topologies obtained by varying γχ between 0.001 and0.1, as reported in Fig. 2 . The results show that the optimal multi-material distribution is very dierent from the single-material optimized topology depicted in Fig. 3 , for values ofγ_χ smaller thanγ_φ. In fact, in this case, the voids present in the single-material struc- ture are replaced by areas of soft material. Contrary, if γ_χ is chosen to be bigger thanγ_φ the solution presents void regions similarly to the single-material case. Fi- nally, we observe that, as expected, when the thick- ness of the diuse interface is too small compared to the element size, the solution does not converge any- more ( Fig. 2a ). Table 1 reports the values of the compliance and of the material fraction index m_χ for dierent values ofγ_χ. From the values in Table 1 , we can see that employing a softer material will increase the compliance of the body, leading at the same time to lighter structures compared to the homogeneous mate- rial case.

(a)γχ= 0.001 _(b)γχ= 0.005

(c)γχ= 0.010 _(d)γχ= 0.020

(e)γχ= 0.050 _(f)γχ= 0.100 Fig. 2: Cantilever beam: Sensitivity study of the multi- material interface parameters γ_χ. If γ_χ ≤ γ_φ = 0.02 dierent values of the graded-material interface param- eters do not aect too much the nal solution, which presents a wide region of soft material lling the voids which, instead, characterize the single-material solution (see Fig. 3 ); whereas, for γ_χ ≥γ_φ, the nal solution of the graded-material case tends to a single-material conguration with multiple holes in the structure. Fi- nally, if we choose γ_χ too small with respect to our element size, the solution does not converge anymore as in Fig. 2a .

4.1.2 Sensitivity study of the slenderness of the structure

On the cantilever beam, we perform a second sensitiv- ity study varying the slenderness ratio s = a/b (i.e., the ratio between the length and the height of the can- tilever beam), for a xed value of the graded-material interface parameter (γ_χ= 0.02). Fig. 4 shows the nal

Fig. 3: Cantilever beam: Sensitivity study of the graded- material interface parametersγ_χ. Single-material opti- mized structure.

Table 1: Cantilever beam: Sensitivity study of the graded-material interface parameters γ_χ. Compliance and material index values for dierent choices ofγ_χ.

γχ compliance mχ convergence

0.001 105.3 0.380 _NO

0.005 122.9 0.265 _YES

0.01 133.0 0.245 _YES

0.02 141.9 0.230 _YES

0.05 154.0 0.225 _YES

0.1 165.4 0.201 _YES

full dense material 52.3 mφ= 0.45 _YES

topologies for three dierent slenderness ratios, where all the resulting structures are characterized by internal regions of softer material and external support of stier material.

4.2 Simply-supported beam

In this second example, we choose Acrylonitrile Buta- diene Styrene (ABS), which is a common thermoplas- tic polymer widely used in 3D printing applications, as material to obtain an optimized simply-supported beam structure. The problem is symmetric and thus we decide to solve only half of the domain as depicted in Fig. 5 , whereh= 50mm anL/2 = 100mm, with a distributed external loadgequal to(0,−50)N/mm applied on the top edge of the structure. The Young modulus and the Poisson coecient of ABS plastic are 2.3GPa and 0.35, respectively. We set m = 0.4, κ_φ = κ_χ = 1 and γφ = γχ = 0.01, while we choose a pseudo-time step

∆τ = 1.0×10⁻⁶ and an initial solutionφ₀= 0.5. 4.2.1 Sensitivity study of the softening factor β Fig. 6 presents the results of a sensitivity analysis per- formed varying the softening factorβ from 1 to 4. The resulting optimized structures show that, introducing

(a)s= 1 _(b)s= 2

(c)s= 4

Fig. 4: Cantilever beam: Sensitivity study of the slan- derness of the structures. Varying the slanderness ratio s = a/b we obtain optimized structures characterized by an outer frame of sti material lled with regions of soft materials.

Table 2: Simply-supported beam: Sensitivity study of the softening factorβ. Compliance and material index values for dierent choices ofβ.

β compliance mχ convergence 1 20.5 0.40 _YES 2 37.3 0.32 _YES 3 46.4 0.24 _YES 4 58.6 0.18 _YES

grey-scale regions in the structure, the optimal design is modied, replacing the typical voids of SIMP ap- proach with areas of soft material. Again we observe in Table 2 that introducing a soft material within the algorithm leads to structure with a smaller material in- dex but higher compliance. The values of the softening factor strongly inuences the nal results and give us the possibility to obtain intermediate structure such as the one in Fig. 6b . Moreover, we can notice that for high values ofβthe results are very similar to each other (see Fig. 6c and Fig. 6d ), thus the higher values of this parameter would depend only on the technological boundaries of the AM process.

x y

^κD

◦ ◦

◦ ◦

◦ ◦ Ω

L/2

h

g

κN

Fig. 5: Simply-supported beam: Initial conguration and problem domain.

(a)β= 1 _(b)β= 2

(c)β= 3 _(d)β= 4

Fig. 6: Simply-supported beam: Sensitivity study of the softening factorβ. Increasing the values of the softening factor, i.e., employing a softer material, the optimized structure does not present anymore the typical holes resulting from a single-material optimization 6a. Voids are now replaced by a region of soft material.

4.2.2 Sensitivity study of the distributed load

On the simply-supported beam, we conduct a second sensitivity study xing β = 3 and increasing the dis- tributed loadgby a factor of 2 and 3, respectively. The resulting structures are reported in Fig. 7 . As we ex-

pected employing a heavier load reduces the areas of soft material, increasing at the same time the number of columnar structures in the nal topology. We want to remark here that the structure of Fig. 7b did not converge even after 1000 iterations. Since the mesh is not modied, this behavior is due to the choice of the stier material, which is in this case too soft for such a heavy load.

(a)g×2 _(b)g×3

Fig. 7: Simply-supported beam: Sensitivity study of the distributed load g. Increasing the load on the upper edge of the structure we observe an increment in the region occupied by full dense material and the presence of holes and columnar structures.

4.3 From numerics to AM products

In order to demonstrate that the numerical results pre- sented so far can be actually realized in practice, we decided to print a sample of the optimized structure depicted in Fig. 6d . To obtain this result we used the Fused Deposition Modeling (FDM) 3D printer present in our PROTOtyping LABoratory (PROTOLAB) at the University of Pavia. This sample is made in ABS plastic material and realized extruding by an oset of 5mm the structure of Fig. 6d . The values ofχare then mapped onto a manufacturing grid where the dimension of each cell depends on technological constraints given by the machine. Finally, by means of a simple boolean operation, we subtract a quadratic region proportional to the average value of χ in each cell. Once the cor- responding CAD model is completed we can directly print it, obtaining the structure of Fig. 8 . Since, the assumption of linear elasticity is not valid for a soft material such as plastic (metal alloys would be a more appropriate choice in this case) we leave the experi- mental validation of the proposed numerical algorithm to future research.

Fig. 8: Simply-supported beam with a distributed load of Fig. 6d printed using FDM 3D printer.

5 Conclusions

In the present work, we have introduced a novel phase- eld topology optimization algorithm based on a graded material denition.

The numerical results show that the additional control parameterχ, introduced in our phase-eld formulation, allows increasing the number of possible optimal de- signs delivered by the topology optimization process.

In particular, we have introduced the possibility to con- trol the distribution of the material density within our structure in a continuous fashion. Such a feature can be in many cases highly desirable, in particular if we consider additive manufacturing applications.

Moreover the algorithm allows to easily control the number of regions with graded material distribution, delivering results which can be in between a fully black- and-white approach and a purely graded-material dis- tribution.

In the near future we aim at investigating mechanical properties of 3D printed structures designed using the graded-material phase-eld topology optimization algo- rithm.

Acknowledgements This work was partially supported by Re- gione Lombardia through the project "TPro.SL - Tech Proles for Smart Living" (No. 379384) within the Smart Living pro- gram, and through the project "MADE4LO - Metal ADditivE for LOmbardy" (No. 240963) within the POR FESR 2014-2020 program. MC and AR have been partially supported by Fon- dazione Cariplo - Regione Lombardia through the project Verso nuovi strumenti di simulazione super veloci ed accurati basati sull'analisi isogeometrica, within the program RST - raorza- mento. This research has been performed in the framework of the project Fondazione Cariplo-Regione Lombardia MEGAsTAR Matematica d'Eccellenza in biologia ed ingegneria come accel- eratore di una nuova strateGia per l'ATtRattività dell'ateneo

pavese. The present paper also benets from the support of the GNAMPA (Gruppo Nazionale per l'Analisi Matematica, la Probabilità e le loro Applicazioni) of INdAM (Istituto Nazionale di Alta Matematica) for ER. A grateful acknowledgment goes to Dr. Ing. Gianluca Alaimo for his support and precious suggestions on additive manufacturing technology.

References

1. M. P. Bendsøe and N. Kikuchi. Generating optimal topolo- gies in structural design using a homogenization method.

Comput. Methods Appl. Mech. Engrg., 71:197224, 1988.

2. M. P. Bendsøe. On obtaining a solution to optimization prob- lems for solid, elastic plates by restriction of the design space.

J. Struct. Mech., 11(4):501521, 1983.

3. O. Sigmund and J. Petersson. Numerical instabilities in topology optimization: A survey on procedures dealing with cheackboards, mesh-dependencies and local minima. Struc- tural Optimization, 16:6875, 1998.

4. G Allaire, F. Jouve, and H. Maillot. Topology optimization and optimal shape design using homogenization. Struct. Mul- tidisc. Optim., 28:8798, 2004.

5. Suzuki. K and N. Kikuchi. A homogenization method for shape and topology optimization. Comput. Methods Appl.

Mech. Engrg., 93(3):291318, 1991.

6. M. Zhou and G. I. N. Rozvany. The coc algorithm, part ii: Topological geometry and generalized shape optimization.

Comp. Meth. Appl. Mech. Engng., 89:197224, 1991.

7. M. P. Bendsøe and O. Sigmund. Material interpolation schemes in topology optimization. Archive of Applied Me- chanics 6, 65:635654, 1999.

8. Martin P. Bendsøe and Ole Sigmund. Topology Optimiza- tion - Theory, Methods, and Applications. Springer Verlag, Germany, 2003. ISBN 3-540-42992-1.

9. A. Gersborg-Hansen, M. P. Bendsøe, and O. Sigmund. Topol- ogy optimization of channel ow problems. Struct. Multidis- cip. O., 30(3):181192, 2005.

10. G. H. Yoon. Topology optimization for stationary uid- structure interaction problems using a new monolithic for- mulation. Int. J. Numer. Methods Engrg., 82(5):591616, 2010.

11. G. H. Yoon, J. Jensen, and Sigmund O. Topology optimiza- tion of acousticstructure interaction problems using a mixed nite element formulation. Int. J. Numer. Methods Engrg., 70(9):10491075, 2007.

12. A. Gersborg-Hansen, M. P. Bendsøe, and O. Sigmund. Topol- ogy optimization of heat conduction using the nite volume method. Struct. Multidiscip. O., 31(4):251259, 2006.

13. C. S. Andreasen and O. Sigmund. Topology optimization of uid-structure-interaction problems in poroelasticity. Com- put. Methods Appl. Mech. Engrg., 31(4):5562, 2013.

14. O. Sigmund and S. Torquato. Composites with extremal thermal expansion coecients. Appl. Phys. Lett., 69(21):

32033205, 1996.

15. B. Bourdin and A. Chambolle. Design-dependent loads in topology optimization. ESAIM Contr. Optim.Calc. Var., 9:

1948, 2003.

16. M. Burger and R. Stainko. Phase-eld relaxation of topology optimization with local stress constraints. SIAM Journal on Control and Optimization, 45(4):14471466, 2006.

17. A. Takezawa, S. Nishiwaki, and M. Kitamura. Shape and topology optimization based on the phase eld method and sensitivity analysis. J. Comp. Phys., 229(7):26972718, 2010.

18. P. Penzler, M. Rumpf, and B. Wirth. A phase-eld model and minimal compliance shape optimization in nonlinear elastic- ity. ESAIM Control Optim. Calc. Var., 2012:229258, 2012.

19. Luca Dedè, Micheal J. Borden, and Thomas JR Hughes. Iso- geometric analysis for topology optimization with a phase eld model. Archives of Computational Methods in Engi- neering, 19(3):427465, 2012.

20. Luise Blank, M. Hassan Farshbaf-Shaker, Harald Garcke, Christoph Rupprecht, and Vanessa Styles. Multi-material Phase Field Approach to Structural Topology Optimization.

In Günter Leugering, Peter Benner, Sebastian Engell, An- dreas Griewank, Helmut Harbrecht, Michael Hinze, Rolf Ran- nacher, and Stefan Ulbrich, editors, Trends in PDE Con- strained Optimization, volume 165, pages 231246. Springer International Publishing, Cham, 2014. ISBN 978-3-319- 05082-9 978-3-319-05083-6.

21. Wentao Yan, Wenjun Ge, Jacob Smith, Stephen Lin, Orion L.

Kafka, Feng Lin, and Wing Kam Liu. Multi-scale modeling of electron beam melting of functionally graded materials.

Acta Materialia, 115:403412, 2016.

22. Zhengtao Gan, Hao Liu, Shaoxia Li, Xiuli He, and Gang Yu.

Modeling of thermal behavior and mass transport in multi- layer laser additive manufacturing of Ni-based alloy on cast iron. International Journal of Heat and Mass Transfer, 111:

709722, 2017.

23. Zhengtao Gan, Gang Yu, Xiuli He, and Shaoxia Li. Numeri- cal simulation of thermal behavior and multicomponent mass transfer in direct laser deposition of Co-base alloy on steel.

International Journal of Heat and Mass Transfer, 104:2838, 2017.

24. Kai Ke Yang, Ji Hong Zhu, Chuang Wang, Dong Sheng Jia, Long Long Song, and Wei Hong Zhang. Experimental valida- tion of 3D printed material behaviors and their inuence on the structural topology design. Computational Mechanics, 61:581598, 2018.

25. Sarah J. Wol, Zhengtao Gan, Stephen Lin, Jennifer L. Ben- nett, Wentao Yan, Gregory Hyatt, Kornel F. Ehmann, Gre- gory J. Wagner, Wing Kam Liu, and Jian Cao. Experimen- tally validated predictions of thermal history and microhard- ness in laser-deposited Inconel 718 on carbon steel. Additive Manufacturing, 27:540551, 2019.

26. Gregoire Allaire and Lukas Jakabcin. Taking into account thermal residual stresses in topology optimization of struc- tures built by additive manufacturing. Mathematical Models and Methods in Applied Sciences, 28:23132366, 2018.

27. D. Brackett, I. Ashcroft, and R. Hague. Topology optimiza- tion for additive manufacturing. Solid Freeform Fabrication Symposium (SFF), Austin, 2014.

28. Lin Cheng, Pu Zhang, Emre Biyikli, Jiaxi Bai, Steve Pilz, and Albert C. To. Integration of topology optimization with ecient design of additive manufactured cellular structures.

Solid Freeform Fabrication Symposium (SFF), Austin, 2015.

29. D. Hickman A. Panesar, M. Abdi and I. Ashcroft. Strate- gies for functionally graded lattice structures derived using topology optimisation for additive manufacturing. Additive Manufacturing, 19:8194, 2018.

30. L. Blank, H. Garcke, M.H. Farshbaf-Shaker, and V. Styles.

Relating phase eld and sharp interface approaches to struc- tural topology optimization. ESAIM Control Optim. Calc.

Var., 20:10251058, 2014.

31. Samuel M. Allen and John W. Cahn. A microscopic theory for antiphase boundary motion and its application to an- tiphase domain coarsening. Acta Metallurgica, 27(6):1085 1095, June 1979. ISSN 00016160.

32. F. Auricchio, E. Bonetti, M. Carraturo, D. Hömberg, A. Re- ali, and E. Rocca. Structural multiscale topology optimiza-

tion with stress constraint for additive manufacturing. Work in progress, 2018.

33. Lin Cheng, Jiaxi Bai, and Albert C. To. Functionally graded lattice structure topology optimization for the design of addi- tive manufactured components with stress constraints. Com- puter Methods in Applied Mechanics and Engineering, 344:

334359, 2019.