On equal-width length-scale control in topology optimization

This paper deals with equal-width length-scale control in topology optimization. To realize this aim, we first review different notions of minimum and maximum length-scale control and highlight some perhaps counterintuitive consequences of the various definitions. Here, we implement equal-width control within the moving morphable components (MMC) framework by imposing the same upper and lower bounds on the width of the components. To avoid partially overlapping beams and nearly parallel beams, as well as beams crossing at small angles, we introduce penalty functions of the angle difference and the minimum distance between any two beams. A penalized optimization formulation of compliance minimization is established and studied in several numerical examples with different load cases and boundary conditions. The numerical results show that equal-width length-scale control can be obtained by using the proposed penalty function in combination with a continuation approach for the amount of penalization.


Introduction
Design optimization aims at finding the optimal placement of materials within a given design domain ⊂ R N . In this paper, we limit the discussion to the case of two material types, for simplicity henceforth referred to as material and void, which occupy the regions M ⊂ and V ⊂ , respectively. Design optimization aims to find the "best" material region M * ⊂ , or more precisely the region that extremizes a given performance measure subject to given constraints.
The task of obtaining size control in some sense has attracted increasing interest over the last two decades. Early attempts to obtain size control in the context of material distribution topology optimization focused on limiting the spatial variation of the density field, for example, by using linear density filtering (Bruns and Tortorelli 2001;Bourdin 2001) or by adding slope constraints (Petersson and Sigmund 1998). Linear density filtering has the disadvantage of not being able to provide sharp transitions between material and void regions. During the last 15 years, there has been an intense interest in developing new non-linear filtering strategies (Guest et al. 2004;Sigmund 2007;Xu et al. 2010;Wang et al. 2011;Svanberg and Svärd 2013) that have the advantage of allowing (almost) sharp transitions between material and void. Recently, Wadbro and Hägg (2015) presented a general filtering framework that includes the vast majority of available filters, which enables a unified analysis of a large class of material distribution methods. Lazarov et al. (2016) comprehensively review different techniques based on density filtering for imposing lengthscale and manufacturability in topology optimization. The problem of controlling the minimum and maximum lengthscales of the design in topology optimization has attracted much interest. In some cases, optimized designs, obtained by using filters similar to the open-close or close-open filters introduced by Sigmund (2007), possess a minimum length-scale on both phases (Hägg and Wadbro 2017;Sigmund 2007). However, the existence of a minimum lengthscale on both phases cannot be guaranteed (Schevenels and Sigmund 2016). The maximum and minimum length-scales have been controlled by Guest (2009), who introduced local constraints for maximum length-scale control together with some Heaviside based filtering for minimum length-scale control. Zhang et al. (2014) introduced constraints, essentially based on length-scale definition (4) below, to the problem formulation and were able to control both the minimum and maximum length-scale of one phase.
2 On length-scale of sets in R N Throughout this paper, we assume that the material region M is a nonempty N-dimensional open set, so M ⊂ R N . Here, we use another set B ⊂ R N , called the structuring element, to define and obtain information about the minimum and maximum length-scale of M, respectively. For r ∈ R + and y ∈ R N , we let and That is, rB denotes the set obtained by scaling B by r and y + rB denotes the set obtained by translating rB by y.
Here, we restrict our attention to the case where the structuring element contains 0 and is open, bounded, symmetric, and convex. Hence, B is the open unit ball for some metric on R N ; more precisely, the Minkowski functional p B (x) = inf{r > 0 | x ∈ rB} is a norm on R N .
Following Hägg and Wadbro (2018), we start by defining the local length-scale relative B at x ∈ M,R B (M; x), as the radius of the largest ball (scaled translate of B) that contains x and is contained within M. That is, By using this local length-scale, the minimum length-scale of M relative B is defined as Similarly, the maximum length-scale of M relative B is Note that, definition (3) specifies a length-scale for each point x ∈ M, while definitions (4) and (5) establish a notion of minimum and maximum length-scale, respectively, for domain M. In words, definitions (4) and (5)

