A framework for plasticity-based topology optimization of continuum structures

In this paper, a framework is proposed for topology optimization of continuum structures considering plasticity. The method merges the rigid-plastic analysis and the density-based topology optimization. To obtain a clean black-and-white design, the density in the objective function is penalized using an exponential function. The solution of the final plasticity-based topology optimization problem exhibits as a sequence of second-order cone programming (SOCP) problems that can be resolved efficiently using the advanced primal-dual interior point method. Compared to the conventional stress-constrained topology optimization techniques, the developed method accounts for plasticity and the finite element analysis of structures does not need to be carried out separately. Further-more, the proposed method requires no relaxation techniques for imposing local stress-constrain and possesses good computational efficiency for large-scale problems.


INTRODUCTION
Since the landmark work on the homogenization method for topology optimization of structures, 1 the research field of topology design has been attracting extensive attention from both academia and industry. To date, numerous approaches have been developed to solve the problem of topology optimization such as the density-based method, 2,3 the level set method, 4,5 the evolutionary structural optimisation method 6 and its variants, 7,8 the phase field method, 9 and the moving morphable components 10,11 to name a few. A comprehensive review of these methods was provided in reference 12 for comparison. The density-based method is maybe the most widely used one among these topology optimization approaches. Despite vast efforts devoted to the algorithm development for the density-based topology design, the focuses of earlier works are majorly placed on the conventional stiffness-based design of elastic continua. For instance, the problem commonly presents as a minimization of elastic strain energy of a structure subjected to given external loads and a volume constraint of materials. 13 The material is considered as elastic, and the strain energy is a function of the displacement and a newly introduced continuous 'density', ranging from 0 to 1, to indicate whether a point of space is occupied by materials. The formulation of topology optimization for elastic materials is well established and can be resolved efficiently using mathematical programming. However, the optimal design led from this approach does not guarantee the feasibility of the stress states with respect to material strength. In other words, it is likely that the stress states of the designed structure are above the yield limit of the material when subjected to the considered external load. Remarkably, the strength of structures is among the most important concerns in practical applications. Hence, in real-world applications, structures from the conventional stiffness-based topology optimization are subjected to sequential modifications and improvement in a later stage to ensure the strength criteria of materials. 14 Further inclusion of stress criteria in the optimization routine is regarded stress constrained topology optimization. 12,15 One way to achieve this goal is by including the yield criterion as additional constraints in the conventional stiffness-based topology optimization. Despite the forthright extension in terms of the mathematical formulation, the resulting problem is much more difficult to address. Alternatively, the stress constrained topology optimization can be formulated as the minimization of the volume subjected to equilibrium equations, which also satisfy the stress constraint. Although this strategy is widely used for considering stress constraints, several challenges burden its application. 12 A pronounced issue is the so-called stress singularity. 16 When the density approaches zero, the low stiffness may cause high deformation and sequentially large local stresses. The local constraints thereby saturate which prevents the removal of materials and results in a solution of the substantial gray region whereas a crisp solid/void result is desired. 17 To alleviate the singularity problem, relaxation techniques such as the -relaxation approach 16 and the qp-approach 18 have to be applied. Another issue is the high computational demand when the stress constraints are enforced at the local points of each element, for example, the numerical integration points of each element. To reduce the computational cost of design with local stress constraints, a single global stress constraint is enforced in the topology optimization using aggregating functions such as the p-norm or the Kreisselmeier-Steinhauser (KS) function 19,20 which, however, leads to a weaker control of the local stress. A compromise approach is to group the elements into blocks based on which regional constraints on stresses are enforced. 21 This strategy reduces the number of constraints dramatically compared to the local-constraints approach while retaining control of the stress behavior.
More recently, a topology optimization approach was proposed for plane strain problems accounting for material strength. 22 Neglecting the elastic behavior, this approach targets the final plastic limit state of structures by combining the direct limit analysis and the density-based topology optimization in the same framework. Specifically, the optimal design problem is exhibited as the minimization of material volume subjected to a stress field which is both statically admissible (i.e., fulfilling the equilibrium equations, stress vector continuity and stress boundary conditions) and plastically admissible (i.e., satisfying f (σ) ≤ 0 where f (σ) is the plastic yield criterion of the material). The material density and stress are the design variables and the material strength is proportional to the material density. 23 The stress constraints are enforced at the element level. Later, Herfelt et al. 24 proposed a more general formulation for this approach that both the lower bound and upper bound limit analyses can be carried out conveniently in the optimal design by using adequate elements. Consequently, the optimal design can be bounded by using the upper bound and relaxed lower bound elements. In the form of a standard convex optimization problem, the final optimization problem is resolved forthrightly using the primal-dual interior point method, which also ensures the solution is the global optimum. Despite its novelty, the topology optimization method developed in reference 24 only leads to a gray-scale design, which dramatically reduces its attraction.
In this paper, a plasticity-based topology optimization framework is developed based on reference 24. Similar to reference 24, the developed method combines limit analysis and density-based topology optimization. However, the plasticity-based topology optimization method developed in this study leads to a black-and-white design, rather than a gray scale design as in reference 24, which is more realizable in manufacturing. The method is developed by first merging the density-based technique into the rigid plastic analysis, resulting in a plasticity-based topology optimization formulation. Further embracing the Solid Isotropic Microstructure with the Penalization for intermediate densities (SIMP) method to steer the intermediate density, the formulation leads to black-and-white layouts instead of gray designs via iterations. Because of the discontinuity of the stress field between elements in the proposed formulation, the filtering operation is introduced for improving the optimal design. Compared to the conventional stress-constrained topology optimization with SIMP, the topology optimization problem presented in this study is resolved straightforward using the advanced primal-dual interior point method with high computational efficiency. It is also found that black-and-white designs can be obtained using this approach without employing any relaxation techniques. Further, the layout from the developed approach is more economical given the concern of plasticity in the developed approach. The correctness and robustness of the proposed method are demonstrated by simulating some typical topology optimization problems.

