An Aggregation-Free Local Volume Fraction Formulation for Topological Design of Porous Structure

Cellular structure can possess superior mechanical properties and low density simultaneously. Additive manufacturing has experienced substantial progress in the past decades, which promotes the popularity of such bone-like structure. This paper proposes a methodology on the topological design of porous structure. For the typical technologies such as the p-norm aggregation and implicit porosity control, the violation of the maximum local volume constraint is inevitable. To this end, the primary optimization problem with bounds of local volume constraints is transformed into unconstrained programming by setting up a sequence of minimization sub-problems in terms of the augmented Lagrangian method. The approximation and algorithm using the concept of moving asymptotes is employed as the optimizer. Several numerical tests are provided to illustrate the effectiveness of the proposed approach in comparison with existing approaches. The effects of the global and local volume percentage, influence radius and mesh discretization on the final designs are investigated. In comparison to existing methods, the proposed method is capable of accurately limiting the upper bound of global and local volume fractions, which opens up new possibilities for additive manufacturing.


Introduction
Porous structures such as bones, coral and sponges have been widely discovered in nature. In recent years, the additive manufacturing (AM) technique provides tremendous opportunities to fabricate cellular structure with complicated shapes. As a powerful tool for innovative design, topology optimization brings the possibility of realizing novel properties beyond traditional materials with porous structure [1], e.g., architected material with programmable Poisson's ratio [2,3], chiral metamaterial [4], bi-mode artificial metamaterial [5] and lattice sandwich composite [6].
Along with the development of AM, much attention has been paid to the optimal design of multi-scale structures in the past few years [7]. In such a multi-scale design, also named concurrent design and hierarchical design, the macrostructure and material constructed by porous structures evolve simultaneously. Since a hierarchical optimization formulation was proposed by Rodrigues et al. [8], there is considerable literature in this aspect [9][10][11][12]. To facilitate the manufacturability, Liu et al. presented a mathematical model for concurrent design on the assumption that the identical microstructure is periodically distributed over the macrostructure [13]. Two separated subsystems, i.e., macrostructure and microstructure, are integrated into one unity by means of numerical homogenization. To realize the resultant design, the connectivity of microstructure needs to be taken into account. Li et al. suggested that the immediate density within a certain range can be classified into the same group [14]. Wang et al. proposed a shape metamorphosis concept with the purpose of generating a sequence of graded microstructures to ensure the connectability [15]. To address the connectivity issues, Liu et al. predefined connective areas in the periodic micro-scale units [16]. Li et al. developed a topology optimization approach for quasi-periodic microstructures by erode-dilate operations [17]. In the framework of the level set method, the variable cutting technique was introduced to obtain the cellular structures [18][19][20], which guarantees connections between the adjacent elements without extra constraints. The lattice structures can be also achieved by the moving morphable components approach [21,22]. Besides the homogenization theory, the substructure and a promising post-process named as de-homogenization can also be the candidate to bridge the macrostructure and porous material [23][24][25][26][27]. For the extensive topic, readers can be referenced to the insightful review [28].
Another way to generate the bone-like structures can be resorting to the porosity control, which may stem from the maximum size constraint in the topology optimization community. Guest and Prévost initially proposed a maximum length dimension restriction by requiring a minimum void volume in a neighborhood surrounding each voxel, in which a penalty term for the objective function and barrier function was introduced [29][30][31]. The maximum size was also integrated into the bi-directional evolutionary structural optimization method [32,33]. Wu et al. proposed infill formulation by imposing the local volume fraction (LVC) on the given region [34]. To handle with vast constraints, the p-norm aggregation function is employed to condense local constraints. Furthermore, this efficient scheme is subsequently extended to the structure containing interior infill wrapped with coated surface, self-supporting structure [35,36], heat transfer structure [37] and multiphase material structure [38]. Liang and Cheng proposed the sequential integer programming to directly solve bounds of linear infill constraints [39]. In the near term, Chen et al. presented a novel parameterized level set method for designing functionally graded cellular structure [40]. Lazarov and Wang presented a band-pass filter, which can result in a maximum size feature [41]. Dou provided an implicit control for local volume percentage by projection manipulation, including two filters followed by two projections and one multiplication [42]. In the case of the explicit and implicit control methods, an apparent defect is that the predefined local constraint cannot be exactly satisfied. Fernández et al. investigated the influence of the p-norm and p-mean function on the distribution and amount of data in the maximum size control [43]. Clausen et al. experimentally proved that the resistance to buckling for the infill structure can be substantially enhanced compared to structures from the traditional minimum compliance formulation [44].
Despite the typical aggregation formulation for porous structure being successfully performed in various engineering applications, the violation of the maximum LVC in existing literature is still a practical problem. In addition, the appropriate value for the aggregation parameter is also a challenge. More recently, the augmented Lagrangian (AL) method was successfully applied to handle millions of stress constraints [45][46][47][48]. Inspired by these facts, we reformulated the LVC constraint by virtue of the AL approach. In comparison with the existing method, the infill structures can be efficiently achieved with the accurate satisfaction of both global and local volume constraints.
The remainder of this paper is indicated as below. Section 2 describes topology optimization formulation for porous structure and the augmented Lagrangian method. Section 3 presents several numerical examples and discussions. Section 4 summarizes the conclusion.

