Load step reduction for adjoint sensitivity analysis of finite strain elastoplasticity

In this paper, load step reduction techniques are investigated for adjoint sensitivity analysis of path-dependent nonlinear finite element systems. In particular, the focus is on finite strain elastoplasticity with typical hardening models. The aim is to reduce the computational cost in the adjoint sensitivity implementation. The adjoint sensitivity formulation is derived with the multiplicative decomposition of deformation gradient, which is applicable to finite strain elastoplasticity. Two properties of adjoint variables are investigated and theoretically proved under certain prerequisites. Based on these properties, load step reduction rules in the sensitivity analysis are discussed. The efficiency of the load step reduction and the applicability to isotropic hardening and kinematic hardening models are numerically demonstrated. Examples include a small-scale cantilever beam structure and a large-scale conrod structure under huge plastic deformations.


Introduction
Shape optimization plays an important role in industrial structure design. Non-parametric shape optimization, which dates back to the work by Zienkiewicz and Campbell (1973), employs nodal coordinates in a finite element system as design variables. Due to a detailed description of structure shape and enlarged design space, non-parametric shape optimization shows great success in engineering applications (Clausen and Pedersen 2006;Furbatto et al. 2009;Böhm and Clausen 2012;Shimoda et al. 2019).
Both gradient-based and gradientless methods have been investigated to address the non-parametric shape optimization with linear finite element systems (Schnack 1979;Meske et al. 2005;Meske 2007; Le et al. 2011). In the recent decade, the focus is more on optimization of path-dependent nonlinear problems (Pedersen et al. 2017). So far, in almost all the literature about structural optimization considering elastoplasticity, the investigations are based on small strain theory where an additive decomposition of the total strain applies Bogomolny and Amir 2012;Amir 2017;Li and Khandelwal 2017;.
Exceptions are only few in recent years (Wallin et al. 2016;Ivarsson et al. 2018).
One major challenge of gradient-based optimizations for path-dependent finite element systems is the significantly increased computational cost in the sensitivity analysis. Due to path dependency, mechanical responses and hence sensitivities of these responses at one load step are determined not only by the current load condition but also by previous load histories (Park and Choi 1999;Choi and Kim 2005). Direct sensitivity analysis, also called direct differentiation method, must follow an incremental sensitivity analysis procedure (Spivey and Tortorelli 1994;Chattopadhyay and Guo 1995;Kleiber and Kowalczyk 1996;Kim et al. 2000;Wisniewski et al. 2003;Gu et al. 2009). A unified framework has been presented on how to formulate sensitivity with adjoint variable method for a wide range of path-dependent system behaviors (Michaleris et al. 1994;Alberdi et al. 2018). It shows that the adjoint variable method should follow a backward solution 1 3 19 Page 2 of 23 procedure (Michaleris et al. 1994;Lee 1999;Maute et al. 2003;Chung et al. 2003;Alberdi et al. 2018). No matter the direct approach or the adjoint approach is employed, the computational effort of pathdependent sensitivity analysis grows in proportional to the number of load steps. Therefore, techniques to reduce the computational cost in the sensitivity analysis are highly demanded.
For this purpose, sensitivity reanalysis in the frame of independent coefficients strategy has been suggested, which employs local modification of a large-scale structure to avoid repeated solutions of full finite element analysis (Liu and Wang 2008). However, this has been demonstrated to be effective only for linear structural systems. Various strategies have been investigated to reduce the number of load steps in the small strain elastoplasticity sensitivities analysis. Numerical study of a structure with single tetrahedral element and a two-bar truss structure under bilinear isotropic hardening model show that sensitivities may only need to be calculated at the last load step given a monotonic load history (Köbler 2015). The sensitivities are only necessary to calculate at the unloading and reloading points given a cyclic load history (Cardoso 2005). It has been theoretically proved that, for small strain elastoplasticity, intermediate elastic load steps could be skipped in the sensitivity analysis, and the former ones among consecutive plastic steps can also be skipped if directions of flow vectors are unchanged (Wang et al. 2017).
The present publication extends the load step reduction strategies from small strain case to finite strain elastoplasticity. It also demonstrates the applicability to various hardening models. This extension is not straightforward for several reasons. Firstly, the multiplicative decomposition of deformation gradient in finite strain theory leads to a more complicated formulation of adjoint sensitivity. Secondly, as will be shown, the prerequisites for load step reduction in finite strain problems are stricter than in small strain cases. The extent to which the load step reduction is possible should be investigated. Thirdly, the accuracy of sensitivities with reduced load steps must be demonstrated for large strain problems.
The paper is organized as follows: Sect. 2 introduces the finite strain elastoplastic analysis procedure which is employed in this study. In Sect. 3, adjoint sensitivity formulation is presented. In Sect. 4, the properties of adjoint variables and load step reduction method for elastic steps and plastic steps are proposed and theoretically investigated. In Sect. 5, the efficiency and accuracy of the load step reduction method is demonstrated with a solid beam structure under severe deformations and a large-scale connecting rod example. In Sect. 2 to Sect. 5, an isotropic hardening model is assumed to reduce complexity of the discussion. In Sect. 6, the load step reduction method is extended to combined hardening and kinematic hardening cases. Additionally, a short discussion on the applicability to nonlinear elasticity and multilinear elastoplasticity are presented. Finally, the conclusions are drawn in Sect. 7.

