Multi-Material Topology Optimization on Separate Tetrahedral Meshes with Explicit Design Resolution by means of Remeshing

1. Introduction

Topology optimization is a method that can support the product developer in the lightweight design of structures and is therefore a central tool in product development to reduce $C{O}_2$ emissions from products. The most widespread method for topology optimization is density-based topology optimization, the development of which began with the publications by Bendsøe and Kikuchi [1] and Bendsøe [2]. In addition to the density-based topology optimization methods, there are many other topology optimization methods. These include, for example, the Optimization Criteria-based methods (), the Phase Field methods [3] or the nature-mimicking Evolutionary Structural Optimization [4,5]. Another topology optimization method is based on the level set method and was primarily developed by Sethian and Wiegmann [6], Wang et al. [7] and Allaire et al. [8]. In this method, the geometry is described during optimization using a discretized level set function. The geometry is evolved by applying the classical shape derivative to determine the front velocity and then by solving a Hamilton-Jacobi equation. [8]

In addition to the level set-based topology optimization approach just described, there are other variants of using the level set method for topology optimization. An overview is given in [9]. One such method is the Velocity Field Level Set method (VFLS) [10]. It allows the use of mathematical programming methods, such as the Method of Moving Asymptotes (MMA) [11] in combination with the level set method. Therby enabling the straightforward integration of several constraints into level set-based topology optimization. [10]

Wang and Wang [12] developed the so-called Color Level Set method for the optimization of multi-material systems. For this purpose, they refined a description of multi-material systems known from image processing for level set topology optimization. It allows the description of $2^n$ materials using n level set functions.

An alternative approach to optimizing multi-material systems is the Reconciled Level Set method (RCLS) [13], in which each material is described by its own level set function. The individual level set functions are evolved independently of each other and any overlaps that arise are subsequently resolved using a correction process [13]. Chen et al. [14] applied the RCLS in the context of multi-material topology optimization of smart energy harvesters and Vogiatzis et al. [15] used it for topology optimization of multi-material negative Poisson’s ratio metamaterials. Tian et al. [16] combined the RCLS with the extended level set method to design ferromagnetic soft active structures using topology optimization.

There are three different approaches to represent the current design for the evaluation of the system response during level set-based topology optimization. In the first approach, a density field of a constant background Finite Element (FE) mesh is adjusted based on the values of the level set function. A smoothed heaviside function is usually used for the adjustment. This variant is similar to density-based topology optimization methods [8]. The second approach applies immersed finite element methods, such as XFEM [17] or CutFEM [18] methods. These use a constant background mesh that is cut by the level set functions. The third approach is to adapt a mesh or remesh the level set function in each iteration [19]. With this approach, the standard Finite Element Method (FEM) can be used.

Allaire et al. [20] developed a level set topology optimization with iterative mesh adaptation. The unstructured mesh is iteratively adapted using classical mesh modification techniques such as edge collapse, vertex smoothing and edge swap. Li et al. [21] combined this mesh adaptation algorithm with the reaction-diffusion equation based level-set method. It uses the topological derivative and a reaction-diffusion equation instead of the shape derivative and the Hamilton-Jacobi equation [22]. Feppon [23] optimized a heat exchanger using also this iterative mesh adaptation approach.

In this article, a multi-material topology optimization method based on RCLS and VFLS with iterative remeshing of the multi-material design is developed. The developed method is validated by two cases. In the developed method, the edge length of the FE mesh can be adjusted. The influence of the edge length and an additional filter radius parameter on the optimized design is investigated in two studies.

2. Methods

In the following, the boundary value problem of elastostatics, which is used in both single and multi-material cases, is introduced first. For the validation of the developed method, the single material optimization problem and the VFLS are then presented. In the next section, the multi-material optimization problem and the RCLS are introduced. Finally, an overview of the developed optimization method is given and its individual components are explained in detail.

2.1. Boundary Value Problem of Elastostatics

This article considers the quasi-static boundary value problem of elastostatics:

$$\begin{array}{ccc}\hfill \mathbf{\nabla }·\mathit{\sigma }-{\mathit{f}}_b& =\mathbf{0}\hfill & \hfill \mathrm{in}\mathrm{\Omega },\\ \hfill \mathit{u}& ={\mathit{u}}_D\hfill & \hfill \mathrm{on}{\xi }_D,\\ \hfill \mathit{\sigma }\mathit{n}& ={\mathit{f}}_t\hfill & \hfill \mathrm{on}{\xi }_N.\end{array}$$

(1)

$\mathit{u}$ is the displacement field to be identified and $\mathit{\sigma }$ is the stress field of the domain $\mathrm{\Omega }$. For the Neumann boundary ${\xi }_n$ with normal vector $\mathit{n}$, the force ${\mathit{f}}_{\mathit{t}}$ is specified as a boundary condition. In addition, in this article the volume forces are assumed to be ${\mathit{f}}_b=\mathbf{0}$ and the displacements ${\mathit{u}}_D=\mathbf{0}$ on the Dirichlet boundary ${\xi }_D$.

