∎ [table]capposition=above [table]captionskip=0cm

11institutetext: Authors 22institutetext: Department of Aerospace and
Mechanical Engineering,
University of Liège, 4000 Liège, Belgium.

Denis Trillet
22email<EMAIL_ADDRESS>
# Note on the minimum length scale and its defining parameters

Analytical relationships for Topology Optimization based on uniform
manufacturing uncertainties.

Denis Trillet  Pierre Duysinx  Eduardo Fernández

(Received: date / Accepted: date)

###### Abstract

The robust topology optimization formulation that introduces the eroded and
dilated versions of the design has gained increasing popularity in recent
years, mainly because of its ability to produce designs satisfying a minimum
length scale. Despite its success in various topology optimization fields, the
robust formulation presents some drawbacks. This paper addresses one in
particular, which concerns the imposition of the minimum length scale. In the
density framework, the minimum size of the solid and void phases must be
imposed implicitly through the parameters that define the density filter and
the smoothed Heaviside projection. Finding these parameters can be time
consuming and cumbersome, hindering a general code implementation of the
robust formulation. Motivated by this issue, in this article we provide
analytical expressions that explicitly relate the minimum length scale and the
parameters that define it. The expressions are validated on a density-based
framework. To facilitate the reproduction of results, MATLAB codes are
provided.

As a side finding, this paper shows that to obtain simultaneous control over
the minimum size of the solid and void phases, it is necessary to involve the
3 fields (eroded, intermediate and dilated) in the topology optimization
problem. Therefore, for the compliance minimization problem subject to a
volume restriction, the intermediate and dilated designs can be excluded from
the objective function, but the volume restriction has to be applied to the
dilated design in order to involve all 3 designs in the formulation.

###### Keywords:

Length Scale Robust Design SIMP

††journal: PREPRINT

##  1 Introduction

Since the seminal work of Bendsøe and Kikuchi (1988), topology optimization
has experienced huge advances and nowadays is being massively adopted in the
industry (Pedersen and Allinger, 2006; Zhou et al, 2011; Zhu et al, 2016).
Among the successful advancements, one can mention the famous density method
under the SIMP interpolation scheme (Bendsøe, 1989). The well-known
shortcomings of SIMP led to a succession of improvements seeking to avoid the
mesh dependency, the checkerboard patterns and the presence of intermediate
densities. To date, one of the most effective approaches to dealing with the
ill-effects of SIMP is the robust design approach that brings the eroded and
dilated versions of the design (Sigmund, 2009). This method considers
manufacturing errors that may incur in a uniformly thinner or uniformly
thicker component compared to the blueprint layout. Sigmund (2009) proposed a
robust formulation that maximizes the performance of the worst performing
design among the eroded, dilated and reference (intermediate) designs. This
means, the formulation guarantees a design with good performance even if it is
eventually eroded or dilated during the manufacturing process. Interestingly,
the robust formulation yields a reference (intermediate) design that features
minimum member size and minimum cavity size (Wang et al, 2011).

In addition to imposing minimum length scale in topology optimization, it has
been seen that the robust formulation provides a more stable convergence than
other projection or filtering strategies known so far, which allows to reach
almost discrete solutions in the density method (Wang et al, 2011). The
formulation has been recently applied in combination with existing methods
that allow to impose stress limits (da Silva et al, 2019; Silva et al, 2020),
overhang angle constraints (Pellens et al, 2018), maximum size restrictions
(Fernández et al, 2020), and geometric non-linearities (Lazarov and Sigmund,
2011; Silva et al, 2020), not only in the density method but also in the level
set method (Chen and Chen, 2011; Andreasen et al, 2020).

Despite proving its effectiveness in several fields of application (Wang et
al, 2011b; Christiansen et al, 2015), the robust design approach presents some
drawbacks with regard to its implementation. For example, the method requires
boundary treatments with respect to the density filter, since the filtering
region can be split at the edges of the design domain affecting the imposed
minimum size (Clausen and Andreassen, 2017; Kumar and Fernández, 2021). In
addition, due to the erosion and dilation distances with respect to the
intermediate design, boundary conditions can be disconnected from the designs
involved in the formulation (Clausen and Andreassen, 2017). Another difficulty
posed by the method is that the minimum length scale must be implicitly
imposed through the parameters defining the density filter and the Heaviside
projection. This drawback has been addressed in the literature using numerical
(Wang et al, 2011) and analytical (Qian and Sigmund, 2013) approaches. The
numerical approach consists of applying the density filter and the Heaviside
projection to a 1D design, measuring the resulting length scale, and repeating
the process several times with different filter and projection parameters to
subsequently construct a graph that relates the involved parameters. The
analytical method consists of applying the density filter as a convolution
integral in a continuous 1D design. The solution of the integration yields
explicit relationships of the projection/filtering parameters with the
obtained length scale (Qian and Sigmund, 2013). As the three field scheme is a
scalar function, the relationships obtained from the 1D designs remain valid
for 2D and 3D (Wang et al, 2011).

The numerical and analytical relationships reported by Wang et al (2011) and
Qian and Sigmund (2013) are intended for the particular case where the minimum
size of the void phase is set equal to that of the solid phase. This
simplifies and facilitates the procedures to explicitly relate the desired
minimum length scale to the filter and projection parameters, but in this way,
the minimum size of the solid phase cannot be different from the minimum size
of the void phase. Nonetheless, due to the rapid popularization of the robust
formulation, it becomes necessary to find a method that allows to impose
different minimum sizes for each phase, in addition to provide access to other
geometric information. For example, the erosion and dilation distances with
respect to the intermediate design have been required to impose a maximum
member size (Fernández et al, 2020) and to obtain designs tailored to the size
of a deposition nozzle (Fernández et al, 2021). To date, the erosion and
dilation distances have been reported for specific minimum length scales,
which hinders widespread applicability of the above methods.

The aim of this work is to provide a method to obtain the filter and
projection parameters that impose user-defined minimum length scales, where
the size of the solid phase can be defined differently from that of the void
phase. In addition, for the desired minimum length scales, the erosion and
dilation distances are provided. To this end, we extend the applicability of
the analytical method proposed by Qian and Sigmund (2013) and the resulting
relationships are validated using the numerical method proposed by Wang et al
(2011) and a set of 2D topology-optimized designs. To facilitate the
replication of results and the application of the methods discussed in this
paper, MATLAB codes are provided.

The remainder of the document is developed as follows. Section 2 introduces
the robust topology optimization formulation and the case studies used to
asses the analytical expressions that relate the minimum len-gth scale to the
parameters that impose it. Section 3 presents the analytical method proposed
by Qian and Sigmund (2013) and describes the main contribution of our work,
which is the extension of the analytical approach to allow choosing
independent minimum sizes for the solid phase and the void phase. Section 4
validates the analytical expressions using the numerical method proposed by
Wang et al (2011). Section 5 discusses the sources of error that are inherent
to the analytical method. Section 6 assesses the analytical expressions on 2D
topology-optimized designs. Section 7 gathers the final conclusions of this
work while Section 8 provides the codes that allow for the replication of
results.

##  2 Problem Definition

The analytical expressions that relate the minimum length scale to the
parameters that define it are developed for the robust topology optimization
formulation based on the eroded, intermediate and dilated designs. This work
considers the density approach based on the SIMP interpolation scheme
(Bendsøe, 1989), even though the proposed methodology can be applied to other
topology optimization approaches with little efforts.

Like most works in the literature, the eroded, intermediate and dilated
designs that constitute the robust formulation are built using a three-field
scheme (Sigmund and Maute, 2013). The first field, denoted by $\bm{\rho}$,
corresponds to the design variables. The second field, denoted by
$\bm{\tilde{\rho}}$, is obtained by a weighted average of the design variables
within a circle of radius $r_{\mathrm{fil}}$. The third field, denoted by
$\bm{\bar{\rho}}$, is obtained by projecting the components of the filtered
field towards 0 or 1. The filter and projection functions are identical to
those provided by Wang et al (2011), however, these are reminded herein for
the sake of clarity as they define the minimum length scale obtained in the
optimized design.

The filter of design variables, or the density filter (Bruns and Tortorelli,
2001; Bourdin, 2001), is defined as follows:

$\tilde{\rho}_{i}=\frac{\displaystyle\sum_{j=1}^{N}\rho_{j}\mathrm{v}_{j}w(\mathbf{x}_{i}-\mathbf{x}_{j})}{\displaystyle\sum_{j=1}^{N}{\mathrm{v}}_{j}w(\mathbf{x}_{i}-\mathbf{x}_{j})}\;\;,$
(1)

where $\rho_{j}$ is the design variable associated to the element $j$ and
$\tilde{\rho}_{i}$ is the filtered variable associated to the element $i$.
$\mathrm{v}_{j}$ is the volume of the element $j$ and
$w(\mathbf{x}_{i}-\mathbf{x}_{j})$ is the weigh of $\rho_{j}$ in the
definition of $\tilde{\rho}_{i}$. As is common practice in the literature, the
weighting function $w(\mathbf{x}_{i}-\mathbf{x}_{j})$ is defined as a linear
and decreasing function with respect to the distance between the elements $i$
and $j$, as follows:

$w(\mathbf{x}_{i}-\mathbf{x}_{j})=\mathrm{max}\left(0,\>1-\frac{\|\mathbf{x}_{i}-\mathbf{x}_{j}\|}{r_{\mathrm{fil}}}\right)\;\;,$
(2)

where $\mathbf{x}_{i}$ and $\mathbf{x}_{j}$ represent the centroid of the
elements $i$ and $j$, respectively. It is recalled that $r_{\mathrm{fil}}$ is
the radius of the density filter.

To reduce the amount of intermediate densities pre-sent in the filtered field
and to build the designs that constitute the robust formulation, the projected
field is obtained with the following smoothed Heaviside function (Wang et al,
2011):

$\bar{\rho}_{i}=\frac{\tanh{(\beta\eta)}\tanh{(\beta(\tilde{\rho}_{i}-\eta))}}{\tanh{(\beta\eta)}\tanh{(\beta(1-\eta))}}$
(3)

where $\beta$ and $\eta$ control the steepness and the threshold of the
projection, respectively. The eroded, intermediate and dilated designs,
denoted by ${\bm{\bar{\rho}}}_{\mathrm{ero}}$,
${\bm{\bar{\rho}}}_{\mathrm{int}}$ and ${\bm{\bar{\rho}}}_{\mathrm{dil}}$, are
obtained from the smoothed Heaviside function of Eq. (3) by using the same
$\beta$ but with different thresholds ${\eta_{\mathrm{ero}}}$,
$\eta_{\mathrm{int}}$ and $\eta_{\mathrm{dil}}$, thus leading to
${{\bm{\bar{\rho}}}_{\mathrm{ero}}}(\bm{\tilde{\rho}},\beta,{\eta_{\mathrm{ero}}})$,
${{\bm{\bar{\rho}}}_{\mathrm{int}}}(\bm{\tilde{\rho}},\beta,{\eta_{\mathrm{int}}})$
and
${{\bm{\bar{\rho}}}_{\mathrm{dil}}}(\bm{\tilde{\rho}},\beta,{\eta_{\mathrm{dil}}})$.