Finite strain elastoplastic analysis
In this section, the finite strain elastoplastic analysis is introduced based on the total Lagrangian formulation and logarithmic strain measure. Following Lee's multiplicative decomposition (Lee and Liu 1967;Lee 1969), the total deformation gradient t 0 X is expressed as the multiplication of an elastic deformation gradient t i X e and a plastic deformation gradient i 0 X p : where 0 refers to the undeformed configuration and i refers to the intermediate stress-free configuration. The elastic deformation gradient can be further decomposed by right polar decomposition: where R e is a unitary matrix, called the elastic rotation tensor, and U e is a symmetric positive-definite matrix, called the elastic right stretch tensor. The logarithmic strain, also known as Hencky strain, is often employed as a proper strain measure for finite strain problem. The elastic logarithmic strain tensor is defined by: The rotated Kirchhoff stress tensor , which is the spatial Kirchhoff stress rotated to the intermediate stress-free configuration by R e (Gabriel and Bathe 1995;Caminero et al. 2011;Neff et al. 2016), is the work-conjugate stress measure to the elastic logarithmic strain: where D e is the elastic constitutive relation matrix. The spatial Kirchhoff stress is eventually obtained by back rotating : The Kirchhoff stress is related to the Cauchy stress through the Jacobian determinant: The internal force can be obtained by integration under deformed volume V or initial volume V 0 Newton-Raphson method is employed in the nonlinear finite element analysis. The flowchart of the method is depicted in Fig. 1. t U is the nodal displacement, t p eqv is equivalent plastic strain, t X p is plastic deformation gradient, t F is external force. The upper left superscript denotes the load step.
The return mapping algorithm is employed to obtain the plastic deformation gradient, equivalent plastic strain and stress tensor at step t + 1 inside each Newton-Raphson iteration. Given the deformation gradient 0 t+1 X at step t + 1, the plastic deformation gradient and equivalent plastic strain at previous step t, the workflow of return mapping algorithm (Eterovic and Bathe 1990;Dvorkin et al. 1994) is summarized in Table 1. It should be noted that, the algorithm works with the rotated stress tensor. The spatial stress is obtained thereafter by back rotating.
As it has been pointed out (Montáns and Bathe 2005;Caminero et al. 2011), Eqs. (16) and (17) can be approximated by This approximation holds for moderately large elastic strains, which is typically fulfilled in metal plasticity. Eterovic and Bathe (1990) also claim that Eq. (19) is exact for isotropic hardening plasticity with associated flow rule or for combined isotropic-kinematic hardening cases where the stress and back stress tensors commute. In view of this, Eq. (18) can be written as When the stopping condition of the return mapping algorithm is satisfied, the following quantities at step t + 1 can be obtained by In Eq. (21), the trial elastic rotation tensor t+1 R e * is used to back rotate the stress tensor to the undeformed configuration. Under the associated flow rule, the incremental plastic stretch Δ p and the trial elastic stress tensor t+1 U e * have the same eigenvectors. Therefore, it can be verified the trial elastic rotation tensor equals the real elastic rotation tensor t+1 R e . Following Eq. (7), the consistent tangent stiffness matrix is where the constitutive relation matrix D is (Simo and Taylor 1985;Crisfield 2000) where The formulation of adjoint sensitivity analysis for finite strain elastoplasticity is derived in this section. The solution procedure for the adjoint variables is also explained. There are two key points in deriving the adjoint formulation. One is the selection of a set of state variables, from which all quantities at each finite element analysis load step can be reconstructed. The selection of state variables is not unique. A proper selection of them will reduce the complexity of the sensitivity formulation. The other key issue is a set of governing equations of the state variables. These equations should be identically equal to zero and have the same number as state variables.
Displacement, stress, strain, plastic strain, equivalent plastic strain and flow vector are typical quantities that to be determined in an elastoplastic analysis. Some of these variables are not independent and can be derived from others. After investigating different combinations, the following four quantities are selected as state variables: displacement, rotated stress tensor, equivalent plastic strain and inverse of plastic deformation gradient The reason of using inversed plastic deformation gradient is it leads to a simpler evaluation of derivatives of the trial rotated stress tensor with Eqs. (8) to (11).
One natural governing equation of the state variables is the equilibrium condition. The residual force is identical to zero at each load step: where s denotes the design variables, and the elastic rotation tensor t R e is a function of t U and t−1 X p−1 .
Other governing equations are found on the element level. According to Eqs. (15) and (20), the rotated stress tensor and the elastic right stretch tensor should follow c. Calculate the trial elastic logarithmic strain and trial elastic stress e * = ln t+1 U e * = 1 2 ln t+1 X e * T t+1 X e * (10) * = D e e * (11) d. Obtain the trial von-Mises equivalent stress where S * is the deviatoric stress of the trial stress * .
2. Check for yield condition. If s * ≤ t Y + E p ⋅ Δ p eqv , then stop. Otherwise, enter return mapping. In the formula, E p is the plastic modulus 3. Return mapping (a) Obtain associated flow vector : The yield and consistency condition is According to Eq. (23), the plastic deformation gradients of two consecutive steps follow The governing Eqs. (32) to (34) on element level define the so-called dependent residual t H (Michaleris et al. 1994), which identically equals zero at each load step: For a system response, it can be expressed as a function of all state variables and design variables where the superscript N is the total number of load steps.
The direct differentiation of the response with respect to a design variable s is To avoid time-consuming evaluations of d t U∕ds and d t V∕ds , two adjoint variable vectors t and t are introduced. They are in the same size as t R and t H respectively. Since t R in Eq. (31) and t H in Eq. (35) are identical to zero at all load steps, adding dot product of t and t R and dot product of t and t H to the response function will not change the value of it, i.e. From Eqs. (31) and (35), the derivatives of t R and t H with respect to the design variable are By taking the total derivative of the response in Eq. (38) with respect to the design variable s, and substituting Eqs. where t K tan is the consistent tangent stiffness matrix following Eq. (24). Since the adjoint variables at step t is dependent on the adjoint variables at the next step t + 1, the solution for the adjoint variables must be carried out backwards from the last load step to the first load step.
19 Page 6 of 23 4 Load step reduction in the adjoint sensitivity analysis According to Eqs. (44) to (46), the major computational effort in the adjoint sensitivity analysis is the backwards solution of adjoint variables. For each load step, a system of linear equations, which is in the same size as the number of degrees of freedom of the underlying finite element model, needs to be solved. Therefore, the computational cost increases in proportional to the total number of load steps. In this section, strategies to reduce the number of load steps that are used in the sensitivity analysis are investigated. Besides the computational cost, the memory cost in implementation may also be challenging if it is not properly handled. The memory cost is measure by number of non-zero values that must be stored at the same time. Following Eqs. (45) and (46), tangent stiffness matrices at the equilibrium point, partial derivatives of residual force and partial derivatives of dependent residual are required in the solutions of adjoint variables. These quantities will be gradually collected after the nonlinear analysis at each step. Due to the backward solution procedure, they must be kept in memory until the last step finite element analysis is finished. Although the stiffness matrices are sparse, the direct storage of all these quantities at all load steps is not the most cost-efficient way.
By looking into the formulation of these quantities in Appendices A and B, the tangent stiffness matrices are fully determined by state variables in Eq. (30). This is also the case for the partial derivative quantities. Therefore, to minimize the memory cost, it is suggested to store only state variables at all load steps. During the sensitivity analysis, intermediate quantities are then retrieved from state variables. The regeneration happens only on element level, where the computational effort is neglectable. However, it is worth mentioning that the storage space is anyway in proportional to the number of load steps in sensitivity analysis. Therefore, strategies to reduce load steps in the sensitivity analysis also contribute to save memory cost in implementation.
Before presenting the load step reduction method, one assumption regarding the responses should be made clear. Given a sequence of load steps L = {1, 2,…, t − 1, t, t + 1,…, N}, the system responses that are investigated in this paper are assumed to be functions of quantities of only the last load step, i.e. it assumes that the system response f satisfies Many responses fulfill these requirements, such as maximum equivalent stress or final displacement. Besides that, many other responses can also be expressed as a composition of such functions, i.e.
where function f k is dependent only on quantities at step k. If this is the case, the sensitivity of individual f k could be calculate first, where the step k is treated as the last load step. And then the sensitivity of g can be obtained eventually by the chain rule.

