Overhang Penalization in Additive Manufacturing via Phase Field Structural Topology Optimization with Anisotropic Energies

A phase field approach for structural topology optimization with application to additive manufacturing is analyzed. The main novelty is the penalization of overhangs (regions of the design that require underlying support structures during construction) with anisotropic energy functionals. Convex and non-convex examples are provided, with the latter showcasing oscillatory behavior along the object boundary termed the dripping effect in the literature. We provide a rigorous mathematical analysis for the structural topology optimization problem with convex and non-continuously-differentiable anisotropies, deriving the first order necessary optimality condition using subdifferential calculus. Via formally matched asymptotic expansions we connect our approach with previous works in the literature based on a sharp interface shape optimization description. Finally, we present several numerical results to demonstrate the advantages of our proposed approach in penalizing overhang developments.


Introduction
Additive manufacturing (AM) is an innovative building technique that produces objects in a layer-by-layer fashion through fusing or binding raw materials in powder and resin forms.Since its introduction in the 1970s, it has shown great versatility in allowing for the creation of highly complex geometries, immediate modifications and redesign, thus making it an ideal process for rapid prototyping and testing.But despite such advantages over traditional manufacturing technologies, AM still has not seen widespread integration in serial production, and there are still many limitations that have yet to be overcome.For an overview of the technologies involved and the challenges encountered by AM, we refer the reader to the review article [1] as well as the report [74].
In this work, we focus on a recurring issue encountered when practitioners employ AM to construct objects.Some categories of AM technologies, such as Fused Deposition Modeling and Laser Metal Deposition, construct objects in an upwards direction (hereafter referred to as the build direction) by repeatedly depositing and then fusing a new layer of raw materials on the surface of the object.Overhangs are regions of the constructed object that when placed in a certain orientation extend outwards without any underlying support.For example, the top horizontal bar of a T-shaped structure placed vertically will be classified as an overhang.The main issue with overhangs is that they can deform under their own weight or by thermal residual stresses from the construction process, and if not supported from below, there is a risk that the object itself can display unintended deformations and the entire printing process can fail.
One natural solution is to identify certain orientations of the object whose overhang regions are minimized prior to printing, and choose the build direction to be one of these orientations.Going back to the T-shape structure example, we can simply rotate the shape by 180 • so that there are no overhang regions (see, e.g., [54,Fig. 5]).However, in general, there is no guarantee that orientations with no overhang regions exist.In this direction we mention the works [70,91] that incorporate various geometric information such as contact surface area and support volume expressed as functions of the build direction within an optimization procedure.On the other hand, simultaneous optimization of build direction and topology has been considered in [88].
Another remedy is to employ support structures that are built concurrently with the object whose purpose is to reduce the potential deformation of overhang regions through mechanical loads or thermal residual stresses.These supplementary structures act like scaffolding to overhangs, and after a successful print are then removed in a post-processing step.Along with increased material costs and printing time, there is a risk of damaging delicate features of the objects during the removal step.Nevertheless, for certain AM technologies such as the aforementioned Fused Deposition Modeling, some form of support structures will always be necessary in order to mitigate the potential undesirable deformations of the finished object.Thus, recent mathematical research has concentrated on the optimizing support structures in an effort to reduce material waste and minimize contact area with the object.Many such works explore optimal support in the framework of shape and topology optimization [3,48,56,60,67], as well as proposing sparse cellular designs that have low solid volume and contact area in the form of lattices [53], honeycomb [64], and trees [87].Further details can be found in the review article [54].
A related approach would be to allow some modifications to the object, as long as the altered design retains the intended functionality of the original, leading to creation of self-supporting objects that do not require support structures at all [30,61,62,63].Our present work falls roughly into this category, where we employ a well-known phase field methodology in structural topology optimization with the aim of identifying optimal designs fulfilling the so-called overhang angle constraint.In mathematical terms, let us consider an object Ω 1 ⊂ R d , d = 2, 3, being built within a hold-all rectangular domain Ω = [0, 1] d in a layer-by-layer fashion with build direction e d = (0, . . ., 1) .The base plate is the region B := [0, 1] d−1 × {0} where Ω 1 is assembled on.For any p ∈ ∂Ω 1 \ B, the overhang angle α is defined as the angle from the base plate to the tangent line at p (denoted as α 1 and α 2 in Figure 1).Then, there exists a critical threshold angle ψ where the portion of Ω 1 is self-supporting (resp.requires support structures) if the overhang angle there is greater (resp.smaller) than ψ.Conventional wisdom from practitioners puts the critical angle ψ at 45 • [34, 60, 82] with the reasoning that every layer then has approximately 50% contact with the previous layer and thus it is well-supported from below, although there are some authors that propose smaller values for ψ (e.g., 40 • in [61,62]).Exact values of this critical angle also depend on the setting of the 3D printer as well as physical properties of the raw materials used.We mention that other works [67] consider an alternate definition, which is given as the angle from the vertical build direction to the outward unit normal of ∂Ω 1 .Denoting this as β, a simple calculation shows that β = 180 • − α.Then, regions where β > 180 • − ψ (typically 135 • if ψ = 45 • ) will require support structures.
In this work we adopt the latter definition, but choose to define the overhang angle as the angle measured from the negative build direction −e d to the outer unit normal, see Figure 2 and also [4,30].Then, given a critical angle threshold ψ (expressed now in radians), the overhang angle constraint in the current context of shape and topology optimization would demand that an optimal geometry for the object Ω 1 has minimal regions where the overhang angles there do not lie in the interval [ψ, 2π−ψ] (i.e., Ω 1 should be self-supporting as much as possible), in addition to other mechanical considerations such as minimal compliance.For structural topology optimization we employ the phase field methodology proposed in [28], later popularized by many authors to other applications such as multi-material structural topology optimization [22,89], compliance optimization [24,78], topology optimization with local stress constraints [29], nonlinear elasticity [73], elastoplasticity [6], eigenfrequency maximization [45,78], graded-material design [32], shape optimization in fluid flow [42,43,44], as well as resolution strategies for some inverse identification problems [21,57].The key idea is to cast the structural topology optimization problem as a constrained minimization problem for a phase field variable ϕ, where the geometry of an optimal design Ω 1 for the object can be realized as a certain level-set of ϕ.In the language of optimal control theory, ϕ acts as a control and influences a response function, e.g., the elastic displacement u of the object, that is the solution to a system of partial differential equations (see (2.3) below).With an appropriate objective functional, for instance, a weighted sum of the mean compliance and a phase field formulation of the overhang constraint (see (2.7) below), we analyze the PDE-constrained minimization problem for an optimal design variable ϕ.
In the above references, previous authors employ an isotropic Ginzburg-Landau functional that has the effect of perimeter penalization (see Section 3 for more details), and in the current context this means no particular directions are preferred/discouraged in the optimization process.In this work we follow the ideas first proposed in [4,34] that use an anisotropic perimeter functional as a geometric constraint for overhangs (see also similar ideas in [2] for support structures), and introduce a suitable anisotropic Ginzburg-Landau functional for the phase field optimization problem to enforce the overhang angle constraint.In one sense, our first novelty is that we generalize earlier works in phase field structural topology optimization by considering anisotropic energies.The precise formulation and details of the problem are given in the next sections.For more details regarding the use of anisotropy in phase field models we refer the reader to [15,16,46,81].
Our second novelty is the analysis of the anisotropic phase field optimization problem.We establish the existence of minimizers (i.e., optimal designs) and derive first-order necessary optimality conditions.It turns out that our proposed approach to the overhang angle constraint requires an anisotropy function γ that is not continuously differentiable.In turn, the derivation of the associated necessary optimality conditions becomes nonstandard.We overcome this difficulty with subdifferential calculus, and provide a characterization result for the subdifferential of convex functionals whose arguments are weak gradients.Then, we perform a formally matched asymptotic analysis for the optimization problem with a differentiable γ to infer the corresponding sharp interface limit.
Lastly, let us provide a non-exhaustive overview on related works that integrate the overhang angle constraint within a topology optimization procedure.Our closest counterpart is the work [4] which employs an anisotropic perimeter functional in a shape optimization framework implemented numerically with the level-set method.The authors observed that oscillations along the object boundary would develop even though the design complied with the angle constraint.This so-called dripping effect is attributed to instability effects of the anisotropic perimeter functional, brought about by the non-convexity of the anisotropy used in [4], see also Example 3.3 and Figure 10 below and the discussion in [52,Sec. 7.3].An alternative mechanical constraint based on modeling the layer-by-layer construction process is then proposed to provide a different treatment of the overhangs and seems to suppress the dripping effect (see also [10] for similar ideas), but is computationally much more demanding.These boundary oscillations have also been observed earlier in [75], which used a Heaviside projection based integral to encode the overhang angle constraint within a density based topology optimization framework, and can be suppressed by means of an adaptive anisotropic filter introduced in [66].Similar projection techniques for density filtering are used in [49], which was combined with an edge detection algorithm to evaluate feasible and non-feasible contours during optimization iterations.In [50], control of overhang angles is achieved by means of filtering with a wedge-shaped support, while in [92] an explicit quadratic function with respect to the design density is used to formulate the self-support constraint.A different approach was proposed in [58,59] using a spatial filtering which checks element densities row-by-row and remove all elements not supported by the previous row.This idea is then extended in [60,72,83] to construct new filtering schemes that include other relevant manufacturing constraints.Another method was proposed in [51] using the frameworks of moving morphable components and moving morphable voids to provide a more explicit and geometric treatment of the problem, and is capable of simultaneously optimizing structural topology and build orientation.Finally, let us mention [85,86], where a filter is developed based on front propagation in unstructured meshes that has a flexibility in enforcing arbitrary overhang angles.Focusing only on overhang angles in enclosed voids, [65] combined a nonlinear virtual temperature method for the identification of enclosed voids with a logarithmic-type function to constraint the area of overhang regions to zero.
The rest of this paper is organized as follows: in Section 2 we formulate the phase field structural optimization problem to be studied, and in Section 3 we describe our main idea of introducing anisotropy, along with some examples and relevant choices, as well as a useful characterization of the subdifferential of the anisotropic functional.The analysis of the structural optimization problem is carried out in Section 4, where we establish analytical results concerning minimizers and optimality conditions.The connection between our work and that of [4] is explored in Section 5 where we look into the sharp interface limit, and, finally, in Section 6 we present the numerical discretization and several simulations of our approach.