Relation to morphological operators
In mathematical morphology, the aim is to gain information about M by probing it by B. The basic morphological operators are erosion, dilation, closing, and opening. Following Heijmans (1995), we define the dilation and erosion of M by B as and respectively. Figure 1   The minimum and maximum length-scales R B and L B can also be defined by using the morphological operators and the convention sup ∅ = 0 as and respectively; recall that rB is the unit ball scaled by r.
A proof of the equivalence between minimum length-scale definitions (4) and (10) is given by Hägg and Wadbro (2018, Theorem 3). The equivalence of definitions (5) and (11) can be shown similarly.
The equivalence between definitions (4) and (10) as well as definitions (5) and (11) provides us with alternative interpretations of length-scales R B and L B . Minimum length-scale R B is the largest radius r so that M can be completely pained with a brush with shape rB without touching any point outside M. Maximum length-scale L B is the smallest radius r so that M is completely removed by an eraser with shape rB keeping the eraser's center outside M.

Relation to the medial axis transform
The medial axis transform, which is extensively used to describe shapes in the field of computer graphics and pattern recognition (Loncaric 1998). Blum (1967) introduced the concept of describing the shape of an object occupying a domain M by its medial axis or skeleton S and a distance map d : S → R + . The skeleton consists of all points in M for which the closest point in R N \ M is not unique, that is Metrics defined via distance functions have been used extensively as a practical tool in the field of computer visualization, while these metrics primarily have been used as a theoretical tool for shape calculus, see for example Delfour & Zolesio's book (Delfour and Zolésio 2011). Note that given M and its skeleton S, there may be multiple distance maps d that satisfy condition (13). To strictly define minimum and maximum length-scales based on the medial axis transform, we need to add a requirement on the distance map d. More precisely, we will, in addition to require that d satisfies condition (13), also stipulate that d(x) is maximal for each x ∈ S. Provided that d is selected as stipulated above, minimum and maximum length-scales can also be defined by using the distance map d associated with the medial axis transform, as suggested in the context of topology optimization by Zhang et al. (2014). The minimum length-scale R B (M) can also be defined as and the maximum length-scale L B (M) can also be defined as That is, the minimum and maximum length scales of domain M may be defined as the smallest and largest value of the distance map for all points on S, respectively. Allaire et al. (2016) proposed an alternative minimum length-scale definition based on the thickness of the design. Consider a point x on the boundary ∂M and assume that its outward directed unit normal n(x) at x is well-defined. Then, the local thicknessT (M; x) of M at x is defined aŝ

Alternative minimum length-scale definition based on thickness
In other words,T (M; x) is the largest t such that for any s ∈ (0, t), it holds that x − sn(x) ∈ M. The minimum thickness T (M) of M is then defined as Later in this manuscript we will also refer to T (M) as the minimum width of the domain M.
We remark that, as also Allaire et al. (2016) acknowledged, the behavior of the minimum thickness T (M) and the minimum length-scale R B (M) can be rather different. The minimum thickness of a square with side length 1 and rounded corners is 1, while R B for this domain is the rounding radius of the corners provided that B is the unit sphere.

The equal length-scale paradox
The end-goal of this work is to find designs with equal member width, or in other words, we want to find designs whose minimum and maximum length-scale coincide. However, this requirement together with length-scale definitions (4) and (5) implies that we cannot have designs that contain equal-width bars that cross. Figure 2 illustrates two bars of width 2 and infinite length that cross. At the crossing point, the local length-scale is √ 2; away from the crossing region, the local length-scale is one. More precisely, in this case, R B = 1 < √ 2 = L B . To avoid this paradoxical consequence of definitions (4) and (5), we choose to use an alternative definition of the maximum length-scale.