Formulation for Porous Structure Design
It is assumed that the structural design domain is partitioned into total Ne finite elements. The SIMP (Simplified Isotropic Material with Penalization) method is used in this study, with three fields [49]: the design variable field ρ, the filtered field ρ and the projected field ρ The existence and null of the element are denoted by the projected field ρ, which is linked to the design variable field ρ by where and θ control the steepness and the threshold of the projection, respectively. Here ρ e is obtained through the density filtering process [50,51].
where N e,i denotes the neighborhood set of elements that fall inside the ith element's filter radius r min . The weighting factor ω ei is given as here ∆(e, i) represents the distance between the center the eth element and ith element. The Yong's modulus E e is parameterized by ρ e according to the SIMP interpolation [52]: where E 0 is the elastic modulus of solid material and E min (=10 −6 E 0 ) is a minor stiffness to prevent singularity in finite element analysis. η denotes the power exponent of the Young's modulus interpolation law, which is typically set as 3.
The local volume fraction can be defined as following [35] where M e,i is the set of adjacent elements that fall inside the ith element's circle or sphere within the influence radius r max . Aiming to generate the porous structure automatically, the optimization problem can be mathematically formulated as pursuing the minimum compliance subject to local or/and global volume fraction constraints, i.e., where K, u and f are the global stiffness matrix, displacement vector and external force vector, respectively. The Objective function is the static compliance c, which is adopted as the performance of structural stiffness. V f represents the global volume faction while ψ is its corresponding threshold. V e is the eth elemental volume of the solid material. The parameter ϕ denotes the local volume fraction for all voxels. It should be highlighted that the restriction on overall material usage is unnecessary since the local volume fraction is already specified. Under normal circumstances, the value of ψ needs to be satisfied with the condition ψ ≤ ϕ when the total volume fraction constraint is active.

Augmented Lagrangian Method
Currently, a variety of optimizers are available in the topology optimization community. For instance, the compliance minimization problems involving sole volume constraint can be effectively resolved using heuristic optimality criteria and rigorous mathematical programming, e.g., the method of moving asymptotes (MMA) [53] and sequential quadratic programming (SQP) [54,55]. The introduction of the Lagrange multiplier to satisfy the constraint equation is nothing new in the level set method (LSM) [56] and the Bi-directional Evolutionary Structural Optimization (BESO) method [57,58]. LSM and BESO do not account for a quadratic term in the constraint function introduced by the AL method. Furthermore, the AL method introduces considerably more constraint functions than the LSM and BESO methods. However, the AL method becomes appealing when a large scale of local stress constraints can be addressed successfully, according to a recent report [46]. In this work, the Equation (6) can be solved as the solution of the unconstrained program, by establishing a sequence of AL functions in each k-th step [47]: where N Φ and N Ψ are the constants related to number of the local and global volume constraints. The penalization terms Φ (k) (ρ) and Ψ (k) (ρ) are defined as follows [47].
where the vector β (k) = β can be rewritten as In this work, we define g = [g 0 , g 1 , ..., g N e ] T as the constraint vector containing both global and local volume constraint. The Lagrange multiplier vector α and penalty factor γ (k) can be updated as following: where κ(>1) and γ max are the update factor and the upper bound for γ respectively, with the adoption to ensure numerical stability.

