Mixed projection- and density-based topology optimization with applications to structural assemblies

TL;DR

This paper introduces a hybrid topology optimization method combining explicit geometric control with density-based flexibility, specifically tailored for optimizing structural assemblies and their interface shapes.

Contribution

It presents a novel mixed projection- and density-based approach that unifies shape and topology optimization for complex structural assemblies.

Findings

Effective optimization of assembly interfaces using geometric profiles.

Ability to impose local constraints and material property variations.

Numerical examples demonstrate design flexibility and complex shape manipulations.

Abstract

In this paper we present a mixed projection- and density-based topology optimization approach. The aim is to combine the benefits of both parametrizations: the explicit geometric representation provides specific controls on certain design regions while the implicit density representation provides the ultimate design freedom elsewhere. This approach is particularly suited for structural assemblies, where the optimization of the structural topology is coupled with the optimization of the shape of the interface between the sub-components in a unified formulation. The interface between the assemblies is defined by a segmented profile made of linear geometric entities. The geometric coordinates of the nodes connecting the profile segments are used as shape variables in the problem, together with density variables as in conventional topology optimization. The variable profile is used to locally impose specific geometric constraints or to project particular material properties. Examples of the properties considered herein are a local volume constraint, a local maximum length scale control, a variable Young's modulus for the distributed solid material, and spatially variable minimum and maximum length scale. The resulting optimization approach is general and various geometric entities can be used. The potential for complex design manipulations is demonstrated through several numerical examples.

Tables

Table 1. Table 1 : Continuation schemes for the numerical examples. Each step consists of 50 50 50 optimization iterations

	Reference			Sec. 6.1, 6.2			Sec. 6.3			Sec. 6.4.1			Sec. 6.4.2,6.4.3
Step	$p_E$	${\mu }_{fil}$	${\beta }_{HS}$	$p_E$	${\mu }_{fil}$	${\beta }_{HS}$	$p_E$	${\mu }_{fil}$	${\beta }_{HS}$	$p_E$	${\mu }_{fil}$	${\beta }_{HS}$	$p_E$	${\mu }_{fil}$	${\beta }_{HS}$
1	1.00	1.00	1.00	1.00	1.25	1.00	1.00	1.50	2.00	1.00	1.25	1.50	1.25	1.25	1.00
2	1.00	1.41	1.46	1.50	1.77	1.46	1.50	2.12	2.93	1.50	1.77	2.52	1.50	1.77	1.41
3	1.00	2.00	2.14	2.00	2.50	2.14	2.00	3.00	4.29	2.00	2.50	4.24	2.00	2.50	2.00
4	1.00	2.83	3.14	2.50	3.54	3.14	2.50	4.24	6.28	2.50	3.54	7.14	2.50	3.54	2.83
5	1.00	4.00	4.59	3.00	5.00	4.59	3.00	5.00	9.19	3.00	5.00	12.00	3.00	5.00	4.00
6	1.00	5.00	6.73	3.00	5.00	6.73	3.00	5.00	13.45	3.00	5.00	20.18	3.00	5.00	5.66
7	1.00	5.00	9.85	3.00	5.00	9.85	3.00	5.00	19.70	3.00	5.00	33.94	3.00	5.00	8.00
8	1.00	5.00	14.42	3.00	5.00	14.42	3.00	5.00	28.84	3.00	5.00	57.08	3.00	5.00	11.31
9	1.00	5.00	21.11	3.00	5.00	21.11	3.00	5.00	32.00	3.00	5.00	96.00	3.00	5.00	16.00
10	1.00	5.00	30.91	3.00	5.00	30.91	3.00	5.00	32.00	3.00	5.00	100.00	3.00	5.00	22.63
11	1.00	5.00	32.00	3.00	5.00	32.00	3.00	5.00	32.00	3.00	5.00	100.00	3.00	5.00	32.00
12										3.00	5.00	100.00

Table 2. Table 2 : Settings adopted in the optimization algorithms of the numerical examples

Sec.	6.1	6.2	6.3.1	6.3.2	6.4.1	6.4.2	6.4.3
$r_{min}$	10	10	3	3	2-10	3	3-6
$r_{\varphi }$	10	10	3	3	2	3	3
$r_{max}$	n/a	n/a	5	5	n/a	7-14	20
${\beta }_{fil}$	20	20	10-5	10-5	15	20	20
$g_0^{\ast }$	0.4	0.4	0.35	0.35	0.4	0.4	0.4
$\alpha$	n/a	n/a	0.5	0.5	n/a	0.6	0.6
$\gamma$	n/a	n/a	n/a	n/a	4	1	1

Equations

$${\tilde{\rho }}_i=\frac{{\sum }_{j\in N_i}w(\Delta {\text{x}}_{ij}){\rho }_j}{{\sum }_{j\in N_i}w(\Delta {\text{x}}_{ij})}$$

$$w(\textbf{x}_{j})=\max\bigl{(}r_{min}-\|\textbf{x}_{j}-\textbf{x}_{i}\|,0\bigr{)}$$

$$\tilde{\mathbf{\rho }}=\mathcal{D}{\left({\text{H}}_s\right)}^{-1}\left(\text{H}\mathbf{\rho }\right)$$

$$\begin{array}{rl}& H_{ij}=w(\Delta {\text{x}}_{ij}),H_{s,i}=\sum _j{H}_{ij}\end{array}$$

$$\begin{array}{rl}& {\bar{\rho }}_i^{ero}=\frac{tanh\left({\beta }_{HS}{\eta }_{ero}\right)+tanh\left({\beta }_{HS}({\tilde{\rho }}_i-{\eta }_{ero})\right)}{tanh\left({\beta }_{HS}{\eta }_{ero}\right)+tanh\left({\beta }_{HS}(1-{\eta }_{ero})\right)},\\ & {\bar{\rho }}_i^{int}=\frac{tanh\left({\beta }_{HS}{\eta }_{ero}\right)+tanh\left({\beta }_{HS}({\tilde{\rho }}_i-{\eta }_{ero})\right)}{tanh\left({\beta }_{HS}{\eta }_{ero}\right)+tanh\left({\beta }_{HS}(1-{\eta }_{ero})\right)},\\ & {\bar{\rho }}_i^{dil}=\frac{tanh\left({\beta }_{HS}{\eta }_{dil}\right)+tanh\left({\beta }_{HS}({\tilde{\rho }}_i-{\eta }_{dil})\right)}{tanh\left({\beta }_{HS}{\eta }_{dil}\right)+tanh\left({\beta }_{HS}(1-{\eta }_{dil})\right)},\\ & \text{with e.g. }{\eta }_{ero}=0.6,{\eta }_{int}=0.5,{\eta }_{dil}=0.4\end{array}$$

$${\varphi }_{i,P}=e^{-\frac{1}{2}{\left(\frac{d_{i,P}^2}{{\beta }_{\varphi }^2}\right)}^{{\mu }_{\varphi }}},\forall P$$

$${\varphi }_P=e^{-\frac{1}{2}{\left(\frac{d_{i,P}^2}{{\beta }_{\varphi }^2}\right)}^{{\mu }_{\varphi }}}+e^{-\frac{1}{2}{\left(\frac{d_{j,P}^2}{{\beta }_{\varphi }^2}\right)}^{{\mu }_{\varphi }}},\forall P$$

$$\begin{array}{rl}& {\varphi }_i=e^{-\frac{1}{2}{\left(\frac{{\bar{d}}_i^2}{{\beta }_{\varphi }^2}\right)}^{{\mu }_{\varphi }}}\text{ for }i=1,\dots ,N_{ele}\\ & \text{where }\forall i:\\ & {\bar{d}}_i^2=\min \left([d_{i,1}^2,d_{i,2}^2]\right)=\\ & =d_{max}^2-\max \left([d_{max}^2-d_{i,1}^2,d_{max}^2-d_{i,2}^2]\right)\end{array}$$

$$\begin{array}{rl}& {\varphi }_i=e^{-\frac{1}{2}{\left(\frac{{\bar{d}}_i^2}{{\beta }_{\varphi }^2}\right)}^{{\mu }_{\varphi }}}\text{ for }i=1,\dots ,N_{ele}\\ & \text{where }\forall i:\\ & {\bar{d}}_i^2=d_{max}^2-\frac{{\left(d_{max}^2-d_{i,1}^2\right)}^{q+1}+{\left(d_{max}^2-d_{i,2}^2\right)}^{q+1}}{{\left(d_{max}^2-d_{i,1}^2\right)}^q+{\left(d_{max}^2-d_{i,2}^2\right)}^q}\end{array}$$

$${\tilde{d}}_i^2=\frac{{\sum }_{j\in N_i}w(\Delta {\text{x}}_{ij}){\bar{d}}_j^2}{{\sum }_{j\in N_i}w(x_j)}$$

