Topology Optimization Using Parabolic Aggregation Function with Independent-Continuous-Mapping Method

Abstract

This paper concentrates on finding the optimal distribution for continuum structure such that the structural weight with stress constraints is minimized where the physical design domain is discretized by finite elements. A novel Independent-Continuous-Mapping (ICM) method is proposed to convert equivalently the binary design variables which is used to indicate material or void in the various elements to independent continuous design variables. Moreover, three smooth mappings about weight, stiffness, and stress of the structural elements are introduced to formulate the objective function based on the so-called concepts of polish function and weighting filter function. A new general continuous approach for topology optimization is given which can eliminate the stress singularity phenomena more efficiently than the traditional 𝜀-relaxation method, and an alternative strain energy method for the stress constraints is proposed to overcome the difficulty in stress sensitivity analyses. Mathematically, by means of a generalized aggregation KS-like function defined as the parabolic aggregation function, a topology optimization model is formulated with the weight objective and single parabolic global strain energy constraints. The numerical examples demonstrate that the proposed methods effectively remove the stress concentrations and generate black-and-white designs for practically sized problems.

1. Introduction

Topology optimization is one of the most challenging research topics but also a research focus in the current structural optimization design field, which has been rapidly developed with a lot of fruitful research work in the past two decades, and the development of topology optimization can be found in the monographs [1–3] and references therein. However, the representative methods to structural optimization design for variable topologies are homogenization method, SIMP method, ESO method and Level-set method, and so forth.

Homogenization method is a very important continuous nonlinear optimization method for topology optimization in structural design [4–13]. This method assigns the orientation and geometry size of the microstructure as the design variables, in which the microstructure is deleted according to the changes in the unicellular size, and leads to anisotropic materials with the middle-sized microstructures. By transforming the difficult structural topology optimization model into a relatively easier sizing optimization model based on the microstructural parameters of composite materials, the approach of homogenization is capable of finding the optimum topologies of the structures. Since this formulation has two kinds of design variables to each element, the sensitivity analysis and the solution process become complicated. However, homogenization method often produces designs with infinitesimal pores in the materials that make the structure unnecessary strong and nonmanufacturing.

After the original homogenization-based method, the solid isotropic material with penalization (SIMP) approach has been developed and well established as its conceptual and practical simplicity [6, 14–16] with applications into many areas such as fluids, multiphysics and multidisciplinary applications, and much more [1, 17–20]. As a variant of the homogenization method, the SIMP method introduces artificially the material with variable pseudodensity between 0-1 by allowing the material to take intermediate densities. However, a general nonlinear interpolation optimization model for the material properties is formulated, and the density of each finite element is chosen as the design variable in which its relationship with the Young’s modulus is expressed by an empirical formula. It was proved that the physical microstructures for SIMP scheme would exist for the satisfaction of some simple conditions on the penalty exponent [16]. Some additional schemes proposed by Sigmund [21] and Petersson and Sigmund [22] should be considered so as to eliminate numerical instabilities such as the checkerboards and mesh-dependency.

Evolutionary structural optimization method referred to as ESO has been presented by Xie and Steven [23, 24] which is successfully applied to optimum material distribution problems for continuum structures including stress considerations, frequency optimization, and stiffness constraints so as to gradually remove material and achieve an optimal design. However, the ESO method is based on the evolutionary strategy in the optimization process by gradually removing ineffective or inefficient material from the ground structure to get the topology plots through each iteration. In order to improve and perfect the traditional ESO method, its later development bidirectional ESO (BESO) is proposed by Querin et al. [25], Yang et al. [26], and Rozvany [19], which allows elements to be added in the locations next to those elements with highest sensitivity values as well as to be eliminated in the region with lowest sensitivity values. For the hard-kill and soft-kill methods in the ESO method, Huang and Xie [27] proposed a new soft-kill BESO method based on the material interpolation scheme with penalization which is the general case of the original hard-kill BESO method [28], when the exponent is not tended to infinity.

Another class of approaches is essentially based on the stress distribution on the current boundaries, which is called as the level set method proposed by Sethian and Wiegmann [29] so as to capture design boundary movements on a fixed Eulerian mesh. In the standard level set method, the structural boundary is described as the zero level set of a higher-dimensional scalar function by an implicit level model (Hamilton-Jacobi PDE) and the rigorous shape derivative [13, 30–32]. Luo et al. [33, 34] proposed the alternative full parameterization level set method for structural shape and topology optimization problems using radial basis functions with compact support, in which the structural boundary is still represented implicitly as the zero level set of the level set function by taking into account the merit of the implicit-free boundary representation so as to bridge the standard level set method with many optimization algorithms. However, the level set method has a good flexibility to describe the topology of complex structures and boundary change.

The structural topology optimization may come down to solve the optimal distribution problems in which subsystems of the overall structure that should be filled with material and should be void, where the skeleton structures including truss structures and frame structures are optimized by the elimination or restoration of some elements and the continuum structures are optimized by making some subdomain into empty or restoration. In practice, the design domain is usually discretized into finite elements in order to obtain a finite dimensional problem with zero-one binary variables, indicating presence or absence of material in the various finite elements, but this leads to integer nonlinear programs which are very hard to solve a global optimum even for relatively small problems [35–37]. Since the large-scale integer programming for solving the current problems has no efficient algorithms, people usually use the low-level design variables such as size or materials instead of the traditional zero-one binary topological design variables. Homogenization method attaches the topological variables on the size of material microstructure, and SIMP method attaches the topological variables on the Young’s modulus.

However, each above-mentioned scheme has its pros and cons; and thus, there is still room left for improvement. Hence, the ideal scheme for topology optimization still has to be invented.

The paper is organized as follows: in Section 2, the independent-continuous-mapping (ICM) method for topology optimization problems is described. The ICM method abstracts the topology variables from low-level design variables such as the cross-section, thickness, and density so as to be independent of the traditional physical quantities and uses continuous variables in [0,1] instead of 0 and 1 binary discrete variables to avoid the difficulties of discrete optimization. Moreover, three smooth mappings about weight, stiffness, and stress of the structural elements are introduced to formulate the objective function based on the so-called concepts of polish function and weighting filter function to complete the transformation of topological variables from discrete to continuous and to the discrete by the approximation of continuous variables with discrete topology variables. In Section 3, an alternative strain energy method for the topology optimization problem with local stress constraints is presented based on the von Mises criteria in the theory of elastic failure and the theory of strain energy under tridirectional stress state to overcome the difficulty in stress sensitivity analyses. In Section 4, a new generalized aggregation KS-like function defined as the parabolic aggregation function is derived, which is different from the two common global functions such as 𝑝-norm and KS function, and an optimization model is formulated with the weight objective and a single global constraints. In Section 5, the serial Lagrange multiplier method (SLMM) is applied to deal with the structural reanalysis of the nonlinear topology optimization problems, and implemental issues are discussed. In Section 6, three widely studied numerical examples are presented to demonstrate the effectiveness of the proposed topology optimization method, where conclusions are drawn in the end of the paper.

2. Independent-Continuous-Mapping (ICM) Method

In this section, we give a brief introduction to the independent-continuous-mapping structural optimization method which is based on the idea presented by Sui [38], Sui and Yang [39], Sui [40] and Sui et al. [41, 42].

2.1. Domain Setting and the Independent Topological Design Variables

Let us consider a fixed design domain Ω in ℝ³ in topology optimization of load carrying continuum structures. The infinite dimensional topology optimization problem is to find the subdomain Ω_solid  filled with material (or the subdomain Ω_void  occupied by void) such that a structure of optimum performance is obtained while constraints are satisfied, where 𝜒(𝑥) is an indicator function that defines the distribution of material [1, 6] 𝜒(𝑥)=1,if𝑥∈Ω_solid,0,if𝑥∈Ω_void,Ω_solid=(𝑥),Ω𝑥∣𝜒=1_void=𝑥∣𝜒(𝑥)=0. (1)(2)

