A continuous model for connectivity constraints in topology optimization

The aim of this work is to present a continuos mathematical model that characterizes and enforces connectivity in a topology optimization problem. That goal is accomplished by constraining the second eigenvalue of an auxiliary eigenproblem, solved together with the governing state law in each step of the iterative process. Our density-based approach is illustrated with 2d and 3d numerical examples in the context of structural design.


Introduction
Connectivity is an important issue in topology optimization (TO) that affects in a different way depending on the physical meaning of the phases involved in the problem. In structural design, for instance, while optimized designs are connected, it is quite often that several enclosed holes appear in minimum compliance structures (Sigmund 2002;Challis et al. 2012;Sigmund et al. 2016;. Those holes are desirable from a stiffness perspective, but lead to a distribution of the void phase which is not connected, and that definitely complicates the fabrication of the final structure.
Something similar happens at the nanoscale when tailoring dispersion properties in photonic crystal waveguides (PhCWs). Sometimes solid regions appear as free-floating members in the final layouts that are impossible to realize (Wang et al. 2011;Christiansen and Sigmund 2021). Additionally, enclosed air-voids appear in the design of 3d bandgap structures. Those pores are advantageous from optimization point of view, but not easy to be manufactured (Men et al. 2014;Swartz et al. 2021Swartz et al. , 2022. Another scenario where lack of connectivity becomes relevant takes place when performing electrode design in piezo modal transducers. Typical electrode patterns correspond to polarization profiles which take on two values only, i.e. either positive or negative polarity, which are the phases here. In general, electrodes hardly ever present profiles of both phases connected, and they typically exhibit isolated features of like-polarity (Ruiz et al. 2018;Donoso and Bellido 2018), which makes it difficult the wiring schemes.
As far as the authors' knowledge, the first strategy to tackle connectivity in structural design was the virtual temperature method (VTM) (Liu et al. 2015;. The idea consists in solving an auxiliary linear thermal problem where void is treated as a conductive material with an internal heat source, whereas solid is treated as an insulator material. Also, some parts of the boundary behave as heat sinks in the form of null Dirichlet boundary conditions. There, void connectivity is imposed by constraining the maximum temperature in the void phase. In a similar way, a thermal problem has been used to enforce connectivity for material design in Andreasen et al. (2014).
VTM was also extended to the molding constraint (Li et al. 2018), and to the electrode connectivity problem in Donoso and Guest (2019) by solving two auxiliary thermal problems, one for each phase. And in Luo et al. (2020), VTM was modified by introducing a nonlinear heat source. That improvement lets generate a uniformly distributed temperature field in all regions covered by enclosed voids, so that the temperature threshold value is independent from the problem. In the last years, other appealing strategies, most of them supported by physical arguments, have been proposed to address the connectivity issue in structural design (Luo et al. 2020;Zhou and Zhang 2019;Wang et al. 2020;Gaynor and Johnson 2020;Sabiston and Kim 2020;Xiong et al. 2020). Furthermore, above reference Swartz et al. (2022) presents different types of constraints (some of them are extensions of Wang et al. (2011) andLiu et al. (2015)) for the design of topology optimized periodic structures like 3d photonic crystals without isolated material components or enclosed pores.
Recently, the authors have developed a new method for imposing connectivity constraints which are based on known results of spectral graph theory. It has been successfully applied to structural design (Donoso et al. 2022) and to electrode design (Donoso et al. 2023). Inspired by those works, we propose here a continuous mathematical model that lets us both characterize and enforce connectivity in any of the two phases (or even in both) involved in a topology optimization problem, and more importantly, regardless of the physical situation. In fact, we also think that it could be extended to ensure connectivity in multi-scale problems (see Wu et al. 2021 for a recent review).
In some ways, our approach presents some similarities with the technique used in Wang et al. (2011) to avoid the appearance of isolated components in PhCWs, and with the VTM itself, but we lay particular emphasis on the fact that our model rests on mathematical basis from spectral theory. The similarities and differences between those approaches will be further discussed below.
The paper is organized as follows. Section 2 presents a continuous mathematical model based on the eigenvalues of the Neumann-Laplacian operator that succeeds in detecting whether a phase in a domain where two phases coexist is connected or not. That can be achieved by solving an appropriate eigenvalue problem. Section 3 shows how the model can be slightly modified to detect and avoid enclosed holes. Section 4 provides a TO-formulation that imposes connectivity over the void phase in structural design. Several numerical examples that corroborate our method are included here. Finally, some conclusions and future work are commented in the last section.

A model for connectivity
It is known (see Courant et al. 1953, Chapter VI, §1.3) that if is a bounded open Lipschitzian set in ℝ N , the first eigenvalue of the Laplacian operator with Neumann boundary condition, the so-called Neumann-Laplacian, is zero; and if is connected, the second eigenvalue is strictly positive. Also if is not connected, we obtain its eigenvalues by collecting and reordering the eigenvalues of each connected component. Therefore, the connectivity of a set can be characterized examining the second eigenvalue of the Neumann-Laplacian: the set is connected if and only if the second eigenvalue is positive. As an example, the first two eigenvalues of the set in Fig. 1a are zero.
Our idea here is to identify the connectivity of a set defined by a density function which is a solution of a minimum compliance problem in TO. Ideally, the solution of such a problem would be a characteristic function of the set where the material is placed, but typically what we obtain is a density function ∈ [0, 1] defined in the reference domain . So, if a sequence of density functions n (illustrated at Fig. 1c) converges (pointwise) to a characteristic function ( Fig. 1b), we would like to establish a relationship between the second eigenvalue of the Neumann-Laplacian in (Fig. 1a) and the limit of the second eigenvalue of the problems where > 0 is fixed. Notice that we need to include this parameter, , in order to have a well posed problem.
First, using (Henrot 2006, Theorem 2.3.3) we have the convergence of the eigenvalues of (1) to the eigenvalues of the problem which is approximated by a density n in c. Gray color representing has been exaggerated to identify the whole domain (a) A set with two connected components.
The following result relates the second eigenvalue of the Neumann-Laplacian in with the second eigenvalue of (2):

Theorem 1 Let be the second eigenvalue of the problem
If denotes the second eigenvalue of (2), then a constant C > 0 exists such that − ≤ C .
Proof Using the Rayleigh quotient (see Henrot 2006), we have the characterization of the second eigenvalue of problem (3): Let * such that the minimum is attained, Let us consider an open bounded domain ′ such that ⊂⊂ ′ . The extension theorem in Sobolev spaces (see Evans 2010, §5.4, Theorem 1) guarantees the existence of u * ∈ H 1 (ℝ N ) such that u * = * in and u * has support in ′ . This extension also satisfies Then u ∈ H 1 ( ) with zero mean value, and therefore, Notice that, as u * = * in then ∇u * = ∇ * in , and, it is clear that ∇u = ∇u * , so the first two terms in the numerator of (5) are equal to ∫ |∇ * | 2 dx. The third term in this numerator and using (4) On the other hand, we can bound the denominator of (5): (note that u * = * in ). So finally, substituting at (5) This finishes the proof. ◻ This result guarantees that the second eigenvalue of a domain and the second eigenvalue of problem (2) is close enough, what can be used to capture the connectivity of a domain.
It is quite common to introduce a shifting method in the eigenvalues' computation, so we will change to − 1 in problems (1-3), and instead of focusing on 0 eigenvalue, we will focus on 1 eigenvalue.
As an example, let us consider be the unit square, and the subset defined by the two black rectangles of Fig. 2. The exact eigenvalues of (3) are known, 1 and the first five (with shifting) correspond to 1, 1, 1 Note that there are two eigenvalues equal to 1 because is not connected. Table 1 shows the numerical approximation of shifted problem (2) for different values of , where a finite element (FE) discretization over two different meshes ( 100 × 100 and 200 × 200 ) have been considered.
While the first eigenvalues are always 1, the second eigenvalues are strictly greater than 1 because of problem (2) is 1 The eigenvalues of the Neumann-Laplacian for a rectangle (0, L) × (0, l) without shifting are m,n = 2 m 2 L 2 + n 2 l 2 , m, n ≥ 0.
set in , which is connected. However, we can see that, for small , the second eigenvalues are close to 1, so we can derive from here that is not connected. However, problem (2) is an ideal situation. In practise, what we have is a problem like (1), with a density with intermediate values between 0 and 1. So, in the same spirit as the SIMP (Solid Isotropic Method with Penalization) does, we introduce the problem where and r is a power that penalizes intermediate densities.
So, the idea of our method relies on imposing that the second eigenvalue of (6) must be strictly greater than 1. Therefore, if is close to a characteristic function of some subset , then has to be connected. In some way, a high value of r introduces an additional requirement for connectivity which favours convergence, but at the same time the rest of eigenvalues are distorted. Table 2 shows the numerical results for problem (6) with different values of r, where is a density obtained by filtering (with a regular conical filter of radius 0.03) the characteristic function of Fig. 2. Note that the second eigenvalue still detects the no connectivity of .

Optimal design without internal voids
The classical compliance minimization problem in topology optimization leads to connected designs because of physical considerations: forces need to be propagated along the structure and isolated features of the material do not contribute to increase stiffness. At the same time, the stiffness requirement typically results in the formation of several enclosed holes that can complicate the fabrication of the final design (specially for 3d problems). Those internal holes in a structure can be seen as different connected components of the set where the material is not placed, or in terms of the density function , the set where = 0 . Therefore, a simple way to avoid the formation of internal holes is imposing the connectivity of that set, or, in terms of the eigenvalue problem (6), forcing the second eigenvalue to be strictly greater than 1, where Notice that we have changed by 1 − as we are now focused on the connectivity of the void phase. However, it is still possible to have a disconnected void phase without any internal hole, as we can see at   where there are no internal holes in the structure, but void phase is not connected. This situation can be easily avoided if we consider a frame around where = 0 . We denote by ♯ to this extended domain (see Fig. 3), and still denote by the extension of by null values outside of . Now, it is clear that the number of enclosed holes is equal to the number of connected components of the void phase minus one. So having no internal holes is equivalent to a connected void phase in ♯ . After discretizing (6) by FE, we consider the discrete problem instead where K and M are the global stiffness and mass matrices, respectively, corresponds to the vector of discretized densities in the extended domain (extended with null values), and ( 2 , 2 ) is the pair eigenvalue and its associated M -orthonormal eigenvector, respectively. Table 3 shows how internal holes are detected in an example of a real minimum compliance structure, and compares the effect of considering or not the extended domain (a small outside frame of void phase of three elements size). In this case, three isolated void areas are identified in the structure delimited by the dashed-line (i.e., without void frame), but they collapse in only two when considering the void frame. Here, as in the rest of numerical experiments that will appear below we have used r = 8 and = 10 −5 .
It is worth mentioning that discrete eigenproblem (8) reminds us to some extent of the one developed by the authors in Donoso et al. (2022) based on graph theory. In that paper, the strategy consisted in considering the centroids of the finite elements in a mesh as nodes, those conforming a graph according to the connections of such elements in the mesh. Void connections were modeled by the product of void densities of adjacent elements. And it was concluded that the multiplicity of the unit eigenvalue lets us know the number of enclosed holes in a finite structure. With the aim of capturing this same behavior from a continuous perspective, we have explored the concept of connected domain through the eigenvalues of the Neumann-Laplacian problem. Undoubtedly, previous work (Donoso et al. 2022) has definitely inspired us to develop the continuous mathematical model proposed here for characterizing void connectivity. In fact, we suspect that the eigenvalue spectrum of both problems (especially the second eigenvalue) should be related, even though the discrete problems are distinct in each case. Regardless, the conclusion is the same in both studies as the second smallest eigenvalue gives us the key to enforce connectivity constraints in TO-problems. In particular, to ensure void connectivity, we must impose that 2 > 1 in ♯ .
From numerics, we need to introduce a high value of r to weaken, in general, the connections in the structures, which in turn let internal holes connect each other. That has to be necessarily complemented with the eigenvalue constraint applied over the extended domain, which forces that no holes remain enclosed in the structure.
On the other hand, as mentioned in the introduction, that above connectivity constraint is similar in appearance to the one used in Wang et al. (2011) to avoid the presence of isolated components in PhCWs. Basically, that goal was achieved by forcing the fundamental free mechanical vibration frequency to stay above a given parameter. And that eigenfrequency was the second eigenvalue of an artificial mechanical eigenproblem. Regarding VTM, even though it is a parabolic problem rather than an eigenproblem, the temperature constraint there may be related to our second eigenvalue, but we dare not make such a claim in a rigorous way.
So, in our opinion, both VTM and Wang et al. (2011) are excellent approaches supported by physical arguments to address the connectivity issue. However, we firmly believe that the Neumann-Laplacian problem must be behind any method that tries to tackle connectivity from a mathematical continuous framework. In this work, we have given the mathematical arguments that bear our proposal.

Formulation and numerical examples
Bearing in mind the considerations of the previous section, the problem formulation for minimum compliance that incorporates void-phase connectivity may be written as where F are the external force vectors and U are the global displacements, so that the objective function is the usual compliance. U is obtained solving the (discrete) elasticity system and K is the global stiffness matrix. The first constraint corresponds to the volume constraint, where v is a vector containing the measure of the elements, V 0 is the (given) volume fraction and | | is the measure of the design domain. And the second constraint corresponds to the second eigenvalue of the problem K(̂E)U = F, that is, the void connectivity constraint.
As usual in TO, the density is filtered with a regular conical filter (Bruns and Tortorelli 2001;Bourdin 2001) that adopts the expression ̃= H R ( ) , where H R is the matrix of the filter for a given filter radius R. With respect to the projection, we have use the one proposed in Guest et al. (2004), with different X , X = E, C . While for the elasticity problem and volume constraint, the values of E are gradually increased as usual to force 0/1 designs, the strategy for the computation of eigenvalues in (10) is to fix C = 16 from the beginning in order to better identify holes. This is why we have different notations for the projected densities: ̂E is the density used in the computation of (9) and volume constraint; and ̂C is the density used in (10), which has been previously extended to the domain ♯ .
In particular, E is gradually increased from 2 to 64 at every 20 iterations, beginning when the iterative process starts to stabilize, which means approximately after 400 iterations. The connectivity constraint is implemented as 2 > 2,min , using a continuation strategy in the lower bound from the very beginning, starting with 2,min = 1.05 and ending with 2,min = 1.15 , increasing that value in 0.01 for each 50 iterations. Also, depending on the continuation strategies considered, we have observed that higher values of r may be required (up to r = 12 in some cases), but in the strategy adopted here, a value of r = 8 is sufficient in most cases. That value remains fixed along the numerical process.
With regard to the sensitivity analysis, only the derivative of an eigenvalue requires special mention, something that is well-known, given by Note that e = 0 on the elements of the extension, so the partial derivatives do not change. Once all the derivatives have been computed, the discretized problem is numerically solved by MMA (the Method of Moving Asymptotes, Svanberg 1987).
An important observation is the fact that 2 never switches to the first (always of value one) eigenvalue due to the parameter , which always forces 2 to stay above 1. Of course, both values may be really close during some iterations, but in such cases, we have corroborated that the derivative formulas considering they are repeated eigenvalues are (a) (b) Fig. 4 Example of a cantilevered beam with a load in the middle point. a minimum compliance design; b the same but forcing that void phase to be connected in the extended domain virtually the same as they are not. It seems to indicate that from a practical point of view they can be treated as single eigenvalues, and that 2 may switch to the third eigenvalue only. The latter should be considered in case we are interested in controlling (not avoiding) the number of holes, a different problem that is in fact a work in progress. Two examples in 2d, a cantilevered beam with a load applied in the middle point and a supported bridge with a vertical load also in the middle point, are used to initially corroborate our method. Optimized structures in both case studies with inner holes and without them are depicted in Figs. 4 and 5, respectively. The mesh sizes used for the 2d examples are 210 × 90 and 300 × 50 , and for both the filter radius was R = 3.6 and the volume fraction V 0 = 0.4.
In Fig. 6, are displayed the snapshots of the topology evolution corresponding to the first example. As in the rest of examples, 2,min is varied from 1.05 to 1.15. Sometimes, the constraint is active in some cases such as in this example, as the value 1.15 is overpassed. The choice of this parameter may negatively affect the robustness of the algorithm judging from our numerical experience. In particular, when increasing 2,min , the set in which optimized solutions are looked for becomes smaller, so, to the best of our knowledge, it is preferable to use values close to 1. It is less restrictive and it satisfies the constraint at the same time. On the other hand, the mathematical meaning of that value is not clear. We suspect that it may be related to the geometry of the set, but we do not currently know concrete results to support this. We consider there are still some pending issues to fully understand the second eigenvalue which are, however, far from the scope of this work. For this reason we will continue to deepen in that invariant in future work.
A more interesting example takes place in 3d, previously studied in Xiong et al. (2020). There, a platform type-structure with null displacements in the red region of the bottom surface is subjected to uniform pressure on the top surface. A quarter structure discretized into 40 × 40 × 30 elements is simulated by symmetry. Three layers of elements on the top of the structure are chosen as a non-design domain. Optimized designs with and without connectivity constraints for V 0 = 0.3 and R = 2 are shown in Fig. 7.

Final remarks
This work proposes a continuous model that imposes void connectivity in structural design, thus avoiding the formation of inner holes in topology optimized structures. The idea behind the method can be used to enforce connectivity in other physical contexts of interest in a straightforward manner, and we plan to explore it in a near future.