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

Authors: Tie Jun; Sui Yun-kang
Journal: Mathematical Problems in Engineering
 (2013)

---

## 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.


---

## Body

## 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 Section2, 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 p-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 ℝ3 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 χ(x) is an indicator function that defines the distribution of material [1, 6]
(1)χ(x)={1,ifx∈Ωsolid,0,ifx∈Ωvoid,(2)Ωsolid={x∣χ(x)=1},Ωvoid={x∣χ(x)=0}.Suppose that the fixed design domainΩ is discretized into a ground structure with n finite elements, we can define the binary design variables ti∈{0,1}  as ti=1 when the ith element is filled with material and ti=0 when the ith element is void such that the problem is then to decide the binary design variable vector t=(t1,t2,…,tn)T∈{0,1}n.Now, we consider the general binary discrete optimization problem with some property constraints as follows.Find topological variablest∈[0,1]n,
(3)minimizeΦ=Φ(t)=∑i=1nϕi,subjecttogj(τ1,…,τn)-g-j≤0,j=1,…,J,0≤ti≤1,i=1,…,n,
where Φ=Φ(t)=∑i=1nϕi  is the objective function, τi is the physical or geometric parameter of the ith element used in the finite analysis, gj(τ1,…,τn) and g-j are the jth constraint response and its upper bound, respectively, J is the total number of constraints, ti is the ith independent topological design variable, t=(t1,…,tn)T is the independent topological design vector, and the total number of design variables is n. In fact, τi may be regarded as ti  function with  τi=0  as ti=0 and τi≠0 as ti≠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,
(4)ti≈ρi,0≤ρi≤1,ti≈di,0≤di≤1,ti≈Ai,0≤Ai≤1,
where ρi, di, and Ai are the ith element normalized density, thickness, and area, respectively, and p 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 variablesti∈{0,1} in (1) can be relaxed to a new step-up function and be defined as (see Figure 1)
(5)ti=s(τiτi0)={1,τiτi0∈(0,1],0,τiτi0=0,
where τi is the physical or geometric parameter of the ith element used in the finite analysis, τi0 is the original and true physical or geometric parameter of the material used of the ith element.Figure 1
The step-up function and traditional topology design variable.As an appropriate relaxation,ti=s(τi/τi0) is allowed to take all possible values between 0 and 1, and it is clear that ti=s(τi/τi0) is an ideal independent topological design variable which has clear physical meaning. However, the step-up function (5) is discontinuous in the point τi/τi0=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 ith element τi is not equal to 0, then ti=s(τi/τi0) is equal to 1, and when τi 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 ti and τi/τi0.


### 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 functionp(τi/τi0) to approximate the step-up function (5) such that ti=s(τi/τi0) becomes a continuous and differentiable function (see Figure 2); that is,
(6)ti=s(τiτi0)≈p(τiτi0),
where the following properties are satisfied.(1)
Continuity and differentiability:p(τi/τi0) is a continuous and differentiable function in [0,1].(2)
Strict increasing monotonicity:p(τi/τi0) is a strictly monotonically increasing function.(3)
Convexity:p(τi/τi0) is a convex function in [0,1].(4)
Boundedness:p(τi/τi0)∈[0,1], p(0)=0, and p(1)=1.(5)
Approximation:∫01|1-p(x)|dx<1/2 or x≤p(x)≤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,ti=p(τi/τi0) can take values in the entire range [0,1] instead of the discrete set {0,1} and make the functional relationship ti=s(τi/τi0) between the topology design variables ti, and the various physical or geometric parameters τi 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 ti=p(τi/τi0) 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 variablesti=p(τi/τi0) take intermediate values between 0-1, ti=p(τi/τi0) not only reflect the proximity to corresponding elements “have” and “no”, but also identify the size of the various physical values τi 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)f:τiτi0⟼p-1(ti);
that is,
(8)τiτi0=f(ti)=p-1(ti)orτi=f(ti)τi0=p-1(ti)τi0.Moreover, we define the inverse function of the polish functionf(ti)=τi/τi0=p-1(ti) as the weighting filter function, where ti∈[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(9)x=tα,α=3,10.32;(b) composite exponential function(10)x=(et/γ-1)(e1/γ-1),γ=0.0621;(c) modified sigmoid function(11)x=122ln1+t/β1-t/β,β=1.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 functionf(ti)=τi/τi0=p-1(ti) plays a very important role as follows.(a) Let the inverse mapping of the step-up function in (5) denote as τi/τi0=s-1(ti), which is a 0-1 binary multivalued function that cannot indicate the size of the physical or geometric parameter τi when the independent topological design variables ti=s(τi/τi0) take intermediate values between 0-1. To overcome this difficulty, we introduce the concept of weighting filter function f(ti)=τi/τi0=p-1(ti) which continuously approximates the inverse mapping of the step-up function τi/τi0=s-1(ti), and it points out the proximity degree that τi closes to τi0 and 0.(b) According to the ICM method, ifτi=wi or Ki, then the variable weight and stiffness matrix of the ith element, wi and Ki are, respectively, recognized by the weighting filter functions as
(12)wi=fw(ti)wi0,(13)Ki=fK(ti)Ki0,
where wi0 and Ki0 are the original weight and stiffness matrix of the ith element for structural material, respectively.Generally, whenτi in (8) represents the different variable physical or geometric parameter, the weighting filter functions could have different forms; that is,
(14)σ-i=fσ(ti)σ-i0,(15)mi=fm(ti)mi0,(16)vi=fv(ti)vi0,(17)Ai=fA(ti)Ai0,(18)di=fd(ti)di0,
where  σ-i, mi, vi, Ai, and di are, respectively, the variable allowable stress, mass, volume, area, and thickness of the ith element for structural material, and σ-i0, mi0, vi0, Ai0, and di0 are the original allowable stress, mass, volume, area, and thickness of the ith 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)ρi=fρ(ti)ρi0,(20)Ei=fE(ti)Ei0,
where ρi and Ei are the variable density and Young’s modulus of the ith element for structural material and ρi0  and Ei0 are the original density and Young’s modulus of the ith 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 variableti=p(τi/τi0) instead of the relative density ρi 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ϕi be the objective function of the ith element used in the finite analysis and ϕi0 the original objective function of the ith element for structural material, then we have that
(21)ϕi=fϕ(ti)ϕi0
and the topology optimization model (3) could be approximated by the new continuous nonlinear topology optimization model; as follows.Find topological variablest∈[0,1]n(22)minimizeΦ=∑i=1nϕi0fϕ(ti)subjecttogj(τ10fτ(t1),…,τn0fτ(tn))-g-j≤0,j=1,…,J,0≤ti≤1,i=1,…,n.(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 variablesti 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σ-i of the ith element, which is recognized by the allowable stress weighting filter functions as (14)
(23)σ-i=fσ(ti)σ-i0
such that
(24)σij≤fσ(ti)σ-i0,
where σ-i0 is the original allowable stress of the ith element for structural material and σij is the von Mises stress for the ith element under the jth load case. That is, the inequality (24) guarantees that σij→0, when ti 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 ρpσl, where they apply a stress constraint for the SIMP model, that is, expressed as a constraint of the von Mises equivalent stress σVM:
(25)σVM≤ρpσl,ifρ>0.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)σij≤σ-i,
then we have that
(27)σij≤σ-i=fσ(ti)σ-i0,(28)σij<ε+fσ(ti)σ-i0,
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)(σij-σ-i0)fσ(ti)=σijfσ(ti)-fσ(ti)σ-i0≤σij-fσ(ti)σ-i0<ε;
that is,
(30)(σij-σ-i0)fσ(ti)<ε.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.