Rigid-perfectly-plastic theory
For a rigid-perfectly-plastic body with volume Ω and surface Γ = Γ u ∪ Γ t where Γ u and Γ t are the kinematic and traction boundaries, respectively, with Γ u ∪ Γ t = ∅, the governing equations consist of • Equilibrium equation • The strain-displacement relation • The constitutive relation • The boundary conditions where. is the Cauchy stress; b is the body force; is the strain; u is the displacement; is the plastic multiplier; f ( ) is the yield function; u is the prescribed displacements (i.e., u = 0); t is the prescribed traction; is the collapse load factor; N consists of components of the outward normal to the boundary Γ t ; and T is the transposed gradient operator, in a plane-stress case, taking the form of According to references 25,26, the rigid-perfectly-plastic analysis can be formulated as the following min-max optimization problem In the above, the maximization part renders the principle of maximum plastic dissipation. The minimization part, on the other hand, concerns the total potential energy and corresponds to equilibrium enforcement. The upper and lower bound theorems follow as special cases of the optimization problem (7). The equivalence between the optimization problem (7) and the governing equations listed in (1)- (5) has been demonstrated in references 25,26 where the Karush-Kuhn-Tucker (KKT) conditions associated with (7) are derived. Using mixed finite elements, the interpolation for the stress and displacement fields are Substituting the above interpolations into the min-max problem (7) results in where N G is the total number of interpolation points for the stress field, and The minimization part of principle (9) with respect to the displacementû can be resolved analytically leading to a maximization problem According to reference 26, problem (13) results in a rigorous upper bound when the mixed finite element shown in Figure 1 is employed. It should be stressed that a solution to the limit analysis problem is a pair of fields ( , u). The solution can be sought via either the statical (or lower bound) method, which involves only stresses as variables, or the kinematical (or upper bound) method, which involves only displacement. However, the mixed approach, involving both stresses and displacements, in some particular cases can reproduce rigorous upper bound solutions as indicated in references 26-28. The mixed element used in this study is from reference 26 that the corner nodes are used for approximating the stress field, which are also numerical integration points, whereas both the corner nodes and the nodes at the middle point of edges are for the displacement field. As a result, the stress field varies linearly within the element and is discontinuous between elements. The displacement field on the other hand is quadratic within the element and continuous between F I G U R E 1 An illustration of the mixed finite element. 26 elements. Such a mixed finite element is indicated as an upper bound element in reference 26, which is also the one used for upper bound limit analysis in the commercial software OptumG2. 29

