A boundary strip indicator for material distribution‑based topology optimization

This article proposes a boundary strip indicator for density-based topology optimization that can be used to estimate the design’s surface area (perimeter in 2D) or identify a coating layer. We investigate the theoretical properties of the proposed boundary strip indicator and propose a differentiable approximation that preserves key properties, such as non-negativity. Finally, we use the boundary strip indicator in a heat conduction design optimization problem for a coated structure. The resulting designs show a strong dependence on the properties of the coating.


Introduction
Topology optimization represents a set of methods used to solve design optimization problems (Bendsøe and Sigmund 2003).The goal of topology optimization is to generate a design that maximizes performance while adhering to specifications (Allaire et al. 1997; Bendsøe and Kikuchi 1988).Typically, the performance is measured using an objective function, which is to be extremized while accounting for specifications encoded in constraints and penalty functions.
A convenient way to get topological or geometrical information about designs is through mathematical morphology, a branch of image analysis stemming from the study of images as mathematical sets.In mathematical morphology, a binary image is examined by probing it with a smaller set, called the structuring element, at every position in the image (Heijmans 1995).
The introduction of mathematical morphology dates back to studies on spatial set algebra (Minkowski 1903) and topology (Matheron 1974).The research on this topic reached one of its applications in digital image processing (Pratt 1991), with the scope of analyzing pictures and restoring noisy ones, among others.In this article, we consider the dilation and erosion operators, first used for topology optimization by Sigmund (2007) to control the minimum feature and minimum inclusion size of the designs.
In a binary setup, a 2D design may be represented as a black-and-white image where black pixels indicate the presence of material and white pixels absence of material.The erosion operator removes a strip of width given by the operator radius from the border of the considered design.In contrast, the dilation operator has the opposite effect, adding material in a strip along the boundary of the design.As a consequence, the dilation operator fills any hole in the design that is smaller than the operator's diameter.Figure 1 shows (i) the structuring element S r , (ii) an example design , (iii) the effect of the erosion operator, and (iv) the effect of the dilation operator.
To enable gradient-based optimization and a fast adjoint-based evaluation of the associated sensitivities, we relax the binary problem by allowing the design variables to take values in the continuum [0, 1].Moreover, we employ fW-mean filters (Wadbro and Hägg 2015) as differentiable approximations of morphological operators.Linear density filters were first introduced in density-based topology optimization by Bruns and Tortorelli (2001) to combat mesh-dependence of the designs and the technique has since been refined using a range of non-linear filters (Guest et al. 2004;Sigmund 2007;Svanberg and Svärd 2013).We note that most filters used in densitybased topology optimization may be analyzed within the fW-mean filter framework (Wadbro and Hägg 2015).
During the last decade, there has been a growing interest in developing methods that enable the identification of a strip around the boundary of the material region.For level set methods, Wang and Kang (2018) proposed a double threshold in the level set function to describe the substrate and the coating.Moreover, Bai and Zuo (2022) upgraded the latter method using multiple level set functions, similarly to the "color" level sets (Wang and Wang 2004), to describe coating problems in a multi-material setting.Another example, is the work by Fu et al. (2019), which uses a levelset-based method, which can easily identify and handle the design's boundary during the design process to design shell-infill structures.However, also in density-based topology optimization there are several methods able to describe material interfaces.Clausen et al. (2015Clausen et al. ( , 2016) ) identified the boundary strip using the spatial gradient of the physical design.Luo et al. (2019) defined the boundary strip as the difference between the original and eroded structure.Yoon and Yi (2019) identified the boundary strip using a combination (similar to the xor operator) of the physical design and a non-linearly filtered counterpart.Wadbro and Niu (2019); Niu and Wadbro (2021) identified the coating using the difference between two dilate-erode operators using different radii.In this article, we propose a boundary strip indicator defined as the difference between a dilation and an erosion filter.The key point separating this boundary indicator from indicators based on spatial gradients lies in its simplicity.While, for example, the method developed by Clausen et al. (2015) involves several different steps (involving projection and smoothing), this novel indicator employs only a standard filtering step.Furthermore, the definition ensures that the boundary strip defined by our indicator has uniform thickness.The proposed indicator has two tuning parameters, namely the radius of the filters, which controls the thickness of the indicated boundary strip, and a non-linearity parameter, which controls the non-linearity of the filters.
The main result of this article is to show that the proposed boundary strip indicator possesses crucial properties of an indicator, such as non-negativity.
The fact that these properties hold sets the proposed boundary strip indicator apart from previous boundary strip indicators.Moreover, using non-negativity and known properties of fW-mean filters, we prove that, for (almost) binary designs, there is (almost) no overlap between the eroded design and the region identified by the boundary strip indicator.
Finally, we demonstrate how the boundary strip indicator can be utilized when optimizing the design of composite materials involving a substrate and a coating layer with different properties.

