Abstract

Optimal geometries extracted from traditional element-based topology optimization outcomes usually have zigzag boundaries, leading to being difficult to fabricate. In this study, a fairly accurate and efficient topology description function method (TDFM) for topology optimization of linear elastic structures is developed. By employing the modified sigmoid function, a simple yet efficient strategy is presented to tackle the computational difficulties because of the nonsmoothness of Heaviside function in topology optimization problem. The optimization problem is to minimize the structural compliance, with highest stiffness, while satisfying the volume constraint. The design problem is solved by a Sequential Linear Programming method. Convergent, crisp, and smooth final layouts are obtained, which can be fabricated without postprocessing, demonstrated by a series of numerical examples. Further, the proposed method has a rather higher accuracy and efficiency compared with traditional TDFM, when the classical topology optimization methods, such as bidirectional evolutionary structural optimization (BESO) and solid isotropic material with penalization (SIMP) method, are taken as benchmark.

1. Introduction

In the field of structural optimization, topology optimization has gained more and more attention. Opposed to the sole variation of the design region boundary in classical shape optimization [1, 2], it aims at searching for the optimal designs by determining the best locations and geometries of cavities in the design domain; in other words, topology optimization allows producing new holes within the design domain.

Topology optimization of continuum structures has its roots in the landmark work of Bendsoe and Kikuchi [3]; it has been developed extensively in the past three decades. Among the various numerical methods, most of them are based on finite element analysis (FEA), where the design domain is discretized by a fine mesh of elements. These topology optimization approaches may be classified into two categories, that is, element-based topology optimization method ([46]) and node-based topology optimization method ([7]).

In the element-based topology optimization method, the optimization procedure is to find the topology of a structure by determining for every point in the design domain whether there should be material (solid element) or not (void element). This type of method including homogenization method ([3]), solid isotropic material with penalization (SIMP) method ([4, 8]), evolutionary structural optimization (ESO) method ([6]), and its improved algorithm C bidirectional evolutionary structural optimization (BESO) method ([912]) has tackled various problems, including uncertain design ([1317]), dynamic problems ([1820]), and designing metamaterials ([21]). The aforementioned element-based topology optimization methods have applications successfully in a number of fields. However, they still have some undesirable shortcomings. From the manufacturing viewpoint, final layouts produced from element-based topology optimization method always have zigzag boundaries even the design domain is discretized into fine elements, so these solutions require some postprocessing to smooth the boundaries before they can be manufactured. This may limit their application.

Node-based topology optimization method is developed to overcome the disadvantages of element-based topology optimization, particularly for the zigzag boundaries issues. It can produce solutions with crisp and smooth edges that require little postprocessing effort to interpret results before they can be manufactured, which is due to the fact that it employs an implicit function to describe the shape and topology of the structure. The implicit function is called topology description function (TDF) in node-based topology optimization method, which is first proposed by Sethian and Wiegmann [22]. In their article, they first apply an explicit jump immersed interface method to compute the structural stresses for a given design domain discretized by finite differences, and then the level set method, which represents the design structure through an embedded topology implicit function (TDF), is used to alter the structural shape, with velocities depending on the stresses in the current design; criteria are also provided for advancing the shape and introducing holes in an appropriate direction. Since this pioneering work, various node-based topology optimization methods have been reported. Reiter and Keulen [23] developed a topology optimization method, where Radial Basis Functions are used to generate the TDF. Genetic algorithm is applied but with large computational costs. Later, some researchers put the idea of level-set-based structural optimization into the shape-sensitivity-based optimization framework. However, the aforementioned approaches are associated with the solution of a Hamilton-Jacobi-type equation, which may considerably slow down the speed of optimization convergence. New strategies have been implemented in the topology optimization method to solve the Hamilton-Jacobi-type equation [2427], aiming at improving the computational efficiency and enhancing the numerical stability [2830]. Recently, Guos group made great progress in parametric level set topology optimization methods, which significantly reduced the number of the design variables and in turn tremendously increase the computational efficiency [3134].

In TDF-based topology optimization method, as discussed in the previous paragraph, Heaviside function is usually employed to define structural boundaries [35]. Using these approaches to design structures may run into computational difficulties because of the nonsmoothness of the Heaviside function. The exact Heaviside function cannot be directly differentiated and is, therefore, often replaced with a smooth approximation of the Heaviside function. To the authors’ knowledge, piecewise polynomials and trigonometric functions are the two main kinds of functions to smooth the Heaviside function. Although, in the traditional TDFM, using the existent smooth functions can successfully obtain the convergent optimum solutions, they may have low accuracy due to the relatively poor coincidence degree of the existent smooth functions and the Heaviside function.

