∎

¹¹institutetext: J. Wang and R. Westermann ²²institutetext: Computer Graphics and Visualzation, Technical University of Munich, Munich, Germany
²²email: {junpeng.wang, westermann}@tum.de ³³institutetext: J. Wu ⁴⁴institutetext: Department of Sustainable Design Engineering, Delft University of Technology, Delft, The Netherlands
Corresponding Author
⁴⁴email: [email protected]

Stress Topology Analysis for Porous Infill Optimization

Junpeng Wang Jun Wu Rüdiger Westermann

(Received: date / Accepted: date)

Abstract

The optimization of porous infill structures via local volume constraints has become a popular approach in topology optimization. In some design settings, however, the iterative optimization process converges only slowly, or not at all even after several hundreds or thousands of iterations. This leads to regions in which a distinct binary design is difficult to achieve. Interpreting intermediate density values by applying a threshold results in large solid or void regions, leading to sub-optimal structures. We find that this convergence issue relates to the topology of the stress tensor field that is simulated when applying the same external forces on the solid design domain. In particular, low convergence is observed in regions around so-called trisector degenerate points. Based on this observation, we propose an automatic initialization process that prescribes the topological skeleton of the stress field into the material field as solid simulation elements. These elements guide the material deposition around the degenerate points, but can also be remodelled or removed during the optimization. We demonstrate significantly improved convergence rates in a number of use cases with complex stress topologies. The improved convergence is demonstrated for infill optimization under homogeneous as well as spatially varying local volume constraints.

Keywords:

Topology optimization Porous infill Stress tensor

^†^†journal: Structural and Multidisciplinary Optimization

1 Introduction

Topology optimization aims at finding the optimal structural layout under relevant design specifications. Topology optimization of multi-scale structures, which dates back to the seminal paper by Bendsøe and Kikuchi (1988), has been a topic of great interest in recent years. The rapid development in this field is partially stimulated by the possibility to fabricate complex structures using additive manufacturing. For an overview of topology optimization approaches for designing multi-scale structures, we refer readers to a recent review article by Wu et al. (2021a).

It has been shown that density-based topology optimization for compliance minimization, under local volume constraints, creates porous infill structures similar to those found in bone (Wu et al., 2018). These bone-mimicking porous structures are lightweight, robust regarding material damages and loading variations, and stable with respect to buckling. The local volume constraints function similarly to maximum length scale control (Guest, 2009). They prevent the forming of large solid regions and, consequently, create porous structures distributed more evenly over the design domain. This approach has been extended, in conjunction with a coating approach proposed by Clausen et al. (2015), to design concurrently structures and porous sub-structures therein, referred to as shell-infill composites (Wu et al., 2017). It has also been applied to design porous shell structures (Träff et al., 2021). Other notable extensions include the design of porous structures with gradation in the porosity and pore size (Schmidt et al., 2019; Das and Sutradhar, 2020), use of multiple materials (Li et al., 2020; Zhao and Zhang, 2021), stress constraints (Kranz et al., 2021), as well as design of self-supporting infill (Liu et al., 2021) and fiber-reinforced structures (Li et al., 2021). Dou proposed a projection operator to implicitly incorporate local volume constraints (Dou, 2020). Besides by density-based approaches, porous infill structures have been designed using an evolutionary design approach (Qiu et al., 2020) and machine learning (Cang et al., 2019).

In this paper, we investigate the convergence behavior of density-based topology optimization with local volume constraints under a single load case. In density-based topology optimization, an important convergence criterion is that the optimized density field converges to a binary or so-called black-white design, i.e., the pseudo density is close to 1 or 0. A few hundred iterations or even more are not uncommon to achieve black-white designs. To improve the convergence rate, a typical solution is to apply a continuation scheme where parameters are updated after a certain number of iterations. However, in some optimization scenarios under local volume constraints, we have observed that certain regions fail to converge to a binary design even after thousands of iterations (see Fig. 1). Interpreting these intermediate density values by applying a threshold results in large solid or void regions, leading to sub-optimal structures.

To analyze the regions where low convergence is observed, we investigate the stress distribution in these regions via trajectory-based visualization (Wang et al., 2020). In particular, we shed light on the relationship between the convergence behavior and the principal stress directions that occurs when simulating on the solid design domain. This approach is inspired by previous work on infill optimization, where uniformly seeded tensor glyphs have confirmed good agreement between the optimized porous infill and the principal stress directions in the solid under load (Wu et al., 2018). In this paper, we exploit advanced mechanisms to perform a topology-based analysis of the stress field, including the use of degenerate points and topological skeletons. At a degenerate point, the principal stress directions cannot be decided, yet a set of hyperbolic and parabolic sectors exist in its surrounding, in which similar patterns of neighboring trajectories are observed (Delmarcelle and Hesselink, 1994). The topological skeleton consists of the boundaries between adjacent sectors—so-called separatrices—and indicate pathways along which the forces are steered towards the degenerate points. In topology optimization, degenerate points have been used to indicate locations where integrability conditions are violated and consistent domain parameterizations cannot be computed (Stutz et al., 2020).

