Efficient adjoint sensitivity analysis of isotropic hardening elastoplasticity via load steps reduction approximation

Accepted Manuscript Efficient adjoint sensitivity analysis of isotropic hardening elastoplasticity via load steps reduction approximation Wenjia Wang, Peter M. Clausen, Kai-Uwe Bletzinger

PII: DOI: Reference:

S0045-7825(16)31119-7 http://dx.doi.org/10.1016/j.cma.2017.07.020 CMA 11522

To appear in:

Comput. Methods Appl. Mech. Engrg.

Received date : 11 September 2016 Revised date : 14 July 2017 Accepted date : 18 July 2017 Please cite this article as: W. Wang, P.M. Clausen, K. Bletzinger, Efficient adjoint sensitivity analysis of isotropic hardening elastoplasticity via load steps reduction approximation, Comput. Methods Appl. Mech. Engrg. (2017), http://dx.doi.org/10.1016/j.cma.2017.07.020 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Efficient adjoint sensitivity analysis of isotropic hardening elastoplasticity via load steps reduction approximation Wenjia Wanga,*, Peter M. Clausenb, Kai-Uwe Bletzingera a

Lehrstuhl für Statik, TU München, Arcisstrasse 21 80333, München, Germany, [email protected], [email protected] b Dassault Systèmes Deutschland GmbH, Albert-Nestler-Str.17, 76131 Karlsruhe, Germany, [email protected]

Abstract: Sensitivity analysis plays a key role in gradient-based shape optimization. In this paper, a highly efficient procedure for design sensitivity analysis of elastoplastic material with isotropic hardening rule is presented. Here, an adjoint variable method is employed, where the adjoint sensitivity formulation is derived from the discrete master equation, yield condition, and consistency condition. Due to history dependence of plasticity, adjoint variables must be solved via a stepwise backward procedure. Therefore, the computational cost increases in proportion to the number of load steps. To address this cost and improve efficiency, condensation rules are proposed to reduce the number of load steps in the sensitivity analysis based on properties of adjoint variables. The properties of adjoint variables and proposed condensation rules are all theoretically proven in this paper. The accuracy and efficiency of the proposed techniques are demonstrated using bar truss structures and solid structures in various complex load cases. The results show that the techniques significantly reduce the number of load steps required in the sensitivity analysis while keeping satisfying accuracy. Keywords: sensitivity analysis, adjoint variable method, elastoplasticity, isotropic hardening, load steps reduction, gradient-based shape optimization

1. Introduction Methods employed to solve shape optimization problems fall into two categories, i.e., gradient-based and non-gradient-based methods. Starting from an initial point, gradient-based algorithms iteratively search for an optimal solution according to sensitivity information, i.e. derivatives of responses with respect to design variables. While design parameters naturally translate to design variables in the parametric shape optimization approach, non-parametric shape optimization instead typically employs nodal coordinates in a finite element model as its design variables. This provides a detailed description of structural shape, thus significantly enlarging the design space much more so than in parametric optimization. Sensitivity analysis for linear structural systems has been widely investigated, with numerous relevant successes in industrial applications. However, in practice, nonlinear structural behaviors occur equally if not more often than linear problems. The three types of nonlinearities are: geometric, material, and boundary condition nonlinearities. These nonlinearities may be encountered simultaneously or individually in large-scale engineering problems. Among the challenges in nonlinear sensitivity analysis, efficiency cannot be over emphasized. Firstly, given the characteristics of nonlinear problems, each sensitivity analysis maybe more complicated than corresponding linear cases, thus requiring more computation time and storage space. Secondly, the total number of sensitivity analysis in a nonlinear optimization procedure is greater than that in a linear optimization due to the increase in the number of iterations and possible need for 1

additional starting points. In the meantime, the accuracy of sensitivity is also important because it influences not only the search direction toward optima, but also the number of iterations required in the optimization procedure. Given the path dependence of general nonlinear material models that describe plasticity or viscosity, mechanical quantities of the current step are characterized not only by current loading conditions, but also by previous loading histories, as are the sensitivities [1, 2]. Therefore, sensitivities are often computed following each load step or time step in an incremental procedure [3, 4]. At each incremental step, a linearized design sensitivity equation must be solved [5]. By properly decomposing the total stress quantities, the derivative of the total stress rather than the incremental stress can also be obtained directly [6]. It has been pointed out that sensitivity analysis for path-dependent problems must be based on the consistent tangent operator, which specifies the variation of stresses with respect to changes in strains in the underlying finite element analysis [7, 8]. In material nonlinear sensitivity analysis, the non-differentiability problem occurs at transition points, such as points that transit from elasticity to plasticity and vice versa [9]. However, it is unlikely to hit such a critical point exactly, and taking the derivatives of adjacent points around the critical point has been shown to not lead to substantial errors [10]. Therefore, it is unnecessary to compute the sensitivities at the initial yield point to obtain sensitivities at subsequent steps. The direct differentiation method is a natural choice for forward marching problems and has been applied to various material nonlinear sensitivity analysis in an incremental process. Successfully, sensitivity computations have been identified for Norton–Soderberg power-law creep materials [6], spread plasticity with geometric nonlinearity under seismic loading [11], plane stress elasto-viscoplasticity [12], anisotropic elastoplastic shell [13], J2 plasticity with multi-yield-surface or nonlinear hardening [14, 15], Prandtl–Reuss elastoplasticity with geometric nonlinearity [10, 16], beam elements with plastic hinges [17], laminate composites with bilinear elastoplastic materials [18], nonlinear elasticity with frictionless contact [19], and elastoplasticity with frictional contact [20, 21]. When the number of design variables is larger than the number of responses, which is always the case for the non-parametric shape optimization, the adjoint variable method is preferred over the direct method in the sensitivity analysis. However, the usefulness of the adjoint variable method for step-by-step incremental problems was once questionable [22, 23]. The argument is based on the fact that each adjoint variable corresponds to only one performance function at a single step rather than the entire history. When sensitivities with respect to state variables at intermediate steps are required, these intermediate step quantities must be treated as additional responses to the system, which would enormously increase the total number of responses, thus weakening the advantages gained by using the adjoint variable method [8, 10, 24-26]. This bias against adjoint variable method, however, has been corrected. The key to solve the problem is to introduce additional adjoint variables to avoid derivatives evaluation with respect to intermediate state variables. Using this approach, shape sensitivity analysis for non-steady forming process is carried out by employing adjoint variables corresponding to the nodal balance equations and additionally to boundary conditions of material flow velocity [27]. By introducing adjoint variables that corresponding to boundary conditions on average residence time field in addition to those corresponding to equilibrium equations, it is able to perform sensitivity analysis 2

and optimization on polymer sheet extrusions and molding processes [28]. Based on the equilibrium-governing equations and a computational fluid grid induced by the elastic deformation of the structure and the fluid system, adjoint sensitivity information for coupled fluid-structure interaction problems is obtained [29]. By introducing adjoint variables for the equilibrium and constitutive relations, sensitivity for transient nonlinear coupled systems are obtained. Further, the sensitivity analysis with elastoplasticity is achieved as a direct extension [30]. It is shown that both direct and adjoint methods maintain their relative advantages as in the linear case for path-dependent problems [31]. Regardless of whether a direct or adjoint method is employed in the path-dependent sensitivity analysis, computational cost increases proportionally to the number of load steps. By empirically testing on single tetrahedral element and a two-bar truss structure comprised of bilinear elastoplastic material, sensitivities may only need to be calculated at the last step for monotonic loading procedure. For cyclical loading histories, design sensitivities must be computed at the unloading starting point and actual load point [31, 32]. However, there is still a need for systematic and theoretical investigations on reducing the number load steps properly, which is especially important in application to large-scale problems. In this paper, an efficient adjoint variable method for nonlinear sensitivity analysis with simultaneous geometric nonlinearity and elastoplasticity is investigated, where the hardening rule is assumed to be isotropic. In addition to this introduction, the paper is organized as follows. In Section 2, the iterative nonlinear finite element solution procedure and return mapping algorithm for elastoplasticity is introduced. Next, adjoint sensitivity formulations and solution procedures for adjoint variables are presented in Section 3. In Section 4, properties of adjoint variables and proposed rules to condense load steps are theoretically proven. Thereafter, the adjoint sensitivity analysis procedure and reduction rules are demonstrated via multiple examples in Section 5. Finally, in Section 6, the conclusions are presented. 2. Elastoplastic and geometric nonlinear finite element analysis The solution of a geometric and elastoplastic material nonlinear problem is the foundation based on which a sensitivity analysis procedure could be established. Therefore, in this section, the classical iterative nonlinear finite element solution procedure and return mapping algorithm for elastoplasticity is introduced. Path-dependence in two typical situations, i.e. two consecutive elastic load steps and consecutive plastic load steps with constant flow vector, are discussed from both physical point of view and finite element analysis point of view. It should be noted that, in this section and throughout the paper, matrix notations are employed in all formulas and deductions. 2.1 Elastoplastic material model and iterative solution procedure Plasticity is a material property in which a material undergoes irreversible deformation in response to applied loads. Elastoplastic material models depict materials that possess both elastic and plastic behavior. Most ductile metals, which are often used in engineering, fall into this category. The typical one-dimensional elastoplastic material behavior in a loading-unloading process is illustrated in Figure 1.

3

𝜎

E tan 𝜎0

𝐸

𝐸

𝜀𝑒

𝜀𝑃

𝜀

Figure1. Elastoplastic material behavior in a loading-unloading process A given elastoplastic material behaves in an elastic manner initially. Once the increasing stress exceeds initial yield strength 𝜎0 , plastic strain 𝜺𝑝 accumulates gradually, which is the irreversible strain after the material is fully unloaded. Total strain 𝜺 is additively decomposed into elastic strain 𝜺𝑒 and plastic strain 𝜺𝑝 as 𝜺 = 𝜺𝑒 + 𝜺𝑝

(1)

Plastic strain is often accompanied by large deformation that may lead to geometric nonlinearities. Therefore, to better fit general cases, elastoplasticity and geometric nonlinearity are simultaneously considered in this study by adopting Green–Lagrange strain and second Piola– Kirchhoff stress. Components of nonlinear Green–Lagrange strain εij in Cartesian coordinates are defined by 1 1 𝜀ij = (𝑢i,j + 𝑢j,i ) + (𝑢k,i ∙ 𝑢k,j ) 2 2

(2)

which can be reformulated in matrix form [33] as 𝟏 𝜺 = 𝑩 L 𝑼 + 𝑼T 𝑩 N 𝑼 𝟐

(3)

where 𝑩L represents the strain-displacement matrix for the linear part, 𝑩N is a symmetric matrix and 𝜺 is defined in Voigt notation. The relation between total stress and total strain is 𝝈 = 𝑫𝑒 𝜺𝑒 = 𝑫𝑒 (𝜺 − 𝜺𝑝 )

(4)

where 𝑫𝑒 represents the elastic constitutive relation and is equal to Young’s modulus in a onedimensional case. In the plastic stage, the stress satisfies the yield condition 𝑦(𝝈): 𝑦(𝝈) = |𝝈| − 𝜎Y (𝜺𝑝 ) = 0

(5)

where |∙| represents an equivalent stress measure, while 𝜎Y (𝜺𝑝 ) represents yield strength that depends on the material properties and accumulated plastic strains.

4

Three types of plasticity are classified according to the development of yield strength, i.e., hardening, perfect plasticity, and softening. Hardening describes the increase in yield strength as plastic strain accumulates, whereas softening describes the opposite. Perfect plasticity is an ideal model in which the yield strength remains the same. In Figure 1, 𝐸 tan denotes the tangent modulus which is related to Young’s modulus 𝐸 and plastic modulus 𝐸 𝑝 [34], i.e., 𝐸 tan =

𝐸𝐸 𝑝 𝐸 + 𝐸𝑝

(6)