Test problem
Given a design optimization problem on a set Ω ⊂ ℝ d , we define a material distribution function in which = 0 indicates absence of material and = 1 pres- ence of material.
To test the boundary strip indicator, we solve a modified version of the standard minimum heat compliance problem over a square domain Ω , as depicted in Fig. 2. We impose ther- mal insulation along the boundary Γ N and set the temperature to zero along the boundary portion Γ D .Moreover, we assume to have material on the left side of the domain.The square domain Ω is exposed to uniform heating.
The set of admissible temperature fields is then For every v ∈ U , we define the linear functional Here, the boundary consists of two parts: Γ D , the portion employing Dirichlet boundary conditions, and Γ N , the portion with Neumann boundary conditions.We assume to have material (in grey) on the left side of the domain representing a uniform source of heat, and the bilinear form where P( ) is the physical thermal conductivity based on the SIMP method (Bendsøe and Sigmund 2003).More precisely, given a minimum physical thermal conductivity  > 0 and a penalization parameter p > 1 , the physical thermal conductivity is given by The minimum physical thermal conductivity ensures that the equilibrium temperature is uniquely defined by the following variational problem: find u ∈ U such that Imposing an upper bound on the volume fraction of Ω con- taining material, V max , we obtain the following formulation of the standard heat compliance problem, Note that penalized material distribution problems, such as the one above, are typically ill-posed-the solutions to the discretized problem are mesh-dependent.Over the years, several strategies to resolve the ill-posedness have been proposed.Borrvall (2001) provides a comprehensive review of classic restriction methods.Here, we are interested in solving a problem where the design includes a coating layer.We briefly describe the implemented regularization when introducing the formulation of the problem involving the composite material below.

Boundary strip indicator
Originally, the morphological operators were defined for sets.The definitions naturally generalize to the case of binary and gray-scale images.The erosion and dilation of a set M by a structuring element B are given by the Minkowski difference M ⊖ B and sum M ⊕ B , respec- tively.Below, we use the extensions of these operators for functions on a bounded domain.For more information about these operators, we refer the reader to the work by Hägg and Wadbro (2018).Let S r (x) be a structuring element having origin in x and radius r.For designs ∶ Ω → [0, 1] , the erosion and dilation operators with respect to the set Ω are defined for x ∈ Ω as and respectively.
Given a material distribution , we introduce the following definition.

Theorem 1 The boundary strip indicator is binary for binary
We consider a discretization of Ω , consisting of n ele- ments Ω 1 , Ω 2 , … , Ω n , and approximate the design variable with an element-wise constant function h .We represent the discretization of the design variable as the vector ∈ [0, 1] n with components i = h (x i ) , where x i ∈ Ω denotes the centroid of element Ω i .
To define the differentiable boundary strip indicator B ( ) on the design vector , where  > 0 is a parameter, we proceed as follows.We use harmonic filters (Svanberg and Svärd 2013), which are examples of fW-mean filters (Wadbro and Hägg 2015), as regular approximations of the erosion and dilation.More precisely, the harmonic erode is defined as where and its inverse are applied componentwise.The weight matrix W = [w ij ] ∈ ℝ n×n satisfies where |N i | denotes the cardinality of the neighborhood of x i , which is defined via the structuring element.
Denoting the all-ones vector as 1 = (1, … , 1) T ∈ [0, 1] n , we define the harmonic dilate filter as In the end, we are led to the following definition.Definition 2 For a vector ∈ [0, 1] n , the harmonic boundary strip indicator is defined as Remark 1 Analogously, we may define boundary strip indicators using other fW-mean filters, such as the geometric, exponential, and power erode and dilate filters.Svanberg and Svärd (2013) proved that D ( ) is con- vex and E ( ) is concave, hence the following proposition holds.
As a corollary to Proposition 2, we have the following theorem.
Theorem 3 The harmonic boundary strip indicator maps into [0, 1] n .In particular, it is non-negative, Proof We start by proving non-negativity.Consider any two material distributions , ∈ [0, 1] n and let B i ( ) denote the ith element of the boundary strip indicator.Since B i ( ) is convex and differentiable, we have Thus, if 12) and ( 21) imply that As shown by Wadbro and Hägg (2015), the ijth entry of the Jacobian of E is Therefore, Eq. ( 24) becomes Since E is an fW-mean filter, we have that E (c1) = c1 for any c ∈ [0, 1] (Wadbro and Hägg 2015), which together with Eq. ( 26) and Definition 2 implies that (20) Thus, inequality (23) implies that, for any ∈ [0, 1] n , Finally, as a consequence of E ( ) ≥ 0, ∀ ∈ [0, 1] n , we have which concludes the proof.◻

