Continuum structural topological optimizations with stress constraints based on an active constraint technique

Stress‐related problems have not been given the same attention as the minimum compliance topological optimization problem in the literature. Continuum structural topological optimization with stress constraints is of wide engineering application prospect, in which there still are many problems to solve, such as the stress concentration, an equivalent approximate optimization model and etc. A new and effective topological optimization method of continuum structures with the stress constraints and the objective function being the structural volume has been presented in this paper. To solve the stress concentration issue, an approximate stress gradient evaluation for any element is introduced, and a total aggregation normalized stress gradient constraint is constructed for the optimized structure under the r−th load case. To obtain stable convergent series solutions and enhance the control on the stress level, two p‐norm global stress constraint functions with different indexes are adopted, and some weighting p‐norm global stress constraint functions are introduced for any load case. And an equivalent topological optimization model with reduced stress constraints is constructed,being incorporated with the rational approximation for material properties, an active constraint technique, a trust region scheme, and an effective local stress approach like the qp approach to resolve the stress singularity phenomenon. Hence, a set of stress quadratic explicit approximations are constructed, based on stress sensitivities and the method of moving asymptotes. A set of algorithm for the one level optimization problem with artificial variables and many possible non‐active design variables is proposed by adopting an inequality constrained nonlinear programming method with simple trust regions, based on the primal‐dual theory, in which the non‐smooth expressions of the design variable solutions are reformulated as smoothing functions of the Lagrange multipliers by using a novel smoothing function. Finally, a two‐level optimization design scheme with active constraint technique, i.e. varied constraint limits, is proposed to deal with the aggregation constraints that always are of loose constraint (non active constraint) features in the conventional structural optimization method. A novel structural topological optimization method with stress constraints and its algorithm are formed, and examples are provided to demonstrate that the proposed method is feasible and very effective. © 2016 The Authors. International Journal for Numerical Methods in Engineering published by John Wiley & Sons Ltd.


SUMMARY
Stress-related problems have not been given the same attention as the minimum compliance topological optimization problem in the literature. Continuum structural topological optimization with stress constraints is of wide engineering application prospect, in which there still are many problems to solve, such as the stress concentration, an equivalent approximate optimization model and etc. A new and effective topological optimization method of continuum structures with the stress constraints and the objective function being the structural volume has been presented in this paper. To solve the stress concentration issue, an approximate stress gradient evaluation for any element is introduced, and a total aggregation normalized stress gradient constraint is constructed for the optimized structure under the r th load case. To obtain stable convergent series solutions and enhance the control on the stress level, two p-norm global stress constraint functions with different indexes are adopted, and some weighting p-norm global stress constraint functions are introduced for any load case. And an equivalent topological optimization model with reduced stress constraints is constructed,being incorporated with the rational approximation for material properties, an active constraint technique, a trust region scheme, and an effective local stress approach like the qp approach to resolve the stress singularity phenomenon. Hence, a set of stress quadratic explicit approximations are constructed, based on stress sensitivities and the method of moving asymptotes. A set of algorithm for the one level optimization problem with artificial variables and many possible non-active design variables is proposed by adopting an inequality constrained nonlinear programming method with simple trust regions, based on the primal-dual theory, in which the non-smooth expressions of the design variable solutions are reformulated as smoothing functions of the Lagrange multipliers by using a novel smoothing function. Finally, a twolevel optimization design scheme with active constraint technique, i.e. varied constraint limits, is proposed to deal with the aggregation constraints that always are of loose constraint (non active constraint) features in the conventional structural optimization method. A novel structural topological optimization method with stress constraints and its algorithm are formed, and examples are provided to demonstrate that the proposed method is feasible and very effective.