Therefore, in this article, we try to use a more accurate function to smooth the Heaviside function in order to improve the accuracy of the TDFM. Modified sigmoid function, thanks to its inherent properties, is employed to smooth the Heaviside function. The TDFM proposed in this paper is activated by the work of Liu et al. [35], but our method has higher accuracy and efficiency compared with their work. Numerical tests demonstrate that, by employing the modified sigmoid function to approximate the Heaviside function, the convergent optimum solutions with crisp and smooth edges are obtained. The solutions obtained require little postprocessing effort to interpret results before they can be manufactured. This characteristic is preferable for topology optimization to the practical engineering application.

The outline of this article is as follows. In Section 2, the strategy of smoothing the Heaviside function is proposed. In Section 3, a topology-description-function-based method for structural optimization using modified sigmoid function is proposed. Remarks on the characteristics of the proposed method are discussed in Section 4 and several numerical examples are shown in Section 5. Finally, the important conclusions are presented in Section 6.

2. Smooth the Heaviside Function

Topology optimization approaches such as the one based on implicit functions to define structural boundaries often employ Heaviside function to vividly represent the absence or presence of the material in the design domain. However, nonsmoothness is the inherent property of Heaviside function, and the topology optimization problem will hardly be solved if we use the Heaviside function directly. It is because discrete values, used as the intermediate design variables during the topology optimization process, are obtained from the calculation using Heaviside function, and it is difficult to solve the design problem utilizing the existing famous numerical computation methods such as Sequential Linear Programming (SLP) [36], Optimality Criteria (OC) method [37], or the method of moving asymptotes (MMA) [38]. Thus, Heaviside function should be smoothed in order to make the design problem easily solved and further improve the practicability of TDFM.

As discussed above, there are two main approaches to smooth the Heaviside function, and both use smooth piecewise function to approximate the Heaviside function. However, the relatively poor coincidence degree of the existent smooth functions and the Heaviside function may lead to low accuracy of the traditional TDFM. This section therefore introduces an alternative approach, where modified sigmoid function is employed to approximate the Heaviside function, aiming at improving the accuracy of the TDFM.

2.1. Sigmoid Function

Sigmoid function [39] is a continuous nonlinear activation function. The origin of the name, sigmoid, is from the fact that the function is -shaped. This function is called logistic function by the statisticians. Using as input, as output, and which is a positive number as a contrast factor term, the sigmoid function can be expressed as

The sigmoid function has the following features: it is a smooth continuous function with its input variable varying from minus infinity to positive infinity and its output lies in the range from 0 to 1. Figure 1 depicts the plot of function for various value of . Note that Figure 1 only shows part of the plot within the input ranges from −1 to 1 for the sake of easy description.

Based on the characteristics of the sigmoid function discussed above, we try to deeply analyze the relationship between the sigmoid function and the Heaviside function and then use modified sigmoid function to smooth the Heaviside function. They will be discussed in the next subsection.

2.2. Smooth the Heaviside Function

The Heaviside function is a noncontinuous function whose value is 0 for negative argument and 1 for positive argument. The function is originally used in the mathematics of control theory to represent a signal that switches on at a specified time and stays switched on indefinitely. It is the integral of the Dirac delta function, which is sometimes written as

Alternatively, the Heaviside function can be defined aswhere is a real number.

Now we try to find the relationship between the sigmoid function and the Heaviside function. It can be clearly seen from Figure 2 that the sigmoid function gradually tends to overlap the Heaviside function with the increase of the value of α. Note that the red lines represent the plot of the Heaviside function. Naturally, we can use the sigmoid function to approximate and smooth the Heaviside function when is a large number, i.e., .

In this paper, modified sigmoid function is employed to smooth the Heaviside function in order to make the topology optimization problem easy to be solved and with high accuracy when topology-description-function-based method is utilized. This will be proposed in the next section.

3. Topology-Description-Function-Based Method

In this section, the basic idea of TDFM is first proposed as a prerequisite. Then modified sigmoid function is introduced for smoothing the Heaviside function in order to make the design problem easy to be solved and improve the accuracy of the TDFM.

3.1. Method Description

TDFM method is based on the main idea that it is to map a function on the reference design domain into geometry. The function is determined by the design parameters, which can be done in several ways. It is worth pointing out that the function is implicit, and it results in one of the most important merits of TDFM, which will be highlighted in Section 4. Without loss of generality, one of the simplest ways to do this is proposed here for two-dimensional space, which can easily be extended to three-dimensional space.