(a)
(b)

Figure 1: (a) Heat sink and (b) force inverter design domains considered in
this study.

To assess the scope of the equations that provide the necessary parameters to
impose the minimum length scale, different topology optimization problems are
sol-ved with variations in the desired minimum length scale. The minimum size
of the optimized designs is measured and compared against the intended values.
In this work, the problems that are chosen to illustrate the developments are
the heat conduction and the non-linear compliant mechanism design. The former
aims at minimizing the thermal compliance subject to a volume restriction
(Wang et al, 2011), whilst in the latter, the formulation of a force inverter
is considered where an output displacement is maximized for a given input
force (Sigmund, 1997). The design domains are shown in Fig. 1. According to
the robust design approach, the topology optimization problems can be written
as follows:

$\begin{split}\text{min}\quad&\mathrm{max}\left(c({{\bm{\bar{\rho}}}_{\mathrm{ero}}}),\;c({{\bm{\bar{\rho}}}_{\mathrm{int}}})\;,c({{\bm{\bar{\rho}}}_{\mathrm{dil}}})\right)\\\
\text{s.t.:}\quad&\mathbf{v}^{\intercal}\>{{\bm{\bar{\rho}}}_{\mathrm{dil}}}\leq
V^{*}_{\mathrm{dil}}(V^{*}_{\mathrm{int}})\\\ &0\leq\rho_{i}\leq
1\;\;,\;\;i=1,...,N\;\;,\end{split}$ (4)

where $c$ represents the thermal compliance of the heat sink or the output
displacement of the force inverter, $V^{*}_{\mathrm{dil}}$ is the upper bound
of the volume restriction and N is the total amount of design variables. As
the design intended for manufacturing is the intermediate design, the upper
bound of the volume constraint is scaled according to the user-defined limit
$V^{*}_{\mathrm{int}}$. The optimization problem and the scaling process of
the volume restriction are the same as described by Wang et al (2011) for the
heat conduction problem. In the nonlinear force inverter, the method to deal
with the mesh distorsion appearing in low density elements is the same as
explained by Wang et al (2014). To avoid overextending the contents of the
manuscript, the interested reader is refered to the cited articles.

##  3 Analytical minimum length scale

Comprehensive numerical tests have shown that the robust formulation in Eq.
(4) yields intermediate designs that have the same topology as their eroded
and dilated projections (Wang et al, 2011). In other words, the erosion
process does not destroy the solid members of the optimized topology
(intermediate), it just makes them thinner. This implies that a minimum size
of the solid members in the intermediate design is imposed according to the
erosion distance. Similarly, the dilation process does not close the cavities
of the optimized topology (intermediate), it only makes them narrower. This
implies that the minimum size of the void phase is imposed according to the
dilation distance.

Considering a two-dimensional design domain, the minimum size of the solid
phase is defined by the radius $r_{\mathrm{min.Solid}}$ of the largest circle
that can be circumscribed into the smallest solid member of the topology,
while the minimum size of the void phase is defined by the radius
$r_{\mathrm{min.Void}}$ of the largest circle that can be circumscribed into
the smallest cavity of the topology. It has been shown that the minimum length
scale, $r_{\mathrm{min.Solid}}$ and $r_{\mathrm{min.Void}}$, are defined by
the filter radius $r_{\mathrm{fil}}$ and the projection thresholds
$\eta_{\mathrm{ero}}$, $\eta_{\mathrm{int}}$ and $\eta_{\mathrm{dil}}$.
Therefore, if specific values are to be imposed for the minimum length scale,
these must be implicitly imposed through the 4 parameters that define the
density filter and the Heaviside projection.

The dependence of the minimum length scale on the projection and filtering
parameters hurdles a general implementation of the robust formulation,
specially in commercial codes, as there is no prompt approach to find the
relationship between the involved parameters and the desired length scale.
Encouraged by this shortcoming, in this section we present and extend the
applicability of the analytical procedure proposed by Qian and Sigmund (2013),
which is described below.

To obtain an explicit relation between the minimum length scale and the
parameters that define it, Qian and Sigmund (2013) proposed to apply the
three-field scheme over a uni-dimensional and continuous design domain. For
instance, to obtain the minimum size of the solid phase, a design domain
$\bm{\rho}(x)$ centered at the coordinate $x_{m}$ and containing solid
elements on a stretch of size $h$ is assumed, as shown in Fig. 2(a). The
continuous form of the density filter reads as follows:

$\begin{matrix}{\tilde{\rho}}(x_{i})&=\frac{\displaystyle\int_{x_{i}-r_{\mathrm{fil}}}^{x_{i}+r_{\mathrm{fil}}}\bm{\rho}(x)\left(1-\frac{\left|x_{i}-x\right|}{r_{\mathrm{fil}}}\right)dx}{\displaystyle\int_{x_{i}-r_{\mathrm{fil}}}^{x_{i}+r_{\mathrm{fil}}}\left(1-\frac{\left|x_{i}-x\right|}{r_{\mathrm{fil}}}\right)dx}\vspace{4mm}\\\
&=\displaystyle\int_{x_{i}-r_{\mathrm{fil}}}^{x_{i}+r_{\mathrm{fil}}}\frac{\bm{\rho}(x)}{r_{\mathrm{fil}}}\left(1-\frac{\left|x_{i}-x\right|}{r_{\mathrm{fil}}}\right)dx\end{matrix}$
(5)

where $x_{i}$ is any coordinate of $x$. For example, by choosing a filter
radius greater than $h/2$, the one-dimensional filtered field in Fig. 2(b) is
obtained. The filtered field can be projected using the smoothed Heaviside
function of Eq. (3). To simplify the analysis, an infinite steepness parameter
($\beta\to\infty$) is assumed. For example, for an Heaviside threshold
$\eta_{i}=0.2$, the design of Fig. 2(c) is obtained. The size of the solid
phase in the projected design is defined by the length $L_{i}$, which can be
obtained by finding $x_{i}$ from the equation $\tilde{\rho}(x_{i})=\eta_{i}$.
However, to this extent, the size of the solid phase is a function
$L_{i}(r_{\mathrm{fil}},\eta_{i},h)$ that depends on the assumed $h$ value. To
determine the value of $h$ and to discover the explicit relation between the
size $L$ of the projected field and the filter and projection parameters, it
is necessary to refer to the foundation of the robust formulation.

(a) 1D domain

(b) Filtered Field

(c) Projected Field

(d) Filtered field

Figure 2: The three field scheme applied to a one-dimensional design domain in
order to obtain the analytical minimum length scale.

An erosion projection removes solid material from the surface of the reference
design. This operation enlarges the cavities and thins the structural members
present in the intermediate reference design (Sigmund, 2009). When robustness
with respect to erosion is desired, the structural members must be present in
both the eroded and the reference designs. It has been shown that a sufficient
condition to preserve robustness is that the solid member in the eroded design
has to be projected by at least an infinitesimal size, i.e.
$L({\eta_{\mathrm{ero}}})\approx 0$, as shown in Fig 2(d). This condition
allows determining the assumed size $h$ by solving the following equation:

$\tilde{\rho}(x_{\mathrm{ero}})=\int_{-h/2}^{h/2}\frac{1}{r_{\mathrm{fil}}}\left(1-\frac{|x|}{r_{\mathrm{fil}}}\right)dx={\eta_{\mathrm{ero}}}$
(6)

where the reference system is placed at $x_{\mathrm{ero}}$ for integration. By
solving the expression in Eq. (6), the size $h$ is obtained as:

$h=2r_{\mathrm{fil}}(1-\sqrt{1-{\eta_{\mathrm{ero}}}})$ (7)

By obtaining the distance $h$ that produces an eroded projection of
infinitesimal size, it is possible to relate the minimum size to the filter
radius for any projection whose threshold meets
$\eta_{i}<{\eta_{\mathrm{ero}}}$. A particular case is to choose
$\eta_{i}=0.5$, which is the value chosen by Qian and Sigmund (2013) to define
the intermediate design. However, to extend the analytical method we use an
arbitrary $\eta_{i}$ such that ${\eta_{\mathrm{i}}}<{\eta_{\mathrm{ero}}}$.

As mentioned above, the size of the solid phase $L$ in the projected field can
be obtained by equating $\tilde{\rho}(x_{i})$ to $\eta_{i}$. Since the
condition of robustness is imposed on $h$ (Eq. 7), the length corresponds to
the minimum size, i.e. $L=2r_{\mathrm{min.Solid}}$. The expression that allows
to relate the minimum size to the filter radius is as follows:

${\tilde{\rho}}(x_{i})=\displaystyle\int_{x_{i}-r_{\mathrm{fil}}}^{x_{i}+r_{\mathrm{fil}}}\frac{\bm{\rho}(x)}{r_{\mathrm{fil}}}\left(1-\frac{\left|x_{i}-x\right|}{r_{\mathrm{fil}}}\right)dx=\eta_{i}$
(8)

For the integration, it is convenient to place the origin of the reference
system at $x_{i}$. To this end, 4 situations must be considered to define the
integration limits, which are summarized in Fig. 3. For example, if the length
of the projected field $L$ is greater than $h$ and the filter radius is
greater than $(L+h)/2$ (Fig. 3(a)), then Eq. (8) becomes:

${\tilde{\rho}}(x_{i})=\displaystyle\int_{\frac{L-h}{2}}^{L-\frac{L-h}{2}}\frac{1}{r_{\mathrm{fil}}}\left(1-\frac{\left|x\right|}{r_{\mathrm{fil}}}\right)dx=\eta_{i}$
(9)

(a) $L\geq h$ and $r_{\mathrm{fil}}\geq\frac{L+h}{2}$.

(b) $L<h$ and $r_{\mathrm{fil}}\geq\frac{L+h}{2}$.

(c) $L\geq h$ and $r_{\mathrm{fil}}<\frac{L+h}{2}$.

(d) $L<h$ and $r_{\mathrm{fil}}<\frac{L+h}{2}$.