Load step reduction for elastic steps
In the following, an elastic load step describes a step of finite element analysis, in which all elements behave elastically. Otherwise, if any element behaves plastically in the load step, then the step is called a plastic load step.
The following property of adjoint variables at elastic steps was first found in additive decomposition based small strain case (Wang et al. 2017). This property also holds for multiplicative decomposition based finite strain elastoplasticity.

Property 1. If load step t is an intermediate elastic step (t is not the last step), then the adjoint variable t at this step is zero.
Due to limited space, the proof of this property is presented in Appendix C. The property eventually leads to the following load step reduction rule.
Elastic load step reduction rule For a sequence of load steps L = {1, 2,…, t − 1, t, t + 1,…, N}, if step t is an intermediate elastic load step, then the exact same sensitivity results can be obtained by skipping step t as the load steps contain only S = {1, 2,…, t − 1, t + 1,…, N}.
Mathematically, it means where (49) The subscript S and L on the left denote the step set from which a quantity is calculated. In Eq. (52), the partial derivatives of R and H are based on load steps set S, where the original step t − 1 and step t + 1 become adjacent load steps. The theoretical proof of this rule is presented in Appendix D. It follows that all intermediate elastic load steps can be skipped in the sensitivity analysis. It should be noted that this rule applies to all types of elements, and no accuracy will be lost after the load step reduction.