Sensitivity Analysis
We will solve the Equation (7) based on the gradient information, and the sensitivity expressions for each term will be derived in this section. The sensitivity of the compliance c with respect to physical density ρ i can be given as [53]: where u i and k 0 are the displacement vector for ith element and local stiffness matrix when physical density ρ i = 1. The sensitivity of the penalization terms Φ (k) (ρ) and Ψ (k) (ρ) with respect to physical density ρ i can be computed as where . The rest of the cases satisfy the following equation: According to the chain rule, the sensitivity of the function ρ e to design variable ρ e can be written as where (·) can represent the individual function c, Φ and Ψ in Equation (7).

MMA Algorithm
In each optimization, the MMA approximation of the Equation (7) at the design variable vector ρ (k) reads [48]: where each individual term is expressed as: during the optimization iterations, the lower asymptote and upper asymptote can be determined as where the parameter s (k) e will be updated dynamically during the optimization iteration; The symbols ρ e and ρ e are the maximum and minimum value of design variable ρ e . When a positive move-limit m is predefined, it yields: in the kth iteration, the design variables are made to satisfy the box constraint [48]: explicitly solving Equation (17) can be accomplished as follows.
here, the e is found as

Numerical Examples
In this section, several two-dimensional numerical tests are presented to prove the effectiveness of the proposed method. The material of the structure has Young's modulus and Poisson's ratio of 1 and 0.3. In all examples, the design domain was composed of the bilinear four-node elements. The projection parameter θ was fixed as 0.5 while started at 1 and was doubled per 50 steps until it reached the maximum 256. The move limit m was set as 0.05 to stabilize iterations. The examples in this paper were built on Matlab, and the program structure was based on this literature [59]. All examples were run on a PC with an Inter i7-9700K processor and 32 GB RAM.