When applying topology analysis to the stress tensor field, it reveals that low convergence occurs around a special type of degenerate points, known as trisectors. Notably, such degenerate points do not always appear, but if so, low convergence is often observed in their surrounding. Due to the isotropy of the stress tensor close to a trisector, the principal stress directions and, thus, a locally consistent binary material layout cannot be decided by the optimizer. In our work, we propose an automatic pre-process that supports the optimizer in finding such a layout, resulting in significantly improved convergence rates in settings where trisectors are paramount. In particular, we build upon the efficient computation of degenerate point locations and separatrices, and prescribe an initial density field where elements along the separatrices are solid and all other elements take an intermediate value (i.e., the local volume upper bound). The solid elements along the separatrices can be changed over the course of the optimization to improve the structural performance, yet we observe that the optimization keeps them more or less unchanged and changes element densities in other parts accordingly. In the vicinity of trisectors, the initialization guides the optimizer towards a stable binary design and enables the optimization process to quickly converge towards a sound global layout. Interestingly, even though one might expect that the imposed initialization biases the optimizer towards a less stiff local optimum, the resulting binary designs show the same or even improved compliance compared to the designs generated by the original approach which exhibit unresolved intermediate density values in the presence of trisectors.

The remainder of this paper is organized as follows. In Section 2, we first review the problem formulation underlying porous infill optimization. Then, in Section 3 we analyse the convergence of porous infill optimization, elaborate on the relationships between optimization convergence and the existence of degenerate points in the stress field, and propose topology-based material initialization to counteract low optimization convergence. The implementation details of performing the topology analysis to the stress tensor field are discussed in Section 4. We demonstrate the effectiveness of our approach in a variety of experiments in Section 5. Section 6 concludes the paper with a discussion of the proposed approach as well as future research directions.

2 Porous Infill Optimization

The low convergence in some design tasks is observed while using the infill optimization approach (Wu et al., 2018), on which and some of its extensions the effectiveness of our method will be demonstrated. For the sake of completeness, we briefly review the formulation of the density-based infill optimization with local volume constraints.

2.1 Local volume constraints

In a discretized design domain, the local volume ( $\bar{\rho}_{e}$ ) of a circular region centered at the centroid of an element, $x_{e}$ , is computed by

\bar{\rho}_{e}=\frac{\sum_{i\in N_{e}}\rho_{i}}{\sum_{i\in N_{e}}1},\quad N_{e}=\{i|\parallel x_{i}-x_{e}\parallel_{2}\leq R_{e}\},\>\>\forall e,

(1)

where $\rho_{i}\in[0,1]$ is the pseudo density for the $i$ -th element. $R_{e}$ denotes the radius of the region on which the local volume is measured.

An upper bound ( $\alpha_{e}$ , $0<\alpha_{e}<1$ ) is imposed on the local volume of each element in the design domain, i.e.,

\bar{\rho}_{e}\leq\alpha_{e}.

(2)

Thus, the local volume constraint involves two parameters, $R_{e}$ and $\alpha_{e}$ . $R_{e}$ indirectly controls the spacing between sub-structures, and $\alpha_{e}$ effectively controls the porosity (Wu et al., 2018). In the original approach, both input fields are prescribed to be homogeneous. Recent developments have demonstrated the use of heterogeneous fields to generate gradations of the porosity and pore size of the optimized porous structures.

Assigning a local volume constraint to each element results in a large number of constraints that need to be considered by the optimizer. Dividing both sides of Eq. 2 by $\alpha_{e}$ , these constraints are aggregated by the $p-$ mean function,

\left(\frac{1}{n}\sum\limits_{e}\left(\frac{\bar{\rho}_{e}}{\alpha_{e}}\right)^{p}\right)^{\frac{1}{p}}\leq 1,

(3)

where $n$ is the number of elements. $p=16$ is found to give a good approximation, and is used in this paper.

2.2 Optimization problem

With the local volume constraint defined, the optimization problem is given by

$\displaystyle\displaystyle\min\limits_{\bm{\phi}}\quad$	$\displaystyle c=\frac{1}{2}\bm{U}^{T}\bm{K}\bm{U},$	(4)
$\displaystyle\displaystyle\mathrm{s.t.}\quad$	$\displaystyle\bm{K}\bm{U}=\bm{F},$	(5)
	$\displaystyle g\left(\bm{\phi}\right)=\left(\frac{1}{n}\sum\limits_{e}\left(\frac{\bar{\rho}_{e}}{\alpha_{e}}\right)^{p}\right)^{\frac{1}{p}}-1\leq 0,$	(6)
	$\displaystyle\phi_{e}\in[0.0,\>\>1.0],\>\>\forall e.$	(7)