Load step reduction for plastic steps
Under certain conditions, the adjoint variable equals zero even at plastic steps. The following property is proved in Appendix E.
Physically, the four prerequisites mean that the flow vectors at adjacent load steps should be in the same direction, the elastic rotation tensor should be constant with respect to the plastic deformation gradient, and the incremental plastic strain should have the same principal directions as the accumulated plastic strain. If all these conditions are met, then the following load step reduction rule for plastic steps could be proved. Due to limited space, the proof is presented in Appendix F.

Plastic load step reduction rule
For a sequence of load steps L = {1, 2,…, t − 1, t, t + 1,…, N}, if the prerequisites in Property 2 are all fulfilled, additionally t Δ p and t−1 X p also commute and t R e ∕ t−1 X p−1 = 0 , then the same sensitivity results can be obtained by skipping load step t as the load steps contain only S = {1, 2,…, t − 1, t + 1,…, N}.
In comparison with small strain case (Wang et al. 2017), the condition on incremental plastic flow is stricter. Not only the consecutive two steps should have the same flow vector, but the incremental quantities should also be consistent with the accumulated plastic quantities. The prerequisites of plastic reduction rule could be satisfied only by 1D bar elements when there is no switch between tension and compression in two consecutive plastic steps. For general 2D and 3D elements, it is impossible that these requirements are all fulfilled. Therefore, the following empirical rule is suggested for finite strain elastoplasticity.
Empirical rule If two consecutive plastic steps are in a monotonic loading procedure, and the incremental plastic flow is in close direction to the accumulated plastic flow, then the former step can be skipped in the sensitivity analysis.
Load steps reduction following the empirical rule will not lead to exact sensitivity results for 2D and 3D elements. By doing so, the accuracy of sensitivities needs to be verified. Its applicability in practice and influence on the sensitivity results will be investigated through numerical examples in the next sections.
It must be made clear that, the proposed reduction rules are derived on the structural level. It means that, as defined at the beginning of Sect. 4.1, the whole structure is treated as plastic even if only just one element behaves plastically. If one plastic element is not in monotonic loading, then the load step for the whole structure cannot be skipped. This strategy applies well for shape optimization, where it is rare that abnormality occurs only in one or a few local elements.
To apply in topology optimization, the load step reduction may be limited by the complicated behaviors of a few or even just one single intermediate element. This doesn't mean the proposed scheme cannot by applied to topology optimization. In such cases, the proposed reduction rules provide guidance on which load steps could be skipped in the sensitivity analysis. A further reduction, however, may still be possible by investigating adjoint variables on the element level. The element level load step reduction is an open issue to investigate. On the other hand, there are already studies on how to avoid abnormality of intermediate element plasticity in topology optimization. One typical solution is to choose separate penalization exponents and lower bounds for stiffness and yield properties in the SIMP interpolation (Maute et al. 1998 (Amir 2017;Zhang et al. 2017). The combination of these techniques and load step reduction rules in topology optimization may be interesting for further investigation.