Suppose that the fixed design domain Ω is discretized into a ground structure with 𝑛 finite elements, we can define the binary design variables 𝑡_𝑖∈{0,1}  as 𝑡_𝑖=1 when the 𝑖th element is filled with material and 𝑡_𝑖=0 when the 𝑖th element is void such that the problem is then to decide the binary design variable vector 𝑡=(𝑡₁,𝑡₂,…,𝑡_𝑛)^𝑇∈{0,1}^𝑛.

Now, we consider the general binary discrete optimization problem with some property constraints as follows.

Find topological variables  𝐭∈[0,1]^𝑛, minimizeΦ=Φ(𝑡)=^𝑛_𝑖=1𝜙_𝑖,subjectto𝑔_𝑗𝜏₁,…,𝜏_𝑛−𝑔_𝑗≤0,𝑗=1,…,𝐽,0≤𝑡_𝑖≤1,𝑖=1,…,𝑛, (3) where Φ=Φ(𝑡)=∑^𝑛_𝑖=1𝜙_𝑖  is the objective function, 𝜏_𝑖 is the physical or geometric parameter of the 𝑖th element used in the finite analysis, 𝑔_𝑗(𝜏₁,…,𝜏_𝑛) and 𝑔_𝑗 are the 𝑗th constraint response and its upper bound, respectively, 𝐽 is the total number of constraints, 𝑡_𝑖 is the 𝑖th independent topological design variable, 𝑡=(𝑡₁,…,𝑡_𝑛)^𝑇 is the independent topological design vector, and the total number of design variables is 𝑛. In fact, 𝜏_𝑖 may be regarded as 𝑡_𝑖  function with  𝜏_𝑖=0  as 𝑡_𝑖=0 and 𝜏_𝑖≠0 as 𝑡_𝑖≠0.

Due to the large dimensionality of the design space, this discrete problem is often relaxed to nonlinear optimization problems with continuous design variables instead of the binary design variables. Several methods have been proposed for topology optimization based on the physical or geometric values such as the variable element thickness [43, 44], cross-sectional area [45], or pseudodensity as the design variables, that is, 𝑡_𝑖≈𝜌_𝑖,0≤𝜌_𝑖𝑡≤1,_𝑖≈𝑑_𝑖,0≤𝑑_𝑖𝑡≤1,_𝑖≈𝐴_𝑖,0≤𝐴_𝑖≤1, (4) where 𝜌_𝑖, 𝑑_𝑖, and 𝐴_𝑖 are the 𝑖th element normalized density, thickness, and area, respectively, and 𝑝 is a penalization parameter in the SIMP method [16, 46].

To obtain the independent topological design variables by abstracting from the above physical parameters and geometrical parameters, the binary design variables 𝑡_𝑖∈{0,1} in (1) can be relaxed to a new step-up function and be defined as (see Figure 1) 𝑡_𝑖𝜏=𝑠^𝑖𝜏⁰_𝑖=𝜏1,^𝑖𝜏⁰_𝑖∈(],𝜏0,10,^𝑖𝜏⁰_𝑖=0, (5) where 𝜏_𝑖 is the physical or geometric parameter of the 𝑖th element used in the finite analysis, 𝜏⁰_𝑖 is the original and true physical or geometric parameter of the material used of the 𝑖th element.

Figure 1 

The step-up function and traditional topology design variable.

As an appropriate relaxation, 𝑡_𝑖=𝑠(𝜏_𝑖/𝜏⁰_𝑖) is allowed to take all possible values between 0 and 1, and it is clear that 𝑡_𝑖=𝑠(𝜏_𝑖/𝜏⁰_𝑖) is an ideal independent topological design variable which has clear physical meaning. However, the step-up function (5) is discontinuous in the point 𝜏_𝑖/𝜏⁰_𝑖=0, but it is reasonable. In the conventional structural topology optimization, no matter how small, as long as the physical or geometric parameter of the 𝑖th element 𝜏_𝑖 is not equal to 0, then 𝑡_𝑖=𝑠(𝜏_𝑖/𝜏⁰_𝑖) is equal to 1, and when 𝜏_𝑖 is equal to 0, then the topological design variable suddenly becomes 0. In fact, the step-up function (5) describes the above-mentioned discontinuous relation between 𝑡_𝑖 and 𝜏_𝑖/𝜏⁰_𝑖.

2.2. The Polish Function and the Continuous Approximation of the Step-Up Function

Since the optimization problem (3) includes finite discontinuities that express the existence or nonexistence of elements, it is, therefore, not easily treated numerically and not solved by classical optimization algorithms based on continuous variables. In order to overcome this difficulty, we introduce the polish function 𝑝(𝜏_𝑖/𝜏⁰_𝑖) to approximate the step-up function (5) such that 𝑡_𝑖=𝑠(𝜏_𝑖/𝜏⁰_𝑖) becomes a continuous and differentiable function (see Figure 2); that is, 𝑡_𝑖𝜏=𝑠^𝑖𝜏⁰_𝑖𝜏≈𝑝^𝑖𝜏⁰_𝑖, (6) where the following properties are satisfied.(1)Continuity and differentiability: 𝑝(𝜏_𝑖/𝜏⁰_𝑖) is a continuous and differentiable function in [0,1].(2)Strict increasing monotonicity: 𝑝(𝜏_𝑖/𝜏⁰_𝑖) is a strictly monotonically increasing function.(3)Convexity: 𝑝(𝜏_𝑖/𝜏⁰_𝑖) is a convex function in [0,1].(4)Boundedness: 𝑝(𝜏_𝑖/𝜏⁰_𝑖)∈[0,1], 𝑝(0)=0, and 𝑝(1)=1.(5)Approximation: ∫¹₀|1−𝑝(𝑥)|𝑑𝑥<1/2 or 𝑥≤𝑝(𝑥)≤1 in [0,1].

Figure 2 

The polish function and independent continuous topological design variable.

Based on the concepts of the step-up function and polish function, 𝑡_𝑖=𝑝(𝜏_𝑖/𝜏⁰_𝑖) can take values in the entire range [0,1] instead of the discrete set {0,1} and make the functional relationship 𝑡_𝑖=𝑠(𝜏_𝑖/𝜏⁰_𝑖) between the topology design variables 𝑡_𝑖, and the various physical or geometric parameters 𝜏_𝑖 become continuous and differentiable from discontinuous and nondifferentiable. This “discrete-relaxation-continuous” mapping process is called as polishing approximation mapping. In this paper, we take 𝑡_𝑖=𝑝(𝜏_𝑖/𝜏⁰_𝑖) as the independent continuous topological design variables. Due to this extension of the design space, the existence of solutions is achieved (see Figure 7).

From the mechanical point of view, when the independent continuous topological design variables 𝑡_𝑖=𝑝(𝜏_𝑖/𝜏⁰_𝑖) take intermediate values between 0-1, 𝑡_𝑖=𝑝(𝜏_𝑖/𝜏⁰_𝑖) not only reflect the proximity to corresponding elements “have” and “no”, but also identify the size of the various physical values 𝜏_𝑖 in the elements.

2.3. The Filter Mappings and the Weighting Filter Functions

To formulate the explicit continuous nonlinear topology optimization model for the optimization problem (3) and obtain the final topological graphs, we define the inverse mapping of the polishing approximation mapping as the filter mapping, which is denoted as following: 𝜏𝑓∶^𝑖𝜏⁰_𝑖⟼𝑝⁻¹𝑡_𝑖; (7) that is, 𝜏_𝑖𝜏⁰_𝑖=𝑓𝑡_𝑖=𝑝⁻¹𝑡_𝑖or𝜏_𝑖=𝑓𝑡_𝑖𝜏⁰_𝑖=𝑝⁻¹𝑡_𝑖𝜏⁰_𝑖. (8)