Example 1
The first example optimizes a short cantilever as the benchmark in the existing literature [34]. Figure 1 shows a rectangular design domain, boundary and loading conditions. The overall dimensionless size of planar structure was: length L = 400 and height H = 200. The left side was fully supported, while a concentrated force was applied at the central point on the right edge. The design domain was subdivided into 400 × 200 elements. The local volume fraction and influence radius was set as ϕ = 0.6 and r max = 6, respectively. For comparison purposes, the optimized topologies procured from the proposed method and the p-norm aggregation scheme and the distributions of the ρ field are drawn in Figure 2, respectively. Additionally, the histograms of physical density ρ as well as local volume fraction ρ are displayed in Figure 3.
From Figure 2, we can observe that the similar spongy-like structures can be obtained automatically via two approaches. The maximum of the local volume percentage 0.598 from the proposed method is slightly below the threshold value. Nevertheless, a serious violation of the local volume fraction occurs in the p-norm aggregation formulation. The visual differences in optimized topologies exist at two corners on the right side. Through the suggested method, the material is prone to be spread in the entire design area, which sounds reasonable in the absence of the total material usage. From the above phenomena, it is recommended that the global volume fraction is included in the proposed formulation. at the central point on the right edge. The design domain was subdivided into 400 × 200 elements. The local volume fraction and influence radius was set as φ = 0.6 and rmax = 6, respectively. For comparison purposes, the optimized topologies procured from the proposed method and the p-norm aggregation scheme and the distributions of the ρ field are drawn in Figure 2, respectively. Additionally, the histograms of physical density ρ as well as local volume fraction ρ are displayed in Figure 3.  height H = 200. The left side was fully supported, while a concentrated force was applied at the central point on the right edge. The design domain was subdivided into 400 × 200 elements. The local volume fraction and influence radius was set as φ = 0.6 and rmax = 6, respectively. For comparison purposes, the optimized topologies procured from the proposed method and the p-norm aggregation scheme and the distributions of the ρ field are drawn in Figure 2, respectively. Additionally, the histograms of physical density ρ as well as local volume fraction ρ are displayed in Figure 3.   From Figure 2, we can observe that the similar spongy-like structures can be obtained automatically via two approaches. The maximum of the local volume percentage 0.598 from the proposed method is slightly below the threshold value. Nevertheless, a serious violation of the local volume fraction occurs in the p-norm aggregation formulation. The visual differences in optimized topologies exist at two corners on the right side. Through the suggested method, the material is prone to be spread in the entire design area, which sounds reasonable in the absence of the total material usage. From the above phenomena, it is recommended that the global volume fraction is included in the proposed formulation.
From Figure 2, the local volume fraction distributions demonstrate that the proposed method is capable of achieving a uniform distribution of local volume fraction for each element, whereas the comparison method fails to control the local volume fraction in the structure's critical-support region and load-bearing region. Furthermore, a more thorough analysis as shown in Figure 3a reveals that more than 80% of the elements in the proposed approach have a local volume fraction between 0.5 and 0.6, and none surpass the maximum limit of 0.6. Approximately 24% of elements have a local volume fraction that exceeds the upper limit as illustrated in Figure 3a. The primary reason for this discrepancy is that the suggested approach can precisely and directly restrict each constraint function, while the p-norm aggregation function roughly estimates the maximum value. Moreover, the physical density tends to a distinct 0-1 distribution as illustrated in Figure 3b.

Example 2
The influence of the global volume fraction on final design was investigated in this example. All parameters were identical with those in Example 1. Four upper bounds of the global volume fraction 0.35, 0.40, 0.45 and 0.50 were used, and the corresponding topologies were plotted in Figure 4. Figure 5 depicted the corresponding iterative curves of compliance, global volume fraction and maximum local volume fraction.
From Figure 4, we can see that a group of porous structures are generated by limiting the total material usages. Consistent with engineering intuition, the structure becomes much stiffer when more material is allowed in the design domain. Alongside this, we observed that the iteration curves exhibited a similar trend using various upper bounds From Figure 2, the local volume fraction distributions demonstrate that the proposed method is capable of achieving a uniform distribution of local volume fraction for each element, whereas the comparison method fails to control the local volume fraction in the structure's critical-support region and load-bearing region. Furthermore, a more thorough analysis as shown in Figure 3a reveals that more than 80% of the elements in the proposed approach have a local volume fraction between 0.5 and 0.6, and none surpass the maximum limit of 0.6. Approximately 24% of elements have a local volume fraction that exceeds the upper limit as illustrated in Figure 3a. The primary reason for this discrepancy is that the suggested approach can precisely and directly restrict each constraint function, while the p-norm aggregation function roughly estimates the maximum value. Moreover, the physical density tends to a distinct 0-1 distribution as illustrated in Figure 3b.

Example 2
The influence of the global volume fraction on final design was investigated in this example. All parameters were identical with those in Example 1. Four upper bounds of the global volume fraction 0.35, 0.40, 0.45 and 0.50 were used, and the corresponding topologies were plotted in Figure 4. Figure 5 depicted the corresponding iterative curves of compliance, global volume fraction and maximum local volume fraction.
From Figure 4, we can see that a group of porous structures are generated by limiting the total material usages. Consistent with engineering intuition, the structure becomes much stiffer when more material is allowed in the design domain. Alongside this, we observed that the iteration curves exhibited a similar trend using various upper bounds of the global volume fraction as shown in Figure 5. The compliance decreased dramatically within the first 50 iterations while two curves of the global volume fraction and the maximum local volume fraction upsurged sharply, accompanied by the value of max ρ e reaching to 1. Then, all curves of the compliance underwent mild fluctuations and eventually reached the plateau. In summary, the material can be automatically and efficiently distributed to form a spongy structure by the proposed method meeting the requirement of the imposed constraints.
of the global volume fraction as shown in Figure 5. The compliance decreased dramatically within the first 50 iterations while two curves of the global volume fraction and the maximum local volume fraction upsurged sharply, accompanied by the value of ( ) max e ρ reaching to 1. Then, all curves of the compliance underwent mild fluctuations and eventually reached the plateau. In summary, the material can be automatically and efficiently distributed to form a spongy structure by the proposed method meeting the requirement of the imposed constraints.