## 2.1. Domain Setting and the Independent Topological Design Variables

Let us consider a fixed design domainΩ in ℝ3 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 χ(x) is an indicator function that defines the distribution of material [1, 6]
(1)χ(x)={1,ifx∈Ωsolid,0,ifx∈Ωvoid,(2)Ωsolid={x∣χ(x)=1},Ωvoid={x∣χ(x)=0}.Suppose that the fixed design domainΩ is discretized into a ground structure with n finite elements, we can define the binary design variables ti∈{0,1}  as ti=1 when the ith element is filled with material and ti=0 when the ith element is void such that the problem is then to decide the binary design variable vector t=(t1,t2,…,tn)T∈{0,1}n.Now, we consider the general binary discrete optimization problem with some property constraints as follows.Find topological variablest∈[0,1]n,
(3)minimizeΦ=Φ(t)=∑i=1nϕi,subjecttogj(τ1,…,τn)-g-j≤0,j=1,…,J,0≤ti≤1,i=1,…,n,
where Φ=Φ(t)=∑i=1nϕi  is the objective function, τi is the physical or geometric parameter of the ith element used in the finite analysis, gj(τ1,…,τn) and g-j are the jth constraint response and its upper bound, respectively, J is the total number of constraints, ti is the ith independent topological design variable, t=(t1,…,tn)T is the independent topological design vector, and the total number of design variables is n. In fact, τi may be regarded as ti  function with  τi=0  as ti=0 and τi≠0 as ti≠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,
(4)ti≈ρi,0≤ρi≤1,ti≈di,0≤di≤1,ti≈Ai,0≤Ai≤1,
where ρi, di, and Ai are the ith element normalized density, thickness, and area, respectively, and p 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 variablesti∈{0,1} in (1) can be relaxed to a new step-up function and be defined as (see Figure 1)
(5)ti=s(τiτi0)={1,τiτi0∈(0,1],0,τiτi0=0,
where τi is the physical or geometric parameter of the ith element used in the finite analysis, τi0 is the original and true physical or geometric parameter of the material used of the ith element.Figure 1
The step-up function and traditional topology design variable.As an appropriate relaxation,ti=s(τi/τi0) is allowed to take all possible values between 0 and 1, and it is clear that ti=s(τi/τi0) is an ideal independent topological design variable which has clear physical meaning. However, the step-up function (5) is discontinuous in the point τi/τi0=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 ith element τi is not equal to 0, then ti=s(τi/τi0) is equal to 1, and when τi 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 ti and τi/τi0.


## 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 functionp(τi/τi0) to approximate the step-up function (5) such that ti=s(τi/τi0) becomes a continuous and differentiable function (see Figure 2); that is,
(6)ti=s(τiτi0)≈p(τiτi0),
where the following properties are satisfied.(1)
Continuity and differentiability:p(τi/τi0) is a continuous and differentiable function in [0,1].(2)
Strict increasing monotonicity:p(τi/τi0) is a strictly monotonically increasing function.(3)
Convexity:p(τi/τi0) is a convex function in [0,1].(4)
Boundedness:p(τi/τi0)∈[0,1], p(0)=0, and p(1)=1.(5)
Approximation:∫01|1-p(x)|dx<1/2 or x≤p(x)≤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,ti=p(τi/τi0) can take values in the entire range [0,1] instead of the discrete set {0,1} and make the functional relationship ti=s(τi/τi0) between the topology design variables ti, and the various physical or geometric parameters τi 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 ti=p(τi/τi0) 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 variablesti=p(τi/τi0) take intermediate values between 0-1, ti=p(τi/τi0) not only reflect the proximity to corresponding elements “have” and “no”, but also identify the size of the various physical values τi 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)f:τiτi0⟼p-1(ti);
that is,
(8)τiτi0=f(ti)=p-1(ti)orτi=f(ti)τi0=p-1(ti)τi0.Moreover, we define the inverse function of the polish functionf(ti)=τi/τi0=p-1(ti) as the weighting filter function, where ti∈[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(9)x=tα,α=3,10.32;(b) composite exponential function(10)x=(et/γ-1)(e1/γ-1),γ=0.0621;(c) modified sigmoid function(11)x=122ln1+t/β1-t/β,β=1.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 functionf(ti)=τi/τi0=p-1(ti) plays a very important role as follows.(a) Let the inverse mapping of the step-up function in (5) denote as τi/τi0=s-1(ti), which is a 0-1 binary multivalued function that cannot indicate the size of the physical or geometric parameter τi when the independent topological design variables ti=s(τi/τi0) take intermediate values between 0-1. To overcome this difficulty, we introduce the concept of weighting filter function f(ti)=τi/τi0=p-1(ti) which continuously approximates the inverse mapping of the step-up function τi/τi0=s-1(ti), and it points out the proximity degree that τi closes to τi0 and 0.(b) According to the ICM method, ifτi=wi or Ki, then the variable weight and stiffness matrix of the ith element, wi and Ki are, respectively, recognized by the weighting filter functions as
(12)wi=fw(ti)wi0,(13)Ki=fK(ti)Ki0,
where wi0 and Ki0 are the original weight and stiffness matrix of the ith element for structural material, respectively.Generally, whenτi in (8) represents the different variable physical or geometric parameter, the weighting filter functions could have different forms; that is,
(14)σ-i=fσ(ti)σ-i0,(15)mi=fm(ti)mi0,(16)vi=fv(ti)vi0,(17)Ai=fA(ti)Ai0,(18)di=fd(ti)di0,
where  σ-i, mi, vi, Ai, and di are, respectively, the variable allowable stress, mass, volume, area, and thickness of the ith element for structural material, and σ-i0, mi0, vi0, Ai0, and di0 are the original allowable stress, mass, volume, area, and thickness of the ith 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)ρi=fρ(ti)ρi0,(20)Ei=fE(ti)Ei0,
where ρi and Ei are the variable density and Young’s modulus of the ith element for structural material and ρi0  and Ei0 are the original density and Young’s modulus of the ith 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 variableti=p(τi/τi0) instead of the relative density ρi 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ϕi be the objective function of the ith element used in the finite analysis and ϕi0 the original objective function of the ith element for structural material, then we have that
(21)ϕi=fϕ(ti)ϕi0
and the topology optimization model (3) could be approximated by the new continuous nonlinear topology optimization model; as follows.Find topological variablest∈[0,1]n(22)minimizeΦ=∑i=1nϕi0fϕ(ti)subjecttogj(τ10fτ(t1),…,τn0fτ(tn))-g-j≤0,j=1,…,J,0≤ti≤1,i=1,…,n.(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 variablesti 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σ-i of the ith element, which is recognized by the allowable stress weighting filter functions as (14)
(23)σ-i=fσ(ti)σ-i0
such that
(24)σij≤fσ(ti)σ-i0,
where σ-i0 is the original allowable stress of the ith element for structural material and σij is the von Mises stress for the ith element under the jth load case. That is, the inequality (24) guarantees that σij→0, when ti 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 ρpσl, where they apply a stress constraint for the SIMP model, that is, expressed as a constraint of the von Mises equivalent stress σVM:
(25)σVM≤ρpσl,ifρ>0.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)σij≤σ-i,
then we have that
(27)σij≤σ-i=fσ(ti)σ-i0,(28)σij<ε+fσ(ti)σ-i0,
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)(σij-σ-i0)fσ(ti)=σijfσ(ti)-fσ(ti)σ-i0≤σij-fσ(ti)σ-i0<ε;
that is,
(30)(σij-σ-i0)fσ(ti)<ε.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 σij be the von Mises stress for the ith element under the jth 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 dΩ by fw(ti), and the objective function Φ=∑i=1nϕi0fΦ(ti) may be defined as
(31)W(t)=∫ΩsolidγmatdΩ=∫Ωfw(t)γmatdΩ=∑i=1nwi=∑i=1nfw(ti)wi0.
Then, the original topology optimization problem with stress constraints can be formulated as follows.Find topological variablest∈[0,1]n;
(32)minimizeW(t)=∑i=1nfw(ti)wi0,subjecttoσij≤σ-i=fσ(ti)σ-i0,0<tmin≤ti≤1,i=1,…,n,j=1,…,J,
where ti is the ith independent continuous topological design variable following (8), t=(t1,…,tn) is the independent continuous topological design vector, the total number of design variables is n,  and the total number of external loads is J. As it is well known, a lower limit tmin of the the ith independent continuous topological design variable is imposed to avoid the singularity of the global stiffness matrix which is given as tmin=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 withn elements will be associated with n×J stress when there is J 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(33)σVM=12[(σ1-σ2)2+(σ2-σ3)2+(σ3-σ1)2],
where σ1, σ2, and σ3 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
(34)uf≤1+ν3E(σ-)2,(35)uf=1+ν6E[(σ1-σ2)2+(σ2-σ3)2+(σ3-σ1)2],
where E, ν 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(36)uf=uγ2,γ≥1,uij=eijvi,
where uij, eij, and vi are the strain energy density, the strain energy, and the volume of the ith element under the jth 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.LetEi be the variable Young’s modulus of the ith element for structural material, substituting (36) into (35), we have that
(37)eij≤1+ν3Ei(γσ-i)2vi=e-ij,
where we remark that e-ij is defined as the variable allowable strain energy of the ith element which is unique under multiload case and may be denoted by e-i=e-ij.It is evident that the structure is safe, and the stress constraintsσij≤σ-i0 are satisfied if the formulation (37) is established.Moreover, substituting (14) and (20) into e-i=e-ij in (37), we have that
(38)fE(ti)e-i=fσ2(ti)e-i0,e-i0=1+ν3Ei0(γσ-i0)2vi,
where e-i0 is defined as the original allowable strain energy of the ith 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 eij≤e-i in (37) instead of the stress constraints σij≤σ-i0 in (32) as following.Find topological variablest∈[0,1]n;
(39)minimizeW(t)=∑i=1nfw(ti)wi0,subjecttoeij≤e-i,0<tmin≤ti≤1,i∈In={1,…,n},j∈Jm={1,…,m}.By the properties of the parabolic aggregation function (see the next section), we have that(40)limp→∞G(eij,p)=maxi,jeij,limp→∞G(e-i,p)=maxie-i.Substituting (40) into (39), the optimization model (32) is equivalent with the following when the parameter p is enough big.Find topological variablest∈[0,1]n;
(41)minimizeW(t)=∑i=1nfw(ti)wi0,subjecttoG(eij,p)≤G(e-i,p)=G-,0<tmin≤ti≤1,i∈In={1,…,n},j∈Jm={1,…,J}.