$${\tilde{\text{d}}}^2=\mathcal{D}{\left({\text{H}}_s^{\varphi }\right)}^{-1}\left({\text{H}}^{\varphi }{\bar{\text{d}}}^2\right)$$

$$\begin{array}{rl}& H_{ij}^{\varphi }=w(\Delta {\text{x}}_{ij}),H_{s,i}^{\varphi }=\sum _j{H}_{ij}.\end{array}$$

$${\varphi }_i=e^{-\frac{1}{2}{\left(\frac{{\tilde{d}}_i^2}{{\beta }_{\varphi }^2}\right)}^{{\mu }_{\varphi }}}\text{ for }i=1,\dots ,N_{ele}$$

$$\begin{array}{rl}\underset{\mathbf{\rho },\text{x}}{\text{minimize :}}& f(\mathbf{\rho },\text{x})\\ \text{subject to :}& g_k(\mathbf{\rho },\text{x})\le 0,k=0,...,m\\ \text{ :}& 0\le {\rho }_i\le 1,i=1,...,N_{ele}\\ \text{ :}& x_{lb}\le x_j\le x_{ub},j=1,...,N_{node}\\ \text{with :}& \text{K}(\mathbf{\rho },\text{x})\text{u}=\text{f}\end{array}$$

$$f(\mathbf{\rho },\text{x})={\text{f}}^T\text{u}({\bar{\mathbf{\rho }}}^{ero},\text{x}).$$

$$\text{K}({\bar{\mathbf{\rho }}}^{ero},\text{x})\text{u}=\text{f}.$$

$$E_i=E_{min}+(E_{max}-E_{min}){\left({\bar{\rho }}_i^{ero}\right)}^{p_E}$$

$$g_0(\mathbf{\rho })=\frac{{\sum }_{i=1}^{N_{ele}}{\bar{\rho }}_i^{dil}{v}_i}{{\sum }_{i=1}^{N_{ele}}{v}_i}-g_{0,dil}^{\ast }\le 0$$

$$g_1(\mathbf{\rho },\text{x})=\frac{{\sum }_{i=1}^{N_{ele}}{\bar{\rho }}_i^{dil}{v}_i{\varphi }_i(\text{x})}{{\sum }_{i=1}^{N_{ele}}{v}_i{\varphi }_i(\text{x})}-g_{1,dil}^{\ast }\le 0$$

$${\hat{\rho }}_i=\frac{{\sum }_{j\in {\hat{N}}_i}{\bar{\rho }}_j^{dil}}{{\sum }_{j\in {\hat{N}}_i}1}$$

$$\hat{\mathbf{\rho }}=\mathcal{D}{\left({\text{H}}_s^{ls}\right)}^{-1}\left({\text{H}}^{ls}{\bar{\mathbf{\rho }}}^{dil}\right)$$

$$\begin{array}{rl}& H_{ij}^{ls}=w_{ij}^{ls},H_{s,i}^{ls}=\sum _j{H}_{ij}^{ls}\end{array}$$