Example 3
In this numerical test, the effects of the local volume fraction constraint and influence radius rmax on optimized results are discussed. The global volume fraction was fixed as ψ = 0.45. The upper bound of the local volume fraction ranged from 0.55 to 0.70. Meanwhile two different influence radii rmax = 6 and 12 were adopted. The resulting topologies by various combinations of the parameter φ and rmax are listed in Figure 6. Figure 7 depicts the relationship between compliance, influence radius, and the local volume fraction constraint.

Example 3
In this numerical test, the effects of the local volume fraction constraint and influence radius r max on optimized results are discussed. The global volume fraction was fixed as ψ = 0.45. The upper bound of the local volume fraction ranged from 0.55 to 0.70. Meanwhile two different influence radii r max = 6 and 12 were adopted. The resulting topologies by various combinations of the parameter ϕ and r max are listed in Figure 6. Figure 7 depicts the relationship between compliance, influence radius, and the local volume fraction constraint.   According to Figure 6, when the local volume fraction limit was relaxed, the number of pores in the porous structure decreased. Additionally, the larger influence radius rmax resulted in the following three effects: larger diameters for holes, fewer holes and stronger rods. As observed from Figure 7, the resultant structure became stiffer as the influence radius rmax and the upper bound of the local volume fraction φ grew. In conclusion, the local volume fraction constraint determined the local porosity, while the influence radius rmax played an important role on the void size of topology. These phenomena can be referenced as the design law for porous structures.

Example 4
In addition to the work described in Example 3, which addressed the impact of the local volume fraction restriction and influence radius rmax, a mesh-independence study was performed in Example 4. Four distinct meshes were utilized, including 400 × 200, 800 × 400, 1200 × 600 and 1600 × 800 elements. Two filtering radii and two volume fractions were fixed as rmin = 3, rmax = 6, φ = 0.6 and ψ = 0.45, respectively. The optimized topologies with the resulting compliance, volume fraction are displayed in Figure 8. According to Figure 6, when the local volume fraction limit was relaxed, the number of pores in the porous structure decreased. Additionally, the larger influence radius r max resulted in the following three effects: larger diameters for holes, fewer holes and stronger rods. As observed from Figure 7, the resultant structure became stiffer as the influence radius r max and the upper bound of the local volume fraction ϕ grew. In conclusion, the local volume fraction constraint determined the local porosity, while the influence radius r max played an important role on the void size of topology. These phenomena can be referenced as the design law for porous structures.

Example 4
In addition to the work described in Example 3, which addressed the impact of the local volume fraction restriction and influence radius r max , a mesh-independence study was performed in Example 4. Four distinct meshes were utilized, including 400 × 200, 800 × 400, 1200 × 600 and 1600 × 800 elements. Two filtering radii and two volume fractions were fixed as r min = 3, r max = 6, ϕ = 0.6 and ψ = 0.45, respectively. The optimized topologies with the resulting compliance, volume fraction are displayed in Figure 8.
As seen from Figure 8, the analogous topologies were generated by different levels of discretization, albeit small discrepancies were noticed. For the densest elements case, more than one million local constraints were involved. Hence, it can be inferred that the proposed method can allow for efficient solution of the computationally demanding problem with mesh refinement.

