This work is licensed under the Creative Commons Attribution-NonCommercial-NoDerivatives 3.0 License.
Introduction
The additive manufacturing technologies have recently demonstrated a unique potential in constructing structures with a high degree of complexity, thereby allowing to process almost directly the designs predicted by topology optimization algorithms. These breakthroughs however come along with new challenges. One of them is to overcome the difficulty of building shapes showing large overhangs, i.e. regions hanging over void without sufficient support from the lower structure.
To give a hint of why such regions are problematic from the point of view of additive manufacturing, a few generalities about these technologies are in order; see [10] or the introduction in [6] for further details. Additive manufacturing is a common label for a whole range of processes, which all begin by a decomposition (or slicing) of the produced shape into a series of two-dimensional layers; these are then constructed successively, one on top of the other; see Figure 1.
Fig. 1
Sketch of the slicing procedure, at the beginning of all additive manufacturing processes.
The construction of each individual layer may be achieved owing to quite different technologies; for instance:
Material extrusion methods, such as the well-known Fused deposition modeling (FDM) technique, proceed by extruding a molten filament from a nozzle, which is deposed under the form of rasters;
Powder bed fusion methods, such as Electron Beam Melting (EBM) or Selective Laster Melting (SLM), rely on a laser to bind the grains of a (metallic) powder together.
The physical origin of the problems caused by overhangs during the construction process depends on the particular technology: in the case of material extrusion methods, assembling a large overhang feature implies that the machine tool has to depose material over void; in the case of powder-bed fusion technologies, overhangs induce large residual stresses, which eventually entail warpage.
The prevailing remedy in the literature to the issue of overhangs consists in erecting a scaffold structure alongside with the assembly of the desired shape; see e.g. [8]. This auxiliary, sacrificial structure has to be removed at the end of the construction process, which is costly and cumbersome.
A different idea is to impose since the design optimization investigations that the shape should be self-supporting - that is, free of overhangs. In this direction, ad hoc criteria, based on a minimum angle between the structural boundary and the horizontal directions, have been used hitherto in order to tackle this issue [11, 13, 14, 119].
In the recent journal articles [2, 3], we have introduced new formulations to deal with overhang constraints in shape and topology optimization algorithms. The purpose of the present note is to give an introduction to this model; as such, it does not contain any original result. The reader is referred to the aforementionned articles for a complete mathematical exposition, and for the discussion of efficient numerical strategies dedicated to the treatment of overhang constraints.
The key idea of our work is to introduce a new constraint functional Psw(Ω) for shape optimization problems, which appraises the constructibility of the shape Ω at each stage of its layer by layer assembly; in particular, overhang constraints are naturally addressed by this formulation. To achieve this, in the setting of the optimization problem, we distinguish the mechanical situation where the final (completed) shape Ω is utilized, on which the optimization criterion is based, and that where Ω, and all the associated intermediate shapes Ωh, are under construction, which guides the definition of our constraint functional.
This note is organized as follows. In Section 2, we introduce the shape optimization problem at stake. In Section 3, we describe the intuitive attempt to rely on geometric functionals to deal with overhang constraints, and explain why it proves insufficient for this purpose. Then, in Section 4, we present our mechanical model of the context in which Ω is constructed, and we formulate our manufacturing constraint functional Psw(Ω) accordingly. Several details of its mathematical analysis are outlined, notably the calculation of its shape derivative, and an algorithm taking advantage of its intrinsic structure is proposed, which allows to accelerate significantly the (costly) calculations entailed by the evaluation of Psw(Ω) and its derivative. Eventually, a numerical validation and several examples are provided in Section 5.
Presentation of the shape optimization problem
A shape is a bounded, regular domain Ω ⊂ ℝd, d = 2, 3, filled with a linear elastic material with Hooke’s law A. In the context of its final utilization, Ω is clamped on a subset ΓD⊂∂Ω, and it is submitted to surface loads f ∈ L2(ΓN)d applied on a region ΓN of ∂Ω disjoint from ΓD. The remaining part $\Gamma := \partial \Omega \setminus \overline{\Gamma_D \cup \Gamma_N}$ is traction-free. The elastic displacement uΩ is the unique solution in $H^1_{\Gamma_D}(\Omega)^d := \left\{u \in H^1(\Omega)^d, \: u =0 \text{on} \Gamma_D \right\}$ to the mechanical system:
Uad is a set of (smooth) admissible shapes Ω, whose boundaries enclose the non optimizable regions ΓD, ΓN, and another one, Γ0, which may overlap ΓD or ΓN, and which is defined in Section 4 below;
P (Ω) is a constraint functional, meant to enforce the constructibility of shapes by additive processes, whose device is the central issue of the present note;
α is a tolerance threshold.
Obviously, other constraints could be added to (3), e.g. constraints on the volume Vol(Ω) ≔ ∫Ωdx, or the perimeter Per(Ω) ≔ ∫∂Ωds of shapes.
Most popular optimization algorithms for the numerical resolution of (3) rely on the derivatives of J (Ω) and P (Ω) with respect to the domain; these are understood in the framework of Hadamard’s method (see e.g. [1, 15, 18, 21): a function F (Ω) of the domain is shape differentiable if the underlying mapping
from W1, ∞(ℝd, ℝd) into ℝ, is Fréchet differentiable at θ = 0; the corresponding derivative is denoted by F′(Ω)(θ). Often, the deformations θ in (4) are restrained to a subset of W1, ∞(ℝd, ℝd); in the following, we shall consider the sets
where k≥ 1 and lk, ∞(ℝd, ℝd) is the set of k times continuously differentiable functions from ℝd into itself, whose derivatives up to order k are uniformly bounded.
For example, the shape derivative of (2) reads, for deformations θ ∈ Θk, k≥1 (see e.g. [5]):
The geometric attempt to deal with overhang constraints
Definition of the geometric functionals
The perhaps most intuitive way to account for overhang constraints - and the most common one in the literature - relies on criteria involving the angle between the normal vector n ≡ nΩ to the structural boundary and the build direction of the additive machine tool. In this section, for simplicity and without loss of generality, this direction is chosen to be the d-th vector ed of the canonical basis of ℝd. In other terms, this geometric description of overhang features brings into play anisotropic perimeter functionals of the form:
and ν is a threshold angle, penalizes the regions of ∂Ω where the angle between the normal vector n and the negative vertical direction −ed is smaller than ν.
where the ψi : ℝd → ℝ, i = 1, ..., m are given pattern functions, and $n_{\psi_i} := \frac{\nabla \psi_i}{\lvert \nabla \psi_i \lvert}$ are the corresponding normalized gradients, compels n to be close to at least one of the directions nψi.
When it comes to the shape derivative of (5), the result of interest is the following [7].
Proposition 1
The functionPg(Ω) defined by (5) is shape differentiable at any admissible shape Ω ∈uadwhen deformationsθare in Θk, k = 1. Its shape derivative reads:
where κ : ∂Ω → ℝ is the mean curvature of ∂Ω and ∇∂Ωψ ≔ ∇ψ − (∇ψη n) nis the tangential gradient of a smooth enough functionψ : ∂Ω → ℝ..
Insufficiency of the geometric constraint functionals: the ‘dripping effect’
In spite of their simplicity, purely geometric criteria of the form (5) generally fail to prevent the appearance of overhang features.
To illustrate this point, let us anticipate a little on Section 5 where our numerical setting is presented in more details, and consider the benchmark two-dimensional MBB Beam example, as depicted in Figure 2. In a computational domain D with size 6 × 1, the structure is anchored at its bottom-right corner, and the vertical displacement is set to 0 at its bottom-left corner. A unit vertical load f = (0, −1) is applied at the middle of its upper side. From the manufacturing point of view, the shape is assembled from bottom to top, i.e. Γ0 coincides with the lower side of D. Taking advantage of the symmetry of the mechanical problem, only half the working domain D is considered during the optimization process; it is meshed by using 300 × 100 ℚ1 elements.
Starting from the initial shape Ω0 of Figure 3 (top), we solve the compliance minimization problem
(Top) initial and (bottom) optimized shapes for Problem (8) in the two-dimensional MBB Beam test-case of Section 3.2.
with the threshold αv = 0.3 for the volume constraint. The resulting optimized design Ω* is depicted in Figure 3 (bottom) and contains large overhanging regions. These overhangs are of great physical significance for the performance of Ω*; hence, it is expected that their removal will prove difficult.
The results displayed in Figure 4 are typical of the ‘optimized’ shapes resulting from the use of geometric functionals such as Pa(Ω) or Pp(Ω) to penalize the overhangs formed by members of such great structural significance. They are obtained by solving the new problem
in other terms, the constraint functional Pp(Ω) imposes that the orientation of the boundary ∂Ω should comply with that of one of the straight lines with normal vector nψi. These directions nψi are uniformly sampled among the set of those making an angle with the negative vertical direction comprised in $[\frac{\pi}{4},\frac{5\pi}{4}]$ .
One first observes that several parts in the resulting designs do comply with the desired orientation, but these are the parts whose structural significance is negligible.
More importantly, oscillations arose on the boundaries of some members, in particular those which are close to horizontal and bear a great part of the loading. This dripping effect has already been observed in the literature; see for instance [19]. In the present situation, the structure contains regions where the minimization of the compliance J (Ω) urges the formation of large horizontal features; the oscillatory patterns have little impact on the mechanical performance (at least when it is measured in terms of the compliance), but they lead to a dramatic decrease in the values of the geometric functionals Pa(Ω) and Pp(Ω). Obviously, the algorithm prefers creating oscillating boundaries than rearranging the large horizontal bars, which would undermine significantly the structural performance of the shape. In this sense, the geometric criteria are not ‘strong enough’ to steer the algorithm to a different optimization path.
Fig. 4
Optimized shapes resulting from Problem (9) in the two-dimensional MBB Beam example, using (top) φ ≡ φa and the threshold angle ν = 45°, and (bottom) φ ≡ φp and the pattern functions ψi defined in (10).
Description and analysis of the mechanical constraint functionals
Assessing the drawbacks of the geometric functionals of Section 3.2, we now introduce a new mechanical constraint functional Psw(Ω), which relies on a simplified model for the additive construction process.
Definition of the self-weight manufacturing compliance Psw(Ω)
The constraint Psw(Ω) relies on the mechanical situation of Ω in the course of the manufacturing process: Ω is enclosed in a box D = S × (0, H), S ⊂ ℝd − 1, representing the build chamber with a vertical build direction ed. For h ∈ (0, H),
is the intermediate shape describing the stage where Ω is assembled up to height h. The boundary ∂Ωh is here decomposed in a somewhat different fashion from that of Section 2:
$\Gamma_h^l = \partial \Omega_h \setminus \overline{\Gamma_0 \cup \Gamma^u_h}$ is the lateral surface.
Eventually, we define lh := {x ∈ ∂ Ω, xd = h}, the part of the boundary ∂Ω that lies at height h see Figure 5. Each intermediate shape Ωh is clamped on Γ0, and is subjected to gravity effects, accounted for by a body force g ∈ H(ℝd). Its elastic displacement $u_{\Omega_h}^c \in H^1_{\Gamma_0}(\Omega_h)^d$ satisfies:
Fig. 5
Intermediate shape Ωh at the height h during the construction of the final structure Ω: the red zone is the lower boundary Γ0and the blue zone is the upper boundary $\Gamma^u_h.$
Our constraint Psw(Ω) of the total structure Ω aggregates the self-weights of all the intermediate shapes, and therefore deserves the name of self-weight manufacturing compliance:
where j : ℝ → ℝ is a smooth function (in practice, we use j(s) = s).
Remark 1
Let us point out an important bias in this model for the manufacturing process: it implicitly relies on the assumption that every layer of material is constructed instantaneously (and so, it does not take into account the way each such layer is assembled). This feature has important consequences on the interpretation of the numerical results, as we shall see in Section 5.
Shape derivative of Psw(Ω)
The statement of our main result about the shape differentiability of Psw requires additional notations. In this section, Ω ∈ Uad is a given admissible shape. We introduce two open sets O1 ⋐ O2 in ℝd and a smooth function χ: ℝd → such that:
In other words O1 is a neighborhood of the ‘flat regions’ of $\partial \Omega \setminus \overline{\Gamma_0}$, and χ is a cutoff function whereby these regions will be ignored in the analysis.
In this context, the relevant sets for perturbations of Ω are the Banach spaces
In particular, vector fields θ ∈ Xk vanish at the points of $\partial \Omega \setminus \overline{\Gamma_0}$ where the normal vector n is parallel to ed
Theorem 2
The functional Psw (Ω) gievn by (13) is shape differentiable at Ω in the sense that the mapping$\theta \mapsto P_{\text{\rm sw}}(\Omega_\theta)$, from Xk into ℝ is differentiable for k ≥ 1. Its derivative is:
The study of the dependence of Psw(Ω) with respect to the domain Ω is not standard since it involves a continuum of domains Ωh, h ∈ (0, H), obtained from Ω by truncation at fixed heights h. The proof proceeds in three steps:
At first, we prove that $\theta \mapsto P_{\text{\rm sw}}(\Omega_\theta)$ is differentiable when variations θ are horizontal, i.e.
and we prove that (15) - (16) hold in this case. This is possible since, in this case, the truncation at level h (Ωθ)h of the deformed shape Ωθ coincides exactly with the deformation (Id + θ)(Ωh) of the truncated domain Ωh. Hence, in this particular case, Theorem 2 follows almost directly from classical techniques from shape optimization, as in e.g. [12, 21].
We prove that an arbitrary deformation$ θ ∈ Θk can be equivalently represented by a horizontal one ξ(θ) ∈ $\xi(\theta) \in X^k_H$, in the sense that
and we analyze the properties of the mapping θ ↦ ξ(θ) (in particular, we calculate its derivative). This relies on the implicit function theorem after rewriting an identity of the form Ωθ = Ωξ under the form F(θ, ξ) = 0, for θ, ξ ∈ Xk and an appropriate mapping F.
Theorem 2 follows from the conclusions of Steps (i) and (ii) by an application of the chain rule.□
Remark 2
Formulas (15) and (16) have a quite intuitive structure: the shape gradient DΩ(x) of Psw(Ω) at some point $x \in \partial \Omega \setminus \overline{\Gamma_0}$ involves the shape gradients of all the self-weights cΩh of the intermediate structures Ωh lying above x, i.e. such that h > xd
Practical calculation of the mechanical constraint and its derivative
The numerical evaluations of Psw(Ω) and P′sw(Ω)(θ), or equivalently 𝒟Ω, rely on a discretization of the height interval (0,H) with a sequence 0 = h0 < h1 <... < hN = H.The intuitive, ‘0th-order’ method to calculate approximations $P_N^0$ and ${\mathcal D}_N^0$ of Psw(Ω) and 𝒟Ω consists in replacing $c_{\Omega_h}$ and $u_{\Omega_h}^c$ by piecewise constant quantities on each interval Ii := (hi,hi+1) before applying (13) or (15):
Because this procedure is low-order, the subdivision {hi}i = 0,...,N of (0,H) has to be quite fine so that the accuracy of this approximation process is guaranteed; this brings about many numerical resolutions of the elasticity system (11) for the $u_{\Omega_{h_i}}^c$ which is very costly.
The computational efficiency of this ‘0th-order’ method can be improved by constructing a higher-order approximation of the mappings $h \mapsto c_{\Omega_h}$ and $h \mapsto u_{\Omega_h}^c$ on each interval Ii; this is achieved by an interpolation procedure based on the values and the derivatives of these mappings at the hi, which have then to be calculated.
Derivatives of $h \mapsto c_{\Omega_h}$ and $h \mapsto u_{\Omega_h}^c$
Consider a fixed shape Ω ∊ 𝒰ad, and a height h ∈ (0,H satisfying:
$$
\begin{equation}\label{eq.hypnnotver}
\text{For any } x \in \ell_h, \text{ the normal vector } n(x) \text{ is not parallel to } e_d.
\end{equation}
$$
Our first result is about the derivative of $h \mapsto c_{\Omega_h}$ ; see [2] for the proof.
Proposition 3
In the above context, the mapping$h \mapsto c_{\Omega_h}$ is differentiable athand:
We now turn to the issue of differentiating the mapping $h \mapsto u_{\Omega_h}^c$ , an operation which has yet to be given an adequate meaning. In this direction, we prove in [2] that (see also [1, 16] for related notions):
There exists t0 > 0 and a mapping $(-t_0,t_0) \ni t \mapsto \phi_t$ satisfying the following properties:
For t ∈ (−t0,t0), φ_t; is a diffeomorphism of ℝd, mapping Ωh onto Ωh−t such that φ;t(Γ0) = Γ0,
$$\text{The mapping }t \mapsto (\phi_t - \text{Id})\ \text{from}\ (-t_0,\ t_0) \text{ into } W^{1,\infty}(\mathbb{R}^d,\ \mathbb{R}^d) \text{ is of class } \mathcal{C}^1,$$
Introducing $V(x) := \frac{d \phi_t(x)}{dt}\lvert_{t=0} \in W^{1,\infty}(\mathbb{R}^d,\mathbb{R}^d)$ , one has,
$$\text{For } x \in \Gamma_h^u, \:\: V(x)\cdot e_d = -1, \text{ and for } x \in \Gamma_h^l, \:\: V(x) \cdot n(x) = 0.$$
The mapping $t\mapsto u_{\Omega_{h-t}^c} \circ \phi_t$ is differentiable from (-t0,t0) into $H^1_{\Gamma_0}(\Omega_h)^d$ . Its derivative at t = 0 is called the Lagrangian derivative $Y_{\Omega_h}$ of $h \mapsto u_{\Omega_h^c}$ . In general, $Y_{\Omega_h}$ depends on the mapping φt used in its definition.
The quantity $U_{\Omega_h}:= Y_{\Omega_h} - \nabla u_{\Omega_h^c} V$ , which we identify as the Eulerian derivative of $h \mapsto u_{\Omega_h}^c$ is the natural candidate for defining its `derivative'. $U_{\Omega_h}$ is the solution in $H^1_{\Gamma_0}(\Omega_h)^d$ to the system:
In particular, UΩk is independent of the mapping ϕt as long as it fulfills (19).
Practical interpolation algorithm for the calculation of Psw(Ω) and 𝒟Ω
The considerations of Section 4.3.1 suggest the following interpolation procedure for calculating first-order approximations, say $P_N^1$ and ${\mathcal D}_N^1$ , of Psw(Ω) and 𝒟Ω respectively. This allows for an accurate and computationally efficient calculation of these quantities, which uses a coarser subdivision {hi}i = 1,..., N of (0, H) than what is needed for the calculation of the 0th-order approximate values $P_N^0$ and ${\mathcal D}_N^0$ .
For i = 0,..., N, calculate the compliances $c_{\Omega_{h_i}}$ as (12) and the displacements $u_{\Omega_{h_i}}^c$ by solving (11).
For i = 0,..., N, calculate the derivative $\frac{d}{dh}(c_{\Omega_h})|_{h=h_i}$ of the compliance owing to Proposition 3.
For i = 1,..., N, calculate the Eulerian derivative $U_{\Omega_{h_i}}$ at hi by using (20).
On each interval Ii, i = 0, ..., N – 1, the compliance cΩh is approximated by a cubic spline ${\tilde c_i}\left( h \right)$ which is uniquely (and explicitly) determined by the data:
notice that the above relation does make sense for x ∈ Ωh regardless of the height h ∈ (hi, hi + 1) since $u_{\Omega_{h_{i+1}}}^c$ and UΩhi+1 are both well-defined on Ωh ⊂ Ωhh + 1.
Alternative mechanical models for the behavior of intermediate shapes
Our definition of the constraint Psw(Ω) in Section 4.1 is guided by the description (11) of the physical behavior of the intermediate shapes Ωh during the construction process: they are only subjected to gravity effects.
Interestingly, different mechanical models could be considered instead of (11), including fictitious ones, thereby giving different focuses to the constraint functional. For instance, one may imagine replacing the uniform gravity load g in (11) by an artificial force, applied only on the upper region of each intermediate shape Ωh: Ωh is now subjected to a body force gh defined by:
where δ > 0 is a small parameter, and g ∈ H1(ℝd)d is a given function. The elastic displacement $u_{\Omega_h}^a$ of Ωh in this context is the unique solution in $H^1_{\Gamma_0}(\Omega_h)^d$ to the system:
As we shall see in Section 5, this formulation is well-suited when it comes to penalizing more specifically the upper region of each intermediate shape Ωh.
As far as the shape derivative of Puw(Ω) is concerned, the exact same proof as that of Theorem 2 can be worked out, taking advantage of the definition (23) of gh, and the conclusions of this Theorem extend verbatim to this new case. Also, using a similar analysis to that of Section 4.3 (and working out similar calculations as in [2]), it is possible to calculate the ‘derivatives’ of the mappings $h \mapsto c_{\Omega_h}^a$ and $h \mapsto u_{\Omega_h}^a$ , which results in expressions similar to (18) and (20), up to additional terms accounting for the dependence of the load gh on h.
Let us eventually emphasize that our choice of the linearized elasticity systems (11) or (24) for the description of the behavior of the intermediate shapes is merely incidental. One could for instance rely on a similar construction, involving instead the heat equation, so as to model cooling effects within the intermediate shapes, and thereby residual stresses; see [4] about this idea.
Numerical illustrations
We end this note with several illustrations of the above theoretical considerations; in Section 5.2, we test the algorithms of Section 4.3 for the evaluation of the constraint Puw(Ω) and its derivative; the subsequent examples (Section 5.3 and below) discuss several numerical strategies for preventing the appearance of overhang features during the shape optimization process. Before entering the core of the matter, we briefly outline our numerical setting.
A few words about the numerical implementation
In all our examples, we rely on the level set method on a fixed computational mesh to represent shapes and their deformations; see [17, 20] about the level set method, and [5, 22] about its use in shape and topology optimization. In a nutshell, one shape Ω, enclosed in a larger, fixed working domain D, is regarded as the negative subdomain of a scalar ‘level set’ function ϕ: D → ℝ, that is:
The evolution in time of a shape Ω(t), driven by a velocity field with normal component V(t, x), can be modeled by the following Hamilton-Jacobi equation, the solution ϕ(t, x) of which, is a level set function for Ω(t):
In our applications, the (scalar) velocity field V(t, x) stems from the resolution of Problem (3) by means of an SLP-type optimization algorithm similar to that presented in [9], and the (pseudo-) time t represents the descent step.
From a numerical point of view, the working domain D is a box in two or three space dimensions; it is discretized by means of a Cartesian mesh 𝒢, i.e. 𝒢 is composed of square elements in 2d and cubes in 3d. The level set function ϕ is discretized at the vertices of 𝒢 and the Hamilton-Jacobi equation (27) is solved by using an explicit second-order upwind scheme on 𝒢, as presented e.g. in [20].
Since the computational mesh 𝒢 is fixed throughout the optimization process of the shape Ω, no mesh of Ω is available for the Finite Element resolution of linearized elasticity systems of the form (1). To circumvent this difficulty, we rely on the so-called ‘ersatz-material’ approximation which consists in filling the void D/Ω̄ with a very soft material with Hooke's tensor εA (in practice ε = 10‒3), thus transferring a system posed on Ω into an approximate one posed on D; see for instance [5].
In all examples, the Young’s modulus of the considered elastic material is normalized as E = 1 and the Poisson’s ratio is set to v = 0.33.
Our first numerical example aims to assess the computational efficiency of the first-order interpolation algorithm proposed in Section 4.3 for the calculation of the constraint Psw(Ω) and its derivative 𝒟Ω.
Let us consider again the two-dimensional MBB Beam described in Section 3.2. At first, we calculate the functional Psw(Ω) and its shape derivative 𝒟Ω in the particular case where Ω = Ω0 (the shape displayed in Figure 3 (top)), by using a uniform subdivision of (0, H) made of 100 layers and the 0th-order approximation scheme, i.e. we evaluate $P_{100}^0$ and ${\mathcal D}_{100}^0$ , which serve as reference values for the comparisons in this section. We then calculate the 0th- and 1st-order approximations $P_N^i$ and ${\mathcal D}_N^i$ , i = 0, 1 associated (0, H) made of N intervals with equal length. We are interested in the behavior of the relative errors:
The results are displayed on Figure 6 (bottom): while the 1st-order approximation method does not bring a lot of improvement when it comes to evaluating the constraint functional P(Ω), it allows for a substantial gain (i.e. a faster convergence with respect to the number N of subdivisions) in the evaluation of its derivative.
Fig. 6
Relative errors of the 0th- and 1st-order approximations of (top) Psw(Ω0) and (bottom) its derivative 𝒟Ω0.
Still in the context of the two-dimensional MBB Beam described in Section 3.2, we turn to the compliance minimization problem under volume constraint:
Optimization of the two-dimensional MBB Beam using a constraint on the self-weight manufacturing compliance Psw(Ω)
Still in the context of the two-dimensional MBB Beam described in Section 3.2, we turn to the compliance minimization problem under volume constraint:
$$\begin{array}{l}\mathop {\min }\limits_\Omega {\text{ }}J(\Omega )\\
{\text{s}}{\text{.t}}{\text{. Vol}}(\Omega ) \leqslant {\alpha _{{v}}}{\text{Vol}}({{D}})\end{array},\:\:\:{\text{ where }}\,J(\Omega ){\text{is the compliance}}(2).$$
We first solve (28), starting from the initial design Ω0 of Figure 3 (top), with threshold value αv =0.3. For the reader’s convenience, the resulting design Ω* is reprinted from Figure 3 (bottom) in Figure 7 (a).
We then add our mechanical constraint Psw(Ω) to the problem (28), and now solve:
for different values of the parameter αc; the results are represented on Figure 7 (b)-(d). Understandingly, decreasing values of αc result in structures which are more rigid from a manufacturing point of view (i.e. with lower values of Psw(Ω)) and more compliant from a structural perspective (i.e. with higher values of J(Ω).
Let us observe however that the optimized shapes still contain overhangs, which are mainly concentrated on their upper regions. This phenomenon arises for mainly two reasons:
As mentioned in Remark 1, the fact that each layer is assumed to be deposited instantaneously has an impact on the values of the self-weights cΩh given by (12): the rigidity of completely flat overhangs (as those appearing in figures 7 (b)-(d)) is actually overestimated; roughly speaking, our modeling only involves the intermediate stages of the construction process where each layer is completely assembled, and where these flat regions are thereby connected to the lower structure. Hence, all the situations where these regions are hanging over void (for instance, at the moments when they are only half-constructed, and when the self-weight would typically take large values) are overlooked.
Secondly, considering a uniform self-weight loading for each intermediate shape Ωh leads to a concentration of the shape gradient of J(Ω) at regions closer to Γ0. Intuitively, this comes from the fact that such regions appear in a greater number of intermediate shapes, which is reflected on the shape derivative of Psw(Ω); see Theorem 2. Thus, the algorithm favours the elimination of the overhangs located in these regions, while upper parts are comparatively less influenced by this requirement.
Fig. 7
Optimized shapes for the two-dimensional MBB Beam example of Section 5.3: (a) optimized shape Ω*for Problem (8) (i.e. without additive manufacturing constraints), and optimized shapes for Problem (29) using parameters (b)αc = 0:50, (c) αc = 0:30, and (d) αc = 0:10.
One remedy for the former issue would consist in using a more precise, ‘pixel by pixel’ model of the construction process introduced in Section 4.1: one could indeed consider the trajectory of the machine tool during the assembly of each layer of the shape, and take as ‘intermediate shapes’ several snapshots in the assembly of each layer. Unfortunately, beyond difficult modeling issues, relying on this idea would cause the computational cost (which is already significant) of the constraint functional to increase dramatically. Therefore, this proposition is not examined in the present article.
A second remedy amounts to changing the mechanical problem (11) characterizing the behavior of the intermediate shapes. More precisely, as proposed in Section 4.4, we consider fictitious body forces applied on the upper region of these shapes, thus better penalizing thin horizontal patterns in the intermediate shapes Ωh. This is the purpose of the next section.
Improvement brought by the upper-weight manufacturing compliance Puw(Ω)
We now take on the two-dimensional MBB Beam example of Section 5.3, replacing the self-weight manufacturing compliance Psw(Ω) with the upper-weight manufacturing compliance Puw(Ω) in the formulation (29) of the optimization problem, with the expectation that it shows a better ability when it comes to removing overhanging features.
The thickness δ of the regions where g is applied in the definition (23) of the body force gh is of the order of the mesh size: δ = Δx.
One observes that the algorithm tends to add more features oriented along the build direction and to connect them together by creating small archs, which have optimal rigidity for self-weight loadings. The results are much more satisfying than those obtained by using the functional Psw(Ω) when it comes to removing overhanging features.
Fig. 8
Optimized shapes for the two-dimensional MBB Beam example of Section 5.4, solving Problem (29) with the upper-weight manufacturing compliancePsw(Ω) and parameters (a)αc = 0:30, (b) αc = 0:10, (c) αc = 0:05, and (d) αc = 0:03.
A three-dimensional example
Our last example is about the optimization of a three-dimensional bridge, as depicted in Figure 9: the dimensions of the working domain D are 6 × 1 × 1. The structure is clamped at its lower-right corners, while at its lower-left corners, the vertical displacement is prevented. A uniform load f = (0, 0, -1) is applied on the upper side of the structure. From the manufacturing point of view, the shape is constructed from top to bottom, i.e. Γ0 coincides with the upper side of D (i.e. the deck of the bridge). Due to the double symmetry of the mechanical problem, only one quarter of D is meshed and the cubic grid 𝒢 is composed of 90 × 15 × 30 elements.
We start by minimizing the structural compliance J(Ω) of the bridge under volume constraint, taking the full working domain D as initial shape; i.e. we solve (8) with αv = 0.10 Figure 10 (left column) shows the resulting optimized shape Ω*, which presents several overhangs of significant structural importance.
Fig. 9
Setting of the three-dimensional bridge test-case.
We then add a constraint on the upper-weight manufacturing compliance Puw(Ω), but instead of solving (29) we rather consider the following optimization problem:
where we look for a shape having at least the same rigidity as Ω*, but which is also more rigid from a manufacturing perspective. Figures 10 (right column) and 11 (left column) show the optimized shapes for αc = 0.70 and αc = 0.10 respectively. For αc = 0.70 only slight changes are observed: the upper bar becomes thinner and the central bars become thicker and are relocated closer to Γ0. For αc = 0.10 the changes are drastic and the same trend as in the two-dimensional MBB beam example is observed.
Fig. 10
Optimized designs for the three-dimensional bridge example of Section 5.5, (left) without manufacturing constraints, (right) solving Problem (30) withαc = 0:7.
Fig. 11
(Left) Different views of the optimized shape for the three-dimensional bridge example of Section 5.5, solving Problem (30) withαc = 0:1; (right) another view on the three-dimensional bridges for Problem (30) with (top) no manufacturing constraint, (middle)αc = 0:7 and (bottom) ac = 0:1.