Here the objective is to minimize the compliance, measured by the strain energy $c$ . $\bm{K}$ is the stiffness matrix in finite element analysis. $\bm{U}$ is the displacement vector, obtained by solving the static elasticity equation (Eq. 5), where $\bm{F}$ is the loading vector. $g\left(\bm{\phi}\right)$ represents the aggregated local volume constraint.

The formulation takes $\phi_{e}$ as the design variable. The pseudo density field ( $\bm{\rho}$ ) is computed from $\bm{\phi}$ by a density filter ( $\bm{\phi}\to\tilde{\bm{\phi}}$ ), followed by a smoothed Heaviside projection ( $\tilde{\bm{\phi}}\to\bm{\rho}$ ). The density filter, with a filter radius $r_{e}$ smaller than $R_{e}$ (Eq. 1), avoids checkerboard patterns resulting from numerical instabilities. The associated equations of this standard operator are omitted here but can be found in e.g., (Wang et al., 2011; Wu et al., 2018). The purpose of the projection $\tilde{{\bm{\phi}}}\to\bm{\rho}$ is to promote a 0-1 solution, by thresholding at the value of $\frac{1}{2}$ ,

\rho_{e}(\tilde{\phi}_{e})=\frac{\tanh(\frac{\beta}{2})+\tanh(\beta(\tilde{\phi}_{e}-\frac{1}{2}))}{2\,\tanh(\frac{\beta}{2})}.

(8)

The smoothed Heaviside function has a parameter, $\beta$ , to control its sharpness. For improving convergence behaviour, a continuation scheme is applied to gradually increase its sharpness, i.e., we start with $\beta=1$ and double its value every 40 iterations until it reaches 128.

To interpolate the Young’s modulus for intermediate densities, we use the modified SIMP (Solid Isotropic Material with Penalization) model,

E_{e}(\rho_{e})=E_{min}+\rho_{e}^{\gamma}(E_{0}-E_{min}),

(9)

where $E_{0}$ is the Young’s Modulus of a fully solid element. $E_{min}$ is a minimum Young’s modulus ( $E_{min}=1.0e^{-6}E_{0}$ in our test), introduced to avoid the singularity of the global stiffness matrix. $\gamma$ is the penalization factor, which is typically set to 3. $E_{e}(\rho_{e})$ is the interpolated Young’s Modulus of the element with density $\rho_{e}$ .

The optimization problem is solved using the method of moving asymptotes (MMA) (Svanberg, 1987). In all experiments performed in this work, the move limit of design variables is set to 0.01 unless specified otherwise.

Refer to caption — Figure 1: (a) Illustration of the design domain (500x250 simulation elements) and boundary conditions. (b, c) The density distributions after 250 and 1000 iterations, respectively, from topology optimization under local volume constraints. Design parameters are $\alpha_{e}=0.6$ , $R_{e}=18$ and $r_{e}=4.5$ . (d, e, f) The sensitivities at 1000 iterations, of the objective $\frac{\partial c}{\partial{\rho}}$ , the constraint $\frac{\partial g}{\partial{\rho}}$ , and $-\frac{\partial c}{\partial{\rho}}/\frac{\partial g}{\partial{\rho}}$ . (g, h, i) show the optimized density fields under different parameter settings: (g) $R_{e}=12$ . (h) $r_{e}=2.6$ , $R_{e}$ varies linearly from $8$ to $24$ , from the left to right side of the design domain. (i) $\alpha_{e}$ varies linearly from $0.4$ to $0.7$ , from the left to right side of the design domain. All other parameters are kept the same as in (b) and (c).

3 Convergence Analysis and Improvement

When applying topology optimization using local volume constraints, in some scenarios it is observed that the iterative optimization process converges very slowly. When inspecting such scenarios in more detail, for instance, by visualizing the material distribution of the intermediate designs, it turns out that in some regions even after several hundreds or thousands of iterations a distinct binary design cannot be achieved by the optimizer. One of such scenarios is shown in Fig. 1. The rectangular design domain is fixed on its left edge. A unit load is applied on the right, while another unit load on the bottom, both in the middle of the edges. In Fig. 1(b), i.e., optimization after 250 iterations, two large grey regions can be observed. While the grey region on the left converges to a binary design after another 250 iterations, the grey region on the right does not result in a binary design even after a few thousands iterations. Applying a threshold to the intermediate densities to set them to either 0 or 1 results in large void or solid regions with sub-optimal mechanical properties or use of material.