Problem formulation
In a bounded domain Ω ⊂ R d with Lipschitz boundary Γ := ∂Ω that exhibits a decomposition Γ = Γ D ∪ Γ g ∪ Γ 0 , we consider a linear elastic material that does not fully occupy Ω.We describe the material location with the help of a phase field variable ϕ : Ω → [−1, 1].In the phase field methodology, we use the level set {ϕ = −1} = {x ∈ Ω : ϕ(x) = −1} to denote the region occupied by the elastic material, and {ϕ = 1} to denote the complementary void region.These two regions are separated by a diffuse interface layer {|ϕ| < 1} whose thickness is proportional to a small parameter ε > 0. Since ϕ describes the material distribution within Ω, complete knowledge of ϕ allows us to determine the shape and topology of the elastic material.
Notation.For a Banach space X, we denote its topological dual by X * and the corresponding duality pairing by •, • X .For any p ∈ [1, ∞] and k > 0, the standard Lebesgue and Sobolev spaces over Ω are denoted by L p := L p (Ω) and W k,p := W k,p (Ω) with the corresponding norms • L p and • W k,p .In the special case p = 2, these become Hilbert spaces and we employ the notation H k := H k (Ω) = W k,2 (Ω) with the corresponding norm • H k .For convenience, the norm and inner product of L 2 (Ω) are simply denoted by • and (•, •), respectively.For our subsequent analysis, we introduce the space Let us remark that, despite we employ bold symbols to denote vectors, matrices, and vector-or matrix-valued functions, we do not introduce a special notation to indicate the corresponding Lebesgue and Sobolev spaces.Thus, when those terms occur in the estimates, the corresponding norm is to be intended in its natural setting.