## 3.1. Minimum Weight with Stress Constraints in Topology Optimization

For the general explicit continuous nonlinear topology optimization model (22), let σij be the von Mises stress for the ith element under the jth 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 dΩ by fw(ti), and the objective function Φ=∑i=1nϕi0fΦ(ti) may be defined as
(31)W(t)=∫ΩsolidγmatdΩ=∫Ωfw(t)γmatdΩ=∑i=1nwi=∑i=1nfw(ti)wi0.
Then, the original topology optimization problem with stress constraints can be formulated as follows.Find topological variablest∈[0,1]n;
(32)minimizeW(t)=∑i=1nfw(ti)wi0,subjecttoσij≤σ-i=fσ(ti)σ-i0,0<tmin≤ti≤1,i=1,…,n,j=1,…,J,
where ti is the ith independent continuous topological design variable following (8), t=(t1,…,tn) is the independent continuous topological design vector, the total number of design variables is n,  and the total number of external loads is J. As it is well known, a lower limit tmin of the the ith independent continuous topological design variable is imposed to avoid the singularity of the global stiffness matrix which is given as tmin=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 withn elements will be associated with n×J stress when there is J 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(33)σVM=12[(σ1-σ2)2+(σ2-σ3)2+(σ3-σ1)2],
where σ1, σ2, and σ3 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
(34)uf≤1+ν3E(σ-)2,(35)uf=1+ν6E[(σ1-σ2)2+(σ2-σ3)2+(σ3-σ1)2],
where E, ν 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(36)uf=uγ2,γ≥1,uij=eijvi,
where uij, eij, and vi are the strain energy density, the strain energy, and the volume of the ith element under the jth 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.LetEi be the variable Young’s modulus of the ith element for structural material, substituting (36) into (35), we have that
(37)eij≤1+ν3Ei(γσ-i)2vi=e-ij,
where we remark that e-ij is defined as the variable allowable strain energy of the ith element which is unique under multiload case and may be denoted by e-i=e-ij.It is evident that the structure is safe, and the stress constraintsσij≤σ-i0 are satisfied if the formulation (37) is established.Moreover, substituting (14) and (20) into e-i=e-ij in (37), we have that
(38)fE(ti)e-i=fσ2(ti)e-i0,e-i0=1+ν3Ei0(γσ-i0)2vi,
where e-i0 is defined as the original allowable strain energy of the ith 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 eij≤e-i in (37) instead of the stress constraints σij≤σ-i0 in (32) as following.Find topological variablest∈[0,1]n;
(39)minimizeW(t)=∑i=1nfw(ti)wi0,subjecttoeij≤e-i,0<tmin≤ti≤1,i∈In={1,…,n},j∈Jm={1,…,m}.By the properties of the parabolic aggregation function (see the next section), we have that(40)limp→∞G(eij,p)=maxi,jeij,limp→∞G(e-i,p)=maxie-i.Substituting (40) into (39), the optimization model (32) is equivalent with the following when the parameter p is enough big.Find topological variablest∈[0,1]n;
(41)minimizeW(t)=∑i=1nfw(ti)wi0,subjecttoG(eij,p)≤G(e-i,p)=G-,0<tmin≤ti≤1,i∈In={1,…,n},j∈Jm={1,…,J}.


