Overhang penalization in additive manufacturing via phase field structural topology optimization with anisotropic energies

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^{\circ}$ 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 $\Omega_{1}\subset\mathbb{R}^{d}$, $d=2,3,$ being built within a hold-all rectangular domain $\Omega=[0,1]^{d}$ in a layer-by-layer fashion with build direction $\bm{e}_{d}=(0,\dots,1)^{\top}$. The base plate is the region $\mathcal{B}:=[0,1]^{d-1}\times\{0\}$ where $\Omega_{1}$ is assembled on. For any $\bm{p}\in\partial\Omega_{1}\setminus\mathcal{B}$, the overhang angle $\alpha$ is defined as the angle from the base plate to the tangent line at $\bm{p}$ (denoted as $\alpha_{1}$ and $\alpha_{2}$ in Figure 1).

Then, there exists a critical threshold angle $\psi$ where the portion of $\Omega_{1}$ is self-supporting (resp. requires support structures) if the overhang angle there is greater (resp. smaller) than $\psi$. Conventional wisdom from practitioners puts the critical angle $\psi$ at $45^{\circ}$ [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 $\psi$ (e.g., $40^{\circ}$ 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 $\partial\Omega_{1}$. Denoting this as $\beta$, a simple calculation shows that $\beta=180^{\circ}-\alpha$. Then, regions where $\beta>180^{\circ}-\psi$ (typically $135^{\circ}$ if $\psi=45^{\circ}$) 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 $-\bm{e}_{d}$ to the outer unit normal, see Figure 2 and also [4, 30]. Then, given a critical angle threshold $\psi$ (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 $\Omega_{1}$ has minimal regions where the overhang angles there do not lie in the interval $[\psi,2\pi-\psi]$ (i.e., $\Omega_{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 $\varphi$, where the geometry of an optimal design $\Omega_{1}$ for the object can be realized as a certain level-set of $\varphi$. In the language of optimal control theory, $\varphi$ acts as a control and influences a response function, e.g., the elastic displacement $\bm{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 $\varphi$.

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 $\gamma$ that is not continuously differentiable. In turn, the derivation of the associated necessary optimality conditions becomes non-standard. 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 $\gamma$ 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.

In a bounded domain $\Omega\subset\mathbb{R}^{d}$ with Lipschitz boundary $\Gamma:=\partial\Omega$ that exhibits a decomposition $\Gamma=\Gamma_{D}\,\cup\,\Gamma_{g}\,\cup\,\Gamma_{0}$, we consider a linear elastic material that does not fully occupy $\Omega$. We describe the material location with the help of a phase field variable $\varphi:\Omega\to[-1,1]$. In the phase field methodology, we use the level set $\{\varphi=-1\}=\{{\bm{x}\in\Omega:\varphi(\bm{x})=-1}\}$ to denote the region occupied by the elastic material, and $\{\varphi=1\}$ to denote the complementary void region. These two regions are separated by a diffuse interface layer $\{|\varphi|<1\}$ whose thickness is proportional to a small parameter $\varepsilon>0$. Since $\varphi$ describes the material distribution within $\Omega$, complete knowledge of $\varphi$ allows us to determine the shape and topology of the elastic material.

In the phase field methodology, the (isotropic) Ginzburg–Landau functional reads as

where $0<\varepsilon\ll 1$ and $\Psi$ is a non-negative double-well potential that has $\pm 1$ as its minima, i.e., $\{s\in\mathbb{R}\,:\,\Psi(s)=0\}=\{\pm 1\}$. Heuristically, the minimization of (3.1) in $H^{1}(\Omega)$ results in minimizers that take near constant values close to $\pm 1$ in large regions of the domain $\Omega$, which are separated by thin interfacial regions with thickness scaling with $\varepsilon$ over which the functions transit smoothly from one value to the other. Formally in the limit $\varepsilon\to 0$ these minimizers converge to functions that only take values in $\{\pm 1\}$. This is made rigorous in the framework of $\Gamma$-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\in L^{1}(\Omega)$ is a function of bounded variation in $\Omega$ if its distributional gradient $\mathrm{D}u$ is a finite Radon measure. The space of all such functions is denoted as $\mathrm{BV}(\Omega)$ and its endowed with the norm $\|\cdot\|_{\mathrm{BV}(\Omega)}=\|\cdot\|_{L^{1}(\Omega)}+\mathrm{TV}(\cdot)$, where for $u\in\mathrm{BV}(\Omega)$, its total variation $\mathrm{TV}(u)$ is defined as

The space $\mathrm{BV}(\Omega,\{a,b\})$ denotes the space of all $\mathrm{BV}(\Omega)$ functions taking values in $\{a,b\}$. We say that a set $U\subset\Omega$ is a set of finite perimeter, or a Caccioppoli set, if its characteristic function $\chi_{U}$, where $\chi_{U}(\bm{x})=1$ if $\bm{x}\in U$ and $\chi_{U}(\bm{x})=0$ if $\bm{x}\notin U$, belongs to $\mathrm{BV}(\Omega,\{0,1\})$. The perimeter of a set of finite perimeter $U$ in $\Omega$ is defined as

while its reduced boundary $\partial^{*}U$ is the set of all points ${{\bf y}}\in\mathbb{R}^{d}$ such that $|\mathrm{D}\chi_{U}|(B_{r}({{\bf y}}))>0$ for all $r>0$, with $B_{r}({{\bf y}})$ denoting the ball of radius $r$ centered at ${{\bf y}}$, and

The unit vector $\bm{\nu}_{U}({{\bf y}})$ is called the measure theoretical unit inner normal to $U$ at ${{\bf y}}$, a theorem by De Giorgi yields the connection $\mathcal{P}_{\Omega}(U)=\mathcal{H}^{d-1}(\partial^{*}U)$, see, e.g., [8]. Then, the result of Modica and Mortola [68, 69] can be expressed as follows: The $\Gamma$-limit of the extended functional

is equal to

with the constant $c_{\Psi}:=\int_{-1}^{1}\sqrt{2\Psi(s)}ds$.

Of particular interest is the following property of $\Gamma$-convergence, which states that if (i) $\mathcal{E}_{0}$ is the $\Gamma$-limit of $\mathcal{E}_{\varepsilon}$, (ii) $u_{\varepsilon}$ is a minimizer to $\mathcal{E}_{\varepsilon}$ for every $\varepsilon>0$, (iii) $\{u_{\varepsilon}\}_{\varepsilon>0}$ is a precompact sequence, then every limit of a subsequence of $\{u_{\varepsilon}\}_{\varepsilon>0}$ is a minimizer for $\mathcal{E}_{0}$. This provides a methodology to construct minimizers of $\mathcal{E}_{0}$ as limits of minimizers to $\mathcal{E}_{\varepsilon}$, provided the associated $\Gamma$-limit is precisely $\mathcal{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 $\bm{\nu}_{U}$ as the measure theoretical inward unit normal, these functionals take the form

with a $C^{1}$ function $\gamma:\mathbb{R}^{d}\to\mathbb{R}$. Two choices were suggested in [4]:

where $\psi$ is a fixed angle threshold, $\bm{e}_{d}$ denotes the build direction, and $\psi_{i}:\mathbb{R}^{d}\to\mathbb{R}$, for $i=1,\dots,m$, are given pattern functions with $\bm{\nu}_{\psi_{i}}:=\nabla\psi_{i}/|\nabla\psi_{i}|$. The first choice $\gamma_{a}$ penalizes the regions of the boundary $\partial^{*}U$ where the angle between the outward normal $(-\bm{\nu}_{U})$ and the negative build direction $(-\bm{e}_{d})$ is smaller than $\psi$, while the second choice $\gamma_{b}$ compels the unit normal $\bm{\nu}_{U}$ to be close to at least one of the directions $\bm{\nu}_{\psi_{i}}$.

A phase field approximation of the anisotropic perimeter functional (3.2) is the following anisotropic Ginzburg–Landau functional

where as before, $\Psi$ is a double well potential with $\pm 1$ as its minima. Then, for a convex $\gamma:\mathbb{R}^{d}\to\mathbb{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 $\mathcal{E}_{\gamma,\varepsilon}$ and $\mathcal{E}_{\gamma,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 $(\varepsilon\to 0)$ of the structural optimization problem as:

where $\mathrm{BV}_{m}(\Omega,\{-1,1\})=\{f\in\mathrm{BV}(\Omega,\{-1,1\})\,:\,\int_{\Omega}f\,\mathrm{dx}=m|\Omega|\}$ and

Note that by Lemma 2.1, the solution operator $\bm{S}:\varphi\mapsto\bm{u}$ is well-defined for $\varphi\in\mathrm{BV}(\Omega,\{-1,1\})$. The connection between (P) and ($\mathbf{P}_{0}$) will be explored in Section 5.