State system
To obtain a mathematical formulation that can be further analyzed, we employ the ersatz material approach, see [5], and model the complementary void region as a very soft elastic material.This allows us to consider a notion of elastic displacement on the entirety of Ω, leading to the displacement vector u : Ω → R d and the associated linearized strain tensor Let C 0 and C 1 be the fourth order elasticity tensors corresponding to the ersatz soft material and the linear elastic material, respectively, that satisfy the standard symmetric conditions and assume there exist positive constants θ and Λ such that for any non-zero symmetric matrices A, B ∈ R d×d : where A : B = d i,j=1 A ij B ij denotes the Frobenius inner product between two matrices.With the help of the phase field variable ϕ, we define an interpolation fourth order elasticity tensor C(ϕ) as where g : R → [−1, 1] is a monotone function satisfying g(−1) = −1 and g(1) = 1.Then, it is clear that in the region {ϕ = −1} (resp.{ϕ = 1}) we obtain the elasticity tensor C 1 (resp.C 0 ) by substituting ϕ = −1 (resp.ϕ = 1) into the definition of C(ϕ).As an example, we can consider where 0 < ε 1 is a constant and for any symmetric matrix A ∈ R d×d , with identity matrix I and Lamé constants λ 1 and µ 1 to obtain a linear interpolating elasticity tensor.Another example for C(ϕ) uses a quadratic interpolation function Let f : Ω → R d denote a body force and g : Γ g → R d denote a surface traction force.On Γ D ⊂ Γ we assign a zero Dirichlet boundary condition for the displacement u and on Γ 0 we assign a traction-free boundary condition.This leads to the following elasticity system for the displacement u, representing the state system of the problem: where n denotes the outward unit normal to Γ = Γ D ∪ Γ g ∪ Γ 0 , and h : is a function introduced so that the body force f only acts on the elastic material {ϕ = −1} (where h(ϕ) = 1) and not on the ersatz material {ϕ = 1} (where h(ϕ) = 0).We extend h(ϕ) as 1 if ϕ < −1 and as 0 if ϕ > 1.
The well-posedness of (2.3) is a simple consequence of [22,Thms. 3.1,3.2],which we summarize as follows: where H d−1 denotes the (d − 1)-dimensional Hausdorff measure.Moreover, there exists a positive constant C, independent of ϕ, such that In addition, let M > 0 and , and u i be the associated solution to (2.4).Then, there also exists a positive constant C, depending on the data of the system and M , but independent of the difference ϕ 1 − ϕ 2 , such that The above result provides a notion of a solution operator, also referred to as the control-to-state operator, S : ϕ → u where u ∈ H 1 D (Ω, R d ) is the unique solution to (2.4) corresponding to ϕ ∈ L ∞ (Ω).As such, we may view the phase field variable ϕ as a design variable that encodes the elastic response of the associated material distribution through the operator S and seek an optimal material distribution ϕ that fulfills suitable constraints and minimizes some cost functional.

Cost functional and design space
We define the design space, i.e., the set of admissible design variables, as (2.6) where m ∈ (−1, 1) is a fixed constant for the mass constraint, and motivated by the context of additive manufacturing, we propose the following cost functional to be minimized: In (2.7), the second term premultiplied by β > 0 is the mean compliance functional, while the first term premultiplied by α > 0 is an anisotropic Ginzburg-Landau functional with anisotropy function γ.We delay the detailed discussion on γ to the next section and remark here that for the isotropic case γ(x) = |x|, x ∈ R d , it is well-known that the Ginzburg-Landau functional is an approximation of the perimeter functional.Therefore, (2.7) can be viewed as a weighted sum between (anisotropic) perimeter penalization and mean compliance.
Lastly, for the non-negative potential function Ψ in (2.7), we require that it has ±1 as its global minima.While many choices are available, in light of the design space V m , we consider Ψ as the double obstacle potential It is worth pointing out that it holds Ψ(ϕ) = Ψ 0 (ϕ) for ϕ ∈ V m .Then, the structural optimization problem we study can be expressed as the following: Minimize J(ϕ, u) subject to ϕ ∈ V m and u solving (2.3). (P) 3 Anisotropic Ginzburg-Landau functional In the phase field methodology, the (isotropic) Ginzburg-Landau functional reads as where 0 < ε 1 and Ψ is a non-negative double-well potential that has ±1 as its minima, i.e., {s ∈ R : Ψ(s) = 0} = {±1}.Heuristically, the minimization of (3.1) in H 1 (Ω) results in minimizers that take near constant values close to ±1 in large regions of the domain Ω, which are separated by thin interfacial regions with thickness scaling with ε over which the functions transit smoothly from one value to the other.Formally in the limit ε → 0 these minimizers converge to functions that only take values in {±1}.This is made rigorous in the framework of Γ-convergence by the seminal work of Modica and Mortola [68,69].
To facilitate the forthcoming discussion, we review some basic properties for functions of bounded variations.For a more detailed introduction we refer the reader to [8,36].A function u ∈ L 1 (Ω) is a function of bounded variation in Ω if its distributional gradient Du is a finite Radon measure.The space of all such functions is denoted as BV(Ω) and its endowed with the norm , where for u ∈ BV(Ω), its total variation TV(u) is defined as The space BV(Ω, {a, b}) denotes the space of all BV(Ω) functions taking values in {a, b}.We say that a set U ⊂ Ω is a set of finite perimeter, or a Caccioppoli set, if its characteristic function χ U , where ).The perimeter of a set of finite perimeter U in Ω is defined as while its reduced boundary ∂ * U is the set of all points y ∈ R d such that |Dχ U |(B r (y)) > 0 for all r > 0, with B r (y) denoting the ball of radius r centered at y, and exists and |ν U (y)| = 1.
The unit vector ν U (y) is called the measure theoretical unit inner normal to U at y, a theorem by De Giorgi yields the connection P Ω (U ) = H d−1 (∂ * U ), see, e.g., [8].Then, the result of Modica and Mortola [68,69] can be expressed as follows: The Γ-limit of the extended functional Of particular interest is the following property of Γ-convergence, which states that if (i) E 0 is the Γ-limit of E ε , (ii) u ε is a minimizer to E ε for every ε > 0, (iii) {u ε } ε>0 is a precompact sequence, then every limit of a subsequence of {u ε } ε>0 is a minimizer for E 0 .This provides a methodology to construct minimizers of E 0 as limits of minimizers to E ε , provided the associated Γ-limit is precisely E 0 .
Returning to our discussion and problem in additive manufacturing, in [4] it was proposed to use anisotropic perimeter functionals to model the overhang angle constraint.In our notation, for a set U of bounded variation with ν U as the measure theoretical inward unit normal, these functionals take the form with a C 1 function γ : R d → R. Two choices were suggested in [4]: where ψ is a fixed angle threshold, e d denotes the build direction, and ψ i : R d → R, for i = 1, . . ., m, are given pattern functions with ν ψ i := ∇ψ i /|∇ψ i |.The first choice γ a penalizes the regions of the boundary ∂ * U where the angle between the outward normal (−ν U ) and the negative build direction (−e d ) is smaller than ψ, while the second choice γ b compels the unit normal ν U to be close to at least one of the directions ν ψ i .
A phase field approximation of the anisotropic perimeter functional (3.2) is the following anisotropic Ginzburg-Landau functional where as before, Ψ is a double well potential with ±1 as its minima.Then, for a convex γ : R d → R that is positively homogeneous of degree one (see (A1)), one has the analogue of the result by Modica and Mortola for anisotropic energies (see [19,20,27,71], and also Lemma 5.1 below), that is where the extended functionals E γ,ε and E γ,0 are defined as This motivates our consideration of the objective functional (2.7) and of the study of the related minimization problem for the overhang angle constraint.Furthermore, let us formally state the corresponding sharp interface limit (ε → 0) of the structural optimization problem as: Note that by Lemma 2.1, the solution operator S : ϕ → u is well-defined for ϕ ∈ BV(Ω, {−1, 1}).The connection between (P) and (P 0 ) will be explored in Section 5.