Moreover, we define the inverse function of the polish function 𝑓(𝑡_𝑖)=𝜏_𝑖/𝜏⁰_𝑖=𝑝⁻¹(𝑡_𝑖) as the weighting filter function, where 𝑡_𝑖∈[0,1].

It is well known that the polish function and the weighting filter function may have a variety of options. Sui [38] suggested three kinds of simple filter functions as following (see Figure 3):

Figure 3 

Three kinds of simple filter functions.

(a) power function 𝑥=𝑡^𝛼,𝛼=3,10.32; (9)

(b) composite exponential function 𝑒𝑥=^𝑡/𝛾−1𝑒^1/𝛾−1,𝛾=0.0621; (10)

(c) modified sigmoid function 1𝑥=22ln1+𝑡/𝛽1−𝑡/𝛽,𝛽=1. (11)

In fact, the inverse functions of the above three kinds of simple weighting filter functions (9), (10), and (11) satisfy continuity and differentiability, strict increasing monotonicity, convexity, boundedness, and approximation. However, the weighting filter function may be considered as the penalty function in the SIMP approach.

In the ICM method, the weighting filter function 𝑓(𝑡_𝑖)=𝜏_𝑖/𝜏⁰_𝑖=𝑝⁻¹(𝑡_𝑖) plays a very important role as follows.

(a) Let the inverse mapping of the step-up function in (5) denote as 𝜏_𝑖/𝜏⁰_𝑖=𝑠⁻¹(𝑡_𝑖), which is a 0-1 binary multivalued function that cannot indicate the size of the physical or geometric parameter 𝜏_𝑖 when the independent topological design variables 𝑡_𝑖=𝑠(𝜏_𝑖/𝜏⁰_𝑖) take intermediate values between 0-1. To overcome this difficulty, we introduce the concept of weighting filter function 𝑓(𝑡_𝑖)=𝜏_𝑖/𝜏⁰_𝑖=𝑝⁻¹(𝑡_𝑖) which continuously approximates the inverse mapping of the step-up function 𝜏_𝑖/𝜏⁰_𝑖=𝑠⁻¹(𝑡_𝑖), and it points out the proximity degree that 𝜏_𝑖 closes to 𝜏⁰_𝑖 and 0.

(b) According to the ICM method, if 𝜏_𝑖=𝑤_𝑖 or 𝐾_𝑖, then the variable weight and stiffness matrix of the 𝑖th element, 𝑤_𝑖 and 𝐾_𝑖 are, respectively, recognized by the weighting filter functions as 𝑤_𝑖=𝑓_𝑤𝑡_𝑖𝑤⁰_𝑖,𝐾_𝑖=𝑓_𝐾𝑡_𝑖𝐾⁰_𝑖, (12)(13) where 𝑤⁰_𝑖 and 𝐾⁰_𝑖 are the original weight and stiffness matrix of the 𝑖th element for structural material, respectively.

Generally, when 𝜏_𝑖 in (8) represents the different variable physical or geometric parameter, the weighting filter functions could have different forms; that is, 𝜎_𝑖=𝑓_𝜎𝑡_𝑖𝜎⁰_𝑖,𝑚_𝑖=𝑓_𝑚𝑡_𝑖𝑚⁰_𝑖,𝑣_𝑖=𝑓_𝑣𝑡_𝑖𝑣⁰_𝑖,𝐴_𝑖=𝑓_𝐴𝑡_𝑖𝐴⁰_𝑖,𝑑_𝑖=𝑓_𝑑𝑡_𝑖𝑑⁰_𝑖, (14)(15)(16)(17)(18) where  𝜎_𝑖, 𝑚_𝑖, 𝑣_𝑖, 𝐴_𝑖, and 𝑑_𝑖 are, respectively, the variable allowable stress, mass, volume, area, and thickness of the 𝑖th element for structural material, and 𝜎⁰_𝑖, 𝑚⁰_𝑖, 𝑣⁰_𝑖, 𝐴⁰_𝑖, and 𝑑⁰_𝑖 are the original allowable stress, mass, volume, area, and thickness of the 𝑖th element for structural material, respectively.

Moreover, we can give the solid isotropic material with penalization (SIMP) approach a new explanation from the view of the ICM method, that is, 𝜌_𝑖=𝑓_𝜌𝑡_𝑖𝜌⁰_𝑖,𝐸_𝑖=𝑓_𝐸𝑡_𝑖𝐸⁰_𝑖, (19)(20) where 𝜌_𝑖 and 𝐸_𝑖 are the variable density and Young’s modulus of the 𝑖th element for structural material and 𝜌⁰_𝑖  and 𝐸⁰_𝑖 are the original density and Young’s modulus of the 𝑖th element, respectively.

(c) Based on the concept of the weighting filter function, the original objective and constraint functions are replaced by relatively simple, approximating functions which are fundamentally different from the penalization scheme and interpolation scheme in the SIMP method. Firstly, the ICM method takes the independent continuous topological design variable 𝑡_𝑖=𝑝(𝜏_𝑖/𝜏⁰_𝑖) instead of the relative density 𝜌_𝑖 in the SIMP method, which changes conventional continuous design variable taking the various physical or geometric parameters. Secondly, from the ICM method point of view, as long as the physical or geometric parameters change in the structural topology optimization process, the related weighting filter function (or interpolation scheme in terms of the SIMP method) should be considered as not limiting to the stiffness matrix. However, the weighting filter functions play the important roles as they eliminate grey scale transitions between solid and void regions.

(d) The weighting filter function can formulate the explicit continuous nonlinear topology optimization model and obtain the final topological graphs.

Let 𝜙_𝑖 be the objective function of the 𝑖th element used in the finite analysis and 𝜙⁰_𝑖 the original objective function of the 𝑖th element for structural material, then we have that 𝜙_𝑖=𝑓_𝜙𝑡_𝑖𝜙⁰_𝑖 (21) and the topology optimization model (3) could be approximated by the new continuous nonlinear topology optimization model; as follows.

Find topological variables  𝐭∈[0,1]^𝑛minimizeΦ=^𝑛_𝑖=1𝜙⁰_𝑖𝑓_𝜙𝑡_𝑖subjectto𝑔_𝑗𝜏⁰₁𝑓_𝜏𝑡₁,…,𝜏⁰_𝑛𝑓_𝜏𝑡_𝑛−𝑔_𝑗≤0,𝑗=1,…,𝐽,0≤𝑡_𝑖≤1,𝑖=1,…,𝑛. (22)

(e) The problems considered in this paper involve (von Mises) stresses, which can exhibit the so-called “singularity phenomena”, that is, due to the discontinuous nature of the stress when the topology design variables 𝑡_𝑖 tend to zero in (14). This stress singularity phenomena has been dealt with by the smooth envelope functions method [47–49] or the𝜀-relaxation approach [36, 50–56].

According to the ICM method, we can also introduce the concept of the variable allowable stress 𝜎_𝑖 of the 𝑖th element, which is recognized by the allowable stress weighting filter functions as (14) 𝜎_𝑖=𝑓_𝜎𝑡_𝑖𝜎⁰_𝑖 (23) such that 𝜎_𝑖𝑗≤𝑓_𝜎𝑡_𝑖𝜎⁰_𝑖, (24) where 𝜎⁰_𝑖 is the original allowable stress of the 𝑖th element for structural material and 𝜎_𝑖𝑗 is the von Mises stress for the 𝑖th element under the 𝑗th load case. That is, the inequality (24) guarantees that 𝜎_𝑖𝑗→0, when 𝑡_𝑖 tends to zero without violating the stress constraint.

In fact, (23) and (24) eliminate the above-mentioned singularity phenomena, since a stress will be exactly zero in a void element, and they also remove the need for taking the complex 𝜀-value in the 𝜀-relaxation approach.

By exploiting the modeling of a simple microstructure [1, 51], they state that the similar variable allowable stress as stress threshold 𝜌^𝑝𝜎_𝑙, where they apply a stress constraint for the SIMP model, that is, expressed as a constraint of the von Mises equivalent stress 𝜎_VM: 𝜎_VM≤𝜌^𝑝𝜎_𝑙,if𝜌>0. (25)