A function in the initial design domain is used to describe the structural topology as follows. Figure 3 illustrates the concept expressed in (4) geometrically.

Then, we combine the properties of the Heaviside function and TDF to describe geometries in a discrete fashion, that is, without intermediate densities (BESO method). Thus, herein, the element densities can be defined in the following:where is considered as the design variable and is the Heaviside function. Equation (5) illustrates the relationship between element densities and the Heaviside function.

In this article, the topology optimization problem, with the objective of minimizing the compliance, aims at finding the structure with highest stiffness, and the constraint on the structural volume is formulated as follows:where denotes the structural compliance. and represent the applied load and displacement vectors, respectively. is the one-dimensional Dirac function. refers to the prescribed objective volume.

3.2. Discrete Processing of the Optimization Problem

TDFM is a kind of topology optimization approach based on the finite element method (FEM). Thus, FEM is employed to discretize the design domain for the optimization problem in (6) before the design problem can be solved. The initial design domain is first divided into homogeneous linear rectangular elements. And then the displacement vectors and are interpolated by a bilinear function, respectively.

where is the total number of nodes in the finite element mesh. denotes the interpolation shape function, which can be obtained in unit isoparametric coordinates as

where and describe horizontal and vertical coordinates, respectively. refers to the node number in the local coordinate system, which is shown in Figure 4.

After the finite element discretization, the optimization problem in (6) can be reformulated as follows:

The displacement interpolation in (10) can be obtained by solving the following equation:where represents the trial function in a confined space. denotes the projection of the surface force on the confined space.

In order to prevent the TDF from presenting fast shock distributions during the process of discretization, the elasticity tensor is regarded as constant in each element and its value is related to the values of on its four nodes. In addition, the idea of SIMP is also used to penalize the element densities:where denotes penalization exponent (typically ).

Then, the stiffness matrix of the element is given in the following:where and refer to the strain matrix and elasticity tensor, respectively. is the elemental stiffness matrix of the solid element.

Although the discretization of the optimization problem has been done successfully, as we discussed in Section 2, it is still hard to solve the problem using mathematical methods due to the nonsmoothness of the Heaviside function. Therefore, the Heaviside function should be regularized in order to make the design problem easily solved. Now, sigmoid function discussed in Section 2 will be modified and later employed to smooth the Heaviside function, which will be detailed below.

3.3. Dealing with Nonsmoothness

The Heaviside function, with its nonsmooth characteristics, is approximated by the modified sigmoid function as follows:where both the smoothing parameters and are small positive numbers with the purpose of avoiding singularity; is the bandwidth with a value of 0.3 and in this paper. The superscript represents the modification for the Heaviside function.

Consequently, the corresponding Dirac function is obtained in the following equation:

Now, the Heaviside function is smoothed successfully, and the design problem in (10) can be reformulated as

In the next subsection, we try to solve the design problem using gradient-based optimization method, where sensitivity analysis is a necessity.

3.4. Sensitivity Analysis

The optimization problem in (16) could be solved using several different approaches such as SLP, OC method, or MMA. For simplicity, an SLP is used in this paper. In mathematical programming framework, sensitivity information is necessary for numerical solution. The sensitivities of the objective and constraint functions with respect to the design variables can be, respectively, obtained by (17) and (18) as follows:

where represents the total number of elements.

3.5. Convergence Criterion

There should be a convergence criterion to decide when the final layout is obtained and stop the optimization procedure. The cycle of finite element analysis and element removal continues until the limited volume is reached and the following convergence criterion (defined in terms of the change in the objective function) is satisfied.

where is the current iteration number, is an allowable convergence tolerance whose value is 0.001 in this article, and is an integer number. Normally, is set to be 5 which implies that the change in the mean compliance over the last 10 iterations is acceptably small.

4. Remarks on the Merits of Proposed TDFM

In this section, the virtues of the present TDFM are highlighted. As mentioned above, the present TDFM is based on implicit functions to define structural boundaries, so crisp and smooth edges are obtained, requiring little postprocessing effort to interpret results before they can be manufactured. This is one of the most important characters of the proposed method compared to the classical element-based topology optimization methods, such SIMP and ESO/BESO, and this merit may enormously improve the practicability of TDFM. In addition, the proposed TDFM has higher accuracy compared to the traditional TDFM due to the fact that the smooth function employed in this article approximating the Heaviside function has higher accuracy. All of this can be verified by several examples in Section 5.