To further analyze the cause of slow convergence, we examine the sensitivities at the 1000 iterations. In Fig. 1(d), the plot of $\frac{\partial c}{\partial{\rho}}$ , the region where low convergence is observed has a high absolute sensitivity, meaning that an increase in density shall be favored for reducing the objective. However, an increase of density in this region will greatly violate the aggregated local volume constraint, as can be seen in Fig. 1(e), the plot of $\frac{\partial g}{\partial{\rho}}$ . Shown in Fig. 1(f) is $-\frac{\partial c}{\partial{\rho}}/\frac{\partial g}{\partial{\rho}}$ , a metric similarly used for deriving a fix-point type update scheme with the optimality criteria (Sigmund, 2001). It can be seen that in the low convergence region the ratio is rather homogeneous and does not indicate a clear material update strategy.

3.1 Relationship between Convergence and Stress

Prior work in infill optimization has shown that the optimized porous structure is in many regions according to the principal stress directions that occur in the solid design domain under equal boundary conditions and external loads. At each point in a 2D solid under load, the stress state is fully described by the stress vectors for two mutually orthogonal orientations. The second-order stress tensor

S(x,y)=\begin{bmatrix}\sigma_{xx}&\tau_{xy}\\ \tau_{xy}&\sigma_{yy}\end{bmatrix}_{(x,y)}

(10)

contains these vectors for the axes of a Cartesian coordinate system. $\sigma_{xx}$ and $\sigma_{yy}$ are the normal stress components along the $x$ and $y$ directions, respectively, $\tau_{xy}$ is the shear stress component.

$S$ is symmetric since the shear stresses given by the off-diagonal elements in $S$ are equal on mutually orthogonal lines. The principal stress directions of the stress tensor indicate the two mutually orthogonal directions along which the shear stresses vanish. These directions are given by the eigenvectors of $S$ , with magnitudes given by the corresponding eigenvalues $\sigma_{1}$ and $\sigma_{2}$ of $S$ . For $\sigma_{1}\geq\sigma_{2}$ , $\sigma_{1}$ is called the major principal stress, and $\sigma_{2}$ the minor principal stress. Accordingly, the corresponding eigenvectors $v_{1}$ and $v_{2}$ are called major and minor principal stress directions. The signs of the principal stress magnitudes classify the stresses into tension (positive sign) or compression (negative sign). However, since there are two principal stresses acting at each point, the classification is with respect to a specific direction.

Fig. 2(a) shows a tensor glyph-based visualization of the stress field, corresponding to the scenario in Fig. 1. Here, the stress tensors are represented by oriented ellipses. The axes of the ellipses are oriented according to the eigenvectors of the stress tensor, and the lengths of their radii are determined by the eigenvalues. The colors of axes indicate the sign of the principal stresses, red for positive and green for negative values. Fig. 1(b), i.e., optimization after 250 iterations, has two large grey regions. Comparing it with Fig. 2(a), it can be seen that the grey region on the left corresponds to $\sigma_{1}=-\sigma_{2}$ , and the one on the right corresponds to $\sigma_{1}=\sigma_{2}$ . While after a few hundred more iterations (see Fig. 1(c)) the grey region on the left converges to a binary design, the grey region on the right shrinks but remains visible. In these regions, the optimizer can favour material growths either along the major or the minor principal stress direction, and it seems that because no preferential direction is present the optimizer has problems to decide for any of them. In the regions where the optimizer doesn’t converge, however, another specific property can be perceived in addition to isotropic stress. As indicated by principal stress lines (PSLs), which are computed by performing numerical integration along the major and minor principal stress directions (see Fig. 2(b)), these regions seem to cover locations where the PSLs indicate directional discontinuities in the tensor field. This observation gives rise to a stress topology-based analysis of the optimization convergence, which we provide in the following.

3.2 Stress Topology-based Analysis

Topology analysis of 2D symmetric second-order tensor fields (e.g., stress tensor fields) has been introduced in the seminal work of Delmarcelle and Hesselink (1994). The topology of a 2D stress tensor field is composed of its degenerate points and the corresponding topological skeleton. At a degenerate point, the stress tensor has repeating eigenvalues, i.e., $\sigma_{1}=\sigma_{2}$ , meaning that the major and minor stress directions cannot be decided. The topological skeleton is given by principal stress lines—so-called separatrices—that start from degenerate points.

An isolated degenerate point can be classified by the winding number of one of the eigenvector fields on a loop surrounding the degenerate point, i.e., the tensor index. The most fundamental types of degenerate points are wedges and trisectors, with a tensor index of $1/2$ and $-1/2$ , respectively. In Section 4, we describe how the degenerate points are computed and classified for a stress tensor field given at the vertices of a Cartesian grid. A major/minor separatrix is a principal stress line starting at a degenerate point and following the major/minor eigenvector field. Let us also refer to Sec. 4 for a discussion of how to determine these directions. Fig. 2c illustrates the major and minor separatrices in the stress field corresponding to the used test scenario in Fig. 1.