Plasticity-based topology optimization
The plasticity-based topology optimization can be constructed in the framework of rigid-perfectly-plastic analysis. By introducing a new design variable-'density' ∈ [0, 1], the optimal problem presents as a minimization of the material volume subjected to force balance equations and yield criteria as below where both f b and the yield function depend on the density . In the topology optimization problem (14), is a known factor since the given external load is denoted by t 0 .
The von Mises yield criterion which is prevalent in the stress constrained topology optimization 15 is adopted in this study. The corresponding yield function is where f y is the yield stress and J 2 is the second invariant of the deviatoric stress. In plane stress cases, it is expressed as When a point has = 0 denoting a void point, all stress components are null due to (15), which mean this point cannot sustain any stresses. Because the yield criteria are enforced on all the stress interpolation points, the density field is approximated in the same way as the stress field, which is Substituting Equation (17) into (14), the optimal design problem is where The derived topology optimization problem (18) is similar to the one presented in reference 24. By introducing a new set of variableŝ= the von Mises yield criterion is reformulated as Thus, the topology optimization problem (18) turns to which is a standard second-order cone programming (SOCP) problem and can be resolved efficiently using the interior-point method available in advanced optimization engines. It is remarkable that, theoretically, the density can be treated as an integer (i.e., either 0 or 1) in mathematical programming. The resulting mixed-integer SOCP, however, is difficult and extremely slow to solve. Thus, the density is treated as continuous for the sake of computational efficiency in reference 24, which leads to a 'gray' optimal design. In this paper, the derived topology optimization problem will be further modified via the penalization of intermediate density to result in a black-and-white solution, which is more desirable.