Figure 3: Situations to be considered when defining the limits of integration. Phase | Conditions |  |  Minimum size
---|---|---|---
Solid | $\eta_{i}\leq 2{\eta_{\mathrm{ero}}}-1$ | and | $\eta_{i}\geq 0.5$ |  | $\displaystyle\frac{2r_{\mathrm{min.Solid}}}{r_{\mathrm{fil}}}=2\sqrt{2-2\eta_{i}}-2\sqrt{1-{\eta_{\mathrm{ero}}}}$
$\eta_{i}\geq 2{\eta_{\mathrm{ero}}}-2+2\sqrt{1-{\eta_{\mathrm{ero}}}}$ | and | $\eta_{i}>2{\eta_{\mathrm{ero}}}-1$ |  | $\displaystyle\frac{2r_{\mathrm{min.Solid}}}{r_{\mathrm{fil}}}=2\sqrt{{\eta_{\mathrm{ero}}}-\eta_{i}}$
$\eta_{i}<4-4\sqrt{1-{\eta_{\mathrm{ero}}}}-2{\eta_{\mathrm{ero}}}$ | and | $\eta_{i}<0.5$ |  | $\displaystyle\frac{2r_{\mathrm{min.Solid}}}{r_{\mathrm{fil}}}=4-2\sqrt{1-{\eta_{\mathrm{ero}}}}-2\sqrt{2\eta_{i}}$
$\eta_{i}\geq 4-4\sqrt{1-{\eta_{\mathrm{ero}}}}-2{\eta_{\mathrm{ero}}}$ | and | $\eta_{i}<2\eta_{\mathrm{ero}}-2+2\sqrt{1-{\eta_{\mathrm{ero}}}}$ |  | $\displaystyle\frac{2r_{\mathrm{min.Solid}}}{r_{\mathrm{fil}}}=2-\frac{\eta_{i}}{1-\sqrt{1-{\eta_{\mathrm{ero}}}}}$
Void | $\eta_{i}\geq 2{\eta_{\mathrm{dil}}}$ | and | $\eta_{i}\leq 0.5$ |  | $\displaystyle\frac{2r_{\mathrm{min.Void}}}{r_{\mathrm{fil}}}=2\sqrt{2\eta_{i}}-2\sqrt{{\eta_{\mathrm{dil}}}}$
$\eta_{i}\leq 2\eta_{d}+1-2\sqrt{{\eta_{\mathrm{dil}}}}$ | and | $\eta_{i}<2{\eta_{\mathrm{dil}}}$ |  | $\displaystyle\frac{2r_{\mathrm{min.Void}}}{r_{\mathrm{fil}}}=2\sqrt{\eta_{i}-{\eta_{\mathrm{dil}}}}$
$\eta_{i}\geq 4\sqrt{{\eta_{\mathrm{dil}}}}-2{\eta_{\mathrm{dil}}}-1$ | and | $\eta_{i}>0.5$ |  | $\displaystyle\frac{2r_{\mathrm{min.Void}}}{r_{\mathrm{fil}}}=4-2\sqrt{{\eta_{\mathrm{dil}}}}-2\sqrt{2-2\eta_{i}}$
$\eta_{i}<4\sqrt{{\eta_{\mathrm{dil}}}}-2{\eta_{\mathrm{dil}}}-1$ | and | $\eta_{i}>2{\eta_{\mathrm{dil}}}+1-2\sqrt{{\eta_{\mathrm{dil}}}}$ |  | $\displaystyle\frac{2r_{\mathrm{min.Void}}}{r_{\mathrm{fil}}}=2-\frac{1-\eta_{i}}{1-\sqrt{{\eta_{\mathrm{dil}}}}}$
Table 1: The explicit relationship between the minimum length scale and the
filter and projection parameters.

Solving the integral of Eq. (9) leads to the following expression:

${\tilde{\rho}}(x_{i})=\frac{h}{r_{\mathrm{fil}}}(1-\frac{L}{2r_{\mathrm{fil}}})=\eta_{i}$
(10)

Finally, replacing Eq. (7) in (10):

$\frac{2r_{\mathrm{min.Solid}}}{r_{\mathrm{fil}}}=2-\frac{\eta_{i}}{1-\sqrt{1-{\eta_{\mathrm{ero}}}}}$
(11)

Eq. (11) explicitly relates the minimum size of the solid phase
($r_{\mathrm{min.Solid}}$), of a projected field defined by $\eta_{i}$, with
the filter radius ($r_{\mathrm{fil}}$) and the erosion threshold
($\eta_{\mathrm{ero}}$). However, Eq. (11) is only valid for $L>h$ and
$r_{\mathrm{fil}}\geq(L+h)/2$. For implementation purposes, it is more
convenient to express the range of application in terms of the projection
thresholds. To this end, $h$ can be replaced from Eq. (7) and $L$ from (11),
which leads to conditions depending only on $\eta_{i}$ and
${\eta_{\mathrm{ero}}}$, as follows:

$\begin{matrix}L>h&\implies\eta_{i}<2\eta_{\mathrm{ero}}-2+2\sqrt{1-{\eta_{\mathrm{ero}}}}\\\\[4.30554pt]
r_{\mathrm{fil}}\geq(L+h)/2&\implies\eta_{i}\geq
4-4\sqrt{1-{\eta_{\mathrm{ero}}}}-2{\eta_{\mathrm{ero}}}\end{matrix}$ (12)

By repeating the procedure from Eq. (9) to (12) for the 4 integration
conditions shown in Fig. 3, a set of equations is obtained which relate the
filter and projection parameters with the minimum size for any projection
threshold $\eta_{i}$, provided that $\eta_{i}<{\eta_{\mathrm{ero}}}$. The set
of equations are summarized in the four first rows of Table 1.

To obtain the relationships that define the minimum size of the void phase,
the same procedure must be used as for the solid phase, however, now starting
from a one-dimensional design domain containing a cavity of size $h$. To avoid
overextending the document with redundant information, this section is limited
to presenting the final equations that define the minimum size of the void
phase. The expressions are summarized in the last 4 rows of Table 1.

Having delivered the set of equations that expand the scope of the method
proposed by Qian and Sigmund (2013), the following section presents a
methodology to use these equations.

##  4 Imposing the desired minimum length scale

In structural design, the minimum length scale control is usually desired
because of design requirements or manufacturing limitations, hence in most
cases, the minimum size of the solid and void phases are known values
established for the intermediate design. Therefore, for the set of equations
presented in Table 1, the radii $r_{\mathrm{min.Solid}}^{\mathrm{int}}$ and
$r_{\mathrm{min.Void}}^{\mathrm{int}}$ are assumed user-defined input values.
In this case, the projection threshold $\eta_{i}$ corresponds to the
projection threshold ${\eta_{\mathrm{int}}}$, hence the desired length scale
for the intermediate design is a function of the projection thresholds and of
the size of the filter, namely,
$r_{\mathrm{min.Solid}}^{\mathrm{int}}({\eta_{\mathrm{int}}},{\eta_{\mathrm{ero}}},r_{\mathrm{fil}})$
and
$r_{\mathrm{min.Void}}^{\mathrm{int}}({\eta_{\mathrm{int}}},{\eta_{\mathrm{dil}}},r_{\mathrm{fil}})$.
Given the number of unknowns
$(r_{\mathrm{fil}},{\eta_{\mathrm{ero}}},\eta_{\mathrm{int}},\eta_{\mathrm{dil}})$,
the system of equations in Table 1 becomes indeterminate and the desired
minimum length scale can be imposed through multiple combinations of
parameters. Nevertheless, such freedom of parameters selection can be reduced
by considering the following three recommendations.

Firstly, a number of advantages have been observed when defining the
intermediate design with a threshold ${\eta_{\mathrm{int}}}=0.5$. For
instance, the projection features lower amounts of intermediate densities
compared to those projections that use a threshold other than 0.5 (Xu et al,
2010; Wang et al, 2011; da Silva et al, 2019), and a threshold
${\eta_{\mathrm{int}}}$ set to 0.5 provides the same size ranges (0.5) for the
erosion and dilatation thresholds, which is convenient for reducing rounding
errors, since small differences in projection thresholds could be insensitive
to the minimum size when using a coarse discretization of the design domain
(Qian and Sigmund, 2013). Secondly, for a particular combination of
thresholds, the filter size ($r_{\mathrm{fil}}$) can become considerably
larger than the desired minimum length scale, which could significantly
increase computational requirements (Lazarov and Sigmund, 2011). This can be
seen in Figs. 4 and 6(b). These figures show graphs that relate the filter
size ($r_{\mathrm{fil}}$) to the minimum size of solid or void phase and to
the erosion or dilation threshold. These graphs show that the closer
${\eta_{\mathrm{ero}}}$ and ${\eta_{\mathrm{dil}}}$ are to
${\eta_{\mathrm{int}}}$, the larger the filter radius, which inevitably
increases computational requirements. Under this observation, we recommend
choosing ${\eta_{\mathrm{ero}}}\geq 0.75$ and ${\eta_{\mathrm{dil}}}\leq
0.25$, thus it is ensured that $r_{\mathrm{fil}}\leq
2r_{\mathrm{min.Solid}}^{\mathrm{int}}$ and $r_{\mathrm{fil}}\leq
2r_{\mathrm{min.Void}}^{\mathrm{int}}$. Thirdly, erosion and dilation
thresholds too distant or too close to the intermediate threshold increases
oscillations of design variables during the optimization process. In general,
a good compromise is to choose $0.10\leq{\eta_{\mathrm{dil}}}\leq 0.4$ and
$0.60\leq{\eta_{\mathrm{ero}}}\leq 0.9$.

Figure 4:  Graphical relationship between the minimum size of the void phase,
the filter radius and the dilation projection.

The first observation removes an unknown from the system of equations, since
${\eta_{\mathrm{int}}}$ is set to 0.5, while the second and third observations
limit the range of the erosion and dilation thresholds. Thus, for a user-
defined minimum length scale, it is possible to develop an algorithm that
solves the system of equations considering the three observations. Here we
propose an algorithm based on graphic relationships, so that the reader can
easily find the desired parameters without the need to resort to a
computational algorithm. Nonetheless, we also provide as supplementary
material a code written in MATLAB named SizeSolution.m that performs the
procedure described below. It is important to note that the above observations
are based on numerical tests considered for specific optimization problems
formulated in the density approach. Therefore, it is possible that under other
topology optimization approaches or formulations the above observations are no
longer valid. However, the proposed procedure can be applied for any other
value of ${\eta_{\mathrm{int}}}$, or any other combination of parameters that
the user may consider convenient.

Given that the ranges of application of the equations in Table 1 are defined
as a function of the projection thresholds ($\eta_{i}$,
${\eta_{\mathrm{ero}}}$ and ${\eta_{\mathrm{dil}}}$), it is rather simple to
construct graphs with respect to them. For instance, as shown in Fig. 5, to
find the range of application of the equations defining the minimum size of
the solid phase, it is simply necessary to know in which region of Fig. 5 the
projection $\eta_{i}$ falls.