At a trisector there are three separatrices in the major (and three in the minor) principal direction field that divide the neighborhood into three sectors sharing this point. Around a wedge there can be either one sector or three sectors. In a sector, similar stress trajectories in the major and minor principal stress direction fields are observed. Fig. 2c, as well as all other experiments we have performed, indicate that the regions where convergence cannot be achieved are always centered around a trisector, while wedges seem to have no influence on the convergence rate. Furthermore, regions where the convergence rate is low but the optimizer can eventually arrive at a stable binary design do not contain any degenerate point (Fig. 1b, left). These regions are characterized by isotropic stress, yet they are not affected by topological changes of the stress field. In Fig. 3, we give two more examples which verify our observations, and clearly indicate the relationships between convergence rate, stress isotropy, and existence of trisectors.

The topological skeletons can be perceived as limits of the principal stress lines close to the boundary between different stress regions. We hypothesize that the porous structures in regions where convergence is not achieved, if a stable binary design should be enforced, follow these skeletons, just as porous structures in other regions follow the principal stress directions. To validate this hypothesis, our idea is to guide the material deposition along the topological skeleton via a skeleton-based initialization of the density field. The initialization sets the optimizer to a state in which a stable binary design in the regions around trisectors is prescribed. The design can be changed during the course of the optimization, yet our experiments demonstrate that the optimizer maintains this design and builds additional support structures around it. These results empirically proof the validity of our hypothesis, and they indicate that the deviation from the prescribed skeleton is not favorable for the objective function.

3.3 Stress Topology-guided Initialization

Typically in density-based topology optimization, the density field is initialized with a constant value. For infill optimization, the constant is chosen as the local volume upper bound. In accordance to the observation that the material layout in porous infill optimization is guided by the principal stress directions, we propose to augment the initialization by setting the densities of elements close to the topological skeleton to a high value. This strategy is fully automatic, since the computation of neither the degenerate points nor the topological skeleton does involve any user intervention.

To generate the initial material layout, the elements which are near the topological skeleton are identified first. In the current implementation, all elements that are touched by any of the PSLs belonging to the skeleton are identified. Then, the initial volume fraction of these elements is set to solid at the beginning of the optimization process. In this way, the initial material field around a trisector degenerate point becomes inhomogeneous, giving rise to sensitivities favoring a unique topology layout. It is worth noting that these pre-embedded solid elements are not passive elements but still belong to the design space, i.e., the material density at these elements can be adjusted by the optimizer if a stiffer design can be achieved. Fig. 4 shows the initial material fields that are used in the test cases ’Cantilever’, ’Bracket’ and ’Bearing’.

The proposed initialization process can be integrated into porous infill optimization in a fully automated way. Once the design domain, material parameters, fixations and external load conditions are given, the following steps are performed:

1)

Finite element analysis to compute the stress field in the fully solid design domain.
2)

Topology analysis including the computation of all trisector degenerate points and the topological skeleton containing these points.
3)

Initialization of the material field according to the topological skeleton.
4)

Topology optimization using local volume constraints for porous infill optimization.

4 Implementation Details

In the following, we discuss the computation of the locations of degenerate points in a given Cartesian simulation grids, as well as the computation of the topological skeleton that is required to initialize the material field. Given the definition of degenerate points, a degenerate point can be located by solving the following system of equations:

\begin{array}[]{l}\sigma_{xx}(x^{*},y^{*})-\sigma_{yy}(x^{*},y^{*})=0,\\ \tau_{xy}(x^{*},y^{*})=0,\end{array}

(11)

where $(x^{*},y^{*})$ denotes the coordinates of the point to be solved for. Here we consider the general situation in topology optimization, i.e., the finite element analysis is performed using axis-aligned quadrilateral finite elements with bilinear shape functions. Thus, each element has four nodes that coincide with the element’s vertices, and the values at the nodes are bilinearly interpolated within the element. Then, Eq. 11 becomes a non-linear system of equations, which can be solved by the Newton-Raphson method.

Since degenerate points usually appear only in a few elements, an efficient way is required to test whether a cell can contain such a point and needs to be further analysed, or can be excluded right away. Therefore, each element is first classified according to the following conditions:

\begin{array}[]{l}\sigma_{xx}(x_{i},y_{i})-\sigma_{yy}(x_{i},y_{i})>0,\quad i=1:4\quad\mathrm{or}\\ \sigma_{xx}(x_{i},y_{i})-\sigma_{yy}(x_{i},y_{i})<0,\quad i=1:4\quad\mathrm{or}\\ \tau_{xy}(x_{i},y_{i})>0,\quad i=1:4\quad\mathrm{or}\\ \tau_{xy}(x_{i},y_{i})<0,\quad i=1:4\end{array}

(12)

