Plastic work constrained elastoplastic topology optimization

An elastoplastic topology optimization framework for limiting plastic work generation while maximizing stiffness is presented. The kinematics and constitutive model are based on finite strain linear isotropic hardening plasticity, and the balance laws are solved using a total Lagrangian finite element formulation. Aggregation of the specific plastic work combined with an adaptive normalization scheme efficiently constrains the maximum specific plastic work. The optimization problem is regularized using an augmented partial differential equation filter, and is solved by the method of moving asymptotes where path‐dependent sensitivities are derived using the adjoint method. The numerical examples show a clear dependence on the optimized maximum stiffness structures for different levels of constrained specific plastic work. It is also shown that due to the history dependency of the plasticity, the load path significantly influences the structural performance and optimized topology.


INTRODUCTION
Topology optimization is a powerful computational method which has gained much attention in the past three decades. The method finds an optimal material layout within a prescribed design space to minimize a cost function and satisfy constraints; it can be applied to a variety of physical problems. 1,2 Among the developed topology optimization schemes, most are restricted to the design of linear elastic structures, where the common objective is to find an optimal compromise between stiffness and material usage. Less attention has been given to the topology optimization of structures involving inelastic material response.
Topology optimization with inelastic material modeling has been applied for design objectives such as maximum ductility 3 and maximum stiffness for reinforced concrete structures. 4 However, the objective of most research on topology optimization with inelasticity is to maximize energy absorption. This is an important design consideration when designing protective systems, for example, in automotive crashworthiness for impact and dampers or mechanical systems for seismic resistance. 5 Recent studies on optimizing energy absorbing structures has shown great promise. Effective energy management structures were designed by Nakshatrala and Tortorelli using a transient small strain elastoplastic formulation. 6 An alternative adjoint sensitivity scheme was presented in the work of Kato et al. 7 Finite strain effects were incorporated for quasi-static structures by Wallin et al., 8 and for dynamically loaded structures by Ivarsson et al., 9 and Alberdi and Khandelwahl. 10 Energy absorption was maximized for small strain elastoplastic microstructures in the work of Alberdi and Khandelwahl, 11 and later for finite strain viscoplastic microstructures in the work of Ivarsson et al. 12 Li et al. 13 studied the effects of kinematic versus isotropic hardening materials for the design of maximum energy absorbing structures under cyclic loads; the same authors also considered anisotropic elastoplastic materials. 14 The design for maximum energy absorption has also been combined with damage constraints by Li et al. 15 using a coupled elastoplastic damage model and an aggregated damage constraint. The aggregation technique was later used for designing structures with pressure sensitive materials subject to damage constraints by Alberdi and Khandelwahl, 16 and elastoplastic materials subject to fracture by Li and Khandelwahl. 17 Both energy absorption and stiffness are seldom optimized together as they are closely related, that is, high initial stiffness commonly gives relatively high energy absorption. However, both of these measures need to be considered when optimizing structures for maximal stiffness subject to plastic dissipation constraints. Amir 18 proposed an elastoplastic topology optimization algorithm for constraining the total plastic strain and volume while maximizing stiffness. The optimization scheme was proposed as an alternative to stress-based optimization. Although more computationally demanding than linear elastic frameworks where the von Mises stress is constrained, this algorithm constrains the plastic strain norm to nearly zero, thus guaranteeing that the maximum stress is close to the yield stress. The study was limited to proportional loading, that is, the influence of the load path was not investigated.
In low-cycle fatigue applications, structures should be designed to be stiff without incurring much plasticity to maximize the fatigue life as predicted by, for example, the Coffin-Manson law (see, e.g., Lemaitre and Chaboche 19 ). However, to the best of our knowledge, no topology optimization design study that treats such considerations exists in the literature. In addition, the majority of topology optimization studies that do consider elastoplastic material response employ proportional loading and thus do not fully investigate the history dependence of the response. Motivated by this gap in the literature, we propose a topology optimization framework for designing structures for maximum stiffness subject to maximum specific plastic work constraints during a single load cycle, thereby ensuring adequate fatigue life. We also investigate the influence of the load path through nonproportional loading.
Our topology optimization utilizes the solid isotropic material with penalization (SIMP) material model 20 for the material interpolation. Efficient schemes for interpolating stiffness (i.e., Young's modulus), mass and yield stress have been presented previously for plastic work maximization. Unfortunately, interpolation schemes for constraining the specific plastic work are not fully understood. In the works of Amir 18 and Alberdi and Khandelwahl, 16 the yield stress and the Young's modulus were interpolated using separate SIMP exponents, similar to the qp-relaxation for stress-based topology optimization 21 and the separate penalization approach proposed by Maute et al. 3 for elastoplastic topology optimization. Alberdi and Khandelwahl 16 also used an additional interpolation scheme to model damage, such that high damage in regions devoid of material is effectively ignored. Notably, in the works of Amir 18 and Alberdi and Khandelwahl, 16 the yield strain is dependent on both the load magnitude and the volume fraction; thus, the response can oscillate between elastic and inelastic during the optimization. Enforcing a delay on the onset of plasticity in intermediate solid-void regions via separate penalization exponents alleviates convergence problems in the solution of the nonlinear state equations, but it also lessens the sensitivity information in these regions which hinders the convergence of the optimization. To increase this sensitivity information and consequently improve the optimization process, the interpolation for the yield stress should mimic the interpolation of the stiffness. This technique was employed by Li and Khandelwal,17 where the yield stress was modeled such that plastic evolution can occur in regions nearly devoid of material.
To capture the plastic work, we use a kinematic and material model based on finite strain elastoplasticity, which has previously been used by Wallin et al. 8 In particular, we apply a material interpolation scheme where both the elastic and elastoplastic response follows an identical power-law. The peak value of the specific plastic work is approximated by spatial aggregation using the p-norm, and subsequently constrained in the optimization problem. The sensitivities required for the method of moving asymptotes (MMA) nonlinear programming algorithm are derived by the adjoint scheme presented in the work of Michaleris et al. 22 The optimization framework is used to solve several example problems in which we consider moderate load magnitudes, as is the normal case in low-cycle fatigue applications, resulting in strains that are not necessarily large. That said, the finite strain model within our formulation comes with little extra computational cost, as compared with a small strain formulation, which increases its applicability. This article starts with the presentation of finite strain elastoplasticity theory in Section 2, which is followed by the numerical solution procedure in Section 3. The plastic work constraint and optimization problem are presented in Section 4 which is followed by the presentations of restriction and thresholding methods in Section 5. Section 6 details the material interpolation schemes. The path-dependent sensitivity scheme is presented in Section 7, and supplemented by the Appendix. In Section 8, we present examples, and finally, conclusions are drawn.

BALANCE LAWS AND CONSTITUTIVE MODEL
We consider finite deformation and therefore the deformed configuration Ω is distinguished from the reference configuration Ω o . The motion of a material particle, which at time t = 0 is located at position X ∈ Ω o , is described by the smooth mapping such that its position at time t > 0 is x = (X, t) = X + u(X, t), where u is the displacement field. The deformation gradient is defined as F = o where o is the gradient with respect to X. We assume a multiplicative split of the deformation gradient, i.e., F = F e F p , where F e and F p are the elastic and plastic components. We also find it useful to define J e = det(F e ) and the elastic Finger strain tensor b e = F e F eT and its isochoric part b e = (J e ) −2∕3 b e . Neglecting body force and inertia, the mechanical balance law and associated boundary conditions are given as where S is the second Piola-Kirchhoff stress, related to the Cauchy stress via = J −1 FSF T , and the boundary of the body Ω o is comprised of the complimentary subsurfaces Ω ou and Ω ot , over which displacement u and traction t are applied. The weak form of (1)-(3) is the basis for the finite element method, and it is stated as: find the kinematically admissible u such that for all kinematically admissible virtual displacements u. In the weak form (4), we introduced the Lagrangian virtual . We use the linear hardening finite strain plasticity constitutive model presented in Simo and Miehe. 23 A viscous extension of this model can be found in our earlier work. 9,12 The hyperelastic component of the constitutive model is defined by the isotropic strain energy potential where and are the bulk and shear modulii. To distinguish elastic from elastoplastic response, the von Mises-type yield function is employed. Elastic response occurs when f < 0 and plastic response occurs when f = 0. In (6), dev is the deviatoric part of the Kirchhoff stress, = J , which is obtained from (5) via = 2 ( w e ∕ b e ) b e , and is the current yield stress, where yo is the initial yield stress, K = H is the drag stress and H is the linear hardening modulus. The evolutions of the internal variable and deformation b e are given by maximum dissipation theory, cf. for example, the work of Ottosen and Ristinmaa, 24 that is, where is the plastic multiplier that is obtained via the KKT-conditions and the consistency conditioṅf = 0. In the evolution Equation (8), the Lie derivative of b e is defined as and n = dev ∕ ‖ ‖ dev‖ ‖ is the normal to the f = 0 yield surface. The internal variables are stored in the vector w = For further details of the evolution laws, we refer to the works of Simo and Miehe 23 and Simo and Hughes. 25

NUMERICAL SOLUTION PROCEDURE
The weak form (4) of the equilibrium equations (1) is discretized using the Galerkin finite element approach, giving the discretized residual equation where u and w are the nodal displacement and integration point internal variable vectors, and the internal force vector at time step n is where B (n) is the strain-displacement matrix, 26 and we use Voigt notation to represent the stress and strain. Subsequently, we drop the superscript (n) and assume that all quantities are evaluated at time step n. The integration of our constitutive evolution equations, that is, the evolutions of b e and , follows the procedures laid out by Simo and Hughes 25 and Wallin et al., 8 and is therefore omitted here. Ultimately, N gp local Gauss point-i residual equations are enforced for all n = 1, 2, … , N.
We use displacement controlled loading wherein the nodal displacement vector is decomposed as u = where u p contains the displacements prescribed on the boundary Ω ou , and the free displacements u f are to be determined. In the loading program, u p is updated and then (u f , w) is evaluated by Newton's method, in which the residuals (12) and (14) are linearized in a Schur complement approach to obtain where the consistent tangent operator K is partitioned 27-29 as This tangent matrix is also utilized in the sensitivity analysis. For details of the linearization, see, for example, the works of Simo and Miehe, 23 Wallin et al., 8 Ivarsson et al. 9

OPTIMIZATION PROBLEM DEFINITION
The partitioning of u and a similar partitioning of F int enable the computation of the compliance C (n) = (F (n) int,p ) T u (n) p at load step n, and the end-compliance In the optimization, we maximize C (N) to generate structures with high stiffness. 28 This contrasts the usual compliance optimization problem that minimizes the compliance. The difference being that here we prescribe nonhomogeneous displacement boundary conditions rather than the usual traction boundary conditions.
To limit the peak plastic work, the specific plastic powerẇ p = y ( )̇is integrated over the load duration to give the specific plastic work which is discretized in time as where the "hatted" quantities( ⋅) are interpolated using the physical volume fraction (discussed in Section 6). Next, we approximate the maximum value ofŵ p over the spatial domain via the p-norm to formulate the constraint where g 1 > 0 is the maximum allowable specific plastic work. In (20), we use P = 8, and refer to Holmberg et al. 30 for a discussion on the effect of the p-norm exponent. The volume V, cf. (30), is constrained by the maximum value V via Finally, we are in position to formulate the optimization problem. Using the constraints (21) and (22), and the end-compliance (17), the optimization problem is stated as solving for the discretized design parameters . Note that the residuals (12) and (14) should hold for all n = 1, 2, … , N and are therefore included as implicit constraints in this reduced space optimization problem. In solving the optimization problem above, we use the normalization technique presented in the work of Le et al. 31 for scalingŵ p PN to improve the approximation of the maximum specific plastic work, that is, we replace the g 1 -constraint (21) by where c is calculated in each optimization iteration I as We follow Le et al. 31 and set the parameter = 0.5 if c oscillates between optimization iterations, otherwise we set = 1. The first optimization iteration uses 1 = 1.

RESTRICTION AND THRESHOLDING
We use restriction to formulate a well-posed topology optimization problem. The design field that is piecewise uniform over the elements is used to define the continuous filtered design field̃through the Helmholtz' partial differential equation (PDE) filter. This filter was originally proposed by Lazarov and Sigmund 32 and recently modified by Wallin et al. 33 to control the boundary effects. Ultimately, we use the design field to evaluate the continuous filtered design field̃by solving subject to where l o and l s are the filter bulk and surface length scales, respectively; N o is the outward unit normal vector to the boundary in the reference configuration, and Δ o is the Laplacian with respect to the material coordinates X. In accordance with Wallin et al., 33 we adopt l s = l o which mimics filter padding wherein the domain edges are padded with void elements, resulting iñ≈ 1 1+(l s ∕l o ) = 0.5 at the design domain boundaries for a solid design variable field. The PDE (26) is solved with a finite element formulation; cf. the work of Wallin et al. 33 for details.
Unfortunately, the filtered desigñfield will contain intermediate regions rather than regions that are either completely full or devoid of material, i.e., regions where 0 <̃< 1 rather thañ= 0, 1. The stress and consequently the specific plastic work evaluated in these intermediate regions is nonphysical, and can severely affect the optimization algorithm. To reduce the size of the intermediate material regions, we utilize Heaviside thresholding, 34,35 and threshold the filtered field̃using a smoothed Heaviside step function H, which approximates the unit step function u s according to where the parameters H and control the unit step approximation such that lim H →∞ H(̃) = u s (̃− ). The filtered-thresholded density ∈ [0, 1] is the physical volume fraction field that is used to evaluate the material constitutive parameters.
In the optimization, we apply a continuation method for increasing the thresholding "steepness" parameter H . Previous studies have presented modified thresholding techniques in which the thresholding parameter is modified in each optimization iteration to conserve volume. 36,37 Such techniques purportedly alleviate oscillatory behavior of the optimization convergence that otherwise often appear during thresholding parameter continuation, cf. the work of Xu et al. 36 We have implemented and investigated such techniques, however in our experience we found that a constant = 0.5 is the most robust with the present optimization formulation.

INTERPOLATION OF MATERIAL PROPERTIES
After the design field has been smoothed and thresholded, we can finally interpolate the material properties using the physical volume fraction field = (̃( )). We interpolate the stiffness, mass and plastic work parameters using the modified version 38 of the SIMP scheme 20 through where p is the interpolation exponent, o is a small residual value and o is the "unscaled," i.e., nominal material parameter of the interpolated parameter , for example, the elastic parameters. Note that (0) = o o and (1) = o . For p ≠ 1, the interpolated material parameter is either artificially high or low in relation to a linear interpolation.

Material parameter interpolations
The interpolation for the bulk and shear modulii parameters o = , use p = p E > 1, and the initial yield stress and hardening modulii o = yo , H use p yo . A common strategy for topology optimization involving inelastic material response is to choose separate penalization exponents, p E and p yo for the stiffness and the initial yield stress, respectively, with p yo < p E . 3 This is done to avoid inelastic strain in intermediate material regions and to completely avoid plasticity as → 0. The particular choice of SIMP exponent p yo = p E − 0.5 has been successfully employed in many studies, see, for example, the works of Amir 18 and Alberdi and Khandelwal. 16 The approach of separate penalization exponents has also been proposed for stress-constrained topology optimization of linear elastic structures as the "qp-approach" 21 to relax the stress in intermediate regions. For topology optimization of maximum energy absorbing structures, this approach is efficient as the yield strain for intermediate regions is artificially higher thereby suppressing plasticity, which improves the convergence when solving the nonlinear state equations. However, different penalty parameters can lead to widely different optimized designs. 39 Furthermore, the artificially high yield strain renders near zero sensitivity of the plastic work constraint which hinders the convergence of the optimization.
To allow plasticity in intermediate material regions, and consequently nonzero plastic work sensitivity in these regions, we assign This interpolation scheme has been used for interpolating the elastic parameters , , initial yield stress yo and hardening parameter H, with equal material constituent residual values o , in the work of Li and Khandelwal, 17 where fracture resistant energy absorbing structures were designed. We follow this approach and choose o = 10 −6 . With these choices, the uniaxial yield strain is independent of volume fraction , and plasticity in material devoid regions can occur; in fact, when the modified SIMP interpolation scheme (29) is used for both the elastic parameters and initial yield stress with equal o , the uniaxial yield strain for = 0 and = 1 is equal regardless of the choice of penalty exponents.
We have also investigated different relaxation techniques for avoiding plasticity in material devoid regions, specifically the -relaxation technique 40 previously used in the works of Luo and Kang 41,42 for relaxing the initial yield stress yo . No significant differences for the optimized designs were observed, however, computational cost for the simulation can be saved via the -relaxation by suppressing plasticity and thus simplifying the local constitutive equations in materially devoid regions. To avoid introducing additional parameters of the optimization, we do not use such relaxation.
In the numerical examples, we use the exponent p E = p yo > 1, and linear interpolation for the mass, that is, the volume constraint (22) is evaluated with This penalization choice results in a low stiffness to volume ratio for intermediate regions, which encourages the optimizer to reduce the size of them.

SENSITIVITY
The optimization problem (23) is solved using the MMA-scheme 43 which requires gradients of the cost and constraint functions. The sensitivity for the objective and the plastic work constraint functions are summarized below and detailed in the Appendix. The sensitivity of the volume constraint is trivial and is therefore omitted. Since the number of design variables greatly exceeds the number of cost and constraint functions in the optimization problem, we use the adjoint sensitivity scheme presented in the work of Michaleris et al. 22 To derive the sensitivity of g o and g 1 of (23), we consider the sensitivity of the generalized function . In the above (31), it is understood that u = u(f ( )) and w = w(f ( )), and thatf ( ) = (̃(̃( ))) .
In the following, we first calculate the sensitivity with respect to the filtered nodal volume fractionsD after which we can evaluate the sensitivity DΘ∕D as described in Reference 32. In (33), the operator D(⋅)/D(⋆) denotes the total derivative of (⋅) with respect to (⋆). Details of the adjoint sensitivity scheme of (33) are provided in the Appendix.

NUMERICAL EXAMPLES
In this section, we present numerical examples where structures are designed for maximum stiffness while constraining their volume and maximum specific plastic work. All optimization problems are solved with the maximum allowed volume V = 0.35Vol(Ω o ) and are initiated with the uniform design = 0.35. The design is updated via the MMA using the standard parameter setting as presented in the work of Svanberg, 44 but to improve the optimization convergence, we follow Verbart et al. 45 and explicitly control the moving asymptotes, that is, the MMA approximation. Specifically, we restrict the maximum distance between an asymptote and the design variable to 0.1. To obtain good scaling of the MMA subproblem, the objective function is normalized by its initial design value, that is, such that g I=1 o = 1, unless otherwise stated.
To decrease the risk of obtaining local minima, we apply continuation schemes for the penalty and thresholding parameters. With increased penalty and thresholding parameters, the nonconvexity of the optimization problem increases, and as such, it is preferable to start from low parameter values and gradually increase them as the optimization progresses. In line with recent performance studies of continuation schemes by Li and Khandelwahl, 37 we apply penalty continuation and subsequently thresholding continuation (see also, e.g., the work of Amir 18 ). The exponent p E is initially equated to 1 and is increased every 10th optimization iteration in increments of 0.1 until the final value p E = 3. The Heaviside parameter H is likewise initially equated to 1, and once the penalty continuation terminates at p E = 3, it is increased such that H ← 1.1 ⋅ H every 10th optimization iteration until it first exceeds 10, whereupon it is equated to max H = 10 and the continuation process terminates. In the subsequent iterations the minimum distance between an asymptote and the design variable is decreased from the default value 10 −2 (cf. the work of Svanberg 44 ) to 10 −4 to avoid oscillations of the p-norm constraint. The optimization terminates after 500 MMA iterations.
For the discretization of the designed plate structure, we use fully integrated 8-node tri-linear cubic elements with 1 mm thickness. We model plane stress conditions such that the plate is traction free on the lateral x − y faces and apply displacement and traction boundary conditions over the x − z and y − z faces.
The material parameters for the solid material are denoted in Table 1. In the limit of small strains, these material parameters give the Young's modulus E ≈ 100 GPa and Poisson's ratio ≈ 0.3. Uniaxial Cauchy stress versus engineering strain curves for this material are given in Figure 1.
The filter bulk length scale is l o = r∕(2 √ 3) mm with r = 4l e , where l e is the element side length. Following Wallin et al., 33 we specify l s = 0 on clamped faces, symmetry faces and on surfaces subjected to traction boundary conditions, and additionally on regions within a distance of 2r from the surfaces subjected to nonzero prescribed displacement. For the remainder of the surface, we specify l s = l o .

L-bracket
We first demonstrate the proposed optimization algorithm by designing the well-known L-bracket, cf. Figure 2. The L-bracket side length is L = 150 mm, and it is subject to a prescribed uniform displacementū y = −1.5 mm and zero traction t x = 0 over the upper 6 mm of the right surface, zero displacement over the top face, and zero traction over the remaining faces. The domain is discretized with 14400 elements of 1 × 1 × 1 mm 3 size. We limit the maximum specific plastic work to be close to zero, by equating g 1 = yo ⋅ 10 −4 . Three equally sized displacement load increments are applied to the bracket. The optimized design presented in Figure 3(A) clearly shows how the material in the vicinity of the reentrant corner has been removed to alleviate the stress concentration and hence plastic work. The von Mises stress distribution v = √ 3∕2 ‖ ‖ dev‖ ‖ in this corner region is uniform, cf. Figure 3(B). This optimization example also shows how the augmented PDE filter and the Heaviside thresholding generate crisp designs. For comparison purposes, the optimization is also performed using the conventional PDE filter with l s = 0. The optimized l s = 0 design is shown in Figure 4(A) and is similar to the l s = l o Figure 3 design. The effect of the l s choice on the domain boundaries is shown in Figure 5 where zoomed-in contours of = 0.5 for the two optimized L-brackets are shown. The figure shows less artificial "sticking" to the domain boundary for the augmented PDE filter (cf. Figure 5A) versus the conventional PDE filter (cf. Figure 5B). At the reentrant corner (x, y) = (60, 60), the filtered volume fraction is̃≈ 0.6 for the l s = l o Figure 3 design and hence the physical volume fraction is greater than 0.5, cf. The convergence histories for the objective and maximum specific plastic work constraint functions for the Figure 3 design are plotted in Figure 6. After approximately 30 optimization iterations, a stable convergence towards a minimum objective function value (blue curve) with active g 1 -constraint (red curve) is seen. The spikes in the objective function history are due to the continuation schemes. It is also observed that the optimization has not completely converged after 500 iterations, however, as shown by the evolution of the topology (cf. Figure 6), no more topological changes occur after iteration 300, rather only shape changes and reductions in the size of the intermediate region occur.
Next, we demonstrate the influence of the maximum allowed specific plastic work on the optimized L-bracket designs. Figure 7 shows four different designs, optimized using the constraint tolerances g 1 = 5.0 yo ⋅ 10 −4 ( Figure 7A),  Figure 7B), g 1 = 1.0 yo ⋅ 10 −4 ( Figures 7C and 3), and g 1 = 0.5 yo ⋅ 10 −4 ( Figure 7D). In all designs, the reentrant corner is replaced by a "fillet". The size of the fillet is increased as g 1 is decreased, as expected. Perhaps less expected is the fact that the topologies of the designs are identical. The objective function and true maximum specific plastic work values for the Figure 7 designs show expected decreases of stiffness and peak plastic work as the constraint is more strictly enforced, cf. Table 2.
Finally, we demonstrate the effect of the final penalty exponent value p E on the optimized L-bracket designs by performing four different optimizations. The resulting optimized designs are shown in Figure 8