$$w_{ij}^{ls}=\{\begin{array}{ll}1,& \text{if }j\in {\hat{N}}_i\\ 0,& \text{otherwise}\end{array}$$

$${\hat{\rho }}_i{\varphi }_i\le \alpha \forall i=1,\dots ,N_{ele}.$$

$$\underset{i}{\max }({\hat{\rho }}_i{\varphi }_i)\le \alpha \text{for }i=1,\dots ,N_{ele}$$

$${\left(\sum _{i=1}^{N_{ele}}{\hat{\rho }}_i^p{\varphi }_i\right)}^{1/p}\le {\left(\sum _{i=1}^{N_{ele}}{\alpha }^p{\varphi }_i\right)}^{1/p}$$

$$g_2(\hat{\rho },\text{x})={\left(\frac{{\sum }_{i=1}^{N_{ele}}{\hat{\rho }}_i^p{\varphi }_i}{{\sum }_{i=1}^{N_{ele}}{\varphi }_i}\right)}^{1/p}-\alpha \le 0$$

$$\Delta x_{max}=\frac{\Delta y}{tan({\theta }_{max})}$$

$$g_3(\text{x})=\max \left(\Delta x_i^2\right)\le \Delta x_{max}^2,\text{for }i=1,\dots ,N_{node}-1$$

$$g_3(\text{x})=\frac{1}{\Delta x_{max}^2}{\left(\sum _{i=1}^{N_{node}-1}{\left(\Delta x_i^2\right)}^p\right)}^{\frac{1}{p}}-1\le 0.$$

Full text

Keywords— topology optimization, shape optimization, projection methods, structural assembly, robust approach

1 Introduction

Topology optimization has experienced an incredible development since its introduction in the seminal paper by Bendsøe and Kikuchi in 1988 (Bendsøe and Kikuchi, 1988) where homogenization of porous micro-structures was used as the underlying parametrization. In general terms, topology optimization addresses the engineering question on how to distribute material within a given domain in order to obtain the best performance of the system considered. The method allows for a great design freedom and often leads to unexpected material layouts that are far from engineering intuition. Among the different directions along which topology optimization has developed, the density-based (Bendsøe, 1989; Zhou and Rozvany, 1991) and level-set (Allaire et al., 2002, 2004; Wang et al., 2003) parametrizations have so far received the most significant attention. An extensive review of various aspects of density-based and homogenization-based topology optimization can be found in the monograph Bendsøe and Sigmund (2003). For recent thorough comparative reviews on topology optimization approaches, the interested reader is referred also to Sigmund and Maute (2013); Deaton and Grandhi (2014).

Even though topology optimization has proved to be able to identify innovative design solutions with a high performance compared to more traditional solutions, it has been always an open question how to actually manufacture these designs. In order to comply with the various manufacturing technologies, certain geometric limitations need to be imposed on the design outcome. While density-based and implicit level-set parametrizations offer significant design freedom, they also suffer from a certain drawback—the absence of direct geometrical control. Indeed, some geometric properties of the optimized structure can be controlled implicitly. Minimum length scale or thickness is the property that attracted most attention of researchers, leading to the formulation of various techniques (Bruns and Tortorelli, 2001; Bourdin, 2001; Guest et al., 2004; Wang et al., 2011; Allaire et al., 2016, to mention a few). It has been shown that also maximum length scale or thickness can be controlled: in density-based procedures using filters and projections (Guest, 2009; Lazarov et al., 2016; Wu et al., 2018) and in level-set approaches using the signed distance function (Allaire et al., 2016). Another more recent geometric limitation that has been addressed by various researchers is the overhang constraint in additive manufacturing (e.g. Gaynor and Guest, 2016; Langelaar, 2017; Allaire et al., 2017; Qian, 2017). Despite these tremendous achievements, the common density-based and level-set approaches do not offer complete and direct control over the resulting geometry, that can be necessary for adapting the computational design procedure to certain manufacturing scenarios.

As an alternative to the density-based and level-set approaches, another class of topology optimization procedures has emerged recently that is based on explicit parametrization of the design. We refer to this class as “geometric projection” as defined in Norato et al. (2015) and include under this definition a variety of parametrizations that see the topological design as a collection of explicit geometric entities that are projected onto a fixed continuum grid for the purpose of finite element analysis. Among those methods one can find the following examples:

1. Early contributions on spline-based representations projected onto fixed grids (Lee et al., 2004; Edwards et al., 2007);
2. A Heaviside projection approach for optimizing the layout of discrete objects (Guest, 2015);
3. Direct geometric projection of discrete elements, in particular planar structures composed of bars (Norato et al., 2015) and three-dimensional structures made of plates (Zhang et al., 2016); and
4. The Methods of Moving Morphable Components (MMC) and Moving Morphable Voids (MMV) that define the topology by projecting explicit geometric entities such as beams and closed splines onto the continuum domain (Guo et al., 2014; Zhang et al., 2017). From a manufacturing-oriented perspective, these methodologies allow for a more direct control over the geometry compared to traditional topology optimization approaches. In fact, most of the procedures mentioned above are based on explicit geometric information thus they provide a direct link between topology optimization and CAD representation. This opens possibilities for direct geometric control—some examples are overhang limitations (Guo et al., 2017) and hole area and boundary curvatures (Yoely et al., 2018). At the same time, traditional topology optimization approaches offer a truly free-form design parametrization that results in rich design spaces. It is not straightforward to obtain this abundance of design possibilities using explicit design parametrizations (Seo et al., 2010). Hence the following question is raised: Can one combine the richness of traditional topology optimization with the geometric control of projection methods?

Thus, this paper presents a design parametrization that enables to get the best of both worlds by mixing projection- and density-based topology optimization. The underlying principle is based on three pillars:

1. An explicit geometric representation provides control in particular regions or for particular design purposes, using shape variables;
2. An implicit density-based representation is used otherwise, using topological variables;
3. Projection functions are used to tie the two parts together, thus coupling the shape and topological variables. A particular application that benefits from such shape-topology coupling is the optimization of pre-stressed concrete beams (Amir and Shakour, 2018). In the current contribution, our main goal is to generalize the approach of Amir and Shakour (2018) and to demonstrate several (out of many) geometric controls that can be obtained. Nevertheless, we are motivated by a particular class of applications, namely structural assemblies—or in other words, structures that are manufactured in parts and then joined together by e.g. hardware assembly, welding, riveting etc. One specific example arose recently with the development of additive manufacturing (AM): The current sizes of printed parts are limited by the sizes of the machines. This is particularly critical in metals, where most AM facilities can manufacture parts that are not longer than a few hundred millimeters. Consequently, for AM of larger structural components, the design should be divided into parts that are manufactured separately and assembled subsequently.

The fact that the ultimate structural component is an assembly of several separate parts may impose certain limitations regarding the geometry and the material properties at the interfaces between the parts. Hence it is not expected that optimizing the component as a whole and subsequently defining interface regions for manufacturing as separate parts will lead to optimal performance. Furthermore, one cannot impose limitations on interface regions without knowing their location, meaning that significant post-processing may be required. The mixed method proposed herein aims to remove these obstacles: We optimize the location and shape of the interfaces between parts using an explicit geometric representation and employ standard density-based topology optimization elsewhere. The outcomes are structures that are optimized with respect to the performance as an assembly while considering limitations at interfaces, without defining their shape and location a-priori.

The remainder of the article is organized as follows: In Sec. 2 we present the coupled parametrization that is based on density and geometric variables. In Sec. 3 we briefly discuss the general type of topology optimization problems considered. Further details on how the mixed formulation is used to control specific entities and to impose limitations related to interfaces of components are provided in Sec. 4. In Sec. 5 we discuss the sensitivity analysis of the various constraints and in Sec. 6 we apply the proposed approach to several example applications. Final conclusions are drawn in Sec. 7.

2 Mixed parametrization of the design problem

In this section we present the parametrization of the topology optimization problem considered in this work. We simultaneously optimize a fictitious density field of the distributed material assigned to each finite element, and the shape of the projection profile used to impose local controls (e.g.,  geometric constraints) or material properties. Hence, the design variables are both fictitious density values used to interpolate the material properties between solid and void regions of the domain, and the geometric parameters that define the projection profile. In the current study, the latter are simply node positions of line segments, however many other geometric representations can be used.

2.1 Density design variables

We follow the popular density-based approach meaning that the material density is expressed in terms of design variables defined at the finite element level that are collected in the vector $\bm{\rho}$. That is, a density variable $\rho_{i}$ is associated to each finite element. Therefore, we have $N_{ele}$ density design variables, where $N_{ele}$ is the number of finite elements adopted for the numerical approximation of the problem. The density variables $\rho_{i}$ can assume $0-1$ values (i.e. $\rho_{i}=\{0,1\}$), defining the presence of void (i.e. $\rho_{i}=0$) or solid material (i.e. $\rho_{i}=1$) in the $i$-th element. However, since for optimization we will use a gradient-based procedure all the variables need to be continuous. Hence, we relax the density variables’ definition and we allow them to assume also the intermediate values between their upper and lower bounds: $0\leq\rho_{i}\leq 1$. In Sec. 3 we will provide further details on the methodology adopted to converge to nearly discrete optimized designs starting from a continuous problem formulation.

2.2 Shape variables of the projection profiles

The projection of selected properties is performed considering one or several evolving 1-D profiles whose location in space is defined through additional geometric variables. Each 1-D profile is a piece-wise linear geometric entity, composed of a chain of piece-wise linear segments. The influence of each profile extends to a sub-domain identified by a strip along its development. For example, in Fig. 1 we show a portion of a horizontal profile made of two linear segments, and three nodes. In this case, we shall consider the coordinates $x_{i}$ of the nodes fixed, and the coordinates $y_{i}$ as additional variables of the problem. In the eventuality of a profile with vertical direction of development, we shall consider the coordinates $y_{i}$ of the nodes fixed, and the coordinates $x_{i}$ as additional variables of the problems. In this work, the geometric description of the profile used for projection is quite simple. In fact, it is made of linear segments composing a chain, and we allow for two cases: horizontal profile (i.e. the $x_{i}$ coordinates of the profile nodes are fixed), and vertical profile (i.e. the $y_{i}$ coordinates of the profile nodes are fixed). Nevertheless, the approach is general and can accommodate more advanced line descriptions (such as splines) and in principle any explicit geometric entity. Furthermore, it can be used in combination with state-of-the-art projection techniques (Guo et al., 2014; Norato et al., 2015).

2.3 Density-based parametrization

We consider a three-field density representation, i.e $\bm{\rho}\rightarrow\tilde{\bm{\rho}}\rightarrow\bar{\bm{\rho}}$ . First, the density variable field $\bm{\rho}$ is regularized to avoid checkerboard patterns and mesh-dependent solutions though the well-known density filter (Bruns and Tortorelli, 2001; Bourdin, 2001). The resulting filtered density field is $\tilde{\bm{\rho}}$ whose components are:

[TABLE]

where $w(\Delta\textbf{x}_{ij})$ is a linear weight function:

[TABLE]

where $r_{min}$ is the specified filter radius, $\textbf{x}_{i}$ is the position of the centroid of the element $i$, and $\textbf{x}_{j}$ is the position of the centroid of the element $j$. If the element $j$ is in the neighborhood $N_{i}$ of the element $i$ (defined by the radius $r_{min}$) the weight $w(\textbf{x}_{j})$ has a positive value different from zero. In a matrix form Eq. (1) can be stated as follows:

[TABLE]

where H is a $[N_{ele}\times N_{ele}]$ matrix, $\textbf{H}_{s}$ is a $[N_{ele}\times 1]$ vector, and $\mathcal{D}()$ is an operator that transforms a vector into a diagonal matrix and vice versa, similarly to the diag() MATLAB function. The entries of H and $\textbf{H}_{s}$ are defined as follows:

[TABLE]

Additionally, we project the filtered densities $\tilde{\bm{\rho}}$ into $\bar{\bm{\rho}}$ with smooth Heaviside functions to improve the convergence of the algorithm towards crisp solid-void material distributions in the final optimized designs (Guest et al., 2004; Xu et al., 2010). In particular, we introduce an eroded density vector $\bar{\bm{\rho}}^{ero}$ and a dilated density vector $\bar{\bm{\rho}}^{dil}$, in addition to the intermediate density vector $\bar{\bm{\rho}}^{int}$ (Sigmund, 2009; Wang et al., 2011). It should be noted that $\bar{\bm{\rho}}^{ero}$, $\bar{\bm{\rho}}^{int}$, and $\bar{\bm{\rho}}^{dil}$ depend explicitly on $\tilde{\bm{\rho}}$, which depends explicitly on $\bm{\rho}$. In problems involving only stiffness and volume it is possible to consider the eroded density field to quantify the worst case for stiffness, and the dilated field to quantify the worst case for volume (Lazarov et al., 2016). The above mentioned projected density vectors are defined as follows:

[TABLE]

where $\beta_{HS}$ control the sharpness of the projection in the transition zone; and $\eta_{ero}$, $\eta_{int}$, and $\eta_{dil}$ are the projection threshold parameters. It should be noted that the final design solution is represented by $\bar{\bm{\rho}}^{int}$, and in the current context it is used during optimization only for updating the volume fraction for the volume constraint. The field $\bar{\bm{\rho}}^{dil}$ is used to evaluate the volume constraint during the optimization process, and $\bar{\bm{\rho}}^{ero}$ is used to calculate the structural stiffness and to evaluate the compliance functional. More details regarding the optimization problem formulation will be given in Sec. 3.

2.4 Coupling between geometric and density variables

We consider piece-wise linear profiles similar to the one shown in Fig. 2 in red. In this case the profile is made of four segments, and the nodes of the segments have variable coordinates. The profile of Fig. 2 develops along the vertical $y$ direction. The $y$ coordinates of the nodes are fixed, the $x$ coordinates of the nodes are variable. Each point $P$ in space has a distance $d_{i,P}$ from the segment $i$ of the profile, measured perpendicularly from the segment. We define also a maximum distance $\beta_{\phi}$, that delimits the region of influence of the profile colored in grey in Fig. 2. Points whose distance is smaller than $\beta_{\phi}$ will be affected by the projection process while points outside this area will not be affected. We consider also the case of a profile that develops along the horizontal direction $x$. In that case, the $x$ coordinates of the nodes are fixed, the $y$ coordinates of the nodes are variable.

For a given configuration of the $i-th$ profile, and once each point $P$ has a distance parameter defined, the projection is defined through a Super-Gaussian function:

[TABLE]

where $\mu_{\phi}$ controls the sharpness of the Super-Gaussian projection.

In the case of two profiles $i$ and $j$ that are either parallel or perpendicular in terms of their direction of development, each point in space will have two distances with respect to two segments associated to each profile: $d_{i,P}$, $d_{j,P}$. If we perform the projection summing up the two contributions we have the following:

[TABLE]

The projection approach based on the sum of the projections of each profile in Eq. (7) results in the graphical behavior shown in Fig. 3.

As it can be seen, when the two profiles are close to each other the points that lay in the vicinity of the two profiles are associated to projected variables $\phi$ whose values are bigger than one. This is an undesired behavior that we wish to avoid. We want the points close to two segments of two different profiles to have a projected variable $\phi$ at most equal to one. The key for avoiding an overestimation of the projection, is to assign to each $i$-th point only the minimum of the two distances to the two profiles:

[TABLE]

The formulation of Eq. (8) allows to select only the minimum distance based on the maximum distance definition. Thus, in a differentiable form Eq. (8) becomes:

[TABLE]

In Eq. (9) $q$ is a large number, e.g. $10^{6}$.

In Fig. 5 it can now be observed that when the two profiles are close to each other the points that lay in the vicinity of the two profiles are associated to projected variables $\phi$ whose values are at most equal to one.

The projected quantity $\phi(\bar{d})$ of Eq. (9) is characterized by sharp transitions between consecutive segments because of the piece-wise linear definition of the projection profile. This can be seen in Fig. 6(b), for the given profile shown in Fig. 6(a).

To smooth the corners of the projected strip, we filter the distance field $\bar{d}^{2}_{e}$,

[TABLE]

or in a matrix form:

[TABLE]

where $\textbf{H}^{\phi}$ is a $[N_{ele}\times N_{ele}]$ matrix, $\textbf{H}^{\phi}_{s}$ is a $[N_{ele}\times 1]$ vector. The entries of $\textbf{H}^{\phi}$ and $\textbf{H}^{\phi}_{s}$ are defined as follows:

[TABLE]

The benefit of the filtering process applied to $\bar{d}$ can be observed in Fig. 6(c). The figure highlights the smoother transitions of the projected quantity between consecutive profile segments. Consequently, the final Super-Gaussian projection is:

[TABLE]

3 Problem formulation

In this section we present the formulation of the optimization problem based on the design variables and the parametrization described in the previous section. We provide also details regarding the governing equations of the problem and additional computational considerations. In particular, the optimization problem considered herein can be formulated in general terms as follows:

[TABLE]

where $\bm{\rho}$ and x are the vectors of the design variables, corresponding to densities and nodal coordinates, respectively; $f(\bm{\rho},\textbf{x})$ is the objective function; $g_{k}(\bm{\rho},\textbf{x})$ are the $m$ inequality constraints of the problem; $N_{ele}$ is the number of finite elements used to discretize the problem; $N_{node}$ is the number of nodes defining the projection profile; and $x_{ub}$ and $x_{lb}$ are predefined upper and lower bounds for the variable geometric coordinates of the projection profiles’ nodes. It should be noted that the shape variables $x_{i}$ are mapped from a normalized domain $[0-1]$ to the actual domain as follows: $X_{i}=x_{i}\,N_{ele,x}\,a$, where $N_{ele,x}$ is the number of elements in the $x$ direction and $a$ is the size of each element. In case of a vertical profile, the same applies to $y_{i}$, $Y_{i}$, and $N_{ele,y}$. Finally, $\textbf{K}(\bm{\rho},\textbf{x})$ is the stiffness matrix of the structural system considered and u are the displacements of the deformed structure due to the external loads f.

As it has been mentioned in Sec. 2.3, we apply a ‘robust’ topology optimization approach (Wang et al., 2011). This improves the control over the minimum length scale, and for problems that consider only stiffness and volume, the robust approach results in considering the eroded design to compute the stiffness and the dilated design to compute the structural volume (Lazarov et al., 2016). The objective function minimized in the topology optimization problem is the structural compliance evaluated considering the eroded density field, i.e. $\bar{\bm{\rho}}^{ero}$. In particular, the compliance is defined as follows:

[TABLE]

The displacement field is computed by solving the linear system of equations for the structural equilibrium. The structural stiffness matrix is computed based on the current value of the design variables, therefore:

[TABLE]

The definition of Young’s modulus $E$ of each element $i$ in the eroded layout is based on the Modifed SIMP interpolation scheme (Sigmund and Torquato, 1997):

[TABLE]

where $E_{min}$ is associated to an ersatz material that represents void and it has a relatively small value to avoid singularities in the stiffness matrix (e.g. $10^{-6}$); $E_{max}$ is the actual Young’s modulus of the distributed material; $p_{E}$ is a penalization factor that for values bigger than one penalizes the intermediate values of $\bar{\rho}^{ero}_{i}$, making them uneconomical and thus driving the optimizer towards near discrete optimized designs; last, $\bar{\rho}^{ero}_{i}$ is defined according to Eq. (5).

We consider also a volume constraint on the full structural domain:

[TABLE]

where $v_{i}$ is the volume of the $i$-th finite element; and $g_{0,dil}^{*}$ is the available solid volume fraction in the dilated layout.

4 Control over projected entities

The topology optimization approach discussed herein is quite flexible as it can control several design features by combining projection- and density-based design variables in the problem formulation. In this section we discuss several of the features that can be controlled. These are a subset of the many possibilities and they have been selected to show the potential of the methodology presented. The control of selected features is expressed through optimization constraints $g_{k}(\bm{\rho},\textbf{x})$ in the optimization problem formulation (14), or by projecting specific mechanical properties.

The first example of a feature that can be controlled through a localized constraint is the total volume occupied by the structural material in the projection area. The possibility of controlling the amount of material in a specific area can be important in the manufacturing process of structural assemblies. In these cases, it may be desirable to reduce the amount of material that needs to be welded or connected at the interface of different sub-components of a structural assembly. It could also be relevant for structural assemblies where one wishes to reduce the structural complexity at the connecting interface of sub-components. This type of control can be defined through a local volume constraint, affecting only the portion of the structure that is included in the projection area. As the shape of the projection profile is also optimized together with the structural topology, this local constraint is imposed on a sub domain that changes during the optimization process depending on the shape of the projection profile. In particular, we define the following constraint:

[TABLE]

where $\bar{\rho}^{dil}_{i}$ is the projected density of the $i$-th element in the dilated layout; $v_{i}$ is the volume of the $i$-th element; $\phi_{i}(\textbf{x})$ is the projection function evaluated in correspondence to the $i$-th element, expressed explicitly as a function of the projection profile coordinates variables x; and $g_{1,dil}^{*}$ is the solid volume fraction allowed in the projection area in the dilated layout. Essentially, Eq. (19) represents a volume constraint imposed only on a sub-domain identified by those elements for which $\phi=1$ (i.e. the gray area of Fig. 2).

Next, we consider an additional constraint that controls the maximum length scale of the solid phase (Guest, 2009; Wu et al., 2018), and we impose it only in the projection area . To this end, we first introduce the vector $\hat{\bm{\rho}}$ whose components $\hat{\rho}_{i}$ describe the local material distribution in a neighborhood $\hat{N}_{i}$ of a element $i$:

[TABLE]

or in a matrix form:

[TABLE]

where $\textbf{H}^{ls}$ is a $[N_{ele}\times N_{ele}]$ matrix, $\textbf{H}^{ls}_{s}$ is a $[N_{ele}\times 1]$ vector. The entries of $\textbf{H}^{ls}$ and $\textbf{H}^{ls}_{s}$ are defined as follows:

[TABLE]

with

[TABLE]

and we require that:

[TABLE]

Eq. (24) represents $N_{ele}$ constraints on the averaged percentage of solid material in a neighborhood of each element $i$ (i.e. $\hat{\rho}_{i}$), where the maximum allowed percentage of material is $\alpha$. These constraints are equivalent to:

[TABLE]

The $\max$ function in Eq. (25) is not differentiable and therefore not suitable for numerical optimization approaches based on first-order information. To reduce the number of constraints in (24) from $N_{ele}$ to one in a differentiable manner, we aggregate them with a $p$-norm into a single constraint that approximates the maximum value of the set considered similarly to Wu et al. (2018):

[TABLE]

In Eq. (26) the quantities $\phi_{i}$ are not raised to the power of $p$. This is done to account also for the elements on the perimeter of the projection area (where $0<\phi_{i}<1$) in the constraint (26) during the early stages of the optimization. However, as the optimization progresses and the projected quantities $\phi_{i}$ assume near discrete [math]-$1$ values, the fact that the quantities $\phi_{i}$ are not raised to the power of $p$ becomes less relevant. After rearranging the terms we obtain the following constraint formulation:

[TABLE]

The next constraint that we wish to discuss is related to the maximum allowed slope of the segments composing the projection profile. As it has been already mentioned, the geometric description of the projection profiles discussed herein relies on simple piece-wise linear segmented profiles. Large variations of the inclination angles between contiguous segments can result in segmented projection profiles with undesired sharp corners (see Fig. 6). This issue has been mitigated with the introduction of a filtered distance of the elements with respect to their associated profile segment. Here, we discuss an additional precaution that can be taken to further regularize the projection profile: We consider a constraint on the slope of the profile segments. If we consider a generic segment belonging to a vertical profile as shown in Fig. 7, $\Delta x$ is the difference between the $x$ coordinates of the segment nodes, $\Delta y$ is the difference between the $y$ coordinates of the segment nodes, and $\theta$ is the angle of the segment with respect to the horizontal direction. In the case of a vertical profile, $\Delta x$ is variable whereas $\Delta y$ fixed. The opposite would be true for a horizontal profile.

A limit on the maximum slope $\theta_{max}$ is equivalent to imposing a limit on the maximum $\Delta x$ for a fixed $\Delta y$:

[TABLE]

Consequently, the constraint is formulated as follows:

[TABLE]

The constraint (29) is equivalent to $N_{node}$-$1$ linear constraints for the shape variables $x_{i}$. It has been defined as a single aggregated constraint to facilitate its use with the optimization algorithm chosen for this work, the Method of Moving Asymptotes (Svanberg, 1987), which is less suited for liner constraints, and in general large numbers of constraints. It should be mentioned, however, that the constraints of Eq. (29) could have been treated separately with other optimization algorithms more suited for larger number of constraints (e.g. Sequential Linear Programming, Sequential Quadratic Programming (Nocedal and Wright, 2006)). Subsequently, we replace Eq. (29) with an approximated differentiable formulation based on a $p$-norm:

[TABLE]

The constraint $g_{3}(\textbf{x})$ defined in Eq. (30) limits the inclination of the projection profile segments, and leads to more regular optimized projection profiles.

The projection variables $\phi_{i}$ allow also to assign different material properties to the elements included in the projection area. In particular, if we wish to assign a Young’s modulus equal to $E_{\phi}$ to the elements in the projection area, we can do that through a material interpolation scheme based on the Modified SIMP:

[TABLE]

where the parameter $r_{E}$ is introduced to prescribe a different Young’s modulus to the elements inside the projection area, where $\phi_{i}=1$ and $E_{i}=E_{\phi}$. The interpolation scheme defined in (31) allows to assign the Young’s modulus $E_{max}$ to the solid finite elements outside the projection area, and the Young’s modulus $E_{\phi}$ to the solid elements inside the projection area. This could be used for example in the case of a structural assembly composed of separate welded components. In this case, $E_{\phi}$ would define the material property along the welded interface identified by the projection profile.

The last controlled feature that we discuss is the spatial variability of the length scale. We implement a density filter with varying radius, similarly to the one suggested in Amir and Lazarov (2018). This allows to obtain a spatially variable minimum or maximum length scale control. More precisely, we consider a density filter definition whose radius varies according to the current relative position in the domain of the element considered and of the projection profile. In order to obtain a variable minimum length scale control, we modify the well known expression for the density filter as follows:

[TABLE]

where $w(\Delta\textbf{x}_{ij},\phi_{i})$ is the modified weight function:

[TABLE]

In (33), the variable filter radius is defined as:

[TABLE]

where $r_{min}$ is the assigned filter radius, and $\gamma$ the filter radius amplification factor. For example if $\gamma=1$, the filter radius is doubled in the projection area, i.e. where $\phi=1$. In a matrix form the filtering transformation (32) reads as follows:

[TABLE]

where H is a $[N_{ele}\times N_{ele}]$ matrix, $\textbf{H}_{s}$ is a $[N_{ele}\times 1]$ vector. The entries of H and $\textbf{H}_{s}$ are defined as follows:

[TABLE]

Alternatively, it is possible to consider a different filter weight. In fact, we can replace the definition (33) with a Gaussian weight as in Amir and Lazarov (2018):

[TABLE]

In Eq. (37) $exp()$ is the exponential function, and $n$ is a predefined positive number. For a density filter with a similar effect of the filter (32), $n$ can be set to $2$.

It is possible to define also a variable maximum length scale control. In regard to this, previously in Eq. (20) the following filtering procedure has been considered:

[TABLE]

where $w(\Delta\textbf{x}_{ij})$ was equivalent to:

[TABLE]

In Eq. (39), $\mathcal{H}(x)$ is the Heaviside function and it returns $1$ for a non zero argument $x$, otherwise it returns [math]. The filter weights of Eq. (39) can be closely approximated also by the Gaussian weight function. By increasing the value of $n$ in Eq. (37), the weights $w_{G}$ become more homogeneous and approximately equal to $1$. However, in Fig. 8(a) it can also be observed that for increasing values of $n$ the effective filter radius reduces up to approximately half of its initial value. To preserve a meaningful definition of the filter weights with respect to the filter radius $r_{min}$, we propose the following modification of the assigned filter radius $r_{min}$ based on the given value of the exponent $n$:

[TABLE]

such that

[TABLE]

It can be observed in Fig. 8(b) that by replacing $r_{min}$ with $\tilde{r}_{min}$ we preserve the desired length scale control also for increasing values of $n$.

Thus, we can now define a constraint for maximum length scale control with variable maximum feature size allowed:

[TABLE]

where

[TABLE]

and $n$ is a positive number higher than $2$ (e.g. $n=\{4,6\}$).

5 Sensitivity analysis

The optimization approach adopted in this work relies on first-order information. Hence, the gradients of the objective and constraints functions need to be calculated. The dependency of the objective function on the optimization problem variables is implicitly defined through the equilibrium equations. Hence in this section we provide additional details regarding the adjoint sensitivity analysis considered to calculate the gradients of the objective function and some of the constraint functions. In particular, some of the functions involved in the problem formulation depend on both the topological variables $\bm{\rho}$ and on the coordinates of the projection profile x, namely the constraints $g_{1}$, $g_{2}$, $g_{4}$, and the objective function $f$ (if the material interpolation scheme (31) is considered). The dependency on both sets of variables originates from the fact that the projection profile configuration defines the subdomain on which some of the constraints are defined (i.e. $g_{1}$, $g_{2}$, and $g_{4}$), and the subdomain on which the solid elements are given the Young’s modulus $E_{\phi}$, if Eq. (31) is considered. The constraint $g_{3}$ depends only on the shape variables x. The explicit calculations of the gradients are provided in the following.

The gradient of the objective function is calculated with an adjoint method where we first define an augmented function $\hat{f}$:

[TABLE]

It should be noted that when the equilibrium is satisfied (i.e. $\textbf{K}(\bm{\rho}_{ero},\textbf{x})\,\textbf{u}-\textbf{f}$=0) then $f=\hat{f}$. The derivative of $\hat{f}$ with respect to a generic variable $\xi$ is given by:

[TABLE]

where we assumed that the force vector f does not depend on the design variables. This does not imply any loss of generality, see for example design-dependent forces in Amir and Shakour (2018). In order to eliminate the unknown implicit derivatives $\frac{\partial\textbf{u}}{\partial\xi}$ from the gradient calculation, the following equation is solved:

[TABLE]

Once the adjoint vector $\bm{\lambda}$ has been calculated, it is possible to calculate the actual gradient:

[TABLE]

We now specify the derivative $\frac{\partial\hat{f}}{\partial\xi}$ with respect to $\bm{\rho}$:

[TABLE]

where

[TABLE]

In Eq. (49), $\mathop{\vphantom{\sum}\mathchoice{\vbox{\hbox{\huge A}}}{\vbox{\hbox{\Large A}}}{\mathrm{A}}{\mathrm{A}}}\displaylimits$ is the matrix assembly operator, and $\textbf{K}_{e,i}$ is the stiffness matrix of each finite element normalized by its interpolated Young’s modulus $E(\bar{\rho}^{ero}_{i})$.

In cases where the material interpolation scheme of Eq. (31) is considered, the stiffness matrix will depend also on the shape coordinates of the projection profile, and as a consequence also the objective function $f$. Hence:

[TABLE]

where

[TABLE]

In Eq. (54), the matrix $\textbf{H}^{\phi}$ and the vector $\textbf{H}_{s}^{\phi}$ are used to perform the filtering of Eq. (10) that transform the distance field $\bar{d}^{2}$ in $\tilde{d}^{2}$.

The functions that define the constraints $g_{i}$ for $i=0,\dots,m$ are formulated explicitly in terms of the variables of the problem, and for this reason the calculation of their gradient does not require a dedicated sensitivity analysis. In particular, the gradient of $g_{0}$ with respect to $\bm{\rho}$ is calculated as follows:

[TABLE]

The gradient of $g_{1}$ with respect to $\bm{\rho}$ is similar to that of $g_{0}$, and in particular:

[TABLE]

Moreover, $g_{1}$ is also a function of x and the associated gradient is defined as follows:

[TABLE]

The gradient of $g_{2}$ with respect to $\bm{\rho}$ is calculated as follows:

[TABLE]

The constraint $g_{2}$ depends also on the shape variables x and its gradient with respect to them is calculated as follows:

[TABLE]

The constraint $g_{3}$ depends solely on the shape variables x. Hence the gradient $\frac{\partial g_{3}}{\partial x_{i}}$ is calculated directly from Eq. (30).

If a filter with variable radius is considered, similarly for example to Eq. (32), the derivative with respect to $x_{i}$ becomes more articulated. In the numerical examples we will also consider variable minimum and maximum length scale controls. In the following we provide example sensitivities to show the effect of filters with variable radius on the sensitivity. If we assume that a variable minimum length scale control is considered, the derivatives of the objective function with respect to the projection shape coordinates $x_{i}$ becomes:

[TABLE]

In Eq. (57):

[TABLE]

where $w_{G}$ was defined in Eq. (41).

Similarly, the gradient of $g_{4}$ with respect to $\bm{\rho}$ is:

[TABLE]

and the gradient of $g_{4}$ with respect to x is:

[TABLE]

in which the results of Eq. (58) still hold.

6 Numerical examples

In this section we discuss several numerical examples that show the capability of the proposed mixed projection- and density-based topology optimization approach. The optimization algorithm is based on a iterative nested approach, where we use the Method of Moving Asymptotes (Svanberg, 1987). The tight coupling between density and geometric shape variables makes it quite challenging to obtain a smooth convergence towards final optimized designs. Hence, during the optimization process we consider conservative moving limits equal to 0.2 for the density variables (i.e. $\rho$) and to 0.005 for the geometric variables that define the shape of the projection profile (i.e. $x$).

In the following numerical examples we consider the continuation schemes shown in Table 1 with steps of $50$ iterations. In Table 1, $p_{E}$ is the penalization parameter of the SIMP interpolation scheme (Eq. (17)); $\mu_{fil}$ defines the sharpness of the Super-Gaussian function used for the projection (Eq. (13)); and $\beta_{HS}$ defines the sharpness of the projection of the design variables $\tilde{\rho}$ in the dilated, intermediate, and eroded layouts (Eq. (5)). In all the examples, $E_{max}$ is equal to $1$, and $E_{min}$ is equal to $10^{-6}$ to avoid singularities in the stiffness matrix. The Poisson’s ratio considered is $\nu=0.3$. Moreover, $\eta_{ero}=0.6$ and $\eta_{dil}=0.4$. During the optimization runs, every $25$ iterations the allowable volume fractions of the dilated layout are updated such that eventually the volume fraction of the intermediate density meets its allowable value.

In the following numerical applications we consider two structural systems: The MBB beam with height H and length L with a ratio of $H/L=1/3$ shown in Fig. 9, and a short cantilever with a height to length ratio equal to $H/L=2/3$ shown in Fig. 11. We first provide reference results for the two structural systems that are used as a base for a systematic comparison on the different design cases that will be discussed. In both cases we consider a “domain extension approach” discussed in Clausen and Andreassen (2017) to avoid boundary effects. In the figures, the dashed black line represents the outer boundaries of the domain, including the added padding. The blue dashed line represents the blueprint layout, where the actual design is located. To obtain the reference results we set $\eta_{ero}=0.75$ and $\eta_{dil}=0.25$.

The MBB beam is discretized with $100\times 300$ finite elements, considering a density filter radius $r=10$ elements, a constant SIMP penalization $p_{E}=1$, and a volume constraint of $40\%$. The optimization ran for $550$ iterations. The optimized topology is shown in Fig. 10 and its final compliance is $f=196.44$.

The short cantilever beam is discretized with $140\times 210$ finite elements, considering a density filter radius $r=8$ elements, a constant SIMP penalization $p_{E}=1$, and a volume constraint of $35\%$. The optimization ran for $550$ iterations. The optimized topology is shown in Fig. 12 and its final compliance is $f=39.89$. In both cases, the load $P=1$.

Table 2 lists additional settings for the optimization algorithms adopted in the numerical examples. In particular: $r_{min}$ is the density filter radius for the minimum length scale control; $r_{\phi}$ is the filter radius used to filter the distance field of the elements with respect to the projection profile; $r_{max}$ is the radius of the maximum length scale control; $\beta_{fil}$ is the distance from the projection profile in each direction and it defines the area of domain affected by the projection; $g_{0}^{*}$ is the allowable volume fraction for the constraint $g_{0}$, and $g_{0,dil}^{*}$ (Eq. (18)) is equal to $1.05\,g_{0}^{*}$ ; $\alpha$ is the volume fraction allowed in the maximum length scale control in each bubble with radius $r_{max}$; $\gamma$ is the amplification factor for $r_{min}$ or $r_{max}$ used to obtain variable length scale controls.

All the following examples are obtained solving Problem (14) with at least two constraints. The first is the total volume constraint $g_{0}$. The second is $g_{3}$ defined in Eq. (30), and it is used in all cases to obtain more regular shapes of the optimized projection profile, by controlling the slopes of the segments composing the projection profile. In particular, we consider an exponent for the $p$-norm in Eq. (30) equal to $10$, and a maximum allowed inclination angle of $60$ degrees. All the numerical analyses were performed on a Linux machine with $8$ Gb of RAM and a dual-core Intel $i7$ CPU at $2.00$ GHz.

6.1 Example 1: MBB beam with local volume constraint

In this example we optimize the MBB beam solving problem (14) with two volume constraints. The first is the total volume constraint $g_{0}$. The second is a local volume constraint applied on the projection area $g_{1}$, which was discussed in Sec. 4. The allowed volume fraction of the intermediate density layout associated to $g_{0}$ is $40\%$, and the volume fraction associated to $g_{1}$ is $25\%$. As it was already mentioned before, we consider also a third constraint to obtain more regular shapes of the optimized projection profile, namely $g_{3}$ which has been defined in Eq. (30). The projection profile is composed of $5$ segments, each spanning $20$ elements in the vertical directions. The $6$ nodes at the segments ends have fixed $y$ coordinates, and variable $x$ coordinates. The geometric variables $x_{i}$ with $i=1,...,6$ are initially set to $0.5$ and their bounds are $x_{lb}=0.33$ and $x_{ub}=0.67$. These bounds correspond to stretches of $50$ elements in the negative and positive horizontal directions. They are represented by dashed lines colored in magenta in Fig. 13. The motivation for defining this allowable region for the interface stems from manufacturing limitations: we consider a case where the whole beam cannot be manufactured in one part, and must be split into two parts whose maximum length is 2/3 of the total length. Then, one would seek the optimal interface location that minimizes compliance while considering properties or limitations on the interface region.

Regarding the parameters that characterize the projection, in this example we consider a radius of the filter (10) $r_{\phi}=10$ elements, and a distance from the projection profile $\beta_{fil}=20$ elements on each side used to define the area on which the local volume constraint $g_{1}$ is applied. Moreover, we initialize the parameter $\mu_{fil}$ of Eq. (13) to $1.25$, and increase it by steps of $\sqrt{2}$ at each continuation step up to a maximum value of $\mu_{fil}^{max}=5$ as shown in Table 1. The parameters that define the robust approach described in Sec. 2.3 are set as follows: $\eta_{ero}=0.6$ and $\eta_{dil}=0.4$. With this numerical example we want to simulate the search for an optimized topology with a volume constraint modified locally in the projection area. This requirement may translate the need for reducing the structural elements’ size at the interface of different parts manufactured separately and subsequently assembled.

The final optimized design has a compliance $f=202.13$ and a volume of the intermediate layout of $V^{int}=40.21\%$. The optimized topology is shown in Fig. 13, where each red circle represents a node of the profile whose horizontal coordinate $x$ is one of the geometric variables. Fig. 14 shows the optimized shape of the projection profile and its associated projection area defined by the variables $\phi_{i}$. As can be observed in Fig. 13, by imposing a reduced amount of material in a sub-portion of the design domain we obtain an optimized design with thinner members in that portion of the domain, and hence a more flexible structure with poorer performance compared to the reference results. More precisely, compared to the reference design of Fig. 10, the obtained design shows a loss in performance of approximately $2.9\%$, even though both structures have about $40\%$ final solid volume fraction. At the same time, the amount of material in the interface region of the result of Fig. 13 is smaller compared to the reference result. This implies that the effort invested in assembly (for example, total welding energy) is expected to reduce compared to the reference result.

6.2 Example 2: MBB beam with a localized Modified SIMP

In this example we optimize the MBB beam solving problem (14) with the local Modified SIMP defined in Eq. (31), considering $r_{E}=0.5$. This implies that the material distributed in the projection area has a Young’s modulus that is half of $E_{max}$. The two constraints $g_{0}$ and $g_{3}$ are also considered in the optimization analysis for this example. Moreover, the projection profile is made of $5$ segments spanning $20$ elements in the vertical directions. The $6$ nodes at the segments ends have fixed $y$ coordinates, and variable $x$ coordinates. The geometric variables $x_{i}$ with $i=1,...,6$ are initially set to $0.5$ and their bounds are $x_{lb}=0.33$ and $x_{ub}=0.67$. These bounds correspond to stretches of $50$ elements in the negative and positive horizontal directions. All the remaining parameters are set according to Table 1 and Table 2. This example shows the possibility of optimizing a structure with an interface between substructures also simultaneously optimized, taking into account the presence of welded material with different mechanical properties connecting the different substructures. It should be noted that the interface between the substructures is defined by the variable shape of the projection profile.

Fig. 15 shows the optimized topology and the optimized configuration of the projection profile. The optimized structure has a compliance $f=207.13$ which is $5.4\%$ higher than the reference design of Fig. 10. The final solid volume of the intermediate layout is $V=40.09\%$. Fig. 16 shows the plot of the Young’s modulus of each solid element of the optimized structure. The black elements have a Young modulus equal to $E_{max}$. The gray elements are inside the area of influence of the projection profiles, and have a reduced Young’s modulus equal to $E_{\phi}=0.5E_{max}$. It is interesting to observe that the optimizer attempts to place the projection profile in the optimized configuration perpendicularly to the intersected structural elements. We presume that by doing this, the optimizer tries to reduce the extent of poor material assigned to the structural elements. Moreover, in Fig 15 it is possible to also observe that the shape of the projection profile is very close to the upper bound of the variables $x_{i}$, namely $x_{ub}$, and almost straight. To further investigate the solution identified by the optimizer we run two additional optimization analyses where we keep the shape variables $x_{i}$ fixed. In one case we set $x_{i}=x_{lb}$ for $i=1,\dots,6$ and we obtain the final optimized topology shown in Fig. 17(a). In the second case we set $x_{i}=0.5$ for $i=1,\dots,6$ which leads to the final optimized topology shown in Fig. 17(b).

The performances of the optimized designs shown in Fig. 17 are both inferior to the performance of our optimized design shown in Fig. 15. This can explain why the optimizer placed the projection profile close to the upper bound $x_{lb}$ as shown in Fig. 15, where presumably the negative effect of the weak material is minimized.

6.3 Example 3: Short cantilever with a localized maximum length scale control

We consider now the optimization of a short cantilever with two projection profiles. In particular, we consider two cases: one case in which the two profiles are vertical, and another case in which one profile is horizontal and another one is vertical. Besides the usual total volume constraint $g_{0}$ and the constraint on the slope of the profiles’ segments $g_{3}$, in this case we consider also the localized volume constraint for local maximum length scale control $g_{2}$ which has been introduced in Sec. 4. For the numerical example of this section, the moving limit for the density variables is set to $0.1$. The moving limit for the geometric variables is set here to $0.005$. This parameter modification was required because during initial numerical experiments the algorithm showed some difficulties in converging towards near discrete final designs. In fact, in this example the portion of domain on which additional controls are imposed is more extended, and this increases the difficulty in converging towards optimized and feasible design solutions. The allowable volume fraction of the intermediate layout is initially set to $35\%$. The allowable volume fraction for the maximum length scale control (the parameter $\alpha$ in Eq. (42)) is set to $50\%$. In both the following examples, the exponent $q$ of the max approximation of Eq. (9) is set to $10^{6}$. The distance from the projection profile is set initially to $\beta_{fil}=10$, and decreased to $5$ at the first continuation scheme step.

6.3.1 Short cantilever with two vertical profiles

The first case considered is the short cantilever with two vertical profiles. The settings for the continuation scheme and the optimization algorithm are shown in Table 1 and Table 2. The two vertical profiles are made of $7$ segments spanning $20$ elements. The $8$ nodes at the segments ends have fixed $y$ coordinates, and variable $x$ coordinates. The geometric variables of the two profiles $x^{1}_{i}$ and $x^{2}_{i}$ have been initially set to $0.25$ and $0.75$. Their bounds are defined by gaps of $30$ elements with respect to the original positions such that $x^{1}_{lb}=0.11$ and $x^{1}_{ub}=0.40$, and similarly $x^{2}_{lb}=0.61$ and $x^{2}_{ub}=0.90$.

Fig. 18 shows the final optimized topology obtained after $550$ optimization iterations, and the final shape and configuration of the two projection profiles. Fig. 18 also displays the upper and lower bounds of the shape variables that are represented by dashed lines colored in magenta. The optimized structure has a compliance $f=44.84$ ($12.4\%$ higher than the reference design of Fig. 12), and a final solid volume of the intermediate layout $V=35.01\%$. The shape nodal variables are represented also in this case by red circles. In particular, in Fig. 18(a) it is possible to observe smaller structural features in the vicinity of the profiles as a direct consequence of the imposed local maximum length scale control $g_{3}$ defined in Eq. (27). This specific feature control may represent for example the requirement of having smaller elements that need to be connected in correspondence of an interface of different parts assembled. From a visual perspective, this topological layout resembles recent results achieved with graded porosity (Schmidt et al., 2019). However, the underlying formulations are quite different. Primarily, we designate the region for imposing maximum length scale using an explicit geometry, whereas in Schmidt et al. (2019) the authors propose manual control, auxiliary density field control or physics-based control.

Fig. 19 shows four intermediate design stages obtained in the initial $100$ optimization iterations. It is possible to observe that initially the algorithm quickly identifies the topology outside of the projection areas. Within the projection areas, the convergence towards a final near discrete design has been observed to be more challenging and slow due to the interplay between shape and density variables.

6.3.2 Short cantilever with orthogonal projection profiles

The short cantilever with orthogonal projection profiles has been optimized using the same settings for the parameters as in the previous example with two vertical profiles (i.e. Sec. 6.3.1). The vertical profile is made of $7$ segments spanning $20$ elements. The horizontal profile is made of $7$ segments spanning $30$ elements. The $8$ nodes at the vertical segments ends have fixed $y$ coordinates, and variable $x$ coordinates. The $8$ nodes at the horizontal segments ends have fixed $x$ coordinates, and variable $y$ coordinates. The variable geometric coordinates of the nodes of the two profiles $y^{1}_{i}$ and $x^{2}_{i}$ have been both initially set to $0.5$. Their bounds are defined by gaps of $30$ elements with respect to the original positions such that $y^{1}_{lb}=0.29$ and $y^{1}_{ub}=0.71$, and similarly $x^{2}_{lb}=0.36$ and $x^{2}_{ub}=0.64$. Fig. 20 shows the final topology obtained after $550$ optimization iterations, and the corresponding shape of the optimized projection area. The design solution obtained in this case is far from being symmetric. The reason could be that initially the horizontal profile is passing exactly through the point of application of the concentrated load. As a consequence the maximum length scale control is imposed initially in the point of application of the load, which significantly weakens the structure. And in fact, in the numerical experiments it has been observed that the optimizer moved the position of the profile away from the point of application of the load from the very beginning of the optimization analysis, thus converging towards a non symmetric final design.

In the two examples discussed in this section, it is possible to observe that the local volume constraint effectively leads towards final optimized designs with a maximum element length scale imposed on the projection areas. In fact, in these areas the optimized designs are characterized by thinner beam-like elements. This demonstrates the capability of the proposed procedure to control the maximum thickness of members in the interface of assemblies, as could be required by certain manufacturing considerations.

6.4 Example 4: Variable minimum/maximum length scale control

In this section we present three examples with different types of controls that impose variable minimum or maximum length scales. More precisely, in the first example we perform the topology optimization of a short cantilever with a spatially variable minimum length scale. We adopt the density filter of Eq. (32) with weights defined according to Eq. (37). Through this filter formulation, a minimum length scale defined by the filter radius $\bar{r}_{min}$ is imposed in the domain inside the projection area (where $\phi=1$), and a minimum length scale defined by the radius $r_{min}$ otherwise. The ratio between $\bar{r}_{min}$ and $r_{min}$ is $5$. The second example considers a spatially variable maximum length scale. For this purpose, we consider the constraint defined in Eq. (42). This constraint is imposed on all the design domain. However, the radius used for imposing the maximum length scale varies, and it is doubled inside the projection area. In this way we allow for thicker features inside the projection area. The third example shows a different approach for achieving a variable maximum length scale control. In particular, we consider the maximum length scale constraint of Eq. (42) with fixed radius $\hat{r}_{min}$ (i.e. $\gamma=0$). At the same time, we consider also a variable density filter for the minimum length scale given in Eq. (32) with filter weights defined in Eq. (37) for $n=2$ and a filter radius doubled in the projection area. The density filter radius for the minimum length scale is kept smaller than that of the maximum length scale in order to be able to impose the desired feature size within the length scale allowed by the maximum length scale control.

6.4.1 Variable minimum length scale

Here we consider a long cantilever beam as shown in Fig. 11 with a $H/L$ ratio of $1/3$. We considered a longer geometry in the horizontal direction of the design domain to allow for a more distinct transition between the different minimum length scales imposed by the different filter radii considered during the optimization. The structure is discretized with $100\times 300$ finite elements with a density filter radius $r_{min}=2$ and a filter radius for the projection profile $r_{\phi}=2$. For the minimum length scale filter, we adopted the filter weights defined in Eq. (41) with a scaling factor of the filter radius in the projection area $\gamma=4$. As a consequence, for the $i$-th element included in the projection profile (i.e. $\phi_{i}=1$) the radius of the minimum length scale is $\bar{r}_{min}=5\,r_{min}$. We consider an allowable volume fraction for the total volume of the intermediate layout equal to $40\%$. The parameters that define the continuation scheme are set as in the previous examples, with the exception of $\beta_{fil}=15$ elements, and $\beta_{HS}$ initialized to $1.5$ and increased by steps of $2^{0.75}$ up to $\beta_{HS,max}=100$. In this example the final value $\beta_{HS,max}$ is significantly higher than that of the other examples because of the large minimum length scale radius $\bar{r}_{min}$ considered in the projection area. A larger filter radius, in fact, requires a more sharp projection in Eq. (5) in order to transform the filtered density filed into a near discrete topology. The optimization process ran for $600$ iterations with steps of the continuation scheme of $50$ iterations. The projection profile is divided into four segments spanning $25$ elements each. Thus a total of five nodes with variable horizontal coordinate $x$ are also considered as deign variables. The nodes’ $x$ coordinates are initially set to $0.5$, and bounded with gaps of $20$ elements. That is, $x_{ub}=0.56$ and $x_{lb}=0.43$. The final optimized topology and projection profile are shown in Fig. 21. Fig. 22 shows the final projection area in black, for which $\phi=1$.

In Fig. 21 it is possible to observe a clear distinction of length scales between the areas of the domain inside and outside the projection area. Outside we observe thinner structural components in accordance with the regular filter radius, and inside only thick members and thick joints exist in accordance with the enlarged filter radius.

6.4.2 Variable maximum length scale

Next, we consider a short cantilever beam as shown in Fig. 11 with a $H/L$ ratio of $1/2$. The structure is discretized with $120\times 240$ finite elements with a density filter radius $r_{min}=3$ elements and a filter radius for the projection profile set to $r_{\phi}=3$ elements. The radius $r_{max}$ for the maximum length scale is set to $7$ elements outside the projection area. We consider an amplification factor $\gamma=1$, such that the $r_{max}$ is doubled inside the projection area. The allowable volume fraction of the intermediate layout is initially set to $40\%$. The allowable volume fraction for the maximum length scale control (the parameter $\alpha$ in Eq. (42)) is set to $60\%$. The projection area is defined in terms of distance from the profile and is set to $\beta_{fil}=20$ elements. The exponent $p$ of the SIMP material interpolation is initialized to $1.25$ and increased to 1.5 at the first continuation scheme step. After, it is increased up to $3$ with steps of $0.5$. The remaining parameter settings are listed in Tables 1 and 2. We consider a single vertical projection profile divided into six segments spanning $20$ elements each. Thus a total of seven nodes with variable horizontal coordinate $x$ are also considered as deign variables. The nodes $x$ cooridnates are initially set to $0.5$, and bounded with gaps of $30$ elements. Hence, their upper bound is $x_{ub}=0.625$ and the lower bound is $x_{lb}=0.375$. The optimization is run for $550$ iterations with steps of the continuation scheme of $50$ iterations.

The compliance of the final optimized structure is $f=87.56$, and it is associated to a blueprint structural layout that occupies $38.2\%$ of the design domain. The final optimized topology and projection profile are shown in Fig. 23. Fig. 24 shows the final projection area in black, for which $\phi=1$. In the final topology it is possible to observe a clear distinction between length scale outside and inside the projection area. Outside, the structure is characterized by thinner and more branched elements, resembling a porous medium. Inside the projection area instead, the number of branching elements is reduced and the elements are thicker. It should be mentioned that the problem formulation discussed in this section would probably benefit from even higher resolutions of the mesh. In fact the maximum length scale control imposed requires high resolutions in order to clearly identify the small features. This of course would imply a higher computational cost and more time to perform the optimization. From an engineering design standpoint, this example shows how one can impose a maximum length scale throughout the design domain, while enlarging it in the assembly interface region in order to reduce the number of member connections.

6.4.3 Maximum length scale control through a variable minimum length scale

Also in this example we consider the short cantilever beam shown in Fig. 11 which has a $H/L$ ratio of $1/2$. The structure is discretized with $120\times 240$ finite elements with a filter radius for the projection profile set to $r_{\phi}=3$ elements. The radius $r_{max}$ for the maximum length scale is set to $20$ element. The density filter radius is set to $r_{min}=3$ elements outside the projection area. We consider an amplification factor $\gamma=1$, such that $r_{min}$ is doubled inside the projection area. The allowable volume fraction of the intermediate layout is initially set to $40\%$. The allowable volume fraction for the maximum length scale control (the parameter $\alpha$ in Eq. (42)) is set to $60\%$. The remaining parameter settings are listed in Tables 1 and 2. We consider a single vertical projection profile divided into six segments spanning $20$ elements each. Thus a total of seven nodes with variable horizontal coordinate $x$ are also considered as design variables. The nodes $x$ coordinates are initially set to $0.5$, and bounded with gaps of $30$ elements. That is, the upper bound is $x_{ub}=0.625$ and the lower bound is $x_{lb}=0.375$ and they are shown in Fig. 25 with the magenta color. The optimization is run for $550$ iterations with steps of the continuation scheme of $50$ iterations.

The compliance of the final optimized structure is $f=64.93$, and it is associated to a blueprint structural layout that occupies $40.12\%$ of the design domain. The final optimized topology and projection profile are shown in Fig. 25. Also in this case, similarly to the results of Sec. 6.4.1, it is possible to observe an increase of the imposed length scale moving from the domain outside to the domain inside the projection area. As expected, the results fulfill the expectation of elements with bigger size inside or in the vicinity of the projection area. Fig. 26 shows the final projection area in black, for which $\phi=1$. In a practical engineering prospective, this example, together with the examples of Sec. 6.4.2 and Sec. 6.4.2, shows an additional way of imposing a maximum length scale throughout the design domain, which is enlarged in the assembly interface region in order to reduce the number of connections between structural elements.

7 Conclusion

We presented a topology optimization approach that blends together projection- and density-based formulations. The resulting formulation relies on an explicit geometric representation, parametrized by shape variables, that allows to impose selective controls over the design in specific regions, while using an implicit density representation elsewhere in the design domain. The two representations are tied together via projection functions that couple the respective shape and density variables.

The mixed formulation is rather general, for two main reasons: 1) Any type of explicit geometric entity can be utilized to define the domain where specific control is desired; and 2) Various geometric and response constraints can be imposed on the regions defined by the projected geometry. In the current study, we are motivated by the optimal design of structural assemblies in 2-D, hence the geometric entities are piece-wise linear segments that define interfaces (or “cuts”) between parts that are manufactured separately and subsequently assembled or joined together. Furthermore, the constraints represent several design considerations that may rise when optimizing an assembled part. Examples are: degraded material properties along the interface; total material volume along the interface; maximum member thickness in the interface region, while no limitation is imposed elsewhere; spatial variation of length scale, e.g. imposing different minimum or maximum length scales near the interface and elsewhere. The examples in Sec. 6 clearly demonstrate the capability to impose such restrictions in specific regions that are not predefined. In other words, the overall topology and the shape of the interface between parts are found simultaneously, giving the optimization the freedom to choose where to impose restrictions while maintaining the full freedom of a density-based parametrization elsewhere.

Acknowledgements

This work has been carried out as part of AATiD –- Advanced Additive Titanium Development Consortium. The authors wish to thank the Israeli Innovation Authority and the industrial partners for their generous financial support.