## 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 eij 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 p-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(42)G(eij,p)=1p[ln(∑i=1n∑j=1mep2eij2)]1/2,
where G(eij,p) 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 p 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 thatG-=(1/p)[ln(∑i=1n∑j=1Je(pe-i0fσ2(ti))2/fE(ti))]1/2 as the variable global allowable strain energy.Theorem 1.
Letv=eu2p2, if the following conditions are satisfied: (1)
u>0,   p>0;(2)
fj(x)>0,  j=1,2,…,n,x∈Rm.And let G(fi(x),p)=(1/p)[ln(∑i=1nep2fi2(x))]1/2; then one has that
(a)(43)limp→+∞G(fi(x),p)=max1≤i≤nfi(x);
(b)(44)minx∈Rmmax1≤i≤nfi(x)=limp→+∞minx∈RmG(fi(x),p).Proof.
(a) Ifv=eu2p2, then we have that
(45)u=1p(lnv)1/2.
For
(46)fi(x)>0,p>0,v[fi(x),p]=ep2fi2(x).
Let max1≤i≤nfi(x)=f-(x). Since u>0, p>0, then we have that
(47)v[f-(x),p]≤∑i=1nv[fi(x),p]≤nv[f-(x),p];
that is,
(48)ep2(f-(x))2≤∑i=1nep2fi2(x)≤nep2(f-(x))2,(49)f-(x)≤1pln[∑i=1nep2fi2(x)]≤[lnnp2+f-2(x)]1/2.
Taking the limit in (49) as p→+∞, we get that
(50)limp→+∞G(fi(x),p)=f-(x)=max1≤i≤nfi(x).
(b) Now, we prove that∂G(fi(x),p)/∂p<0. Consider the following:
(51)∂G(fi(x),p)∂p=-1p2(ln[∑i=1nep2fi2(x)])1/2+(ln[∑i=1nep2fi2(x)])-1/2·(∑i=1nep2fi2(x))fi2(x)∑i=1nep2fi2(x)=(p2(∑i=1nep2fi2(x))fi2(x)-(∑i=1nep2fi2(x))×(ln[∑i=1nep2fi2(x)]))×(p2(ln[∑i=1nep2fi2(x)])1/2∑i=1nep2fi2(x))-1=(-(∑i=1nep2fi2(x))(ln[∑j=1i≠jnep2fj2(x)-p2+1]))×(p2(ln[∑i=1nep2fi2(x)])1/2∑i=1nep2fi2(x))-1<0.
Since ∂G(x,p)/∂p<0, then the function
(52)G(fi(x),p)=1p[ln(∑i=1nep2fi2(x))]1/2
is descent for variable p.
LetG(f-(x),p)=(1/p)ln[nep2f-2(x)]1/2, then limp→+∞G(f-(x),p)=f-(x).
For an increasing ordered series of numbers(53)0<p1<p2<⋯<pk<⋯<+∞,
we have that
(54)G(f-(x),p1)≥G(fi(x),p1)≥G(fi(x),p2)≥⋯≥G(fi(x),pk)≥⋯≥limp→+∞G(fi(x),p)=max1≤i≤nfi(x),(55)minx∈RmG(f-(x),p1)≥minx∈RmG(fi(x),p1)≥⋯≥minx∈RmG(x,pk)≥⋯≥minx∈Rmlimp→+∞G(fi(x),p).
Taking the limit in (55) as p1→+∞, 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 increasep until the numerical overflow in order to determine the range of p, so that the linear and parabolic aggregation function approximates well the maximum envelope function maxi,jfij(x).Example A.
(56)
Find
x
∈
R
s
.
t
.
min
x
max
1
≤
i
≤
3
f
i
(
x
)
,
where
f
1
(
x
)
=
1
x
,
f
2
(
x
)
=
x
2
,
f
3
(
x
)
=
6
-
x
.
It can be seen from Figure5 that the optimal solution is x*=2 and minxmax1≤i≤3fi(x*)=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 Table1, where the linear aggregatation function is G(fi(x),p)=(1/p)[ln(∑i=1nepfi(x))] and the parabolic aggregatation function is G(fi(x),p)=(1/p)[ln(∑i=1nep2fi2(x))]1/2.Table 1
Iterative history of two aggregation functions in A.
The parameterp of the linear aggregation function
2
20
50
60
100
≥150
x
* the optimal solution
1.8647
1.9861
1.9944
1.9954
1.9973
overflow
G
(
x
*
)
4.1261
4.0251
4.0100
4.0083
4.0050
overflow
The parameterp of the parabolic aggregation function
1
2
3
4
5
≥6
x
* the optimal solution
1.9662
1.9914
1.9962
1.9979
1.9986
overflow
G
(
x
*
)
4.0629
4.0157
4.0070
4.0039
4.0025
overflowFigure5 and Table 1 show that the number of iterations is 17 and the parameter p≥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 p≥3, and the satisfactory solution is available. But the parabolic aggregate function can better approximate to the maximum envelope function maxi,jfij(x) 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, letfw(t)=tαw=t, fK(t)=fE(t)=fσ(t)=tαK=tαE=tασ=t3.44, and 0<tmin≤ti≤1.The element strain energy satisfies the following relationship as(57)eij=pijTKi-1pij2,
where Ki and pij are the variable elemental stiffness and the elemental node force vector, respectively.Substituting the weighting filter functionKi=fK(ti)Ki0 in (13) into (57), then we have that
(58)eij=pijT(Ki0)-1pij2fK(ti)=pijT(Ki0)-1pij2tiαK.Assume in thekth structural analysis that we have that
(59)eij(k)=pijT(Ki0)-1pij2fK(ti(k))=pijT(Ki0)-1pij2(ti(k))αK;
then in the (k+1)th structural analysis, and we have that
(60)eij=pijT(Ki0)-1pij2fK(ti)=pijT(Ki0)-1pij2fK(ti(k))·fK(ti(k))fK(ti)=eij(k)·fK(ti(k))fK(ti)=eij(k)(ti(k)ti)αK.In the context, we define the active set asIa={i∣0<tmin<ti<1,i=1,…,n}, the lower passive set as Ia_={i∣ti=tmin,i=1,…,n}, and the upper passive set as Ia¯={i∣ti=1,i=1,…,n}. Substituting (59) and (60) into G(eij,p)=(1/p)[ln(∑i=1nep2eij2)]1/2=G- in (41), then in the (k+1)th structural analysis we have that
(61)F(t)=∑i=1n[∑j=1Jep2[eij(k)·(fK(ti(k))/fK(ti))]2-Je(pe-i0(fσ2(ti)/fE(ti)))2]=0.Letxi=[fK(ti)]-2=[fE(ti)]-2=(ti)-2αK, aij(k)=p2[eij(k)fK(ti(k))]2=p2[eij(k)(ti(k))αK]2, and q=p2; then (58) is equivalent with
(62)F(t)=g(x)=∑i∈Ia_(∑j=1Jeq(eij(k))2-Jeq(e-i0·tmin2ασ-αE)2)+∑i∈Ia(∑j=1Jeqaij(k)xi-Jeq(e-i0)2·xi-(2ασ-αE)/αk)+∑i∈Ia¯(∑j=1Jeq(eij(k))2-Jeq(e-i0)2).In order to achieve the constraint in (59) explicitly, we demand the first-order partial derivatives of g(x)   on the design variables:
(63)∂g∂xi=∑j=1J[eqaij(k)xiqaij(k)]+Jeq(e-i0)2·xi-(2ασ-αE)/αkq(e-i0)2×[(2ασ-αE)αK]xi(-(2ασ-αE)-αK)/αK.Then, the first-order Taylor expansion formula ofg(x) is
(64)g(x)=g(x(k))+∑i∈Ia[∂g(x(k))∂xi](xi-xi(k)).
Substituting the intermediate variable xi=ti-2αk into (64), then we have that
(65)g(x)=b0(k)+∑i∈Iabi(k)ti-2αk,
where
(66)b0(k)=g(x(k))-∑i∈Ia[∂g(x(k))∂xi]xi(k),bi(k)=∂g(x(k))∂xi.Hence, the approximate model of (41) can be expressed as follows.Find topology variablest∈[0,1]n;
(67)minimizeW(t)=∑i=1nfW(ti)Wi0=∑i∈Ia_wi0tminαw+∑i∈Iawi0tiαw+∑i∈Ia¯wi0,subjecttob0(k)+∑i∈Iabi(k)ti-2αk=0,0<tmin≤ti≤1,i∈In={1,…,n},j∈Jm={1,…,J}.Setting aside consideration of interval constrainttmin≤ti≤1, then the augmented Lagrange function of (67) can be expressed as
(68)L(t,λ)=∑i∈Ia_wi0tminαw+∑i∈Iawi0tiαw+∑i∈Ia¯wi0+λ(b0(k)+∑i∈Iabi(k)ti-2αK).Taking extreme conditions for the saddle point in (68),
(69)∂L∂tm=αwwm0tmαw-1-2λαKbm(k)tm-2αK-1=0,(m∈Ia),(70)∂L∂λ=b0(k)+∑i∈Iabi(k)ti-2αK=0.By (69) we can have that
(71)tm=[2λ·(αKbm(k))(αwwm0)]1/(2αK+αw),(m∈Ia).Substituting (71) into (70), we can have that
(72)λ=[-b0(k)(∑i∈Iabi(k)((2αKbi(k))/(αwwi0))-2αK/(2αK+αw))]-(2αK+αw)/2αK.Substituting (72) into (71), we can have that
(73)t~m=[-b0(k)(∑i∈Iabi(k)((2αKbi(k))/(αwwi0))-2αK/(2αK+αw))]-1/2αK×[(2αKbm(k))(αwwm0)]1/(2αK+αw).With consideration on interval constrainttmin≤ti≤1(i=1,…,n),  we can have that
(74)tm(k+1)={tmin,(t~m≤tmin),t~m,(tmin<t~m<1),1,(t~m≥1).Record and update the active setIa and the passive set Ia_ and Ia¯; then return to (73) to calculate t~m, this process is repeated until the same set of active and passive sets does not change and meets the convergence criteria for iteration:
(75)|W(k+1)-W(k)|W(k)≤ε,
where W(k) and W(k+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 Figure6(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/m3, the Poisson’s radio is ν=0.3, and the elastic module of material is E=3×105 N/m2. In order to avoid stress concentration, the load P=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) Geometry and load condition for the short cantilever beam problem. (b) The optimal topology structure for the short cantilever beam problem.
(a)
(b)(a) The elemental strain energye distribution for the short cantilever beam problem. (b) The independent continuous topological design variables t-distribution for the short cantilever beam problem. (c) Stress σ distribution for the short cantilever beam problem.
(a)
(b)
(c)Thus, the optimal topology structure is shown in Figure6(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 p=2.Table 2
The distribution of the independent continuous topological design variables in the optimal topology structure of the short cantilever beam.
Value of design variables
t
=
t
min
t
min
<
t
<
0.1
0.1
≤
t
<
0.2
0.2
≤
t
<
0.3
0.3
≤
t
<
0.4
0.4
≤
t
<
0.5
Number of the elements
817
3
38
54
54
55
Value of design variables
0.5
≤
t
<
0.6
0.6
≤
t
<
0.7
0.7
≤
t
<
0.8
0.8
≤
t
<
0.9
0.9
≤
t
<
1.0
t
=
1.0
Numberof the elements
75
80
12
0
0
12Figure 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 Figure11 (cf. [63]), in which material constants are E=2.1e8 pa, ν=0.3, and ρ=1.0 kg/m3, where a vertically force P=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 Figure12, 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 p=2.Table 3
The distribution of the independent continuous topological design variables in the optimal topology structure of the MBB-beam.
Value of design variables
t
=
t
min
t
min
<
t
<
0.1
0.1
≤
t
<
0.2
0.2
≤
t
<
0.3
0.3
≤
t
<
0.4
0.4
≤
t
<
0.5
Number of the elements
2850
242
504
502
782
260
Value of design variables
0.5
≤
t
<
0.6
0.6
≤
t
<
0.7
0.7
≤
t
<
0.8
0.8
≤
t
<
0.9
0.9
≤
t
<
1.0
t
=
1.0
Number of the elements
138
30
0
0
0
92Figure 12
The optimal topology structure for the MBB-beam problem.(a) The elemental strain energye distribution for the MBB-beam problem. (b) The independent continuous topological design variables t-distribution for the MBB-beam problem. (c) Stress σ distribution for the MBB-beam problem.
(a)
(b)
(c)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 Figure17 (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 P. The material constants are E=21000 pa, ν=0.3, ρ=7800 kg/m3, and l=6.0 m,  and the uniformly distributed load is P=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 Figure18, 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 p=10.Table 4
The distribution of the independent continuous topological design variables in the optimal topology structure of the L-shape beam.
Value of design variables
t
=
t
min
t
min
<
t
<
0.1
0.1
≤
t
<
0.2
0.2
≤
t
<
0.3
0.3
≤
t
<
0.4
0.4
≤
t
<
0.5
Number of the elements
386
13
71
55
71
83
Value of design variables
0.5
≤
t
<
0.6
0.6
≤
t
<
0.7
0.7
≤
t
<
0.8
0.8
≤
t
<
0.9
0.9
≤
t
<
1.0
t
=
1.0
Numberof the elements
113
67
19
2
0
92Figure 18
The optimal topology structure for the L-shape beam problem.(a) The elemental strain energye distribution for the L-shape beam problem. (b) The independent continuous topological design variables t-distribution for the L-shape beam problem. (c) Stress σ distribution for the L-shape beam problem.
(a)
(b)
(c)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.


## 6.1. The Short Cantilever Beam

The first example is the short cantilever beam shown in Figure6(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/m3, the Poisson’s radio is ν=0.3, and the elastic module of material is E=3×105 N/m2. In order to avoid stress concentration, the load P=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) Geometry and load condition for the short cantilever beam problem. (b) The optimal topology structure for the short cantilever beam problem.
(a)
(b)(a) The elemental strain energye distribution for the short cantilever beam problem. (b) The independent continuous topological design variables t-distribution for the short cantilever beam problem. (c) Stress σ distribution for the short cantilever beam problem.
(a)
(b)
(c)Thus, the optimal topology structure is shown in Figure6(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 p=2.Table 2
The distribution of the independent continuous topological design variables in the optimal topology structure of the short cantilever beam.
Value of design variables
t
=
t
min
t
min
<
t
<
0.1
0.1
≤
t
<
0.2
0.2
≤
t
<
0.3
0.3
≤
t
<
0.4
0.4
≤
t
<
0.5
Number of the elements
817
3
38
54
54
55
Value of design variables
0.5
≤
t
<
0.6
0.6
≤
t
<
0.7
0.7
≤
t
<
0.8
0.8
≤
t
<
0.9
0.9
≤
t
<
1.0
t
=
1.0
Numberof the elements
75
80
12
0
0
12Figure 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 Figure11 (cf. [63]), in which material constants are E=2.1e8 pa, ν=0.3, and ρ=1.0 kg/m3, where a vertically force P=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 Figure12, 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 p=2.Table 3
The distribution of the independent continuous topological design variables in the optimal topology structure of the MBB-beam.
Value of design variables
t
=
t
min
t
min
<
t
<
0.1
0.1
≤
t
<
0.2
0.2
≤
t
<
0.3
0.3
≤
t
<
0.4
0.4
≤
t
<
0.5
Number of the elements
2850
242
504
502
782
260
Value of design variables
0.5
≤
t
<
0.6
0.6
≤
t
<
0.7
0.7
≤
t
<
0.8
0.8
≤
t
<
0.9
0.9
≤
t
<
1.0
t
=
1.0
Number of the elements
138
30
0
0
0
92Figure 12
The optimal topology structure for the MBB-beam problem.(a) The elemental strain energye distribution for the MBB-beam problem. (b) The independent continuous topological design variables t-distribution for the MBB-beam problem. (c) Stress σ distribution for the MBB-beam problem.
(a)
(b)
(c)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 Figure17 (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 P. The material constants are E=21000 pa, ν=0.3, ρ=7800 kg/m3, and l=6.0 m,  and the uniformly distributed load is P=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 Figure18, 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 p=10.Table 4
The distribution of the independent continuous topological design variables in the optimal topology structure of the L-shape beam.
Value of design variables
t
=
t
min
t
min
<
t
<
0.1
0.1
≤
t
<
0.2
0.2
≤
t
<
0.3
0.3
≤
t
<
0.4
0.4
≤
t
<
0.5
Number of the elements
386
13
71
55
71
83
Value of design variables
0.5
≤
t
<
0.6
0.6
≤
t
<
0.7
0.7
≤
t
<
0.8
0.8
≤
t
<
0.9
0.9
≤
t
<
1.0
t
=
1.0
Numberof the elements
113
67
19
2
0
92Figure 18
The optimal topology structure for the L-shape beam problem.(a) The elemental strain energye distribution for the L-shape beam problem. (b) The independent continuous topological design variables t-distribution for the L-shape beam problem. (c) Stress σ distribution for the L-shape beam problem.
(a)
(b)
(c)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 asp-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.



---

*Source: 735063-2013-12-04.xml*