Double-clamped plate
This double clamped plate optimization example further demonstrates the influence of the maximum specific plastic work constraint tolerance value g 1 on the optimized design, cf. Figure 9. The design domain of the plate has the height and width 0.5L and 2L, respectively, with L = 150 mm. A uniform transverse displacementū y = 0.0035L = 0.525 mm and zero traction t x = 0 are prescribed over the center 30 mm region of the plate's top x − z face. Due to symmetry, only the left half of the domain is analyzed; this 150 × 75 × 1 mm 3 domain is discretized by 150 × 75 × 1 = 11250 elements of 1 mm side length.

F I G U R E 8 L-bracket designs for different penalty exponent values p E
F I G U R E 9 Double-clamped plate First, we perform an optimization assuming elastic material response to generate a reference design whereby we increase the initial yield stress tõy o = 10 10 ⋅ yo . In this way, the g 1 constraint is not active and effectively we maximize the end-compliance under the volume constraint g 2 . Figure 10(A) shows both the optimized topology (left subplot) and the von Mises stress distribution (right subplot) for this design. The stress plot uses a colormap ranging from 0 to 200 MPa to indicate the regions of the optimized structures that exceeds the yo = 200 MPa yield stress. The topology consists of straight bars connecting the highly stressed loading and support regions; the highest von Mises stress is located at the bottom of the clamped region.
Next, we perform four different optimizations with elastoplastic material response, that is, with yo = 200 MPa, and with the constraint tolerance values g 1 = 4.0 yo ⋅ 10 −4 ( Figure 10B), g 1 = 2.0 yo ⋅ 10 −4 ( Figure 10C), g 1 = 1.5 yo ⋅ 10 −4  Figure 10D), and g 1 = 0.75 yo ⋅ 10 −4 ( Figure 10E). In all cases, the objective function value is normalized by the g I=1 o value obtained from the elastic reference optimization. The four resulting optimized structures all differ from the Figure 10(A) elastic reference structure, as more material is placed in the bottom center and in the upper support regions. As shown in the right subplots in Figure 10(B-E), this results in smoother stress fields as compared with the elastic reference structure. The optimized topology of the g 1 = 4.0 yo ⋅ 10 −4 design is identical to the elastic reference structure; however, to decrease the peak plastic work, material is redistributed as previously stated. As a result, the stress concentration and consequently plasticity is avoided in the loading and support regions.