Solid beam under severe bending
In this section, the load step reduction rules are verified through a cantilever beam structure meshed with 3D tetrahedral elements. The finite element model is depicted in Fig. 2. One end of the structure is fixed. Two external forces in horizontal x-direction and in vertical y-direction are applied on the free end of the structure simultaneously. A bilinear isotropic hardening material is assumed. The Young's modulus is 210 GPa, plastic modulus is 50 GPa, initial yield stress is 235 MPa and Poisson's ratio equals 0.3. There are 11 load steps as described in Fig. 3. The vertical load pointing downwards increases in the first three steps and then decreases gradually. The horizontal load in x direction increases throughout the procedure. According to the nonlinear finite element analysis, all the load steps are plastic steps.
The maximum equivalent plastic strain at the last step is 25.2%. The contour of the equivalent plastic strain is presented in both initial and deformed configurations in Fig. 4. Areas where equivalent plastic strain is larger than 5% are depicted in red. The deformation shows that the beam structure is severely bended under the given loads.
In the sensitivity analysis, the design variables are the vertical coordinates of center nodes on the bottom surface of the beam. These nodes are indicated by red dots in Fig. 2. The maximum equivalent plastic strain at the fixed end and average vertical displacement at the free end are defined as two system responses.
Since the structure is in an increasingly bending procedure under the given loads, the plastic flows are all in close directions. According to the empirical rule, all the intermediate load steps in monotonic loading stage could be skipped. Therefore, only step 3, where the vertical load turns from increasing to decreasing, and the last load step must be included in the sensitivity analysis. The sensitivity results using only step 3 and step 11 are compared with sensitivities using all load steps in Figs. 5 and 6. In subfigures (a), the values of sensitivities are plotted with respect to the position of design nodes. Subfigures (b) present the percentage error in term of relative value, which is defined by In subfigures (a), sensitivities obtained by central finite differencing scheme with a perturbation size of 10 −5 mm are also presented. Both Figs. 5a and 6a show that, the adjoint sensitivity results match well with finite defencing results. Therefore, the adjoint sensitivity analysis procedure is validated.
Focusing on the adjoint sensitivity results, the results with reduced load steps match well with that using all load steps. The percentage errors at most of the design variables are small. It increases significantly as closing to the free end of the beam. This increase in error is partially because the sensitivity at the free end is close to zero. The relative error is calculated with a small denominator and hence is amplified. Filtering out sensitivities whose absolute value is smaller than 1% of the maximum sensitivity, the sensitivity errors for both responses are summarized in Table 2. The average errors are smaller than 3%, which demonstrate a good match of sensitivities when the load steps are reduced. Hence the empirical rule applies for this example.

Solid beam under severe bending and twisting
One of the prerequisites of the empirical rule is incremental plastic flow should be in close direction to the accumulated plastic flow. Due to the vagueness of this statement, it is worth presenting a case to show when the load steps can't be further reduced even under a monotonic loading procedure.
In this example, the same cantilever beam model as in Sect. 5.1 is employed. The horizontal and vertical load history are plotted in Fig. 7, and there are 20 load steps in total. Besides these two forces, a constant force of 2 kN is applied at the free end in horizontal Z direction. All the load steps are plastic steps.
The deformation of the beam at representative load steps are depicted in Fig. 8 with contour of equivalent plastic strain. The red color represents the area with equivalent plastic strain larger than 5%. The maximum equivalent plastic strain at final step is 36.3%. The deformations show that the beam slightly bends out of the X-Y plane at the beginning. After large enough plastic strain is accumulated near the fixed end of beam, the out-of-plane load in Z direction gradually causes the beam to twist.
According to the empirical rule, the turning points of monotonic load procedures should be included in the sensitivity analysis. These are step 3, step 13, and step 20. Under twisting of the beam, the principal directions of stress tensor will change significantly. It leads to the change in direction of the associated plastic flow. Therefore, the prerequisite of   the empirical rule is violated. Following the empirical rule, from step 15 on all the steps must be included in the sensitivity analysis although they are in a monotonic loading procedure. It means that the load steps in the sensitivity analysis can only be reduced to step 3, step 13, and steps from 15 to 20 as depicted by pentagrams in Fig. 7. Take the vertical displacement at the free end as the response function. The sensitivity results with all load steps and the reduced load steps are compared in Fig. 9. It shows a good match between these two results.
To demonstrate the necessity of step 15 to step 19 in the sensitivity analysis, several try-outs are also presented in Fig. 9. In each of these results, the same reduced load steps are used in the sensitivity analysis except one step between 15 and 19 is additionally skipped. It shows that, if any step between 15 and 19 is skipped in the sensitivity analysis, the result will have significant errors. This example shows that, to follow the empirical rule, good engineering judgment may be needed in practical applications. A solution to avoid the subjectivity could be quantification of the empirical rule. This is also important for implementing the empirical rule into a non-intrusive optimization procedure because there are cases where structural behaviors change essentially between iterations and hence an automatic load steps reselection is required. How to properly quantify the conditions in the empirical rule and set criteria for selection are still open issues to be investigated.