INTRODUCTION
Because most of the materials in real structural design problems may suffer from strength failure, a stress-based topology optimization method (STOM) has also attracted much interest. Earlier efforts on generating optimal configurations of continuum structures subjected to von Mises stress constraints over each individual element can be traced back to late 1990s, see, for example, Cheng and Zhang [1] and Shim and Manoochehri [2]. Since then, the topology optimization with local stress constraints has been studied by Duysinx and Sigmund [3], Duysinx and Bendsøe [4], Pereira et al. [5] and Bruggi [6], Jeong and Park [7], Moon and Yoon [8], James and Waisman [9] París and Colominas [10], Yoon [11], etc.
To the best of our knowledge, there are three theoretical issues related to the STOM, namely, the singularity issue, the element wisely defined stress constraint issue, and the highly nonlinear constraint issue. The singularity issue in the STOM indicates that the stress constraints defined at the centers of the finite elements of the void domain should be removed at every iteration of the optimization [6,[12][13][14][15][16]. Many relaxation methods have been proposed to resolve this issue [6,[12][13][14][15][16], in which the qp-relaxation method has been applied owning to its simplicity and effectiveness.
Furthermore, as the stress constraints are defined at every finite element in the topology optimization of stress constrained continuum structures, the number of stress constraints in the optimization problem becomes very large. Therefore, the large-scale nature of local stress constraints presents a major obstacle to the practical application of stress constrained topology optimization. For tackling this difficulty, two main strategies, namely, the 'active-set method' and the 'global constraint method', for reducing the number of stress constraints have been adopted in existing studies. In the 'active-set method', only the potentially critical stress constraints are taken into account at each iteration step [17]. By discarding a fraction of the stress constraints that are far from active, this method saves the computational cost to some extent. However, the selected active constraint set at the end of the optimization usually vary from 10-40% of the total number of constraints and may still be large when a fine discretization is used. In this sense, the 'active-set method' does not bring the required efficiency. By coupling a global compliance enforcement with local stress constraints, Bruggi and Duysinx [18] proposed a selection method leading to further reduction in the number of active constraints. The 'global constraint method' is more efficient because it aggregates all the local stress constraints into one Kreisselmeier-Steihauser or p-norm global constraint function. Many authors, such as Duysinx and Sigmund [3],Yang and Chen [13], Qiu and Li [19], and París et al. [20], have employed such a strategy. The main disadvantages of the global constraint method are the high nonlinearity of the global function and a weaker control on the stress level. (Generally speaking, the global function constraints often are non-active constraints during the optimization procedure.) The high nonlinearity may cause unstable convergence during the optimization procedure, and the weaker control on the stress level may lead to a serious violation of some local stress constraints. The technique of 'block aggregation' [21] by using more global constraints to reduce the number of constraints involved in a single global function would be helpful in avoiding these disadvantages, but not completely. The adaptive normalization scheme proposed by Le et al. [16] works well in controlling the maximum stress. However, it is difficult to obtain the exact optimal solution because this scheme may prevent the local stresses of some elements to reach the permissible strength limit. Recently, some remedies [22][23][24] addressed the highly nonlinear behavior including stress concentration issue of the aggregation stress constraints.
In the STOM, the highly nonlinear issue of the aggregation stress constraints, such as these aggregation stress constraint approximation problem with certain accuracy and a large design variable region, and an equivalent approximate optimization model issue, such as the control problem of the maximum several von-Mises stresses whose quantities and locations change with an increase in iteration step, remain problematic and require further mathematical and scientific modifications and contributions.
To resolve these issues, we proposed a new and effective topological optimization method of the continuum structure with the objective function being the structural volume and the constraint functions being the stresses, by constructing several active aggregation constraint functions, and building a set of algorithm for the one level optimization problem with artificial variables and a lot of possible non-active design variables, based on the aforementioned related researches, the method of moving asymptotes (MMA), and various numerical approaches. This paper is organized as follows. In Section 2, we describe element topological variables, interpolation functions, a stress relaxation approach, a basic stress-based topology optimization model, and our development of an equivalent topology optimization approximate model with stress aggregation gradient constraints and stress aggregation constraints. In Section 3, we describe one order derivative formula of the aggregation stress gradient functions and smooth aggregation stress functions, and their explicit approximations. In Sections 4 and 5, an optimization strategy, design space adjustments, and an optimization solving numerical approach are presented. And in Section 5, two numerical examples are presented to show the advantages of the proposed method. In Section 6, a conclusion is given.

Element topological variables and interpolation functions.
The discrete topological variable e that takes value 0 or 1, is replaced by the continuous topological variable e that varies between min e and 1. Consequently, the difficulty arising in the discrete optimization is avoided by mathematical programming optimization algorithms. The volume and the elastic modulus of an element are recognized by their interpolation functions: where V e and E e are the volume and the elastic modulus of the e-th element, respectively. And V e 0 and E e 0 are the original volume and the original elastic modulus of the e-th element, respectively. e is the topology variable of the e-th element. By referring to the rational approximation for material properties [25], the interpolation functions g v . e / and g k . e / are selected as follows, where˛v and are two empirical parameters, D 4, and˛v D 1.25 are used in this paper examples.

The qp stress relaxation approach.
It is assumed that the von Mises stress e; r vm of the e th element for the r th load case is approximately measured by the von Mises stress of the geometrical center point of the e th element for the r th load case. Then, the stress vector ¢ e; r of the e th three dimensional element based on its nodal displacement vector u e; r for the r th load case can be approximately expressed as ¢ e;r D ® e;r xx ; e;r yy ; e;ŕ´; e;r xy ; e;r x´; e;r y´¯T D D e B e . r u e / D N S e . r u e / D g k . e / N S 0 e . r u e / (3) where D e is the conventional elastic modulus matrix of the e th element and B e is the conventional geometrical matrix at the geometrical center point of the e th element. N S e and N S 0 e denote the stress matrix and the original stress matrix at the geometrical center point of the e-th element, respectively.
To eliminate the arising of the stress singularity phenomenon, referring to the qp approach idea OE6 , here, Equation (4) for stress constraint relaxations is adopted: where g q . e / D q e =.1 C .1 e //; q < 1 where e Y is the allowable von Mises stress limit of the e-th element for the full density material (failure stress), q is an empirical parameter, and q D 0:7 is used in this paper examples. For the convenience of consequent formula expressions, it is assumed that S e l .l D 1; 2; : : : ; 6/ (for three Figure 1. Element categories in the fixed finite mesh domain: (a) a topology structure only including removed element (white areas) and retained elements (black areas); (b) the transferring topology structure with retained elements (black areas), removed elements (white areas), and artificial material elements (reddish areas), corresponding to Figure 1(a).
dimension structures) denotes the l-th row of g k . e /=.g q . e / e Y / N S 0 e , and Q S e l .l D 1; 2; : : : ; 6/ denotes the l-th row of .1 q/. e / . q/ N S 0 e = e Y .

A basic stress-based topology optimization model.
In this study, a fixed maximum finite element mesh is adopted. All structural characteristic data during the structural optimization procedure are stored and used, based on the fixed maximum finite element mesh. When the structural topology obtained by the proposed procedure gets close to the vicinity of the optimum topology, many elements with small topology variable values are removed, and a smaller size structure within the fixed maximum finite element mesh can be obtained (Figure 1(a)). When a structure being similar to Figure 1(a) occurs at some optimization iteration step, it must be assured that the structure (black area in Figure 1(a)) be non-singular, and some error removed elements can be restored. Here, a layer of new elements is inserted around the cavities and boundaries of the current structure during the structural optimization procedure (Figure 1(b)). And these element topological variables are all set to m j .j D 1; 2; : : : ; m h / Á N 1 m D 0:0001 at the first iteration step after that many elements with small topology variable values are removed. Here, they are called as artificial material elements; m h is the number of the inserted artificial material elements. These new elements are incorporated in structural analyses and optimization designs at consequent several iteration steps so that design space expansions and reductions can be automatically carried out.
Structural elements can be categorized into two kinds: (1) designable elements and (2) nondesignable elements. The non-designable element numbering may be expressed as i p .p D 1; 2; : : : ; P /, its topological variables i p D 1 .p D 1; 2; : : : ; P / do not change during the optimization process, and P is the number of the non-designable elements. Q is the number of the designable elements, whose numbering may be represented by n q .q D 1; 2; : : : ; Q/. The designable elements include the artificial material elements and the designable and retained elements (whose topology variable values are between 0 and 1, and whose numbering is expressed as h i .i D 1; 2; : : : ; R/, and R is the number of the designable elements and retained elements).
The structural topology optimization problem in terms of a minimum volume approach with stress constraints can be stated as Equation ( where the lower limit min n v D 0:00001 of the n v -th element topological variable n v is set so that the structure optimized is always kept to be non-singular during the optimization process. And s i ; r vm denotes the von Mises stress of the s i th element for the r th load case, V denotes the optimized structural volume, and L is the number of the load cases applied on the structure.

An approximation expressions of stress gradients
Generally speaking, when there is a small change of the design variables of the neighborhood around an element, there is a dramatic change of the stress quantity of the element, and this phenomenon is called as the stress concentration. The stress concentration factors of some typical structures, such as components with U-notches and V-notches, are evaluated in structural engineering. There also are the highly nonlinear behavior of the design variables in the aggregation functions, such as the p-norm and the segregated global stress norms [13,[20][21][22][23][24], except this nonlinear behavior caused by the stress concentrations.
To solve the stress concentration issue, an approximate stress gradient evaluation for any element is introduced for any load case. Here, the approximation stress gradient absoluteness i;r T of the i th element for the r th load case is defined as where Y i D ¹e; kx i x e k 6 r 0 ; i ¤ eº denotes an element set related to the stress gradient of the i-th element and r 0 is a pre-defined radius parameter for a small element neighborhood. n i T represents the number of elements included in the set Y i , and x i and x e are the coordinate vectors of the geometrical center points for the i-th element and the e-th element, respectively. And r i e D kx i x e k denotes the distance between two geometrical center points for the i-th element and the e-th element. r 0 D 1:1d mesh has been used in the example problems of this paper, and d mesh is the minimum space distance between the finite element grid lines.
Furthermore, by referring to the stress aggregation functions [19][20][21], a total aggregation normalized stress gradient function is constructed for the optimized structure under the r th load case and is used to approximately form an approximate maximum normalized stress gradient constraint of the structural elements under the r th load case. A total aggregation normalized stress gradient function under the r th load case is expressed by Equation (9).
where one stress gradient aggregation parameter q sg D 5 is adopted in the examples of this paper.

An equivalent topology optimization approximate model with stress constraints
To dramatically reduce the number of stress constraints, and construct an equivalent stress-based structural topology optimization model with only several constraints, here, both aggregation stress constraints are constructed, and aggregation stress gradient constraints are introduced to replace all the stress constraints in Equation (6). It is well known that a threshold parameter for any material stress gradient has not been investigated in engineering structures, and when structural material quantity decreases, the total aggregation stress gradient function of the structure may increase for same load case. Based on the aggregation stress gradient of the current structure, a varied upper limit with an increase in the iteration number is adopted to construct an aggregation stress gradient constraint so that there will not be a big change of the maximum element stress gradient in the structure at an iterative step.
To overcome the unstable convergence caused by the high nonlinearity of the 'global constraint method' during the optimization procedure, two p-norm global constraint functions with different indexes for any load case are adopted. To enhance the control on the stress level, some weighting p-norm global constraint functions for any load case are introduced. Referring to Refs. [17,[19][20][21], the optimization model (6) at the k-th iterative step is temporarily transferred into the approximate model (10).
where q sg , q c ; p˛.c D 1; 2I˛D 1; 2; : : : n w / and n w are empirical parameters and n w D 2; q sg D 5; q 1 D 8; q 2 D 16; p 1 D p 2 D 16 are used in this paper examples. r T indicates an introduced aggregation normalized stress gradient limit for the r th load case at the k th outer loop iteration step. N k;L n v D max.. where Sn D 1 denotes the first optimization adjustment phase given in the consequent content and Sn D 2 denotes the optimization acceleration phase given in the consequent content.This trust scheme can reduce or avoid the big fluctuations of the maximum von Mises stress of the optimized structure and assures certain accuracy of the approximate smooth aggregation functions in Equation (10). And weighting functions r;s i .x; y;´/ in Equation (10)  where n w in Equations (10) and (11b) represents the number of weighting aggregation stress constraint functions with weighting functions r;s i .x; y;´/ under the r th load case and W r;q s i .x; y;´/ is similar to a rectangular window function with a variable d r; q s i .x; y;´/ of the signal processing. In Equation (11c), d r; q s i .x; y;´/; q D 1; 2 ;˛is related to the geometrical center point of the element with a so-called q-th region maximum von Mises stress of the optimized structure under the r th load case. And it is required that the numbering of the center elements of small structural neighborhoods with the maximum to the n w -th region maximum von Mises stresses of the optimized structure under the r th load case, respectively, is determined. The so-called qth maximum von Mises stress of the optimized structure and its element numbering are defined as follows: (1) At first, setting U r 1 D ¹s 1 ; s 2 ; ; s P CR º, the numbering of the element possessing the maximum von Mises stress within the element set U r 1 of the optimized structure under the r th load case , is easily determined and recorded as s r 1 .
(2) Setting z.s r 1 / is the geometrical center point coordinate of the s r 1 -th element, a set U r 2 D ® s i j z.s i / z.s r 1 / > and s i 2 U r 1¯m ay be obtained, in which z.s i /is the geometrical center point coordinate of the s i -th element. Similarly, the numbering of the element possessing the maximum von Mises stress within the element set U r 2 of the optimized structure under the r th load case is determined and recorded as s r 2 . Here, the maximum von Mises stress within the element set U r 2 of the optimized structure under the r th load case is called as the so-called 2-th region maximum von Mises stress of the optimized structure under the r th load case.
(3) Similarly, U r q D°s i j z.s i / z.s r q 1 / > and s i 2 U r q 1 ± ; q D 1; 2; ; n w I r D 1; 2; L, may be determined, and the numbering of the element possessing the maximum von Mises stress within the element set U r q of the optimized structure under the r th load case is determined and recorded as s r q . Here, the maximum von Mises stress within the element set U r q of the optimized structure under the r th load case is called as the so-called q th region maximum von Mises stress of the optimized structure under the r th load case.
In Equation (11c), d r; q s i .x; y;´/; q D 1; 2 ;˛is the distance between the geometrical center point (.x; y;´/ denoting its coordinate ) of the s i -th element and the geometrical center point of the element possessing the q-th region maximum von Mises stress of the optimized structure under the r th load case; is a pre-defined empirical parameter approximately as a radius of a small structural neighborhood, which center point is the geometrical center point of the element possessing the q-th region maximum von Mises stress of the optimized structure under the r th load case. D 10d mesh is used in the first example of this paper.
Generally speaking, W r;q s i .x; y;´/ is correlated with the topology variables of the small structural neighborhood including the element with the q-th region maximum von Mises stress of the optimized structure under the r th load case during the optimization procedure. To simplify the model formula, the change relation of W r;q s i .x; y;´/ with topology variables is omitted here, and it is assumed that W r;q s i .x; y;´/ is not related to topology variables. Smooth aggregation stress functions f r;w ;˛D 1; 2; : : : ; n w in Equation (10), which obviously enhance the actions of the second region maximum von Mises stress, the third region maximum von Mises stress of the optimized structure under the r th load case, etc., can be constructed by using weighting functions r;s i .x; y;´/. It is found that weighting aggregation stress functions f r;w and aggregation stress functions f r;c so , respectively, possess completely different distributions within the structural domain. The equivalence of the model (10) and the model (6) can be enhanced by using these weighting aggregation stress constraints and aggregation stress gradient constraints in order to strongly control structural maximum stress levels during the optimization procedure.
In the optimization model (10), the second equation to the fourth equation may be derived by using some simple multiplications, additions, and powers of the second equation in the model (6). Because each structural element possesses a very different stress level for any load case, if the von-Mises stresses of only several structural elements excess the allowable von Mises stress limit, these aggregation constraints in the model (10) basically are non-active constraints (namely, loose constraints), even if at the optimum solution state of the model (6). When the model (10) is solved by adopting a mathematical programming method, its series solutions are mainly affected by the objective minimization requirements and the trust region, and it is very difficult to control the maximum von Mises stress of the optimized structure. Therefore, the model (10) with reduced constraints still is not equivalent to the model (6). To construct a dramatically reduced optimization model being approximately equivalent to the model (6), a heuristic operation is added to alter these smooth aggregation constraints so that at least one of these smooth aggregation constraints in the model (10) is an active constraint at each iterative step, and its approximate series solutions are closely related to the structural stress and stress sensitivity fields. Being referred to the varied constraint limit approach in Refs. [26][27][28], the original optimization model (6) is transferred into the approximate model (12), which replaces the model (10).

One order derivative formula of the aggregation stress gradient functions and smooth aggregation stress functions
In this paper, e; r vm =g q . e / is called as an effective von Mises stress of the e-th element for the r th load case; the derivative of e;r vm =.g q . e / e Y / of the e-th element with respect to n v is expressed as (for three dimensional structure) re six stress components of the e-th element for the r th load case. The stiffness K n v 0 is the stiffness matrix of the n v -th element when n v D 1. e;j u n v denotes an element displacement vector, which is related only to the n v th element of the structure because of the unit load F j at the j -th degree of freedom of the e-th element in the structure. And r u n q is an element displacement vector, which is related only to the n q -th element of the structure under the r-th load case. And S e l j denotes the component of S e l in Section 2.1.2, and Q S e p is referred to Section 2.1.2. The derivative of aggregation stress functions f r;c so . /; c D 1; 2 in the optimization model (10) with respect to n v at .k 1/ n v ; v D 1; 2; : : : ; Q may be expressed as Equation (15) where r N N u is obtained by solving the adjoint state Equation (16) and N N n v ; r vm;ad is given in Equation (17b). Here, where F s i ;j denotes a unit virtual load at the j -th degree of freedom of the s i -th element in the structure and r N N u n v is an element displacement vector, which is related only to the n v th element of the structure because of the new virtual load N N F r applied on the structure. The derivative of weighting aggregation stress functions f r;w . /;˛D 1; 2; : : : ; n w I r D 1; 2; : : : ; L in the optimization model (10) with respect to n v at .k 1/ n v ; v D 1; 2; : : : ; Q may be expressed as Equation (18) where r N N N u is obtained by solving the adjoint state Equation (19) and N N N n v ; r vm;ad is given in Equation (20b). Here, is an element displacement vector, which is related only to the n v th element of the structure in the solution r N N N u of the Equation (19) because of another new virtual load N N N F r applied on the structure.
The derivative of total aggregation stress gradient functions f r T C . /; r D 1; 2; : : : ; L in the optimization model (10) with respect to n v at .k 1/ n v ; v D 1; 2; : : : ; Q may be given as follows: where r N u is obtained by solving the adjoint state Equation (22). N n v ; r vm;1 and N n v ; r vm;2 are given in Equation (23b) and (23c), respectively. Here, where r N u n v is an element displacement vector, which is related only to the n v th element of the structure in the solution r N u of the Equation (22) because of the new virtual load N F r applied on the structure.

Explicit approximations of aggregation functions
Gradient-based optimization algorithms that rely on nonlinear but convex separable approximation functions have proven to be very effective for large-scale structural optimization [32]. Here, the MMA expression Q f r T C;M .¡/ of f r T C .¡/ with respect to topology variables y v D n v ; v D 1; 2; : : : ; Q at .k 1/ n v ; v D 1; 2; : : : ; Q is adopted and can be given as Equation (24) where where y max v D 1 and y min v D min n v are the upper limit and lower limit of y v , respectively. Here, are the upper limit and lower limit of moving asymptotes, and other notations are referred to the references [31][32][33].
Observing the conventional dual method of optimization problems with several constraints and many design variables, it is known that the quadratic Taylor series expansions of constraint functions and an objective function are very beneficial to obtaining the explicit closed solutions by using Karush-Kuhn-Tucker (KKT) conditions. Moreover, from Equations (7)-(9), it is found that the function f r T C . / is a serious strong nonlinear function of topology variable vector . To enlarge the design variable region of explicit approximations of aggregation stress gradient functions, reciprocal intervening variables x v D 1= n v ; v D 1; 2; : : : ; Q are selected as design variables, and a quadratic Taylor series expansion of f r T C at x .k/ is adopted. The one-order derivatives of f r T C with respect to design variables x v D 1= n v at .k 1/ n v ; v D 1; 2; : : : ; Q may be obtained by Equations (26) and (21).
The second order derivatives of f r T C ; r D 1; 2; : : : ; L with respect to design variables x v D 1= n v at .k 1/ n v ; v D 1; 2; : : : ; Q may be obtained by adopting Equation (24) and its convex approximation treatment, and its expression is given in Equation (27).
In this paper, a convex quadratic Taylor series expansion of f r T C . /; r D 1; 2; : : : ; L at x .k 1/ is approximately expressed as follows.
Similarly, a convex quadratic Taylor series expansion of f r;c so .¡/ ; r D 1; 2; : : : ; LI c D 1; 2 and f r;w .¡/; r D 1; 2; : : : ; LI˛D 1; 2; : : : ; n w at x .k 1/ , respectively, may approximately be expressed as follows  (29) for each load case. There are similar specialties for approximate functions Q f r;c so .x/ and Q f r;w .x/ in Equations (30) and (31). To reduce the errors of these approximate functions, an optimization strategy in the Section 4.4 of this paper will be given. In the optimization strategy, when the many outer optimization iterations are completed, many elements with small topology variables, which do not belong to the optimum structure, will be removed, and the design space reduction can be carried out by a phase transferring step.

A transferred approximate optimization model
Here, the objective function and constraint functions of the approximate model (12) are replaced by their quadratic Taylor series expansions. Following Svanberg's approach [30], artificial variables w D w 1 ; w 2 ; : : : ; w m T ; w > 0, are introduced into the optimization problem (12), so that the following enlarged problem is addressed. Even if there is an empty solution set for the optimization model (12), a solution with the minimum violating the constraints in Equation (12) can be obtained.
Many simulation examples administrate that a solution with the minimum violating the constraints in Equation (12) at many cases still belongs to the feasible solution set of the optimization model (6). Moreover, being referred to Ref. [31,32], the second order term of constraint functions is introduced into the objective function by using the Lagrange multiplier vector obtained at the previous outer loop iteration; the model (12) is transformed into an equivalent approximate quadratic programming model (34) ; c D 1; 2I r D 1; 2; : : : ; L (36c) In Equations (34) and 35, .k 1/ j is the j -th Lagrange multiplier obtained at the previous outer loop iteration. The constants c j must be chosen large enough so that the variables w j are zero at the optimal solution, in case the original problem has a non-empty feasible set and fulfills a constraint qualification. Being referred to Refs. [29][30][31][32], c j D 3000 and d j D 1; j D 1; 2; : : : ; m are selected in the examples of this paper.
In order to deal with the checkerboard problem of topology optimization, the approach in Ref. [33] is adopted to modify Q c j v and one-order and two-order sensitivities of the Q f .k/ obj .x/ for all real material elements. The modified values of Q c j v and and one-order and two- for all real material elements are substituted into Equations (34) and (35).

A smooth dual solving method for the approximate sub-problem
When there are many design variables and several constraints for an optimization problem, the dual solving method can dramatically reduce the computational quantity. However, because there are design variable small trust regions x k;L v 6 x v 6 x k;U v , there may be many non-active design variables in the model (34). If the conventional dual method in Refs. [26][27][28] is adopted to solve the Lagrange multipliers, the second order matrix of the dual programming may become an ill matrix, so that it is very difficult to obtain the Lagrange multipliers. Here, a smooth dual solving method is proposed to solve the approximate sub-problem.
The programming problem of Equation (34) . / and w j . j / be the design variables and artificial variables at the k th outer loop iteration step and the (lC1) th inner loop iteration step. Equation (40) may be obtained by using the KKT condition.
. / ± ± ; v D 1; 2; : : : ; Q where Here, x .k;lC1/ v . / and w j . j / in Equations (40) and (41b) are approximately transformed into the following smoothing functions N x .k;lC1/ v . / and N w j . j / by use of Â p .t / D p ln.1 C e t=p / and by choosing two small parameters p 1 and p 2 far less than 1.
where p 10 and p 20 are two empirical parameters and many numerical simulations show that a value between 0.03 and 0.3 can be selected as p 10 or p 20 :p 10 D 0:095; and p 20 D 0:094 are used in this paper examples.

A criterion of related parameter renewals in the outer iterations
The KKT conditions of the model (34) with varied aggregation constraint limits may approximately be satisfied when the design variables and the Lagrange multipliers obtained achieve at the optimum state for one set of varied displacement limits. To reduce iteration number, a relaxed criterion of related parameter renewals at the outer iterative steps is given. The maximum value ı .x .k/ ; .k/ / Á 2 ( k 2 ¹1; 2; : : : ; 5º ) at the beginning five iteration steps of the optimization process is selected as a reference value. It is considered that the optimum state of the optimization problem (40) has been achieved for one set of varied aggregation constraint limits, whenever x and verify one of Equations (49) and (50), and the varied aggregation constraint limits are renewed by adopting Equation (13).

The proposed optimization strategy and design space adjustments
When there are many elements with small topology variables in the optimized structure, the relative errors of approximate functions Q f r T C .x/, Q f r;c so .x/; Q f r;w .x/ in Equations (29)-(31) with respect to their original nonlinear functions f r T C . /; f r;c so .x/, and f r;w .x/, respectively, become bigger. To improve the accuracy of these approximate functions and efficiency of the structural optimization solving, the whole optimization iteration process is divided into a smooth optimization phase, a phase transferring step, and an accelerating optimization phase with heuristic operations. The solution obtained by the proposed procedure can approach to the vicinity of the optimum topology when the smooth optimization phase are completed. At the phase transferring step, many elements with small topology variables, which do not belong to the optimum structure, will be removed, and the design space reduction can be carried out.Then, at the accelerating optimization phase, heuristic operations will be given to improve the accuracy of the aforementioned approximate functions and make the design topology good black/white. Figure 2 depicts a flowchart of the proposed method. In Figure 2, the approaches of structural response and virtual load response analyses, varied aggregation constraint limit renewals, structural approximate model constructions, and topology variable renewals have been given in Sections 2 to 4 of this paper. Here, heuristic operations at these optimization phases are given as follows.
The beginning optimization iterations belong to the smooth optimization phase. The topology variables of some inefficient material elements can smoothly be changed by applying the smooth dual solving algorithm because of the action of the trust region, for the model (12) with varied aggregation constraint limits. There may be many elements with small topology variables in the optimal topology when the many outer optimization iterations are completed. At this case, there are certain relative errors of approximate aggregation functions with respect to their original nonlinear functions in the optimization model, and the objective reduction quantity becomes smaller with an increase in outer loop iterative number. A criterion for the smooth optimization phase (first phase) end is given as follows. 0:3 where, 0:3 mean denotes the average value of all the topology variables below 0.3. " 4 ; " 5 and " 6 are three prescribed small empirical parameters, k rd is a prescribed empirical parameter, and " 4 D 0:05, " 5 D 10 3 ; " 6 D 1: 10 3 , and k rd D 65 are adopted in the proposed examples.

heuristic operation at a phase transferring step.
To improve the accuracy of aggregation approximate functions in the Section 3.2, a phase transferring step is adopted, and the design space reduction is carried out. Its initial design is the near optimum design obtained when the smooth optimization phase ends. Firstly, after solving Equation (34) by adopting the smooth dual method, two element sets N N T d and N N T a can be obtained as follows If one of Equations (52) and (53) is satisfied, Equation (55) is carried out. Here, the second term in Equation (55) is used to enable that the maximum effective von-Mises stress of the optimized structure will not dramatically increase. Then, a layer of artificial material elements is inserted around the cavities and boundaries of the current structure at the beginning of the consequent iteration step (Section 2.1.2). These new elements are incorporated in structural analyses and optimization designs at consequent iteration steps so that some error removed elements can be restored. Moreover, the elements with topology variable min n q D 0:00001 are removed from the optimized structure as the structure shown in Figure 1(b), and the mesh information of this modified optimized structure, are renewed. At the beginning of the consequent iteration steps, the varied aggregation function constraint limits described earlier are changed according to the Equation 13, for the models (12).
When aforementioned operations are finished, the new topology variable vector is set as .kC1/ , the iteration step number is set as k C 1, and optimization iterations will enter the accelerating optimization phase with heuristic operations.

Heuristic operation at the accelerating optimization phase.
When the phase transferring step ends, the solution obtained by the proposed procedure can get close to the vicinity of the optimum topology. Another heuristic operation is then given to enable the design topology good black/white distribution after solving Equation (34) by adopting the smooth dual method. The heuristic operation is described as follows.
When the structural topology variables n q .q D 1; 2; : : : ; Q/ are renewed by adopting the smooth dual method, the removal threshold del th for topology variables are determined as follows.

The stopping iteration criterion for the minimum volume approach with stress constraints.
The outer loop finishes successfully whenever the following condition is satisfied.
where k st is the number of the outer routines after applying the phase transferring. k st and " 7 are two predefined small empirical parameters, and k st D 15 and " 7 D 0:001 are adopted in the proposed examples.

NUMERICAL EXAMPLES
In this section, two numerical examples are shown to demonstrate the effectiveness of the proposed approach for the solution of stress-related topology optimization problems.

A L-shape beam design
An L-shape beam shown in Figure 3 is used to examine the numerical performance of the proposed approach. The material, load, and geometry data are chosen as dimensionless because the main purpose is to test the numerical performance of the proposed algorithm. The thickness of the structure is 1. The material parameters are Young Module E D 1 and the Poisson ratio D 0.3 for solid material. A plane stress state is assumed and the four-node bilinear rectangular elements are used to interpolate the primary and adjoint displacements. The finite element mesh consists of 5184 same size quadrilateral elements. A vertical force of F D 1 is applied on the middle part of the right vertical edge and is distributed on two elements. And the initial design structural volume is 0.64, and the initial structural maximum von Mises stress is 68.62. It is assumed that an allowable von Mises stress limit is e Y D 70 for the full density material. There are the empirical parameters and , respectively, in Equations (11a) and (11c) of this paper, empirical parametersˇ0;ˇ1, andˇ2 in Section 2.3, empirical parameters " 1 ; " 2 , and " 3 in Section 4.3, empirical parameters " 4 ; " 5 , " 6 , and " 7 in Section 4.4, empirical parameters "; ; p 10 , and p 20 in Section 4.2, and empirical parameters 1 d , 2 d ;˛d el ; da , and add th in Section 4.4.1 and 4.4.3.
Although here, there are 21 empirical parameters to be determined, only and ;ˇ0;ˇ1, andˇ2 are important empirical parameters and are treated as active empirical parameters. Other empirical parameters are easily given and have been given in corresponding sections of this paper. For example, the empirical parameter values "; ; p 10 , and p 20 given in Section 4.2 by many trial calculations for shape and size designs of simple structures, can be adopted in the following examples. These five active empirical parameters are also easily given by the small quantity change requirements of structural performance parameters, the optimization speed requirements in different iteration phases and the aforementioned considerations. For example,ˇ0;ˇ1, andˇ2 only are changed by the results of two optimization phases,ˇ0;ˇ1; andˇ2 must be increased when the optimization solution needs too  The green part denotes elements whose topology variable values all are between 0.85 and 1.0, the yellow part denotes elements whose topology variable values all are between 0.7 and 0.85, the blue part denotes elements whose topology variable values all are between 0.55 and 0.7, the red part denotes elements whose topology variable values all are between 0.3 and 0.55, the dark blue part denotes elements whose topology variable values all are between 0.1 and 0.3, the dark red part denotes elements whose topology variable values all are between 0.01 and 0.1, and the black part denotes elements whose topology variable values all are below 0.01. many iteration numbers, andˇ0;ˇ1, andˇ2 must be reduced when there are some bad topological configurations in the optimization process.
Here, the empirical parameters D The topological optimization history of the L-shape beam structure, which was obtained by the proposed method, is presented in Figure 4; each figure in Figure 4 gives the corresponding topology variable distribution, and the meanings of color parts are referred to the note of Figure 4. Table I lists the structural maximum effective von Mises stresses, the volumes and structural aggregation function values of the topologies in Figures 4(1)-(15), respectively.
It is seen from Figure 4 that the solutions obtained at the outer iteration steps before the 60th iteration step do not possess clear topology configurations because of the applications of the design variable trust regions, and the topology configurations obtained at the outer iteration steps after the 60th iteration step become clearer and clearer because of the active constraint actions brought by the varied constraint limits, and possess many grey elements whose topology variables are below 0.2. Figure 4 (8) is the topology (only including elements whose topology variables are large than 0.2) obtained at the 111th iterative step, and just the phase transferring step ends. And the topology in Table I. the structural characteristic data and aggregation function values corresponding to the historical topologies in Figure 4.    Figure 5 depicts the stress cloud pictures of the structures corresponding to the topologies in Figure 4. The corresponding stress distribution of the topology obtained becomes more uniform   with an increase in outer loop iterative number until the 178th outer loop iterative step. It is seen that the stress concentration problem of this example is solved by the proposed method. The volume optimization history and the maximum effective von Mises stress optimization history of the L-shape beam structure, which were obtained by the proposed optimization algorithm, are presented in Figure 6. The volume curve in Figure 6 shows that the material volume decreases 351 with an increase in outer loop iterative number until the 178th outer loop iterative step. At the 179th outer loop iterative step, the structural volume slightly increases, and just the optimum topology was abstracted.
It can be seen from Figure 6 that the maximum effective von Mises stress of the structure slightly decreases or increases with an increase in outer loop iterative number between the 3th to the 111th Figure 9. Optimal design of the L-shape beam problem under the volume minimum formulation with a compliance constraint and a stress constraint for the full density material with an allowable von Mises stress limit e Y D 70 and the corresponding stress distribution in Ref. [24]. And its optimal design volume is 0.2086. Figure 10. Initial design domain and its load case. Figure 11. The optimization topology history of the structure with two load cases, obtained by the proposed method, and the meanings of color parts are referred to the note of Figure 4.
outer loop iterative step. At the 111th iterative step, there is a large jump of the maximum effective von Mises stress, just the phase transferring step ends, many elements with topology variables below 0.2 are removed, and there is a change of the element location of the maximum effective von Mises stress. Henceforth, the maximum effective von Mises stress gradually decreases with an increase in outer loop iterative number until the 137th outer loop iterative step. At the 137th outer loop iterative step, there is a jump of the maximum effective von Mises stress, because many elements with small topology variables were removed. The maximum effective von Mises stress slightly fluctuates between the 138th and the 162th outer loop iterative step because of the element location changes of the maximum effective von Mises stress caused by structural stress re-distributions. Henceforth, the maximum effective von Mises stress gradually increases with an increase in outer loop iterative number until the optimum topology is obtained. It is found from Figure 6 that a stable convergent results can be obtained and the structural stress level may efffectively be controlled by the proposed method. Figure 12. The stress cloud pictures of the structures corresponding to the topologies in Figure 11, for the second load case. Figure 7 depicts the optimization history of the aggregation stress gradient function f 1 T C . /, the aggregation stress constraint functions f 1;c so . / ; c D 1; 2, and the weighting aggregation stress constraint functions f 1;w . /;˛D 1; 2 for the L-shape beam structure. These aggregation function specialties are similar to the specialty of the maximum effective von Mises stress curve in Figure 6. It can be seen from Figure 7 that the aggregation stress gradient function fluctuates between the 138th and the 162th outer loop iterative step, and the aggregation stress constraint functions f 1;c so . / ; c D 1; 2, and the weighting aggregation stress constraint functions f 1;w . /;˛D 1; 2 all are far less 1 in the whole optimization process. The optimization process results of this example show that most or at least one of the aggregation constraints in the model (12) are active constraints at each iteration step. It is seen from Figure 7 that a stable convergence of the optimization process can be achieved.
X. Guo et al. (2011) [24] gave various L-shape beam designs under a compliance minimum formulation with a volume constraint and a stress constraint for the full density material with an allowable von Mises stress limit e Y D 70, and under the volume minimum formulation with a compliance constraint and a stress constraint for the full density material with e Y D 70. Comparing the L-shape beam design results in Figures 8 and 9 given by X. Guo et al. and Figure 4(15), it is found that the optimum topology obtained by the proposed method is similar to the results in Figures 8 and  9, and possesses a fewer volume than the result in Figure 9. It is very interesting that the structural local stresses can be controlled very well in the optimization process, and the stress concentration problem of this example can be solved by the proposed method without any structural total performance requirements.  .   Table II. The structural characteristic data and aggregation function values corresponding to the historical topologies in Figure 11 for the first load case.

Numbering
The maximum effective von Mises total aggregation total aggregation total aggregation stress /MPa stress gradient stress function 1 stress function 2 Volume /d m 3 The topological optimization history with topology variable distributions of the structure, which was obtained by the proposed method, is presented in Figure 11. Figure 12 gives the stress cloud pictures of the structures corresponding to the topologies in Figure 11, respectively, for the second load case. Table II lists the structural maximum effective von Mises stresses, the volumes, and structural aggregation function values of the topologies in Figure 11(1)- (13), respectively, for the first load case. For the second load case, the volume optimization history and the maximum effective von Mises stress optimization history of the structure with two load cases, which were obtained by the proposed optimization algorithm, are presented in Figure 13. Figure 14 depicts the optimization Figure 13. The volume and maximum effective von-Mises stress optimization histories,for the second load case, of the structure with two load cases, which were obtained by the proposed method. Figure 14. The aggregation function optimization histories,for the second load case, of the structure with two load cases, which were obtained by the proposed method.
histories of the aggregation stress gradient function f 1 T C . /, the aggregation stress constraint functions f 1;c so . / ; c D 1; 2, and the weighting aggregation stress constraint functions f 1;w . /;˛D 1; 2, for the second load case, of the structure with two load cases. The similar specialties in the aformentioned example can be found. And the sharp re-entrant corners have been smoothed in the optimal design shown in Figure 11 (13). The optimal topology obtained by the proposed method is very reasonable.

CONCLUSIONS
An effective STOM based on active constraint technique is developed. Examples are given to demonstrate that feasibility and effectiveness of the proposed method. The following conclusions can be drawn: (1) An equivalent topological optimization model with reduced stress constraints is proposed.
(a) An total aggregation stress gradient constraint, two aggregation stress constraints with different indexes and smooth aggregation weighting stress constraints for any load case are introduced. Being incorporated with a varied constraint limit scheme and a trust region scheme, etc., active aggregation constraints can be formed at each iterative step, a stable convergence can be achieved, and the strong control on the stress level can be carried out. (b) An equivalent topological optimization approximation model proposed in this paper, benefits to the alleviation of stress concentrations and the control of all local stresses.
(2) A set of stress quadratic explicit approximations is given, based on the MMA.
(3) A set of smooth dual algorithm for the optimization sub-problem with artificial variables and many possible non-active design variables, is proposed. A novel structural topological optimization method with stress constraints and its algorithm are formed, (4) Examples are provided to demonstrate that the structural maximum effective von Mises stress is controlled very well during the optimization procedure, and the stress concentration problem is solved by the proposed method. The proposed method is very effective for solving the stress-based topology optimization problem. where the meanings of all notations are referred to Section 2 and 3.1 of this paper. Therefore, Equations (15)-Equation (17)  ; Q is similarly derived as follows.