Anisotropy function, Wulff shape and Frank diagram
Consider an anisotropic density function γ : R d → R satisfying (A1) γ is positively homogeneous of degree one: which immediately implies γ(0) = 0.
(A2) γ is positive for non-zero vectors: (A3) γ is convex: Note that it is sufficient to assign values of γ on the unit sphere ∂B 1 (0) in R d , since by the one-homogeneity property (A1) we can define for any p = 0 with p = p/|p| γ(p) := |p|γ p .
For such anisotropy density functions and smooth hypersurfaces Γ with normal vector field ν, we define the anisotropic interfacial energy as Then, the isoperimetric problem involves finding a hypersurface Γ * that minimizes F γ under a volume constraint.This problem has been well-studied by many authors (see for instance [39,40,79,80] and [52] and the references therein) and the solution is given as the boundary of the region called the Wulff shape [90] where γ * is the dual function of γ defined as The dual function γ * also satisfies (A1)-(A3) and hence we can view the Wulff shape W as the 1-ball of γ * .Besides the Wulff shape, another region of interest that is used to visualize the effects of the anisotropy is the Frank diagram [38], which is defined as the 1-ball of γ: which, due to (A3), is always a convex subset of R d .To see how the shape of the boundary of F determines which directions of unit sphere in R d is preferred by the anisotropic density function γ, let us consider three examples.
Example 3.1 (Isotropic case).Consider γ(q) = |q| for q ∈ R d .Then, the associated Frank diagram is just the unit ball in R d , with boundary {γ(q) = |q| = 1}.As all points on the boundary are equidistant to the origin, all directions of the unit sphere in R d are equally preferable.
Example 3.2 (Convex example).Consider the Frank diagram shown in the left of Figure 3, whose boundary is composed of a circular arc C and a horizontal line L. The black dot denotes the origin in R 2 .For any unit vector in R 2 , we denote by φ ∈ [0, 2π) the angle it makes with the negative y-axis measured anticlockwise (see also Figure 2).Then, there exists θ > 0 such that all unit vectors with angle in [0, θ] ∪ [2π − θ, 2π) are associated with the horizontal line L in Figure 3, while all unit vectors with angle in (θ, 2π − θ) are associated with the circular arc C. Notice, if the origin lies on L, then θ = π 2 , and C is the upper semicircle.
Let p ∈ C and q ∈ L be arbitrary, and set p = p |p| , q = q |q| as their unit vectors.It is clear from the figure that |p| ≥ |q|, and since γ(p) = γ(q) = 1, by (A1) we see that This implies that γ( p) ≤ γ( q), and from the viewpoint of minimizing the interfacial energy F γ in (3.5), directions p are preferable to directions q.Consequently, from the Frank diagram in Figure 3   With such γ, consider two spatial points x 1 and x 2 at the same height, see Figure 4. Connecting them via a horizontal straight line is energetically expensive since this is associated to a direction with angle zero (where on the boundary of the Frank diagram is closest to the origin).An energetically more favorable connection is a zigzag path from x 1 and x 2 whose normal vectors oscillate between r 1 and r 2 .This is similar to a behavior termed "dripping effect" in [4] (see also [75,Fig. 15]), which is the tendency for shapes to develop oscillatory boundaries in order to meet the overhang angle constraints.

Relevant examples of anisotropic density function
Motivated by the above discussion, in this section we provide some examples of γ that achieve the Frank diagrams shown in Figure 3.