with the tangent modulus and plastic modulus ranging from positive to negative, thereby reflecting different types of plasticity >0 hardening 𝐸 tan or 𝐸 𝑝 { = 0 perfect plasticity <0 softening

(7)

The total plastic strain is an accumulation of incremental plastic strains at each load step, i.e., 𝑝

𝜺𝑝 = ∑ ∆𝜺t t

(8)

where the subscript t denotes the load step and the incremental plastic strain is the product of an 𝑝 incremental equivalent plastic strain denoted by ∆𝜀equ and a flow vector. The flow vector is a normalized vector denoted by 𝐚 [35]. It thus has 𝑝

𝑝

∆𝜺t = ∆𝜀equ,t ∙ 𝐚t

(9)

In metal plasticity, it is typically assumed that incremental plastic strain follows an associated-flow rule, which means that the flow vector is in the same direction as the normal of the yield surface, i.e., 𝐚=

𝑑𝑦 𝑑|𝝈| = 𝑑𝝈 𝑑𝝈

(10)

If the direction of the flow vector is governed by a plastic potential other than the yield function, it then follows a non-associated flow rule. The Newton–Raphson method is often employed in nonlinear finite element analysis. Figure 2 depicts the flowchart of this method for elastoplasticity at new load step t+1 from current load step t [34].

5

𝑝

𝑝

Given: 𝑼t , 𝝈t , 𝜺t , 𝜺equ,t , 𝑭t , 𝑭t+1

𝑝

Assume: ∆𝑼 = 𝟎, ∆𝝈 = 𝟎, ∆𝜺𝑝 = 𝟎, ∆𝜀equ = 0 Update: 𝑼t+1 = 𝑼t + ∆𝑼 𝝈t+1 = 𝝈t + ∆𝝈 𝑝 𝑝 𝑝 𝑝 𝑝 𝑝 𝜺t+1 = 𝜺t + ∆𝜺 𝜀equ,t+1 = 𝜀equ,t + ∆𝜀equ

𝑝

𝑝

𝑲(𝑼t+1 , 𝝈t+1 , 𝜺t+1 , 𝜀equ,t+1 ) ∙ ∆𝑼 = 𝑭t+1 − 𝑭int  ∆𝑼

𝑝

𝑝

𝑝

∆𝑼  ∆𝝈, ∆𝜺𝑝 , ∆𝜀𝑒𝑞𝑢  update: 𝑼t+1 , 𝝈t+1 , 𝜺t+1 , 𝜀equ,t+1

𝑭int

𝑹 = 𝑭t+1 − 𝑭int

No

𝑹 <ε

Yes 𝑝

𝑝

return 𝑼t+1 , 𝝈t+1 , 𝜺t+1 , 𝜀equ,t+1 Figure 2. Newton–Raphson iterative solution procedure for elastoplasticity Since the constitutive relation between stress and strain for an elastoplastic material is not 𝑝 constant, stress 𝝈, plastic strain 𝜺𝑝 , and equivalent plastic strain 𝜀equ are also required in addition to displacement so as to obtain the tangent stiffness matrix at each step. These quantities cannot normally be determined explicitly and are usually obtained via a return mapping algorithm in the step highlighted by dashed frame. The return mapping algorithm is introduced in the next subsection together with the formulation of the consistent tangent stiffness matrix. 2.2 Return mapping algorithm and the consistent tangent stiffness matrix Given the quantities of the current step and a specified strain increment, the aim of the return mapping algorithm is to determine stress, plastic strain, and equivalent plastic strain of the new 𝑝 load step at each quadrature point. More specifically, given stress 𝝈t , plastic strain 𝜺t , equivalent 𝑝 plastic strain 𝜀equ,t , and incremental strain ∆𝜺 which is obtained from incremental displacement 𝑝

∆𝑼, trial stress is calculated under the assumption that ∆𝜀equ = 0 via ∆𝝈𝑡𝑟𝑦 = 𝑫𝑒 ∙ ∆𝜺 6

(11)

𝝈𝑡𝑟𝑦 = 𝝈t + ∆𝝈𝑡𝑟𝑦 Next, the equivalent stress of trial stress is compared to the yield strength at step t. Here, yield strength for an isotropic hardening material is calculated as 𝑝

𝑝 𝜎Y,t (𝜀equ,t )

𝜀equ,t

= 𝜎0 + ∫

0

𝑝

𝐸𝑝 ∙ 𝑑𝜀equ

(12)

If the equivalent stress is less than or equal to the yield stress, then the new load step is an elastic step. Therefore, the trial stress is the real stress and the incremental plastic strain is zero. Otherwise, the step contains plastic deformation and the iterative return mapping algorithm defined below is employed [34]. 𝑝

Set ∆𝜀equ = 0, ∆𝝈try = 𝑫𝑒 ∙ ∆𝜺, and 𝝈try = 𝝈t + ∆𝝈try. 1. Calculate trial flow vector 𝐚. 2. Update the increment of plastic strains as follows: ∆𝜆 =

|𝝈try | − 𝜎Y,t − ∫∆𝜀𝑝

equ,t

𝐸𝑝 ∙ 𝑑𝑠

𝐚 ∙ 𝑫𝑒 ∙ 𝐚T + 𝐸 𝑝 𝑝 ∆𝜀equ

𝑝

= ∆𝜀equ + ∆𝜆

∆𝜺𝑝 = ∆𝜺𝑝 + ∆𝜆 ∙ 𝐚 3. Update trial stress via 𝝈try = 𝝈try −𝑫𝑒 ∙ ∆𝜆 ∙ 𝐚. 4. Check whether |𝝈try | − 𝜎Y,t − ∫∆𝜀𝑝

equ,t

𝐸𝑝 ∙ 𝑑𝑠 = 0. If not, go to Step 1; otherwise, stops.

When the stopping criteria of the return mapping algorithm is satisfied, it gives 𝝈t+1 = 𝝈try 𝑝

𝑝

𝜺t+1 = 𝜺t + ∆𝜺𝑝 𝑝

𝑝

(13) 𝑝

𝜀equ,t+1 = 𝜀equ,t + ∆𝜀equ Given the above quantities, the consistent tangent stiffness matrix and internal force are expressed, based on [33], as 𝑲 = ∫ (𝑩L + 𝑼T 𝑩N )T 𝑫(𝑩L + 𝑼T 𝑩N)𝑑𝑉 + ∫ [𝛁U (𝑼T 𝑩N )]T 𝝈𝑑𝑉 𝑉

(14)

𝑉

N T

𝑭int = ∫ (𝑩L + 𝑼T 𝑩 ) 𝝈𝑑𝑉

(15)

𝑉

Note that the constitutive matrix 𝑫 is equal to 𝑫𝑒 in the elastic state. In the plastic state, the constitutive matrix is defined, based on [34, 36], by 𝑫 = 𝑫𝑒𝑝 = 𝑸−1 ∙ 𝑫𝑒 − where 7

𝒅 𝐚 ∙ 𝒓 + 𝐸𝑝

(16)

𝑸 = 𝑰 + 𝑫𝑒 ∙

𝑑𝐚 𝑝 ∆𝜀 𝑑𝝈 equ

(17)

𝒓 = 𝑸−1 ∙ 𝑫𝑒 ∙ 𝐚T 𝒅 = 𝒓 ∙ 𝒓T 2.3 Typical mechanical and computational behavior for elastoplasticity

In the finite element analysis of an elastoplastic and geometric nonlinear system, the loading procedure is often divided into a sequence of load steps. It is not only for the purpose of describing a loading history, but also for solvers to achieve better convergence. In the following, two typical situations in a sequence of load steps are discussed. They reveal that, if the difficulty in convergence can be disregarded, then some of load steps could be skipped without influencing the finite element analysis results. These properties enlighten further discussions on the load steps reduction in the sensitivity analysis. 2.3.1 Two consecutive elastic load steps A load step t is called an elastic load step if all elements behave elastically in that step. Assume step t and step t+1 are two consecutive elastic steps, then the following shows that the equilibrium state at step t+1 is independent of load step t. 𝑝

𝑝

Given quantities at step t-1, i.e., 𝑼t−1 , 𝝈t−1 , 𝜺t−1 and 𝜺equ,t−1, since both step t and step t+1 are elastic, it has 𝝈t+1 = 𝝈t−1 + 𝑫𝑒 (𝜺t+1 − 𝜺t−1 ) 𝑝

𝑝

𝜺t+1 = 𝜺t−1 𝑝

(18)

𝑝

𝜺equ,t+1 = 𝜺equ,t−1 Because the total strain is determined only by displacement via Eq. (3), so the stress 𝝈t+1 is a function of only displacement 𝑼t+1 . Thus the equilibrium displacement 𝑼t+1 could be solved with equilibrium equation N T

𝑭t+1 = 𝑭int = ∫ (𝑩L + 𝑼Tt+1 𝑩 ) ∙ 𝝈t+1 (𝑼t+1 )𝑑𝑉

(19)

𝑉

Therefore, the equilibrium states at time t+1 is determined solely by quantities at time t-1 and external load at t+1. They are independent of any quantity at step t. From physical point of view, it reveals that, if two consecutive elastic steps t and t+1 are encountered, then the former step will not influence the equilibrium quantities at the latter step. Thus the former elastic step t could be skipped in the solution procedure if it will not lead to convergence issues. What’s more, if there is a small perturbation in the structural system, the change of equilibrium quantities at the latter step will also be independent of the former step. Therefore, for a response that is a function of quantities at the latter step t+1, its sensitivities will not be influenced by the former step t, i.e., the former step could be skipped in the sensitivity analysis. 2.3.2 Consecutive plastic load steps with constant flow vectors 8

A load step t is called a plastic load step if any element behaves plastically in that load step. Assume step t-1, step t and step t+1 are consecutive plastic steps, and have constant flow vectors at every Gaussian point, then the followings show that the equilibrium states at step t+1 is independent of step t. 𝑝

𝑝

Given 𝑼t−1 , 𝝈t−1 , 𝜺t−1 , 𝜺equ,t−1 , and denote constant flow vector 𝐚, then according to the return mapping algorithm, it has 𝑝

𝑝

𝝈t = 𝝈t−1 + 𝑫𝑒 (𝜺t − 𝜺t−1 ) − 𝑫𝑒 (𝜺t − 𝜺t−1 ) 𝑝

𝑝

𝑝

𝜺t − 𝜺t−1 = ∆𝜀equ,t ∙ 𝐚 𝑝

∆𝜀equ,t

(20)

|𝝈t | − |𝝈t−1 | = 𝐚 ∙ 𝑫𝑒 ∙ 𝐚T + 𝐸 𝑝

and 𝑝

𝑝

𝝈t+1 = 𝝈t + 𝑫𝑒 (𝜺t+1 − 𝜺t ) − 𝑫𝑒 (𝜺t+1 − 𝜺t ) 𝑝

𝑝

𝑝

𝒑

|𝝈t+1 | − |𝝈t | = 𝐚 ∙ 𝑫𝑒 ∙ 𝐚T + 𝐸 𝑝

𝜺t+1 − 𝜺t = ∆𝜀equ,t+1 ∙ 𝐚 ∆𝜀equ,t+1

(21)

Substitute Eq. (20) into Eq. (21), it yields 𝑝

𝑝

𝝈t+1 = 𝝈t−1 + 𝑫𝑒 (𝜺t+1 − 𝜺t−1 ) − 𝑫𝑒 (𝜺t+1 − 𝜺t−1 ) 𝑝

𝑝

𝑝

𝜺t+1 − 𝜺t−1 = ∆𝜀equ ∙ 𝐚 𝑝

∆𝜀equ =

(22)

|𝝈t+1 | − |𝝈t−1 | 𝐚 ∙ 𝑫𝑒 ∙ 𝐚T + 𝐸 𝑝

Because total strain at step t+1 is determined by displacement 𝑼t+1 , via Eq. (3), so the stress 𝝈t+1 is also a function of only displacement 𝑼t+1 . Thus the equilibrium displacement 𝑼t+1 could be solved implicitly with Eqs. (19) and (22). Therefore, similar to the discussion in subsection 2.3.1, the equilibrium states at time t+1 is independent of step t, and load step t could be skipped in both finite element analysis and sensitivity analysis. 3. Sensitivity analysis with adjoint variable method In this section, the adjoint sensitivity formulation for elastoplasticity is introduced based on a selection of independent state variables. The backward solution procedure of adjoint variables is discussed thereafter. 3.1 Independent state variables Physical quantities such as stresses and strains of an elastic structure system are determined solely by the equilibrium displacement. Plasticity, however, introduces quantities that cannot be determined by displacements, such as stresses, plastic strains and equivalent plastic strains. According to Eq. (10), associated flow vectors are determined by stresses. As shown in Eqs. (14) to (16), the constitutive matrix, tangent stiffness matrix, and internal force vector of an 9

elastoplastic system are fully determined by equilibrium displacements, stresses, equivalent plastic strains and flow vectors. Further, given stresses and strains, plastic strains can be obtained through relation 𝑝

𝑫𝑒 (𝜺t − 𝜺t ) = 𝝈t

(23)

Therefore, all physical quantities at a load step can be explicitly expressed as functions of three variables, i.e., equilibrium displacements, stresses and equivalent plastic strains. As described in the return mapping algorithm, the incremental of these three quantities at a new load step is mutually determined. However, the accumulated total quantities of these variables are indeed independent. Thereby, a state vector 𝑽t that is comprised of stress and equivalent plastic strain is defined at each Gaussian point as 𝝈t 𝑽t = (𝜀 𝑝 ) equ,t

(24)

As discussed above, 𝑼t and 𝑽t at all load steps are two sets of independent state variables that are able to determine all mechanical quantities at any load step. Any system response can be either explicitly or implicitly expressed as a function of these variables and design variables s, i.e. 𝑓 = 𝑓(𝑼1 (𝑠), 𝑼2 (𝑠), … , 𝑼N (𝑠), 𝑽1 (𝑠), 𝑽2 (𝑠), … , 𝑽N (𝑠), 𝑠)

(25)

where N is the total number of load steps. It is worth mentioning that the selection of independent variables is not unique. For example, it is also possible to take plastic strain as an independent variable to replace the stress. The stress can then be obtained through Eq. (23). The difference in the selection of independent variables will influence the expression of equations and deductions in following sections. However, without verification in this paper, it should not influence any result and conclusion presented hereafter. 3.2 Deduction of adjoint sensitivity formulation According to Eq. (25), the sensitivity of a response 𝑓 with respect to a design variable is obtained directly following the chain rule N

N

𝑑𝑓 𝜕𝑓 𝜕𝑓 T 𝑑𝑼t 𝜕𝑓 T 𝑑𝑽t = +∑ +∑ 𝑑𝑠 𝜕𝑠 𝜕𝑼t 𝑑𝑠 𝜕𝑽t 𝑑𝑠 t=1

(26)

t=1

where direct evaluations of derivatives of 𝑼t and 𝑽t with respect to the design variable are required. If a large number of design variables are presented, which is typical in non-parametric shape optimization problems, then these evaluations are time-consuming, and the adjoint variable method is alternatively pursued. The idea of the adjoint variable method is to avoid the derivatives evaluation of 𝑼t and 𝑽t by canceling them out with additional terms that are products of adjoint variables and total derivatives of functions that are identical to zero. The adjoint variables are dependent on responses instead of design variables, which makes it preferred when number of design variables is larger than the number of responses. 10

The following deduction of adjoint sensitivity formulation is analog to that presented in [30], the difference is that here the definition of dependent residual function 𝑯t is simplified. In plasticity, the incremental stress and equivalent plastic strain follow yield and consistency conditions, thus 𝑝

𝑝

𝝈t = 𝝈t−1 + 𝑫𝒆 [(𝜺t − 𝜺t−1 ) − (𝜺t − 𝜺t−1 )] 𝑝 (𝜀equ,t

−

𝑝 𝜀equ,t−1 ) ∙

𝑝 [|𝝈| − 𝜎Y (𝜀equ,t )]

(27)

=0

(28)

∆𝜺t = 𝜺t − 𝜺t−1 = (𝜀equ,t − 𝜀equ,t−1 ) ∙ 𝐚t

(29)

where the increment of plastic strain is 𝑝

𝑝

𝑝

𝑝

𝑝

Eqs. (27) and (28) define a vector function 𝑯t as follows, which is identical to zero at any load step and is named dependent residual in literature [30] 𝑝

𝑯t (𝑼t , 𝑼t−1 , 𝑽t , 𝑽t−1 , 𝑠) = (

𝑝

𝝈t−1 + 𝑫𝒆 [(𝜺t − 𝜺t−1 ) − (𝜺t − 𝜺t−1 )] − 𝝈t 𝑝

𝑝

𝑝

(𝜀equ,t − 𝜀equ,t−1 ) ∙ [|𝝈| − 𝜎Y (𝜀equ,t )]

)≡𝟎

(30)

Here, the derivative of 𝑯t with respect to a design variable s is equal to 𝑑𝑯t 𝜕𝑯t 𝜕𝑯t 𝑑𝑼t 𝜕𝑯t 𝑑𝑼t−1 𝜕𝑯t 𝑑𝑽t 𝜕𝑯t 𝑑𝑽t−1 = + + + + 𝑑𝑠 𝜕𝑠 𝜕𝑼t 𝑑𝑠 𝜕𝑼t−1 𝑑𝑠 𝜕𝑽t 𝑑𝑠 𝜕𝑽t−1 𝑑𝑠

(31)

In addition, according to the equilibrium condition, the residual force function is zero at each load step N T

𝟎 ≡ 𝑹t (𝑼t , 𝑽t , 𝑠) = 𝑭t − 𝑭int = 𝑭t − ∫ (𝑩L + 𝑼T 𝑩 ) 𝝈𝑑𝑉 (t = 1,2, … , N)

(32)

𝑉

Here, the derivative of 𝑹t with respect to a design variable is 𝑑𝑹t ∂𝑹t 𝜕𝑹t 𝑑𝑼t 𝜕𝑹t 𝑑𝑽t = + ∙ + ∙ 𝑑𝑠 ∂s 𝜕𝑼t 𝑑𝑠 𝜕𝑽t 𝑑𝑠

(33)

If pre-multiply adjoint variables 𝝀t and 𝜸t to 𝑹t and 𝑯t respectively, and then subtract them from the response function, the response will not change because these terms are zero, i.e., it has N

𝑓=𝑓−

N

∑ 𝝀Tt 𝑹t t=1

− ∑ 𝜸Tt 𝑯t

(34)

t=1

According to Eqs. (26), (31) and (33), the derivative of Eq. (34) with respect to a design variable is N

N

N

t=1

t=1

t=1

𝑑𝑓 𝜕𝑓 𝜕𝑓 T 𝑑𝑼t 𝜕𝑓 T 𝑑𝑽t 𝜕𝑹t 𝜕𝑹t 𝑑𝑼t 𝜕𝑹t 𝑑𝑽t = +∑ +∑ − ∑ 𝝀Tt ( + + ) 𝑑𝑠 𝜕𝑠 𝜕𝑼t 𝑑𝑠 𝜕𝑽t 𝑑𝑠 𝜕𝑠 𝜕𝑼t 𝑑𝑠 𝜕𝑽t 𝑑𝑠 N

𝜕𝑯t 𝜕𝑯t 𝑑𝑼t 𝜕𝑯t 𝑑𝑼t−1 𝜕𝑯t 𝑑𝑽t 𝜕𝑯t 𝑑𝑽t−1 − ∑ 𝜸Tt ( + + + + ) 𝜕𝑠 𝜕𝑼t 𝑑𝑠 𝜕𝑼t−1 𝑑𝑠 𝜕𝑽t 𝑑𝑠 𝜕𝑽t−1 𝑑𝑠 t=1

11

(35)

By merging similar terms and forcing the coefficients of 𝑑𝑼t ⁄𝑑𝑠 and 𝑑𝑽t ⁄𝑑𝑠 to be zero, it has 𝑑

N

N

t=1

t=1

𝑑𝑓 𝜕𝑓 𝜕𝑹t 𝜕𝑯t = − ∑ 𝝀Tt − ∑ 𝜸Tt 𝑑𝑠 𝜕𝑠 𝜕𝑠 𝜕𝑠 𝜕𝑓 𝑇 𝜕𝑹N 𝜕𝑯N 𝑑𝑼N +( − 𝝀TN − 𝜸TN ) 𝜕𝑼N 𝜕𝑼N 𝑑𝑠 ⏟𝜕𝑼N =0

𝜕𝑓 𝑇 𝜕𝑹N 𝜕𝑯N 𝑑𝑽N +( − 𝝀TN − 𝜸TN ) 𝜕𝑽N 𝜕𝑽N 𝑑𝑠 ⏟𝜕𝑽N =0 N

𝜕𝑹t−1 − ∑ (𝝀Tt−1 ⏟ 𝜕𝑼t−1 t=2

+

𝜸Tt−1

𝜕𝑯t−1 𝜕𝑯t 𝜕𝑓 T 𝑑𝑼t−1 T + 𝜸t − ) 𝜕𝑼t−1 𝜕𝑼t−1 𝜕𝑼t−1 𝑑𝑠

(36)

=0

N

𝜕𝑹t−1 𝜕𝑯t−1 𝜕𝑯t 𝜕𝑓 T 𝑑𝑽t−1 − ∑ (𝝀Tt−1 + 𝜸Tt−1 + 𝜸Tt − ) 𝜕𝑽t−1 𝜕𝑽t−1 𝜕𝑽t−1 𝑑𝑠 ⏟ 𝜕𝑽t−1 t=2

=0

𝜕𝑯1 𝑑𝑼0 𝜕𝑯1 𝑑𝑽0 −𝜸1T ( + ) 𝜕𝑼0 𝑑𝑠 𝜕𝑽0 𝑑𝑠 ⏟ ≡0

where a sequence of linear systems of equations for adjoint variables are obtained 𝜕𝑹N T 𝜕𝑼N 𝜕𝑹N T ( 𝜕𝑽N 𝜕𝑹t−1 T 𝜕𝑼t−1 𝜕𝑹t−1 T ( 𝜕𝑽t−1

𝜕𝑯N T 𝜕𝑓 𝜕𝑼N 𝜕𝑼N 𝝀N ∙( )= 𝜸N 𝜕𝑓 𝜕𝑯N T 𝜕𝑽N ) ( 𝜕𝑽N )

𝜕𝑯t−1 T 𝜕𝑯t 𝑇 𝜕𝑓 𝜕𝑼t−1 𝜕𝑼t−1 𝜕𝑼t−1 𝝀t−1 ∙( )=− ∙ 𝜸t + (t = N, … ,2) 𝜸t−1 𝜕𝑓 𝜕𝑯t 𝑇 𝜕𝑯t−1 T 𝜕𝑽t−1 ) ( 𝜕𝑽t−1 ) ( 𝜕𝑽t−1 )

(37)

(38)

Finally, the sensitivity of the response is equal to N

N

t=1

t=1

𝑑𝑓 𝜕𝑓 𝜕𝑹t 𝜕𝑯t = − ∑ 𝝀Tt − ∑ 𝜸Tt 𝑑𝑠 𝜕𝑠 𝜕𝑠 𝜕𝑠

(39)

Note that the linear systems in Eqs. (37) and (38) can only be solved by going backward, i.e., by traversing from the last step to the first step. 12

3.3 Solutions of the adjoint variables To evaluate sensitivities using Eq. (39), the adjoint variables must first be obtained according to Eqs. (37) and (38). Solving Eq.(37) directly, it yields 𝜕𝑹N T 𝜕𝑯N T 𝜕𝑯N −T 𝜕𝑹N T 𝜕𝑓 𝜕𝑯N T 𝜕𝑯N −T 𝜕𝑓 − − ( ) 𝝀N = 𝜕𝑼N 𝜕𝑼N 𝜕𝑽N 𝜕𝑽N 𝜕𝑼N 𝜕𝑼N 𝜕𝑽N 𝜕𝑽N

(40)

𝜕𝑯N −T 𝜕𝑓 𝜕𝑹N T 𝜸N = − 𝝀 ) ( 𝜕𝑽N 𝜕𝑽N 𝜕𝑽N N Next, it will show that the following equation holds for each load step: T

𝜕𝑹t T 𝜕𝑯t T 𝜕𝑯t −T 𝜕𝑹t T 𝜕𝑹t 𝜕𝑹t 𝜕𝑯t −1 𝜕𝑯t − =( − ) = −𝑲t 𝜕𝑼t 𝜕𝑼t 𝜕𝑽t 𝜕𝑽t 𝜕𝑼t 𝜕𝑽t 𝜕𝑽t 𝜕𝑼t

(41)

To prove this, first, derivatives of residual force and dependent residual are deduced as N T

𝜕 ∫𝑉 (𝑩L + 𝑼T 𝑩 ) 𝝈t 𝑑𝑉 𝜕𝑹t T =− = − ∫ 𝛁𝐔 (𝑼T 𝑩N ) 𝝈t 𝑑𝑉 𝜕𝑼t 𝜕𝑼t 𝑉 𝜕 𝜕𝑹t = −( 𝜕𝝈t 𝜕𝑽t

𝜕

)∫ 𝑝 𝜕𝜀equ,t 𝑉

N T

N T

(𝑩L + 𝑼T 𝑩 ) 𝝈t 𝑑𝑉 = − (∫ (𝑩L + 𝑼T 𝑩 ) 𝑑𝑉

(42)

0)

(43)

𝑉

The formulation for 𝑯t can be treated separately for an elastic and a plastic step. If load step t is an elastic step, then 𝑯t has the form 𝑝

𝑝

𝝈t−1 + 𝑫𝒆 [(𝜺t − 𝜺t−1 ) − (𝜀equ,t − 𝜀equ,t−1 ) ∙ 𝐚t ] − 𝝈t 𝑯t = ( )=0 𝑝 𝑝 𝜀equ,t − 𝜀equ,t−1

(44)

According to Eq (3), the derivatives of 𝑯t are 𝜕𝜺t L T N 𝑒 𝑒 𝜕𝑯t 𝑫 = ( 𝜕𝑼t ) = (𝑫 (𝑩 + 𝑼 𝑩 )) 𝜕𝑼t 0 0

(45)

𝜕𝜺t−1 𝑒 𝜕𝑯t 𝜕𝑯t−1 = (−𝑫 𝜕𝑼t−1 ) = − 𝜕𝑼t−1 𝜕𝑼t−1 0

(46)

𝜕𝜺t 𝜕𝜺t−1 𝜕(𝜺t − 𝜺t−1 ) 𝜕𝑯t 𝑒 𝑒 𝑒 𝑫 𝑫 𝑫 =( )=( 𝜕𝑠 ) − ( 𝜕𝑠 ) 𝜕𝑠 𝜕𝑠 0 0 0

(47)

𝜕 𝜕𝑯t = (𝜕𝝈 𝜕𝑽t t

𝜕

𝜕𝑯t −𝑰 −𝑫𝑒 𝐚t )= ) 𝑯t = ( 𝑝 𝜕𝜀equ,t 0 1 𝜕𝑽t

𝜕 𝜕𝑯t = (𝜕𝝈 𝜕𝑽t−1 t−1

𝜕

) 𝑯t 𝑝 𝜕𝜀equ,t−1 13

𝑰 =( 0

−1

𝜕𝑯t 𝑫𝑒 𝐚t )=− −1 𝜕𝑽t

(48) (49)

where 𝜕𝜺t 𝜕𝑩L 1 𝜕𝑩N = 𝑼t + 𝑼Tt 𝑼 𝜕𝑠 𝜕𝑠 2 𝜕𝑠 t

(50)

Instead, if load step t is a plastic step, then 𝑯t has the form 𝑝

𝑝

𝝈t−1 + 𝑫𝒆 [(𝜺t − 𝜺t−1 ) − (𝜀equ,t − 𝜀equ,t−1 ) ∙ 𝐚t ] − 𝝈t 𝑯t = ( )=𝟎 𝑝 |𝝈| − 𝜎Y (𝜀equ,t )

(51)

which is different from Eq. (44) in the second component. Therefore, derivatives of dependent residual with respect to displacements and design variables are identical to Eqs. (45) to (47), and derivatives with respect to state vectors are 𝜕 𝜕𝑯t = (𝜕𝝈 𝜕𝑽t t

𝜕

−𝑸 −𝑫𝑒 𝐚t =( T ) 𝐚t −𝐸 𝑝

) 𝑯t 𝑝 𝜕𝜀equ,t

𝜕 𝜕𝑯t = (𝜕𝝈 𝜕𝑽t−1 t−1

𝜕

) 𝑯t 𝑝 𝜕𝜀equ,t−1

(52)

𝑫𝑒 𝐚 t ) 0