Demonstration with a connecting rod structure
In this section, the applicability of load step reduction rules is demonstrated through a connecting rod example under cyclic load history. The finite element model of a typical connecting rod (or called conrod) in an internal combustion piston engine is presented in Fig. 10a. The design variables are x-coordinates of 223 nodes which are highlighted by red dots in Fig. 10b. These nodes lie on the outer surface between the small end and big end of the conrod. A bilinear elastoplastic material is assumed. Material properties are listed in Table 3.
The rod bolts of the structure are fixed, and a horizontal force (Fx) and a vertical force (Fy) are applied uniformly on the inner surface of the small end. The load history is   Fig. 11. The forces are generated by the internal pressure on the piston in a four-stroke cycle. The magnitudes of the forces are artificially enlarged to represent an extreme load case during engine failure. In the load cycle, the structure is purely compressed in y-direction at step 1. From step 2 to step 7, the compression load increases. At the same time, a horizontal force in x-direction is gradually applied which causes bending of the rod. From step 8 to step 12, the x-direction force is gradually unloaded to zero and the y-direction force is partially unloaded. The horizontal load increases in negative x-direction at step 13 to step 16 accompanied by a slight increase of the compression load. From step 17 to step 19, the horizontal force is fully unloaded to zero while the y-direction force gradually increases to the same level as in step 1.
The contour of the equivalent plastic strain at the last load step is depicted in both initial and deformed configurations in Fig. 12. Areas where equivalent plastic strain is larger than 10% are depicted in red. The maximum equivalent plastic strain is 22.3%, which is close to strain at break of typical steel material. The maximum value point lies on the outer surface near the big end. The average equivalent plastic strain of 106 elements around this point is defined as one system response. The average x-displacement and average y-displacement of the small end at the last step are defined as other two responses.
According to the nonlinear analysis, steps 2 to step 7 and step 16 are plastic, other load steps are elastic. Following the elastic load step reduction rule, all elastic load steps are skipped in the sensitivity analysis except step 19. Based on the empirical rule, plastic steps 2 to 6 are also skipped since they lie in a monotonic load procedure. The three load steps that must be involved in the sensitivity analysis are highlighted by pentagrams in Fig. 11 and listed in Table 4.
In Fig. 13, the contour of sensitivities calculated with the reduced load steps are compared with those using all load steps. The results match very well for all three responses.
The percentage errors of sensitivities are evaluated and summarized in Table 5. To eliminate the effect of small denominator, sensitivities whose absolute value is smaller than 1% of the maximum sensitivity are ignored. It shows that the average percentage errors of sensitivities for all three responses are less than 10%. Hence, the empirical rule applies well in this example.

Efficiency in terms of computational time
Up to here, the efficiency of load step reduction rules is measured by the number of load steps used in the sensitivity analysis. It is assumed that the computational time of sensitivity analysis will be reduced proportionally to the number of load steps. This assumption is verified explicitly in this section.
Besides the number of load steps, the computational time of sensitivity analysis is also dependent on many other factors including degrees of freedom of the underlying FEmodel, number of responses, number of design variables, computing environment, coding efficiency, etc. To focus on  the factor of number of load steps, the ratio of sensitivity analysis time using reduced load steps to the time using all load steps could be evaluated. In Fig. 14, the ratio of time of sensitivity analysis is presented in the vertical axis. The horizontal axis presents the ratio of number of reduced load steps to the number of all load steps. There are four points in the figure, representing the results of four numerical examples in this paper. For better organization of the paper, the example presented in a later Sect. 6.1 is also included here.
It clearly shows the reduction of sensitivity analysis time is proportional to the reduction of number of load steps with proportionality constant close to one. Although there are not enough data from statistics point of view, the linear relation could still be concluded, which is also in accordance with intuition.

