This paper presents an extended sequential element rejection and admission (SERA) topology optimization method with a region partitioning strategy. Based on the partitioning of a design domain into solid regions and weak regions, the proposed optimization method sequentially implements finite element analysis (FEA) in these regions. After standard FEA in the solid regions, the boundary displacement of the weak regions is constrained using the numerical solution of the solid regions as Dirichlet boundary conditions. This treatment can alleviate the negative effect of the material interpolation model of the topology optimization method in the weak regions, such as the condition number of the structural global stiffness matrix. For optimization, in which the forward problem requires nonlinear structural analysis, a linear solver can be applied in weak regions to avoid numerical singularities caused by the over-deformed mesh. To enhance the robustness of the proposed method, the nonmanifold point and island are identified and handled separately. The performance of the proposed method is verified by three 2D minimum compliance examples.
Topology optimization aims to obtain the optimal material distribution with prescribed constraints in a fixed domain. Since the pioneering work of Bendsoe et al. [
In an early discrete variable method, the evolutionary structural optimization (ESO) method as proposed by Xie et al. [
In hard-kill type methods, the ability to change topology is usually related to the material introduction strategies [
Introducing a weak material has been proven to be very successful in solving static small deformation problems and has become a standard model in topology optimization. To ensure smooth access to the optimal solution, the physical property of the weak material must be controlled within a reasonable range. If the value is too large, the stiffness provided by it cannot be ignored, which may lead to an overly rigid performance evaluation. If too small, the condition number of the structural stiffness matrix becomes worse, and large numerical errors may occur in the results. Therefore, the condition number is usually selected as 10−3. In addition, in some specific optimization problems, a weak material may cause numerical instability or even incorrect results. For example, when a structure has a large deformation, excessive mesh deformation may occur in weak regions, which leads to nonconvergence of the Newton–Raphson iteration [
The hard-kill type methods and some other methods that remove the elements in the holes can avoid the problems caused by weak materials [
From the perspective of model consistency, the physical model and finite element (FE) model are consistent in the hard-kill type methods, but the FE model of the soft-kill type methods is an approximate model from the physical model. If the inconsistencies between the two models are eliminated, it is possible to solve the above problems in the framework of soft-kill type methods.
In SERA, the concepts of real material and virtual material are introduced to distinguish the elements (it should be noted that the real material and virtual material here correspond to the solid and weak material mentioned above. To avoid confusion, only the concepts of solid and weak materials are retained in the next content. The sensitivities of solid elements and weak material elements (weak elements) are ranked, and two thresholds are selected to determine the elements to be removed or to be added. This method indirectly divides the design domain into solid regions and weak regions but does not further identify them. The design domain is still regarded as a whole for FEA in this method. Tong et al. [
Based on the characteristics of the above methods, a topology optimization method with a region partitioning strategy is proposed in this paper. The method is implemented based on SERA. The FEA is carried out regionally. To avoid the influence of weak material on the result of solid regions, FEA of solid regions is carried out first, where the weak elements are eliminated temporarily in this step. When the FEA of the solid regions is completed, FEA of the weak regions is performed with the displacement solution of solid regions as the Dirichlet boundary conditions (
In the implementation, two special structural features, nonmanifold points and islands, may arise in solid regions after region segmentation. These two structural features might cause instability in the numerical solution of the displacement via FEA. We introduce concepts and algorithms in graphics to recognize and solve them: the nonmanifold points are identified by an undirected graph and treated by the minimum volume filling (MVF) algorithm and minimum sensitivity filling (MSF) algorithm [
The rest of the paper is arranged as follows: the topology optimization algorithm with a region partition strategy is detailed in
In SERA, solid regions and weak regions can be distinguished by densities of elements. To avoid the negative effects of weak regions on solid regions, the FEAs for the two types of regions are carried out in sequence. The optimization formula is as follows:
The only difference between the
In the proposed method, the FEA of solid regions is prioritized. Weak elements are temporarily removed in this step. When FEA of the solid regions is completed, FEA of the weak regions is performed with the displacement solution of solid regions as the Dirichlet boundary conditions. The FEA procedure is shown in
Since the element densities can only be one or
In regard to the finite deformation problem, the nonlinear FE formula must be considered:
The most frequently applied scheme for solving nonlinear systems is the Newton–Raphson algorithm. It is based on a Taylor series development of
After the displacement vector
In the FE model of weak regions, Dirichlet boundary conditions are imposed on the interface to maintain the continuity of the displacement field. The implementation process is as follows: first, the components of the displacement vector
The global stiffness matrix, displacement vector and load vector of the weak regions are also divided into two parts according to the noninterface nodes and the interface nodes. The updated FE formula is:
Since
When the load-independent design problem is considered, the weak regions usually have no external load so that
The displacement vector
In the case of small deformation problems, the sensitivity number
When addressing small deformation problems, according to Huang et al. [
It should be noted that the sensitivity numbers in
Two sensitivity thresholds
The above strategies can effectively avoid the adverse effects caused by weak elements and maintain the integrity of sensitivities. However, in the execution process, structural instability may occur due to some special structural features. To ensure the robustness of the algorithm, these structural features must be solved.
After regional partitioning, there are two structural features that may cause numerical instability: nonmanifold points and islands. They must be addressed to prevent system ill conditioning. In this section, we introduce some strategies to identify and address them.
Hinged elements often appear in solid regions when material is almost removed (
It is necessary to introduce the concept of a 2D manifold to describe nonmanifold points. A surface is a 2D manifold if it has the following characteristics: (1) Every point on the surface of an object has a sufficiently small neighborhood isomorphic to the disk on the plane; (2) If the surface is unclosed (it has boundaries), every point on the boundaries must have a sufficiently small neighborhood isomorphic to the half disk in the plane. Since it is impossible to describe each point on a surface by computer, the conditions identifying that a surface is a 2D manifold are represented by vertices, edges and faces of the discrete grid: (1) each edge can only be shared by two faces; (2) the edges and faces adjacent to each vertex can form one ring around it. For example, there are three points in
As the design domain is discretized by a structured rectangular mesh, the nonmanifold point is presented in only one form: adjacent elements of the vertex contain two diagonally distributed solid elements and two diagonally distributed weak elements. The procedure of identifying nonmanifold points is as follows and the flow chart is as shown in Create the vertex set Calculate the vertex set Create the edge set Calculate the vertex set Calculate the nonmanifold point set
The strategy to accommodate nonmanifold point is converting the weak elements adjacent to the nonmanifold point into solid elements, as shown in
Structures consisting of nonmanifold points can be divided into two categories: hinged structures and checkerboard patterns. Since the adjacent weak elements have no intersection in hinged structures, converting any of the surrounding weak elements can address the problem (
To reduce the disturbance to the volume and the objective, we introduced two algorithms to address nonmanifold points: the minimum volume filling algorithm (MVF) and the minimum sensitivity filling algorithm (MSF). The MVF requires filling a minimum number of solid elements to reject all nonmanifold points, and the MSF requires that the elements that are filled have the least impact on the objective.
The concept of an undirected graph is introduced to implement the above two algorithms. The graph
The undirected graph
Minimum Volume Filling (MVF)
The MVF can be transformed into minimum point set coverage, which means that each edge in the undirected graph must be associated with at least one of its vertices. In graph theory, the minimum point set covering is equal to the maximum matching so that the Hungarian algorithm [ Minimum Sensitivity Filling (MSF)
The MSF can be transformed into minimum point set coverage with weight in the undirected graph. In the algorithm, the sensitivity values of the weak elements adjacent to each nonmanifold point are added to the vertices of the undirected graph as weights. The minimum point set coverage with weight requires that every edge in the weighted undirected graph is associated with at least one vertex and the sum of the weights of these vertices is minimum.
As shown in
In
In SERA, the initial structure is usually made entirely of solid elements. In the early optimization process, solid elements are gradually transformed into weak elements. This may lead to a situation in which some solid elements are completely separated from the main structure, forming an island (
To accommodate islands, we need to identify them first. The recognition of the nonmanifold point mentioned above is the recognition of elements, while an island is a region composed of an uncertain number of elements. Although the domain is divided in SERA, the numbers of solid regions and weak regions are still unknown. We introduce the fire-burning method [
To realize the fire-burning method, the adjacent relation of the elements must be calculated first. The concept of a graph is still used to explain here. For the graph
The graph
For a structured grid, the coding rules of elements and their nodes are generally fixed and regular. For example, the simplest coding method is to encode from top to bottom and from left to right. We can evaluate whether the two elements are adjacent by evaluating whether there is a common node between the two elements, and the adjacency matrix of the elements can be established by this information.
The process of the fire-burning method is as follows:
Determine an element as the initial burning point; Search the adjacent untraversed elements with the same density as the current burning points; these elements are regarded as the new burning points. Iterate Step 2 until there are no more qualified elements around the current burning points. All the elements in the current region can be found by the iteration. Determine an untraversed element as a new burning point, and repeat Steps 2 and 3. Once there are no untraversed elements in the design domain, the process ends, and all the regions have been found.
Since all the regions have been found, we can determine whether a region is an island by evaluating whether there are constrained points or load points in the region. Considering little contribution to the objective function and adverse effects on the numerical stability, we group islands into weak regions for FEA.
In this section, the concept of a conditional number is introduced to measure the ill conditioning of a system, and three examples are performed to verify the effectiveness of the proposed method. The first one is the short beam example; the second one is the C-shape plate example [
The condition number of a matrix is a well-known measure of ill conditioning. For an
When considering the static small deformation problem, the design domain is discretized by the FE mesh, and the governing equation can be expressed by a linear equation:
In SERA, the element stiffness matrix is expressed as follows:
A simple example can illustrate this.
Since the condition number of the global stiffness matrix is a direct parameter to measure numerical stability, it can be seen from the above that two cases can lead to ill-conditioned systems:
An extremely small ratio of the minimum and maximum densities. Local unstable structures.
The method proposed in
The load and support conditions of the short beam example are shown in
In SERA, the optimal structure can usually be obtained faster by deleting materials than by adding materials. Therefore, the initial structures in the subsequent examples are all composed of solid elements. In the first half of the optimization process, the solid elements are gradually rejected until the volume constraint is met, and then the same amounts of solid elements and weak elements are increased or decreased until convergence.
In this example, four values, 10−9, 10−18, 10−27, and 10−36, are selected as the design values
To verify the computational efficiency of the proposed method, the computational time is counted. The original SERA consumed 9.8 and 18.4 s for computing the two examples in
In this example, we also verify the stability of the proposed method under different mesh densities. Three mesh densities are considered: 1) coarse, consisting of 100 × 50 rectangular elements; 2) medium, consisting of 200 × 100 rectangular elements; 3) fine, consisting of 300 × 150 rectangular elements. The loading position is shown in
The optimal structural topologies are almost the same for the three mesh densities. Only the boundary smoothness and the size of local holes are slightly different. Therefore, the robustness of the proposed method is not affected by the mesh density. However, there is another factor to consider as the mesh density increases: in the proposed method, the nonmanifold point processing is time-consuming. Once the elements increase, the nonmanifold points may increase, which will reduce computational efficiency. To adapt to large-scale computation, we consider how to improve the efficiency of nonmanifold point processing in future work.
To verify the effectiveness of the proposed method for solving excessive mesh deformation under finite deformation, we considered an example proposed by Yoon et al. [
For better comparison, we analyzed the structure for the following three options: (1) geometric linearity or geometric nonlinearity; (2) inclusion or exclusion of weak regions; and (3) with or without a region partition strategy. The results are shown in
Geometric nonlinearity | Weak regions | Region partition | Convergence | |
---|---|---|---|---|
× | × | - | √ | |
√ | × | - | √ | |
× | √ | × | √ | |
√ | √ | × | × | |
√ | √ | √ | √ |
When no weak region is included (
Now, we apply the topology optimization method with a region partitioning strategy to this problem. The design domain and boundary conditions are still consistent with the model in
As shown in
The design domain and boundary condition of the cantilever beam are shown in
In this example, we verify the difference in optimization results under geometric linear and geometric nonlinear assumptions based on a region partitioning strategy. The effects of different loads on the topology optimization results are examined simultaneously. The results are shown in
The third line in
Numerical instability problems caused by a material interpolation model in the weak region are considered in this paper. A topology optimization model with geometrical region partitioning and the corresponding FEA procedure is proposed. Numerical strategies to identify nonmanifold points and island structures are discussed, and numerical instability during the optimization procedure can be alleviated by properly addressing the nonmanifold points. Numerical examples illustrate that the proposed method can solve the ill-conditioned matrix problem in FEA and nonconvergence phenomena caused by excessive mesh deformation in finite deformation. The region partition strategy proposed in this paper can be implemented from SERA to other topology optimization methods, such as the SIMP method and boundary evolution type method.