where $(x_{i},y_{i}),\>\>i=1:4$ refers to the four nodal coordinates of a finite element. It can be easily shown that an element cannot contain a degenerate point if any of the conditions in Eq. 12 is true. If none of the conditions is true, the element needs to be further analyzed to locate a degenerate point in its interior. Fig. 5 shows a possible distribution of the eigenvalues corresponding to the major and minor principal stress directions in a quadrilateral simulation element containing a degenerate point.

In a symmetric second tensor field, two types of stable degenerate points exist: trisectors and wedges. They are indicated by characteristic patterns of the PSLs in their vicinity, and are determined from the so-called tensor gradients (see Delmarcelle and Hesselink (1994) for a comprehensive derivation). First, the partial derivatives of the tensor are introduced as

\begin{array}[]{l}a=\frac{1}{2}\frac{\partial{(\sigma_{xx}-\sigma_{yy})}}{\partial x}\quad b=\frac{1}{2}\frac{\partial{(\sigma_{xx}-\sigma_{yy})}}{\partial y}\\ c=\frac{\partial{\tau_{xy}}}{\partial x}\quad\quad\quad\quad\>\>d=\frac{\partial{\tau_{xy}}}{\partial y}\end{array}

(13)

These derivatives are then used to compute the invariant under rotation

\delta=ad-bc.

(14)

The sign of $\delta$ determines the type of the degenerate point. I.e., a trisector degenerate point is indicated by $\delta<0$ , and a wedge degenerate point is indicated by $\delta>0$ . At a trisector degenerate point, there are three major and three minor separatrices starting from this point. In contrast, two separatrices start from a wedge, one coincides with the major PSL and the other one with the minor PSL (see Fig. 2c). These separatrices are termed the topological skeleton of a stress tensor field, i.e., the topological skeleton is composed of the PSLs starting from the degenerate points. Compared to the PSLs not belonging to the topological skeleton, the tangent of the topological skeleton at the degenerate point is not unique, since there is an infinite set of principal stress directions at such points. To solve this problem, Delmarcelle and Hesselink (1994), propose that the tangents to the topological skeleton at the degenerate points are the real root(s) of the cubic equation

dx^{3}+(c+2b)x^{2}+(2a-d)x-c=0.

(15)

5 Results and Discussions

In this section, we use several examples to demonstrate the effectiveness of the proposed initialization for density-based porous infill optimization. All design domains are discretized by Cartesian finite element grids with unit size simulation elements. The Young’s Modulus and Poisson’s ratio are set to 1.0 and 0.3, respectively. Convergence improvement is quantified by the sharpness measurement

s=\frac{4}{n}\sum\limits_{e}\rho_{e}(1-\rho_{e})

(16)

A small value of $s$ indicates a sharper binary design of the optimized topology.

Figure 6 shows the binary designs that are generated using an initialization by topological skeletons. The same parameter settings as in Fig. 1 are used here. As can be seen, in all cases a binary design is achieved regardless of the area of the region around the degenerate point where convergence is not achieved by the original approach. Table 1 compares the mechanical properties of the designs generated by both approaches, as well as the used material and the sharpness (cf. Eq. 16) of the designs after 1000 optimization iterations. Notably, even after some thousands of iterations convergence cannot be reached via original porous infill optimization. As can be seen from the sharpness values, the proposed initialization strategy improves the convergence behavior of porous infill optimization considerably. In all test cases, a distinct binary design has been reached within the given number of iterations. The difference in compliance and material fraction is rather small.

Table 1: Quality and convergence comparison. Each pair of rows shows the compliance, the fraction of solid material, and the sharpness of the resulting designs when using porous infill optimization without (top) and with (bottom) stress topology-guided material initialization.

Cases	Compliance	Solid fraction	Sharpness
Fig. 1c	27.99	0.479	$1.8\times 10^{-2}$
Fig. 6a	28.26	0.473	$5.4\times 10^{-3}$
Fig. 1g	27.37	0.501	$1.1\times 10^{-2}$
Fig. 6b	27.63	0.497	$1.1\times 10^{-3}$
Fig. 1h	35.49	0.377	$1.0\times 10^{-2}$
Fig. 6c	35.76	0.379	$5.4\times 10^{-3}$
Fig. 1i	27.98	0.482	$1.2\times 10^{-2}$
Fig. 6d	28.47	0.480	$1.3\times 10^{-3}$
Fig. 3c	19.15	0.490	$1.7\times 10^{-2}$
Fig. 7a	19.32	0.482	$5.8\times 10^{-3}$
Fig. 3f	45.67	0.515	$6.6\times 10^{-2}$
Fig. 7b	45.49	0.504	$4.8\times 10^{-3}$