𝑰 =( 0

(53)

where 𝑸 is defined in Eq. (17). Due to the different formulations for an elastic and a plastic step, Eq. (41) is also proven separately. First, if loading step t is an elastic step, substitute Eqs. (42), (43), (45), (48) into Eq. (41) and compare with Eq. (14), it has 𝜕𝑹t 𝜕𝑹t 𝜕𝑯t−1 𝜕𝑯t − 𝜕𝑼t 𝜕𝑽t 𝜕𝑽t 𝜕𝑼t T

N T

= − ∫ 𝛁𝐔 (𝑼 𝑩 ) 𝝈t 𝑑𝑉 𝑉

N T

+ (∫ (𝑩L + 𝑼T 𝑩 ) 𝑑𝑉

N

−𝑰 −𝑫𝑒 𝐚t ) (𝑫𝑒 (𝑩L + 𝑼T𝑩 )) 0 1 0

0) (

𝑉

N T

(54)

N

= − ∫ (𝑩L + 𝑼T 𝑩 ) 𝑫𝑒 (𝑩L + 𝑼T 𝑩 ) 𝑑𝑉 − ∫ [𝛁U (𝑼T 𝑩N )]T 𝝈t 𝑑𝑉 𝑉

𝑉

= −𝑲t

Given the symmetry of the tangent stiffness matrix, it then has T

𝜕𝑹t T 𝜕𝑯t T 𝜕𝑯t −T 𝜕𝑹t T 𝜕𝑹t 𝜕𝑹t 𝜕𝑯t −1 𝜕𝑯t − =( − ) = −𝑲t 𝜕𝑼t 𝜕𝑼t 𝜕𝑽t 𝜕𝑽t 𝜕𝑼t 𝜕𝑽t 𝜕𝑽t 𝜕𝑼t Second, for a plastic load step t, denote 14

(55)

𝜕𝑯t −1 𝑨 =( 𝑪 𝜕𝑽t

𝑩 ) 𝑫

(56)

According to Eq. (52), it can be verified that −1

𝑨=𝑸

𝒅 ∙ 𝑫𝑒 −1 − 𝐚 ∙ 𝒓 + 𝐸𝑝

(57)

Therefore, substitute Eqs. (42), (43), (45), (52) into Eq. (41) and compare with Eqs. (14) and (16), it has 𝜕𝑹t 𝜕𝑹t 𝜕𝑯t −1 𝜕𝑯t − 𝜕𝑼t 𝜕𝑽t 𝜕𝑽t 𝜕𝑼t T

N T

N T

N T

N

𝑩 𝑫𝑒 (𝑩L + 𝑼T 𝑩 ) )( ) 𝑫 0