The first graph proposed in this work is shown in Fig. 6(a) and gathers 4
parameters, the minimum length scale ($r_{\mathrm{min.Solid}}^{\mathrm{int}}$
and $r_{\mathrm{min.Void}}^{\mathrm{int}}$) and the projection thresholds
(${\eta_{\mathrm{ero}}}$ and ${\eta_{\mathrm{dil}}}$), provided that
${\eta_{\mathrm{int}}}=0.5$. In this graph, the user can easily find the set
of erosion and dilation threshold that leads to the desired length scale.
Then, the user can access the graph in Fig. 6(b) to obtain the filter radius.
For example, for the following minimum length scale,
$r_{\mathrm{min.Solid}}^{\mathrm{int}}=3$ elements and
$r_{\mathrm{min.Solid}}^{\mathrm{int}}=3$ elements, the graph in Fig. 6(a) is
accessed with a value of 1.0 for the ordinate. According to the aforementioned
observations, the combination of thresholds [${\eta_{\mathrm{ero}}}$,
${\eta_{\mathrm{dil}}}$] that can be chosen among others are:

$\displaystyle[{\eta_{\mathrm{ero}}},{\eta_{\mathrm{dil}}}]$
$\displaystyle=[0.75,\;0.25]\;,$ (13a)
$\displaystyle[{\eta_{\mathrm{ero}}},{\eta_{\mathrm{dil}}}]$
$\displaystyle=[0.80,\;0.20]\;,$ (13b)
$\displaystyle[{\eta_{\mathrm{ero}}},{\eta_{\mathrm{dil}}}]$
$\displaystyle=[0.85,\;0.15]\;,$ (13c)
$\displaystyle[{\eta_{\mathrm{ero}}},{\eta_{\mathrm{dil}}}]$
$\displaystyle=[0.90,\;0.10]\;.$ (13d)

Arbitrarily, [0.75, 0.25] is selected, and from the graph in Fig. 6(b), it is
obtained that $r_{\mathrm{fil}}=2r_{\mathrm{min.Solid}}^{\mathrm{int}}=6$
elements.

It should be noted that the graph in Fig. 6(a) is constructed considering a
resolution of 0.05 in the projection thresholds. This is due to the fact that
the discretization of the filter radius in topology optimization is generally
coarse, and decimal numbers smaller than 0.05 in the threshold value have
usually a negligible effect on the minimum length scale of the optimized
design. However, if a better resolution is required, the user can resort to
the attached code SizeSolution.m.

Figure 5: Graphic representation of the applicability of the set of equations
in Table 1 proposed for the solid phase. The denomination of the zone
corresponds to the row number of Table 1.

(a)

(b)

(c)

(d)

(a) (b)
(c) (d)

Figure 6: Graphical relationships between the minimum length scale and
$\eta_{\mathrm{ero}}$, $\eta_{\mathrm{int}}$, $\eta_{\mathrm{dil}}$, and
$r_{\mathrm{fil}}$. The minimum length scale is defined for the intermediate
design. The graphs are built using both, the analytical and the numerical
method. Graphs (a) and (b) are designed to obtain the projection thresholds
and the filter radius, respectively. Graphs (c) and (d) are designed to obtain
the dilation and erosion distances, respectively.

As previously mentioned, the erosion and dilation distances could be required,
for instance, when implementing maximum size restrictions (Fernández et al,
2020). These distances can be easily obtained from the equations of Table 1.
As shown in Fig. 2(d), the dilation distance, denoted by $t_{\mathrm{dil}}$,
is the offset distance between the intermediate and dilated designs.
Analogously, the erosion distance, denoted by $t_{\mathrm{ero}}$, is the
offset distance between the intermediate and eroded designs. Therefore:

$\begin{matrix}t_{\mathrm{dil}}&=r_{\mathrm{min.Solid}}^{\mathrm{dil}}{(\eta_{\mathrm{dil}},\eta_{\mathrm{ero}},r_{\mathrm{fil}})}-r_{\mathrm{min.Solid}}^{\mathrm{int}}{(\eta_{\mathrm{int}},\eta_{\mathrm{ero}},r_{\mathrm{fil}})}\vspace{2mm}\\\
t_{\mathrm{ero}}&=r_{\mathrm{min.Void}}^{\mathrm{ero}}{(\eta_{\mathrm{ero}},\eta_{\mathrm{dil}},r_{\mathrm{fil}})}-r_{\mathrm{min.Void}}^{\mathrm{int}}{(\eta_{\mathrm{int}},\eta_{\mathrm{dil}},r_{\mathrm{fil}})}\end{matrix}$
(14)

where the minimum size of the solid phase in the dilated design
($r_{\mathrm{min.Solid}}^{\mathrm{dil}}$) is obtained by choosing
$\eta_{i}={\eta_{\mathrm{dil}}}$ in the equations of Table 1. Analogously, the
minimum size of the void phase in the eroded design
($r_{\mathrm{min.Void}}^{\mathrm{ero}}$) is obtained by choosing
$\eta_{i}={\eta_{\mathrm{ero}}}$. In this way, analytical expressions can be
obtained for the dilation and erosion distances.

It is noted from Eq. (14) that the erosion and dilation depend on the 3
projection thresholds and on the size of the filter, namely,
$t_{\mathrm{ero}}(r_{\mathrm{fil}},{\eta_{\mathrm{ero}}},{\eta_{\mathrm{int}}},{\eta_{\mathrm{dil}}})$
and
$t_{\mathrm{dil}}(r_{\mathrm{fil}},{\eta_{\mathrm{ero}}},{\eta_{\mathrm{int}}},{\eta_{\mathrm{dil}}})$.
Therefore, if for some particular reason the user needs to impose the minimum
length scale, and the erosion and dilation distances, then the following
system of equations must be solved:

$\begin{split}r_{\mathrm{min.Solid}}^{\mathrm{int}}&=r_{\mathrm{min.Solid}}^{\mathrm{int}}(r_{\mathrm{fil}},{\eta_{\mathrm{int}}},{\eta_{\mathrm{ero}}})\vspace{1mm}\\\
r_{\mathrm{min.Void}}^{\mathrm{int}}&=r_{\mathrm{min.Solid}}^{\mathrm{int}}(r_{\mathrm{fil}},{\eta_{\mathrm{int}}},{\eta_{\mathrm{dil}}})\vspace{2mm}\\\
t_{\mathrm{ero}}&=r_{\mathrm{min.Solid}}^{\mathrm{int}}(r_{\mathrm{fil}},{\eta_{\mathrm{ero}}},{\eta_{\mathrm{int}}},{\eta_{\mathrm{dil}}})\vspace{1mm}\\\
t_{\mathrm{dil}}&=r_{\mathrm{min.Solid}}^{\mathrm{int}}(r_{\mathrm{fil}},{\eta_{\mathrm{ero}}},{\eta_{\mathrm{int}}},{\eta_{\mathrm{dil}}})\\\
\end{split}$ (15)

The system in Eq. (15) is determined and there is only one combination of
parameters [$r_{\mathrm{fil}}$, ${\eta_{\mathrm{ero}}}$,
${\eta_{\mathrm{int}}}$, ${\eta_{\mathrm{dil}}}$] that leads to the desired
length scale [$r_{\mathrm{min.Solid}}^{\mathrm{int}}$,
$r_{\mathrm{min.Void}}^{\mathrm{int}}$, $t_{\mathrm{ero}}$,
$t_{\mathrm{dil}}$]. This can be illustrated with the following example. For
each combination of thresholds [$\eta_{\mathrm{ero}}$, $\eta_{\mathrm{dil}}$]
given in Eq. (13), the erosion and dilation distances [${t_{\mathrm{ero}}}$ ,
${t_{\mathrm{dil}}}$ ] are provided, as follows:

$\displaystyle[{\eta_{\mathrm{ero}}},{\eta_{\mathrm{dil}}}]$
$\displaystyle=[0.75,\;0.25]\;,[{t_{\mathrm{ero}}},{t_{\mathrm{dil}}}]=[1.76,\;1.76]$
(16a) $\displaystyle[{\eta_{\mathrm{ero}}},{\eta_{\mathrm{dil}}}]$
$\displaystyle=[0.80,\;0.20]\;,[{t_{\mathrm{ero}}},{t_{\mathrm{dil}}}]=[1.99,\;1.99]$
(16b) $\displaystyle[{\eta_{\mathrm{ero}}},{\eta_{\mathrm{dil}}}]$
$\displaystyle=[0.85,\;0.15]\;,[{t_{\mathrm{ero}}},{t_{\mathrm{dil}}}]=[2.21,\;2.21]$
(16c) $\displaystyle[{\eta_{\mathrm{ero}}},{\eta_{\mathrm{dil}}}]$
$\displaystyle=[0.90,\;0.10]\;,[{t_{\mathrm{ero}}},{t_{\mathrm{dil}}}]=[2.43,\;2.43]$
(16d)

It is recalled that all combinations of thresholds [$\eta_{\mathrm{ero}}$,
$\eta_{\mathrm{dil}}$] in Eq. (16) lead to the minimum length scale
$r_{\mathrm{min.Solid}}^{\mathrm{int}}=r_{\mathrm{min.Void}}^{\mathrm{int}}=3$
elements. However, they all result in different erosion and dilation
distances.

To the authors’ knowledge, the need to impose specific values for the erosion
and dilation distances has not yet been claimed, so this manuscript is limited
to providing these values for a given combination of parameters. The proposed
graphs are shown in Fig. 6(c) and 6(d), which depend on the erosion and
dilation thresholds and on the corresponding minimum size. For instance,
considering $[{\eta_{\mathrm{ero}}},{\eta_{\mathrm{dil}}}]$ = [0.75, 0.25],
the erosion and dilation distances are $t_{\mathrm{ero}}=0.58\times 3$
elements and $t_{\mathrm{dil}}=0.58\times 3$ elements, where 0.58 is obtained
from the graphs.

##  5 Sources of error

The analytical procedure developed in this article allows to quickly obtain
the set of filtering and projection parameters that imposes the desired
minimum length scale. However, in practice, the length scale of the optimized
design often differs from the desired values (Qian and Sigmund, 2013). This is
mainly due to the fact that the analytical method assumes (i) a continuous
domain, (ii) a perfect Heaviside projection ($\beta\to\infty$), and (iii) that
the eroded and dilated fields project an infinitesimal minimum size in the
solid and void phases, respectively. The assumptions (i), (ii), and (iii) are
not met in topology optimization mainly because the design domains are
discretized into finite elements and because the Heaviside projection is
smoothed. To assess the error introduced by these assumptions, in the
following we compare the analytical method with the numerical method proposed
by Wang et al (2011), which considers a discrete domain and a smoothed
Heaviside function.