In the following, we prove that (23) and (24) can satisfy the𝜀-relaxation conditions.

Without loss of generality, assuming the stress constraints in the continuum structure design [36, 51–56] as 𝜎_𝑖𝑗≤𝜎_𝑖, (26) then we have that 𝜎_𝑖𝑗≤𝜎_𝑖=𝑓_𝜎𝑡_𝑖𝜎⁰_𝑖,𝜎_𝑖𝑗<𝜀+𝑓_𝜎𝑡_𝑖𝜎⁰_𝑖, (27)(28) where 𝜀>0, 𝜀→0, and 𝜀 is an arbitrary small quantity in which its value determines the extent of relaxation and the inequality.

Moreover, the inequality (24) and (28) can imply 𝜎_𝑖𝑗−𝜎⁰_𝑖𝑓_𝜎𝑡_𝑖=𝜎_𝑖𝑗𝑓_𝜎𝑡_𝑖−𝑓_𝜎𝑡_𝑖𝜎⁰_𝑖≤𝜎_𝑖𝑗−𝑓_𝜎𝑡_𝑖𝜎⁰_𝑖<𝜀; (29) that is, 𝜎_𝑖𝑗−𝜎⁰_𝑖𝑓_𝜎𝑡_𝑖<𝜀. (30)

It is well known that the inequality (30) is exactly the ordinarily formula in the 𝜀-relaxation approach.

In conclusion, the key idea of the ICM method is the independent continuous topological design variable, the polish and weighting filter mapping, and the simple discrete or inverse process.

3. An Alternative Strain Energy Method for the Topology Optimization Problem with Stress Constraints

In this section, an alternative strain energy method for the topology optimization design problem with stress constraints is proposed. The method is implemented in the topology optimization procedure to overcome the difficulty in stress sensitivity analyses.

3.1. Minimum Weight with Stress Constraints in Topology Optimization

For the general explicit continuous nonlinear topology optimization model (22), let 𝜎_𝑖𝑗 be the von Mises stress for the 𝑖th element under the 𝑗th load case, 𝛾_mat=const the material density, and Ω the original reference design domain. However, the integral on the volume of Ω_solid filled with material of the material density 𝛾_mat=const can be written as an integral on the total volume Ω by multiplying 𝑑Ω by 𝑓_𝑤(𝑡_𝑖), and the objective function Φ=∑^𝑛_𝑖=1𝜙⁰_𝑖𝑓_Φ(𝑡_𝑖) may be defined as 𝑊(𝑡)=_Ωsolid𝛾_mat𝑑Ω=_Ω𝑓_𝑤(𝑡)𝛾_mat=𝑑Ω_𝑛𝑖=1𝑤_𝑖=^𝑛_𝑖=1𝑓_𝑤𝑡_𝑖𝑤⁰_𝑖. (31) Then, the original topology optimization problem with stress constraints can be formulated as follows.

Find topological variables  𝐭∈[0,1]^𝑛; minimize𝑊(𝑡)=^𝑛_𝑖=1𝑓_𝑤𝑡_𝑖𝑤⁰_𝑖,subjectto𝜎_𝑖𝑗≤𝜎_𝑖=𝑓_𝜎𝑡_𝑖𝜎⁰_𝑖,0<𝑡_min≤𝑡_𝑖≤1,𝑖=1,…,𝑛,𝑗=1,…,𝐽, (32) where 𝑡_𝑖 is the 𝑖th independent continuous topological design variable following (8), 𝑡=(𝑡₁,…,𝑡_𝑛) is the independent continuous topological design vector, the total number of design variables is 𝑛,  and the total number of external loads is 𝐽. As it is well known, a lower limit 𝑡_min of the the 𝑖th independent continuous topological design variable is imposed to avoid the singularity of the global stiffness matrix which is given as 𝑡_min=0.001.

However, Sui et al. [42] and Tie et al. [57] study topology optimization of continuum structure with globalization of stress constraints using an efficient strain energy method.

In fact, most of papers referred to what was previously mentioned and the monographs ([1–3] and references therein) have been focused on global objectives, such as compliance, displacement, and frequency. It is fewer about the strain energy method in the structural topology optimization.

Since the essence of topology optimization lies in searching for the optimal path of transferring loads, the stress of the structural elements in the structural topology optimization has to be considered to obtain a reliable design. Clearly the stress of the structural elements is a local quantity, there are three difficulties for solving the local problem as follows: (a) the dimension number of the design variable is too extensive, and it is too difficult to be solved. However, every structural element in the structural topology optimization is just with the corresponding one design variable. Even for one stress at a point in each element region, it is only restricted as a stress constraint of the region; the whole structure with 𝑛 elements will be associated with 𝑛×𝐽 stress when there is 𝐽 load cases. (b) Since the stress of the structural elements is highly nonlinear with respect to the design variables, the computational cost of stress sensitivity analysis is too expensive to be accepted. (c) The stress singularity phenomena must be also considered in the topology optimization problem mentioned above.

Thus it is necessary to introduce an alternative strain energy method to transform the local stress constraints into the strain energy constraints of the structural elements by means of the von Mises criteria in the theory of elastic failure.

3.2. The Derivation of an Alternative Strain Energy Method

It is well known that the von Mises stress at a point for deformation body is calculated by the formula as following 𝜎_VM=12𝜎₁−𝜎₂²+𝜎₂−𝜎₃²+𝜎₃−𝜎₁², (33) where 𝜎₁, 𝜎₂, and 𝜎₃ are the principal stresses.

Let 𝜎_VM≤𝜎, and looking back upon the shear strain energy density from the von Mises strength theory in (33), then we have that 𝑢^𝑓≤1+𝜈(3𝐸𝜎)²,𝑢^𝑓=1+𝜈𝜎6𝐸¹−𝜎₂²+𝜎₂−𝜎₃²+𝜎₃−𝜎₁², (34)(35) where 𝐸, 𝜈 are the Young’s modulus, Poisson’s ratio and 𝜎 is the allowable stress of the structural material, respectively.

Since the shear strain energy density is impossible to surpass the strain energy density, we have that 𝑢^𝑓=𝑢𝛾²𝑢,𝛾≥1,_𝑖𝑗=𝑒_𝑖𝑗𝑣_𝑖, (36) where 𝑢_𝑖𝑗, 𝑒_𝑖𝑗, and 𝑣_𝑖 are the strain energy density, the strain energy, and the volume of the 𝑖th element under the 𝑗th load case, respectively.

Note that the amplification factor 𝛾 can be computed in terms of so-called corrector terms and be given as 𝛾=1.2 in our work. Other values of 𝛾 could be considered, but they give similar results.

Let 𝐸_𝑖 be the variable Young’s modulus of the 𝑖th element for structural material, substituting (36) into (35), we have that 𝑒_𝑖𝑗≤1+𝜈3𝐸_𝑖𝛾𝜎_𝑖²𝑣_𝑖=𝑒_𝑖𝑗, (37) where we remark that 𝑒_𝑖𝑗 is defined as the variable allowable strain energy of the 𝑖th element which is unique under multiload case and may be denoted by 𝑒_𝑖=𝑒_𝑖𝑗.

It is evident that the structure is safe, and the stress constraints 𝜎_𝑖𝑗≤𝜎⁰_𝑖 are satisfied if the formulation (37) is established.

Moreover, substituting (14) and (20) into 𝑒_𝑖=𝑒_𝑖𝑗 in (37), we have that 𝑓_𝐸𝑡_𝑖𝑒_𝑖=𝑓²_𝜎𝑡_𝑖𝑒⁰_𝑖,𝑒⁰_𝑖=1+𝜈3𝐸⁰_𝑖𝛾𝜎⁰_𝑖²𝑣_𝑖, (38) where 𝑒⁰_𝑖 is defined as the original allowable strain energy of the 𝑖th element which is also unique under multiload case.