As discussed above, the proposed TDFM in this article has higher accuracy compared to the traditional TDFM due to the fact that the smooth function employed approximating the Heaviside function has higher accuracy when the value of is large enough. From Figure 2, it seems that, by giving the parameter a larger value, the sigmoid function may simulate the Heaviside function well. However, it also increases the nonlinearity of the problem and may possibly have remarkable and negative influence on the solution of the problem. And the solution of the problem may more easily be trapped into local optima. Thus, in this paper, a relatively large number of will be used.

In traditional TDFM, piecewise polynomials and trigonometric functions are the two main kinds of functions to smooth the Heaviside function. Equations (20) and (21) are two well-known smooth functions to approximate the Heaviside function; they are, respectively, named smooth function 1 and smooth function 2 for the sake of easy description.where and represent smooth function 1 and smooth function 2, respectively. and have the same meanings as those in (14).

Given that a large may lead to many computational problems, the value of is selected as 50 (a relatively large number) in the numerical examples in Section 5. And this selection will be tested using numerical examples. Apparently, it can be seen, from Figure 5, that the sigmoid function has a relatively high accuracy to approximate the Heaviside function compared with smooth function 1 and smooth function 2; it can also visually describe the reason why the proposed TDFM has higher accuracy compared to the traditional TDFM. Note that Figure 5 only shows the images when the range of is from −0.3 to 0.3.

5. Numerical Examples

In this section, we present three numerical examples in order to illustrate the applicability of the proposed method as well as the benefit from using the modified sigmoid function to smooth the Heaviside function. For test purpose, the material properties, Young’s modulus of 1 GPa, and Poisson’s ratio of 0.3 are assumed, and the volume constraint is 50% of the design domain in all the numerical examples. For comparison purpose, traditional TDFM, SIMP, and BESO are also applied to solve the design problems. In SIMP, rmin = 3.0, while in BESO, ER = 2%, ARmax = 50%, and rmin = 3.0. Note that traditional TDFMs using and smoothing the Heaviside function are, respectively, named Traditional 1 and Traditional 2 for convenient description in the numerical examples. The central processing unit (CPU) times presented are based on optimization runs carried out on Lenovo Personal Computer equipped with 1 Intel Core i3-2120 processor with a clock speed of 3.3 GHz, 2.99 GB RAM, and Windows XP 32-bit operating system. Notice that the classical SIMP and BESO methods are taken as benchmarks in the following examples, and the computational time by using these two approaches will not be considered as their efficiency is varied with different unfixed parameters.

5.1. Topology Optimization of a Short Cantilever

A classical topology optimization problem of a corner loaded short cantilever has been extensively studied by many researchers. The cantilever, fixed on the upper surface, is 60 mm long and 40 mm deep with a point load of 3 kN acting at the light corner as illustrated in Figure 6. The design domain is modeled with four-node plane stress elements.

The optimal material distributions obtained from the proposed TDFM, traditional TDFM 1, traditional TDFM 2, BESO method, and SIMP method are compared in Figures 7(a)7(e). Both final layouts show the distribution of smoothed nodal densities. In the figures, the black areas denote the solid regions in the final design outcome; the white areas indicate the void regions. It can be seen that, as discussed in Section 4, the optimum structure produced by the proposed TDFM has crisp and smooth edges when the modified sigmoid function is implemented, which has the benefit that it requires little postprocessing effort to interpret the final layout before it can be manufactured. This merit of the proposed TDFM can be clearly manifested compared to those produced by BESO method and SIMP method. It is worth pointing out that crisp and smooth edges can be also obtained by BESO method and SIMP method when small FEM grid is used but the computational efficiency will be greatly reduced. In addition, the final result produced using the proposed TDFM is almost the same as those obtained by the other topology optimization methods, which demonstrates the validity of the proposed method.

Figure 8 shows the evolution history of the structural compliance of solution produced by the proposed TDFM. It can be seen that the compliance gradually decreases and then is convergent to an almost constant value within less than 21 iterations. For comparison purpose, Table 1 gives the final compliance of the different topologies produced by different topology optimization methods. It can be found that the structural compliance of the proposed TDFM solution is in excellent agreement with that of classical SIMP solution and slightly differs from that of classical BESO solution, whereas compliances of traditional TDFM 1 and traditional TDFM 2 are larger than classical BESO and SIMP solutions. Therefore, now, it is safe to say that the proposed method is much more accurate than the traditional TDFMs. In addition, CPU time is also compared among the proposed TDFM, traditional TDFM 1, and traditional TDFM 2 in Table 1. The proposed approach is found to be up 55% and 41% more efficient in terms of computational time when compared to traditional TDFM 1 and traditional TDFM 2, respectively, due to the implantation of the sigmoid function in TDFM.

5.2. Topology Optimization of Simply Supported Beam Structure