Maximum length-scale revisited
Instead of defining the maximum length-scale of M by using the local length-scales or directly by using the distance map and skeleton associated with the medial axis transform, one could try to mix concepts from the medial axis transform and mathematical morphology. First, note that the structural skeleton S may include more points than C, the set of all center points of the balls defining the local length-scale. Formally C may be defined as For the domain in Fig. 2 where two crossing beams comprise M, the left and right images in Fig. 3 illustrate the sets C and S, respectively. The skeleton S includes all points in M that do not have a unique closest point on the boundary of M, while the set of all center points C holds only the center points of the circles that define the local length-scalê We define the maximum width D B (M) as twice the smallest scaling factor r such that M is a subset of its skeleton dilated by rB. That is, The factor 2 is included to account for that the thickness of a long bar is twice the radius of the largest ball that fits inside the bar. We remark that by using the maximum length-scale constraint proposed by Zhang et al. (2014) (More precisely, when we let parameter 2 → 0 in problem (3.1) in their article.), we get length-scale control according to maximum length-scale definition (19), but not necessarily according to maximum length-scale defintion (5). (Recall that distance map d is not uniquely defined by expression (13).) In the context when two bars cross at a right angle, recall Fig. 2, definition (19) results in a maximum width of 2, which is the same as the minimal thickness according to definition (17). However, in other cases, such as the one detailed below, maximum width definition (19) produces unintuitive results. Consider a sequence of domains M n , where each domain is the union of the unit square and 4 · 2 n balls of radius r n = 2 −n−1 centered at (0, a k ), (1, a k ), (a k , 0), and (a k , 1) for k = 1, 2, . . . , 2 n and where a k = (2k − 1)r n . The domains M 2 , M 3 , and M 4 are illustrated in Fig. 4. For domain M n , the minimum length-scale following definition (4) is R B (M n ) = 2 −n−1 and the maximum length-scale following definition (5) is L B (M n ) = 0.5. (Recall that the side length of the unit square is 1, while the  (18) and (12), respectively. Here, M is the same as in Fig. 2; its boundary is illustrated by the thin dashed lines Fig. 4 A sequence of domains obtained as the union of the unit square and balls of radius 1/4 (left) 1/8 (middle) and 1/16 (right) centered along the boundary of the unit square radius of the unit ball is 1.) However, the maximum lengthscale according to definition (19) is D(M n ) = 2 −n−1 . Note that, the skeleton (illustrated by the black lines in Fig. 4) becomes increasingly dense in each step. In essence, the closest point from a point inside the unit square to the boundary of M n is one of the points where the boundary of the unit square and the boundary of M n intersect. In the limit n → ∞, we have that R B (M n

Discussion
Above, we have discussed various possibilities of quantifying the minimum and maximum length-scales of a domain M. To proceed towards our end-goal of imposing equalwidth length-scale, we need to select which definitions to use. In this study, we use definition (17) for the minimum length-scale and definition (19) for the maximum length-scale. By using these definitions, we avoid the equalwidth paradox discussed in Section 2.4. Therefore, we can use these definitions for the equal-width length-scale control without inconsistency. However, we may still arrive at the case illustrated in Fig. 4, where the maximum length behaves unintuitively. The problem is that, due to the wiggly boundary of the domain, the skeleton includes many lines that are located close together.

Problem formulation and design specification
Here, we study a standard test-problem, compliance minimization of a cantilever beam. We define our design domain D to be a rectangle with height H and length L as illustrated in Fig. 5. The structure is clamped at its left boundary D and subject to a traction force f over the right boundary L . Note that f does neither need to be constant nor nonzero over the full boundary L ; for each experiment below, we will specify the load condition. The material is isotropic with Young's modulus 1 and Poisson's ratio 0.3. Below, we assume that material occupies the region ⊂ D (and that supp{f } ⊂ ∂ , which is Lipschitz). We let u denote the steady state displacement of the domain. The task is to minimize the compliance of the structure, that is to minimize the integral of f · u over L .
In the present study, we use the so-called Moving Morphable Components (MMC) approach by Guo et al. (2014). This method provides an explicit and geometric description of the topological designs. In principle, the domain occupied by material is composed of N c components, that is, = ∪ N c n=1 n , where n denotes the region occupied by the nth component. In MMC, each component is described by its center position (x c , y c ), length L > 0, width (or thickness) t > 0, and orientation (directional sine) η ∈ [−1, 1]. That is, the domain n occupied by the nth component is parameterized by vector d n = [x c n , y c n , L n , t n , η n ] T . Thus, to describe the layout of a structure comprising N c components, we use a design vector The minimum compliance of the structure with the constraint on the volume of material can be formulated as:

Fig. 5
The material region is a subset of the rectangular shaped design domain D . The structure is fixed along the left wall D and a traction force is applied on the boundary segment L where U d is the set of feasible designs and C(d) is evaluated as (u) for the displacements u that solve equilibrium equation, detailed below, for the domain defined by the design d.
Instead of describing accurately the boundary of the components using X-FEM approach or adaptive meshing along boundary as proposed by Guo et al. (2014), we use a fictitious domain approach (also called ersatz material model) by adopting the MMC implementation by Zhang et al. (2016). That is, we use a fixed and regular finite element mesh, and the element stiffness is interpolated based on the values of the topological description functions (TDF) of the nodes of the corresponding element. There is one topological description function for each component, e.g., the region n occupied by the nth component can be described by the TDF φ n , defined so that where x is an arbitrary point in the design domain D and ∂ n is the boundary of n .
The structure is fixed at D ; thus, the set of admissible kinematic displacements are The state problem is to find the displacements u ∈ U such that Here, the bilinear for a and the linear form are defined as where E is a constant fourth-order elasticity tensor corresponding to the given isotropic material, (u) is the second order strain tensor of u, the colon ":" denotes the contraction of two tensors.
To realize equal-width control with a pre-specified width t eq , the upper and lower bounds of the design variables corresponding to component width must all be the same; that is, t n = t eq for n = 1, . . . , N c . However, this condition is not sufficient to guarantee that the full structure has equal structural width overall. Fig. 6 illustrates two such cases. Hence, to avoid partial overlapping beams as well as beams crossing at very small angles as well as nearly parallel beams, we choose to penalize such configurations as detailed in the next section.

Penalty formulation
As discussed above, simply ensuring that the width of all components is the same is not sufficient to ensure equal-width length-scale control. To obtain size control, we need to ensure that unwanted configurations, such as those illustrated in Figs. 4 and 6, are suppressed.
In the discussion in this section, we first consider two beams, both of width t eq , parameterized by d i and d j as described above, respectively. That two beams partially overlap is by itself not a problem, since the crossing angle may be large. Similarly, the presence of two parallel beams is by itself unproblematic, as they may be non-overlapping. However, we would like to avoid any situation where two beams are nearly parallel and partially overlap, except if the overlap occurs close to both beams'. Thus, we penalize the product of a function depending on the angle difference between the two beams and another function that depends on the minimum distance between the two beams. Below, we detail how we have selected these functions for the current application.

Penalty for angle difference
We let θ i = arcsin η i and θ j = arcsin η j ; recall that η i and η j are the variables that determine the direction of the ith and j th component of , respectively. First, let us start by considering the penalty with respect to angle difference. First, we want this penalty to be symmetric in that the numbering of the two beams should not matter, that is, only the absolute angle difference should matter. To avoid the non-differentiability of the absolute values, we use a penalty function that is based on the squared angle difference. More precisely, we let the penalty function depend on the approximate angle difference where ε is a small positive regularization parameter. The function α is differentiable in θ i and θ j and approximates |θ i −θ j | with a maximum error dictated by the regularization parameter ε. Moreover, we want the penalty function to have its maximum value when the bars are parallel, that is, when α(θ i , θ j ) = 0 and its minimum for α(θ i , θ j ) = π/2 and we require that this function be symmetric around π/2. One example of such a function, which we use in our numerical experiments, is There are many choices of possible penalty function that satisfy the above conditions as well as smooth approximations of the absolute value function. The functions used here are merely used to serve as a proof of concept for the proposed approach to obtain equal-width length-scale control. Figure 7 graphs the angle penalty p a and its exponentiation (p a ) q with power q = 6 as functions of the absolute angle difference |θ i − θ j |. Here, as well as in the numerical experiments, we use = 10 −4 to compute α(θ i , θ j ) following definition (25). Figure 7 shows the p a and p 6 a as functions of the absolute angle difference |θ i − θ j | for two crossing bars. The angle penalty is scaled so that it is 1 when the bars Fig. 7 Penalty function p a versus absolute angle difference |θ i − θ j | for parameter = 10 −4 are parallel, or fully overlapping, while the angle penalty is negligible when the crossing angle is close to 90 degrees. For the basic angle penalty, p a the amount of penalization is small only for a small range of angles around 90 degrees, and as a consequence, many of the bars that cross in the optimized designs will do so with a crossing angle around 90 degrees. By using the exponentiation (p a ) q , for example, with power q = 6, as shown in Fig. 7, the angle penalty will have relatively small value in a wider range around 90 degrees. Thus, by using p 6 a instead of p a , bars may cross in a large range of angles with only a minor influence on the value of the angle penalty. In a practical situation, it is likely that one wishes to have more control, for example, be able to select at which angle the penalization "starts". To obtain such control, one can use some projection-based technique, such techniques have been successfully used in density-based topology optimization (Xu et al. 2010;Wang et al. 2011).

Penalty for inter-component distance
We start by noticing that each point on the centerline of beam i can be written on the form p i (s i ) := (x c i , y c i ) + s i (cos θ i , η i ), where s i is a parameter satisfying |s i | ≤ l i /2. Here, we are interested in computing the distance between the two beams. To do so, we study the minimum distance between a point on the centerline of beam i and a point on beam j 's centerline. In principle, if this distance is larger than t eq , the beams do not overlap. However, the beams need not overlap if the distance is smaller, for example, if the endpoints of the two-beams are the points of minimum distance. However, if we only consider points p i (s i ) := (x c i , y c i ) + s i (cos θ i , η i ) with |s i | ≤ q i = (l i − w sp )/2, then the two beams will (at least partially) overlap if the minimum distance between their centerlines' mid-sections is smaller than t eq . The symbol w sp is defined as the predetermined space along the centerlines at the two ends of the beam. In cases where the length of the component is very short, that is if l i is smaller than w sp , then we will consider all points along the centerline by setting q i = l i /2. Thus, we compute Next, we expand the objective function above to find that So, evaluating the minimum distance between the two segments boils down to solving a quadratic programming problem. If the lines are not parallel, this problem has a unique solution that can be computed very quickly, for example, by using an active set algorithm. When evaluating derivatives, we can capitalize on that the same constraints will be active in the perturbed case as in the original when considering sufficiently small perturbations. (In some cases, we may need to settle for a one-sided difference.) Note that, if no constraints are active, then the segments cross and the distance between them is zero-after any small perturbation, the segments will still cross and thus the gradient will be zero.
On the other hand, if the lines are parallel, the quadratic problem may no longer have a unique solution. The case when the problem does not have a unique solution is the worst possible case; in this case, the two centerline segments are parallel and there are multiple pairs of closest points between these segments. In this case, one could argue that the gradient with respect to the length as well as rotation of the segments can be considered to be 0.
We note that if the distance between centerlines is greater than the imposed component thickness, then the segments do not overlap. Here, to also penalize configurations with nearly overlapping components, we use the following tanh inspired penalty function: where k is a slope parameter. Figure 8 shows penalty function p d as a function of inter-component centerline distance β for three values of parameter k. In the numerical experiments, we use k = 5.

Summing up the penalties
As mentioned already at the beginning of this section, we would like to penalize structures that include beams that are nearly parallel and partially overlap. To do so, we form a total penalty function P by summing the product of the angle penalty and distance penalty for all pairs of beams in the structure; that is, We add a scaled version of this total penalty function to the objective function and obtain our penalized optimization problem: where γ is a penalty parameter.

Center loading
Here, the design domain D , illustrated in Fig. 5, has length L = 1.5, height H = 1, and is meshed into 300 × 200 elements. A concentrated traction force is applied at the center of L and at most 40 % of the design domain is allowed to be filled with material. We first solve the problem without considering our endgoal of obtaining a structure with equal-width control. More precisely, we solve problem (31) with parameter γ = 0 to find a design with minimum compliance by the basic MMC method with the lower and upper thickness bounds set to 0.02 and 0.1, respectively. For all runs with central loading, a symmetry condition on the design along the midline in the horizontal direction is enforced throughout the optimization.
The initial layout, illustrated in Fig. 9, consists of 24 components: 4 × 3 crosses uniformly spaced in the design domain. Figure 10 shows the optimized design from this experiment. The compliance of the optimized beam is 40.00. As a second reference study, we set the lower and upper bounds for the thickness to t eq = 0.08, but we do not include any penalty to promote design with components of equalwidth, so γ = 0 also in this experiment. The compliance of the optimized beam shown in Fig. 11 is 40.90. These reference designs in Figs. 10 and 11 do not have the desired equal-width length-scale control. Each component in the design in Fig. 11 is restricted to have a specified thickness. However, the complete designs include multiple components that are nearly parallel and partially overlap; so, the full structure has a maximum width that is much larger than t eq . For parameter q = 1, the values of total penalty function (30) for the designs in Fig. 10 and Fig. 11  Having established our reference solutions, we turn to the problem of obtaining a design with equal-width length-scale control. To that aim, we solve optimization problem (31) by MMC, impose lower and upper bounds for the thickness by setting t eq = 0.08, and use a so-called continuation strategy to increase penalty parameter γ and gradually increase γ from 0 to 60 during the optimization. The space w sp is set to 0.01; and the distance w used to define distance penalty p d is set to 0.2. Figures 12 and 13 show the optimized designs obtained when using total penalty (30) with parameter q = 1 and q = 6, respectively. The compliance corresponding to the optimized designs in Figs. 12 and 13 is 81.16 and 55.22, respectively.
From Figs. 12 and 13, it can be seen that all bars have the uniform width and no obvious partially overlapped or paralleled bars exist, which is also reflected by the relatively small value of the total penalty. More precisely, the total Fig. 10 Cantilever beam optimized without equal-width length-scale control using the layout from Fig. 9 as initial design Fig. 11 Cantilever beam optimized with an imposed width of all individual components, but without equal-width length-scale control. The layout from Fig. 9 was used as initial design for the optimization penalty P is 0.22 and 0.02 for the results in Figs. 12 and 13, respectively. Note that the total penalty values above are computed with q = 1 and q = 6, respectively, so they are not directly comparable to each other. However, these numbers can be directly compared to and are magnitudes smaller than the corresponding value for the optimized designs (Figs. 10 and 11) for the unpenalized problem.
Next, we investigate the influence of the initial guess for the MMC on the final result. More precisely, we study the same problem as above, with the load at the right-center point; but instead of starting from the design in Fig. 9, we start from the design in Fig. 14. This design consists of 20 components: 2 × 5 crosses uniformly spaced in the design domain. Also, in this case, we first solve the problem by MMC with the lower and upper thickness bounds of the individual components set to 0.02 and 0.1, respectively. The optimized design without equal-width control is shown in Fig. 15 and its compliance is 43.51. By just changing Fig. 12 The optimized design with the equal length control t eq = 0.08-optimized compliance 81.16 and total penalty 0.22, where the total penalty function is applied and computed using the angle penalty p a Fig. 13 The optimized design with the equal length control t eq = 0.08-optimized compliance 55.22 and total penalty 0.02, where the total penalty function is applied and computed using the angle penalty p 6 a the lower and upper thickness bounds of the individual components to t eq = 0.08, we obtain an optimized design that is similar in structure to the design presented in Fig. 15 and has compliance 44.81. Figures 16 and 17 show the optimized designs obtained when using the layout in Fig. 14 as an initial guess; total penalty (30) with parameter q = 1 and q = 6, respectively; w sp = 0.01; w = 0.2; and a continuation approach gradually increasing penalty parameter γ from 0 to 60. The compliances corresponding to the optimized designs in

Corner loading
As our second experiment, we aim for equal-width lengthscale control while minimizing the compliance of a cantilever beam loaded by a concentrated traction force applied at the lower right corner of the computational   Fig. 18 and consists of 12 components: 2 × 3 partially overlapping crosses.
Also, in this case, we present reference results obtained without equal-width control. Just as in the previous experiments, we first solve the problem of finding a design with minimum compliance with the lower and upper thickness bounds set to 0.04 and with penalty parameter γ = 0 throughout the optimization. Figure 19 shows final design, whose compliance is 89.02. Similarly as in the experiments for the center loading, simply setting the lower and upper thickness bounds to t eq only yields minor changes in the overall layout for the optimized beam.
Next, we consider the penalized problem. Figures 20  and 21 show the optimized designs obtained when using the layout in Fig. 18 as initial guess; total penalty (30) with Fig. 16 The optimized design with the equal length control t eq = 0.08-optimized compliance 84.95 and total penalty 2.29, where the total penalty function is applied and computed using the angle penalty p a Fig. 17 The optimized design with the equal length control t eq = 0.08-optimized compliance 52.62 and total penalty 0.02, where the total penalty function is applied and computed using the angle penalty p 6 a .
parameter q = 1 and q = 6, respectively; w sp = 0.01; w = 0.1; and gradually increasing penalty parameter γ from 0 to 20. The compliance corresponding to these optimized designs is 112.76 and 106.74, respectively. Next, we try to optimize the beam with corner loading and setting the upper and lower thickness bounds to 0.02. That is, the bars have half the width compared to the earlier experiments. We mesh the design domain D of length L = 1.0 height H = 0.5 by 300 × 150 elements and set the maximum volume V = 0.3| D |. The initial design for this test case is shown in Fig. 22 and consists of 16 components: 2 × 4 partially overlapping crosses. We solve problem (31) using total penalty (30) with parameters q = 1, w sp = 0.01, w = 0.1, and gradually increasing penalty parameter γ from 0 to 20. Figure 23 shows the final design, which has compliance 154.91 and total penalty 2.69.

Other boundary conditions
In our final experiment, we use another set of boundary conditions. Here, we consider the setup illustrated in Fig. 18 The initial layout of the components for the experiments with corner loading. We use the layout in Fig. 25 as initial guess. Figures 26  and 27 show optimized beams obtained when using total penalty (30) with parameter q = 1 and q = 6, respectively; w sp = 0.01; w = 0.1; and gradually increasing penalty parameter γ from 0 to 3. Here, the initial design consists of 15 bars, one horizontal bar placed so that the distributed force is applied over its top side; the remaining bars are placed in a pattern with four plus three crosses on two rows. The compliance corresponding to these optimized designs is 9.996 and 11.36, respectively. Figure 28 shows the iteration history for the optimization process that resulted in the design in Fig. 27. The scaling factor for the penalty is chosen so that both functions are of the same order. Missing data for the compliance graph corresponds to iterations where the compliance is greater than 25.

Fig. 20
Optimized design for the corner loading with the equal length control t eq = 0.04-optimized compliance 112.76 and total penalty 1.09, where the total penalty function is applied and computed using the angle penalty p a Fig. 21 Optimized design for the corner loading with the equal length control t eq = 0.04-optimized compliance 106.74 and total penalty 0.06, where the total penalty function is applied and computed using the angle penalty p 6 a Fig. 22 The initial layout of the components for the experiments with corner loading and thinner components   Optimized design for the distributed top loaded case equal length control t eq = 0.04-optimized compliance 9.996 and total penalty 1.07, where the total penalty function is applied and computed using the angle penalty p a Fig. 27 Optimized design for the distributed top loaded case equal length control t eq = 0.04-optimized compliance 11.36 and total penalty 0.08, where the total penalty function is applied and computed using the angle penalty p 6 a 6 Concluding discussion Without any additional treatment, constraining the lower and upper bounds of the components' widths in the Moving Morphable Components (MMC) approach is not sufficient to get a resulting design that satisfies the desired lengthscale constraints. This is due to that the components comprising the design may cross at very small angles as well as overlap. In this study, we added a penalty function that considers the angle difference and the distance between components to the optimization in order to promote designs with desired properties. By using a penalized optimization formulation, equal-width length-scale control has been successfully realized in several numerical examples with different loading cases.
To obtain good results with equal-width length-scale control, penalty parameter γ should be increased gradually. For example, for the results in Fig. 27, γ is zero until iteration 50. At iteration 50, γ is set to 0.025 and is then increased by 0.0005 per iteration during iterations 51-149. At iteration 150, γ is set to 0.15 and is then increased by 0.001 per iteration during iterations 151-299. During iterations 300-499, γ = 2; and after iteration 500, γ = 3. Our numerical experiments suggest that the compliance should be Fig. 28 Iteration history illustrating how the compliance C(d), the total penalty P (d), as well as penalty parameter γ change during the optimization process that resulted in the design in Fig. 27 significantly larger than γ times the total penalty function at the beginning of the optimization. This ensures that the compliance drives the evolution of the design at the beginning of the optimization process. For the experiments in this paper, it is possible to keep the magnitude of the compliance larger than γ times the total penalty function throughout the optimization. Increasing the penalty parameter too fast, so that it dominates the objective function, may drive the design into an, from compliance minimization viewpoint, undesirable state. On the other hand, if the penalty contribution to the objective function is too small, then nearly parallel components will not be properly separated.
The effects of the imposed equal-width, the minimum distance between components, and the initial guesses on the optimization results are compared and discussed. Moreover, we tested two different powers, q = 1 and q = 6, for the angle penalty (p a ) q . For q = 1, the penalty function drives bars to cross at an angle close to 90 degrees. Promoting bars to cross at a right angle may conflict with the aim of minimizing the compliance; this may also result in slow convergence and an "oscillatory" behavior of the intermediate designs. For q = 6, the penalty function is relatively small over a wider angle-interval; and the resulting designs have bars crossing at a larger variety of angles. Overall, in our numerical experiments, the choice q = 6 typically yielded better designs with a smaller value of the objective function as compared to the results obtained with q = 1. The settings for penalty parameter γ in problem (31) may have a large influence on the progress of the optimization and hence the optimized design. Here, we used a continuation approach, to gradually increase the amount of penalty applied during the optimization.
We remark that the figures showing the optimized designs show the layout of the components exactly as the MMC parameterization specifies. Some of these design includes imperfect joints, and in particular components that are nearly or only weakly connected; see, for example, Figs. 13 and 23. However, components that do not overlap may partially overlap one or more elements and thus be considered as connected from the finite element viewpoint as the ersatz material model only considers the relative area included for each element in the mesh. The incomplete connection among components could be improved by introducing a joint connection constraint for moving morphable components; see, for example, the work by Hoang and Jang (2017) who introduced masks to realize versatile thickness control and prevent overlap of components. One possible disadvantage of the penalized formulation is that the objective function may be dominated by the penalty term, which might result in that the optimized design deviates somewhat from the result of compliance minimization, for example, some hanging parts are present in Figs. 20 and 21. Another drawback of the proposed method is that the number of terms included in penalty formulation (30) grows with the number of component pairs, that is as O(N 2 c ). Hence, the problem of obtaining equal-width length-scale control with a very large number of MMCs and very thin size is still very challenging due to large analysis cost and complexity of optimization.
Equal-width length-scale control is realized in the MMC framework, but a similar implementation may be realized in some other optimization parameterization. The designs obtained by MMC combined with our newly proposed penalty depend on the initial layout of components, which is expected since the designs obtained by the basic MMC also depend on the initial component layout. Moreover, as discussed earlier, the final design also depends on how the penalty parameter is adjusted during the optimization. Some possible extensions include to consider the design with multi-material, such as to use the ordered multimaterial interpolation (Zuo and Saitou 2017) and consider the strength at the intersection.