To solve the boundary value problem 1, an additional constitutive law is required. For this purpose, Hook’s linear-elastic law is used:

$$\begin{array}{c}\hfill \mathit{\sigma }=\tilde{\mathbb{C}}\left[\mathit{\epsilon }\left(\mathit{u}\right)\right].\end{array}$$

(2)

The strain tensor $\mathit{\epsilon }\left(\mathit{u}\right)$ is given by $\mathit{\epsilon }\left(\mathit{u}\right)={\displaystyle \frac{1}{2}}(\nabla \mathit{u}+{(\nabla \mathit{u})}^T)$. In Equation (2), the elasticity tensor $\tilde{\mathbb{C}}$ is a function that can be described in the single-material case by $\tilde{\mathbb{C}}={\chi }_1\left(\mathit{x}\right){\mathbb{C}}_1$ and in the multi-material case by $\tilde{\mathbb{C}}={\sum }_i{\chi }_i\left(\mathit{x}\right){\mathbb{C}}_i$. The indicator functions ${\chi }_i\left(\mathit{x}\right)$ of the i-th material each have the value 1 if the i-th material is present and 0 otherwise.

Using the constitutive law and the above simplification, the weak form of the boundary value problem can be derived:

$${\int }_{\mathrm{\Omega }}\mathit{\epsilon }\left(\mathit{u}\right)·\tilde{\mathbb{C}}\left[\mathit{\epsilon }\left(\mathit{\psi }\right)\right]d\mathrm{\Omega }={\int }_{{\xi }_N}{\mathit{f}}_{\mathit{t}}·\mathit{\psi }d\xi \forall \mathit{\psi }\in U_{ad}.$$

(3)

Here $U_{ad}$ denotes the space of the kinematically admissible displacement fields.

If one also introduces the two abbreviations

$$\begin{array}{cc}\hfill \mathit{a}(\mathit{u},\mathit{\psi })& ={\int }_{\mathrm{\Omega }}\mathit{\epsilon }\left(\mathit{u}\right)·\tilde{\mathbb{C}}\left[\mathit{\epsilon }\left(\mathit{\psi }\right)\right]d\mathrm{\Omega },\hfill \\ \hfill \mathit{l}(\mathit{u},\mathit{\psi })& ={\int }_{{\xi }_N}{\mathit{f}}_{\mathit{t}}·\mathit{\psi }d\xi ,\hfill \end{array}$$

(4)

Equation (3) can be written as $\mathit{a}(\mathit{u},\mathit{\psi })=\mathit{l}\left(\mathit{\psi }\right)$.

2.2. Single-Material

2.2.1. Optimization Problem

The optimization problem for the single material case is:

$$\begin{array}{cc}\hfill & \underset{\eta }{min}{w}_c{\int }_{\mathrm{\Omega }}{\displaystyle \frac{1}{2}}\mathit{\epsilon }\left(\mathit{u}\right)·{\mathbb{C}}_1\left[\mathit{\epsilon }\left(\mathit{u}\right)\right]d\mathrm{\Omega }+w_r{\int }_{\partial \mathrm{\Omega }}d\xi \hfill \\ \hfill \mathrm{s}.\mathrm{t}.& a(\mathit{u},\mathit{\psi })-l\left(\mathit{\psi }\right)=0,\forall \mathit{\psi }\in U_{ad}\hfill \\ \hfill & f_{V1}\le {\overline{f}}_{V1}\hfill \\ \hfill & v^{min}\le {\eta }_i\le v^{max},i=1,...,n_d\hfill \end{array}$$

(5)

In Equation (5), the weighted objective function of compliance and regularization term is considered, with the two weighting factors $w_c$ and $w_r$. The latter ensures that a strut-like structure is created and also regularizes the optimization problem [8]. The boundary value problem 3 in weak form exists as an equality constraint and the displacements $\mathit{u}\left(\mathit{\eta }\right)$ must satisfy this constraint at all times. The target volume fraction ${\overline{f}}_{V1}$ of the material in relation to the volume of the design space $V_D$ is specified as an inequality constraint. The design variables ${\eta }_i$ are also constrained by the lower and upper bounds $v^{min}$ and $v^{max}$. In the case of this article, the bounds are normal velocities, as the design variables in the VFLS are the normal velocities.

2.2.2. Velocity Field Level Set Method

The description of the part surface during optimization is done using the level set function $\phi :{\mathbb{R}}^3\to \mathbb{R}$. An example of a level set function is depicted in Figure 1. The level set function is positive inside the part, zero on the part surface and negative outside the part (in the rest of the design domain). In summary, the following applies:

$$\left\{\begin{array}{cc}\phi \left(\mathit{x}\right)\ge 0,\hfill & \forall \mathit{x}\in \mathrm{\Omega }\left(\mathrm{part}\right)\hfill \\ \phi \left(\mathit{x}\right)=0,\hfill & \forall \mathit{x}\in \partial \mathrm{\Omega }(\mathrm{part}\mathrm{surface})\hfill \\ \phi \left(\mathit{x}\right)<0,\hfill & \forall \mathit{x}\in D\setminus \mathrm{\Omega }\hfill \end{array}\right.$$

(6)

In this article, a special form of the level set function, the so-called Signed Distance Function (SDF), is used to numerically stabilize the evolution of the surface described below [8]. The function values of the SDF are the signed distances to the surface inside and outside the part. For example, in the part, $\phi \left(\mathit{x}\right)=d(\mathit{x},\partial \mathrm{\Omega })$ applies, whereby the function d returns the distance between $\mathit{x}$ and $\partial \mathrm{\Omega }$. SDFs also have the property that $\parallel \nabla \phi \left(\mathit{x}\right)=1\parallel$. [24]

The evolution of the level set function and thus the surface ($\phi \left(\mathit{x}\right)=0$) in a previously given or calculated normal velocity field $v^n$ is described by the Hamilton-Jacobi equation [24]:

$${\displaystyle \frac{\partial \phi }{\partial t}}+v^n\parallel \nabla \phi \parallel =0.$$

(7)

The normal velocity field $v^n$ is linked to the design variables $\mathit{\eta }$ by linear interpolation in the topology optimization using VFLS [10]:

$$v^n\left(\mathit{x}\right)=\sum _{i=1}^{n_d}{b}_i\left(\mathit{x}\right){\eta }_i.$$

(8)

Here, $b_i\left(\mathit{x}\right)$ is the global basis function of the linear FEM associated with node i and ${\eta }_i$ is the respective design variable. When discretizing the normal velocity field, other basis functions can be used in addition to the basis functions from the linear FEM. An overview of this is given in [10].

The design variables are varied by the optimization algorithm based on previously calculated sensitivities. In the VFLS method, the derivation of the sensitivities is based on the shape derivative [10]. The shape derivative of the compliance component of the objective function is based on the simplification described above (no volume forces and no variation of the Neumann and Dirichlet boundary) [8,10]:

$$D{J}_c\left(\mathrm{\Omega }\right)\left(\mathit{\theta }\right)={\int }_{\partial \mathrm{\Omega }}\mathit{\epsilon }\left(\mathit{u}\right)·{\mathbb{C}}_1\left[\mathit{\epsilon }\left(\mathit{u}\right)\right]\mathit{\theta }·\mathit{n}d\xi ,$$

(9)

and the shape derivative of the regularization component [8]:

$$D{J}_r\left(\mathrm{\Omega }\right)\left(\mathit{\theta }\right)={\int }_{\partial \mathrm{\Omega }}\kappa \mathit{\theta }·\mathit{n}d\xi .$$

(10)

The mean curvature $\kappa$ of the surface $\partial \mathrm{\Omega }$ is defined as $\kappa =\nabla ·\mathit{n}$ [8]. Using $v^n=\mathit{\theta }·\mathit{n}$, the derivative of the objective function with respect to the design variables ${\eta }_i$ is given by [10]:

$${\displaystyle \frac{\partial J}{\partial {\eta }_i}}=w_c{\int }_{\partial \mathrm{\Omega }}\mathit{\epsilon }\left(\mathit{u}\right)·{\mathbb{C}}_1\left[\mathit{\epsilon }\left(\mathit{u}\right)\right]b_{i}d\xi +w_r{\int }_{\partial \mathrm{\Omega }}\kappa b_{i}d\xi ,i=1,...,n_d.$$

(11)

Similarly, the following derivative for the volume constraint can be derived [10]:

$${\displaystyle \frac{\partial f_{V1}}{\partial {\eta }_i}}=-{\displaystyle \frac{1}{V_D}}{\int }_{\partial \mathrm{\Omega }}{b}_{i}d\xi ,i=1,...,n_d.$$

(12)

2.3. Multi-Material

2.3.1. Optimization Problem

The optimization problem for the multi-material case with two materials is as follows:

$$\begin{array}{cc}\hfill & \underset{\mathit{\eta }}{min}{w}_c{\int }_{\mathrm{\Omega }}{\displaystyle \frac{1}{2}}\mathit{\epsilon }\left(\mathit{u}\right)·\tilde{\mathbb{C}}\left[\mathit{\epsilon }\left(\mathit{u}\right)\right]d\mathrm{\Omega }+w_{ri}{\int }_{\mathrm{\Gamma }}d\mathrm{\Gamma }+w_{rs}{\int }_{\partial \mathrm{\Omega }}d\xi \hfill \\ \hfill \mathrm{s}.\mathrm{t}.& a(\mathit{u},\mathit{\psi })-l\left(\mathit{\psi }\right)=0,\forall \mathit{\psi }\in U_{ad}\hfill \\ \hfill & f_{V1}\le {\overline{f}}_{V1}\hfill \\ \hfill & f_{V2}\le {\overline{f}}_{V2}\hfill \\ \hfill & v^{min}\le {\eta }_i\le v^{max},i=1,...,n_d\hfill \end{array}$$

(13)

For the regularization in the multi-material case, the area of the interface $\mathrm{\Gamma }$ between the two materials and their surface $\partial \mathrm{\Omega }$ is used. The former regularization ensures that the optimization problem is solvable and the latter ensures strut-like structures in the result. The weighting is performed analogously to the single-material case with the weighting factors $w_c$, $w_{ri}$ and $w_{rs}$. In addition, there is another volume constraint for the volume fraction of the second material $f_{V2}$ in the multi-material case compared to the single-material case.

2.3.2. Reconciled Level Set Method

This article applies the RCLS [13] to describe the multi-material case. Each material is represented with its own level set function in the RCLS. Here, these are the two level set functions ${\phi }_1$ and ${\phi }_2$ (see Figure 2).

The two level set functions are developed separately using the Hamilton-Jacobi Equation (7). The sensitivities from the single material case are used to determine the normal velocity fields required for this (see Equations (11) and (12)). This is feasible because the RCLS is coupled with the VFLS method in this article. However, the elasticity tensor associated with the material (or the level set function) must be used for the compliance component of the sensitivity (e.g. for ${\phi }_2\left(\mathit{x}\right)$${\mathbb{C}}_2$).

In the multi-material case, the design variable vector takes the form $\mathit{\eta }={[{\eta }_1^1,...,{\eta }_{n_d}^1,{\eta }_1^2,...,{\eta }_{n_d}^2]}^T$.

Since the two level set functions are developed independently in RCLS, the individual level set functions may overlap during evolution (see Figure 3). However, this problem can be solved using the Merriman-Bence-Osher (MBO) operator [13]:

$${\phi }_i={\displaystyle \frac{1}{2}}({\hat{\phi }}_i-\underset{i\ne j}{max}{\hat{\phi }}_j).$$

(14)

This means that the RCLS works like a predictor-corrector method. In the first step, the level set functions are further evolved in the respective existing velocity field without considering the overlap (predictor step). This results in the temporary level set functions ${\hat{\phi }}_1$ and ${\hat{\phi }}_2$ (predictor step). The overlaps are then removed using Equation (14) (correction step). The described method is also depicted in Figure 3.

2.4. Optimization Method

The optimization loop of the developed method is shown in Figure 4. At the beginning of each iteration, the current design, represented by the level set functions on the level set mesh, is remeshed with tetrahedra using the method described below in Section 2.4.5. The result of the remeshing is the FE mesh of the current design. This is then used to determine the system response using linear FEM and to calculate the sensitivities (see Section 2.2.2). The calculated sensitivities serve as input for the MMA [11], which updates the design variables (normal velocities). The implementation of the MMA follows [25]. Since the normal velocities are only non-zero in the immediate vicinity (on the surface) of the level set functions, they are extrapolated to the entire design domain (constant level set mesh) using the velocity extension method described in Section 2.4.1. To stabilize the optimization method, the extrapolated normal velocities are then filtered using the Helmholtz equation (see Section 2.4.2). The filtered normal velocities are then used to evolve the level set functions (see Section 2.4.3). The evolved level set functions are then corrected using the MBO operator and subsequently their SDF property is restored using redistancing (see Section 2.4.4).

2.4.1. Velocity Extension

The normal velocity field on the entire level set mesh is required for each level set function to be evolved. However, since the sensitivities are only non-zero in the vicinity of the zero level sets (surfaces) and consequently the normal velocities (design variables) updated by the optimizer are also non-zero, these must be extrapolated to the rest of the design domain (LS mesh).

For this purpose, a method is developed in this article with an analogous approach to the redistancing method described in [26]. Extrapolation is carried out in the developed method in a two-stage process, which is outlined in Figure 5. In the initialization phase, the surface described by the zero level set ($\phi \left(\mathit{x}\right)=0$) is sampled. For this, the surface elements (triangles) belonging to the zero level set are derived from the FE mesh using the material ids. The values of the normal velocities are then determined for each node of the surface mesh using linear interpolation with the FE basis functions from the LS mesh.

The result of the initialization phase is a point cloud of the surface $C_s=\{{\mathit{x}}_1,{\mathit{x}}_2,...,{\mathit{x}}_{n_s}\}\subset {\mathbb{R}}^3$ analogous to [26] with interpolated values and their adjacent triangles. The point cloud is also used to create a search data structure (kd tree).

In the extrapolation phase, the projection of each vertex of the level set mesh (design domain) is identified on the discretized level set surface mesh. The search data structure and the list of adjacency triangles are used for this purpose. At the identified point, the normal velocity is interpolated using linear interpolation and assigned to the original node in the level set mesh.

The two steps described can be parallelized straight-forward by means of multi-threading, since both the interpolation step in the initialization phase for each vertex of the surface mesh and the extrapolation step for each vertex of the LS mesh can be carried out independently of each other.

Figure 5. Schematic representation of the velocity extension from the zero level set (blue) to the target vertex based on [27]. For this purpose, the pole of the target vertex is determined on the zero level set and then the normal velocity of the target vertex (orange) is calculated from the normal velocities of the triangle (green) using linear interpolation.

Figure 5. Schematic representation of the velocity extension from the zero level set (blue) to the target vertex based on [27]. For this purpose, the pole of the target vertex is determined on the zero level set and then the normal velocity of the target vertex (orange) is calculated from the normal velocities of the triangle (green) using linear interpolation.

2.4.2. Velocity Filtering

After the normal velocities have been extrapolated into the entire design domain, a filter step is then performed. The filter is applied separately for the two normal velocity fields. The aim of filtering is to smooth the velocity field in order to stabilize the optimization process. The filter used is the one described in [28], which is based on the Helmholtz-PDE and is well known in the field of homogenization-based topology optimization. The filter equation (Helmholtz PDE) is [28]:

$$-r_f^2{\nabla }^2{\check{v}}^n+{\check{v}}^n=v^n.$$

(15)

${\check{v}}^n$ is the filtered normal velocity field that results as the solution of Equation (15). The filter parameter $r_f$ is a length parameter that describes the size of the filter mask, but not its explicit radius [28].

To solve Equation (15), its weak form is generated. The building blocks from the evaluation of the system response are used for evaluation the weak form. The weak form is solved on the LS mesh with homogeneous Neumann boundary conditions ${\displaystyle \frac{\partial {\check{v}}^n}{\partial \mathit{n}}}=0$ and the unfiltered normal velocity field $v^n$ as the right-hand side [28].

The advantages of this filter approach over explicit filters are the reuse of the building blocks including their parallelization and the automatic consideration of non-convex design spaces (holes, undercuts) [28].

2.4.3. Advection

In this section, the solution of the Hamilton-Jacobi Equation (7) on tetrahedral meshes using the finite difference method is presented. The method described in [29] is used for this purpose.

The Hamilton-Jacobi Equation (7) can be written as [29]:

$${\displaystyle \frac{\partial \phi }{\partial t}}+H(\nabla \phi ,\mathit{x})=0\mathrm{with}H(\nabla \phi ,\mathit{x})=v^n\parallel \nabla \phi \left(\mathit{x}\right)\parallel .$$

(16)

The numerical solution of the equation using finite differences requires a stable and accurate upwind approximation of the Hamiltonian $H(\nabla \phi ,\mathit{x})$ [29]. In the method for tetrahedral meshes developed in [29], the approximation using an implicit mesh dual to the level set mesh (cf. Figure 6) is applied as follows:

$$H_{+\alpha }^u=v_{\alpha }{\displaystyle \frac{{\sum }_{c\in \alpha }{\gamma }_c\parallel \nabla {\phi }_{z\left(c\right)}\parallel }{{\sum }_{c\in \alpha }{\gamma }_c}}.$$

(17)

Here, the upwind approximation $H_{+\alpha }^u$ is formed by weighting the cell-centered values $H_z$ surrounding the node $\alpha$ with the corner weights ${\gamma }_c$. These are calculated as [29]:

$${\gamma }_c=max(0,{\displaystyle \frac{{\mathit{v}}_c}{\parallel {\mathit{v}}_c\parallel }}·{\mathit{a}}_c).$$

(18)

In this case, ${\mathit{a}}_c$ are the cell corner unit normals, which are derived by summation and normalization from the outward corner surface area normals ${\mathit{A}}_c={\sum }_{i\in c}{\mathit{A}}_i=-{\sum }_{i\in c}{\mathit{S}}_i$. The corner velocities in the case of normal advection are given as ${\mathit{v}}_c=v_{\alpha }{\displaystyle \frac{\nabla {\phi }_{z\left(c\right)}}{\parallel \nabla {\phi }_{z\left(c\right)}\parallel }}$. For the calculation of the corner velocities and the upwind approximation, the cell centered gradients are still required:

$$\nabla {\phi }_z={\displaystyle \frac{1}{V_z}}\sum _{i\in z}{\mathit{A}}_{\mathit{i}}{\phi }_{p\left(i\right)}.$$

(19)

Here $V_z$ denotes the volume of the current tetrahedron of the level set mesh, $i\in z$ means that all surfaces of the tetrahedron are summed up and ${\phi }_{p\left(i\right)}$ are the level set values associated with the surfaces ${\mathit{A}}_i$ [29].

The integration in time is carried out using the two-stage Runge-Kutta method. It is referred to as two-stage because the intermediate step $n+{\displaystyle \frac{1}{2}}$ is calculated for the calculation of the $n+1$ time step. [29]

$$\begin{array}{cc}\hfill {\phi }_{\alpha }^{n+{\displaystyle \frac{1}{2}}}& ={\phi }_{\alpha }^n-{\displaystyle \frac{\mathsf{\Delta }t}{2}}{H}_{+\alpha }^n,\hfill \\ \hfill {\phi }_{\alpha }^{n+1}& ={\phi }_{\alpha }^n-\mathsf{\Delta }t{H}_{+\alpha }^{n+{\displaystyle \frac{1}{2}}}.\hfill \end{array}$$

(20)

2.4.4. Redistancing

In the progress of the optimization, the signed distance property ($\parallel \nabla \phi \left(\mathit{x}\right)=1\parallel$) of the level set functions disappears due to the evolution by means of the Hamilton-Jacobi Equation (7) and after application of the MBO operator 14. However, this is essential for the stability of the optimization process, which is why it is restored by means of redistancing [8,14].

This article uses a two-step explicit method for redistancing. In the first step, the level set function is discretized using the Marching Tetrahedra (MT) method. For this purpose, all tetrahedra of the level set mesh are looped. For each tetrahedron, the intersection points of the edges with the zero level set function are determined and the respective MT case is identified based on the edge intersection pattern (see Figure 7) [30].

In the second step, the closest point on the discretized surfaces and the distance to this point are determined for each vertex of the level set mesh. The distance calculation is accelerated using an AABB tree. In addition to the distance, the sign is required, which is determined using the Fast-Windig-Number [31].

2.4.5. Remeshing

This section describes the remeshing method and its implementation for resolving the current design in each iteration. A copy of the LS mesh serves as the starting point. Analogous to the redistancing method, all tetrahedra are looped and split based on the MT cases (see Figure 7) and the individual tetrahedra are assigned the respective material id. The resulting mesh is then improved using the method described in [32]. The method consists of the following four iterative phases [32]:

• split all edges with $l_e>l_{e_{max}}$
• collapse edges $l_e<l_{e_{min}}$
• swap edges of elements with $q<q_{min}$
• smooth all vertices

In the first phase, edges with a length greater than $l_{e_{max}}={\displaystyle \frac{4l_{tar}}{3}}$ are split by creating an additional vertex in the middle of the edge to be split. All tetrahedra adjacent to the updated edge are updated. In the second phase, all edges with a length less than $l_{e_{min}}={\displaystyle \frac{4l_{tar}}{5}}$ are collected and then all permissible collapses are performed. The first two phases of the remeshing method primarily ensure that the edge length converges in towards the target edge length $l_{tar}$. In the third phase, all tetrahedra whose quality q is less than $q_{min}$ are identified. [32] The volume-length measure [33] is applied as quality measure. For the adjacent edges, all possible edge swaps are tested and the one with the best possible quality is performed. The edge swap phase is primarily used to improve the tetrahedral quality. In the fourth phase, all vertices are repositioned in a hierarchical smoothing. For this purpose, the vertices are divided into volume and boundary vertices and successively smoothed. The classification of the vertices is explained in more detail below. [32] Smoothing the boundary vertices involves a projection step onto the surface, where we project onto the level set functions using Newton’s method instead of onto an Moving Least Square based geometric definition as in [32].

In the phases listed above, many operations occur that are not permitted because they are not feature-preserving or because inverted tetrahedra are created during execution. For the former, a concise set of rules was established in [32], which is based on a classification of vertices, edges and tetrahedra as well as associated conditions. The classification uses the number of labels of the adjacent tetrahedra (here: #material ids). For vertices, for example, if there is exactly one material in the adjacent tetrahedra, it is a volume vertex and if there is more than one material, it is a boundary vertex (see Figure 8a). The classification then allows rules to be defined that describe which mesh modifications are permitted. This is illustrated by the example of an edge collapse in Figure 8b,c [32].

The four phases described above are followed by a quality improvement phase, in which the last remaining tetrahedra with poor quality are eliminated [32]. At this point, we also use the method described in [34], which classifies the bad tetrahedra and then applies targeted operations depending on the type of tetrahedron.

For the implementation of the mesh modification (edge collapses, edge swap, ...) we use the unique cavity operator from [35]. This allows the mesh modifications required for the phases described above to be represented using the same operator and therefore the same code. This can significantly reduce the implementation and maintenance effort for the mesh modification [35].

The different operations can be implemented by different initialization and subsequent filling of the cavity [35]. Figure 9 shows an edge collapse and an edge swap using a 2d triangle mesh. For the edge collapse $A\to B$ (Figure 9a), all elements (tetrahedrons or triangles) are added to the cavity adjacent to vertex A for the initialization, after which they are all deleted. The adapted mesh is obtained by connecting all surfaces of the cavity that do not contain vertex B with target vertex B.

In order to obtain a robust remeshing method, the first step is to test whether inverted tetrahedra are created during the filling process of the cavity operator. The filter-based orientation predicate described in [36] is used for this purpose. Secondly, when initially cutting the mesh with the level set functions based on the MT cases, a vertex can be zero, negative or positive instead of only positive or negative.

A disadvantage of the cavity operator is the somewhat slower performance compared to explicitly implemented mesh operators. However, this disadvantage can be remedied by thread-based parallelization. In this paper, the vertex smoothing and the edge-collapse phase are parallelized. For the vertex smoothing phase, a heuristic graph coloring is used to identify independent cavities. The vertices of each color are then processed using #threads. After the vertices of a color have been smoothed, the threads are synchronized [37].

The edge collapse phase is parallelized using the producer consumer pattern. We use #threads - 1 producer and 1 consumer. After the edges with too short a length have been identified, the producers test all the preconditions described above and then attempt to perform the collapse virtually. If a collapse is valid, it is entered in a queue. The consumer removes the pre-selected collapses and executes them with renewed preliminary checks. The latter is necessary, as the mesh in the vicinity of the edge may have changed in the time between insertion into the queue and execution of the collapse due to other collapses executed by the worker, meaning that the collapse may no longer be valid. The parallelization strategy is based on the fact that only about $10\%$ of the edge collapses per iteration are valid and can be executed.

The cavity operator described only works inside the volume and not at the boundary. Since no extra operator is implemented for the boundary region and the method described requires a convex region, a mesh preparation step is necessary [32]. In the first step, the convex hull of the region is determined and an offset is then formed from this using the method described in [38]. The resulting volume is filled with tetrahedra. We use the fTetWild algorithm [39] for this. This preparation now makes it possible to apply the cavity operator described above to the entire area, as the filled offset means that the surface tetrahedra are surrounded by elements (see Figure 10). [32]

3. Numerical Experiments

The following section first presents the two cases that are used in the article to illustrate and validate the developed optimization method. In addition, a summary of the parameters that apply to all numerical examples is then provided. The parameters apply for as long as not directly stated otherwise. This is followed by two single-material optimization examples to validate the single-material case and then two multi-material examples to validated the multi-material case.

3.1. Cases

The first case is a cube with side length 80 with chamfered edges on the top (see Figure 11). The case serves as an example for almost arbitrary design spaces of the optimization method. However, it can be compared at the same time with the classic example of the box without chamfering from the literature. The box is supported at four bearing points on the underside. The displacement of the vertices of one of these supports is restricted in all three spatial directions ($u_x=u_y=u_z=0$). The vertices of the other supports can move freely in the plane, but not in the z-direction ($u_z=0$). The square load application and constraint surfaces are slightly protruding. The first element layer adjacent to the surfaces is always considered as a non-design space in the optimization and cannot be modified.

The second case is a cantilever clamped on one side (see Figure 12), which is subjected to a force F in the negative z-direction on the front surface. Two square surfaces with a side length of 15 length units are used for clamping. The vertices of the upper surface are completely fixed in all three spatial directions ($u_x=u_y=u_z=0$) and the vertices of the lower surface are only fixed in the x and y directions ($u_x=u_y=0$). Only half of the design domain is considered with an additional symmetry boundary condition so that the result of the optimization is a symmetrical design. The dimensions of the design domain can be taken from Figure 12. As with the box, the load application and constraint surfaces are slightly protruding and the first element layer is a non-design space.

For both the single and multi-material cases, $v_{min}=-0.8$ and $v_{max}=0.8$ are used. A normal velocity of $-0.2$ is specified for the non-design spaces described above. The filter parameter is selected as $1.25$ times the average edge length and the target edge length is $l_{tar}=2.0$. The target volume fraction is set to $15\%$ of the design space. In the multi-material case, this is divided equally between the two materials so that $f_{V1}\le 0.075$ and $f_{V2}\le 0.075$ applies.

To determine the convergence of the optimization, a convergence criterion is evaluated after the optimization has been performed. This takes into account the 5 previous iterations in the sense of a moving average and uses this filtered value to calculate the percentage change between two points. In the single material case, the limit for the percentage change is set at $1e^{-4}$ and $1e^{-3}$ in the multi-material case.

3.2. Single-Material

For the two single material cases, we use the parameters described below in addition to those listed above. A Young’s modulus of $E=210,000$ and $\nu =0.3$ is used for the material. The weighting is selected as $w_c=0.6$, $w_r=0.4$ for the box case and $w_c=0.5$, $w_r=0.5$ for the cantilever. The regularization term is activated after iteration 30 to ensure that it does not dominate the optimization at the beginning.. The level set function is redistanced every 5 iterations. The c parameter of the MMA algorithm is increased from a starting value of 1.0 to 100.0 within 150 iterations by means of a linear ramp. The latter results in a slow and smooth volume reduction, which gives the optimization method sufficient time to shape the design according to the objective function.

Figure 13a shows the initial hole structure used to initiate the level set function and Figure 13b shows the corresponding optimization result at convergence at iteration 107 of the box in the single material case. It is evident that four main struts are formed to transfer the loads during the optimization process. These transfer the load directly from the point of application to the support points. In addition, cross struts develop in the form of a cross, which prevent the three floating supports from sliding sideways. The optimization result is in line with the state of the research (see e.g. Shape optimization book). However, it should be noted that depending on the optimization parameter (e.g. the c parameter of the MMA algorithm), the optimization sometimes reaches a different local optimum. This results in a rectangular cross strut structure between the four main struts instead of the cross-shaped cross strut structure.

Figure 14a,b show the filtered history of the objective function and the volume constraint respectively. In order to filter out the numerical fluctuation introduced by the remeshing process, a moving-average filter is applied that takes into account the five past iterations. The history of the volume constraint shows the uniform decrease in volume over time.

Figure 15a shows the initial hole structure of the cantilever in the single-material case and Figure 15b the resulting optimization result at iteration 60. The corresponding curves of the filtered objective function and volume constraint are shown in Figure 16a resp. b.

The optimization result shows the usual truss structure consisting of triangles and agrees with the results from the state of research (see [8]).

In summary, it can be said that the developed optimization method provides plausible results in the single material case, which are also consistent with the results from the current state of research. As a result, the developed method is considered validated for the single material case.

3.3. Multi-Material

Two isotropic materials with Young’s modulus of $E_1=70,000$ and $E_2=35,000$ are used for the multi-material cases. The Poisson’s ratio for the two materials is ${\nu }_1={\nu }_2=0.3$. The weighting is selected as $w_c=0.8$, $w_{ri}=0.15$ and $w_{rs}=0.05$ for both cases. The two regularization terms are activated after iteration 30 to ensure that they do not dominate the optimization at the beginning.. Since the SDF property disappears when the MBO operator is applied, the two level set functions are redistanced in each iteration. The c parameter of the MMA algorithm is increased from a starting value of 0.1 to 400.0 within 250 iterations by means of a linear ramp.

In Figure 17a,b the initial hole structure for the multi-material case of the box as well as the optimization result is shown. The structure of the multi-material result is very similar to the single-material case. In this case, the square cross strut structure has formed. The less stiff material is shifted by the optimization near the load application.

Figure 18a,b show the evolution of the objective function and the evolution of the two volume constraints.

In Figure 19a,b the initial hole structure for the multi-material case of the cantilever and the resulting optimization result are shown. The optimization result is very similar to the single material result. The material with the lower stiffness is shifted to the rear part of the cantilever in particular. It is also clearly visible that the bars have a greater thickness in the areas of the less stiff material (blue) in order to compensate for the disadvantage of the lower Young’s modulus.

Figure 20a,b show the evolution of the objective function and the evolution of the two volume boundary conditions.

Since the developed optimization method for the multi-material case reuses the already validated components from the single-material case and the optimization results for the multi-material case are also plausible, the developed method is also considered validated for the multi-material case.

4. Sensitivity Studies

In the following, the influence of the two parameters element edge length ratio of the two meshes and filter radius on the optimization result is examined. One study is carried out for each parameter. In both cases, the symmetrical cantilever clamped on one side serves as the object of investigation for the studies. The values of the other parameters are selected analogously as above.

4.1. Edge Length Ratio

In Figure 21, the optimization results are shown for different edge lengths of the FE mesh with a constant edge length of the level set mesh (change in ratio). It can be seen that the optimization runs in different local optima, but no systematic pattern can be identified. In addition, the local optima are very close to each other in terms of structure and material distribution.

Further numerical experiments have shown that satisfactory results can always be achieved with approximately the same edge length of the two meshes.

4.2. Filter Radius

In Figure 22, the optimization results for three different filter parameters of the Helmholtz filter are listed. It can be seen that this has almost no influence on the final optimization result in the selected range. If the filter radius is increased further, it appears that the optimizer generates coarser features, but the two materials are separated during the optimization and it is therefore aborted.

5. Conclusions

In this article, a multi-material topology optimization method was developed to support the product developer for the development of structures in multi-material design. The method combines the Velocity Field Level set method with the Reconciled Level set method. The method uses a constant mesh for the development of the level set functions and a second adaptive mesh for the evaluation of the system response. For the validation of the method, it was first validated on two single-material and then on two multi-material examples.

The influence of the edge length ratio of the two meshes and the influence of the filter parameter on the optimization result were examined in two sensitivity studies. It was shown that the influence is low as long as the parameters are within a reasonable range.

The method forms the basis for further optimization problems that require an explicit resolution of the structures in the multi-material design during optimization. This includes the authors’ current research, in which the influence of the interfaces between the materials, as presented for the first time in [27], is taken into account. The authors plan to publish their findings in the near future.