Thus, we can modify (32) to a new formulation according to (37) such that the strain energy constraintis 𝑒_𝑖𝑗≤𝑒_𝑖 in (37) instead of the stress constraints 𝜎_𝑖𝑗≤𝜎⁰_𝑖 in (32) as following.

Find topological variables𝐭∈[0,1]^𝑛; minimize𝑊(𝑡)=^𝑛_𝑖=1𝑓_𝑤𝑡_𝑖𝑤⁰_𝑖,subjectto𝑒_𝑖𝑗≤𝑒_𝑖,0<𝑡_min≤𝑡_𝑖≤1,𝑖∈𝐼_𝑛={1,…,𝑛},𝑗∈𝐽_𝑚={1,…,𝑚}. (39)

By the properties of the parabolic aggregation function (see the next section), we have that lim_𝑝→∞𝐺𝑒_𝑖𝑗,𝑝=max_𝑖,𝑗𝑒_𝑖𝑗,lim_𝑝→∞𝐺𝑒_𝑖,𝑝=max_𝑖𝑒_𝑖. (40)

Substituting (40) into (39), the optimization model (32) is equivalent with the following when the parameter 𝑝 is enough big.

Find topological variables  𝐭∈[0,1]^𝑛; minimize𝑊(𝑡)=^𝑛_𝑖=1𝑓_𝑤𝑡_𝑖𝑤⁰_𝑖,subjectto𝐺𝑒_𝑖𝑗,𝑝≤𝐺𝑒_𝑖,𝑝=𝐺,0<𝑡_min≤𝑡_𝑖≤1,𝑖∈𝐼_𝑛={1,…,𝑛},𝑗∈𝐽_𝑚={1,…,𝐽}. (41)

4. The Parabolic Aggregation Function and the Global Strain Energy Approach

In this section, we present a new strategy to deal with topology optimization of continuum in (39). However, the above strain energy-based approach produces highly nonlinear objective function and constraints so as to give rise to a complicated problem. That is, in finite-element-based topology optimization, the number of strain energy 𝑒_𝑖𝑗 is very large, and it is necessary to take into account that this approach may require large computing resources, and sensitivity computation is prohibitively costly.

A way out for overcoming this difficulty is to apply an aggregation function to deal with the large number strain energy which is frequently used in the topology optimization with stress constraints [52, 53, 56, 58–60]. Now we should explore a new aggregation function which differs from the common 𝑝-norm and KS function [61] and is very suitable for the topology optimization of continuum structures with stress or strain energy inequality constraints.

The proposed aggregation function that aggregates all the elemental strain energy into single global strain energy can be formulated as 𝐺𝑒_𝑖𝑗=1,𝑝𝑝ln^𝑛𝑚_𝑖=1𝑗=1𝑒^{𝑝2𝑒2𝑖𝑗1/2}, (42) where 𝐺(𝑒_𝑖𝑗,𝑝) is a KS-like global function and designated as the parabolic aggregation function ([62]; see Figure 4) at the same time we call KS function as the linear aggregation function. Indeed, the parameter 𝑝 does not depend on the strain energy, and the exponent of the parabolic aggregation function becomes dimensionless.

Figure 4 

The three-dimensional surface graphics of the parabolic aggregation function.

A simple mathematical proof for the global property of the parabolic aggregation function will be given here later. Moreover, in the context we denote that 𝐺=(1/𝑝)[ln(∑^𝑛_𝑖=1∑^𝐽_𝑗=1𝑒^{(𝑝𝑒0𝑖𝑓2𝜎(𝑡𝑖))2/𝑓𝐸(𝑡𝑖)})]^1/2 as the variable global allowable strain energy.

Theorem 1. Let 𝑣=𝑒^𝑢2𝑝2, if the following conditions are satisfied: (1)𝑢>0,   𝑝>0;(2)𝑓_𝑗(𝑥)>0,  𝑗=1,2,…,𝑛,𝑥∈𝑅^𝑚.And let 𝐺(𝑓_𝑖(𝑥),𝑝)=(1/𝑝)[ln(∑^𝑛_𝑖=1𝑒^{𝑝2𝑓2𝑖(𝑥)})]^1/2; then one has that
(a)lim_𝑝→+∞𝐺𝑓_𝑖(𝑥),𝑝=max_{1≤𝑖≤𝑛}𝑓_𝑖(𝑥); (43)
(b)min_{𝑥∈𝑅𝑚}max_{1≤𝑖≤𝑛}𝑓_𝑖(𝑥)=lim_𝑝→+∞min_{𝑥∈𝑅𝑚}𝐺𝑓_𝑖(𝑥),𝑝. (44)

Proof. (a) If 𝑣=𝑒^𝑢2𝑝2, then we have that 1𝑢=𝑝(ln𝑣)^1/2. (45) For 𝑓_𝑖(𝑥)>0,𝑝>0,𝑣𝑓_𝑖(𝑥),𝑝=𝑒^{𝑝2𝑓2𝑖(𝑥)}. (46) Let max_{1≤𝑖≤𝑛}𝑓_𝑖(𝑥)=𝑓(𝑥). Since 𝑢>0, 𝑝>0, then we have that 𝑣𝑓(𝑥),𝑝≤^𝑛_𝑖=1𝑣𝑓_𝑖(𝑥),𝑝≤𝑛𝑣𝑓(𝑥),𝑝; (47) that is, 𝑒^{𝑝2(𝑓(𝑥))2}≤^𝑛_𝑖=1𝑒^{𝑝2𝑓2𝑖(𝑥)}≤𝑛𝑒^{𝑝2(𝑓(𝑥))2},𝑓(𝑥)≤1𝑝ln^𝑛_𝑖=1𝑒^{𝑝2𝑓2𝑖(𝑥)}≤ln𝑛𝑝²+𝑓²(𝑥)^1/2. (48)(49) Taking the limit in (49) as 𝑝→+∞, we get that lim_𝑝→+∞𝐺𝑓_𝑖(𝑥),𝑝=𝑓(𝑥)=max_{1≤𝑖≤𝑛}𝑓_𝑖(𝑥). (50)
(b) Now, we prove that 𝜕𝐺(𝑓_𝑖(𝑥),𝑝)/𝜕𝑝<0. Consider the following: 𝜕𝐺𝑓_𝑖(𝑥),𝑝1𝜕𝑝=−𝑝²ln^𝑛_𝑖=1𝑒^{𝑝2𝑓2𝑖(𝑥)1/2}+ln^𝑛_𝑖=1𝑒^{𝑝2𝑓2𝑖(𝑥)−1/2}⋅∑^𝑛_𝑖=1𝑒^{𝑝2𝑓2𝑖(𝑥)}𝑓²_𝑖(𝑥)∑^𝑛_𝑖=1𝑒^{𝑝2𝑓2𝑖(𝑥)}=𝑝^2𝑛_𝑖=1𝑒^{𝑝2𝑓2𝑖(𝑥)}𝑓²_𝑖(𝑥)−^𝑛_𝑖=1𝑒^{𝑝2𝑓2𝑖(𝑥)}×ln^𝑛_𝑖=1𝑒^{𝑝2𝑓2𝑖(𝑥)}×𝑝²ln^𝑛_𝑖=1𝑒^{𝑝2𝑓2𝑖(𝑥)𝑛1/2}_𝑖=1𝑒^{𝑝2𝑓2𝑖(𝑥)−1}=−^𝑛_𝑖=1𝑒^{𝑝2𝑓2𝑖(𝑥)}ln^𝑛_{𝑗=1𝑖≠𝑗}𝑒^{𝑝2𝑓2𝑗(𝑥)−𝑝2}×𝑝+1₂ln^𝑛_𝑖=1𝑒^{𝑝2𝑓2𝑖(𝑥)𝑛1/2}_𝑖=1𝑒^{𝑝2𝑓2𝑖(𝑥)−1}<0. (51) Since 𝜕𝐺(𝑥,𝑝)/𝜕𝑝<0, then the function 𝐺𝑓_𝑖(𝑥)=1,𝑝𝑝ln^𝑛_𝑖=1𝑒^{𝑝2𝑓2𝑖(𝑥)1/2} (52) is descent for variable 𝑝.
Let 𝐺(𝑓(𝑥),𝑝)=(1/𝑝)ln[𝑛𝑒^{𝑝2𝑓2(𝑥)}]^1/2, then lim_𝑝→+∞𝐺(𝑓(𝑥),𝑝)=𝑓(𝑥).
For an increasing ordered series of numbers 0<𝑝₁<𝑝₂<⋯<𝑝_𝑘<⋯<+∞, (53) we have that 𝐺𝑓(𝑥),𝑝₁≥𝐺𝑓_𝑖(𝑥),𝑝₁≥𝐺𝑓_𝑖(𝑥),𝑝₂≥⋯≥𝐺𝑓_𝑖(𝑥),𝑝_𝑘≥⋯≥lim_𝑝→+∞𝐺𝑓_𝑖(𝑥),𝑝=max_{1≤𝑖≤𝑛}𝑓_𝑖(𝑥),min_{𝑥∈𝑅𝑚}𝐺𝑓(𝑥),𝑝₁≥min_{𝑥∈𝑅𝑚}𝐺𝑓_𝑖(𝑥),𝑝₁≥⋯≥min_{𝑥∈𝑅𝑚}𝐺𝑥,𝑝_𝑘≥⋯≥min_{𝑥∈𝑅𝑚}lim_𝑝→+∞𝐺𝑓_𝑖(𝑥),𝑝. (54)(55) Taking the limit in (55) as 𝑝₁→+∞, then the result (b) follows.