Remark 2
The proof of Theorem 3 relies only on the mapping, convexity, and concavity properties of the fW-mean filters used to approximate the dilate and the erode filters.Thus, everything holds also for the exponential and the geometric filters defined by Svanberg and Svärd (2013).
Noting that E , D ∶ [0, 1] n → [0, 1] n and recalling Defini- tion 2, we now derive some corollaries to Theorem 3. First, we immediately obtain the ordering for any ∈ [0, 1] n of the harmonic erosion with parame- ter  > 0 and the harmonic dilation with parameter  > 0 .Note that this is the same ordering that holds for the (exact) erosion and dilation.Second, we demonstrate that E ( ) and B ( ) are complementary in the sense that if one is large the other must be small.More precisely, for a given vector ∈ [0, 1] n and element i, we define ∈ [0, 1] as = 1 − E i ( ) .For this , we have Similarly, for a given vector ∈ [0, 1] n and element i, we define δ ∈ [0, 1] as δ = 1 − B  i ( ) , which yields the inequalities Let us now consider E ( ) and B ( ) that are (almost) binary; that is, E ( ) and B ( ) are (close to) either 0 or 1 for (almost) all elements i.Then, if E i ( ) = 1 − is (close to) 1 for some i; that is, is (close to) 0, then B i ( ) is (close to) 0 by bound (32).Similarly, if B  i ( ) = 1 − δ is (close to) 1 for some i; that is, δ is (close to) 0, then E i ( ) is (close to) 0 by bound (33).Thus, in the (almost) binary case there is (almost) no overlap between regions where E ( ) is (close to) 1 and regions where B ( ) is (close to) 1.In the follow- ing numerical experiments, it is this complementarity that guarantees that if the final design is (almost) free of intermediate values, then there is (almost) no overlap between the substrate and the coating layer.To experimentally verify this complementary, we suggest a new quality measure in Eq. ( 43).The proposed measure, which quantifies the overlap between the substrate and the coating, is reported for the final designs obtained in the numerical experiments presented below.Furthermore, the proposed indicator relies only on fWmean filters, that are independent from the chosen discretization.Additionally, we do not introduce new tuning parameters, as the proposed indicator relies only on the filters' radius and non-linearity parameter, as in other standard filtering techniques.