An example of a convex anisotropy
To fix the ideas, let us begin with the two-dimensional case.For a fixed constant α ∈ (0, 1), we consider the function where the set V ⊆ R 2 will be determined in the following.To achieve the boundary of the Frank diagram F shown in the left of Figure 3, we notice that ∂F ∩ V is a circular arc of radius 1 centered at the origin, while ∂F ∩ V c implies x 2 = −α which is a horizontal line segment at height −α.Following the convention in Figure 2 where angles are measured anticlockwise from the negative x 2 -axis, we can parameterize the circular arc by (cos φ, sin φ) for φ ∈ [θ, 2π − θ] where θ := cos −1 (α), see Figure 3. Hence, the set V in the definition (3.9) can be characterized as Remark 3.2.Note that in the limiting case α = 1, the above set V is the entire plane R 2 .Hence we have the isotropic case γ(x) = |x| for all x ∈ R 2 , and ∂F is simply the unit circle.
For a parametric characterization of the set V , let c := α √ 1−α 2 and consider the two straight lines {x 2 = cx 1 } and {x 2 = −cx 1 } dividing R 2 into eight regions (see Figure 5), which we label as Region 1, 2, . . ., 8 in an anticlockwise direction starting from the positive x 1 -axis.Then, the horizontal straight line portion ∂F ∩ V c of the Frank diagram at height x 2 = −α is contained in Regions 6 and 7, whose union is described by the set {x 2 < −c|x 1 |}, while the circular arc portion ∂F ∩V is contained in Regions 1, . . ., 5 and 8, whose union is described by the set {x 2 ≥ −c|x 1 |}.Hence, a parameteric characterization of the set V in (3.9) is Generalizing to the d-dimensional case, we obtain From (3.10) we see that γ α ∈ C 0 (R d ) but it is not continuously differentiable at the points x d = −α|x|.Furthermore, it is clear that γ α satisfies the assumptions (A1)-(A3) and hence (A4).

An example of a non-convex anisotropy
For completeness, we provide an example of an anisotropic function γ that yields a nonconvex Frank diagram as seen in the right of Figure 3. Referring to the construction of the previous example, we present the two-dimensional case first.Fix λ ∈ (0, 1) and let (0, −λα) denote the intersection of the two straight line segments (which we call L 1 and L 2 respectively) just below the origin.Denoting by p − the endpoint of the left line . Consider, for some constant b > 0 to be identified, This corresponds to Region 6 in the right of Figure 5 that contains the line segment L 1 , which can be parameterized as Notice that γ satisfies (A1)-(A2) (due to the modulus), and a short calculation shows that for Hence, choosing b = (λα Similarly, a choice of tangent vector for L 2 is τ = ( √ 1 − α 2 , (λ − 1)α) , so that a normal vector for x 1 } which corresponds to Region 7 in Figure 5 that contains the line segment L 2 .Then, a short calculation shows that γ(x) = 1 for x ∈ L 2 .Thus, an example of an anisotropic function γ that give rise to a Frank diagram whose boundary is the right figure in Figure 3 is for α ∈ (0, 1), λ ∈ (0, 1) and x = (x 1 , x 2 ) ∈ R 2 .Notice that in the limit λ → 1, we recover the convex anisotropic function γ α defined in (3.10).
To generalize to the d-dimensional case, we notice that the lines L 1 and L 2 in the above discussion are now replaced by the lateral surface S of a cone with apex (0, 0, . . ., 0, −λα) ∈ R d , which can be parameterized as Then, by similar arguments leading to (3.11), we obtain the function , where we can verify that γ α,λ (x) = 1 for x ∈ S. In Figure 6 we display the Frank diagrams for non-convex anisotropy functions of the form (3.11) with λ = 0.5 and α= 0.7, 0.5, 0.3.

Subdifferential characterization
The analysis of the structural optimization problem (P) follows almost analogously as in [22].The major difference is the anisotropic Ginzburg-Landau functional in (2.7).From the examples of γ discussed in the previous subsection, our interest lies in anisotropy density functions that are convex and continuous, but not necessarily C 1 (R d ), which necessitates a non-trivial modification to the analysis performed in [22].In this subsection we focus only on the gradient part of (2.7) and investigate its subdifferential in preparation for the first-order necessary optimality conditions for (P) (cf.Theorem 4.2).We define A : R d → [0, ∞) as and from (A1)-(A4), it readily follows that A is convex, continuous with A(0) = 0, positive for non-zero vectors, positively homogeneous of degree two: and there exist positive constants c and C such that Associated to such a function A, we consider the integral functional is a maximal monotone operator from H 1 (Ω) to H 1 (Ω) * (see [13,Thm. 2.43]), which is equivalent to the property that, for any f ∈ H 1 (Ω) * , there exists at least one solution ϕ 0 ∈ D(∂F) with ξ 0 ∈ ∂F(ϕ 0 ) such that Let us now provide a useful lemma that characterizes the subdifferential ∂F(ϕ) in terms of elements of ∂A(∇ϕ) (see also [13, p. 146, Problem 2.7]).Heuristically, elements of ∂F(ϕ) are the negative weak divergences of elements of ∂A(∇ϕ).