Now, we compare the linear and parabolic aggregation function by the following numerical Example A. In the calculation, we should gradually increase 𝑝 until the numerical overflow in order to determine the range of 𝑝, so that the linear and parabolic aggregation function approximates well the maximum envelope function max_𝑖,𝑗𝑓_𝑖𝑗(𝑥).

Example A. Find𝑥∈𝑅s.t.min_𝑥max_1≤𝑖≤3𝑓_𝑖(𝑥),where𝑓₁(𝑥)=1𝑥,𝑓₂(𝑥)=𝑥²,𝑓₃(𝑥)=6−𝑥. (56)
It can be seen from Figure 5 that the optimal solution is 𝑥^∗=2 and min_𝑥max_1≤𝑖≤3𝑓_𝑖(𝑥^∗)=4 in Example A.

Figure 5 

Three test functions in Example A.

The comparison results for Example A of the linear and parabolic aggregatation function could be seen in Table 1, where the linear aggregatation function is 𝐺(𝑓_𝑖(𝑥),𝑝)=(1/𝑝)[ln(∑^𝑛_𝑖=1𝑒^{𝑝𝑓𝑖(𝑥)})] and the parabolic aggregatation function is 𝐺(𝑓_𝑖(𝑥),𝑝)=(1/𝑝)[ln(∑^𝑛_𝑖=1𝑒^{𝑝2𝑓2𝑖(𝑥)})]^1/2.

Table 1 

Iterative history of two aggregation functions in A.

Figure 5 and Table 1 show that the number of iterations is 17 and the parameter 𝑝≥60 when the satisfactory solution is available for the linear aggregate function. At the same time, for the parabolic aggregate function, the number of iterations is also 17, and the computation is not increased when the parameter 𝑝≥3, and the satisfactory solution is available. But the parabolic aggregate function can better approximate to the maximum envelope function max_𝑖,𝑗𝑓_𝑖𝑗(𝑥) and results in the significant reduction of the number of structural reanalyses so as to make a smaller computational cost. The latter example will illustrate this point in Section 6.

5. Solution Strategy Based on the ICM Method and the Parabolic Aggregate Function

In this section, the optimization problem (41) can be solved by many continuous-type optimization algorithm. However, the simplest choice of the numerical treatment is still the serial Lagrange multiplier method (SLMM) which could efficiently deal with the structural reanalysis of the nonlinear topology optimization problems and obtain an explicit expression of the independent continuous topological design variables.

Now, we will give the process of the structural reanalysis based on the ICM method and the serial Lagrange multiplier method (SLMM).

For convenience, let 𝑓_𝑤(𝑡)=𝑡^𝛼𝑤=𝑡, 𝑓_𝐾(𝑡)=𝑓_𝐸(𝑡)=𝑓_𝜎(𝑡)=𝑡^𝛼𝐾=𝑡^𝛼𝐸=𝑡^𝛼𝜎=𝑡^3.44, and 0<𝑡_min≤𝑡_𝑖≤1.

The element strain energy satisfies the following relationship as 𝑒_𝑖𝑗=𝑝^𝑇_𝑖𝑗𝐾_𝑖⁻¹𝑝_𝑖𝑗2, (57) where 𝐾_𝑖 and 𝑝_𝑖𝑗 are the variable elemental stiffness and the elemental node force vector, respectively.

Substituting the weighting filter function 𝐾_𝑖=𝑓_𝐾(𝑡_𝑖)𝐾⁰_𝑖 in (13) into (57), then we have that 𝑒_𝑖𝑗=𝑝^𝑇_𝑖𝑗𝐾⁰_𝑖⁻¹𝑝_𝑖𝑗2𝑓_𝐾𝑡_𝑖=𝑝^𝑇_𝑖𝑗𝐾⁰_𝑖⁻¹𝑝_𝑖𝑗2𝑡^𝛼𝐾_𝑖. (58)

Assume in the 𝑘th structural analysis that we have that 𝑒^(𝑘)_𝑖𝑗=𝑝^𝑇_𝑖𝑗𝐾⁰_𝑖⁻¹𝑝_𝑖𝑗2𝑓_𝐾𝑡_𝑖^(𝑘)=𝑝^𝑇_𝑖𝑗𝐾⁰_𝑖⁻¹𝑝_𝑖𝑗2𝑡_𝑖^{(𝑘)𝛼𝐾}; (59) then in the (𝑘+1)th structural analysis, and we have that 𝑒_𝑖𝑗=𝑝^𝑇_𝑖𝑗𝐾⁰_𝑖⁻¹𝑝_𝑖𝑗2𝑓_𝐾𝑡_𝑖=𝑝^𝑇_𝑖𝑗𝐾⁰_𝑖⁻¹𝑝_𝑖𝑗2𝑓_𝐾𝑡_𝑖^(𝑘)⋅𝑓_𝐾𝑡_𝑖^(𝑘)𝑓_𝐾𝑡_𝑖=𝑒^(𝑘)_𝑖𝑗⋅𝑓_𝐾𝑡_𝑖^(𝑘)𝑓_𝐾𝑡_𝑖=𝑒^(𝑘)_𝑖𝑗𝑡_𝑖^(𝑘)𝑡_𝑖^𝛼𝐾. (60)

In the context, we define the active set as 𝐼_𝑎={𝑖∣0<𝑡_min<𝑡_𝑖<1,𝑖=1,…,𝑛}, the lower passive set as 𝐼_𝑎={𝑖∣𝑡_𝑖=𝑡_min,𝑖=1,…,𝑛}, and the upper passive set as 𝐼_𝑎={𝑖∣𝑡_𝑖=1,𝑖=1,…,𝑛}. Substituting (59) and (60) into 𝐺(𝑒_𝑖𝑗,𝑝)=(1/𝑝)[ln(∑^𝑛_𝑖=1𝑒^{𝑝2𝑒2𝑖𝑗})]^1/2=𝐺 in (41), then in the (𝑘+1)th structural analysis we have that 𝐹(𝑡)=^𝑛_𝑖=1^𝐽_𝑗=1𝑒^{𝑝2[𝑒(𝑘)𝑖𝑗⋅(𝑓𝐾(𝑡𝑖(𝑘))/𝑓𝐾(𝑡𝑖))]2}−𝐽𝑒^{(𝑝𝑒0𝑖(𝑓2𝜎(𝑡𝑖)/𝑓𝐸(𝑡𝑖)))2}=0. (61)