Influence on the computational cost of optimization
As shown in previous examples, the accuracy of sensitivity will be lost to some extent while reducing the number of plastic load steps. It naturally leads to the question, how the loss of sensitivity accuracy will influence the optimization? Will the gain in the sensitivity analysis be consumed Fig. 13 Comparison of sensitivity results of three system responses. In each subfigure, the result with reduced load steps is depicted on the left and the result using all load steps is presented on the right  These are open question to be investigated. The following three aspects should be paid attention to while addressing the questions. They show the answers are not straightforward. Firstly, will and how the slight loss of accuracy jeopardizes information provided to an optimization algorithm? In gradient-based optimizations, sensitivities are usually used to generate a search direction in which to improve system responses. As demonstrated, following the empirical rule, the degree of inaccuracy of the sensitivities are very limited. There is no obvious change of sensitivity magnitudes, and the spatial distributions are also the same. Therefore, it is sufficient to believe that the directional information is still well preserved after proper load step reduction. The influence on the optimization process is thus minimized.
Secondly, on the nature of optimization algorithm, in which framework the sensitivities are utilized leads to algorithm-dependent answers to aforementioned questions. In the literature, efforts are taken to reduce computational cost by improving optimization procedures and algorithms. An accelerated gradient algorithm has been proposed (Arouri and Sayyafzadeh 2020), which is less sensitive to the gradient approximation accuracy than the steepest descent algorithm. Switching back and forth between accurate sensitivity and Broyden approximation (Li et al. 2007), both trust region and SQP based minmax algorithm may take even fewer optimization iterations to converge to an optimal point. The method of moving asymptotes (MMA) is a popular algorithm for shape and topology optimization (Svanberg 1987). An interesting study presented by Amir (2021) shows that, aiming at delivering viable engineering designs within limited number of function evaluations, inexact design sensitivities do not lead MMA algorithm to a significantly inferior solution. In short, the choice of algorithms will lead to different requirement on sensitivity accuracy.
Last but not least, how to use sensitivities in shape optimization procedure is also a topic. Two general problems in gradient-based shape optimization are the mesh-dependency of sensitivities and non-smooth shapes resulted by noisy sensitivity fields. Successful techniques to address these issues are sensitivity filters (Le et al. 2011;Stück and Rung 2011;Sigmund and Maute 2012;Bletzinger 2014), sensitivity weighting (Kiendl et al. 2014) and vertex-morphing (Hojjat et al. 2014). The common point of these methods is that the shape sensitivities are smoothed before they are used for design update. These intentional modifications will unavoidably lead to discrepancy between post-processed sensitivities and raw sensitivities. The requirement on the accuracy of raw sensitivities is loosened to some extent under these circumstances.
Generally speaking, how the loss of raw sensitivity accuracy influences the optimization performance is a case-by-case question. This paper focuses on the reduction of computational cost in sensitivity analysis with least loss of raw accuracy. In practical optimization applications, the requirement on the sensitivity accuracy should be considered in a comprehensive way.
6 Extension to more general constitutive models

Extension to kinematic and combined hardening model
In this section, the load step reduction technique is extended to kinematic hardening and mixed hardening case. The adjoint sensitivity analysis and the reduction rules are briefly discussed under these models. The applicability of the empirical rule is demonstrated through a conrod example.
In the adjoint sensitivity analysis with kinematic or mixed hardening, the major differences to isotropic hardening case lie in the state variables and dependent residual. They are caused by the introduction of the back stress and the hardening ratio. Denote the back stress under stress-free configuration by t b . For bilinear mixed hardening elastoplasticity, the yield surface is where is the hardening ratio lies in the range from 0 to 1. If is equal to 1, it describes a pure isotropic hardening behavior, and with equals 0, it describes a pure kinematic hardening model. The back stress of two consecutive steps follows where Except the state variables in Eq. (30), the back stress is taken as an additional one, i.e.
Correspondingly, Eq. (55) is an additional governing equation in the dependent residual: There are no other changes in the adjoint sensitivity analysis procedure. The two properties and two load steps reduction rules presented in Sect. 4.1 and 4.2 can be proved following the same procedures as in Appendices C to F. Therefore, the same empirical rule as in Sect. 4.2 is also proposed to reduce plastic load steps in the sensitivity analysis.
The following example demonstrates the applicability of the empirical rule for kinematic hardening elastoplasticity. In this example, the same conrod model as in Sect. 5.3 is used, including the structure, boundary conditions and design nodes. The material properties are also the same as in Table 3, except the hardening model is pure kinematic, i.e. = 0 . The load history on the small end of the conrod is depicted in Fig. 15, which is the same as in Fig. 11. However, it should be noted that, all the load steps behave plastically with kinematic hardening material. The contour of the equivalent plastic strain at the turning steps of monotonic load procedures are depicted in Fig. 16. Areas where equivalent plastic strain is larger than 20% are depicted in red. The maximum equivalent plastic strain at the final step is 73%. It shows that the structure experiences severe back and forth bending under the given load history.
According to the empirical rule, all the turning steps which are plastic should be involved in the sensitivity analysis. These are step 7, step 12, step 16 and step 19. In comparison with isotropic hardening case, the step 12 is additional. This is because step 12 is a plastic step under kinematic hardening model. The sensitivity results of average equivalent plastic strain response, which is defined in Sect. 5.3, is presented in Fig. 17. It shows that the sensitivity analysis with reduced load steps match well with results using all load steps. Excluding sensitivities whose absolute value is smaller than 1% of the largest value, the maximum percentage error is 41% with an average error of only 3.6%. The point with maximum error is also identified in Fig. 17. It shows that the large error is due to the small denominator in calculating with Eq. (53). Therefore, the applicability of empirical rule to kinematic hardening model is demonstrated.