Boundary strip as coating layer
In this section, we show one of the possible applications of the boundary strip indicator.Following Yoon and Yi (2019), we introduce a coating layer in the design.Here, for a generic material distribution vector , E ( ) repre- sents the substrate of the coated structure, while B ( ) describes the coating layer.Accordingly, the outline of the complete design, consisting of both substrate and coating layer, corresponds to E ( ) + B ( ) = D ( ).
To achieve minimum length scale control on the final design, we employ a filtering technique inspired by the open-close strategy from Sigmund (2007).To avoid the extra cost associated with robust formulations (Sigmund 2009;Wang et al. 2011), we employ the scheme suggested by Hägg and Wadbro (2018).That is, the material distribution entering the physical conductivity is defined using a harmonic opening of the design vector while a harmonic closing is employed in the volume constraint to enforce minimum length scale control on the layout of the composite structure.This filtering technique promotes a minimum size control on the substrate, thus impeding the appearance of thin features consisting of coating only.For the sake of notational simplicity, given the design vector , we let be the design vectors after the application of the open and the close filters, respectively.
Before defining the physical thermal conductivity and the volume, we introduce two gain parameters that express the difference in properties of the coating and the substrate.The conductivity gain,  gain > 0 , is defined by the ratio of the thermal conductivity of the coating, c , and that of the substrate, s ; that is, Without loss of generality and to match the presentation in Sect.2, we let s = 1 .Similarly, the cost gain, C gain > 0 , expresses the relative cost of the coating, c c , compared to that of the substrate, c s ; that is, Using a SIMP-inspired material interpolation scheme, we express the physical thermal conductivity of the composite structure as where the erode filter identifies the substrate and the boundary strip indicator the coating layer.Introducing an upper bound on the total cost, leads us to a weighted volume constraint of the form Considering the discretization from Sect. 3, we employ the finite element method for solving the governing equation.In this case, we rewrite the bilinear form in Eq. ( 4), substituting the SIMP expression for the physical thermal conductivity of the homogeneous material with the physical thermal conductivity of the composite structure P c ( O ) , as where u is the discrete representation of the temperature field, K c ( O ) represents the stiffness matrix, and f the source term.
Finally, we use a compliance type objective function, J(u) = f T u .Thus, the discretized optimization problem we study is To solve the optimization problem, we use an optimality criteria method (Bendsøe and Sigmund 2003).To evaluate the effect of the filters on the sensitivities, we adapt the methodology of Sigmund (2009) and Hägg and Wadbro (2017).In this case, following the chain rule, the gradient of the objective function is calculated as

Numerical experiments
This section presents and discusses numerical results on optimization problem (40).The simulations aim to examine the effects on the final design when changing the thickness of the coating layer or the gain parameters gain and C gain .
The numerical experiments are performed in MATLAB R2022b.
We discretize the computational square domain of side L by 256 × 256 square elements.A continuation method similar to the one used by Tovar and Khandelwal (2011) is applied to gradually adjust the penalization and filter nonlinearity parameters p and .To promote a grey-scale-free final design, we increase the penalty parameter and decrease the non-linearity parameter until p = 3 and = 10 −8 .To evaluate whether a design has gray areas or not, we use the measure of non-discreteness inspired by the one suggested by Sigmund (2007).Measure (42) assumes the value 0% only when the complete design, defined as D ( C ) = E ( C ) + B ( C ) is binary.Note that the complete design corresponds to the (complete) layout of material entering the volume constraint (38).
Moreover, to evaluate the overlap between the substrate and the coating layer, we define the measure of overlapping fraction as The latter indicates the cardinality of the intersection between the sets of elements where the substrate and the coating layer assume values greater than a cutoff value of 1/2, divided by the number of elements.
In the figures presented below, the layout of the substrate, E ( C ) , is presented in a white-cyan color scale and the layout of the coating, B ( C ) , in a white-black color scale, where white corresponds to void in both color scales.Thus, the complete layout of the material in the figures corresponds to the layout of the complete design used to measure the non-discreteness.
To study the effects that different coating layer thicknesses and costs may have on the final design, we analyze two cases.In the first case, we vary the thickness of the ( 41) coating layer, and in the second case, we vary the gain parameters.In both cases, we keep the gain parameters equal, so gain = C gain .
To vary the thickness of the coating layer in the first experiment, we vary the filter radius.Specifically, if we consider the boundary strip indicator with a filter radius r defined as a fraction of the domain's size, the width of the coating layer will be 2r.Figure 4 shows how the final design changes when varying the filter radius.We display the results for a filter radius equal to L/128, L/64, L/32, and L/16, respectively.The gain parameters are gain = C gain = 2.
In our second experiment, we investigate how different thermal conductivities and costs of the coating affect the final design.In these numerical experiments, we keep a fixed filter radius of L/64 and solve the optimization problem for gain = C gain = 2, 4, 8, and 16. Figure 5 displays the result- ing designs.
Figures 6 and 7 show the evolution of the simulation for r = L∕64 and r = L∕32 , respectively, with gain = C gain = 2 .Here, we notice the growth of the coating along the simulation.This effect is a direct consequence of the choice of a density-based method.Accordingly, around the initial iterations, the M nd is still too high to identify a sharp boundary strip around the design.Therefore, we appreciate the effect of the coating layer on the design the more we get closer to the end of the simulation.Nonetheless, by comparing Figs. 6 and 7, we can observe that the coating influences the evolution of the designs from the first steps of the simulations.
Note that in all cases, the measure of non-discreteness of the final designs ranges between 0.0001% and 0.0016% .Conse- quently, by computing the measure of overlapping fraction, we get a value of M of = 0 for all the presented results.Fig- ures 6 and 7 also illustrate the history of the objective during the optimization process.In contrast to a standard problem, the objective value appears to grow during the iterations.This is because we are solving a sequence of problems, each with increasing penalty and non-linearity parameters, which makes the coating more pronounced.One possible reason for the growth of the objective between these problems is that it is best to focus the high-conductivity material close to the boundary Γ D .However, the coating is at the boundary of the design, which is typically located far away from Γ D .
Finally, we explore the effects of a coating layer that is cheaper and less conductive than the substrate.Figure 8 shows the resulting designs using a filter radius equal to L/64 and gain parameters gain = C gain = 1∕2, and 1/8, respec- tively.However, to achieve a low M nd in the final designs for  gain , C gain < 1 we need a higher final penalty parameter, namely p = 10.
We find that increasing the coating's cost and conductivity or the coating layer thickness yields smoother optimized designs.Following the intuition, the design's perimeter becomes shorter as the coating gets more conductive and expensive.Similarly, a thicker layer is more conductive and more expensive than a thin one, leading to smoother structures requiring less coating material.Recall that the lack of coating on the substrate located at the left boundary of the design domain is a consequence of our modeling choice of placing fixed material to the left of this boundary, as seen in Figs. 2 and 3. Other cases can be implemented by defining different material distributions outside the design domain.