In this numerical example, the well-known simply supported beam problem is considered. The support is located on both left-hand side and right-hand side of the initial design domain (see Figure 9). The design domain, 100 mm in length and 40 mm in height, is divided into four-node plane stress elements. The load is applied at the bottom middle of the design domain with a value of 2 kN.

In Figure 10(a), the topology obtained by the proposed TDFM is presented. It is compared with those produced by the other methods discussed in the previous numerical example (see Figures 10(b)10(e)). As can be seen in Figure 10, the final layout produced by the proposed TDFM has crisp and smooth edges when the modified sigmoid function is implemented compared with those obtained from the classical BESO method and SIMP method. This merit may considerately improve the practical engineering application of the proposed TDFM, as it requires little postprocessing effort to interpret the optimum before it can be manufactured. It should be noted that the topology proposed from the proposed TDFM is in excellent agreement with that of classical SIMP solution and is different from that of classical BESO solution in which two more thin bars appear. In addition, the proposed TDFM solution slightly differs from those of traditional TDFM solutions as the latter ones have one or two additive very small holes. It is because different smoothed Heaviside functions are utilized.

The evolution history of the structural compliance of the proposed TDFM solution is depicted in Figure 11. It indicates that it has good convergence such that the compliance gradually decreases and then is convergent to an almost constant value within less than 24 iterations.

The final structural compliances and CUP time of different topologies are compared in Table 2. Analyzing Table 2, we can see that the compliance of the proposed TDFM solution is almost the same as that of the classical SIMP solution, whereas the compliances of the traditional TDFM solutions are larger than that of the classical SIMP solution. It can be concluded that the proposed TDFM has much higher accuracy compared with traditional TDFMs. It can be also found from Table 2 that the proposed approach yields 56% and 65% saving, respectively, in terms of computational time compared to traditional TDFM 1 and traditional TDFM 2. Thus, in this specific example, we can conclude that the proposed TDFM is very efficient. The designer can benefit greatly from this merit when designing actual engineering structures, which is always very time-consuming.

5.3. Topology Optimization of a Beam Structure Subjected to Two Applied Loads

In the last example, topology optimization of a beam structure with consideration of two concentrated forces applied at its upper surface is studied. As shown in Figure 12, the support is located on both left-hand side and right-hand side of the initial design domain. Using an element size of 1 mm, the design domain is discretized by four-node plane stress elements.

The compared results of various topology optimization methods are illustrated in Figures 13(a)13(e). It follows from the comparison among Figures 13(a), 13(d), and 13(e) that the proposed TDFM solution is very similar to the classical BESO and SIMP solutions. So, the effectiveness of the proposed approach can be verified in this numerical example. It should be noticed that the optimum topologies generated from the traditional TDFMs have two relatively big holes in the left bottom and right bottom, respectively. We may say that the traditional TDFMs are less effective than our proposed approach in this specific example. But all the traditional TDFMs and the proposed TDFM can produce crisp and smooth final layouts, which makes them more advantageous than BESO and SIMP in this point.

Figure 14 demonstrates the good convergence of the proposed approach. Evolution history of the structural compliance presented in Figure 14 shows that the compliance gradually decreases and then is convergent to an almost constant value within less than 20 iterations.

In order to investigate the accuracy and efficiency of the proposed approach, Table 3 summarizes the final structural compliances of various topologies and the CPU time to produce these final layouts. Similar to the previous two numerical examples, our approach proposed in this article has higher accuracy compared to traditional TDFMs.

6. Concluding Remarks

In this paper, a relatively accurate TDFM is reported for linear elastic problems. Modified sigmoid function instead of the usual piecewise polynomials and trigonometric functions is employed to smooth the Heaviside function. In doing this, the inherent nonsmoothness property of the Heaviside function is addressed. The computational difficulties due to the nonsmoothness of the Heaviside function in the topology optimization problem are successfully tackled.

The major contribution of this article comes from the improved efficiency of the proposed algorithm when measured in terms of computational effort involved in the FEA runs required to solve the optimum solution compared with the traditional TDFMs. In addition, the accuracy is also greatly enhanced when the classical SIMP and BESO methods are taken as benchmarks. Moreover, numerical examples demonstrate the effectiveness of the proposed approach and also highlight the merit such that the proposed TDFM solutions have crisp and smooth edges compared with the classical element-based topology optimization methods, that is, BESO and SIMP.

Data Availability

No data were used to support this study.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This work was supported by NSFC Project (11671342 and 11672104), Hunan Education Department Key Project (17A210), and Hunan Province Natural Science Fund (2018JJ2374).