𝑨 0) ( 𝑪

= − ∫𝑉 𝛁𝐔 (𝑼 𝑩 ) 𝝈t 𝑑𝑉 + (∫ (𝑩L + 𝑼T 𝑩 ) 𝑑𝑉 𝑉 N

T

N T

= − ∫ (𝑩L + 𝑼T 𝑩 ) 𝑑𝑉 ∙ 𝑨 ∙ 𝑫𝒆 ∙ (𝑩L + 𝑼T 𝑩 ) − ∫ 𝛁𝐔 (𝑼 𝑩 ) 𝝈t 𝑑𝑉

(58)

𝑉

𝑉

N T

N

T

N T

= − ∫ (𝑩L + 𝑼T 𝑩 ) 𝑫𝑒𝑝 (𝑩L + 𝑼T 𝑩 ) 𝑑𝑉 − ∫ 𝛁𝐔 (𝑼 𝑩 ) 𝝈t 𝑑𝑉 𝑉

𝑉

= −𝑲t

Using Eq. (41), adjoint variable 𝝀N and 𝜸N are obtained by solving the linear system of equations at the last load step, i.e., 𝑲 N ∙ 𝝀N = − (

𝜕𝑓 𝜕𝑯N T 𝜕𝑯N −T 𝜕𝑓 − ) 𝜕𝑼N 𝜕𝑼N 𝜕𝑽N 𝜕𝑽N

𝜕𝑯N −𝑇 𝜕𝑓 𝜕𝑹N 𝑇 𝜸N = − 𝝀 ) ( 𝜕𝑽N 𝜕𝑽N 𝜕𝑽N N According to Eq. (38), going backward 1 are obtained by solving the linear systems

(59)

, adjoint variables 𝝀t and 𝜸t for step N-1 through

𝜕𝑯t T 𝜕𝑯t −T 𝜕𝑯t+1 T 𝜕𝑯t+1 T 𝜕𝑓 𝑲 t ∙ 𝝀t = − ( − ) 𝜸t+1 + 𝜕𝑼t 𝜕𝑽t 𝜕𝑽t 𝜕𝑼t 𝜕𝑼t 𝜕𝑯t −T 𝜕𝑯t+1 T 𝜕𝑹t 𝑇 𝜕𝑓 𝜸t = − 𝜸t+1 + 𝝀 )+ ( 𝜕𝑽t 𝜕𝑽t 𝜕𝑽t t 𝜕𝑽t

(60)

4. Reducing computational cost by condensing the load steps From the previous section, it can be seen that the major computational cost in elastoplastic sensitivity analysis is the solution of adjoint variables. Since matrix 𝜕𝑯t ⁄𝜕𝑽t is block-wise diagonalized, the linear systems to obtain 𝜸t are automatically decoupled, and thus allowing us to 15

ignore the solution cost of 𝜸t . To obtain 𝝀t at all load steps, it requires solutions of N linear systems each with a size equal to the total degrees of freedom of the underlying finite element model. Therefore, the total computational cost is roughly proportional to the number of load steps. In this section, rules for reducing the number of load steps required in the sensitivity analysis are investigated. As discussed in section 2.3.1, the former one of two consecutive elastic load steps could be skipped in the sensitivity analysis. As discussed in section 2.3.2, intermediate steps in a sequence of plastic load steps with constant flow vectors could also be skipped in the sensitivity analysis. In this section, above two observations are extended to cases of a single elastic load step and two consecutive plastic load steps. Beforehand, two properties regarding the adjoint variables are proposed and proven. Based on these properties, two condensation rules, which show that a majority of load steps can actually be skipped during the sensitivity analysis, are then proposed and theoretically proven. For a sequence of load steps L={1, 2,..., t-2, t-1, t,…, N}, system responses that are investigated in this section are assumed to be functions only of the last step quantities. Mathematically speaking, this assumes that a system response 𝑓 satisfies 𝜕𝑓 =0 𝜕𝑼t 𝜕𝑓 =0 𝜕𝑽t

for all t < N

(61) for all t < N

On one hand, a number of practical mechanical responses fall into this category, including but not limited to residual displacements, maximum equivalent plastic strain, and maximum equivalent stress in the entire loading history. On the other hand, a general response g can often be expressed as a composition of a set of such functions, i.e., g(𝑼1 , 𝑼2 , … , 𝑼N , 𝑽1 , 𝑽2 , … , 𝑽N ) = ℎ°(𝑓1 (𝑼1 , 𝑽1 ), 𝑓2 (𝑼2 , 𝑽2 ), … , 𝑓N (𝑼N , 𝑽N ))

(62)

where 𝑓k (𝑼k , 𝑽k ) is dependent only on 𝑼k and 𝑽k. The sensitivity of g can therefore be obtained by the chain rule with sensitivities of 𝑓k. To calculate the sensitivity of 𝑓k, it only needs to treat step k as the last load step, thus 𝑓k is a function that satisfies the assumption in Eq. (61). 4.1 Two properties of adjoint variables In this subsection, a load step is called an elastic step if all elements in the finite element system behave elastically. Otherwise, the load step is called a plastic load step. Given this convention, we propose and prove two properties described below regarding adjoint variables. Property 1. If load step t−1 (t−1
(63)

According to Eqs. (49) and (53), 𝜕𝑯t 𝑰 𝑫𝑒 𝐚t =( ) 𝜕𝑽t−1 𝟎 𝑐 𝑐={

−1 0

if step t is elastic if step t is plastic

(64)

Eqs. (45), (46), (63) and (64) together lead to 𝜕𝑯t−1 T 𝜕𝑯t−1 −T 𝜕𝑯t T 𝜕𝑯t T − 𝜕𝑼t−1 𝜕𝑽t−1 𝜕𝑽t−1 𝜕𝑼t−1 =

𝜕𝑯t−1 T 𝜕𝑯t−1 −T 𝜕𝑯t T 𝜕𝑯t−1 T + 𝜕𝑼t−1 𝜕𝑽t−1 𝜕𝑽t−1 𝜕𝑼t−1

= ((𝑩L + 𝑼t−1 𝑩N )𝐓 𝑫𝒆 = ((𝑩L + 𝑼t−1 𝑩N )𝐓 𝑫𝒆

−𝑰 𝟎 𝑰 𝟎 )( ) + 𝑰) 0) (( T −𝐚t−1 𝑫𝑒 1 𝐚Tt 𝑫𝑒 𝑐 𝟎 𝟎 ) 0) ( T (𝐚t − 𝐚Tt−1 )𝑫𝑒 𝑐 + 1

(65)

=𝟎

According to Eq. (60) and the assumption described in Eq. (61), it has 𝑲t−1 ∙ 𝝀t−1

𝜕𝑯t−1 T 𝜕𝑯t−1 −T 𝜕𝑯t T 𝜕𝑯t T 𝜕𝑓 = −( − =𝟎 ) 𝜸t + 𝜕𝑼t−1 𝜕𝑽t−1 𝜕𝑽t−1 𝜕𝑼t−1 𝜕𝑼t−1

(66)

which leads to 𝝀t−1 = 𝟎

(67)

□ Property 2. If load steps t−1 and t are both plastic, and flow vectors satisfy 𝜸𝝈,t−1 𝜸 d𝐚 𝐚t−1 = 𝐚t , t−1 = 𝟎, then 𝝀t−1 = 0, 𝜸t−1 = (𝜸 𝑝 ) = ( 𝝈,t ). d𝝈t−1 𝜀 ,t−1 0 equ

Proof: Denote 𝜕𝑯t−1 −1 𝑨 𝑩 =( ) 𝑪 𝑫 𝜕𝑽t−1

(68)

Since step t is plastic, according to Eq. (53) and the assumption that 𝐚t−1 = 𝐚t , it has 𝜕𝑯t 𝑰 𝑫𝑒 𝐚 t 𝑰 𝑫𝑒 𝐚t−1 =( )=( ) 𝜕𝑽t−1 𝟎 0 𝟎 0 Thus, with Eqs. (45) and (46), it obtains 17

(69)

𝜕𝑯t−1T 𝜕𝑯t−1 −T 𝜕𝑯t T 𝜕𝑯t T 𝜕𝑯t−1 T 𝑨T 𝑪T 𝑰 𝟎 − = [( T )( T ) + 𝑰] 𝑒 T 𝐚 𝑫 0 𝜕𝑼t−1 𝜕𝑽t−1 𝜕𝑽t−1 𝜕𝑼t−1 𝜕𝑼t−1 t−1 𝑩 𝑫

(70)

Since load step t-1 is also plastic, according to Eq. (52), 𝜕𝑯t−1 −𝑸 =( T 𝐚t−1 𝜕𝑽t−1

−𝑫𝑒 𝐚t−1 ) −𝐸 𝑝

(71)

where 𝑸 = 𝑰 + 𝑫𝑒 ∙

𝑑𝐚t−1 𝑝 ∙ ∆𝜀equ,t−1 𝑑𝝈t−1

(72)

Note it assumes that 𝑑𝐚t−1 =𝟎 𝑑𝝈t−1

(73)

𝑸=𝑰

(74)

Therefore,

and further 𝑰=

T 𝜕𝑯t−1 −T 𝜕𝑯t−1 T = (𝑨T 𝜕𝑽t−1 𝜕𝑽t−1 𝑩

𝑪T ) ( −𝑰 𝑫T −𝐚Tt−1 𝑫𝑒

𝐚t−1 ) −𝐸 𝑝

(75)

Therefore T

(𝑨 T

𝑩

0 −𝑰 𝟎 𝑪T ) ( 𝑰 )=( ) T 𝑒 T 𝐚t−1 𝑫 0 𝟎 0 𝑫

(76)

Substitute Eqs. (45) and (76) into Eq. (70), it yields 𝜕𝑯𝑡−1 T 𝜕𝑯𝑡−1 −T 𝜕𝑯𝑡 T 𝜕𝑯𝑡 T T − = ((𝑩L + 𝑼t−1 𝑩N ) 𝑫𝒆 𝜕𝑼𝑡−1 𝜕𝑽𝑡−1 𝜕𝑽𝑡−1 𝜕𝑼𝑡−1

𝟎 𝟎 0) (𝟎 1) = 𝟎

(77)

Therefore, according to Eq. (60) and the assumption described in Eq. (61), it has 𝑲t−1 ∙ 𝝀t−1

𝜕𝑯𝑡−1 T 𝜕𝑯𝑡−1 −T 𝜕𝑯𝑡 T 𝜕𝑯𝑡 T 𝜕𝑓 = −( − =𝟎 ) 𝜸t + 𝜕𝑼𝑡−1 𝜕𝑽𝑡−1 𝜕𝑽𝑡−1 𝜕𝑼𝑡−1 𝜕𝑼t−1

(78)

which leads to 𝝀t−1 = 0 Further,

18

(79)

𝜕𝑯𝑡−1 −T 𝜕𝑯𝑡 T 𝜕𝑽𝑡−1

𝜕𝑽𝑡−1

T

= (𝑨T 𝑩

𝑪T ) ( 𝑰 T 𝑒 𝑫T 𝐚t−1 𝑫

0 𝑰 ) = −( 0 0

0 ) 0

(80)

Therefore, 𝜸t−1 = −

𝜕𝑯𝑡−1 −T 𝜕𝑽𝑡−1

(

𝜕𝑯𝑡 T 𝜕𝑽𝑡−1

𝜸t +

𝜕𝑹𝑡−1 T 𝜕𝑽𝑡−1

𝝀t−1 ) +

𝜕𝑓 𝑰 =( 0 𝜕𝑽t−1

𝜸 0 ) 𝜸 = ( 𝝈,t ) 0 t 0

(81)

□ 4.2 Computational cost reduction via condensing the load steps Based on above properties of adjoint variables, following two rules are proposed and theoretically proven in this subsection. These rules show that, under certain conditions, some load steps can be skipped during sensitivity analysis without influencing the results. Rule 1. Given a sequence of loading step L={1, 2,..., t−2, t−1, t,…, N}, if step t-1 (t-1
(82)

where the superscript S and L denote quantities obtained with load steps sequence S and L respectively. Lastly, since the finite element analysis procedure are not changed, it is also obvious that, quantities and partial derivatives which are determined by displacements, stresses, strains, equivalent plastic strains at steps other than t-1 are also the same for load steps sequence S and L. Specifically, it has 𝜕𝑹Sn 𝜕𝑹Ln (for n ≠ t − 1) = 𝜕𝑠 𝜕𝑠 𝜕𝑯Sn 𝜕𝑯Ln = (for n ≠ t − 1 and n ≠ t) 𝜕𝑠 𝜕𝑠 and 19

(83)

𝜕𝑯Sn 𝜕𝑯Ln = (for n ≠ t − 1) 𝜕𝑼n 𝜕𝑼n 𝜕𝑯Sn 𝜕𝑯Ln = (for n ≠ t − 1) 𝜕𝑽n 𝜕𝑽n 𝜕𝑯Sn 𝜕𝑯Ln = (for n ≠ t − 1 and n ≠ t) 𝜕𝑽n−1 𝜕𝑽n−1

(84)

𝜕𝑯Sn 𝜕𝑯Ln = (for n ≠ t − 1) 𝜕𝑼n−1 𝜕𝑼n−1 𝑲Sn = 𝑲Ln (for n ≠ t − 1) Next, we enter the proof of the first condensation rule. Proof: Mathematically, it needs to prove that 𝑑𝑓 L 𝑑𝑓 S = 𝑑𝑠 𝑑𝑠

(85)

According to Eq. (39), both sensitivity formulations are composed of three terms, i.e., 𝑑𝑓 L 𝜕𝑓 = − ΛL − Γ L 𝑑𝑠 𝜕𝑠

(86)

𝑑𝑓 S 𝜕𝑓 = − ΛS − Γ S 𝑑𝑠 𝜕𝑠 where t−2 L

Λ =

T ∑ 𝝀Ln n=1

N

L L 𝜕𝑹Ln T 𝜕𝑹t−1 T 𝜕𝑹n + 𝝀Lt−1 + ∑ 𝝀Ln 𝜕𝑠 𝜕𝑠 𝜕𝑠 n=t

t−2

N

n=1

n=t

(87)

S S T 𝜕𝑹n T 𝜕𝑹n ΛS = ∑ 𝝀Sn + ∑ 𝝀Sn 𝜕𝑠 𝜕𝑠

and t−2 L

Γ =

T ∑ 𝜸Ln n=1

N

L L L 𝜕𝑯Ln T 𝜕𝑯t−1 T 𝜕𝑯t T 𝜕𝑯n + 𝜸Lt−1 + 𝜸Lt + ∑ 𝜸Ln 𝜕𝑠 𝜕𝑠 𝜕𝑠 𝜕𝑠 n=t+1

t−2

ΓS =

S T 𝜕𝑯n ∑ 𝜸Sn 𝜕𝑠 n=1

+

S T 𝜕𝑯t 𝜸St

𝜕𝑠

N

+ ∑ n=t+1

S T 𝜕𝑯n 𝜸Sn

(88)

𝜕𝑠

By comparing both formulas in Eqs. (87) and (88), and taking Eqs. (83) and (84) into account, Eq. (85) holds if the following three items are proven: 20

1. 𝝀Lt−1 = 0 2. 𝝀Sn = 𝝀Ln and 𝜸Sn = 𝜸Ln (for n ≤ t − 2) 3. 𝜸Lt−1

T 𝜕𝑯L t−1 𝜕𝑠

+ 𝜸Lt

T 𝜕𝑯L t

= 𝜸St

𝜕𝑠

T 𝜕𝑯St 𝜕𝑠

It is rather straightforward to obtain the first item from Property 1 in section 4.1, since step t−1 is elastic. For the second item, we first prove that it holds at step t-2. According to Eqs. (60), (61) and (84), it has T

𝝀Lt−2 = 𝑲Lt−2

=

−1 𝑲St−2

−1

𝜕𝑯Lt−2 𝜕𝑯Lt−2 ∙ [( 𝜕𝑼t−2 𝜕𝑽t−2 T

S

S

𝜕 𝑯t−2 𝜕𝑯t−2 ∙( 𝜕 𝑼t−2 𝜕 𝑽t−2

−T

−T

𝜕𝑯Lt−1

T

−

𝜕𝑽t−2

𝜕𝑯Lt−1 𝜕𝑼t−2

T

L

L

T

) 𝜸Lt−1 ]

(89)

T

𝜕𝑯t−1 L 𝜕𝑯t−1 L 𝜸 − 𝜸 ) 𝜕 𝑽t−2 t−1 𝜕 𝑼t−2 t−1

According to item 1 and assumption in Eq. (61) 𝜸Lt−1

=−

𝜕𝑯Lt−1

(

𝜕𝑽t−1 L

=−

−T

𝜕 𝑯t−1 𝜕 𝑽t−1

𝜕𝑯Lt

𝜕𝑽t−1

−T

T

𝜸Lt

+

𝜕𝑹Lt−1

T

𝝀Lt−1 ) +

𝜕𝑽t−1

L T

𝜕 𝑯t 𝑰 𝜸L = ( 𝐚t−1 T 𝑫𝑒 𝜕𝑽t−1 t

𝑰

0

𝜕𝑓 𝜕𝑽𝑡−1

0 𝑰 )( −1 𝐚t T 𝑫𝑒

0 S )𝜸 𝑐 t

(90)

S

𝐚t−1 T − 𝐚t T )𝑫𝑒 −𝑐 ) 𝜸t

= ((

Therefore, L

T

𝜕𝑯t−1 L 𝑰 𝜸 =( 𝐚t−1 T 𝑫𝑒 𝜕 𝑽t−2 t−1

𝑰 0 𝑰 0 S 0 )( ) 𝜸S = ( T 𝒆 )𝜸 −1 (𝐚t−1 T − 𝐚t T )𝑫𝑒 −𝑐 t 𝐚t 𝑫 𝑐 t

S T

𝜕 𝑯t = 𝜸S 𝜕𝑽t−2 t

(91)

According to Eqs. (45), (46) and (84) −

𝜕𝑯Lt−1 𝜕𝑯Lt−2 𝜕𝑯St−2 𝑫𝑒 (𝑩L + 𝑼t−2 T 𝑩N )) = = =( 𝜕𝑼t−2 𝜕𝑼t−2 𝜕𝑼t−2 0

(92)

Therefore, L

−

T

𝜕𝑯t−1 L T 𝜸 = ((𝑩L + 𝑼t−2 T 𝑩N ) 𝑫𝑒 𝜕 𝑼t−2 t−1 L

T

N T

0) ((𝐚

= ((𝑩 + 𝑼t−2 𝑩 ) 𝑫

𝑒

21

t−1

0) 𝜸St

T

𝑰 − 𝐚t T )𝑫𝑒 T

𝜕𝑯St−2 S = 𝜸 𝜕𝑼t−2 t

0 ) 𝜸S −𝑐 t

(93)

As discussed above, step t-2 and step t are two adjacent steps in load sequence S. Therefore, according to Eq. (46), it has 𝜕𝑯St−2 𝜕𝑯St =− 𝜕𝑼t−2 𝜕𝑼t−2

(94)

Substitute Eqs. (91), (93) and (94) into Eq. (89), it yields T

𝝀Lt−2

=

−1 𝑲St−2

𝜕𝑯St−2 𝜕𝑯St−2

∙(

−T

𝜕𝑼t−2 𝜕𝑽t−2

𝜕𝑯St 𝜕𝑽t−2

T

T

𝜕𝑯St − ) 𝜸St = 𝝀St−2 𝜕𝑼t−2

(95)

Further, it follows 𝜸Lt−2 = −

=−

𝜕𝑯Lt−2

−T

𝜕𝑽t−2

𝜕𝑯St−2

−T

𝜕𝑽t−2

𝜕𝑯Lt−1

(

(

𝜕𝑽t−2

𝜕𝑯St

𝜕𝑽t−2

T

𝜸Lt−1 +

T

𝜸St

+

𝜕𝑹St−2 𝜕𝑽t−2

𝜕𝑹Lt−2

T

𝜕𝑽t−2

𝝀Lt−2 )

(96)

T

𝝀St−2 )

=

𝜸St−2

Since 𝝀S is equal to 𝝀L at t-2 and 𝜸S is also equal to 𝜸L at t-2, according to Eq. (84), the second item is obtained obviously 𝝀Sn = 𝝀Ln and 𝛄Sn = 𝛄Ln (for n ≤ t − 2)

(97)

For the third item, according to Eqs. (47) and (90), it has 𝜸Lt−1

L T 𝜕𝑯t−1

𝜕𝑠

T 𝑰 = 𝜸St ( 𝟎

+ 𝜸Lt

L T 𝜕𝑯t

𝜕𝑠

𝜕𝜺 𝜕𝜺 𝜕𝜺t 𝜕𝜺t−1 𝒆 𝑫𝑒 (𝐚t−1 − 𝐚t ) 𝑫𝒆 ( t−1 − t−2 ) − ) S T (𝑫 ( )( ) + 𝜸 𝜕𝑠 𝜕𝑠 𝜕𝑠 𝜕𝑠 ) t −𝑐 0 0

𝜕𝜺t 𝜕𝜺t−2 𝒆 = 𝛄S,t (𝑫 ( 𝜕𝑠 − 𝜕𝑠 )) 0 = 𝛄S,t

𝜕𝑯St 𝜕𝑠

(98)

□

Rule 2. Given a sequence of load steps L={1, 2,…, t−2, t−1, t,…, N}. If steps t−1 and step t both are plastic load steps, 𝐚t−1 = 𝐚t , and d𝐚t−1⁄d𝝈t−1 = 0, then skipping step t−1, i.e. the load steps contain only S={1, 2,..., t−2, t,…, N} will not change the sensitivity result. Proof: With the same discussion at the beginning of the proof for the previous rule, it is sufficient to prove the same three items: 22

1. 𝝀Lt−1 = 0 2. 𝝀Sn = 𝝀Ln and 𝜸Sn = 𝜸Ln (for n ≤ t − 2) 3. 𝜸Lt−1

T 𝜕𝑯L t−1 𝜕𝑠

+ 𝜸Lt

T 𝜕𝑯L t 𝜕𝑠

= 𝜸St

T 𝜕𝑯St 𝜕𝑠

The first item is obtained from Property 2 in Section 4.1. The second item will be proven in a similar way as in the previous proof. Firstly, Eqs. (89) and (92) still hold, and substitute Eq. (92) into Eq. (89) yields T

𝝀Lt−2

=

−1 𝑲St−2

𝜕𝑯St−2 𝜕𝑯St−2 ∙( 𝜕𝑼t−2 𝜕𝑽t−2

−T

𝜕𝑯Lt−1 𝜕𝑽t−2

T

T

𝜸Lt−1

𝜕𝑯St−2 L + 𝜸 ) 𝜕𝑼t−2 t−1

(99)

Since T N 𝑒 L 𝜕𝑯St 𝜕𝑯St−2 =− = (𝑫 (𝑩 + 𝑼t−2 𝑩 )) 𝜕𝑼t−2 𝜕𝑼t−2 0

(100)

According to Property 2 in section 4.1, it has T T 𝜸S𝛔,t 𝜕𝑯St−2 L 𝜕𝑯St T 𝑒 𝜸S𝛔,t L T N 𝜸 = ((𝑩 + 𝑼t−2 𝑩 ) 𝑫 0) ( ) = − ( ) 𝜕𝑼t−2 t−1 𝜕𝑼t−2 𝜸S𝜀𝑝 ,t 0 equ

=

T 𝜕𝑯St − 𝜸S 𝜕𝑼t−2 t

(101)

Additionally, L

T

𝜕 𝑯t−1 L 𝑰 𝜸 =( 𝐚t−1 T 𝑫𝒆 𝜕 𝑽t−2 t−1

T 𝜸S𝛔,t 𝜕𝑯St 0 𝜸S𝛔,t 𝑰 0 )( )=( T 𝒆 )( )= 𝜸S 0 𝐚t 𝑫 0 𝜸S𝜀𝑝 ,t 𝜕𝑽t−2 t 0

(102)

equ

Substituting Eqs.(101) and (102) into Eq.(99), it derives T

𝝀Lt−2

=

−1 𝑲St−2

𝜕𝑯St−2 𝜕𝑯St−2 ∙( 𝜕𝑼t−2 𝜕𝑽t−2

−T

T

T

𝜕𝑯St 𝜕𝑯St − ) 𝜸St = 𝝀St−2 𝜕𝑽t−2 𝜕𝑼t−2

(103)

The same as Eq. (96), it further leads to 𝜸Lt−2 = 𝜸St−2

(104)

𝝀Sn = 𝝀Ln and 𝛄Sn = 𝛄Ln (for n ≤ t − 2)

(105)

and it is followed by

23

According to Property 2, Eqs. (47) and (84), the third item is obtained via the following deduction 𝜸Lt−1

L T 𝜕𝑯t−1

𝜕𝑠

+ 𝜸Lt

L T 𝜕𝑯t

𝜕𝑠

T 𝑆 𝜕𝜺t−2 𝜕𝜺t 𝜕𝜺t−1 𝜸𝛔,t 𝑆 𝒆 𝜸𝛔,t 𝑫 ( − ) 𝑫𝒆 ( − ) =( ) ( ) + ( ) ( 𝑆 𝜕𝑠 𝜕𝑠 𝜕𝑠 𝜕𝑠 ) 𝜸 𝑝 0 𝜀equ ,t 0 0 T 𝑆 𝜕 𝜺 𝜕 𝜺 𝜸𝛔,t t t−2 𝒆 T 𝜕𝑯S =( 𝑆 ) (𝑫 ( 𝜕𝑠 − 𝜕𝑠 )) = 𝜸St 𝜕𝑠t 𝜸𝜀𝑝 ,t 0

𝜕𝜺t−1

T

(106)

equ

□ The first rule above is applicable to all types of elements. It states that among a sequence of load steps, all elastic load steps except for the last step can be ignored in the sensitivity analysis procedure. The second rule above states that for two consecutive plastic steps t-1 and t, as long as certain conditions are satisfied, the former plastic step can be skipped. These conditions include 𝐚t−1 = 𝐚t , and d𝐚t−1⁄d𝝈t−1 = 0. For bar elements that undergo only tension or compression, the flow vector is 1 in tension 𝐚={ −1 in compression

(107)

Further, as long as there is no switch from tension to compression or vice versa in two consecutive plastic steps, the derivative of the flow vector is zero, i.e., d𝐚 d𝛔

=𝟎

(108)

Therefore, the prerequisites of Rule 2 are strictly satisfied with bar truss structures when all elements keep their tension or compression status in two adjacent plastic steps. Thus, sensitivity analysis with condensed load steps via Rule 1 and Rule 2 should generate exactly the same results as using all load steps. However, for general two- and three-dimensional elements, these two requirements are difficult or impossible to be fulfilled. Take a general three-dimensional element with the von Mises yield criterion and associated flow rule as an example. The associated flow vector is [35] 𝐚=

3 [s , s , s , s , s , s ]T 2|𝝈| 11 22 33 12 13 23

where sij represent deviatoric stress. The derivative of flow vector with respect to stress is

24

(109)

d𝐚 d𝛔

=

1 −0.5 −0.5 1 −0.5 −0.5

1

| 𝝈|

−0.5 −0.5 1 3

([

− 𝐚 ∙ 𝐚T 3

3]

(110) )

which is always non-zero. Therefore, Rule 2 cannot be directly applied to general structures. The two prerequisites physically ask for one thing, i.e., that the flow vectors in two consecutive plastic steps do not change much if at all. Based on this principle, the following empirical rule is proposed, which relaxes the requirements for condensing plastic load steps. Empirical rule: Given a sequence of load steps L = {1, 2,., t−2, t−1, t,…, N}. If steps t−1 and t are both plastic steps and in a monotonic loading procedure, then step t−1 can be skipped in the sensitivity analysis. This empirical rule applies to all element types. The aim of this rule is to reduce the number of load steps to the extent possible in the sensitivity analysis. By doing so, some accuracy will be sacrificed in exchange for efficiency. In the next section, the effectiveness of this empirical rule and its influence on the accuracy of the sensitivity results will be investigated through numerical examples. 5. Numerical examples of sensitivity analysis with condensed load steps In this section, firstly, the adjoint sensitivity analysis procedure is verified through a bar truss example. Next, the load step condensation rules are demonstrated with the same structure using multiple loading conditions. Lastly, the efficiency and accuracy of the condensation techniques, especially the empirical rule, are investigated using two solid structures with three-dimensional tetrahedral elements under complex load cases. 5.1 Verification of adjoint sensitivity analysis

Y position (mm)

In this subsection, a bar truss structure is employed to verify the sensitivity analysis with proposed adjoint variable method. As depicted in Figure 3, the beam-shaped truss structure is composed of 100 bar elements. One end of the structure is fixed while an external force is applied to the other end. Here, a bilinear elastoplastic material with isotropic hardening and associated flow rule is assumed. The Young’s modulus of the material is 210GPa, with initial yield stress of 235MPa and the plastic modulus of 100GPa.

Truss beam structure 10 0 -10 0

F 20

40 60 X position (mm) 25

80

100

Figure 3. Configuration of a 100-bar truss structure with design nodes indicated by red dots A loading-unloading history applied to the structure is presented in Figure 4(a), while the finite element analysis results, specifically the plastic strain history of an element at the fixed end, is shown in Figure 4(b). During the first few load steps, the structure behaves elastically. Then, plastic strain increases until the maximum load is reached at step 40. Finally, the load is gradually released. The deformation of the structure at the fully loaded step (i.e., step 40) and fully unloaded step (i.e., step 80) are depicted in Figure 4(c) and 4(d), respectively. Loading History

Plastic strain at fixed end

40

0.014 0.012

Plastic strain

Force (N)

30

20

0.01 0.008 0.006 0.004

10

0.002 0 0

20

40 Load step

60

0 0

80

(a) Loading history on the structure

20

2D bar truss, loading end

80

2D bar truss, fully unloading 60

50

50

40

40

Y position (mm)

Y position (mm)

60

(b) Plastic strain with respect to the loading history

60

30 20 10 0 -10

40 Load step

30 20 10 0

20

40 60 X position (mm)

-10

80

(c) Undeformed and deformed structure at the fully loaded step (i.e. load step 40)

20

40 60 X position (mm)

80

(d) Undeformed and deformed structure at the fully unloaded step (i.e. load step 80)

Figure 4. Loading history and finite element analysis results of the bar truss structure: (a) the loading-unloading history on the structure, (b) plastic strain history of an element at the fixed end, (c) deformation of the structure at the fully loaded step and (d) remaining plastic deformation after the fully unloaded step. Design nodes of the structure locate in the lower side of the structure are indicated by red dots in Figure 3. The vertical coordinates of these nodes are taken as design variables. Here, the residual vertical displacement after unloading is selected as the system response. Figure 5(a) presents sensitivity results obtained with the adjoint variable method and the global finite difference method, which is employed as reference results. 26

In the global finite difference method, the central scheme finite differencing is employed to improve the accuracy of the results. The perturbation size used in the finite differencing also influences the results. Therefore, several perturbation sizes ranging from 10-3mm to 10-10mm are tested. It is shown that converged results are achieved for this example with perturbation sizes from 10-4mm to 10-9mm. In Figure 5, the finite difference results with a perturbation size of 5×10-5mm are presented. Sensitivity comparison

adjoint sensitivity/global FD

3

1.005 adjoint global FD

adjoint/global FD

relative sensitivity

Sensitivity (mm/mm)

4

2 1

1

0 -1 0

20 40 60 80 x-coordinate of design node (mm)

100

0.995 0

(a) Sensitivity comparison of relative value

20 40 60 80 x-coordinate of design node (mm)

100

(b) Sensitivity comparison of relative value

Figure 5. A comparison of sensitivity results Further, the relative sensitivity is defined as relative sensitivty =

sensitivity with adjoint variable method sensitivity with global finite differencing

(111)

If the relative sensitivity equals one, then it means that the two results are exactly the same. The relative sensitivities for this example are shown in Figure 5(b). It shows that the sensitivities of adjoint variable method match perfectly with the reference results, from both real value and relative value points of view. 5.2 Demonstration of condensation rules with a bar truss structure In this subsection, to thoroughly demonstrate the efficiency and accuracy of the condensation rules, three cases are presented, in which three different loading histories are applied to the structure in Figure 3. In the first case, a vertical force load is applied to the free end of the structure following a total number of 35 load steps, as depicted in Figure 6. The loading history is composed of the following three stages: (1) an initial loading stage from step 1 to step 12; (2) a partially unloading stage from step 13 to step 22; and (3) a reloading stage from step 23 to step 35. Note that the final load level in step 35 is larger than the intermediate level in step 12. Based on finite element analysis, elastic and plastic steps are distinguished in Figure 6 by using green discs to represent elastic load steps and red diamonds to illustrate plastic steps.

27

15

35

Force (N)

12 10

5 Loading History elastic step plastic step

0 0

5

10

15 20 Load step

25

30

35

Figure 6. Load history of the first load case, with a pentagram depicting the condensed load step Following condensation Rule 1, all elastic load steps can be skipped in the sensitivity analysis. For the plastic steps, since the force load is monotonically increasing during the loading and reloading stages, the state of each bar remains the same, either in tension or compression. Therefore, the former step of each pair of consecutive plastic steps can be neglected according to Rule 2, and thus only step 12 and step 35 are left. Note that step 12 and step 35 are now two consecutive plastic steps, and they are in the same direction and there is no change in tension and compression for any bar in these two steps. Therefore, step 12 can also be skipped. Given this example, only the last load step is required for the sensitivity analysis as listed in Table 1 and indicated with a pentagram in Figure 6. Table 1. Condensed load step for the first case Condensed load step Corresponding load step 1 35

Load 15N

The sensitivities of two system responses, i.e., the vertical displacement at the loading node and the equivalent plastic strain of an element at the fixed end are evaluated. In Figure 7, sensitivity results calculated by the condensed load step are compared with those of using full load steps. Both the real values and relative values are presented. The relative value is defined as the division of the sensitivity using the condensed load steps by the sensitivity using full load steps, i.e. relative value =

sensitivitycondensed load steps sensitivityfull load steps

(112)

As expected, results show that sensitivities with the condensed load step are exactly the same as those using all load steps for both responses.

28

3.5

1.001 full load steps condensed LS

Condensed:L S F ull:L S 1.0005

2.5

relative value

Sensitivity (mm/mm)

3

2 1.5 1

1

0.9995

0.5 0 0

20 40 60 80 x-coordinate of design node (mm)

0.999 0

100

(a) Sensitivity of the displacement response

20 40 60 80 x-coordinate of design node (mm)

100

(b) Relative value of sensitivities with different load steps for the displacement response

-5

0

x 10

-1

Condensed:L S F ull:L S 1.0005

-2

relative value

Sensitivity (/mm)

1.001

full loading steps condensed LS

-3 -4

1

0.9995

-5 -6 0

20 40 60 80 x-coordinate of design node (mm)

(c) Sensitivity of the plastic strain response

0.999 0

100

20 40 60 80 x-coordinate of design node (mm)

100

(d) Relative value of sensitivities with different load steps for the plastic strain response

Figure 7. A comparison of sensitivities for the first load case In the second load case, the vertical force load follows a total number of 50 load steps, as depicted in Figure 8. Here, the loading history consists of the following four stages: (1) an initial loading stage from step 1 to step 10; (2) a full unloading stage from step 11 to step 20; (3) a reverse loading stage from step 21 to step 35; and (4) an unloading stage from step 36 to step 50. In this case, the maximal force of reverse loading is larger than that of the initial loading stage.

29

10

10

Loading History elastic step plastic step

Force (N)

5 0

50

-5 -10 -15 0

35 10

20 30 Load step

40

50

Figure 8. Loading history of the second load case, with pentagrams depicting the condensed load steps According to Rule 1, all elastic load steps except for the last step can be skipped in the sensitivity analysis. Since the force load is monotonically increasing during the loading and reverse loading stages, the state of tension or compression remains the same for each bar in either stage. Therefore, the former one in each pair of consecutive plastic steps can be ignored. As such, in this case, only plastic step 10, plastic step 35, and elastic step 50 are left. Step 10 and step 35 are now two consecutive plastic steps. In contrast to the first example, the external forces in these two steps are in the opposite direction, thus some of the bars change from tension to compression, while some change from compression to tension. Therefore, step 10 cannot be skipped. In summary, step 10, step 35, and step 50 are needed in the sensitivity analysis in this case. These steps are listed in Table 2, and indicated in Figure 8 with pentagrams. Table 2. Condensed load steps of the second load case Condensed load steps Corresponding load step 1 10 2 35 3 50

Load 10N −15N 0N

As in the first load case, the two system responses here are the same. Sensitivity results calculated by the condensed load steps are compared with those obtained using all load steps. Again, both the real values and relative values are presented in Figure 9. It is shown that the sensitivities with the condensed load steps perfectly match those with all load steps.

30

1.001

Condensed:L S F ull:L S

full LS condensed LS

0.4

1.0005

relative value

Sensitivity (mm/mm)

0.6

0.2 0

1

0.9995

-0.2 -0.4 0

20 40 60 80 x-coordinate of design node (mm)

0.999 0

100

(a) Sensitivity of displacement response

20 40 60 80 x-coordinate of design node (mm)

100

(b) Relative value of sensitivities with different load steps for displacement response

-5

1

x 10

0.5

Condensed:L S F ull:L S 1.0005

relative value

Sensitivity (/mm)

1.001 full LS condensed LS

0 -0.5

1

0.9995

-1 -1.5 0

20 40 60 80 x-coordinate of design node (mm)

(c) Sensitivity of plastic strain response

0.999 0

100

20 40 60 80 x-coordinate of design node (mm)

100

(d) Relative value of sensitivities with different load steps for plastic strain response

Figure 9. A comparison of sensitivities for the second load case In the third case, aside from the vertical load, another load in the horizontal x-direction is simultaneously applied on the free end of the structure. As shown in Figure 10, both loads follow 50 load steps. In spite of the load history of each individual force, they together define the following three stages: (1) an initial plastic loading stage from step 1 to step 5; (2) an elastic stage from step 6 to step 38; and (3) a vertical plastic loading stage from step 39 to step 50.

31

5

5000

40 Loading History elastic step plastic step

4000

Loading History elastic step plastic step

35

50

3000

Force (N)

Force (N)

30

2000

25 20 15

5

10 1000

0 0

10

20 30 Load step

40

50

5

50

0 0

(a) Loading history in horizontal x-direction

10

20 30 Load step

40

50

(b) Loading history in vertical y-direction

Figure 10. Loading history of the third load case, with pentagrams depicting the condensed load steps In the initial loading stage, bars are primarily dominated by the monotonically increasing extension force in the x-direction, thus there is no switch of their states. According to Rule 2, the first four steps can be skipped here. In the vertical plastic loading stage, tension or compression is determined by the bending moment caused by the vertical force. Because the vertical force is monotonic, all steps except for the final step can be skipped in the sensitivity analysis. For the two steps that are left, the horizontal force leads to a tensile status for most of the bars at step 5. At step 50, the vertical force leads to a bending deformation of the structure, thus some bars change to the compression states. Therefore, both step 5 and step 50 should remain in the sensitivities analysis. Table 3 lists the condensed load steps, and these two steps are indicated with pentagrams in Figure 10. Table 3. Condensed load steps for the third load case Condensed load steps Corresponding load step Load 1 5 Fx = 5000N, Fy = 5N 2 50 Fx = 0N, Fy = 40N In this third load case, the two system responses are the same as in the previous two cases; sensitivity results are presented in Figure 11 with real values and relative values. As expected, the sensitivities with the condensed load steps perfectly match those using full load steps. 0.5

Condensed:L S F ull:L S 1.0005

-0.5

relative value

Sensitivity (mm/mm)

0

1.001 full LS condensed LS

-1 -1.5 -2

1

0.9995

-2.5 -3 0

20 40 60 80 x-coordinate of design node (mm)

0.999 0

100

32

20 40 60 80 x-coordinate of design node (mm)

100

(a) Sensitivity of displacement response

(b) Relative value of sensitivities with different load steps for displacement response

-4

5

x 10

0

Condensed:L S F ull:L S 1.0005

relative value

Sensitivity (/mm)

1.001 full LS condensed LS

-5 -10

1

0.9995

-15 -20 0

20 40 60 80 x-coordinate of design node (mm)

(c) Sensitivity of plastic strain response

0.999 0

100

20 40 60 80 x-coordinate of design node (mm)

100

(d) Relative value of sensitivities with different load steps for plastic strain response

Figure 11. A comparison of sensitivities for the third load case 5.3 Demonstration of condensation rules with solid structures In this subsection, the condensation rules, especially the empirical rule, are demonstrated using two solid structure examples, where the prerequisites of Rule 2 are in no case satisfied. Firstly, Figure 12 illustrates a solid cantilever beam structure with 286 tetrahedral elements. The material model and properties are the same as those described in Section 5.1, with Poisson’s ratio equals to 0.3. The vertical coordinates of the design nodes, which are indicated by red dots in the figure, are chosen as design variables. The total number of design variables is 31.

F

X Z

Y Figure 12. Solid cantilever beam structure with a length of 300 mm, a width of 15 mm, and a height of 15 mm Two force loads in the vertical y-direction and horizontal z-direction are applied simultaneously to the structure. The vertical force follows a constant loading-partially unloadingholding constant history, whereas the horizontal force follows a pure loading procedure that starts from the load step when the vertical force begins to keep constant. They together define the following three stages: (1) an initial plastic loading stage from step 1 to step 5; (2) an elastic stage from step 6 to step 8; and (3) a plastic loading stage from step 9 to step 18. The loading history is 33

shown in Figure 13(a), the deformation and contour of the equivalent plastic strain at the last step are depicted in Figure 13(b). 10000

Force (N)

8000

Fy Fz elastic step plastic step

6000

4000

2000

0 0

5 2

4

18 6

8 10 Load step

12

14

16

18

(a) Loading history in the vertical y-direction and horizontal z-direction

(b) Undeformed structure and contour plot of equivalent plastic strain on deformed structure at the last load step

Figure 13. Loading history of the solid beam example and finite element analysis results In the sensitivity analysis, the first four steps are skipped according to the empirical rule. The three elastic steps are skipped due to Rule 1. In the second plastic loading stage, the resultant force of both loads changes its direction at each step; however, the force in the y-direction remains constant and the force in the z-direction monotonically increases. Therefore, the empirical rule can still apply, and only the last step is preserved for sensitivity analysis. In summary, the condensed load steps contain two steps as highlighted in Figure 13(a) and listed in Table 4. Table 4. Condensed load steps for the first solid structure example Condensed load steps Corresponding load step Force 1 5 Fy = 5000N, Fz = 0N 2 18 Fy = 2000N, Fz = 10000N The sensitivities of two responses are investigated. One is the average displacement in the vertical y-direction at the loading end, which is influenced by both loads in the y and z direction. The other is the average equivalent plastic strain of elements at the fixed end. Sensitivities calculated using the condensed load steps are compared with those using full load steps in Figure 14 from both the real value and relative value points of view.

34

0.3 0.25

Condensed:L S F ul l:L S 1.05

0.2

relative value

Sensitivity (mm/mm)

1.1

full loading steps condensed LS

0.15 0.1 0.05

1

0.95

0 -0.05 0

50

100 150 200 250 x-coordinate of design node (mm)

0.9 0

300

(a) Sensitivity of average vertical y-displacement at the loading end

50 100 150 200 250 x-coordinate of design node (mm)

300

(b) Relative value of sensitivities for displacement response

-4

2

x 10

1.3

Condensed:L S F ul l:L S 1.2

relative value

Sensitivity (/mm)

0

-2

-4

1 0.9

-6

-8 0

1.1

full loading steps condensed LS 50 100 150 200 250 x-coordinate of design node (mm)

0.8 0

300

(c) Sensitivity of average equivalent plastic strain at the fixed end

50 100 150 200 250 x-coordinate of design node (mm)

300

(d) Relative value of sensitivities for equivalent plastic strain response

Figure 14. Comparison of sensitivities of the solid beam structure example As per usual, the sensitivities of the two strategies match quite well with one another from the real value point of view as shown in Figure 14(a) and 14(c). The prerequisites regarding flow vectors are not satisfied with solid elements, and the violations are more severe in the second loading stage since the resultant load and thus flow vectors change their directions at every load step of the second stage. Therefore, deviations of results using the two strategies are observed as expected through the relative value results shown in Figure 14(b) and 14(d). Further, absolute percentage errors of sensitivities are calculate from the relative values following Eq. (113). sensitivitycondensed load steps absolute percentage error = | − 1| × 100% sensitivityfull load steps

(113)

The results show that the maximal error of displacement sensitivities is 8.4% with an average error of 2.7%. The maximal error of the equivalent plastic strain sensitivities is 16.3% with an average error of 5.9%.

35

It should be noted that, by quantifying the influence of design variables on responses, sensitivities are usually employed to find vectors along which design variables should be adjusted so as to improve responses in gradient-based optimization algorithms. These vectors are determined more by sensitivities with large values than very small sensitivities. Therefore, in practical applications, more attention should be paid to the accuracy of large sensitivities, i.e. sensitivities that are in the same order of magnitude as the maximal sensitivity value. Additionally, the percentage error is biased larger when the numerator (here refers to the difference between sensitivity with condensed load steps and sensitivity with full load steps) keeps the same but the denominator (here refers to the value of sensitivity with all load steps) gets smaller. Due to these two reasons, it is more meaningful to analyze the percentage errors together with their corresponding real sensitivity values. In Table 5, the total of 31 design variables are divided into three groups according to their displacement sensitivity values as shown in the first column. The first group contains design variables whose absolute displacement sensitivity is no less than one tenth of the maximal absolute sensitivities of all design variables. The second group contains design variables whose absolute sensitivity lies in one hundredth to one tenth of the maximal absolute sensitivity. The last group contains design variables whose absolute sensitivity is less than one hundredth of the maximum. The number of design variables in each group is listed in the second column of Table 5. The maximal absolute percentage error of sensitivities in each group is listed in the third column, and the average of absolute percentage errors of sensitivities in each group is presented in the last column. Table 5. Absolute percentage errors of displacement sensitivities for solid beam example Range of absolute sensitivity ( × max(|Sensitivity|))

Number of design variables

Maximal absolute percentage error

Average of absolute percentage error

[0.1,1] [0.01,0.1) [0,0.01)

18 7 6

3.8% 5.3% 8.4%

1.7% 3.9% 4.2%

As shown in Table 5, the maximal absolute percentage error in the first group is 3.8% with the average error of 1.7%. Both the maximal error and average error increase in the second and the third group with smaller sensitivity values. The maximal error in the second group is 5.3% with the average error of 3.9%. The absolute percentage errors of sensitivities for the equivalent plastic strain response are analyzed in the same way and presented in Table 6. It shows that the maximal absolute percentage error of sensitivities in the first group is 6.6%, and the average error is 4.1%. The maximal error in the second group is 11.0% with the average error of 5.3%. Table 6. Absolute percentage errors of equivalent plastic strain sensitivities for solid beam example Range of absolute sensitivity ( × max(|Sensitivity|))

Number of design variables

Maximal absolute percentage error

Average of absolute percentage error

[0.1,1] [0.01,0.1) [0,0.01)

3 7 21

6.6% 11.0% 16.3%

4.1% 5.3% 6.4%

Secondly, the applicability of load step condensation for solid elements is discussed through an example of control arm structure in an automotive suspension system. The finite element model of the control arm is presented in Figure 15(a). Due to the symmetry of the structure about 36

the XY plane, only one half of the model is established. The model consists of 8925 4-node tetrahedral elements and 2744 nodes. Except for the symmetric boundary condition, upper left end and upper right end are both fixed. Forces are applied uniformly on inner surface nodes at lower bearing joint. The design variables are the z-coordinates of 612 nodes that are indicated by red dots in Figure 15(b). These nodes lie on the surface between upper ends and the lower bearing joint. A bilinear elastoplastic material follows isotropic hardening and the associated flow rule is assumed. Material parameters are listed in Table 7.

Upper ends, fixed Symmetric about XY plane

Fy

Lower bearing joint

Fx (a) Model description and boundary conditions

(b) Design nodes (highlighted with red dots)

Figure 15. Finite element model of a control arm structure Table 7. Material parameters of the control arm example Young’s modulus (GPa) Plastic modulus (GPa) Initial yield stress (MPa) 210 50 350

Poisson’s ratio 0.3

A total number of 15 load steps as depicted in Figure 16 are applied on inner surface nodes in XY plane. The structure is gradually loaded in both horizontal X and vertical Y directions from step 1 to step 6, the structure behaves plastically in this stage. Both forces are partially unloaded in step 7 and step 8. From step 9 to step 12, the structure is horizontally loaded in the opposite direction. In the meantime, the vertical load increases in step 10 and remains the same in step 10 and step 11. Both forces are elastically unloaded from step 13 to step 15.

37

Fy (kN)

15

Loading history elastic step plastic step

0

1

2

3

14

4

-50

13

9

8

5 7

-100

12

-250

-200

-150

6

10

11 -100

-50 Fx (kN)

0

50

100

150

Figure 16. Loading history on the control arm structure The contour of equivalent plastic strain at the last step is depicted in Figure 17(a). The maximum value point lies on the side surface between lower bear joint and upper right end. The average of equivalent plastic strains of 93 elements in this area is defined as the first response. The contour of displacement at the last step is depicted in Figure 17(b). The average residual displacement of inner surface nodes of lower bearing joint after unloading are defined as the second response.

max. equivalent plastic strain

(a) Contour of equivalent plastic strain at final step

(b) Contour of displacement (mm) on both original and deformed structure at the last step

Figure 17. Finite elements analysis results of the control arm example According to Rule 1, all elastic load steps are skipped in the sensitivity analysis. Further, plastic steps 1 to 5 and steps 9 to 11 are also skipped following the empirical rule since they lie in monotonically loading procedures with respect to neighbouring load steps. Finally, the three condensed load steps in sensitivity analysis are indicated by pentagrams in Figure 16 and listed in Table 8. Table 8. Condensed load steps in sensitivity analysis for the control arm example Condensed load steps Corresponding load step Force 38

1 2 3

Fx = 100kN, Fy = −100kN Fx = −200kN, Fy = −100kN Fx = 0kN, Fy = 0kN

6 12 15

In Figure 18 and Figure 19, the contour plots of sensitivities for both responses are depicted on the design surface respectively. In each of these two figures, the sensitivity results obtained with condensed load steps are presented on the left and results obtained with all load steps are presented on the right. The negative value of sensitivity means that if a design node move in +Z direction, i.e., the structure increases its thickness in Z-direction at that node point, then the response will decrease in value. Oppositely, positive value of sensitivities mean responses will increase in value if structure thickness is increased. The sensitivities of the first response in Figure 18 show that, if thickness is increased close to the maximal equivalent plastic strain point, then the average value around that point can be decreased. This is caused by the decrease of maximal equivalent plastic strain due to the increase in thickness. Since equivalent plastic strain is directly related to equivalent stress, therefore, similar to stress responses, the change of structure far away from the critical area will have little influence on the equivalent plastic strain of critical area. The almost zero sensitivities far away from the critical area are in consistent with this intuition.

(a) Sensitivities of equivalent plastic strain response with condensed load steps

(b) Sensitivities of equivalent plastic strain response with all load steps

Figure 18. Contour plot of sensitivities of equivalent plastic strain response. Sensitivity results with condensed load steps are depicted on the left and results with all load steps are presented on the right. The sensitivity plots of the second response in Figure 19 show that the residual displacement of lower bearing joint could be reduced through the increase of thickness in the area close to side faces where plastic deformation is large. This is caused by the decrease of plastic strains near side faces due to the increase in thickness.

39

(a) Sensitivities of displacement response (b) Sensitivities of displacement response with condensed load steps with all load steps Figure 19. Contour plot of sensitivities of displacement response. Sensitivity results with condensed load steps are depicted on the left and results with all load steps are presented on the right.

By comparing the results of Figure 18(a) to Figure 18(b), and Figure 19(a) to Figure 19(b), it can be seen that the sensitivities calculated with the condensed load steps match well with those with all load steps. The same discussion about absolute percentage errors as in the solid beam example is again carried out here. Table 9 and Table 10 present the sensitivity errors for the equivalent plastic strain response and the displacement response, respectively. The results of equivalent plastic strain response show that the maximal absolute percentage error of sensitivities that are no less than one tenth of maximum is 4.8% with the average error of only 0.7%. The maximal absolute percentage error of displacement sensitivities in the first group is 13.2% with the average error of 4.4%. The maximal error and average error increase in the second and the third group with smaller sensitivity values. The maximal error in the third group of both responses are over 100%. However, as discussed, the inaccuracy of sensitivities with small values is insignificant when they are employed in an optimization procedure. Table 9. Absolute percentage errors of equivalent plastic strain sensitivities for control arm example Range of absolute sensitivity ( × max(|Sensitivity|))

Number of design variables

Maximal absolute percentage error

Average of absolute percentage error

[0.1,1] [0.01,0.1) [0,0.01)

64 79 469

4.8% 12.6% >100%

0.7% 2.8% 12.7%

Table 10. Absolute percentage errors of displacement sensitivities for control arm example Range of absolute sensitivity ( × max(|Sensitivity|))

Number of design variables

Maximal absolute percentage error

Average of absolute percentage error

[0.1,1] [0.01,0.1) [0,0.01)

305 258 49

13.2% 23.3% >100%

4.4% 6.8% 17.6%

40

The results of both the solid beam example and the control arm example show that the absolute percentage errors of large sensitivities (no less than one tenth of maximal sensitivity value) are significantly smaller than the errors of sensitivities with small values. According to Table 5, Table 6, Table 9 and Table 10, the average absolute percentage errors of large sensitivities are less than 5% in all cases. The average absolute percentage errors of sensitivities lie in one hundredth to one tenth of maximal value are less than 7% in all cases. They demonstrate that the condensation rules especially the empirical rule provide a reliable procedure for load steps reduction in adjoint sensitivity analysis. Following these rules, satisfying sensitivities results for optimization could be obtained with significantly reduced number of load steps. 6. Conclusions In this paper, efficient sensitivity analysis with adjoint variable method for simultaneous elastoplasticity and geometric nonlinearities is investigated. The adjoint sensitivity formulation introduces a sequence of adjoint variables, which are obtained via a backward procedure due to the path dependence of material nonlinearity. As a result, computational costs grow in proportion to the total number of load steps, which does not scale well. To address this cost and improve efficiency, several properties regarding adjoint variables are proposed; based on these properties, it is theoretically proven that elastic load steps except for the last step can be skipped in the sensitivity analysis. It is also theoretically proven that, under certain conditions, the former of two adjacent plastic load steps can also be skipped. When the conditions are not satisfied, an empirical rule is proposed in the case of general types of elements. Following these condensation rules, the total number of load steps required for sensitivity analysis can be significantly reduced. Both a theoretical discussion and specific numerical examples of a bar truss structure with several typical loading histories show that these strategies are exact for bar elements. Two solid structure examples with solid tetrahedral elements and complicated loading histories also demonstrate that, these load steps condensation strategies can highly efficiently obtain sensitivities while keeping high accuracy for general problems. 7. Acknowledgements This research is supported by the European Commission under the 7th Framework Program through the Marie Curie Action: “Industry-Academia Partnerships and Pathways (IAPP)” of the “Capacities” Program. Project: LaSciSO, Grant number 285782. 8. References [1] Y. H. Park and K. K. Choi, Shape design sensitivity analysis of nonlinear 2-D solids with elasto-plastic material, Structural Optimization, Vol. 18(4), pp.236-246, 1999. [2] K. K. Choi and N. H. Kim, Structural Sensitivity Analysis and Optimization 2: Nonlinear Systems and Applications, New York, Springer, 2005. [3] C. Wu and J. S. Arora, Design sensitivity analysis and optimization of nonlinear structural response using incremental procedure, AIAA Journal, Vol. 25(8), pp.1118-1125, 1987. [4] S. Schwarz, K. Maute, and E. Ramm, Topology and shape optimization for elastoplastic structural response, Computer Methods in Applied Mechanics and Engineering, Vol. 190(15-17), pp.2135-2155, 2001. [5] Y. H. Park and K. K. Choi, Design Sensitivity Analysis of Truss Structures with Elastoplastic Material, Mechanics of Structures and Machines: An International Journal, Vol. 24(2), pp.189-216, 1996. [6] J. Kato, H. Hoshiba, S. Takase, K. Terada and T. Kyoya, Analytical sensitivity in topology optimization for elastoplastic composites, Structural and Multidisciplinary Optimization, Vol. 52(3), pp.507-526, 2015. 41

[7] C. A. Vidal, H. S. Lee, and R. B. Haber, The consistent tangent operator for design sensitivity analysis of historydependent response, Computing Systems in Engineering, Vol. 2(5-6), pp.509-523, 1991. [8] C. A. Vidal and R. B. Haber, Design sensitivity analysis for rate independent elastoplasticity, Computer Methods in Applied Mechanics and Engineering, Vol. 107(3), pp.393-431, 1993. [9] M. Ohsaki and J. S. Arora, Design Sensitivity Analysis of Elastoplastic Structures, International Journal for Numerical Methods In Engineering, Vol. 37(5), 737-762, 1994. [10] S. Schwarz and E. Ramm, Sensitivity analysis and optimization for non-linear structural response, Engineering Computations, Vol. 18(3/4), pp.610-641, 2001. [11] A. Habibi and H. Moharrami, Nonlinear sensitivity analysis of reinforced concrete frames, Finite Elements in Analysis and Design, Vol. 46(7), pp.571-585, 2010. [12] M. Kleiber and P. Kowalczyk, Sensitivity analysis in plane stress elasto-plasticity and elasto-viscoplasticity, Computer Methods in Applied Mechanics and Engineering, Vol. 137(3-4), pp.395-409, 1996. [13] S. Holopainen, Parameter sensitivity of anisotropic elasto-plastic shell, 6th World Congress on Structural and Multidisciplinary Optimization, Rio de Janeiro, 2005, Brazil. [14] K. Wisniewski, P. Kowalczyk, and E. Turska, On the computation of design derivatives for Huber-Mises plasticity with non-linear hardening, International Journal for Numerical Methods in Engineering, Vol. 57(2), pp.271-300, 2003. [15] Q. Gu, J. P. Conte, A. Elgamal, and Z. Yang, Finite element response sensitivity analysis of multi-yield-surface J2 plasticity model by direct differentiation method, Computer Methods in Applied Mechanics and Engineering, Vol. 198(30-32), pp.2272-2285, 2009. [16] S. Schwarz, Sensitivitätsanalyse und Optimierung bei nichtlinearem Strukturverhalten, Doktor Dissertation, Institut für Baustatik, Universität Stuttgart, 2001. [17] C. B. W. Pedersen, Topology optimization of 2D-frame structures with path-dependent response, International Journal for Numerical Methods in Engineering, Vol. 57(10), pp.1471-1501, 2003. [18] A. Chattopadhyay and R. Guo, Structural design sensitivity analysis for composites undergoing elastoplastic deformation, Mathematical and Computer Modelling, Vol. 22(2), pp.83-105, 1995. [19] C. O. Spivey and D. A. Tortorelli, Tangent operators, sensitivity expressions, and optimal design of non-linear elastica in contact with applications to beams, International Journal for Numerical Methods in Engineering, Vol. 37(1), pp.49-73, 1994. [20] N. H. Kim, K. K. Choi, and J. S. Chen, Shape Design Sensitivity Analysis and Optimization of Elasto–Plasticity with Frictional Contact, AIAA JOURNAL, Vol. 38(9), pp.1742-1753, 2000. [21] S. Stupkiewicz, J. Korelc, M. Dutko and T. Rodi, Shape sensitivity analysis of large deformation frictional contact problems, Computer Methods in Applied Mechanics and Engineering, Vol. 191(33), pp.3555-3581, 2002. [22] Y. S. Ryu, M. Haririan, C. C. Wu, and J. S. Arora, Structural design sensitivity analysis of nonlinear response, Computers & Structures, Vol. 21(1-2), pp.245-255, 1985. [23] M. Kleiber, H. Antunez, T. D. Hien, and P. Kowalczyk, Parameter Sensitivity in Nonlinear Mechanics, Theory and Finite Element Computations. Wiley: New York, 1997. [24] J. J. Tsay and J. S. Arora, Nonlinear structural design sensitivity analysis for path dependent problems. Part 1: General theory, Computer Methods in Applied Mechanics and Engineering, Vol. 81(2), pp.183-208, 1990. [25] P. Kowalczyk and M. Kleiber, Shape sensitivity in elasto-plastic computations, Computer Methods in Applied Mechanics and Engineering, Vol. 171(3), pp.371-386, 1999. [26] T. Hisada, Recent Progress in Sensitivity Nonlinear FEM-Based Sensitivity Analysis, JSME International Journal, Series A, Vol. 38(3), pp.301-310, 1995. [27] S. H. Chung, L. Fourment, J. L. Chenot, and S. M. Hwang, Adjoint state method for shape sensitivity analysis in non-steady forming applications, International Journal for Numerical Methods in Engineering, Vol. 57(10), pp.14311444, 2003. [28] D. E. Smith, Design sensitivity analysis and optimization for polymer sheet extrusion and mold filling processes, International Journal for Numerical Methods in Engineering, Vol. 57(10), pp.1381-1411, 2003. [29] K. Maute, M. Nikbay, and C. Farhat, Sensitivity analysis and design optimization of three-dimensional non-linear aeroelastic systems by the adjoint method, International Journal for Numerical Methods in Engineering, Vol. 56(6), pp.911-933, 2003. [30] P. Michaleris, D. A. Tortorelli, and C. A. Vidal, Tangent operators and design sensitivity formulations for transient non-linear coupled problems with applications to elastoplasticity, International Journal for Numerical Methods in Engineering, Vol. 37(14), pp.2471-2499, 1994. 42

[31] J. B. Cardoso, Structural Design Sensitivity Analysis of Elasti-Plastic History-Dependent Response, 6th World Congress of Structural and Multidisciplinary Optimisation, Rio de Janeiro, 30 May-03 June 2005, Brazil. [32] J. Köbler, Adjungierte Sensitivitätsberechnung und Gestaltoptimierung von Materialien mit nicht-linearem und plastischem Verhalten, Master Thesis, Institut für Technische Mechanik, Karlsruhe Institut für Technology, 2015. [33] P. Pedersen, Analytical stiffness matrices with Green–Lagrange strain measure, International Journal for Numerical Methods in Engineering, Vol. 62(3), pp.334-352, 2005. [34] J. C. Simo, T. J. R. Hughes, Computational Inelasticity, New York, Springer, pp. 12-13, 35-39, 49-51, 127-128, 1998. [35] M. A. Crisfield, Non-linear finite element analysis of solids and structures Vol. 1: Essentials, John Wiley & Sons Ltd, pp. 156-158, 162-168, 178-181, 2000. [36] F. Dunne, N. Petrinic, Introduction to Computational Plasticity, Oxford University Press, pp. 151-159, 2005.

43

Highlights  The adjoint variable method for design sensitivity analysis of isotropic hardening elastoplasticity and geometric nonlinearity is presented.  Two special numerical properties of adjoint variables are proposed and theoretically proven.  Two theoretical rules for load steps reduction in sensitivity analysis are proposed and theoretically proven based on above two properties.  An empirical rule for load steps condensation in practical sensitivity analysis are proposed.  Numerical examples, where bar truss structures and solid structures under complicated multistep load cases are employed. Results show that the method significantly reduces the computational cost in the sensitivity analysis while keeping high accuracy.