Discussion and outlook
In this paper, we have introduced a novel boundary strip indicator and proved it is well-defined.The main difference between this indicator and other ones lies in its simplicity.While other methods might require several regularization steps (like projections and smoothing processes), this novel indicator utilizes only a standard filtering step.Furthermore, its definition ensures that the boundary strip defined by our indicator has uniform thickness.Moreover, we showed the applicability of the introduced indicator for optimizing designs with a coating layer.We have conducted numerical experiments, investigating the possibility of altering the properties and cost of the coating.
Although the numerical experiments are performed in a two-dimensional setup, the boundary strip indicator is defined and examined for a general dimension d.Therefore, it is possible to extend the treatment to a three-dimensional setup.In such a case, the indicator identifies a tubular neighborhood of the design's surface.
Moreover, we limited our attention to problems using a structured mesh, which is standard in topology optimization.The theory behind the approach is applicable to unstructured meshes as well.However, implementing the approach efficiently on unstructured meshes may require additional care and considerations.
Another application of the boundary strip indicator would be the introduction of a perimeter or surface constraint in problem (7).Further development of the presented method would make it possible to solve material distribution problems for composite materials with more than two layers.For the implementation of multi-layered materials, a new filter is needed to identify inner layers.However, the latter could be implemented as a combination of existing filters with the introduced boundary strip indicator.

Fig. 1 Fig. 2
Fig.1Effects of the erosion operator (iii) and of the dilation operator (iv) in a binary case, applied to the design in (ii).Here, the radius r of the structuring element S r in (i) indicates the radius of the operators

Figure 3 ,
Figure3, shows (i) a structuring element S r , (ii) a binary test-case design , and (iii) its boundary strip B( ) com- puted with the boundary strip indicator, assuming there is material (in black) outside of the design domain on the left side.Note that the diameter (2r) of the structuring element characterizes the uniform width of the boundary strip.The boundary strip indicator can, for example, give us a surface area (perimeter in 2D) estimate of the object defined by the material distribution function .

Fig. 3
Fig.3Binary test-case design (ii) and its respective boundary strip indicator B( ) (iii), assuming there is material outside of the design domain on the left side (in black).Here, the radius r of the structuring element S r in (i) indicates the radius of the operators

Fig. 4
Fig. 4 Here, we vary the filter radius, keeping a fixed gain = C gain = 2 .Cyan represents the substrate and black represents the coating layer

Fig. 5
Fig. 5 Final design having substrate in cyan and coating layer in black.Here, we variate the gain parameters, namely gain = C gain = 2, 4, 8, 16 , keeping a fixed filter radius of L/64

Fig. 6 Fig. 8
Fig. 6 Evolution of the simulation fixing the parameters to r = L∕64 and gain = C gain = 2 and respective objective function value