Ambiguous phase assignment of discretized 3D geometries in topology optimization
AAmbiguous phase assignment of discretized 3D geometries intopology optimization
Jorge L. Barrera, Kurt MauteAnn and H. J. Smead Department of Aerospace Engineering Sciences,University of Colorado at Boulder, 3775 Discovery Dr, Boulder, CO 80303, USACorresponding author: [email protected] 25, 2020
Abstract
Level set-based immersed boundary techniques operate on nonconforming meshes while providing acrisp definition of interface and external boundaries. In such techniques, an isocontour of a level set fieldinterpolated from nodal level set values defines a problem’s geometry. If the interface is explicitly tracked,the intersected elements are typically divided into sub-elements to which a phase needs to be assigned.Due to loss of information in the discretization of the level set field, certain geometrical configurationsallow for ambiguous phase assignment of sub-elements, and thus ambiguous definition of the interface.The study presented here focuses on analyzing these topological ambiguities in embedded geometriesconstructed from discretized level set fields on hexahedral meshes. The analysis is performed on three-dimensional problems where several intersection configurations can significantly affect the problem’stopology. This is in contrast to two-dimensional problems where ambiguous topological features existonly in one intersection configuration and identifying and resolving them is straightforward. A setof rules that resolve these ambiguities for two-phase problems is proposed, and algorithms for theirimplementations are provided. The influence of these rules on the evolution of the geometry in theoptimization process is investigated with linear elastic topology optimization problems. These problemsare solved by an explicit level set topology optimization framework that uses the extended finite elementmethod to predict physical responses. This study shows that the choice of a rule to resolve topologicalfeatures can result in drastically different final geometries. However, for the problems studied in thispaper, the performances of the optimized design do not differ.
In topology optimization, a parametrized field that describes a problem’s geometry can freely evolve insearch of an optimal design. The level set method ([1, 2, 3]) provides a convenient framework for thisprocess. Traditionally, immersed analysis techniques smear the interface over elements (e.g., the Ersatzmaterial method [4]). However, a rapidly growing class of these techniques are able to capture and trackthe interface explicitly (e.g., cut finite element method [5], extended finite element method [6]). In suchtechniques, embedding the geometry in a computational domain via a discretized level set field (LSF) mayallow for ambiguous definitions of the discretized geometry at any stage of the design optimization process.Therefore, methods for detecting and resolving these ambiguities are needed. Yet, with only a few exceptions([7, 8, 9]), work on level set topology optimization omits mentioning and addressing this issue.Ambiguities of discretized geometries can be resolved via user-defined rules. The effect of a given ruleon the optimization process and the optimized designs, i.e. whether some user-defined rules can lead tosuboptimal designs, is still an open question. This scenario is illustrated in the 2D schematics of Fig. 1,where two different rules, R and R , are used to resolve geometric ambiguities. The sketch shows that,starting from the same initial design, not only different intermediate designs are possible, but also thatthe optimized designs may differ in geometry and performance. Hence, in this work we investigate the1 a r X i v : . [ c s . C E ] F e b nfluence of these rules used to determine ambiguous geometrical configurations in the context of topologyoptimization.Figure 1: Probable effect of phase assignment using either rule R or R to resolve ambiguous geometricalfeatures in the optimization process.In this paper, we focus on level set topology optimization approaches that use fixed background meshes todiscretize the geometry and/or the state variable fields. These approaches typically split elements in whichthe interface is embedded ([3, 10]). The intersected elements are decomposed into sub-elements such that aface in 3D, or an edge in 2D, aligns with the interface defined by the isocontour of the LSF. A solid phaseneeds to be assigned to each of the generated sub-elements. However, certain geometrical configurationssuffer from ambiguous definitions of the interface, and thus the phase cannot be uniquely determined for allsub-elements in these cases. This scenario is displayed in the zoomed-in region of the 2D embedded topologyin Fig. 2. The level set values at the corner nodes of the center background element allow for two differentinterpretations of the geometry. Although both configurations are consistent with the discretized LSF, oneof them creates connections between the void inclusions.In a 2D domain discretized by quadrilateral elements, ambiguous topological features exist only if anelement is intersected by an isocontour more than once. In 3D domains composed of hexahedral backgroundelements, this issue is also found in simpler intersection configurations that describe a single interface, andis aggravated if multiple interfaces are discretized within an element ([11, 12, 13]). Furthermore, specialcare may be needed for neighboring background elements that contain sub-elements with ambiguous phasesassignment. Inadequate handling of sub-elements that share common faces across adjacent backgroundelements may introduce inconsistencies in the geometry description ([8]).Although ambiguous geometrical configurations are avoided if a LSF is discretized using tetrahedralelements, hexahedral elements remain the preferred choice ([14]). This is in part because hexahedral elementsprovide a good compromise between low computational costs and accurate approximations. Furthermore, LSF
Discretized domain
Solid
Void
Topology Ambiguities in Level Set Methods
The layout of a two-phase solid-void problem in a design domain, Ω, is described by a LSF, φ , as follows: φ ( X ) > , ∀ X ∈ Ω S ,< , ∀ X ∈ Ω V , = 0 , ∀ X ∈ Γ S,V . (1)The design domain is subdivided into the solid and void phases, Ω S and Ω V , respectively, such that Ω S ∪ Ω V =Ω. The interface, denoted by Γ S,V , is identified by the zero level set isocontour.The LSF is discretized by tri-linear shape function on a structured mesh, H . This mesh is composedof regular eight-node hexahedra that are classified into two groups, H U and H I , with H U ∪ H I = H . Thesubdomain H U contains all the un-intersected hexahedra. A hexahedron, H i , is part of this set if all theirnodal level set values are either positive or negative. Furthermore, the subset H I consists of intersectedhexahedra defined as all hexahedra for which at least one of their eight nodal level set values differs in signfrom the remaining nodal level set values.All intersected hexahedral elements, i.e. H i ∈ H I , are split into tetrahedral sub-elements via a tetra-hedralization algorithm. The vertices of the tetrahedra are constructed from vertices of the backgroundelement and intersection points along the edges. Hence, some of the faces of these generated tetrahedracontain intersection points only. Tetrahedralization algorithms employed for this subdivision receive thenodal level set values and the intersection points of a hexahedron H i as input and return the geometricalinformation of the set of generated tetrahedra, T H i . In this work, all the generated tetrahedra are containedin the set T = { T H i : H i ∈ H I } . The Delaunay method ([31, 32]) is used for tetrahedralizing the intersectedhexahedral elements in this study.The intersection points along the edges are defined by the points with a zero level set value. The levelset values along an edge are linearly interpolated from the level set values of the edge’s corner nodes. Thus,the linear interpolation of the LSF limits the possible configurations to one intersection per edge. Since eachnode can have a positive or negative level set value, a total of 2 = 256 nodal phase assignment configurations(NPACs) are possible in an eight-node hexahedral element. After accounting for reflective and rotationalsymmetries, the number of unique NPACs for intersected elements is reduced to the 14 cases shown in Fig. 3.The reader is referred to [11, 12, 13] for a case-by-case analysis, as well as detailed visual schematics.Figure 3: All unique NPACs using a discretized LSF onto an intersected eight-node hexahedron.In this work, the discretized geometry is constructed by connecting intersection points along intersectededges. Thus, a small discrepancy between the discretized geometry and the approximated level set field usinglinear interpolation functions is expected. The left side of Fig. 4 illustrates this in a simple 2D schematic.It can be seen that the zero isocontour of the LSF, described by a blue dashed curve, is approximated by astraight blue line between two intersected edges. 4 A
LSF isocontours
Intersection
Figure 4: Two contiguous hexahedra which discretized geometry present ambiguous phase assignment (cen-ter). The highlighted faces show (left) the difference between approximated level set field and discretizedgeometry, and (right) level set field contours on such face and its saddle point.Once intersected hexahedra have been split into tetrahedra, one or multiple interfaces must be definedto approximate the problem’s geometry. This implies that either a solid or a void phase must be assigned toevery tetrahedron, T i ∈ T . In some cases, however, one or more tetrahedra can be defined as solid or voidwhile keeping consistency with the associated NPAC. The lack of uniqueness inherent in such cases leads tomultiple possible interface definitions. In other words, different discretized domains can be generated fromthe same discretized LSF depending on the phase assigned to these ambiguous tetrahedra, as shown in Fig. 2in 2D. An ambiguous tetrahedron, AT, is defined as any tetrahedron whose vertices are all intersection points.Hence, a minimum of four intersection points is needed for an AT to exist. For this reason, all but theNPAC-1 configuration can contain ambiguities; see Fig. 3. Despite the number and shape of the generatedtetrahedra may change depending on the tetrahedralization algorithm employed ([32, 33]), ATs cannot beavoided.To simplify the analysis of topology ambiguities, the NPACs described in Fig. 3 are categorized as shownin Fig. 5. Examples of continuous LSFs to be embedded in a discretized domain are shown on the left column.The discretized geometries are shown in dark gray, together with the ATs highlighted in light gray in thecenter column. The right column shows only the ATs for better visualization. In the first case, illustrated inFig. 5(a), no ATs are present since only three intersection nodes exist. Fig. 5(b) shows the most innocuousNPAC with ATs, which occurs in hexahedra with a single interface. In this configuration, assigning differentphases to the ATs results in different shapes but not topologies. If the background element is intersectedmore than once, as shown in Fig. 5(c), then ATs can connect/disconnect interfaces, and thus alter thetopology. Note that in both Figs. 5(b) and (c), none of ATs faces is part of the background hexahedronfaces. We classify them as Internal Ambiguous Tetrahedra, herein labelled as IATs. However, NPACs canalso include ATs with one of their faces overlapped with one of the hexahedron’s faces, as is depicted inFig. 5(d). We term these Boundary Ambiguous Tetrahedra as BATs.
A key characteristic of this type of ambiguities is that at least one of the faces of an IAT lays on aninterface. If the NPAC leads to only one interface, all ambiguous sub-elements are IATs, as shown inFig. 5(b). Furthermore, the difference between the interface geometry approximated by the IATs and theone described by the linearly interpolated LSF is typically negligible.5igure 5: Illustrations of the four cases of intersection configurations for a tetrahedralized hexahedron: (a)no ambiguous sub-elements; (b) IATs with one interface; (c) IATs with multiple interfaces; (d) IATs andBATs with multiple interfaces.
This type of AT can only exist in intersected hexahedra that are able to describe multiple interfaces. Fur-thermore, BATs are always accompanied by IATs. A well-known issue of NPACs with BATs is guaranteeingconsistency in phase assignment across two adjacent hexahedral background elements that share a face onwhich a BAT is defined; see [25, 8]. The sketch at the center of Fig. 4 shows this scenario. To achieve aconsistent phase assignment, BATs sharing a hexahedron face must have the same phase.
Internal ambiguities are ubiquitous and, in most cases, harmless since they do not alter significantly thegeometry. However, hexahedral background elements with multiple interfaces can have IATs which phaseassignment may have a strong influence on the topology of a given design. In most of those scenarios, theanalysis has to include BATs.Inappropriate treatment of ambiguous configurations that combine IATs and BATs can result in undesiredgeometries, especially when thin features are present in the design. Shell-like structures, which are notuncommon in engineering design problems, can be greatly affected since the appearance of dents may be6avored. This is illustrated in the example in Fig. 6. Here, the influence of ATs on the geometry representationof a spherical shell is shown for three different wall thicknesses. The shell has a constant inner radius andits thickness is reduced by varying the external radius. An eighth of the 3 × × × ×
10 hexahedral background mesh and is highlighted in dark gray for visualization purposes. Inthis example, void phase is assigned to all ATs. The discretized geometry, together with the ratio betweenvolumes of AT and solid subdomain, are provided for the three configurations on the right side of Fig. 6. Itcan be observed that once the difference between the inner and outer radii is smaller than the edge length of abackground element, the volume of the ambiguous subdomain increases up to 26.71% of the solid subdomain.Note that the extreme misrepresentation of the spherical shell geometry can be avoided by assigning a solidphase to all ATs. However, assigning the same phase to ATs is not necessarily the most beneficial option forall intersection configurations, as is shown in Section 5.To assess whether it is more beneficial in the optimization process to use a simple criterion, such as theone used for the spherical shell, or geometrical information, a group of rules that systematically assigns aphase to AT is presented in next section. Ω V
The ambiguities described in Section 2.2 can be resolved by applying phase assignment rules either locally(i.e., separately for each hexahedral background element) or globally (i.e., on groups of adjacent hexahedra).The local and global rules explored in this paper are denoted by L and G , respectively. These rules collecta subset of techniques that allow for the free evolution of the geometry in a shape or topology optimizationprocess, but this subset is by no means complete. The reader is referred to [14, 25, 27, 30] for more informationon rules used to this end. L ) The local rules described in this section resolve IATs and BATs separately, using one criterion for eachambiguity type. Consistent phase assignment across hexahedra (i.e., BAT phase assignment) is achievedby enforcing the asymptotic decider [25]. This technique uses the nodal level set values of the hexahedronface containing a face of a BAT to compute the level set value at the saddle point, φ sp , assuming a bi-linerinterpolation of the level set field over the hexahedron face. The right side of Fig. 4 shows a schematic of thelocation of φ sp on a face that contains BATs. The level set value at the saddle point is computed as follows: φ sp = φ A φ C − φ B φ D φ A + φ B + φ C + φ D , (2)where φ A , φ B , φ C , and φ D correspond to the nodal level set values of the sharing hexahedron face. The signof φ sp determines the phase to be assigned. If φ sp ≥
0, then all BATs in the adjacent hexahedra sharing the7ace are assigned the solid phase. Conversely, if φ sp <
0, BATs are assigned the void phase. The reader isreferred to [25] for details of this technique.After BATs have been resolved via the asymptotic decider, the phase of the remaining IATs is determinedusing one of the rules described below (in order of increasing complexity in terms of implementation). L Ω S / L Ω V ) In this rule, all IATs are assigned the same phase throughout the optimization process. If ambiguous sub-elements are assigned the solid phase, the rule is denoted as L Ω S ; conversely, if void phase is assigned, it isnamed L Ω V . This approach requires no additional computations for phase determination, as can be seen inAlgorithm 1. Algorithm 1:
Resolving IATs using the L scheme// Loop over all intersected hexahedra of the background mesh for H i ∈ H I do // Check if any of the generated tetrahedral sub-elements are ambiguous if IAT (cid:54) = {} then // Assign predefined phase to all IAT p = Ω S / Ω V for L Ω S / L Ω V L ) An easy-to-implement alternative that promotes a smoother evolution of the design during the optimizationprocess is choosing the phase assigned to the IATs in the current design iteration, D it , based on an elementalphase, p H i , function of the previous design iteration, D it − . If, for example, a hexahedron with the configu-ration NPAC-4 shown in Fig. 5(b) is part of the void phase at D it − , this rule prevents connecting the twointerfaces by setting the IATs to void phase. In this scheme, all internal ambiguities are assumed to be solidphase in the first design iteration (i.e., at D it = 1) to promote a topology preserving evolution of the designat the initial stage of the optimization process. An implementation of this logic is detailed in Algorithm 2. Algorithm 2:
Resolving IATs using the L scheme// Loop over all hexahedra of the background mesh for H i ∈ H do // Check if element is not intersected to initialize elemental phase p H i if H i ∈ H U then // Store phase of un-intersected hexahedron p H i = Ω S ∨ p H i = Ω V else ( H i ∈ H I ) if D it = 1 then // If this is the first design iteration, set elemental phase to solid phase p H i = Ω S // Assign phase p H i to the ambiguous sub-elements if IAT (cid:54) = {} then IATs phase = p H i L ) Expanding the idea of the asymptotic decider to 3D, IATs can be resolved based on the sign of a singleinterpolated level set value in H i . The decider point (equivalent to the saddle point of the asymptoticdecider) is defined as the average of the centroids of the IAT within a background element. This is illustrated8n 2D in Fig. 7, where the centroid of the ambiguous subdomain, ¯ c IAT , is used to determine a phase. Notethat information of the previous optimization step is not needed in this rule. Also, in hexahedra with asingle interface, the phase is assigned in accordance to the unambiguous LSF used to construct the interface.A scheme that implements this rule is shown in Algorithm 3.Figure 7: 2D schematics of the L rule. Algorithm 3:
Resolving IATs using the L scheme// Loop over all intersected hexahedra of the background mesh for H i ∈ H I do // Check if hexahedron contains ambiguous sub-elements if IAT (cid:54) = {} then // Compute centroids of every IAT using the coordinates of their vertices, X IAT vi c IAT = N v =4 (cid:88) i =1 X IAT vi /N v // Compute the average of the centroids of the N IAT ambiguous sub-elements¯ c IAT = N IAT (cid:88) i =1 c IAT i /N IAT // Interpolate the level set value at centroid using linear shape functions, N , and the nodallevel set values, ˆ φ H i , of the containing hexahedron φ ¯ c IAT = N (¯ c IAT ) ˆ φ H i // Assign elemental phase based on the sign of the level set value if ( φ ¯ c IAT ≥ p H i = Ω S else p H i = Ω V // Assign a phase to the ambiguous sub-elementsIATs phase = p H i L A + / L A − ) To avoid interpolating the LSF to resolve the IATs, a rule based on a geometric indicator is proposed. In thisrule, the shared surface area between each IAT and its unambiguous neighboring tetrahedra are computedper phase. The shared surface areas of two neighboring IATs are omitted from this calculation. The phaseof the unambiguous tetrahedra with whom the IATs share the maximum ( L A + ) or minimum ( L A − ) surfaceis chosen for all IAT. Figure 8 shows a 2D schematic of this rule.Note that implementing this scheme requires knowledge of the phase of adjacent tetrahedra. Hence, theconnectivity of the generated tetrahedra must be constructed per intersected hexahedron. An approach thatemploys this rule is summarized in Algorithm 4. G ) In this set of rules, the asymptotic decider is no longer used to enforce consistent phase assignment to BATsacross intersected hexahedra. Here both IAT and BAT are resolved using the same criterion. Consequently,9 T
Resolving IATs using the L scheme// Loop over all intersected hexahedra of the background mesh for H i ∈ H I do // Check if hexahedron contains ambiguous sub-elements if IAT (cid:54) = {} then // Initialize areas for both phases A Ω S = A Ω V = 0// Loop over all ambiguous tetrahedra in intersected hexahedron for T H i j ∈ IAT do // Loop over the faces, F k , of current ambiguous tetrahedron for F k ∈ T H i j do // Get tetrahedra sharing current face T H i F k and T H i F k sharing F k // Skip face computations if both tetrahedra are ambiguous. if ( T H i F k & T H i F k ) ∈ IAT s then continue else // Compute face area, A F i // Update area of shared tetrahedra based on the phase of unambiguoustetrahedron( A Ω S = A Ω S + A F i ) ∨ ( A Ω V = A Ω V + A F i )// Assign elemental phase based on shared area if A Ω S > A Ω V then p H i = Ω S for L A + , or p H i = Ω V for L A − else p H i = Ω V for L A + , or p H i = Ω S for L A − // Assign a phase to the ambiguous sub-elementsIATs phase = p H i chains of connected ATs across multiple elements can be formed. The criteria described below are applied ei-ther on the entire domain ( G ) or per cluster of connected ambiguous tetrahedra connected across intersectedtetrahedra ( G ). G Ω S / G Ω V ) This rule enforces either solid ( G Ω S ), or void ( G Ω V ) phase to all ambiguous tetrahedra. Similar to L , noadditional implementation is needed once the ATs have been identified. However, in contrast to L , hereBATs are also included. Note that for this rule the choice of using G Ω S or G Ω V may have a major impact onthe resulting geometry compared to the other rules presented in this paper.10 .2.2 Shared area between ambiguous and uniquely defined tetrahedra per phase ( G A + / G A − ) This global rule is an extension of the L local rule and is enforced on clusters of connected tetrahedra.These clusters can exist across intersected hexahedra instead of being contained within a single hexahedron;see Section 3.1.4. Hence, knowledge of the connectivity of all tetrahedra, which changes for every designiteration, is needed to identify these clusters. A graph search scheme to generate this connectivity thatavoids issues associated with memory and/or computational time consumption is provided in Algorithm5. Note that implementing this rule is laborious, especially in parallelized software. Existing third partylibraries, such as MOAB ([34]) or the STK package that is part of Trilinos ([35]), provide robust frameworksto manage mesh databases and can be used to ease the search of connected ATs. In terms of implementation, local rules L , L and L , and global rule G , are more attractive as they requirelittle effort. Additionally, these rules have a low computational overhead. However, rules L and G providea more consistent/smooth geometry description at the cost of higher computational effort. The advantage ofthese last two rules comes at the cost of involved implementations of problem specific geometric indicatorsto resolve ambiguities.Local rules L , L and L allow for imperfections in the geometry description due to disparity in phaseassignment between BAT and IAT. Either solid floating pieces embedded in a mostly void subdomain, orsmall holes in a mostly solid subdomain, are likely to occur in hexahedra with multiple interfaces, see [28]. Inthe former case, the numerical technique used for the physical analysis must include appropriate treatment ofthese disconnected pieces, such as consistent generalized enrichment when using XFEM and adding stiffnessto avoid zero energy modes; see Section 4.3. The local rule L does not suffer from the issues describedabove since the phase of the interior clusters is determined from the phase of the surrounding unambiguoustetrahedra.Global rules have a more pronounced effect on the resulting discretized geometry compared to local rules.Although only rarely observed, assigning different phases to large clusters of connected ATs can significantlyalter the topology. Differences in topologies are less evident if local rules are used because of the asymptoticdecider. The influence (if any) of using the rules described in this section on the smoothness of the evolutionof the geometry in a topology optimization process, as well as the performance of the optimized designs, hasbeen unexplored to date. In this paper, a study of these rules is presented in Section 5. The design, i.e. the geometry, is controlled by a vector of design variables, s := { s ∈ IR N s | φ low ≤ s i ≤ φ up , i = 1 , ..., N s } bounded between user-defined lower, φ low , and upper, φ up , bounds. A design variable isassigned to each node of the structured hexahedral background mesh; thus, N s equals the number of nodesof the background mesh; see [3, 10] for details.To increase numerical stability and enhance convergence of the optimization problem, the explicit re-lationship between the design variables and the nodal coefficients of the discretized LSF is defined by thelinear filter described in [36]. This distance-based filter is defined as: φ i = 1 N rf (cid:88) j =1 w ij N rf (cid:88) j =1 w ij s j , w ij = max (0 , r f − | X i − X j | ) , (3)where N r f is the number of nodes within the filter radius r f ; and | X i − X j | is the Euclidean distance betweennodes i and j. Within each hexahedral background element, the LSF is interpolated tri-linearly. Furthermore,the regularization scheme described below in Section 4.1 is adopted to promote a uniform spatial gradientof the LSF near the solid-void interface while converging to either a positive or negative target value awayfrom the interface. 11 lgorithm 5: Resolving ATs using the G scheme// Initialize container of clusters, C = []// Initialize flag for N T ∈ H I tetrahedra of intersected background elements, P : P i = false, i = 1 , ..., N T ∈ H I // Loop over all intersected hexahedra of the background mesh for All tetrahedra faces, F T ∈ H I do // Get tetrahedra sharing current face, F i ,// Check if any tetrahedra sharing F i is ambiguous and has not been processed if ( T F i ∈ IAT ) & ( P T F i = false) then // Initialize a face list with current face, C F = [ F i ]// Initialize an empty cluster of AT, C AT = []// Process tetrahedra in face list while C F (cid:54) = [] do // Loop over tetrahedra sharing current face, T j ∈ F i for T j ∈ F i do // If tetrahedron is ambiguous, process it if T j ∈ AT then // Add tetrahedron to current cluster C AT = [ C AT , T j ]// Loop over faces in current tetrahedron, F k ∈ T j for F k ∈ T j do // Add faces to list C F = [ C F , F k ]// Flag tetrahedron as processed to avoid repeating searches in first loop P T j = true// Remove current face in C F list// If a new cluster was generated from the processed tetrahedra if C AT (cid:54) = [] then // Add new cluster to container of clusters C = [ C , C AT ]// After all clusters of connected ATs have been generated, loop over them for C AT i ∈ C do // Initialize areas for both phases: A Ω S = A Ω V = 0// Loop over faces of current cluster for F j ∈ C AT i do // Get tetrahedra sharing current face, T k ∈ F j // Check if only one of the shared tetrahedra is ambiguous if Only one F j ∈ AT then // Compute face area, A F j // Update area of phase of shared tetrahedra (either p = Ω S or p = Ω V ) A p = A p + A F j // Assign phase based on shared area if A Ω S > A Ω V then p = Ω S for G A + , or p = Ω V for G A − else p = Ω V for G A + , or p = Ω S for G A − // Loop over all ambiguous sub-elements in current cluster for T j in C ATi do Assign phase p .1 Optimization problem formulation The optimization problems considered in this work are formulated as:min s z ( s , u ) = w F ( s , u ) + w P P er ( s ) + w P Reg ( s ) s.t. : g i ( s , u ) ≤ , i = 1 , ..., N g . u ∈ IR N u , s ∈ Π = { IR N s | φ low ≤ s i ≤ φ up , i = 1 , ..., N s } (4)The first component of the objective, z , represents the quantity of interest to be minimized, F , such as strainenergy, target displacement, and mass. The second term is the normalized perimeter control penalty, addedto avoid the emergence of irregular geometric features: P P er = (cid:82) Γ S,V dA (cid:82) Γ S,V dA , (5)which is normalized by the perimeter of the initial design, Γ S,V . The last component is a penalty thataims at regularizing the LSF. This penalty ensures (i) an approximately uniform spatial gradient near thesolid-void interface; and (ii) that the LSF away from the interface convergence to target upper and lowerbounds. The normalized formulation of P Reg reads: P Reg = 1 (cid:90) Γ S,V dA (cid:20) (cid:90) Ω D α (1 − w ) (cid:16) φ ( s ) / ˜ φ − sign ( φ ) (cid:17) dV + (cid:90) Ω D α w (cid:16) ||∇ φ ( s ) || / ||∇ ˜ φ || − (cid:17) dV + (cid:90) Ω D α (1 − w )( (cid:107)∇ φ ( s ) (cid:107) ) dV (cid:21) , (6)with α , α , and α as weighting factors. The first term penalizes the difference between the LSF, φ ( s ),and a target value, ˜ φ , away from the interface. The second term penalizes the difference between the normof the spatial gradient of the LSF, ∇ φ ( s ), and a target gradient norm, ∇ ˜ φ , in the vicinity of the interface.The third term enforces a gradient of zero away from the interface. The parameter w is a measure of thedistance of a node to the interface defined as: w = e − γ w ( I ( X ) /I max − , (7)where the parameter γ w determines the proximity of the interface. In this regularization scheme, the distanceof a node to the interface is approximated by a neighborhood-level field I ( X ), as described in [15], andillustrated in Fig. 9. The I max parameter represents a measure of the of the region along the interface overwhich the approximation I ( X ) is evaluated. The higher I max , more background elements surrounding theinterface are considered. In this work, we set ˜ φ = 3 h ( sign ( φ )), ∇ ˜ φ = 1, I max = 1, and the weighting factors α i in the range of [0 . − . w i in the objective function are selected such that thecontributions of the regularization and perimeter control penalty components are significantly lower than F .The design needs to satisfy a set of N g problem dependent inequality constrains, [ g , ..., g N g ] that considerlimits on, e.g., mass, maximum allowable stress, and minimum eigenvalue. The number of design variablesare denoted by N s and the number of state variables is N u ; see Eq. 4. The state variables u are governedby a set of discretized partial differential equations, which are satisfied at each design in the optimizationprocess.The parameter optimization problems are solved using a gradient-based nonlinear programming scheme,namely the Globally Convergent Method of Moving Asymptotes (GCMMA), with no inner iterations; see[37]. The design sensitivities of the objective function with respect to the design variables are defined as dzd s = ∂z∂ s + ∂z∂ u d u d s . (8)13
The state variable vector, ˆ u i ( X ), at node i of a two-material problem (phase I, phase II) is approximated as:ˆ u i ( X ) = L (cid:88) l =1 (cid:18) H ( − φ ( X )) N eN (cid:88) k =1 N k ( X ) δ klq u k, Ω S il + H ( φ ( X )) N eN (cid:88) k =1 N k ( X ) δ klq u k, Ω V il (cid:19) , (11)where H is the Heaviside function depending on the LSF: H ( φ ) = (cid:40) , ∀ φ ( X ) > , ∀ φ ( X ) < . (12)In Eq. 11, L denotes the maximum number of enrichment levels, N k ( X ) represents the nodal shapefunction, and δ klq is the Kronecker delta which selects the active enrichment level q for node k . The Kroneckerdelta ensures that displacements at node k are only interpolated by a single set of degrees of freedom definedat node position X such that the partition of unity principle is satisfied ([41]). The number of nodes perelement is given by N eN . The reader is referred to [42, 43, 44] for an in-depth explanation of the generalizedHeaviside enrichment strategy employed.Note that other immersed boundary techniques with a crisp interface definition, such as CutFem [5] andHIFEM [45], are equally applicable without loss of generality.14 .3.2 Governing equations Static equilibrium is described by the following residual equation augmented with stabilization terms: R = R U + R N Γ + R G Γ + R S , (13)in which the weak form of the linear elastic governing equation, R U , is defined as: R U = (cid:90) Ω S δ (cid:15) : σ dV − (cid:90) Γ Ω SN δ u T ν dA. (14)The displacement field and the test function are denoted by u and δ u , respectively. Traction forces, T ν , areapplied along the boundary Γ Ω S ν . The material tensor for isotropic linear elasticity, D , together with theinfinitesimal strain tensor, (cid:15) = ( ∇ u + ∇ u T ), define the Cauchy stress σ = D : (cid:15) .The unsymmetric version of Nitsche’s method is employed to weakly enforce Dirichlet boundary condi-tions ([46, 47]). These conditions are applied through the following residual component: R N Γ = − (cid:90) Γ (cid:74) δ u (cid:75) σ · ν dA + (cid:90) Γ δ ( σ · ν ) (cid:74) u (cid:75) dA + γ N E/h (cid:90) Γ (cid:74) δ u (cid:75)(cid:74) u (cid:75) dA. (15)The normal vector on a domain boundary is denoted by ν . The jump operator (cid:74) · (cid:75) is defined as (cid:74) u (cid:75) = u − ¯ u , (cid:74) δ u (cid:75) = δ u − δ ¯ u . (16)The integrals of Eq. 15 correspond to the standard consistency, adjoint consistency, and the Nitsche penaltyterms, respectively. The third term is scaled by the Young’s modulus, E , the element size, h , and the penaltyfactor γ N . This last term provides additional control over the accuracy at which a boundary condition (BC)is enforced. Our numerical experiments show that choosing a penalty factor within a range of γ N = [10 , u and δ u by tri-linear shape functions, ill-conditioning is mitigated by applyingthe following penalty: R G Γ = Ehγ G (cid:88)(cid:124)(cid:123)(cid:122)(cid:125) F ∈ F cut (cid:90) F (cid:74) ∇ δ u · ν e (cid:75)(cid:74) ∇ u · ν e (cid:75) dA. (17)Element faces in the vicinity of the interface for which at least one of the two adjacent elements is intersectedare included in F cut , as explained in [49]. Normals of these element faces are denoted by ν e . The influenceof the ghost penalty term presented above is controlled by the penalty factor γ G , which is usually chosen tobe within γ G = [0 . , .
01] based on numerical experiments; see [50, 51].To suppress rigid body motion of disconnected solid subdomains that may emerge and develop, selectivestructural springs are added to the governing equations; see [50, 52]. An additional stiffness term is appliedonly to solid subdomains that are disconnected from the support via the following residual component: R S = γ S E/h (cid:90) Ω D δ u · u dV. (18)The parameter γ S denotes the spring stiffness scaling and is non-zero only for the free-floating pieces ofsolid material. These pieces are identified by an auxiliary indicator field constructed from the solution of adiffusion problem, as described in [52].In the final example, a gradient stabilized scalar stress field, τ ( X ), is post-processed via the XFEMinformed smoothing procedure described in [53]; see Section 5.3. The additional set of state variables iscomputed by solving the residual equation: R = R τ = (cid:90) Ω I δτ ( τ − σ V M ) dV + h γ τ (cid:88)(cid:124)(cid:123)(cid:122)(cid:125) F ∈ F cut (cid:90) F (cid:74) δ ∇ τ (cid:75)(cid:74) ∇ τ (cid:75) dA. (19)The field σ V M ( X ) and the parameter γ τ represent the von Mises stress and the ghost penalty weight,respectively. Spatial oscillations of τ are mitigated by penalizing the jump in spatial stress gradients acrosselemental faces (second term in Eq. 19). 15 Numerical Examples
The local and global rules presented in Section 3 are studied in this section to assess their influence on theevolution of the topology, and in the optimized design. Furthermore, the influence of the initial seedingpattern of holes on the optimization process is investigated for different rules. These rules are tested onstructural linear elastic single material, solid-void problems solved by the explicit level set XFEM topologyoptimization framework detailed in Section 4. A strain energy optimization problem is formulated in Example1, minimizing the displacement in a portion of the structure is investigated in Example 2, and in thelast example, the objective consists of minimizing the total mass while satisfying a stress constraint. Theoptimization problem is considered converged if the constraints are satisfied and the relative change inobjective between two consecutive design iterations is less than 1 × − . The initial seeding of the designdomains with holes is constructed using patterns of cuboid or spherical primitives.The governing equations are discretized using the XFEM approach outlined in Section 4.3. Relevantparameters used for the numerical examples 1 and 2 are listed in Table 1 in self-consistent units. Table 5summarizes the parameters used in the last example in SI units. The example problems consist of a one-waycoupled set of governing equations; i.e., the diffusion problem describing the auxiliary indicator field, thestabilized linear elasticity equations (Eq. 13) and, for the third problem, the stabilized stress projectionequations (Eq. 19). The systems of equations of the first and second examples are solved using an ILUpreconditioner ([54]). The third example was solved iteratively via a Trilinos algebraic multigrid solver([55]). Table 1: Parameters common to all numerical examples function of the element size h .Parameter ValueTarget LSF ˜ φ = 2 . h LSF regularization control γ P Reg = − log (0 . f r = 1 . h Nitsche penalty factor γ N = 100 . γ G = 0 . γ S = 1x10 − Stress ghost penalty factor (Ex. 3) γ τ = 0 . A beam problem is used to investigate the effect of the topology consistency rules detailed in Section 3as the mesh is refined. In this example, we focus on analyzing differences in the optimized designs thattypically favor thin-walled shear-webs when using a level set topology optimization approach. Three levelsof refinement are considered.The problem setup is shown in Fig. 10. The design domain of size 240 × ×
40 is simply supported at allfour corners of the bottom face, as shown on the left side of Fig. 10. A traction load T X = − . ×
10 area. Only one quarter of the design domain is simulated byapplying symmetry boundary conditions along the X − X and X − X planes at the center of the designdomain. The loading and support regions highlighted in dark gray are excluded from the design domain.Structured hexahedral meshes with element sizes of h = [4 . , . , .
0] are considered. The initial design forthis problem is constructed using cuboids as void inclusions. The right side of Fig. 10 shows the hole seedingarrangement for the finest mesh used in this example.The optimization problem seeks to minimize the strain energy, Ψ, subject to a mass, M , constraint; andis formulated as follows: min s z ( s , u ) = w Ψ( s , u ) / Ψ + w P P er ( s ) + w P Reg ( s ) s.t. : g ( s ) = M ( s ) / M T − γ m ≤ S
20 is enforced using the continuation strategy described in Table2 to avoid premature removal of mass resulting in suboptimal designs. This constraint is reduced every 12design iterations in a total of 8 steps (i.e., for D it ≤ M T . The optimization problem starts feasible for all mesh sizes, i.e., the initialhole seeding satisfies the initial mass constraint. The weights of the objective are w i =[ 0.85, 0.05, 0.10].Table 2: Continuation parameters of beam problem for γ m = 0 .
20 and three mesh sizes.Parameter D it = 1 D it = 13 D it = 25 D it = 37 D it = 49 D it = 61 D it = 73 D it = 85 γ m | h =4 . γ m | h =2 . γ m | h =1 . The evolution of the strain energy, Ψ, and mass constraint, g , for the results using the finest mesh (i.e., h = 1 .
0) are shown in Fig. 11. The sudden violations of the constraint, and to a lesser extent in the objective,are consequence of the mass constraint being updated at every continuation step. Overall, a smooth behavioris observed in all cases. The strain energies of the optimized designs are summarized in Table 3 and showthat the relatively small differences in performance are reduced as mesh is refined. In this example, the rulesused to resolve ambiguities have a negligible effect on both the smoothness of the evolution of the design, aswell as the performance of the optimized designs.Table 3: Final strain energy of beam problem using three different mesh sizes.Mesh size L Ω S L Ω V L L L A + L A − G Ω S G Ω V G A + G A − h = 4 9 . . . . . . . . . . h = 2 9 . . . . . . . . . . h = 1 9 . . . . . . . . . .
1% in all cases. The plot on theright of Fig. 13 shows the ratio of the volume of ATs for which a solid phase was assigned with respect tothe total solid volume, V AT Ω S /V Ω S . As expected, global rules G Ω S and G Ω V provide upper and lower limitsfor the cases considered; see Section 3.2.1. Furthermore, it can be seen that rules L , L , L , and G tendto distribute the phase assignment more evenly between the void and solid subdomains. In all cases, thevolume of the ATs is, at all times, considerably smaller than the uniquely defined domain. To demonstrate that the avoidance of ambiguities along the shear webs are not specific for the particularmass constraint studied above, lower target mass constraints, i.e., γ m = [0 . , . G Ω S and G Ω V and a mesh size of h = 1 . γ m = [0 . , . D it = 1 D it = 13 D it = 25 D it = 37 D it = 49 D it = 61 D it = 73 D it = 85 γ m | h =1 . γ m | h =1 . G Ω S and G Ω V phase assignment rules, for decreasingmass constraint from left to right and highlighting ATs in red. It can be seen that the reduction in massrequirements is satisfied by generating new holes rather than decreasing the thickness of the vertical thinwalls. A small influence on the topology consistency rule employed is observed on the shape and size of theseholes; see onsets of right column of Fig. 14. In this example, all ATs are internal, i.e., the configuration shownin Fig. 5(b). Hence, determine ambiguities as solid or void in the optimized designs can only generate shapechanges while keeping the same topology. However, optimized designs can also include BATs resulting fromhexahedra with multiple intersections. This is shown in Example 3 using a different objective formulation. In the second problem, we investigate the evolution of a design when the initial topology is significantlyaffected by the rule chosen to resolve ambiguities. The design domain of size 80 × ×
15 is subject to theboundary conditions specified on the left side of Fig. 15. All four top corners of the domain are clampedin all three directions within the volumes highlighted in Fig. 15. A prescribed traction of T X = − . × X − X and X − X planes at the center of the domain are considered. The domainis discretized using a mesh with an element edge length of h = 0 . G Ω S and G Ω V rules, and for γ m = [0 . , . , . X − X plane symmetry
15 is enforced through a continuation scheme. The continuation step size is 12, andthe maximum mass allowed decreases by setting γ m = [0 . , . , . , . , . , . , . , . , . G Ω S and G Ω V since these rules provide the greatest difference in theresulting geometries. To better visualize this, the additional volume defined as void when using the G Ω V ruleis highlighted in red on the right side of Fig. 16. The void subdomain is shown in this case. It can be seenthat rule G Ω V generates connections between the void spheres. These connections reduce the solid volumeby ≈ G Ω S and G Ω V rules on the top and bottom halves, respectively. Both designs converge to four legs spreading fromthe top corner clamped subdomains to the region at the bottom where the distributed load is applied.Different arrangements of truss-like connections between the legs and the center of the domain on the topare observed. Overall, the optimized structures differ only locally, and similar performances are achieveddespite starting from significantly different initial objectives. Hence, in this example, the optimizationframework compensates for the phase assignment rule and develop similar discretized geometries.The evolution of the V AT /V Ω S ratio is plotted in Fig. 18. Note the drastic reduction in the total volumeof ATs at the initial stages of the optimization process. The ATs in the optimized designs represent approxi-mately 1% in the extreme cases examined. Similar to the previous example, the majority of these ambiguitiesare originated from intersected hexahedra with a single interface. However, unlike in the previous example,here the ATs on the surface produce more pronounced differences in smoothness.Note that in problems with a strict mass requirement such as the one shown here, a large number ofholes may need to be seeded relatively close to each other since they have to cover a large portion of thedesign domain. In such scenarios, rules that favor void phase assignment, i.e. L Ω V and G Ω V , can make theinitial design unfeasible if the relative difference in volume of ATs is significant. For this reason, using suchrules can result in premature merging of holes and convergence to sub-optimal designs. In the final example, all local and global rules described in this paper are studied with a realistic engineeringdesign problem. This problem is characterized by a non-standard 3D geometry of the design domain, andcomplexities associated with multiple objective components and a stress constraint. The objective of thisoptimization problem consists of finding a structure that supports a payload given a set of supports andbolts for attaching the structure to the payload.The problem design and non-design subdomains are illustrated in Fig. 19. The design subdomain (Ω dI )is colored in light gray, while the non-design subdomains for the payload box (Ω pII ), the supports (Ω sII ),and bolts (Ω bII ) are colored in dark gray. A 3D view together with side and bottom views are provided fora better visualization of the constrained subdomains. A uniform pressure load acts on the top surface ofthe payload box. The entire structure is subject to a body force in X direction, representing an equivalentshock loading. The material properties and load conditions for this static problem can be found in Table 5.The goal of this optimization problem is to minimize the mass, M , of the structure that supports theFigure 16: Difference in the initial topology of the structural suspender problem (void phase) as a result ofthe phase assignment rule used. 21igure 17: Optimized designs of suspender problem using the G global rules. G Ω V
01, and w P τ = 1 . × . The objective scaling parameter is set to z sc = 10 . g , with γ m = 0 .
30) imposes an upper limit on the mass, in addition to consideringthe mass in the objective. The stress constraint, g , controls the stress field such that is does not exceed σ maxV M through the penalty: P τ = (cid:90) Ω d Υ dV. (26)This penalty term measures the volume where the von Mises stress exceeds the stress limit. To this end, theintegrand, Υ, is increased as the stress exceeds the stress limit. To ensure differentiability at Υ = σ maxV M , theintegrand is formulated as follows:Υ = (cid:40)(cid:2) (( τ − σ maxV M ) /σ maxV M ) + ξ (cid:3) / − ξ Υ , ∀ τ − σ maxV M > , , ∀ τ − σ maxV M ≤ , (27)23here the parameter ξ Υ > τ > σ maxV M , and is set to 0.1. The constant σ maxV M is given by the yield stress of Ti-6Al-4V, reduced by a safety factor of 2.0; see Table 5. Note that,although initially inactive, this stress constraint regulates the mass removal effect of the objective at thefinal stages of the optimization process.The design domain, supports and bolts are immersed into a hexahedral computational domain of size15 . × . × . cm . The computational domain is discretized by a hierarchically refined tensor mesh.This mesh is locally refined within the design domain to reduce computational cost while keeping a sufficientresolution in relevant regions. A 3D view and two plane views at the center in the X − X and X − X planes are shown on the left side and center of Fig. 20, respectively. The initial hole seeding pattern is shownon the right side of Fig. 20. Note that, unlike in Example 2, in this setup no connections are created betweenvoid inclusions with any of the rules employed.Figure 20: Bracket problem analysis setup: (left) 3D view of locally refined mesh used together with twoplane views; and (right) initial seeding in design domain.The mass and strain energy of the optimized designs are summarized in Table 6. The maximum differencein the mass of the optimized designs across all design rules is less than 0.3%. The differences in strain energyare negligible as well. Similar to what was observed in the previous examples, the performance is notcompromised by any of the topology consistency rules tested.Table 6: Final mass and strain energy of bracket example using all local and global rules. Z component L Ω S L Ω V L G L A + L A − G Ω S G Ω V G A + G A − M ( × − ) 4 . . . . . . . . . . × ) 5 . . . . . . . . . . L Ω V and G A − rules, additional thinner barscomplement the main frame connecting the box and the front supports.Overall, although the outer frame of the bracket is similar in all cases, secondary features with variousshapes and sizes appear in various regions of the design domain depending on the rule used to resolve ATs. Ageneral trend of these features in terms of location cannot be inferred. Fluctuations in the stress distributionscan also be observed. Regions of maximum stress are concentrated on thinner features in rules L A + , G Ω V ,and G A + . Nevertheless, the objective convergence behavior is comparable in all cases, and the differences inperformance between results using any of the rules is negligible.24igure 21: Bracket problem optimized designs using all local and global rules for resolving topology ambi-guities. 25 Summary and Future Work
The effect of ambiguities in the phase assignment of tetrahedralized intersected hexahedra was analyzed inthe context of level set topology optimization. To ease the analysis of the underlying issue, these ambigu-ities were classified as internal ambiguous tetrahedra (i.e., all faces of an ambiguous tetrahedron inscribedwithin the hexahedron), and boundary ambiguous tetrahedra (i.e., one of the ambiguous tetrahedron facescoinciding with a face of the containing hexahedron). Inconsistencies encountered for specific intersectionpatterns were resolved through rules applied locally (i.e., individually for each intersected hexahedron), orglobally (i.e., across connected intersected hexahedra). Using local rules, consistent phase assignment ofboundary ambiguous tetrahedra across intersected elements was enforced via the asymptotic decider. Inaddition, internal ambiguous tetrahedra were resolved using criteria based on predetermined user-defined orgeometrical indicators. Conversely, in the global rules described in this paper, a single criterion is appliedto both boundary and internal ambiguous tetrahedra.The sensitivity of the design process and performance of optimized designs on the aforementioned phaseassignment rules were investigated using an explicit level set topology optimization framework. Althoughshape and topological changes, such as small dents in the interface or holes in regions with thin features, canappear, they were not observed in the optimized designs. The minor influence of the ambiguous tetrahedraresolution schemes on the optimized design is likely due to the ability of the optimization process to implicitlyaccount for the phase assignment rule. Distinct geometrical configurations developed initially lead to verysimilar discretized geometries in optimized designs. Overall, the numerical experiments showed that theoptimized designs achieved similar performances for various optimization problem formulations regardlessof the rule employed. None of the rules examined appeared to be superior to the others, nor any particularrule seems to adversely affect the optimization process. As long as a rule is used consistently throughoutthe optimization process, a smooth evolution of the objective can be expected regardless of the rule used.Nevertheless, we recommend that studies on 3D level set topology using immersed boundary methods withexplicit interface representations, such as the XFEM, explicitly state the criteria used for resolving topologicalambiguities. This enhances the reproducibility of results and, thus, an informed analysis on features formedin the optimized designs.Additional rules that promote smoothness in the interface while not restricting the evolution of designs,such as curvature measures, should be studied. Furthermore, the numerical efficiency and robustness of therules provided should be tested for solid-solid domains, and various physical models.
Compliance with Ethical Standards
The authors declare that they have no conflict of interest. Upon request, the authors will provide the fullset of input parameters for each topology optimization problem presented in the paper. Any optimizationframework with an implementation of the rules described here to determine a phase for ambiguous tetrahedrashould be able to reproduce the results shown in the numerical examples section.
Acknowledgments
All authors acknowledge support of the Air Force Office of Scientific Research (grant FA9550-16-1-0169).The first author acknowledges support from NSF (Grant 1463287). The second author acknowledges supportfrom the Defense Advanced Research Projects Agency (DARPA) under the TRADES program (AgreementHR0011-17-2-0022). The opinions and conclusions presented in this paper are those of the authors and donot necessarily reflect the views of the sponsoring organizations.