F I G U R E 10
The topology for the g 1 = 2.0 yo ⋅ 10 −4 design differs from the elastic reference structure and the g 1 = 4.0 yo ⋅ 10 −4 design as an additional supporting structural member connects the two bars closest to the supporting region. In addition, this g 1 = 2.0 yo ⋅ 10 −4 design contains more curved bars with smooth edges which presumably reduce the maximum specific plastic work. Finally, the g 1 = 1.5 yo ⋅ 10 −4 and g 1 = 0.75 yo ⋅ 10 −4 designs result in yet other topologies, wherein less material is devoted to the lower support regions and the bars exhibit greater curvatures.
Similar to the L-bracket optimization, the stiffness and the peak specific plastic work for the Figure 10(B-E) designs decreases as the constraint is more strictly enforced, cf. Table 3.
To give a fairer comparison of the objective function values, the optimized elastic reference structure 10(A) is reanalyzed using the elastoplastic material modeling, i.e., using yo = 200 MPa, whereupon the stiffness decreases from 1.147 to 1.032, which is less than that of the g 1 = 4.0 yo ⋅ 10 −4 and g 1 = 2.0 yo ⋅ 10 −4 designs, and moreover, the maximum specific plastic work is much higher than the other designs, cf. Table 3.
In Figure 11, we plot the convergence histories for the objective and maximum plastic work constraint functions for the structures presented in Figure 10. The figures show that the objective function has been maximized for all designs (A-E) (cf. the top figure), and that the constraint g 1 is active for the four Figure 10 (B-E) designs (cf. the bottom figure). Some oscillations occur in the objective function during the initial stage of the optimization for Figure 10 (B-E) designs, which are most likely due to the stricter constraints. It is also seen that for the first 200 optimization iterations, while the penalty exponent continuation is performed, the compliance and consequently −g o decreases every 10th iteration. Once the penalty continuation is completed, stable convergence towards the maximum compliance with an active g 1 constraint is observed for all of the designs.