Example 5
A 2D bone-like structure was optimized to exam the viability of the suggested method dealing with irregular shape and mesh. The finite element model consisted of 35,516 elements and 35,912 nodes, and four layers of elements were kept as the passive region. The bottom edge was entirely fixed while both the top left and top right corners were loaded diagonally downward at θ = 45 • as illustrated in Figure 9. The magnitude of the left and right load was 2545 and 1697, respectively. The upper bounds for two volume fractions were set as ϕ = 0.7 and ψ = 0.6. For comparison, the topology optimization regardless of the local volume fraction was also performed. Two resultant topologies are displayed in Figure 10. As seen from Figure 8, the analogous topologies were generated by different levels of discretization, albeit small discrepancies were noticed. For the densest elements case, more than one million local constraints were involved. Hence, it can be inferred that the proposed method can allow for efficient solution of the computationally demanding problem with mesh refinement.

Example 5
A 2D bone-like structure was optimized to exam the viability of the suggested method dealing with irregular shape and mesh. The finite element model consisted of 35,516 elements and 35,912 nodes, and four layers of elements were kept as the passive region. The bottom edge was entirely fixed while both the top left and top right corners were loaded diagonally downward at θ = 45° as illustrated in Figure 9. The magnitude of the left and right load was 2545 and 1697, respectively. The upper bounds for two volume fractions were set as φ = 0.7 and ψ = 0.6. For comparison, the topology optimization regardless of the local volume fraction was also performed. Two resultant topologies are displayed in Figure 10.     As shown in Figure 10, the optimized topologies reflect visual differences, i.e., the material accumulation regardless of the local volume proportion and real cross-section of the bone by porosity control.
The robustness of the final designs concerned to the force direction was verified. The resultant compliance of the bone-like structure sustaining loads in different directions were computed and are plotted in Figure 11. The compliance fluctuates substantially within the scope of the investigation by employing the traditional stiffness optimization. Comparatively, the curve from the suggested method underwent a relatively gradual shift. The optimized results clearly demonstrate that the bone-like structure is insensitive to the direction of external load. As shown in Figure 10, the optimized topologies reflect visual differences, i.e., the material accumulation regardless of the local volume proportion and real cross-section of the bone by porosity control.
The robustness of the final designs concerned to the force direction was verified. The resultant compliance of the bone-like structure sustaining loads in different directions were computed and are plotted in Figure 11. The compliance fluctuates substantially within the scope of the investigation by employing the traditional stiffness optimization. Comparatively, the curve from the suggested method underwent a relatively gradual shift. The optimized results clearly demonstrate that the bone-like structure is insensitive to the direction of external load.

Discussion and Conclusions
This paper proposed an alternative approach to generate the porous structure. Innovatively, masses of local volume fraction constraints were appended to the objective function. The optimum in the primary problem can be obtained by dividing into an array

Discussion and Conclusions
This paper proposed an alternative approach to generate the porous structure. Innovatively, masses of local volume fraction constraints were appended to the objective function. The optimum in the primary problem can be obtained by dividing into an array of sub-problems in the framework of the augmented Lagrangian function. The MMA approximation and solution at each iteration were conducted. The performance of the proposed method was numerically studied in detail, and its effectiveness to produce the infill structure was demonstrated by five 2D numerical examples.
It was found that the vast majority of the predefined LVCs can be strictly satisfied, which can be viewed as a remarkable advantage superior to previous methods that are unable to control the LVCs precisely. In addition, the upper bound of the local volume fraction and the influence radius synthetically determined the local porosity and void size. The proposed method can also provide solutions independent of discretization and insensitive to variation of the external force direction. In future work, the nonlinear or transient topology optimization problem will be furthered based on the current study.
Author Contributions: Conceptualization and writing-original draft preparation, K.L. and Z.C.; validation and software, C.Z. and X.Y.; writing-review and editing, N.S. All authors have read and agreed to the published version of the manuscript.
Funding: This research was funded by Huaneng Technology Funds, grant number HNKJ20-H88.
Data Availability Statement: Data sharing not applicable.