Let 𝑥_𝑖=[𝑓_𝐾(𝑡_𝑖)]⁻²=[𝑓_𝐸(𝑡_𝑖)]⁻²=(𝑡_𝑖)^−2𝛼𝐾, 𝑎^(𝑘)_𝑖𝑗=𝑝²[𝑒^(𝑘)_𝑖𝑗𝑓_𝐾(𝑡_𝑖^(𝑘))]²=𝑝²[𝑒^(𝑘)_𝑖𝑗(𝑡_𝑖^(𝑘))^𝛼𝐾]², and 𝑞=𝑝²; then (58) is equivalent with 𝐹(𝑡)=𝑔(𝑥)=_{𝑖∈𝐼𝑎}^𝐽_𝑗=1𝑒^{𝑞(𝑒𝑖𝑗(𝑘))2}−𝐽𝑒^{𝑞(𝑒0𝑖⋅𝑡𝜎𝐸2𝛼−𝛼min)2}+_{𝑖∈𝐼𝑎}^𝐽_𝑗=1𝑒^{𝑞𝑎(𝑘)𝑖𝑗𝑥𝑖}−𝐽𝑒^{𝑞(𝑒0𝑖)2⋅𝑥𝜎𝐸𝑘𝑖−(2𝛼−𝛼)/𝛼}+_{𝑖∈𝐼𝑎}^𝐽_𝑗=1𝑒^{𝑞(𝑒𝑖𝑗(𝑘))2}−𝐽𝑒^{𝑞(𝑒0𝑖)2}. (62)

In order to achieve the constraint in (59) explicitly, we demand the first-order partial derivatives of 𝑔(𝑥)   on the design variables: 𝜕𝑔𝜕𝑥_𝑖=^𝐽_𝑗=1𝑒^{𝑞𝑎(𝑘)𝑖𝑗𝑥𝑖}𝑞𝑎^(𝑘)_𝑖𝑗+𝐽𝑒^{𝑞(𝑒0𝑖)2⋅𝑥𝜎𝐸𝑘𝑖−(2𝛼−𝛼)/𝛼}𝑞𝑒⁰_𝑖²×2𝛼_𝜎−𝛼_𝐸𝛼_𝐾𝑥^{(−(2𝛼𝜎−𝛼𝐸)−𝛼𝐾)/𝛼𝐾}_𝑖. (63)

Then, the first-order Taylor expansion formula of 𝑔(𝑥) is 𝑔(𝑥)=𝑔𝑥^(𝑘)+_{𝑖∈𝐼𝑎}𝜕𝑔𝑥^(𝑘)𝜕𝑥_𝑖𝑥_𝑖−𝑥_𝑖^(𝑘). (64) Substituting the intermediate variable 𝑥_𝑖=𝑡^−2𝛼𝑘_𝑖 into (64), then we have that 𝑔(𝑥)=𝑏₀^(𝑘)+_{𝑖∈𝐼𝑎}𝑏_𝑖^(𝑘)𝑡^−2𝛼𝑘_𝑖, (65) where 𝑏₀^(𝑘)=𝑔𝑥^(𝑘)−_{𝑖∈𝐼𝑎}𝜕𝑔𝑥^(𝑘)𝜕𝑥_𝑖𝑥_𝑖^(𝑘),𝑏_𝑖^(𝑘)=𝜕𝑔𝑥^(𝑘)𝜕𝑥_𝑖. (66)

Hence, the approximate model of (41) can be expressed as follows.

Find topology variables 𝐭∈[0,1]^𝑛; minimize𝑊(𝑡)=^𝑛_𝑖=1𝑓_𝑊𝑡_𝑖𝑊⁰_𝑖=_{𝑖∈𝐼𝑎}𝑤⁰_𝑖𝑡^𝛼𝑤_min+_{𝑖∈𝐼𝑎}𝑤⁰_𝑖𝑡^𝛼𝑤_𝑖+_{𝑖∈𝐼𝑎}𝑤⁰_𝑖,subjectto𝑏₀^(𝑘)+_{𝑖∈𝐼𝑎}𝑏_𝑖^(𝑘)𝑡^−2𝛼𝑘_𝑖=0,0<𝑡_min≤𝑡_𝑖≤1,𝑖∈𝐼_𝑛={1,…,𝑛},𝑗∈𝐽_𝑚={1,…,𝐽}. (67)

Setting aside consideration of interval constraint 𝑡_min≤𝑡_𝑖≤1, then the augmented Lagrange function of (67) can be expressed as 𝐿(𝐭,𝜆)=_{𝑖∈𝐼𝑎}𝑤⁰_𝑖𝑡^𝛼𝑤_min+_{𝑖∈𝐼𝑎}𝑤⁰_𝑖𝑡^𝛼𝑤_𝑖+_{𝑖∈𝐼𝑎}𝑤⁰_𝑖+𝜆𝑏₀^(𝑘)+_{𝑖∈𝐼𝑎}𝑏_𝑖^(𝑘)𝑡^−2𝛼𝐾_𝑖. (68)

Taking extreme conditions for the saddle point in (68), 𝜕𝐿𝜕𝑡_𝑚=𝛼_𝑤𝑤⁰_𝑚𝑡^𝛼𝑤_𝑚⁻¹−2𝜆𝛼_𝐾𝑏_𝑚^(𝑘)𝑡^−2𝛼𝐾_𝑚⁻¹=0,𝑚∈𝐼_𝑎,𝜕𝐿𝜕𝜆=𝑏₀^(𝑘)+_{𝑖∈𝐼𝑎}𝑏_𝑖^(𝑘)𝑡^−2𝛼𝐾_𝑖=0. (69)(70)

By (69) we can have that 𝑡_𝑚=𝛼2𝜆⋅^𝐾𝑏_𝑚^(𝑘)𝛼_𝑤𝑤⁰_𝑚^{1/(2𝛼𝐾+𝛼𝑤)},𝑚∈𝐼_𝑎. (71)

Substituting (71) into (70), we can have that 𝜆=−𝑏₀^(𝑘)∑_{𝑖∈𝐼𝑎}𝑏_𝑖^(𝑘)2𝛼_𝐾𝑏_𝑖^(𝑘)/𝛼_𝑤𝑤⁰_𝑖^{−2𝛼𝐾/(2𝛼𝐾+𝛼𝑤)−(2𝛼𝐾+𝛼𝑤)/2𝛼𝐾}. (72)

Substituting (72) into (71), we can have that 𝑡_𝑚̃=−𝑏₀^(𝑘)∑_{𝑖∈𝐼𝑎}𝑏_𝑖^(𝑘)2𝛼_𝐾𝑏_𝑖^(𝑘)/𝛼_𝑤𝑤⁰_𝑖^{−2𝛼𝐾/(2𝛼𝐾+𝛼𝑤)−1/2𝛼𝐾}×2𝛼_𝐾𝑏_𝑚^(𝑘)𝛼_𝑤𝑤⁰_𝑚^{1/(2𝛼𝐾+𝛼𝑤)}. (73)

With consideration on interval constraint 𝑡_min≤𝑡_𝑖≤1(𝑖=1,…,𝑛),  we can have that 𝑡_𝑚^(𝑘+1)=𝑡_min,𝑡_𝑚̃≤𝑡_min,𝑡_𝑚̃,𝑡_min<𝑡_𝑚̃,𝑡̃<11,_𝑚≥1. (74)