Figure 7 shows converged density distributions of the ’Bracket’ and ’Bearing’, obtained using the topology-based initialization strategy. Low convergence regions from the original approach (cf. Fig. 3) are removed. The convergence is again confirmed by a reduction in the sharpness value. In the ’Bearing’ result from the original approach (cf. Fig. 3f) the area of the low convergence regions is large. In this case, the sharpness value is reduced by an order of magnitude by the proposed initialization. It leads to a stiffer structure with less material consumption.

We further test the applicability of the proposed initialization on topology optimization with both local and global volume constraints. The global volume constraint is

\frac{1}{n}\sum\limits_{e}\rho_{e}-\alpha_{total}\leq 0.

(17)

The test is performed on a square where its four corners are loaded (Fig. 8(a)), an example taken from Stutz et al. (2020). This example has two trisectors, as shown in Fig. 8(b). From Fig. 8(c), it can be seen that the central region is largely grey after 1000 optimization iterations. The grey region disappears in the optimized result from the proposed initialization (Fig. 8(d)). The significant improvement in convergence can also be seen from the plot of the sharpness over iterations, shown in Fig. 9(right). As the large grey region is replaced by a binary design, the material consumption reduces from 0.400 to 0.378 and the compliance value decreases marginally from 26.02 to 25.96.

Interestingly, under certain design specifications, the original approach is able to create binary designs at the presence of degenerate points. This happens if the specified local volume bound is small. Fig. 10 shows the cantilever example optimized under $\alpha_{e}=0.4$ , with a homogeneous initialization (a) and with the proposed topology-guided initialization (b). The sharpness values of results from both initialization strategies are very close. The one with the proposed initialization leads to a slightly larger compliance ( $103.8\%$ ), while consuming $3.8\%$ less material. Our experiments also revealed that the local convergence may, to some extend, be alleviated by a more aggressive move limit in the MMA solver. Fig. 10(c) shows the optimized result with a move limit of 0.1 (in contrast to a limit of 0.01 in previous examples), under a homogeneous initialization. While a binary design is obtained, the design has irregular large void regions that do not agree with the intention of creating distributed porous infill structures. The introduction of the topology-guided initialization is able to fill the void (see Fig. 10(d)). The latter design consumes $3.6\%$ more material, and decreases the compliance by $4.3\%$ .

6 Conclusions

In this work, we have analyzed the convergence of porous infill optimization towards a stable binary design. In a number of experiments we have shown that low convergence regions may appear in this variant of topology optimization, prohibiting an automatic generation of a distinct and mechanically sound binary design. By analyzing the topology of the stress field that arises in the solid object, the existence of trisector degenerate points in this field could be determined as the major cause of low convergence. Based on this observation, we have proposed an initialization process for porous infill optimization that quickly guides the optimization towards a stable binary design. This process generates an initial solid material layout along the topological skeleton of the stress field, which is comprised of principal stress lines starting at the trisector degenerate points.

In the future, we intend to shed light on the following extensions of the proposed approach: Firstly, we aim to consider the application of stress topology-based porous infill optimization to three-dimensional (3D) domains. Therefore, the convergence of 3D porous infill optimization first needs to be analyzed, using dedicated visualization techniques for 3D scalar fields. Then, since degenerate points become lines and surfaces in 3D, the relationships between the 3D stress field topology and the local convergence ratio needs to be investigated. Based on these investigations, specific initialization strategies and material growth processes need to be developed. Secondly, we will consider stress topology analysis for homogenization-based infill optimization. In particular, we will address the automatic generation of a 2D quad-dominant mesh where the mesh edges align with the principal stress directions. Porous infill optimization, under a single load case, tends to lay out the material along the mutually orthogonal principal stress lines, and—with our proposed initialization—automatically handles the material layout around degenerate points where quad meshing approach have difficulties to construct a consistent mesh structure (Wu et al., 2021b). We will build upon this observation and combine stress topology-guided porous infill optimization with the enforcement of material deposition along the principal stress lines.

Acknowledgements This work was supported in part by a grant from German Research Foundation (DFG) under grant number WE 2754/10-1.

Conflict of interest The authors declare that they have no conflict of interest.

Replication of results All important details have been presented in the paper. To facilitate re-implementation and reuse of this method, we provide a demonstration program¹¹1https://github.com/Junpeng-Wang-TUM/Infill $\_$ plus in Matlab. This demo works with rectangular design domains (i.e., Fig. 1 and Fig. 8). This demo additionally includes the functionality for topology analysis and visualization of the stress tensor field.