Varying load paths: Double L-bracket
In this double L-bracket example we demonstrate the importance of the load path on the optimized structures, cf. Figure 12. The height of the bracket is L = 150 mm. Two uniform transverse displacement boundary conditions, respectively,ū y1 = 1.25 mm and t x1 = 0 andū y2 = 1.25 mm and t x2 = 0, are applied over the 0.048L wide regions of the outer y − z left and right faces. We consider two different prescribed displacement paths, shown in Figure 13. Path #1 ramps up both prescribed displacementsū y1 andū y2 simultaneously. For load path #2,ū y2 is first incremented until it reaches 75% of its final value, whileū y1 is fixed at 0 mm. Then, bothū y2 andū y1 are increased to their final values 1.25 mm. Notably, both paths prescribe identical final boundary conditions. In an elastic setting, the different paths generate identical design results, this is not the case for optimizations based on elastoplasticity. The domain is discretized by 13750 elements of side length 1.2 mm and thickness 1 mm. We set l s = l o on all surfaces, except the clamped and the outer left and right surfaces, where l s = 0 is specified.
To maximize the compliances C (N) left and C (N) right at the terminal load step N, at the left and right ends of the bracket, the net end-compliance C (N) = C (N) left + C (N) right may be maximized. However, such maximization allows either of the left or right reaction tractions to take nearly zero values. Therefore, we instead minimize to ensure nonvanishing reaction tractions. and maximum specific plastic work constraint histories (bottom) for the optimized topologies shown in Figure 10 F I G U R E 12 Double L-bracket Similar to the previous example, we first perform an optimization using elastic material response to be used as a reference design. For this case, we use load path #1. Figure 14(A) shows that the optimized topology is nearly symmetric with respect to the y-axis. This is expected because the terminal load case is nearly asymmetric about the y-axis, and because the hyperelastic response is similar (but not identical) in tension and compression for small strains.
Next, we constrain the specific plastic work, and perform optimizations for load paths #1 and #2. Based on our numerical experiments, the final value of the penalty exponent p E is increased in these optimizations to the terminal value p E = 4. The continuation for H is still initiated once p E > 3. The maximum specific plastic work constraint tolerance is The double L-bracket optimized designs corresponding to load paths #1 and #2 are shown in Figure 15(A,B), respectively. Both designs exhibit smooth boundaries in the vicinity of the reentrant corners, which is not the case for the elastic reference design, cf. Figure 14. The von Mises stress is clearly reduced in these designs, cf. Figures 14(B) and 16. The maximum specific plastic work values are max(ŵ p ) = 0.99 g 1 and max(ŵ p ) = 1.01 g 1 for load paths #1 and #2, respectively, that is, both optimizations reach the target value within 1%. Due to the different load paths, it is not surprising that the compliance between the two designs varies; the objective function values are g o = 0.536g I=1 o and g o = 0.511g I=1 o for the load path #1 and #2 designs, respectively.
Despite the fact that the terminal load states are identical, and maximum values and their corresponding locations of the specific plastic work are nearly identical, the optimized topologies in Figure 15(A,B) differ. Similar to the elastic reference design, the structure optimized for load path #1 is nearly symmetric (cf. Figure 15A) whereas the structure optimized for load path #2 clearly is not symmetric (cf. Figure 15B). The material distribution in the center region of the designs also differs. For load path #1, the center region contains a V-shaped boundary with a sharp corner at x = 0, whereas for the load path #2, the center region is smooth to avoid high shear strains and consequently high specific plastic work. To show this, and to clearly show the performance differences between the designs, we analyze them with interchanged load paths, that is, the structures optimized for load paths #1 and #2 are subjected to load paths #2 and  Figure 17.
Again we observe that̃≈ 0.6 at the corners (x, y) = (±30, 60) which gives rise to reentrant corners. Figure 17(A) shows that the structure optimized for load path #1 exhibits the maximum specific plastic work at the point of the V with value max(ŵ p ) = 21.97 g 1 when subjected to load path #2, i.e., its value is approximately 22 times higher than the allowable value. On the other hand, the structure optimized for load path #2 also performs well when it subjected to load path #1; the peak value of max(ŵ p ) = 1.17 g 1 is located at the left reentrant corner, cf. Figure 17(B).