Analysis of the structural optimization problem
By Lemma 2.1 the solution operator S : L ∞ (Ω) → H 1 D (Ω, R d ), S(ϕ) = u where u is the unique solution to (2.3) corresponding to the design variable ϕ, is well-defined.This allows us to consider the reduced functional J : in the structural optimization problem (P).Invoking [35, Thm. 1, §8.2.2], the convexity of γ and hence of A(•) = 1  2 |γ(•)| 2 yields that the gradient term in (2.7) is weakly lower semicontinuous in H 1 (Ω).Then, following the proof of [22,Thm. 4.1] we obtain the existence of an optimal design to (P).
Proof.Since the proof is now rather standard, we only sketch some of the essential details.Using the bound (2.5) we infer that the reduced functional J is bounded from below in V m .This allows us to consider a minimizing sequence {φ n } n∈N ⊂ V m such that J (φ n ) → inf ζ∈Vm J (ζ) as n → ∞.Then, again by (2.5) and also (A4) we deduce that {φ n } n∈N is uniformly bounded in H 1 (Ω) ∩ L ∞ (Ω), from which we extract a non-relabeled subsequence converging weakly to ϕ in H 1 (Ω).As V m is a convex and closed set, hence weakly sequentially closed, we infer also ϕ ∈ V m and by passing to the limit infimum of J (φ n ), employing the weak lower semicontinuity of the gradient term in (2.7) and also the weak convergence S(φ n ) → S(ϕ) in H 1 (Ω, R d ) we infer that J (ϕ) = inf ζ∈Vm J (ζ), which implies that ϕ is a solution to (P).
To derive the first-order optimality conditions, we first take note of the following result concerning the Fréchet differentiability of the solution operator S, which can be inferred from [22,Thm. 3.3]: and Then, we introduce the adjoint system for the adjoint variable p associated to ϕ, whose structure is similar to the state equation (2.3) for u: Via a similar argument (see also [22,Thm. 4.3]), the well-posedness of the adjoint system is straightforward: ) is the unique solution to (2.4), then p = βu.
Our main result is the following first-order optimality conditions for the structural optimization problem with anisotropy (P): ) and ϕ ∈ V m be an optimal design variable with associated state u = S(ϕ).Then, there exists ξ ∈ ∂A(∇ϕ) almost everywhere in Ω such that for all φ ∈ V m .
Proof.Recalling the definition of the functional F in (3.14), we introduce so that the reduced functional J can be expressed as J (φ) = αεF(φ) + K(φ).Using the differentiability of S we infer that K is also Fréchet differentiable with derivative at 1] a.e. in Ω} given as where w = S (ϕ) [ζ] is unique solution to the linearized system (4.1) and u = S(ϕ).We can simplify this using the adjoint system by testing (4.3) with v = w and testing (4.1) with η = p to obtain the identity so that, together with the relation p = βu, (4.6) becomes On the other hand, the convexity of F and the optimality of ϕ ∈ V m imply that holds for all t ∈ [0, 1] and arbitrary φ ∈ V m .Dividing by t and passing to the limit t → 0 yields The above inequality allows us to interpret ϕ ∈ V m as a solution to the convex minimization problem min where denotes the indicator function of the set V m .Using the well-known sum rule for subdifferentials of convex functions, see, e.g., [13,Cor. 2.63], the inequality (4.8) can be interpreted as 0 where ∂ denotes the subdifferential mapping in H 1 (Ω).This implies the existence of elements v ∈ ∂F(ϕ) and ψ ∈ ∂I Vm (ϕ) such that 0 = αεv + K (ϕ) + ψ. (4.10) For arbitrary φ ∈ V m , by definition of ∂I Vm (ϕ) we have Then, using Lemma 3.1, there exists ξ ∈ L 2 (Ω, R d ) with ξ ∈ ∂A(∇ϕ) almost everywhere in Ω and ( αεξ, The optimality condition (4.4) is then a consequence of the above and (4. for all ζ ∈ H 1 (Ω).In the above, the subdifferential of the indicator function I [−1,1] has the explicit characterization Indeed, instead of (4.9) we can interpret ϕ ∈ V m as a solution to the minimization problem this yields (4.11) with Lagrange multiplier µ := −|Ω| −1 (K (ϕ) [1] + θ).

Sharp interface asymptotics
Our interest is to study the behavior of solutions under the sharp interface limit ε → 0, which connects our phase field approach with the shape optimization approach of [4].

Γ-convergence of the anisotropic Ginzburg-Landau functional
We begin with the Γ-convergence of the extended anisotropic Ginzburg-Landau functional (3.4a) to the extended anisotropic perimeter functional (3.4b), which is formulated as follows: Lemma Then, there exists a sequence {u ε } ε>0 of Lipschitz continuous functions on Ω such that Then, there exists a non-relabeled subsequence ε → 0 and a limit function u such that The first and second assertions are known as the liminf and limsup inequalities, respectively, while the third assertion is the compactness property.In the following, we outline how to adapt the proofs of [20,Thms. 3.1 and 3.4] for the Γ-convergence result in the multi-phase case.Our present setting corresponds to the case N = 2 in their notation.

Formally matched asymptotic analysis
In this section we consider an anisotropy function γ ∈ C 2 (R d )∩W 1,∞ (R d ) satisfying (A1)-(A3) in the structural optimization problem (P), as well as a more regular body force f ∈ H 1 (Ω, R d ).The differentiability of γ implies the characterization of the subdifferential ∂A(∇ϕ) as the singleton set {[γDγ](∇ϕ)}, where Then, the corresponding optimality condition for a minimizer (ϕ ε , u ε ) to (P) becomes and its strong formulation reads as (see (4.11)) [γDγ](∇ϕ ε ) with Lagrange multiplier µ ε for the mass constraint.Our aim is to perform a formally matched asymptotic analysis as ε → 0, similar as in [25,Sec. 5], in order to infer the sharp interface limit of (5.5).The method proceeds as follows, see [31]: formally we assume that the domain Ω admits a decomposition where Ω I ε is an annular domain and ϕ ε ∈ V m satisfies The zero level set Λ ε := {x ∈ Ω : ϕ ε (x) = 0} is assumed to converge to a smooth hypersurface Λ as ε → 0, and that (ϕ ε , u ε ) admit an outer expansion in regions in Ω ± ε as well as an inner expansion in regions in Ω I ε .We substitute these expansions (in powers of ε) in (2.3) and (5.5), collecting terms of the same order of ε, and with the help of suitable matching conditions connecting these two expansions, we deduce an equation posed on Λ.For an introduction and more detailed discussion of this methodology we refer to [37,47].
To start, let us collect some useful relations.As γ is positively homogeneous of degree one, taking the relation γ(tp) = tγ(p) for t > 0 and p ∈ R d \ {0} and differentiating with respect to t, and also with respect to p leads to with the latter relation showing that Dγ is positively homogeneous of degree zero.Then, it is easy to see γDγ is positively homogeneous of degree one.Next, from (3.12), we see that A is positively homogeneous of degree two.Taking the relation A(tp) = t 2 A(p) and differentiating with respect to t twice leads to (5.9)

Outer expansions
For points x in Ω ± ε , we assume an outer expansion of the form where all functions are sufficiently smooth and the summations converge.From (5.6) we deduce that ϕ 0 (x) ∈ {−1, 1}, ϕ i (x) = 0 for i ≥ 1.

Inner expansions
We assume the outer boundary Γ + ε and inner boundary Γ − ε of the annular region Ω I ε can be parameterized over the smooth hypersurface Λ.Let ν denote the unit normal of Λ pointing from Ω − to Ω + , and we choose a spatial parameter domain U ⊂ R d−1 with a local parameterization r : U → R d of Λ.Let d denote the signed distance function of Λ with d(x) > 0 if x ∈ Ω + , and denote by z = d ε the rescaled signed distance.Then, by the smoothness of Λ, there exists δ 0 > 0 such that for all x ∈ {|d(x)| < δ 0 }, we have the representation where s = (s 1 , . . ., s d−1 ) ∈ U .In particular, r(s) ∈ Λ is the projection of x to Λ along the normal direction.This representation allows us to infer the following expansion for gradients and divergences [47]: for scalar functions b(x) = b(s(x), z(x)) and vector functions j(x) = j(s(x), z(x)).In the above ∇ Λ is the surface gradient on Λ and div Λ is the surface divergence.Analogously, for a vector function j, we find that For points close by Λ in Ω I ε , we assume an expansion of the form

Matching conditions
In a tubular neighborhood of Λ, we assume the outer expansions and the inner expansions hold simultaneously.Since Γ ± ε are assumed to be graphs over Λ, we introduce the functions Furthermore, we assume an expansion of the form ε converge to Λ, in the computations below we will deduce the values for Y ± 0 (s).Then, by comparing these two expansions in this intermediate region we infer matching conditions relating the outer expansions to the inner expansions via boundary conditions for the outer expansions.For a scalar function b(x) admitting an outer expansion ∞ k=0 ε k b k (x) and an inner expansion ∞ k=0 ε k B k (s, z), it holds that (see [31] and [47,Appendix D])

Analysis of the inner expansions
We introduce the notation B sym = 1 2 (B + B ) for a second order tensor B. Plugging in the inner expansions to the state equation (2.3), to leading order O( 1 ε 2 ) we find that Taking the product with U 0 and by the symmetry of C we have the relation Integrating over z and by parts, using ∂ z ν = 0 and by the matching condition lim z→Y ± 0 Coercivity of C yields that ∂ z U 0 ⊗ ν sym = 0 and hence ∂ z U 0 (s, z) = 0.This implies U 0 is constant in z and by the matching conditions we obtain (5.13) Then, using ∂ z U 0 = 0, to the next order O( 1 ε ) we obtain from the state equation (2.3) Integrating over z and using the matching condition we obtain From the optimality condition (5.5a), to leading order O( 1 ε 2 ) we obtain the trivial equality 0 = 0 due to ∂ z U 0 = 0. Recalling (5.11) for µ −1 , to the next order O( 1 ε ) of (5.5a) we have where we again use ∂ z U 0 = 0 to see that no terms from the elastic part contribute at order O( 1 ε ), the expression (5.12) for the divergence in the inner variables, the relation Dγ(ν) • ν = γ(ν) from (5.7), the fact that γDγ is positively homogeneous of degree one, so that (γDγ) , as well as the fact that ν is independent of z.We now construct a solution to (5.16) fulfilling the conditions Φ 0 (s, z) → ±1 as z → Y ± 0 (s) as per the matching conditions, and also satisfies ∂ z Φ 0 (s, z) > 0. Let ψ(z) be the monotone function defined in (5.3) which satisfies the boundary value problem with ψ(0) = 0 and ψ (z) > 0 for all z ∈ R. We define Φ 0 (s, z) = ψ z γ(ν(s)) , and a short calculation shows that

The property ψ(± π
2 ) = ±1 and the matching conditions for Φ 0 provide the identifications Moreover, since ψ is monotone and γ(ν) is positive, we can infer that ∂ z Φ 0 (s, z) > 0. Now, multiplying (5.16) with ∂ z Φ 0 and integrating over z yields the so-called equipartition of energy after using Ψ 0 (−1) = 0 and lim z→Y − 0 ∂ z Φ 0 (s, z) = 0. Performing another integration over z leads to the relation . (5.17)To the next order O(1), we find that (5.18) The aim is to multiply the above with ∂ z Φ 0 and integrate over z, but first, let us derive some useful relations.Using the symmetry C ijkl = C jikl of the elasticity tensor and ∂ z U 0 = 0, we deduce that, with the shorthand Hence, together with ∂ z (C(Φ 0 )E U ν) = 0 from (5.14) and the matching conditions we obtain the relation By (5.8), (5.9), (5.16), (5.17) and the matching conditions of ∂ z Φ 0 , we have and upon multiplying (5.18) with ∂ z Φ 0 and integrate over z, employ the matching conditions leads to the following solvability condition for Φ 1 : (5.19) Note that the assumed higher regularity f ∈ H 1 (Ω, R d ) allows us to define f on Λ.Thus, the sharp interface limit ε → 0 consists of the system (5.10)posed in Ω ± , along with the boundary conditions (5.13), (5.15) and (5.19) on Λ.

BGN anisotropies
For the numerical approximation of an optimal design variable ϕ ∈ V m to the structural optimization problem (P), we first consider a discretization of the variational inequality (5.4) with differentiable γ.In general, the term γDγ is highly nonlinear, and to facilitate the subsequent discussion regarding its discretization we first consider a matrix-type reformulation of the anisotropy density function introduced by the works of the first and third authors, see, e.g., [14,15,16] for more details.Suppose for some L ∈ N, the anisotropy density function γ can be expressed as γ l (q), where γ l (q) := q • G l q for q ∈ R d , ( with symmetric positive definite matrices G l ∈ R d×d for l = 1, . . ., L. Then, a short calculation shows that and recalling that A( For p close to q, we now perform a linearization DA(q) ≈ B(p)q, where so that B(q)q = DA(q) for all q ∈ R d .
Furthermore, as {G l } L l=1 are symmetric positive definite, the matrix B(p) is also symmetric positive definite for all p ∈ R d .We term anisotropies γ that are of the form (6.1) as BGN anisotropies.Notice that for the definition of the linearized matrix B it is sufficient to have anisotropy densities γ ∈ C 0 (R d ).However, due to our examples (3.10) and (3.11) it will turn out that the corresponding G l are no longer constant matrices.Thus, inspired by the works [14,15,16] we extend their approach to the case where G l are piecewise constant and dependent on q, as shown in the following examples.In our numerical approximation to the optimality condition detailed in the next section we simply replace DA(q) = [γDγ](q) by B(q)q, where B(q) now involves piecewise constant matrices G l (q).Example 6.1.The convex anisotropy (3.10) can be expressed as an extended BGN anisotropy γ α (q) = q • G(q)q, where we denote by P 1 (o) the set of all affine linear functions on o.In addition we define the convex subsets where (•, •) denotes the L 2 -inner product on Ω, as well as We now introduce a finite element approximation of the structural optimization problem (P) and the optimality conditions described above.Given Here τ denotes the time step size, (•, •) h is the usual mass lumped L 2 -inner product on Ω, and A, B C = ´Ω CA : B dx for any fourth order tensor C and any matrices A and B. We implemented the scheme (6.5) with the help of the finite element toolbox ALBERTA, see [77].To increase computational efficiency, we employ adaptive meshes, which have a finer mesh size h f within the diffuse interfacial regions and a coarser mesh size h c away from them, see [11,18] for a more detailed description.Clearly, the system (6.5a)decouples, and so we first solve the linear system (6.5a) in order to obtain u n h , and then solve the variational inequality (6.5b) for ϕ n h .We employ the package LDL, see [33], together with the sparse matrix ordering AMD, see [9], in order to solve (6.5a).For the variational inequality (6.5b) we employ a secant method as described in [26] to satisfy the mass constraint (ϕ n h , 1) = m, and use a nonlinear multigrid method similar to [55] for solving the variational inequalities over V h that arise as the inner problems for the secant iterations.The second method always converged in at most five steps.Finally, to increase the efficiency of the numerical computations in this paper, at times we exploit the symmetry of the problem and perform the computations in question only on half of the desired domain Ω.In those cases we prescribe "free-slip" boundary conditions for the displacement field u n h on the symmetry plane Γ D i , that is, we replace (6.4) with All computations performed in this work are for spatial dimension d = 2.For the remainder of the paper we consider the quadratic interpolation function (2.2) in the elasticity tensor C(ϕ), forcing f = 0, objective functional weightings β = 1 and α = 1 unless further specified.

Optimality condition without elasticity
We first investigate the setting with g = 0, so that (6.5a) yields u n h = 0 and (6.5b) reduces to an Allen-Cahn variational inequality on V h m .From the discussion in Section 3, this is a phase field approximation of the volume preserving anisotropic mean curvature flow (due to the mass constraint), and it is expected that the long time behavior of solutions to display the Wulff shape (3.6) corresponding to the anisotropy γ, see Figure 7.In Figure 8 we display snapshots of the solution at times t = 0, 0.001, 0.005 and 0.03 with anisotropy (6.2) for parameter values α = 0.5, δ = 0.1 and ε = 1/(16π).We notice from the bottom plot of the discrete Ginzburg-Landau energy E h γ (ϕ n h ) = (εA(∇ϕ n h ) + 1 ε Ψ(ϕ n h ), 1) h that it is decreasing over time, and on the top right we observe the expected Wulff shape attained near equilibrium.

Dripping effect of a straight interface
Next, we investigate the role of the convexity/non-convexity of γ in the generation of the dripping effect [4,10,75].We solve (6.5) with g = 0 and initial condition taken as some large perturbation of a straight interface.Taking (6.2) as the anisotropy with α = 0.5, δ = 0.1 and ε = 1/(32π), in Figure 9 we display the snapshots of the solution at times t = 0, 0.001, 0.005 and 0.02, where it is clear that the oscillations are dampened over time.Thus, with a convex anisotropy the dripping effect is suppressed.In contrast, taking γ as the non-convex anisotropy (3.11) with λ = 0.5 now yields Figure 10, where the corresponding snapshots of the solution at times t = 0, 0.001, 0.005, 0.02 are displayed.
Here we clearly observe the behavior described in Example 3.3.

Cantilever beam computation
Inspired by [4, Figure 8], for the state system (2.3), we consider the setting Ω = (− 1 2 , 1 2 ) × (− 1 2 , 1 2 ), Γ D = (− 1 4 , 1 4 ) × {− 1 2 }, Γ g = (−0.02,0.02) × { 1 2 }, g = g 0 , (6.7)In Figure 12 we provide a more detailed comparison of the cantilever beam designs for the isotropic case α = 1 and the strongly anisotropic cases α= 0.3, 0.2, 0.1 with loading magnitudes g = 50 and g = 30.In both situations, the anisotropic designs exhibit less interior structures and connecting bridges with steeper slopes than the isotropic designs, which is to be expected from the construction of γ α .It is worth mentioning that some (but not all) of the interfacial regions (colored deep blue) are thicker for smaller values of α.This is attributed to the fact that, from the expression of Φ 0 obtained in the formally matched asymptotic analysis in Section 5.2, the interfacial thickness at a spatial point x = (s, z) in the transformed coordinate system is proportional to γ(ν(s)) with unit normal ν.Recalling the left of Figure 3, the normal vectors with directions associated to the line segment L attain higher values of γ, and thus a thicker interfacial region, compared to directions associated with the circular arc C. Smaller values of α amplify the thickness as the line segment L is closer to the origin, leading to higher values of γ.
In Figure 13 we display refined designs taking as initial conditions the final states in Figure 12 as well as a smaller value of ε = 1/(64π).The overall designs are similar to Figure 12 with some minor changes, particularly for the isotropic case α = 1 with g = 50.We also note that the thicknesses of the interfacial layers in Figure 13 are smaller than those in Figure 12, particularly for α = 0.1, which is due to the smaller values of ε used.

Figure 1 :
Figure 1: Definition of overhang angle measured from the base plate to the tangent line on the boundary.

Figure 2 :
Figure 2: The definition of overhang angle used in this work.(Left) A structure with four angles measured from the negative build direction −e d to the outward unit normal on the boundary.(Right) Visualization of the angles on the unit circle, where the angle φ i is associated to unit normal n i .

Figure 3 :
Figure 3: Illustrations of the Frank diagram for Examples 3.2 and 3.3.(Left) Convex case and (Right) non-convex case.

Example 3 . 3 (
Non-convex example).Consider the Frank diagram shown in the right of Figure 3, whose boundary encloses a non-convex set.Let r 1 and r 2 denote the two unit vectors whose angles, say θ and 2π − θ, respectively, associate to the two endpoints of the circular arc.From previous discussions, the anisotropy density function γ will prefer directions with angles in (θ, 2π − θ).

Figure 4 :
Figure 4: Connection between two spatial points with non-convex anisotropic density function γ whose Frank diagram looks like the right of Figure 3.