References

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, DOI 10.1016/0045-7825(88)90086-2
Cang et al. (2019) Cang R, Yao H, Ren Y (2019) One-shot generation of near-optimal topology through theory-driven machine learning. Computer-Aided Design 109:12–21, DOI 10.1016/j.cad.2018.12.008
Clausen et al. (2015) Clausen A, Aage N, Sigmund O (2015) Topology optimization of coated structures and material interface problems. Computer Methods in Applied Mechanics and Engineering 290:524–541, DOI 10.1016/j.cma.2015.02.011
Das and Sutradhar (2020) Das S, Sutradhar A (2020) Multi-physics topology optimization of functionally graded controllable porous structures: Application to heat dissipating problems. Materials & Design 193:108775, DOI 10.1016/j.matdes.2020.108775
Delmarcelle and Hesselink (1994) Delmarcelle T, Hesselink L (1994) The topology of symmetric, second-order tensor fields. In: Proceedings Visualization’94, IEEE, pp 140–147
Dou (2020) Dou S (2020) A projection approach for topology optimization of porous structures through implicit local volume control. Structural and Multidisciplinary Optimization pp 1–16, DOI 10.1007/s00158-020-02539-x
Guest (2009) Guest JK (2009) Imposing maximum length scale in topology optimization. Structural and Multidisciplinary Optimization 37(5):463–473, DOI 10.1007/s00158-008-0250-7
Kranz et al. (2021) Kranz M, Lüdeker JK, Kriegesmann B (2021) An empirical study on stress-based fail-safe topology optimization and multiple load path design. Structural and Multidisciplinary Optimization pp 1–22, DOI 10.1007/s00158-021-02969-1
Li et al. (2020) Li H, Gao L, Li H, Tong H (2020) Spatial-varying multi-phase infill design using density-based topology optimization. Computer Methods in Applied Mechanics and Engineering 372:113354, DOI 10.1016/j.cma.2020.113354
Li et al. (2021) Li H, Gao L, Li H, Li X, Tong H (2021) Full-scale topology optimization for fiber-reinforced structures with continuous fiber paths. Computer Methods in Applied Mechanics and Engineering 377:113668, DOI 10.1016/j.cma.2021.113668
Liu et al. (2021) Liu Y, Zhou M, Wei C, Lin Z (2021) Topology optimization of self-supporting infill structures. Structural and Multidisciplinary Optimization 63(5):2289–2304, DOI 10.1007/s00158-020-02805-y
Qiu et al. (2020) Qiu W, Jin P, Jin S, Wang C, Xia L, Zhu J, Shi T (2020) An evolutionary design approach to shell-infill structures. Additive Manufacturing 34:101382, DOI 10.1016/j.addma.2020.101382
Schmidt et al. (2019) Schmidt MP, Pedersen CB, Gout C (2019) On structural topology optimization using graded porosity control. Structural and Multidisciplinary Optimization 60(4):1437–1453, DOI 10.1007/s00158-019-02275-x
Sigmund (2001) Sigmund O (2001) A 99 line topology optimization code written in Matlab. Structural and multidisciplinary optimization 21(2):120–127, DOI 10.1007/s001580050176
Stutz et al. (2020) Stutz F, Groen J, Sigmund O, Bærentzen J (2020) Singularity aware de-homogenization for high-resolution topology optimized structures. Structural and Multidisciplinary Optimization 62(5):2279–2295, DOI 10.1007/s00158-020-02681-6
Svanberg (1987) Svanberg K (1987) The method of moving asymptotes—a new method for structural optimization. International journal for numerical methods in engineering 24(2):359–373, DOI 10.1002/nme.1620240207
Träff et al. (2021) Träff EA, Sigmund O, Aage N (2021) Topology optimization of ultra high resolution shell structures. Thin-Walled Structures 160:107349
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, DOI 10.1007/s00158-010-0602-y
Wang et al. (2020) Wang J, Wu J, Westermann R (2020) A globally conforming lattice structure for 2d stress tensor visualization. In: Computer Graphics Forum, Wiley Online Library, vol 39, pp 417–427
Wu et al. (2017) Wu J, Clausen A, Sigmund O (2017) Minimum compliance topology optimization of shell–infill composites for additive manufacturing. Computer Methods in Applied Mechanics and Engineering 326:358–375
Wu et al. (2018) Wu J, Aage N, Westermann R, Sigmund O (2018) Infill optimization for additive manufacturing – approaching bone-like porous structures. IEEE Transactions on Visualization and Computer Graphics 24(2):1127–1140, DOI 10.1109/TVCG.2017.2655523
Wu et al. (2021a) Wu J, Sigmund O, Groen JP (2021a) Topology optimization of multi-scale structures: a review. Structural and Multidisciplinary Optimization pp 1–26, DOI 10.1007/s00158-021-02881-8
Wu et al. (2021b) Wu J, Wang W, Gao X (2021b) Design and optimization of conforming lattice structures. IEEE Transactions on Visualization and Computer Graphics 27(1):43 – 56, DOI 10.1109/TVCG.2019.2938946
Zhao and Zhang (2021) Zhao Z, Zhang XS (2021) Design of graded porous bone-like structures via a multi-material topology optimization approach. Structural and Multidisciplinary Optimization 64(8):677––698, DOI 10.1007/s00158-021-02870-x