Extension to nonlinear elasticity and multilinear plasticity
The deductions and examples in this paper are all based on bilinear elastoplasticity. However, the applicability can be extended to nonlinear elasticity and multilinear plasticity without difficulties.
To extend for nonlinear elasticity, the change is in the dependency of elastic modulus on the total strain. The elastic constitutive relation D e is not constant, but dependent on the strain at each load step: Therefore, the dependent residual is formulated as For multilinear plasticity, the change is in the relation of yield strength to the equivalent plastic strain. The formulation of dependent residual is the same as in Eq. (35). The only change appears when the derivative of yield strength with respect to equivalent plastic strain is calculated, where the result is not a constant E p , but a function of t p eqv : With these changes, all the deduction presented in this paper can still be verified. Therefore, the two properties, two load steps reduction rules and the empirical rule still apply.

Conclusion
In this paper, load step reduction method for adjoint sensitivity analysis of finite strain elastoplasticity is investigated. Two properties regarding adjoint variables are proved theoretically. Based on these properties, it proves that intermediate elastic load steps can be skipped in the sensitivity analysis without loss of accuracy; under certain conditions, the former ones of consecutive plastic load steps can also be skipped.
19 Page 18 of 23 phase. The nonlinear finite element analysis procedure is not influenced by it. Therefore, the partial derivatives of R with respect to design variables are the same for both step sets, i.e.: where the subscript S and L on the left denote the step set, from which R and H are calculated. Furthermore, partial derivatives of R and H are also the same for both step sets if they are independent on quantities at step t. Specifically, it means Obviously, adjoint variables from step N to step t + 1 are the same for both step sets since there is nothing changed yet in the backward solution procedure, i.e., Now we begin the proof.
Proof: Mathematically, it is to prove Eq. (50). Taking out equal terms in Eqs. (51) and (52), then it suffices to prove the following three items: The first one is obtained from Appendix C since step t is elastic.
For the second one, we firstly prove it holds at load step t − 1. Combining Eqs. (46) n s = n L and n s = n L , for for n ≤ t − 1.

Appendix E. Property of adjoint variables at a plastic step
In this appendix, the property of adjoint variables at a plastic step is mathematically proved.

Property 2. If
Since t+1 Δ p and t X p commute, then To c o m p u t e ln t+1 U e * ∕ t X p−1 , we c o m p u t e t+1 U e * ∕ t X p−1 first. Taking derivative of Eq. (9), using Eq. (8) and t+1 R e ∕ t X p−1 = 0 , it yields Premultiplying t+1 R e T on both sides results in For material which follows associated flow rule, from t+1 Δ p and t X p commute follows t+1 U e * and t X p also commute. Therefore using Eq. (E.6) it follows     Funding Open Access funding enabled and organized by Projekt DEAL.

Data availability Not applicable.
Code availability Not applicable.

Conflict of interest
The authors declare that they have no conflict of interest.  s H s ary conditions, mechanical loads, material properties and design nodes are described in full detail in Sect. 5.1. For the second to the fourth example, the finite element models in ABAQUS format are attached to this submission. At the same time, the boundary conditions, mechanical loads, material properties and design nodes are described respectively in Sect. 5.2, Sect. 5.3 and Sect. 6 in detail. To reproduce the results, the finite element analysis procedure is presented in Sect. 2 and the adjoint sensitivity analysis procedure is derived in Sect. 3. The authors wish to withhold the source code for commercialization purposes. This includes Python code implementing the finite strain elastoplastic analysis and adjoint sensitivity analysis.

Replication of results
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.