Record and update the active set 𝐼_𝑎 and the passive set 𝐼_𝑎 and 𝐼_𝑎; then return to (73) to calculate 𝑡_𝑚̃, this process is repeated until the same set of active and passive sets does not change and meets the convergence criteria for iteration: 𝑊^(𝑘+1)−𝑊^(𝑘)𝑊^(𝑘)≤𝜀, (75) where 𝑊^(𝑘) and 𝑊^(𝑘+1) are the total weight of the structure of previous and current iterations, respectively, and 𝜀=0.001 is the convergence precision.

6. Numerical Examples

In this section, we present three numerical results to demonstrate the capabilities of the proposed methods. The structural optimization problem in all examples is to solve the optimization model (67) so as to obtain the optimal topologies of the original topology optimization problem with stress constraints in (32) and its equivalent optimization model with strain energy constraints in (41). Actually, all the examples are two-dimensional. However, the three-dimensional problems would work as well in principle.

6.1. The Short Cantilever Beam

The first example is the short cantilever beam shown in Figure 6(a), with a clamped left hand side, and a vertical load applied at the middle of the right hand side. The rectangular domain (0.1 m thick) containing the structure is discretized into 60 × 20  meshes by using four-node membrane element. In the model, the density of material is used to be 𝜌=1.0 kg/m³, the Poisson’s radio is 𝜈=0.3, and the elastic module of material is 𝐸=3×10⁵ N/m². In order to avoid stress concentration, the load 𝑃=6.0 N is located at the three nodes in the middle of the right side. The allowable stress is 100 N, and the convergence precision is limited to be 0.001.

(a)

(b)

Figure 6 

(a) Geometry and load condition for the short cantilever beam problem. (b) The optimal topology structure for the short cantilever beam problem.

(a)

(b)

(c)

Figure 7 

(a) The elemental strain energy 𝑒 distribution for the short cantilever beam problem. (b) The independent continuous topological design variables 𝑡-distribution for the short cantilever beam problem. (c) Stress 𝜎 distribution for the short cantilever beam problem.

Thus, the optimal topology structure is shown in Figure 6(b), and the historical curves of the elemental maximum stress, the structural strain energy, the elemental maximum strain energy, and the structural weight are shown in Figures 7, 8, 9, and 10, respectively. The optimal design was obtained in 69 iterations, where the number of elements reserved is 383 (see Table 2), the total structure strain energy is 0.0290543 J, the average elemental strain energy of the element reserved is 7.58597 E-5 J, and the elemental maximum strain energy is 0.000246914 J. From strength considerations in the optimal topology structure, the maximum stress of the structural elements is at the node 609, and its value is 81.9 N which meets the strength requirements, and the total weight of the optimal structure is 2.4827 kg, in which the parameter of the parabolic aggregatation function is 𝑝=√2.

Table 2 

The distribution of the independent continuous topological design variables in the optimal topology structure of the short cantilever beam.

Figure 8 

The historical curves of the elemental maximum stress and the structural weight for the short cantilever beam problem.

Figure 9 

The historical curves of the structural strain energy and the structural weight for the short cantilever beam problem.

Figure 10 

The historical curves of the elemental maximum strain energy and the structural weight for the short cantilever beam problem.

6.2. The MBB-Beam

The second example is the MBB-beam shown in Figure 11 (cf. [63]), in which material constants are 𝐸=2.1𝑒8 pa, 𝜈=0.3, and 𝜌=1.0 kg/m³, where a vertically force 𝑃=100 N  is applied at the center point of the upper side, and the bottom-left corner is supported as a roller, while the bottom-right corner is a fixed support. The design domain (1 m thick) is discretized into 180 × 30  meshes by Q4 elements. The allowable stress is 1.0 kN,  and the convergence precision is limited to be 0.001.

Figure 11 

Geometry and load condition for the MBB-beam problem.

Thus, the optimal topology structure is shown in Figure 12, and the historical curves of the elemental maximum stress, the structural strain energy, the elemental maximum strain energy, and the structural weight are shown in Figures 13, 14, 15, and 16, respectively. The optimal design was obtained in 46 iterations, where the number of elements reserved is 2550 (see Table 3), the total structure strain energy is 0.0519843 J, the average elemental strain energy of the element reserved is 2.0386E-5 J, and the elemental maximum strain energy is 21.446 E-5 J. From strength considerations in the optimal topology structure, the maximum stress of the structural elements is at the node 5520, and its value is 975 N which meets the strength requirements, and the total weight of the optimal structure is 88.2662 kg, in which the parameter of the parabolic aggregatation function is 𝑝=2.

Table 3 

The distribution of the independent continuous topological design variables in the optimal topology structure of the MBB-beam.

Figure 12 

The optimal topology structure for the MBB-beam problem.

(a)

(b)

(c)

Figure 13 

(a) The elemental strain energy 𝑒 distribution for the MBB-beam problem. (b) The independent continuous topological design variables 𝑡-distribution for the MBB-beam problem. (c) Stress 𝜎 distribution for the MBB-beam problem.

Figure 14 

The historical curves of the elemental maximum stress and the structural weight for the MBB-beam problem.

Figure 15 

The historical curves of the structural strain energy and the structural weight for the MBB-beam problem.

Figure 16 

The historical curves of the elemental maximum strain energy and the structural weight for the MBB-beam problem.

6.3. The L-Shape Beam

The third example is the L-shape beam shown in Figure 17 (cf. [63]), which is clamped on its top boundary and loaded on the mid 1/3 of the right side by the uniformly distributed load 𝑃. The material constants are 𝐸=21000 pa, 𝜈=0.3, 𝜌=7800 kg/m³, and 𝑙=6.0 m,  and the uniformly distributed load is 𝑃=0.3 N. The L-shape domain (0.027 m thick) containing the structure is discretized into 18 × 18 × 3  meshes by using four-node membrane element with 972 elements. The allowable stress is 70 N, and the convergence precision is limited to be 0.001.

Figure 17 

Geometry and load condition for the L-shape beam problem.

Thus, the optimal topology structure is shown in Figure 18, and the historical curves of the elemental maximum stress, the structural strain energy, the elemental maximum strain energy, and the structural weight are shown in Figures 19, 20, 21, and 22, respectively. The optimal design was obtained in 32 iterations, where the number of elements reserved is 586 (see Table 4), the total structure strain energy is 0.0221408 J, the average elemental strain energy of the element reserved is 3.7783E-5 J, and the elemental maximum strain energy is 0.000343573 J. From strength considerations in the optimal topology structure, the maximum stress of the structural elements is at the node 686, and its value is 160 N which meets the strength requirements, and the total weight of the optimal structure is 1762.9 kg, in which the parameter of the parabolic aggregatation function is 𝑝=√10.

Table 4 

The distribution of the independent continuous topological design variables in the optimal topology structure of the L-shape beam.

Figure 18 

The optimal topology structure for the L-shape beam problem.

(a)

(b)

(c)

Figure 19 

(a) The elemental strain energy 𝑒 distribution for the L-shape beam problem. (b) The independent continuous topological design variables 𝑡-distribution for the L-shape beam problem. (c) Stress 𝜎 distribution for the L-shape beam problem.

Figure 20 

The historical curves of the elemental maximum stress and the structural weight for the L-shape beam problem.

Figure 21 

The historical curves of the structural strain energy and the structural weight for the L-shape beam problem.

Figure 22 

The historical curves of the elemental maximum strain energy and the structural weight for the L-shape beam problem.

7. Conclusions

This paper develops three novel methods for solving stress-based topology optimization problems, which are the independent-continuous-mapping (ICM) method, an alternative strain energy method, and the parabolic aggregation function, that is, different from the two common global functions such as𝑝-norm and KS function. The numerical examples demonstrate that the proposed methods effectively remove the stress concentrations and generate black-and-white designs for practically sized problems.

Acknowledgments

This paper is supported by the National Natural Science Foundation of China (10872012, 70901010, and 11172013), the Foundation of National Key Laboratory for Structural Analysis of Industrial Equipment (Dalian University of Technology) (GZ1008), and the Postdoctoral Science Foundation of Beijing University of Technology (X0001015200803).