The procedure for obtaining the minimum size using the numerical method is
analogous to the analytical method but now using a one-dimensional domain
discretized into $N$ elements. For example, to obtain the minimum size of the
solid phase the three-field scheme is applied to a one-dimensional design
domain $\bm{\rho}$ containing solid elements in a length $h$. Then, the size
of the eroded, intermediate and dilated fields are measured, and the length
$h$ is adjusted so that the resulting eroded field has an infinitesimal size
(1 solid element). This process is repeated several times for different values
of ${\eta_{\mathrm{ero}}}$, and the resulting minimum size is normalized with
respect to the size of the chosen filter. For implementation details regarding
the numerical method, the reader is referred to the works of Wang et al (2011)
and Fernández et al (2020), and to the attached MATLAB code named
NumericalSolution.m.

To validate the analytical and numerical methods, the latter is implemented
using a Heaviside function at $\beta=500$, a filter radius of 1000 elements
and a design domain discretized into 10 thousand elements, so the numerical
method can be considered continuous and under similar assumptions than the
analytical one. The relationships obtained with the numerical method can be
seen in the graphs of Fig. 6. The agreement of results from both methods
allows us to validate the set of equations provided in Table 1.

The three sources of error that affect the analytical method are discussed
below.

###  5.1 Continuous design domain

When the design domain is discretized using finite elements, the radii that
define the minimum sizes ($r_{\mathrm{fil}}$,
$r_{\mathrm{min.Solid}}^{\mathrm{int}}$,
$r_{\mathrm{min.Solid}}^{\mathrm{dil}}$,
$r_{\mathrm{min.Void}}^{\mathrm{int}}$ and
$r_{\mathrm{min.Void}}^{\mathrm{ero}}$) are defined by a discrete number of
elements. In this case, the rounding error is $\pm$ 1 finite element in the
radius that defines the minimum size. Therefore, to reduce this error, it is
sufficient to reduce the size of the elements by mesh refinement. To
illustrate this remark graphically, we consider the following two
discretizations of a one-dimensional domain, one containing 100 elements and
the other containing 200 elements. In our implementation, the size of the
filter is chosen equal to $10\%$ of the domain size. Therefore, for the
discretizations containing 100 and 200 elements, the filter radius contains 10
and 20 elements, respectively. For each discretization, the relationship
between the minimum size of the solid phase, the size of the filter and the
erosion threshold is plotted in Fig. 7, using both the analytical and the
numerical method.

As the analytical method does not depend on the discretization, it provides
identical relationships on both discretizations. However, the rounding error
associated to the analytical method does depend on the discretization and can
be plotted as an offset of the analytical curve, as shown by the dashed curves
in Fig. 7. The rounding error, denoted as $\delta r$, corresponds to 1 element
in the radius, i.e. 2 elements in the diameter. That is, the vertical offset
of the analytical curve due to the rounding error is equal to $2\delta r=2/10$
and to $2\delta r=2/20$ in Figs. 7(a) and 7(b), respectively.

(a) $N=100$, $r_{\mathrm{fil}}=10$.

(b) $N=200$, $r_{\mathrm{fil}}=20$

Figure 7: Relationship between the minimum size of the solid phase, the eroded
threshold, and the filter radius. The dotted line represents the analytical
relationship including the rounding error.

The error curves agree with the results obtained from the numerical method,
which shows that the error of the analytical method coming from the assumption
of a continuous domain can be easily estimated. This estimation allows to know
either the margin of error of the minimum size ($\delta r$) for a given
analytical threshold (${\eta_{\mathrm{ero}}}$ or ${\eta_{\mathrm{dil}}}$), or
the margin of error of the analytical threshold ($\delta\eta$) for a given
minimum size, as shown in Fig. 7(a).

If the user of the robust approach needs to impose precisely the minimum size,
he can choose the projection thresholds (${\eta_{\mathrm{ero}}}$ and
${\eta_{\mathrm{dil}}}$) that lead to a larger filter radius, or he can use
the numerical method (Wang et al, 2011) using a discretization representative
of the design domain to be optimized.

###  5.2 Perfect Heaviside projection

The analytical method developed in this work assumes a perfect Heaviside
function, which results in a projected field $\bm{\bar{\rho}}$ containing
discrete densities (0 and 1). However, in topology optimization, the Heaviside
function is smoothed resulting in regions of intermediate densities that lie
on the surface of the optimized structure. For this reason, in practice, once
the optimized solution ${{\bm{\bar{\rho}}}_{\mathrm{int}}}$ is obtained, a
cut-off value $\varepsilon$ is usually added to the projected densities in
order to get the optimized structure. Therefore, the final design intended for
manufacturing is obtained by projecting the design
${{\bm{\bar{\rho}}}_{\mathrm{int}}}$ as follows:

${\bar{\rho}^{\varepsilon}}_{i}=\left\\{\begin{matrix}1\;\;,&\text{if}\;\;\;\;\bar{\rho}_{\mathrm{int}(i)}\geq\varepsilon\vspace{1mm}\\\
0\;\;,&\text{otherwise}\end{matrix}\right.$ (17)

Since the post-processed design $\bm{\bar{\rho}^{\varepsilon}}$ is the one
intended for manufacturing, the analytical expressions relating the minimum
length scale and the parameters that define it must be elaborated for the
field $\bm{\bar{\rho}^{\varepsilon}}$. As will be explained hereafter, the
equations provided in Table 1 can be easily adapted to consider the cut-off
value $\varepsilon$.

Figure 8:  Illustration of a transition region where projected densities reach
intermediate values.

Consider Fig. 8 illustrating a transition zone between the solid structure and
the void region. In blue it is shown the filtered field $\bm{\tilde{\rho}}$
and in red the projected field $\bm{\bar{\rho}}$ obtained with $\beta=32$. In
this illustration, the density cut-off is defined as $\varepsilon=0.95$. A
perfect Heaviside projection ($\beta\to\infty$) would define the size of the
solid zone at the coordinate $x_{i}$. However, the smoothed Heaviside
projection defines the solid/void transition in the coordinate
$x_{\varepsilon}$. The idea then is to find the projection threshold
$\eta_{\varepsilon}$ for which a perfect Heaviside projection produces the
same solid/void transition coordinate than the one that is obtained with a
smoothed Heaviside projection with a threshold $\eta_{i}$. Thus, the value of
the filtered density present at the coordinate $x_{\varepsilon}$ must be
found. To this end, from Fig. 8 it is observed that:

$\varepsilon=\frac{\tanh{(\beta\eta_{i})}\tanh{(\beta(\eta_{\varepsilon}-\eta_{i}))}}{\tanh{(\beta\eta_{i})}\tanh{(\beta(1-\eta_{i}))}}$
(18)

Assuming $\beta>10$ at the end of the optimization process, which is a common
practice when dealing with Heaviside projection, Eq. (18) can be simplified
to:

$\tanh{(\beta(\eta_{\varepsilon}-\eta_{i}))}=2\varepsilon-1$ (19)

From Eq. (19), the filtered density $\eta_{\varepsilon}$ can be obtained:

$\eta_{\varepsilon}=\eta_{i}+\frac{1}{\beta}\mathrm{atanh}(2\varepsilon-1)$
(20)

We recall that the filtered density $\eta_{\varepsilon}$ represents the value
for which a perfect Heaviside function with threshold $\eta_{\varepsilon}$ and
a smoothed Heaviside function with threshold $\eta_{i}$ result in the same
size for the solid and void phases. Therefore, to take into account the
intermediate densities resulting from a smoothed Heaviside function, the
shifting term $\mathrm{atanh}(2\varepsilon-1)/\beta$ must be added to the
thresholds involved in the analytical equations (Table 1). Note that if
$\beta\to\infty$ or $\varepsilon=0.50$, then the shifting term is zero and
$\eta_{\varepsilon}=\eta_{i}$.

(a)

(b)

Figure 9: (a) Minimum length scale obtained with the analytical method
assuming a perfect Heaviside projection ($\beta\to\infty$), and with the
numerical method using a smoothed Heaviside projection ($\beta=32$) and three
different cut-off values $\varepsilon$. (b) The analytical curves obtained by
adding the shifting term $\mathrm{atanh}(2\varepsilon-1)/\beta$ to the
thresholds.

To assess the error introduced by a smoothed Heaviside projection, the minimum
length scale obtained with the analytical and the numerical methods are
compared. The analytical method is applied with $\beta\to\infty$, so it does
not consider the shifting term on the thresholds. For the numerical method we
use $\beta=30$, a filter radius equal to 200 and a one-dimensional domain
discretized into 2000 elements in order to avoid rounding errors and isolate
the effect of using a smoothed Heaviside projection. In the numerical method,
three cut-off values are used, $\varepsilon=0.01$, $\varepsilon=0.50$ and
$\varepsilon=0.99$. The results are shown in Fig. 9(a).

The results show that the influence of using a smoothed Heaviside function is
observed when $\varepsilon\neq 0.5$ and when ${\eta_{\mathrm{ero}}}>0.75$. For
this reason, we recommend to use a density cut-off value equal to 0.5, as this
avoids the need to add the shifting term to the projection thresholds.
Nevertheless, the user can easily adjust the analytical curves for
$\varepsilon\neq 0.5$, since it is only required to add
$\mathrm{atanh}(2\varepsilon-1)/\beta$ to the projection thresholds, as shown
in the corrected curves in Fig. 9(b).

It is worth mentioning that probably the effect of using a smoothed Heaviside
function is only seen for ${\eta_{\mathrm{ero}}}\geq 0.75$ because the
projected field becomes more discrete as the projection threshold is closer to
0.5 (Wang et al, 2011).

###  5.3 Infinitesimal size

It is recalled that the analytical and numerical methods are developed under
the condition of robustness which ensures that the manufactured design will
feature a good performance even if the blueprint design is uniformly thinned
or thickened during the manufacturing process. This condition dictates that
the structural members/cavities must be projected by at least one solid/void
element in the case of an eroded/dilated design.

To obtain the minimum size, the analytical method assumes that the size of the
structural members/cavities is infinitesimal ($\approx 0$) in the
eroded/dilated design, but in practice this does not occur. In topology
optimization, it is observed that the smallest structural members/cavities in
the eroded/dilated design are composed of one or a few elements, generally
described with intermediate densities. Thus, in practice, the minimum size of
the solid phase in the eroded design ($r_{\mathrm{min.Solid}}^{\mathrm{ero}}$)
results in a discrete number of elements. This can be seen in Fig. 10, which
is build using the numerical method. There, a filtered field and its eroded
projection under the robust condition are shown. It can be seen that the
minimum size of the eroded design is not infinitesimal, and therefore the
minimum size of the solid phase in the intermediate and dilated design would
be larger than the value predicted by the analytical method.

Figure 10: Illustration in which an infinitesimal size is not reached for the
eroded projection.

The minimum size of the solid and void phases in the eroded and dilated
designs ($r_{\mathrm{min.Solid}}^{\mathrm{ero}}$ and
$r_{\mathrm{min.Void}}^{\mathrm{dil}}$) are defined by the size of the
elements that discretize the design space. In addition, the amount of elements
with intermediate densities defining the minimum size depends on the steepness
of the smoothed Heaviside function ($\beta$), therefore, the error introduced
by assuming an infinitesimal size in the robust condition is related to the
two sources of error mentioned previously. To isolate the effect of the
infinitesimal size in the robust condition and illustrate the error introduced
by the assumption of infinitesimal size, we make use of the numerical method.
The numerical method is implemented in a discretized domain containing a large
number of elements ($10^{4}$ elements) and using a large filter size ($10^{3}$
elements) to simulate a continuous domain. The steepness parameter of the
smoothed Heaviside projection is set as $\beta=512$. Thus, the size and
density of the elements are excluded from the analysis. The effect of not
achieving an infinitesimal size in the condition of robustness is
intentionally imposed in the numerical code. For this, the robustness
condition is considered satisfied if
$r_{\mathrm{min.Solid}}^{\mathrm{ero}}=\alpha\>r_{\mathrm{fil}}$. Considering
that $r_{\mathrm{min.Solid}}^{\mathrm{ero}}$ represents one finite element in
topology optimization, and that representative values for $r_{\mathrm{fil}}$
are between 2 and 10 elements, it is reasonable to consider values for
$\alpha$ ranging from 0.1 to 0.5. Taking into consideration the above, we
build the graphical solutions that relate the minimum size in the solid phase
($r_{\mathrm{min.Solid}}^{\mathrm{int}}$), the filter radius and the erosion
threshold. The graph is shown in Fig. 11.

Figure 11: Effect of not reaching an infinitesimal size in the condition of
robustness $r_{\mathrm{min.Solid}}^{\mathrm{ero}}=0$. This condition is
assessed by imposing
$r_{\mathrm{min.Solid}}^{\mathrm{ero}}=\alpha\>r_{\mathrm{fil}}$.

The graph shows that the error of not achieving an infinitesimal size in the
eroded design produces a minimum size of the solid phase bigger than the value
predicted by the analytical method. The error related to the infinitesimal
size is low in comparison to the rounding error, so the latter would be the
most relevant source of error from the analytical method that assumes a
continuous design domain.

In summary, this section discussed the scope of the analytical method (Qian
and Sigmund, 2013) by comparing it with the numerical one (Wang et al, 2011),
both developed for a one-dimensional design domain. To this end, different
sources of error were examined. In general, the errors can be controlled by
mesh refinement or by correcting the projection thresholds according to the
cut-off $\varepsilon$ value. In the following section, the analytical method
is assessed using 2D-topology optimization problems.

##  6 Numerical examples and discussion

This section examines the reliability of the analytical expressions provided
in Table 1. To this end, a set of 2D topology optimization problems are
solved, from which the length scale is measured graphically and compared with
the imposed values. Then, some designs obtained with maximum size constraints
are provided to illustrate the use of the erosion and dilation distances.
Finally, this section provides a remark regarding the simplified robust
formulation, where the intermediate and dilated designs are removed from the
objective function.

###  6.1 Minimum length scale

The minimum length scale is assessed using the heat exchanger design problem
described in Section 2. A set of results is generated from this design
problem, which differ in the desired minimum length scale and in the
discretization used. Specifically, three sets of solutions are obtained by
discretizing the design domain into $100\times 100$, $200\times 200$ and
$400\times 400$ quadrilateral elements. In addition, three different length
scales are prescribed for each discretization, which are reported as the ratio
between the minimum size of the solid phase and the minimum size of the void
phase, i.e.
$r_{\mathrm{min.Solid}}^{\mathrm{int}}/r_{\mathrm{min.Volid}}^{\mathrm{int}}$.
The chosen ratios are $1/2$, $1/1$ and $2/1$. The minimum size in the solid
phase is the same in all scenarios and is defined as a physical dimension. In
number of finite elements, the radius that defines the minimum size of the
solid phase ($r_{\mathrm{min}}^{\mathrm{int}}$) is equal to 1, 2 and 4, for
the discretizations that use $100^{2}$, $200^{2}$ and $400^{2}$ elements,
respectively. It is well known that the initial values of design variables
have a huge influence on the resulting topology when it comes to thermal
compliance minimization (Yan et al, 2018), hence, in order to facilitate the
comparison of results, we impose a base structure as a starting point, which
is shown in Fig. 12. Before presenting the results, the procedure to obtain
the minimum length scale from the optimized designs is detailed.

Figure 12: Initial distribution of design variables considered for the thermal
compliance minimization problem.

The minimum size of the solid phase is measured graphically by counting the
number of finite elements that define the size of the thinnest structural
branch. Similarly, the minimum size of the void phase is measured by counting
the elements in the radius of the largest circumference that can be inscribed
at the re-entrant corners of the design. For example, consider the design of
Fig. 13 where the minimum length scales
$r_{\mathrm{min.Solid}}^{\mathrm{int}}=\mathrm{min.Void}^{\mathrm{int}}=2$
elements are imposed. To determine the real minimum size of the void phase,
the largest circle in Fig. 13 that falls into the re-entrant corners of the
design is identified. The corners analyzed are those that form a sharp angle
between two structural branches. Fig. 13 shows three representative re-entrant
corners of the design depicted in Fig. 13. From there it is observed that the
minimum size is given by a circle of radius 2 elements (zone C). Similarly,
the largest region that fits into the thinnest structural members is
determined, as shown in Fig. 13. There, the minimum size of the solid phase is
given by a circle of radius 1.5 elements (zone D).

After describing the test case, the obtained results are presented. Table 2
contains the nine results generated in this example (3 length scales $\times$
3 discretizations). The imposed minimum length scales are reported graphically
next to each solution. The minimum size of the void phase is indicated in
blue, while the minimum size of the solid phase in magenta. The table also
reports the 3 parameters required to impose the desired minimum length scales,
i.e. $\eta_{\mathrm{ero}}$, $\eta_{\mathrm{dil}}$, and $r_{\mathrm{fil}}$.
These parameters have been obtained using the analytical method implemented in
the MATLAB code provided with this paper (SizeSolution.m). It is recalled that
all the examples assume $\eta_{\mathrm{int}}=0.5$.

(a)

(b)

(c)

(d)

(a)
(b)
(c)
(d)

Figure 13: Illustration of the minimum length scale measurement. In (a) the test regions and their radius sizes given in number of finite elements. In (b) the minimum size of the void phase, (c) the optimized heat exchanger, and (d) the minimum size of the solid phase. The imposed length scales are graphically shown in (c), at the upper left corner of the design. Mesh | $r_{\mathrm{min.Solid}}^{\mathrm{int}}/r_{\mathrm{min.Void}}^{\mathrm{int}}$
---|---
$1/2$ | $1/1$ | $2/1$
using [$\eta_{\mathrm{ero}}$, $\eta_{\mathrm{dil}}$] = [0.65, 0.05] | using [$\eta_{\mathrm{ero}}$, $\eta_{\mathrm{dil}}$] = [0.70, 0.30] | using [$\eta_{\mathrm{ero}}$, $\eta_{\mathrm{dil}}$] = [0.80, 0.42]
$100\times 100$ |  |  |
$r_{\mathrm{min.Solid}}^{\mathrm{int}}=1$, $r_{\mathrm{fil}}=2.58$. |  $r_{\mathrm{min.Solid}}^{\mathrm{int}}=1$, $r_{\mathrm{fil}}=2.24$. |  $r_{\mathrm{min.Solid}}^{\mathrm{int}}=1$, $r_{\mathrm{fil}}=1.81$.
$200\times 200$ |  |  |
$r_{\mathrm{min.Solid}}^{\mathrm{int}}=2$, $r_{\mathrm{fil}}=5.16$. |  $r_{\mathrm{min.Solid}}^{\mathrm{int}}=2$, $r_{\mathrm{fil}}=4.47$. |  $r_{\mathrm{min.Solid}}^{\mathrm{int}}=2$, $r_{\mathrm{fil}}=3.62$.
$400\times 400$ |  |  |
$r_{\mathrm{min.Solid}}^{\mathrm{int}}=4$, $r_{\mathrm{fil}}=10.32$. |  $r_{\mathrm{min.Solid}}^{\mathrm{int}}=4$, $r_{\mathrm{fil}}=9.94$. |  $r_{\mathrm{min.Solid}}^{\mathrm{int}}=4$, $r_{\mathrm{fil}}=7.24$.
Table 2: Optimized designs for the heat exchanger problem. The imposed volume
constraint is $V^{*}_{\mathrm{min}}=20\%$.

The measured minimum size of the solid and void phases are reported on Figs.
14 and 14, respectively. To construct these graphs, the measured minimum sizes
are normalized with respect to the desired minimum sizes, and the values are
placed in the ordinate coordinate. This is done for each discretization (which
determines the abscissa) and for each length scale. This procedure is also
carried out with results obtained from the numerical method (Wang et al,
2011), which is executed using representative discretizations for each case.
The values obtained from the numerical method are labeled as 1D in the graphs
of Fig. 14, while the values measured from the optimized designs are labeled
as 2D. The error bars in Fig. 14 illustrate the rounding error present in the
analytical method (see Section 5.1 for definition). Rounding error lines have
not been plotted in Fig. 14 because the length scale ratios have different
rounding errors.

As a general observation, we can point out from Fig. 14 that mesh refinement
reduces the error between the measured minimum size and the desired minimum
size (which is imposed through the set of parameters provided by the
analytical method). This is consistent with the observations made in the
previous section, where a continuous and uni-dimensional design domain was
used. In addition, the predicted error for the analytical method (the error
bars) are also consistent with the measured sizes, which validates the scope
and limitations of the analytical equations discussed in Section 5.

Regarding the graph in Fig. 14, it can be mentioned that in the nine 2D
designs, the error between the measured and desired minimum sizes is half a
finite element. This error is relatively large in the coarse discretization
(50$\%$ error), and small in the fine discretization (12.5$\%$ error).
However, despite the fact that the analytical method is not exact, it seems to
be accurate enough in the examples examined, since in all cases the error is
the same, half a finite element. It is interesting to note that the numerical
method that assumes a discrete 1D domain, estimates a minimum size that
differs by half a finite element with respect to the imposed value, which
matches the measured error. However, the numerical method provides a minimum
size of the solid phase larger than the measured one. This is probably due to
the fact that the condition of robustness refers to one finite element and not
to an infinitesimal size, as discussed in Section 5.3.

On the other hand, the graph built for the void phase (Fig. 14) shows a
different pattern. The measured minimum size of the void phase is often equal
to the desired one and even bigger for some designs. This could be explained
by the simple fact that the chosen test case does not present a geometric
singularity over the minimum size of the void phase, therefore it is not
possible to guarantee that the smallest re-entrant corner of the design will
indeed correspond to the minimum size imposed by the robust formulation. The
representative case might be that where $r_{\mathrm{min.Void}}^{\mathrm{int}}$
is chosen twice the size $r_{\mathrm{min.Solid}}^{\mathrm{int}}$ (ratio 1/2).
In such a case, the measured error corresponds to half a finite element
smaller than the desired value, that is, the same result as for the solid
phase.

(a)

(b)

(a)
(b)

Figure 14: The graphs summarize the measured minimum size of (a) the solid
phase and (b) the void phase from the designs in Table 2. For each
discretization, two set of results are provided, labeled as $2D$ and $1D$. The
$2D$ represents the graphical measurement normalized with respect to the
imposed value. The $1D$ represents the minimum size obtained from the
numerical method normalized with respect to the imposed value.

###  6.2 Erosion and dilation distances

In the following example we illustrate the use of the erosion and dilation
distances that are provided for a desired minimum length scale. To this end,
we use maximum size constraints, where the erosion and dilation distances are
essential information to impose a consistent length scale in the robust
formulation (Fernández et al, 2020). For the sake of completeness of the
manuscript, the non-linear force inverter is considered in this illustrative
example. The topology optimization problem of the non-linear force inverter
including maximum size constraint reads as follow :

$\begin{split}\text{min}\quad&\text{max}\left(\mathrm{c}({{\bm{\bar{\rho}}}_{\mathrm{ero}}}),\;\mathrm{c}({{\bm{\bar{\rho}}}_{\mathrm{int}}})\;,\mathrm{c}({{\bm{\bar{\rho}}}_{\mathrm{dil}}})\right)\\\
\text{s.t.:}\quad&\mathbf{v}^{\intercal}\>{{\bm{\bar{\rho}}}_{\mathrm{dil}}}\leq
V^{*}_{\mathrm{dil}}(V^{*}_{\mathrm{int}})\\\
&\mathrm{G_{ms}}({\bm{\bar{\rho}}}_{\mathrm{dil}})\leq 0\\\ &0\leq\rho_{i}\leq
1\;\;,\;\;i=1,...,N\;\;,\end{split}$ (21)

where $\mathrm{G_{ms}}$ is the maximum size constraint defined exactly as in
(Fernández et al, 2020). Therefore, $\mathrm{G_{ms}}$ represents a p-mean
aggregation function that gathers local maximum size restrictions. As in the
cited work, the $p$ aggregation exponent is set at 100.

Regarding the optimization parameters, these represent the implementation of
Wang et al (2014), i.e. the Heaviside parameter $\beta$ is initialized at 1
and is increased by 1 every 20 iterations until a value $\beta=16$ is reached.
Then, 20 more iterations are carried out with a $\beta=32$. The SIMP penalty
parameter is set at 3. We have found that this setting of parameters works
well for introducing maximum size restrictions into the non-linear force
inverter formulated under the robust design approach.

To impose minimum and maximum length scales, the maximum size constraint
$\mathrm{G_{ms}}$ is applied on the dilated design (Fernández et al, 2020). To
do so, the regions where the local maximum size constraints are applied must
be scaled according to the dilation distance. For instance, if the desired
maximum size in the intermediate design is defined by a circle of radius
$r_{\mathrm{max.Solid}}^{\mathrm{int}}$, the maximum size constraint should be
formulated for the dilated design using a circle of radius:

$r_{\mathrm{max.Solid}}^{\mathrm{dil}}=r_{\mathrm{max.Solid}}^{\mathrm{int}}+t_{\mathrm{dil}}$
(22)

Equation (22) explains the need of knowing the dilation distance when imposing
maximum size restrictions in the robust formulation. As mentioned previously,
this information can be easily obtained from the graphs generated by the
analytical equations. For example, in the following, the half force inverter
depicted in Fig. 1 is solved. The design domain is discretized into $200\times
100$ quadrilateral finite elements. The minimum size of the solid phase is set
as $r_{\mathrm{min.Solid}}^{\mathrm{int}}=2$ elements, while the maximum size
is set as $r_{\mathrm{max.Solid}}^{\mathrm{int}}=3$ elements. Two different
values for the minimum size of the void phase are chosen,
$r_{\mathrm{min.Void}}^{\mathrm{int}}=2$ and
$r_{\mathrm{min.Void}}^{\mathrm{int}}=3$ elements. The filter and projection
parameters used to impose the desired length scale are reported in Table 3. In
all cases, the volume constraint is set to $25\%$.

$r_{\mathrm{min.Void}}^{\mathrm{int}}$ | $r_{\mathrm{fil}}$ | $\eta_{\mathrm{ero}}$ | $\eta_{\mathrm{dil}}$ | $t_{\mathrm{dil}}$ | $t_{\mathrm{ero}}$
---|---|---|---|---|---
$2$ | $4.47$ | $0.70$ | $0.30$ | $1.03$ | $1.03$
$3$ | $4.47$ | $0.70$ | $0.11$ | $2.41$ | $1.03$
Table 3: Parameters to impose $r_{\mathrm{min.Solid}}^{\mathrm{int}}=2$ on the
design with $\beta=32$.

To obtain the dilation distance $t_{\mathrm{dil}}$, the projection thresholds
that define the minimum length scale must be defined. These values can be
graphically obtained from Fig. 6, or from the attached MATLAB code
(SizeSolution.m). For example, using this code, the sets of parameters in
Table 3 are obtained. The dilation distance is also reported there. For the
minimum sizes $r_{\mathrm{min.Void}}^{\mathrm{int}}=2$ and
$r_{\mathrm{min.Void}}^{\mathrm{int}}=3$ elements, the dilation distances are
$t_{\mathrm{dil}}=1$ and $t_{\mathrm{dil}}=2$ elements, respectively (the
numbers have been rounded to the nearest integer).

The results are shown in Figs. 15 and 16. Each figure reports 2 optimized
designs, which are obtained with and without maximum size restrictions. To
facilitate the visual comparison between results, the designs are placed with
respect to the symmetry axis and are reported in their deformed configuration.
The imposed length scale is also reported graphically next to each design. As
in the previous examples, the blue and magenta circles represent the minimum
size of the void and solid phase, respectively. The black circle represents
the maximum size desired for the intermediate design (which is imposed through
the dilated design)111The maximum size is actually imposed using an annular
region (Fernández et al, 2020), but for illustrative purposes, a circular
region is drawn..

Figure 15: Nonlinear force inverter without (upper half) and with (lower half)
maximum size constraints. The minimum length scales are
$r_{\mathrm{min.Solid}}^{\mathrm{int}}$ =
$r_{\mathrm{min.Void}}^{\mathrm{int}}$ = 2 elements. The maximum size is
$r_{\mathrm{max.Solid}}^{\mathrm{int}}=3$ elements. Figure 16: Nonlinear force
inverter without (upper half) and with (lower half) maximum size constraints.
The minimum length scales are $r_{\mathrm{min.Solid}}^{\mathrm{int}}$ = 2
elements and $r_{\mathrm{min.Void}}^{\mathrm{int}}$ = 3 elements. The
prescribed maximum size is $r_{\mathrm{max.Solid}}^{\mathrm{int}}=3$ elements.

The results are consistent with the imposed length scales. In the 4 designs
reported in Figs. 15 and 16, the minimum size of the solid and void phases are
met with a half-finite element of error (which agrees with the analytical
rounding error). Regarding the maximum size, this also matches the imposed
value. However, due to the inherent drawbacks of the aggregation function and
the strong non-linearity of the force inverter design problem, there are local
regions where the imposed maximum size restriction is not met, such as the
horizontal bar in the design of Fig. 16.

This force inverter test case shows the usefulness of the proposed equations
and provided codes, since they allow to quickly obtain the filter and
projection parameters, and the dilation distance that allow to impose the
desired length scales.

###  6.3 Remark on the simplified robust formulation

Before concluding this manuscript we devote a discussion concerning the robust
formulation. It is widely known that the robust formulation of the
optimization problem can be simplified when the objective function is
monotonically dependent on the volume fraction. The most widespread case in
the literature is the compliance minimization problem, where the intermediate
and dilated fields have less compliance than the eroded design and therefore
can be removed from the objective function without compromising the robustness
of the formulation. In this case, the objective function is evaluated only for
the eroded design, with the intermediate and dilated fields remaining solely
for formulating design constraints. When the volume restriction is the only
constraint included in the optimization problem, the constraint can be
formulated using the intermediate design, assuming that this design field is
the one intended for manufacturing. However, it has been mentioned in the
literature that evaluating the volume restriction in the dilated design
promotes convergence to better optimums since numerical instabilities are
prevented. In this section we add another reason that has not been mentioned
so far (to the best of the authors’ knowledge). If the dilated design is not
included in the volume restriction, then it is not possible to ensure the
control over the minimum size of the void phase. To explain this statement,
consider the following two robust topology optimization formulations for the
thermal compliance minimization problem:

P.IP.II $\overbrace{\hskip 85.35826pt}\hskip 11.38109pt\overbrace{\hskip
105.27519pt}$ $\displaystyle{\min_{\bm{\rho}}}$ $\displaystyle\quad
c({{\bm{\bar{\rho}}}_{\mathrm{ero}}})$ $\displaystyle{\min_{\bm{\rho}}}$
$\displaystyle\quad c({{\bm{\bar{\rho}}}_{\mathrm{ero}}})$
$\displaystyle\mathrm{s.t.:}$
$\displaystyle\quad\mathbf{v}^{\intercal}{{\bm{\bar{\rho}}}_{\mathrm{int}}}\leq
V^{*}_{\mathrm{int}}\>,$ $\displaystyle\mathrm{s.t.:}$
$\displaystyle\quad\mathbf{v}^{\intercal}{{\bm{\bar{\rho}}}_{\mathrm{dil}}}\leq
V^{*}_{\mathrm{dil}}(V^{*}_{\mathrm{int}})$ (23) $\displaystyle\quad
0\leq{\rho_{i}}\leq 1$ $\displaystyle\quad 0\leq{\rho_{i}}\leq 1$

The two optimization problems, P.I and P.II in Eq. (6.3), are formulated under
the robust design approach, but P.I evaluates the volume constraint directly
in the intermediate design while P.II does it through the dilated design,
whose upper bound $V^{*}_{\mathrm{dil}}$ is scaled according to
$V^{*}_{\mathrm{int}}$. From Eq. (6.3) it is clear that
P.I$({{\bm{\bar{\rho}}}_{\mathrm{ero}}},{{\bm{\bar{\rho}}}_{\mathrm{int}}})$
and
P.II$({{\bm{\bar{\rho}}}_{\mathrm{ero}}},{{\bm{\bar{\rho}}}_{\mathrm{int}}},{{\bm{\bar{\rho}}}_{\mathrm{dil}}})$,
or alternatively, P.I$({\eta_{\mathrm{ero}}},{\eta_{\mathrm{int}}})$ and
P.II$({\eta_{\mathrm{ero}}},{\eta_{\mathrm{int}}},{\eta_{\mathrm{dil}}})$.
Therefore, P.I is not influenced by the dilation threshold nor the dilated
design.

We recall that the condition of robustness imposed for the void phase involves
the dilated design, which has to project a cavity with at least one void
element to be present in the 3 fields that constitute the robust formulation.
The influence of the dilated design on the minimum size of the void phase can
be seen graphically in Fig. 6(a) by considering a fixed erosion threshold and
different dilation thresholds. For example, for the set
$[{\eta_{\mathrm{dil}}}$, ${\eta_{\mathrm{ero}}}]$ consider the points
$[0.14\>,\>0.60]$ and $[0.40\;,\;0.60]$. The first point corresponds to a
length scale where the minimum size of the void phase is equal to that of the
solid phase, while the second point imposes the size of the void phase twice
that of the solid phase. In the following, problems P.I and P.II are solved
using the two sets of thresholds indicated previously. To help to see the
differences, the volume fraction ratio has been increased to $0.3$. These sets
of parameters are resumed in Table 4. The results are summarized in Fig. 17.

Clearly, the results obtained using P.II show a length scale consistent with
the expected one, as the design in Fig. 17(a) features bigger reentrant
corners than the design in Fig. 17(b). However, for any value of
${\eta_{\mathrm{dil}}}$, the result from P.I is always the same and
corresponds to that shown in Fig. 17(d). The interpretation that can be made
of the P.I problem is that it uses a dilation threshold equal to the
intermediate one, i.e. ${\eta_{\mathrm{dil}}}={\eta_{\mathrm{int}}}=0.5$. This
can be corroborated by solving P.II with the set
$[{\eta_{\mathrm{dil}}}\>,\>{\eta_{\mathrm{ero}}}]$ equal to $[0.50,0.60]$,
whose result is shown in Fig. 17(c).

$r_{\mathrm{min.Solid}}^{\mathrm{int}}$ | $r_{\mathrm{min.Void}}^{\mathrm{int}}$ | $r_{\mathrm{fil}}$ | $\eta_{e}$ | $\eta_{d}$
---|---|---|---|---
$4$ | $0$ | $12.65$ | $0.60$ | $0.50$
$4$ | $4$ | $12.65$ | $0.60$ | $0.40$
$4$ | $8$ | $12.65$ | $0.60$ | $0.14$
Table 4: Parameters to compare the robust formulation and its simplified form.

(a) P.II, ${\eta_{\mathrm{dil}}}=0.40$.

(b) P.II, ${\eta_{\mathrm{dil}}}=0.14$.

(c) P.II, ${\eta_{\mathrm{dil}}}=0.50$.

(d) P.I for any ${\eta_{\mathrm{dil}}}$.

Figure 17: Heat exchanger design problem using two different variations of the
robust formulation. P.I evaluates the volume constraint directly in the
intermediate design, while P.II does it through the dilated design. Here,
${\eta_{\mathrm{ero}}}=0.60$ and $r_{\mathrm{min.Solid}}^{\mathrm{int}}=4$
elements.

##  7 Conclusion

The robust topology optimization formulation based on uniform manufacturing
errors has gained increasing acceptance in the topology optimization
community. This is mainly due to its ability to control the minimum size of
both the solid and void phases, and its potential to be combined with other
topology optimization approaches. Despite the increasing popularity of the
formulation, no method was yet available to easily obtain the filter and
projection parameters that produce the desired minimum length scales. This
need encouraged us to further develop the analytical method proposed by Qian
and Sigmund (2013). The scope and limitations of this method were assessed
using the numerical method of Wang et al (2011) and a set of 2D design results
from two topology optimization problems, the thermal compliance minimization
problem and the non-linear force inverter.

In addition to providing a fast and effective way to obtain the parameters
that impose the desired minimum length scale, this work shows that to obtain
simultaneous control over the minimum sizes of the solid and void phases, it
is necessary to involve the 3 fields (eroded, intermediate and dilated) in the
robust topology optimization problem. For example, for the compliance
minimization problem subject to a volume restriction, it is known that
intermediate and dilated designs can be excluded from the objective function,
but the volume restriction has to be applied to the dilated design in order to
involve all 3 designs in the formulation.

##  8 Replication of results

This manuscript contains two MATLAB codes as supplementary material. The first
is called SizeSolution.m and provides a list of filter and projection
parameters that impose user defined minimum length scales. The second is
called NumericalSolution.m and builds the graphs in Fig. 6 using the numerical
method proposed by Wang et al (2011).

###### Acknowledgements.

The authors acknowledge the research project FAFIL (Fabrication Additive par
Dépôt de Fil), funded by INTERREG and the European Regional Development Fund
(ERDF).

## Conflict of interest

On behalf of all authors, the corresponding author states that there is no
conflict of interest.

## References

  * Andreasen et al (2020) Andreasen CS, Elingaard MO, Aage N (2020) Level set topology and shape optimization by density methods using cut elements with length scale control. Structural and Multidisciplinary Optimization pp 1–23
  * Bendsøe and Kikuchi (1988) Bendsøe MP, Kikuchi N (1988) Generating optimal topologies in structural design using a homogenization method. Computer methods in applied mechanics and engineering 71(2):197–224
  * Bendsøe (1989) Bendsøe M (1989) Bendsoe, m.p.: Optimal shape design as a material distribution problem. structural optimization 1, 193-202. Structural Optimization 1:193–202
  * Bourdin (2001) Bourdin B (2001) Filters in topology optimization. International journal for numerical methods in engineering 50(9):2143–2158
  * Bruns and Tortorelli (2001) Bruns TE, Tortorelli DA (2001) Topology optimization of non-linear elastic structures and compliant mechanisms. Computer methods in applied mechanics and engineering 190(26-27):3443–3459
  * Chen and Chen (2011) Chen S, Chen W (2011) A new level-set based approach to shape and topology optimization under geometric uncertainty. Structural and Multidisciplinary Optimization 44(1):1–18
  * Christiansen et al (2015) Christiansen R, Lazarov B, Jensen J, Sigmund O (2015) Creating geometrically robust designs for highly sensitive problems using topology optimization: Acoustic cavity design. Structural and Multidisciplinary Optimization 52:737–754
  * Clausen and Andreassen (2017) Clausen A, Andreassen E (2017) On filter boundary conditions in topology optimization. Structural and Multidisciplinary Optimization 56(5):1147–1155
  * da Silva et al (2019) da Silva GA, Beck AT, Sigmund O (2019) Topology optimization of compliant mechanisms with stress constraints and manufacturing error robustness. Computer Methods in Applied Mechanics and Engineering 354:397 – 421
  * Fernández et al (2020) Fernández E, Yang Kk, Koppen S, Alarcón P, Bauduin S, Duysinx P (2020) Imposing minimum and maximum member size, minimum cavity size, and minimum separation distance between solid members in topology optimization. Computer Methods in Applied Mechanics and Engineering 368:113,157
  * Fernández et al (2021) Fernández E, Ayas C, Langelaar M, Duysinx P (2021) Topology optimization for large-scale additive manufacturing: Generating designs tailored to the deposition nozzle size (Under Review.)
  * Kumar and Fernández (2021) Kumar P, Fernández E (2021) A numerical scheme for filter boundary conditions in topology optimization on regular and irregular meshes (Under Review, arXiv:2101.01122v1)
  * Lazarov and Sigmund (2011) Lazarov BS, Sigmund O (2011) Filters in topology optimization based on helmholtz-type differential equations. International Journal for Numerical Methods in Engineering 86(6):765–781
  * Pedersen and Allinger (2006) Pedersen C, Allinger P (2006) Industrial Implementation and Applications of Topology Optimization and Future Needs, vol 137, Springer, pp 229–238
  * Pellens et al (2018) Pellens J, Lombaert G, Lazarov B, Schevenels M (2018) Combined length scale and overhang angle control in minimum compliance topology optimization for additive manufacturing. Structural and Multidisciplinary Optimization
  * Qian and Sigmund (2013) Qian X, Sigmund O (2013) Topological design of electromechanical actuators with robustness toward over-and under-etching. Computer Methods in Applied Mechanics and Engineering 253:237–251
  * Sigmund (1997) Sigmund O (1997) On the design of compliant mechanisms using topology optimization. Journal of Structural Mechanics 25(4):493–524
  * Sigmund (2009) Sigmund O (2009) Manufacturing tolerant topology optimization. Acta Mechanica Sinica 25(2):227–239
  * Sigmund and Maute (2013) Sigmund O, Maute K (2013) Topology optimization approaches. Structural and Multidisciplinary Optimization 48(6):1031–1055
  * Silva et al (2020) Silva G, Beck A, Sigmund O (2020) Topology optimization of compliant mechanisms considering stress constraints, manufacturing uncertainty and geometric nonlinearity. Computer Methods in Applied Mechanics and Engineering 365:112,972
  * Wang et al (2011) Wang F, Lazarov BS, Sigmund O (2011) On projection methods, convergence and robust formulations in topology optimization. Structural and Multidisciplinary Optimization 43(6):767–784
  * Wang et al (2011b) Wang F, Jensen J, Sigmund O (2011b) Robust topology optimization of photonic crystal waveguides with tailored dispersion properties. JOSA B 28:387–397
  * Wang et al (2014) Wang F, Lazarov BS, Sigmund O, Jensen JS (2014) Interpolation scheme for fictitious domain techniques and topology optimization of finite strain elastic problems. Computer Methods in Applied Mechanics and Engineering 276:453 – 472
  * Xu et al (2010) Xu S, Cai Y, Cheng G (2010) Volume preserving nonlinear density filter based on heaviside functions. Structural and Multidisciplinary Optimization 41(4):495–505
  * Yan et al (2018) Yan S, Wang F, Sigmund O (2018) On the non-optimality of tree structures for heat conduction. International Journal of Heat and Mass Transfer 122:660–680
  * Zhou et al (2011) Zhou M, Fleury R, Patten S, Stannard N, Mylett D, Gardner S (2011) Topology optimization-practical aspects for industrial applications. In: 9th World Congress on Structural and Multidisciplinary Optimization
  * Zhu et al (2016) Zhu JH, Zhang WH, Xia L (2016) Topology optimization in aircraft and aerospace structures design. Archives of Computational Methods in Engineering 23:595–622