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

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-ContinuousMapping (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.


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][2][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][5][6][7][8][9][10][11][12][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][15][16] with applications into many areas such as fluids, multiphysics and multidisciplinary applications, and much more [1,[17][18][19][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 Mathematical Problems in Engineering 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 meshdependency.
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][31][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][36][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.

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

Domain Setting and the Independent Topological Design
Variables.Let us consider a fixed design domain Ω in R 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 () is an indicator function that defines the distribution of material [1,6] 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  = ( 1 ,  2 , . . .,   )  ∈ {0, 1}  .Now, we consider the general binary discrete optimization problem with some property constraints as follows. Find where Φ = Φ() = ∑  =1   is the objective function,   is the physical or geometric parameter of the th element used in the finite analysis,   ( 1 , . . .,   ) and   are the th constraint response and its upper bound, respectively,  is the total number of constraints,   is the th independent topological design variable,  = ( 1 , . . .,   )  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], crosssectional area [45], or pseudodensity as the design variables, that is, 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 0 1 1 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) where   is the physical or geometric parameter of the th element used in the finite analysis,  0  is the original and true physical or geometric parameter of the material used of the th element.
As an appropriate relaxation,   = (  / 0  ) is allowed to take all possible values between 0 and 1, and it is clear that   = (  / 0  ) is an ideal independent topological design variable which has clear physical meaning.However, the stepup function (5) is discontinuous in the point   / 0  = 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   = (  / 0  ) 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   / 0  .

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 (  / 0  ) to approximate the step-up function (5) such that   = (  / 0  ) becomes a continuous and differentiable function (see Figure 2); that is, where the following properties are satisfied.
The polish function and independent continuous topological design variable.
(5) Approximation: Based on the concepts of the step-up function and polish function,   = (  / 0  ) can take values in the entire range [0, 1] instead of the discrete set {0, 1} and make the functional relationship   = (  / 0  ) 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   = (  / 0  ) 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   = (  / 0  ) take intermediate values between 0-1,   = (  / 0  ) not only reflect the proximity to corresponding elements "have" and "no", but also identify the size of the various physical values   in the elements.

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: that is, Moreover, we define the inverse function of the polish function (  ) =   / 0  =  −1 (  ) as the weighting filter function, where 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): (a) power function  =   ,  = 3, 10.32; (b) composite exponential function (c) modified sigmoid function 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 (  ) =   / 0  =  −1 (  ) plays a very important role as follows.
(a) Let the inverse mapping of the step-up function in (5) denote as   / 0  =  −1 (  ), 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   = (  / 0  ) take intermediate values between 0-1.To overcome this difficulty, we introduce the concept of weighting filter function (  ) =   / 0  =  −1 (  ) which continuously approximates the inverse mapping of the step-up function   / 0  =  −1 (  ), and it points out the proximity degree that   closes to  0  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 where  0  and  0  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, where   ,   , V  ,   , and   are, respectively, the variable allowable stress, mass, volume, area, and thickness of the th element for structural material, and  0  ,  0  , V 0  ,  0  , and  0  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, where   and   are the variable density and Young's modulus of the th element for structural material and  0  and  0  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   = (  / 0  ) 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.Let   be the objective function of the th element used in the finite analysis and  0  the original objective function of the th element for structural material, then we have that and the topology optimization model (3) could be approximated by the new continuous nonlinear topology optimization model; as follows. Find (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][48][49] or the relaxation approach [36,[50][51][52][53][54][55][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)   =   (  )  0  (23) such that where  0  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 : In the following, we prove that ( 23) and ( 24) can satisfy the -relaxation conditions.

Mathematical Problems in Engineering
Without loss of generality, assuming the stress constraints in the continuum structure design [36,[51][52][53][54][55][56] as then we have that 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 that is, 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.

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.

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  0   Φ (  ) may be defined as Then, the original topology optimization problem with stress constraints can be formulated as follows.
Find topological variables t ∈ [0, 1]  ; minimize where   is the th independent continuous topological design variable following (8),  = ( 1 , . . .,   ) 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.

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 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 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 where   ,   , and V  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 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   ≤  0  are satisfied if the formulation (37) is established.
Moreover, substituting (14) and (20) into   =   in (37), we have that where  0  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   ≤  0  in (32) as following.(41)

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][59][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 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.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  = as the variable global allowable strain energy.
For an increasing ordered series of numbers we have that min Taking the limit in (55) as  1 → +∞, 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 ,   ().( It can be seen from Figure 5 that the optimal solution is  * = 2 and min  max 1≤≤3   ( * ) = 4 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 . 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.

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).
The element strain energy satisfies the following relationship as where   and   are the variable elemental stiffness and the elemental node force vector, respectively.Substituting the weighting filter function   =   (  ) 0  in ( 13) into ( 57), then we have that Assume in the th structural analysis that we have that then in the ( + 1)th structural analysis, and we have that 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) =  in (41), then in the ( + 1)th structural analysis we have that Let , and  =  2 ; then (58) is equivalent with In order to achieve the constraint in ( 59) explicitly, we demand the first-order partial derivatives of () on the design variables: Then, the first-order Taylor expansion formula of () is Substituting the intermediate variable where Hence, the approximate model of ( 41) can be expressed as follows.
Record and update the active set   and the passive set   and   ; then return to (73) to calculate t , this process is repeated until the same set of active and passive sets does not change and meets the convergence criteria for iteration:       (+1) −  ()       ()  ≤ , where  () and  (+1) are the total weight of the structure of previous and current iterations, respectively, and  = 0.001 is the convergence precision.Weight (kg) Figure 9: The historical curves of the structural strain energy and the structural weight for the short cantilever beam problem.

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 3 , the Poisson's radio is ] = 0.3, and the elastic module of material is  = 3 × 10 5 N/m 2 .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.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.    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.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.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.9kg, in which the parameter of the parabolic aggregatation function is  = √ 10.

Figure 1 :
Figure 1: The step-up function and traditional topology design variable.
(d) The weighting filter function can formulate the explicit continuous nonlinear topology optimization model and obtain the final topological graphs.

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

Figure 5 :
Figure 5: Three test functions in Example A.

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

Figure 8 :
Figure 8: The historical curves of the elemental maximum stress and the structural weight for the short cantilever beam problem.0 10 20 30 40 50 60 70

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

Figure 11 :
Figure 11: Geometry and load condition for the MBB-beam problem.

Figure 12 :Figure 13 :
Figure 12: The optimal topology structure for the MBB-beam problem.

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

6. 2 .
The MBB-Beam.The second example is the MBB-beam shown in Figure11(cf.[63]), in which material constants are  = 2.18 pa, ] = 0.3, and  = 1.0 kg/m 3 , 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 16 : 3 N/m l/ 2 Figure 17 :
Figure 16: The historical curves of the elemental maximum strain energy and the structural weight for the MBB-beam problem.

Figure 18 :
Figure 18: The optimal topology structure for the L-shape 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 .The material constants are  = 21000 pa, ] = 0.3,  = 7800 kg/m 3 , and  = 6.0 m, and the uniformly distributed load is  = 0.3N.

Figure 19 :
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.
This paper develops three novel methods for solving stress-based topology optimization problems, which are Mathematical Problems in Engineering

Figure 20 :Figure 21 :Figure 22 :
Figure 20: The historical curves of the elemental maximum stress and the structural weight for the L-shape beam problem.

Table 1 :
Iterative history of two aggregation functions in A.

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

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