CONCLUSIONS
A topology optimization framework based on elastoplasticity for maximizing stiffness while limiting the maximum specific plastic work is presented. clearly demonstrate the influence of different levels of allowed peak specific plastic work on the optimized structures and how the load path significantly influences the optimized structures. The optimization framework is based on a finite strain isotropic hardening kinematics and constitutive model and limits the specific plastic work during a single load cycle. For future work, the constitutive model could be extended to involve kinematic hardening and account for the Bauschinger effect. In this way, the fatigue life may be predicted for cyclic loading with the techniques presented in this article.

DATA AVAILABILITY STATEMENT
Relevant parts of the data that support the findings of this study are available from the corresponding author upon reasonable request.

APPENDIX A. ADJOINT SENSITIVITY
This section presents the adjoint sensitivity scheme for obtaining DΘ∕D . In the adjoint method, we augment the function Θ by adjoint variables f and and the implicit residuals (12) and (14) Dw (n−1) D̃+ where C is the collection of all Gauss point C i residuals, and we recall that F ext = 0, that is, R = F int , and note that D(⋅)∕Du p = 0. We also note that, since the implicit residuals (12) and (14) hold,Θ = Θ and DΘ∕D̃= DΘ∕D̃. Ultimately we assign the, heretofore arbitrary, adjoint variables such that ) T (1) ) T

(N−n)
] . (A5) The annihilation of the implicit derivatives starts by considering (A4): we set in order to annihilate the Du (N) f ∕D̃-term. Next, we set that is, to annihilate the Dw (N) ∕D̃-term of (A4). By inserting (A8) into (A6) we find or equivalently (1) where ) T is the (transpose of the) partitioned tangent operator at time step N and