PENALIZATION AND FILTERING
To penalize the intermediate density, the objective function ∫ Ω dΩ in (14) is replaced by ∫ Ω h( )dΩ in which h( ) = e p(1− ) with p being a factor greater than or equal to 1. Because of the exponential term in the objective function, the discretized problem is a non-SOCP problem and cannot be resolved as for (18). To seek its solution, the function is approximated as h( ) = c n where c n is a known factor calculated based on the density obtained from the last iteration step, for instance where subscript n indicates the corresponding variable at nth iteration step. Then the discretized topology optimization with penalization reads min (̂,ρ)Lŝ Problem (25) is now a standard SOCP problem that can be resolved forthrightly using the primal-dual interior point method, which will be briefly summarized in the next section.
To sum up, the solution algorithm for the proposed plasticity-based topology optimization with penalization is as follows: (i) Assume the density = 1 for the entire computational domain; (ii) Calculate c n at stress interpolation points based on the known density n using (24); (iii) Solve the optimization problem (25) using the primal-dual interior point method to attain the density field and the stress field; (iv) Perform density filtering across the domain; (v) Stop the iteration if convergence criterion is satisfied (i.e., the objective function, Obj =L̂in Equation (25) The density filtering operation in (iv) is carried out by first calculating the density of each element, e , as an average of the density at the three corner nodes of the element. According to reference 30, the filtered density of an element is theñ where v i and e i are the volume and the density of the ith element, respectively, and N e denotes the total number of elements located in the filtering region of the element, which is a circle of radius R. The weighting function is w (x i ), and the Gaussian (bell shape) distribution function is used such that (29) in which x e is the coordinate of the centroid of the element whose density is filtered and x i is the coordinate of the centroid of the ith element within the filtering region. The parameter d is R∕2 with R being 1.5 times the mesh size.
Remarkably, because the proposed topology optimization problem (25) is derived from the rigid-perfectly-plastic analysis problem (13), the layout resulting from (25) fulfills the governing Equations (1)- (5). In other words, the applied force in fact is the maximum force the designed structure can sustain based on the plastic theory. This contrasts with the conventional stress-constrained topology optimization, which provides a more conservative solution. This point will be further discussed in the numerical example section.

PRIMAL-DUAL INTERIOR POINT METHOD FOR SOCP
In general, a second-order cone programming (SOCP) problem can be written in the form where c, x ∈ R n , b ∈ R m , A ∈ R m×n , and K i is one of the following cones: Quadratic cone Rotated quadratic cone Apparently, both the objective function and the equality constraint in the final optimization problem (25) are linear as these in (30). The inequality constraints in (25) are also the standard quadratic cones, for instance, the cone in (31). Thus, the final optimization problem (25) is a standard SOCP problem of the form (30). The advanced primal-dual interior point method for solving the standard SOCP problem (30) has been detailed in reference 31. To seek the solution, the dual problem of (30) is first defined max (y,s) b T y subject to A T y + s = c (33) where K * i is the dual cone of K i such that Noting that the cone K i in (30) is self-dual meaning K i = K * i , the dual problem (33) is then expressed as max (y,s) b T y subject to A T y + s = c (35) According to the duality theory, 32 solving primal problem (30) or dual problem (35) is equivalent to solving the system The above system (i.e., (36)) involves neither the gradients nor Hessians of the nonlinear cone constraints for defining the optimal point of the SOCP problem. Indeed, system (36) can be regarded as a generalization of linear programming as indicated in reference 33. It can be resolved efficiently by employing a generalization of the Goldman-Tucker homogeneous model for linear optimization. We refer readers to reference 31 where the equivalence between the SOCP problem (30) and the system (36) is proven and an efficient algorithm, namely the primal-dual interior point method, for solving (36) is introduced in detail. Note that the algorithm documented in reference 31 also leads to an advanced optimization engine MOSEK, 34 which is adopted in this study for solving the final topology optimization problem (25) which is a SOCP problem.

NUMERICAL EXAMPLES
In

A plate under shear load
The first example concerns a plate subjected to a shear load as shown in Figure 2. This is a well-known stress-constrained topology problem. Despite its simplicity, this problem clearly illustrates the difficulties of topology design with stress constraints and, thus, serves as a classical benchmark for stress-constrained topology optimization schemes. 18,20,[35][36][37] In this study, the setup of the problem is in line with that in reference 18. The size of the plate is 4 × 8 m 2 and the left side is clamped. The applied shear force, F = 1 kN, is distributed along a central portion of length 1.5 m on the right surface. The yield stress of the material is f y = 1 kPa. The black region around the load application zone is enforced to be full material. Figure 3 illustrates the layout and the normalized von Mises stress (i.e., the ratio of von Mises stress and f y ) obtained from the developed approach and the PolyStress code developed in reference 38 a code for local stress-constrained topology optimization using the augmented Lagrangian method. The mesh size (e.g., the edge length of a typical element) in both simulations is h e = 0.05 m. Clearly, similar material layouts are produced by the plasticity-based topology optimization method developed in this study and the conventional local stress-constrained topology optimization method, for instance, PolyStress in reference 38.
However, it should be stressed that a perfect agreement between layouts from these two methods is not essential given the significant differences in the fundamentals underpinning the two approaches. To clarify this point, we consider an elastic perfectly plastic material whose behavior is shown in Figure 4A. The conventional stress-constrained topology optimization solves the elastic equilibrium equation, and the yielding of all material points is suppressed. In other words, it searches for solutions within the range that a continuum structure behaves in elasticity. Thus, in this approach, the F I G U R E 2 A plate under shear load. structure is deemed to fail when any local yielding starts (see Figure 4B) which is too conservative. On the other hand, the developed plasticity-based topology optimization approach accounts for the plastic deformation, implying that local yielding at the material point level is allowed if the global structure is still stable. This targets the real limit load of the whole structure (see Figure 4B). Hence, theoretically, the volume ratio (defined as the volume of the designed layout over the volume of the original domain) from the conventional stress-constrained method will be larger than that from the developed plasticity-based approach. This can be demonstrated in Figure 5 where the convergence history of the two approaches for this problem is illustrated. A converged solution is attained with 6 iterations for the developed approach and around 15 iterations are required in PolyStress. The conventional stress-based approach produces a more conservative solution (volume ratio of 0.276) than the developed plasticity-based approach does (volume ratio of 0.256).
Note worthily, holes are observed in the solutions from both approaches as shown in Figure 3A,B. This is because a very fine mesh is used in the simulation and the stress state in these areas is quite low. Consequently, in the iterations, the corresponding materials are suppressed in these areas leading to small holes. As indicated in reference 13, this is a commonly observed phenomenon in density-based topology optimization and can be alleviated by adopting a relatively larger mesh if no detailed design is required. As shown in Figure 6A, with a larger mesh size (i.e., h e = 0.2) the final layout from the proposed method does not have any holes. Figure 6B,C show the layouts from reference 18 where the conventional stress-constrained topology optimization method is adopted using the same mesh size. As seen, relaxation techniques must be employed otherwise a blurry solution is obtained ( Figure 6B). Note that the PolyStress also employs the relaxation technique to attain those black-and-white layouts shown in Figure 3A. The developed plasticity-based topology optimization method, on the contrary, can result in a black-and-white design without the use of any relaxation technique, which is majorly attributed to the strong convergence properties of the primal-dual interior point method for second-order cone programming problems. 31,34

Clamped beam
The second example is a double-clamped beam as plotted in Figure 7. This problem has been considered using different topology optimization techniques. [39][40][41] In this study, the setup is in line with that reported in reference 39. The length of the beam is 20 m, and the height is 5 m. The yield stress of the material is set to be 300 kPa. A uniformly distributed vertical force of 383.2 kN is applied along 2.5 m at the center of the top surface. Due to symmetry, only the right half of the domain is optimized.
In the topology design, the mesh size is he = 0.1 m. The dark region shown in Figure 7 has a fixed density of 1. Figure 8 illustrates the evolutionary history of the structure in the optimal design using the proposed method. As seen,   if no iteration is carried out meaning the approach proposed in reference 24 is used, the solution is a gray layout. With 5 iterations, a satisfactory black-and-white layout is attained. Further iterations (i.e., iteration number greater than 5) have little influence on the black-and-white layout, although the convergence criterion (i.e., ‖Obj n+1 −Obj n ‖ Obj n+1 ≤ 1 × 10 −4 ) is fulfilled at the 10th iteration for the proposed method ( Figure 9). Figure 9 also shows that the converged volume ratio from the conventional stress-constrained method (i.e., PolyStress) is higher than that from the plasticity-based method developed in this study, which echoes the statement in section 3. Despite that, the final layouts from PolyStress and the developed method are similar as shown in Figure 10. Additionally, it can be seen from Figure 10 that, for areas of low von Mises stress, the corresponding material density is low. This is the expected case for both the PolyStress simulation and the simulation using the developed plasticity-based method.
To investigate the influence of the parameter p of the exponential penalization function on the design, the problem is re-analyzed using p = 0 − 6, and 10. The converged layouts from different simulations are shown in Figure 11. A very gray design is attained for p = 0 since in this case no penalization is enforced. An increase of p leads to a clearer solution. For p = 5, a satisfactory clean black-and-white layout is obtained. Further increase of p has little influence on the black-and-white design. For instance, the layout from p = 5 coincides with these from p = 6 and 10. Thus, p ≥ 5 is recommended for the proposed method.

Bridge design
The last example considers the problem of bridge design. Figure Figure 13 shows the optimization history of the bridge obtained from the proposed design method. It is shown that a converged black-and-white design can be achieved with 10 iterations. In this example, we set the iteration number to be 50 even though the convergence criterion is fulfilled at the tenth iteration. As shown in Figure 13, there is little difference in the layouts after 10 iterations. This verifies the good convergence property of the proposed iteration approach for seeking a black-and-white design.
It is recognized that the computational cost of stress-constrained topology optimization depends heavily on the total number of stress constraints. Thus, efforts are devoted to reducing the total stress constraint number by developing a global stress constraint approach and grouped aggregation approach. 21,36,42 To disclose the relationship between the computational demand of the proposed method and the number of stress constraints, the problem is re-analyzed using elements of sizes h e ranging from 0.  Figure 14. Overall, the cost of the proposed method is linearly proportional to the total number of stress constraints. This is an admirable feature given that computational cost normally increases polynomially with the number of degrees of freedom in the conventional stress-constrained topology optimization method.

CONCLUSIONS
In this paper, a method for black-and-white topology optimization of continuum structures is proposed accounting for plasticity theory. The method is exhibited as a sequence of continuous convex topology optimization problems, in the standard SOCP form, that can be resolved efficiently using the advanced primal-dual interior point method available in modern optimization engines. The penalization of the density is performed in the objective function to steer the intermediate density towards integers 0 and 1 and results in a black-and-white optimal design with the help of the density filtering operation. Compared to the conventional stress constrained topology optimization, the proposed method requires no separate finite element analysis of the continua since the density and the stress field can be solved simultaneously in the used primal-dual interior point method. Additionally, the relaxation techniques commonly required in the conventional density-based topology optimization with stress constraints are not necessary for the proposed method. Because the proposed method is developed based on the rigid-perfectly-plastic analysis, the design targets the limit load of the continuum structure rather than the load under which the structure behaves in pure elasticity. Despite the high nonlinearity of the strength-based topology optimization, the proposed approach shows a good computational efficiency that the computational demand is linearly proportional to the number of used element nodes.