Variational sensitivity analysis of elastoplastic structures applied to optimal shape of specimens

Abstract

The aim of this paper is to improve the shape of specimens for biaxial experiments with respect to optimal stress states, characterized by the stress triaxiality. Gradient-based optimization strategies are used to achieve this goal. Thus, it is crucial to know how the stress state changes if the geometric shape of the specimen is varied. The design sensitivity analysis (DSA) of the stress triaxiality is performed using a variational approach based on an enhanced kinematic concept that offers a rigorous separation of structural and physical quantities. In the present case of elastoplastic material behavior, the deformation history has to be taken into account for the structural analysis as well as for the determination of response sensitivities. The presented method is flexible in terms of the choice of design variables. In a first step, the approach is used to identify material parameters. Thus, material parameters are chosen as design variables. Subsequently, the design parameters are chosen as geometric quantities so as to optimize the specimen shape with the aim to obtain a preferably homogeneous stress triaxiality distribution in the relevant cross section of the specimen.

1 Introduction

In biaxial tests the stress state can be characterized by the stress triaxiality η defined by

$$ \eta = \frac{\sigma_{\mathsf{m}}}{\sigma_{\mathsf{eq}}} = \frac{I_{1}(\boldsymbol{\mathsf{\sigma}})}{3 \sqrt{3 J_{2}(\text{dev}\boldsymbol{\mathsf{\sigma}})}}, $$

(1)

with the mean stress $\sigma _{\mathsf {m}} = \frac {1}{3} I_{1}(\boldsymbol {\mathsf {\sigma }})$ and the equivalent von Mises stress $\sigma _{\mathsf {eq}} = \sqrt {3 J_{2}(\text {dev}\boldsymbol {\mathsf {\sigma }})}$. The stress triaxiality classifies the stress state to be tensile if η > 0, compressive if η < 0, or pure shear if η = 0. New specimens for biaxial tests have been developed in Gerke et al. (2017), which widen the range of stress triaxiality compared to traditional specimens. This is crucial to distinguish between different micro-mechanical damage mechanisms that are responsible for material failure. With the new specimens, it is possible to capture various damage mechanisms with only one type of specimen by varying the loading conditions. These specimens shall be used to identify material parameters for complex damage models. Thus, it is important to be able to separate different damage mechanisms. In Gerke et al. (2019), special attention is drawn to the so-called X0-specimen, which is also subject of investigation in this paper. The optimization goal is to further improve the shape of the chosen X0-specimen so as to achieve a preferably homogeneous stress triaxiality distribution in the area the specimen is predicted to fail. With this at hand, a parameter identification for special damage models that are to capture different damage mechanisms will be easier and more accurate.

For design optimization, gradient-based methods are utilized. The gradient information is determined by means of a variational approach proposed in Barthold (2002); Barthold (2008); Barthold et al. (2016); Barthold and Stein (1996). It is based on an enhanced kinematic concept that offers a rigorous separation of structural and physical quantities and provides exact gradient information efficiently with moderate effort. Additionally, it allows simultaneous computation of stress states and sensitivities within a finite element framework. Thus, the proposed approach is distinct from other variational techniques but can be linked to well-known methods such as the Material Derivative Approach, cf. e.g. Zolésio (1981) and the Domain Parametrization Approach, cf. e.g. Phelan and Haber (1989). A detailed comparison of the mentioned approaches is beyond the scope of this paper.

In contrast to a purely elastic mechanical response, here the deformation history deserves special attention in structural as well as sensitivity analysis. Numerous researchers investigated the theoretical and computation details which are linked to pioneering contributions such as Michaleris et al. (1994). In our research, we follow the preparatory work published in Wiechmann (2000). The shape optimization of specimens has already been treated in literature, and a few contributions are highlighted. The reduction of stress concentration in the transition zone between a straight-sided zone and a wider end zone has been discussed in, e.g. Klemensø et al. (2007). The relation of the geometric shape of the specimens on the parameter identification problem from a mathematical point of view has been investigated in Bauer et al. (2016). Shape optimization of biaxially loaded specimens is more sophisticated as outlined in, e.g. Makris et al. (2010) and Etling and Herzog (2018). This paper pursues our research on design sensitivity analysis applied to elastoplasticity outlined in Liedmann and Barthold (2018a) and represents an extended version of contribution Liedmann and Barthold (2019).

In Section 2, the main optimization problem is stated in terms of the objective function, design variables, and nonlinear constraints. After some preliminaries given in Section 3 about the viewpoint of kinematics and the notation used in this paper, the model equations for the chosen J₂ elastoplastic model are summarized in Section 4 and the solution strategy is described. Section 4.1 focuses on the mechanical response analysis. Based on this, the variational sensitivity analysis is consequently performed, and the concept to compute the mechanical response and stress sensitivities are outlined in Section 4.2. Section 5 is concerned with numerics. The method of $\bar {F}$ is briefly explained in Section 5.1. Furthermore, the concept of the so-called design velocity matrix is described in Section 5.2, followed by some remarks about the computational effort in Section 5.3. In Section 6, the method is applied to two optimization problems, a parameter identification in Section 6.1 and the shape optimization of the chosen X0-specimen in Section 6.2. Section 7 summarizes the paper and draws a conclusion on the findings. Discretization using $\bar {F}$ finite elements and the corresponding explicit discrete operators for structural as well as sensitivity analysis are given in the Appendix.

2 Optimization tasks

In this section, the two optimization problems subject in this paper are introduced, namely the main shape optimization and a preceding parameter fitting of the underlying material model used for the simulations to experimental data. Both problems are solved utilizing gradient based optimization methods and analytical gradients are provided. The computation of the gradients and the underlying structural analysis problem is presented in Section 4.

2.1 Shape optimization

The main goal in this paper is the shape optimization of the X0-specimen presented, e.g. in Gerke et al. (2019). The initial geometry is illustrated in Fig. 1. In the notched area, where the specimen is predicted to fail, the stress triaxiality distribution shall be preferably homogeneous. Therefore, the approach is to maximize the stress triaxiality intensity in that area by varying the geometric design variables defining the shape of the notches. The design vector $\boldsymbol {p} = \left [ R_{1} \ R_{2} \ R_{3} \ d \right ]$ contains the design variables indicated in Fig. 2.

Fig. 1

Initial geometry of X0-specimen

Fig. 2

Design variables in the notched area

Here, l_cs and t_cs denote the length and thickness of the central cross-section, respectively. The objective function is chosen to be

$$ \max_{p} \quad J(\boldsymbol{p}) = || \boldsymbol{\eta}(\boldsymbol{p}) ||, $$

(2)

where the vector η stores the stress triaxiality values at the cross-section. The problem is constrained in terms of the equality constraint of the form

$$ {c}^{\mathsf{eq}} = t_{\mathsf{cs}} \cdot l_{\mathsf{cs}} - \bar{A}_{\mathsf{cs}} = 0, $$

(3)

which forces a constant cross-section area $\bar {A}_{\mathsf {cs}}$ of the specimen. An additional condition is that the penetration depth d has always to be smaller than the radius of the notch in thickness direction R₃

$$ {c}^{\mathsf{in}} = d - R_3 < 0 $$

(4)

to prevent sharp edges in the geometry and the finite element mesh. The results of the shape optimization are presented in Section 6.2.

2.2 Parameter fitting

As in future investigations, the results shall be validated in real biaxial experiments; in a first step, the elastic and plastic material parameters of the underlying constitutive model are fitted to approximate the behavior of the real material (AlCuMgSi) as good as possible. Therefore, the first inverse problem to solve is a parameter identification, in which the global load-displacement curve of a uniaxial tension test is fitted to experimental data. The objective function is of least squares type and takes the form

$$ \min_{m} \quad J(\mathbf{m}) = || \boldsymbol{f}^\mathsf{R}(\mathbf{m}) - \boldsymbol{F}^\mathsf{R} ||, $$

(5)

where the vector F^R represents the reaction forces measured in the experiment and the vector f^R(m) stores the simulated reaction forces depending on the design variables $\boldsymbol {m} = [E \nu \sigma _{0} \sigma _{\infty } H \beta ]$.

Details on the meaning of the constitutive parameters are presented in Section 4.1. Results of the curve fitting procedure are presented in Section 6.1.

3 Preliminaries

3.1 Tensor notation and operators

Symbolic tensor notation is used in this paper. Vectors are bold face italic, second-order tensors are upright bold face, and sans serif and fourth-order tensors are blackboard bold characters

$$ \begin{array}{@{}rcl@{}} \boldsymbol{v} &=& v_i \boldsymbol{e}_i, \quad \boldsymbol{\mathsf{A}} =A_{ij} \boldsymbol{e}_i\otimes\boldsymbol{e}_j, \\ \mathbb{T} &=& T_{ijkl} \boldsymbol{e}_i\otimes\boldsymbol{e}_j\otimes\boldsymbol{e}_k\otimes\boldsymbol{e}_l. \end{array} $$

(6)

Single contraction is either denoted as dot product for vectors or omitted for higher order tensors, double contraction is denoted as colon, e.g.

$$ \begin{array}{@{}rcl@{}} \boldsymbol{\mathsf{A}} \boldsymbol{v} &=& \boldsymbol{\mathsf{A}}\cdot\boldsymbol{v} = A_{ij}v_j \boldsymbol{e}_i, \\ \mathbb{T} : \boldsymbol{\mathsf{A}} &=& T_{ijkl}A_{kl} \boldsymbol{e}_i\otimes\boldsymbol{e}_j. \end{array} $$

(7)

Special transpositions of fourth-order tensors are denoted as superscripts, e.g.

$$ \mathbb{T}^{\overset{23}{\mathsf{T}}} = T_{ikjl} \boldsymbol{e}_i\otimes\boldsymbol{e}_j\otimes\boldsymbol{e}_k\otimes\boldsymbol{e}_l, $$

(8)

which represents the transposition of the second and third bases of a fourth-order tensor.

In some equations, the special product $\overset {21}{\ast }$ is used, which denotes a single contraction of the second basis of a fourth-order tensor with the first basis of a second-order tensor, i.e.,

$$ \mathbb{T} \overset{21}{*} \boldsymbol{\mathsf{A}} = T_{ijlm}A_{jk} \boldsymbol{e}_i\otimes\boldsymbol{e}_k\otimes\boldsymbol{e}_l\otimes\boldsymbol{e}_m. $$

(9)

The special fourth-order tensors $\mathbb {I}_{\mathsf {sym}}$ and $\mathbb {I}_{{\mathsf {dev}}}$ denote the symmetric and deviatoric projection tensors with the coefficients

$$ \begin{array}{@{}rcl@{}} {I}_{ijkl}^{\mathsf{s}} &=& \frac{1}{2} (\delta_{ik} \delta_{jl}+\delta_{il} \delta_{jk}),\\ {I}_{ijkl}^{\mathsf{d}} &=& {I}_{ijkl}^\mathsf{s}-\frac{1}{3} \delta_{ij} \delta_{kl}. \end{array} $$

(10)

3.2 Enhanced kinematics

In the context of structural optimization, the material body is not considered with a fixed reference configuration. Thus, within a general continuum mechanical framework, it is convenient to work with an enhanced kinematic concept, see Fig. 3. Motivated by an improved viewpoint on the material body using arguments from differential geometry, cf., Noll (1974); Bertram (1989); Barthold (2002), the main idea is the rigorous separation of physical and geometrical quantities. By the introduction of a fixed parameter space ${\mathcal{B}}$ with Cartesian basis Z_i and local coordinates Θ, the classical deformation mapping φ, defined as

$$ \boldsymbol{\varphi} : (\boldsymbol{X},t) \mapsto \boldsymbol{x}(\boldsymbol{X},t) $$

(11)

can be decomposed into two independent mappings, i.e., the design-dependent geometry mapping κ(Θ,s) and the time-dependent motion mapping μ(Θ,t)

$$ \boldsymbol{\kappa} : (\boldsymbol{\varTheta},s) \mapsto \boldsymbol{X}(\boldsymbol{\varTheta},s) \text{ and } \boldsymbol{\mu} : (\boldsymbol{\varTheta},t) \mapsto \boldsymbol{x}(\boldsymbol{\varTheta},t). $$

(12)

Fig. 3

Enhanced kinematics

With the corresponding tangent mappings

$$ \boldsymbol{\mathsf{K}} = {\nabla}_{\varTheta} \boldsymbol{\kappa} = \boldsymbol{G}_{i} \otimes \boldsymbol{Z}^{i} \text{ and } \boldsymbol{\mathsf{M}} = {\nabla}_{\varTheta}\boldsymbol{\mu} = \boldsymbol{g}_i \otimes \boldsymbol{Z}^i, $$

(13)

not only the deformation mapping but also the deformation gradient can be decomposed and written as

$$ \boldsymbol{\varphi} = \boldsymbol{\mu} \circ \boldsymbol{\kappa}^{\mathsf{-1}} \text{ and } \boldsymbol{\mathsf{F}} = \boldsymbol{\mathsf{M}} \boldsymbol{\mathsf{K}}^{\mathsf{-1}} = \boldsymbol{\mathsf{g}}_i \otimes \boldsymbol{\mathsf{G}}^{i}. $$

(14)

Consequently, for volume mappings, we obtain

$$ \int\limits_{\mathcal{M}}\mathrm{d}v = \int\limits_{\mathcal{K}} J \mathrm{d}V = \int\limits_{\mathcal{B}} J J_K \mathrm{d}V_{{\varTheta}} = \int\limits_{\mathcal{B}}J_M \mathrm{d}V_{{\varTheta}}, $$

(15)

with the determinants $J = \det {\boldsymbol {\mathsf {F}}}$, $J_{K} = \det \boldsymbol {\mathsf {K}}$, and $J_{M} = \det \boldsymbol {\mathsf {M}}$. The main advantage of this enhanced viewpoint of kinematics is that implicit dependencies do not arise until the definition of global equilibrium. The parameter s is used as a scalar design variable parameterizing the material body in the reference configuration $\mathcal {K} = \mathcal {K}(s)$ and the referential points X = X(s).

Note, due to the choice of the elastoplastic material model, the plastic intermediate configuration $\mathcal {P}$ is defined up to a rigid body rotation, cf., Simo and Hughes (2000). Thus, the tangent mappings F_e and F_p stay undefined during the computational procedure.

3.3 Variations and derivatives

The total variation of any quantity (∙)(u;s;h_n) is given by the sum of the partial variations

$$ \begin{array}{@{}rcl@{}} \delta (\bullet)(\boldsymbol{u};\boldsymbol{s};\boldsymbol{h}_n) &=& \delta_u (\bullet)(\boldsymbol{u};\hat{\boldsymbol{s}};\hat{\boldsymbol{h}}_n) \\ &&+ \delta_s (\bullet)(\hat{\boldsymbol{u}};\boldsymbol{s};\hat{\boldsymbol{h}}_n)\\ &&+ \delta_{h_n} (\bullet)(\hat{\boldsymbol{u}};\hat{\boldsymbol{s}};\boldsymbol{h}_n), \end{array} $$

(16)

where, e.g., the partial variation with respect to u with fixed $\hat {\boldsymbol {s}}$ and $\hat {\boldsymbol {h}}_{n}$ is defined as

$$ \delta_{u}(\bullet) = \frac{\mathrm{d}}{\mathrm{d}\epsilon}(\bullet)(\boldsymbol{u}+\epsilon \delta\boldsymbol{u};\hat{\boldsymbol{s}};\hat{\boldsymbol{h}}_{n})\big|_{\epsilon=0}. $$

(17)

This notation is used throughout the whole paper.

4 Response and sensitivity analysis

Following, the material model, solution strategies for computing the structural response and response sensitivities are outlined. The chosen mechanical model based on a multiplicative split of the deformation gradient F = F_eF_p can be found in the relevant literature, e.g. Elguedj and Hughes (2014); Simo (1988a, b) and captures large deformations and strains.

4.1 Mechanical response

The nonlinear weak form of equilibrium $R : \mathcal {V} \to \mathbb {R}$, commonly called residual, in the reference configuration $\mathcal {K}$ reads

$$ \begin{array}{@{}rcl@{}} R(\boldsymbol{u};\boldsymbol{v}) &=& \int\limits_\mathcal{K}\boldsymbol{\mathsf{P}}^{\mathsf{K}}(\boldsymbol{u}) : \text{Grad}\boldsymbol{v} \mathrm{d}V \\ &&- \int\limits_\mathcal{K} \boldsymbol{b}_0 \cdot \boldsymbol{v} \mathrm{d}V - \int\limits_{\partial\mathcal{K}} \boldsymbol{t}_0 \cdot \boldsymbol{v} \mathrm{d}A = 0, \end{array} $$

(18)

where Gradv is the gradient of the test function $\boldsymbol {v} \in \mathcal {V}$ and $\boldsymbol {u} \in \mathcal {V}$ is the displacement field. Here, $\mathcal {V}$ denotes the usual Sobolev space of all admissible displacements. The chosen strain energy function given by

$$ \begin{array}{@{}rcl@{}} W &=& U(J) + W(\boldsymbol{\widetilde{\mathsf{C}}}, \boldsymbol{\widetilde{\mathsf{C}}}_{\mathsf{p}}^{\mathsf{-1}})\\ &=& \frac{K}{2} \left[\frac{1}{2} (J^2-1)-\ln J\right] + \frac{G}{2} \left[\widetilde{\boldsymbol{\mathsf{C}}} : \widetilde{\boldsymbol{\mathsf{C}}}_{\mathsf{p}}^{\mathsf{-1}} - 3\right] \end{array} $$

(19)

is an extension to the compressible range of the Neo–Hookean model with the advantage of poly-convexity and is suitable for large deformations and strains, cf. Simo (1988a). Here, K and G denote the elastic compression and shear moduli, respectively. Stresses can be determined via a return-mapping algorithm, cf. Elguedj and Hughes (2014); Simo (1988b). Further, $\widetilde {\boldsymbol {\mathsf {C}}}_{\mathsf {p}}^{\mathsf {-1}}$ denotes the isochoric part of the inverse plastic contribution of the left Cauchy–Green deformation tensor, that is, the pull-back of the isochoric elastic left Cauchy–Green deformation tensor $\widetilde {\boldsymbol {\mathsf {b}}}_{\mathsf {e}} = \widetilde {\boldsymbol {\mathsf {F}}}_{\mathsf {e}} \widetilde {\boldsymbol {\mathsf {F}}}_{\mathsf {e}}^{\mathsf {T}}$ using the isochoric deformation gradient $\widetilde {\boldsymbol {\mathsf {F}}}$

$$ \widetilde{\boldsymbol{\mathsf{C}}}_{\mathsf{p}}^{\mathsf{-1}} = \widetilde{\boldsymbol{\mathsf{F}}}^{\mathsf{-1}} \widetilde{\boldsymbol{\mathsf{b}}}_{\mathsf{e}} \widetilde{\boldsymbol{\mathsf{F}}}^{\mathsf{-T}}. $$

(20)

The von Mises–type yielding condition only depends on the deviatoric part of the stresses, and thus is pressure independent

$$ f = || \boldsymbol{\mathsf{s}} || - \sqrt{\frac{2}{3}} k(\alpha) \leq 0, $$

(21)

with s = devτ denoting the deviatoric part of the Kirchhoff stress tensor τ. Hardening during plastic flow is defined by the nonlinear function that depends on the isochoric hardening variable α

$$ k(\alpha) = \sigma_{0} + \sigma_{\infty} [1-\exp(-\beta \alpha)] + H \alpha, $$

(22)

where the material parameters σ₀ and $\sigma _{\infty }$ denote the initial and limit yield stress, respectively. H is the linear hardening modulus and β is a dimensionless parameter. Note that although the deformation process is assumed to be quasi-static, a pseudo-time has to be introduced in order to solve the evolution of the internal hardening variables.

The set of history variables that has to be saved in each load step and integration point is $\boldsymbol {h} = \{\widetilde {\boldsymbol {\mathsf {C}}}_{\mathsf {p}}^{\mathsf {-1}}, \alpha \}$. Evolution of the internal history variables is described by the flow rule

$$ \dot{\widetilde{\boldsymbol{\mathsf{C}}}}^{\mathsf{-1}}_{\mathsf{p}} = -\frac{2}{3} \dot{\gamma}\text{ tr }\boldsymbol{\mathsf{b}}_\mathsf{e} \boldsymbol{\mathsf{F}}^\mathsf{-1} \boldsymbol{\mathsf{n}} \boldsymbol{\mathsf{F}}^\mathsf{-T} $$

(23)

and the hardening law

$$ \dot{\alpha} = \sqrt{\frac{2}{3}} \dot{\gamma}. $$

(24)

Due to the nonlinear hardening function, within an implicit time integration, the increment of the plastic multiplier Δγ has to be solved within a local Newton–Raphson procedure, see Algorithm 2. The solution has to be a Kuhn–Tucker point defined by the conditions

$$ {\varDelta}\gamma =0, f(\boldsymbol{\mathsf{\tau}},\alpha) \leq 0, {\varDelta}\gamma f = 0. $$

(25)

With the increment of the plastic multiplier at hand, the update of the deviatoric stress tensor can be determined as

$$ \boldsymbol{\mathsf{s}} = \boldsymbol{\mathsf{s}}^\mathsf{tr}-2 \bar{\mu} {\varDelta}\gamma \boldsymbol{\mathsf{n}}, \boldsymbol{\mathsf{n}} = \frac{\boldsymbol{\mathsf{s}}^\mathsf{tr}}{|| \boldsymbol{\mathsf{s}}^\mathsf{tr} ||}, $$

(26)

with $\bar {\mu } = \frac {1}{3} G \text { tr } \widetilde {\boldsymbol {\mathsf {b}}}_{\mathsf {e}}^{\mathsf {tr}}$, $\boldsymbol {\mathsf {s}}^{\mathsf {tr}} = G \text {dev}\widetilde {\boldsymbol {\mathsf {b}}}_{\mathsf {e}}^{\mathsf {tr}}$ and the flow direction n. The superscript tr denotes the trial state in the elastic predictor step of the return-mapping procedure, see Algorithm 1. Due to the assumption of incompressible plastic flow, the hydrostatic pressure can easily be computed to

$$ p = p^{\mathsf{tr}} = U^{\prime}(J) = \frac{K}{2 J} \left( J^{2}-1\right). $$

(27)

The total Kirchhoff stress tensor reads

$$ \boldsymbol{\mathsf{\tau}} = J p \boldsymbol{\mathsf{I}} + \boldsymbol{\mathsf{s}} $$

(28)

and the first Piola–Kirchhoff stress tensor is given by

$$ \boldsymbol{\mathsf{P}}^{\mathsf{K}} = \boldsymbol{\mathsf{\tau}} \boldsymbol{\mathsf{F}}^{\mathsf{-T}}. $$

(29)

Consistent linearization of the stress tensor has been derived in Simo (1988a, b), and the numerical solution strategy is outlined. An excellent summary of the solution algorithm can also be found in Elguedj and Hughes (2014) in the context of isogeometric analysis. However, using the $\bar {F}$ method and in view of sensitivity analysis, it is convenient to use the algorithmic material tangent $\mathbb {A} = \frac {\partial \boldsymbol {\mathsf {P}}^{\mathsf {K}}}{\partial \boldsymbol {\mathsf {F}}}$, which can be split into a volumetric and deviatoric part

$$ \mathbb{A} = \frac{\partial\boldsymbol{\mathsf{P}}^{\mathsf{K}}_{\text{vol}}}{\partial\boldsymbol{\mathsf{F}}} + \frac{\partial\boldsymbol{\mathsf{P}}^{\mathsf{K}}_{\text{dev}}}{\partial\boldsymbol{\mathsf{F}}} = \mathbb{A}_{\text{vol}} + \mathbb{A}_{\text{dev}}. $$

(30)

The explicit forms of the volumetric and deviatoric parts read

$$ \mathbb{A}_{\text{vol}} = K J^2 \boldsymbol{\mathsf{F}}^\mathsf{-T} \otimes \boldsymbol{\mathsf{F}}^\mathsf{-T} - J p \left( \boldsymbol{\mathsf{F}}^\mathsf{-T} \otimes \boldsymbol{\mathsf{F}}^\mathsf{-T}\right)^{\overset{24}{\mathsf{T}}} $$

(31)

and

$$ \mathbb{A}_{\text{dev}} = \mathbb{S}_{\text{dev}} {~}^{\overset{21}{*}}\boldsymbol{\mathsf{F}}^\mathsf{-T} - \boldsymbol{\mathsf{s}} \left( \boldsymbol{\mathsf{F}}^\mathsf{-T} \otimes \boldsymbol{\mathsf{F}}^\mathsf{-T}\right)^{\overset{24}{\mathsf{T}}}, $$

(32)

respectively. The fourth order tensor $\mathbb {S}_{\text {dev}}$ is given by

$$ \begin{array}{@{}rcl@{}} \mathbb{S}_{\text{dev}} &=& \frac{\partial\boldsymbol{\mathsf{s}}}{\partial\boldsymbol{\mathsf{F}}} \\ &=& G \left[\beta_0 \mathbb{I}_{\text{dev}} - \frac{2}{3} {\varDelta}\gamma \beta_1 \boldsymbol{\mathsf{n}} \otimes \boldsymbol{\mathsf{I}} + \beta_2 \boldsymbol{\mathsf{n}} \otimes \boldsymbol{\mathsf{n}}\right] : \mathbb{B} \\ & =& \mathbb{H} : \mathbb{B}, \end{array} $$

(33)

with the factors

$$ \beta_0 = 1-2 {\varDelta}\gamma \frac{\overline{\mu}}{|| \boldsymbol{\mathsf{s}}^\mathsf{tr} ||}, \beta_1 = 1+2 \frac{\overline{\mu}}{f^{\prime}}, \beta_2 = \beta_1-\beta_0 $$

(34)

and the identities

$$ \mathbb{B} := \frac{\partial\widetilde{\boldsymbol{\mathsf{b}}}_{\mathsf{e}}^{\mathsf{tr}}}{\partial\boldsymbol{\mathsf{F}}} = \left[\mathbb{F} {~}^{\overset{21}{*}} \left( \widetilde{\boldsymbol{\mathsf{C}}}_{\mathsf{p},n}^{\mathsf{-1}} \widetilde{\boldsymbol{\mathsf{F}}}^{\mathsf{T}}\right) + \widetilde{\boldsymbol{\mathsf{F}}} \widetilde{\boldsymbol{\mathsf{C}}}_{\mathsf{p},n}^{\mathsf{-1}} \mathbb{F}^{\overset{12}{\mathsf{T}}}\right], $$

(35)

$$ \mathbb{F} := \frac{\partial\widetilde{\boldsymbol{\mathsf{F}}}}{\partial\boldsymbol{\mathsf{F}}} = J^{-\frac{1}{3}} \left[\left( \boldsymbol{\mathsf{I}} \otimes \boldsymbol{\mathsf{I}}\right)^{\overset{23}{\mathsf{T}}} - \frac{1}{3} \boldsymbol{\mathsf{F}}\otimes\boldsymbol{\mathsf{F}}^\mathsf{-T}\right]. $$

(36)

The whole computational procedure is summarized in Algorithm 1 in pseudo code format.

4.2 Response sensitivity

Once a global equilibrium point has been computed, that is, a solution of (18), response sensitivities are obtained variationally at continuous level. As (18) has to hold for any design change δs, its total variation has to vanish, cf. Barthold (2002); Liedmann and Barthold (2018a); Wiechmann (2000), and we can write

$$ \begin{array}{@{}rcl@{}} \delta R &=& \delta_u R + \delta_s R + \delta_{h_n} R = 0 \\ &=& k(\boldsymbol{v},\delta\boldsymbol{u}) + p(\boldsymbol{v},\delta\boldsymbol{s}) + h(\boldsymbol{v},\delta\boldsymbol{h}_n), \end{array} $$

(37)

where the two bilinear forms $k : \mathcal {V} \times \mathcal {V} \to \mathbb {R}$ and $p : \mathcal {V} \times \mathcal {S} \to \mathbb {R}$ represent the partial variations of the weak equilibrium condition w.r.t. displacements and design, respectively. The third bilinear form $h : \mathcal {V} \times \mathcal {G} \to \mathbb {R}$ corresponds to the deformation history that has to be captured and considered for all total variations, cf. Liedmann and Barthold (2018a, b) and Wiechmann (2000). The spaces $\mathcal {S}$ and $\mathcal {G}$ denote Sobolev spaces with all admissible design and history functions, respectively. We refer to, e.g. Haug et al. (1986); Sokołowski and Zolésio (1992); Kim and Choi (2004); Choi and Kim (2004) for details on the function spaces and the mentioned notation using linear and bilinear operators.

4.3 Matrix form of sensitivity relations

From the matrix form of (37)

$$ \begin{array}{@{}rcl@{}} \delta\mathbf{R} &=& \left[\frac{\partial\mathbf{R}}{\partial\mathbf{u}}\right] \delta\mathbf{u} + \left[\frac{\partial\mathbf{R}}{\partial\mathbf{s}}\right] \delta\mathbf{s} + \left[\frac{\partial\mathbf{R}}{\partial\mathbf{h}_n}\right] \delta\mathbf{h}_n = \mathbf{0} \\ &=& \mathbf{K} \delta\mathbf{u} + \mathbf{P} \delta\mathbf{s} + \mathbf{H} \delta\mathbf{h}_n, \end{array} $$

(38)

where $\mathbf {R} \in \mathbb {R}^{\mathsf {ndof}}$, the sensitivity matrix $\mathbf {S} \in \mathbb {R}^{\mathsf {ndof \times ndv}}$ can be derived, that is, the total derivative of the structural response u w.r.t. design changes s

$$ \begin{array}{@{}rcl@{}} \delta\mathbf{u} &=& -\mathbf{K}^\mathsf{-1} \left[\mathbf{P} \delta\mathbf{s} + \mathbf{H} \delta\mathbf{h}_n\right]\\ &=& -\mathbf{K}^\mathsf{-1} \left[\mathbf{P}+\mathbf{H} \mathbf{Z}_n\right] \delta\mathbf{s} = \mathbf{S} \delta\mathbf{s}. \end{array} $$

(39)

Here, the matrix $\mathbf {Z}_{n} \in \mathbb {R}^{\mathsf {nhv \times ndv}}$ is the total design derivative of the history variables

$$ \delta\mathbf{h}_{n} = \mathbf{Z}_{n} \delta \mathbf{s}, $$

(40)

see Section 4.4 for details.

$\mathbf {K} \in \mathbb {R}^{\mathsf {ndof \times ndof}}$ denotes the tangent stiffness matrix, $\mathbf {P} \in \mathbb {R}^{\mathsf {ndof \times ndv}}$ denotes the so-called pseudo load matrix and $\mathbf {H} \in \mathbb {R}^{\mathsf {nhv \times nhv}}$ is the history sensitivity matrix. Here, ndof is the number of total degrees of freedom, ndv is the number of design variables, and nhv is the number of history variables. Note that assuming h₀ = 0 and δh₀ = 0 at time t = t₀, from (39) we obtain

$$ \mathbf{S}_{0} = -\mathbf{K}^{\mathsf{-1}} \mathbf{P},\quad \mathbf{Z}_{0} = \mathbf{0} $$

(41)

for the first pseudo-time increment. The sensitivity part corresponding to the deformation history has to be evaluated at the end of each time step that causes plastic yielding and saved for the subsequent step. Note that all values are evaluated in the current pseudo-time step t = t_n+ 1 except those that have been saved in the prior step with subscript n.

4.3.1 Pseudo load matrix

Referring to (38), the pseudo load matrix is the partial derivative of the residual R w.r.t. design s. A pull-back to parameter space yields

$$ \begin{array}{@{}rcl@{}} \delta_{s} R &=& \int\limits_{\partial\mathcal{B}} \delta_s\boldsymbol{\mathsf{P}}^{\mathsf{K}} : \text{Grad}\boldsymbol{v} J_K \mathrm{d}V_{{\varTheta}}\\ && + \int\limits_{\partial\mathcal{B}} \boldsymbol{\mathsf{P}}^{\mathsf{K}} : \delta_s\text{Grad}\boldsymbol{v} J_K \mathrm{d}V_{{\varTheta}} \\ &&+ \int\limits_{\partial\mathcal{B}} \boldsymbol{\mathsf{P}}^{\mathsf{K}} : \text{Grad}\boldsymbol{v} \delta_s J_K \mathrm{d}V_{{\varTheta}}, \end{array} $$

(42)

where $J_{K} = \det \boldsymbol {\mathsf {K}}$. At this point, we have to consider the choice of design variables as this will affect functional dependencies.

Geometric design

For geometric design, the design vector is chosen as the vector of geometry points s := X. With the partial variations of the gradient of the test function

$$ \delta_X \text{Grad}\boldsymbol{v} = -\text{Grad}\boldsymbol{v} \text{Grad}\delta\boldsymbol{X}, $$

(43)

the determinant of the geometry gradient

$$ \delta_X J_K = J_K \text{Div}\delta\boldsymbol{X}, $$

(44)

the first Piola–Kirchhoff stress tensor

$$ \delta_X \boldsymbol{\mathsf{P}}^{\mathsf{K}} = \frac{\partial\boldsymbol{\mathsf{P}}^{\mathsf{K}}}{\partial\boldsymbol{\mathsf{F}}} : \delta_X\boldsymbol{\mathsf{F}} = \mathbb{A} : \delta_X\boldsymbol{\mathsf{F}} $$

(45)

and the deformation gradient

$$ \delta_X\boldsymbol{\mathsf{F}} = -\text{Grad}\boldsymbol{u} \text{Grad}\delta\boldsymbol{X}, $$

(46)

cf. Barthold (2002), Materna (2009), and Kijanski (2018), we can write

$$ \begin{array}{@{}rcl@{}} \delta_{X} R &= & - \int\limits_{\partial\mathcal{K}} \left[\left( \text{Grad}\boldsymbol{u} \text{Grad}\delta\boldsymbol{X}\right) : \mathbb{A}\right] : \text{Grad}\boldsymbol{v} \mathrm{d}V \\ && - \int\limits_{\partial\mathcal{K}} \boldsymbol{\mathsf{P}}^{\mathsf{K}} : \text{Grad}\boldsymbol{v} \text{Grad}\delta\boldsymbol{X} \mathrm{d}V \\ && + \int\limits_{\partial\mathcal{K}} \boldsymbol{\mathsf{P}}^{\mathsf{K}} : \text{Grad}\boldsymbol{v} \text{Div}\delta\boldsymbol{X} \mathrm{d}V, \end{array} $$

(47)

where we have used the identity dV = J_KdV_Θ.

Constitutive parameters

For the sensitivity regarding material parameters, the design vector is chosen as the vector containing the material parameters s := m. As only the stress tensor depends on material parameters, we find

$$ \delta_m\text{Grad}\boldsymbol{v} = \boldsymbol{0}, \delta_m J_K = 0, $$

(48)

but

$$ \delta_m\boldsymbol{\mathsf{P}}^{\mathsf{K}} = \frac{\partial\boldsymbol{\mathsf{P}}^{\mathsf{K}}}{\partial\boldsymbol{m}} \delta\boldsymbol{m} = \underline{\boldsymbol{\mathsf{M}}} \delta\boldsymbol{m}, $$

(49)

where $\underline {\boldsymbol {\mathsf {M}}}$ is the third-order tensor connecting the first Piola–Kirchhoff stress tensor with the material parameters. Thus, we can write

$$ \delta_{m} R = \int\limits_{\mathcal{K}} \text{Grad}\boldsymbol{v}^\mathsf{T} \underline{\boldsymbol{\mathsf{M}}} \delta\boldsymbol{m} \mathrm{d}V. $$

(50)

For reasons of brevity, without detailed derivation, for the vector of material parameters $\boldsymbol {m} = [E \nu \sigma _{0} \sigma _{\infty } H \beta ]$, the explicit form of the tensor $\underline {\boldsymbol {\mathsf {M}}}$ reads

$$ \begin{array}{@{}rcl@{}} \underline{\boldsymbol{\mathsf{M}}} &=&\frac{1}{2}\left( J^{2} - 1\right) \boldsymbol{\mathsf{F}}^{\mathsf{-T}} \otimes \frac{\partial K}{\partial\boldsymbol{m}} \\ && + \boldsymbol{\mathsf{H}}_{\mathsf{m}}\boldsymbol{\mathsf{F}}^{\mathsf{-T}} \otimes \frac{\partial G}{\partial\boldsymbol{m}} \\ && - 2 \sqrt{\frac{2}{3}} \frac{\bar{\mu}}{f^{\prime}}\boldsymbol{\mathsf{n}} \boldsymbol{\mathsf{F}}^{\mathsf{-T}} \otimes \frac{\partial k}{\partial\boldsymbol{m}}, \end{array} $$

(51)

Here, the second-order tensor H_m reads

$$ \boldsymbol{\mathsf{H}}_{\mathsf{m}} = \left( \beta_0\mathbb{I}_{\text{sym}} + \beta_2\boldsymbol{\mathsf{n}} \otimes \boldsymbol{\mathsf{n}}\right):\text{dev}\widetilde{\boldsymbol{\mathsf{b}}}_{\mathsf{e}}^{\mathsf{tr}} - 2{\varDelta}\gamma\frac{\bar{\mu}}{G}\beta_1\boldsymbol{\mathsf{n}}. $$

(52)

The partial derivatives of the bulk modulus K, the shear modulus G, and the nonlinear hardening function k are given by

$$ \frac{\partial K}{\partial\boldsymbol{m}} = [a_1 a_2 0 0 0 0], $$

(53)

with $a_{1} = \frac {1}{3-6 \nu }, \quad a_{2} = \frac {2 E}{3 (1-2 \nu )^{2}}$,

$$ \frac{\partial G}{\partial\boldsymbol{m}} = [b_1 b_2 0 0 0 0], $$

(54)

with $b_{1} = \frac {1}{2 (1+\nu )}, \quad b_{2} = -\frac {E}{2 (1+\nu )^{2}}$ and

$$ \frac{\partial k}{\partial\boldsymbol{m}} = [ 0 0 1 c_1 \alpha c_2], $$

(55)

with $c_{1} = 1-\exp (-\beta \alpha ), \quad c_{2} = \sigma _{\infty } \alpha \exp (-\beta \alpha )$, respectively.

4.3.2 History sensitivity matrix

The history sensitivity matrix is the partial derivative of the residual R w.r.t. the internal history variables of the prior pseudo-time step h_n and reads

$$ \delta_{h_{n}} R = \int\limits_{\partial\mathcal{B}} \delta_{h_n}\boldsymbol{\mathsf{P}}^{\mathsf{K}} : \text{Grad}\boldsymbol{\eta} J_K \mathrm{d}V_{{\varTheta}}, $$

(56)

as only the stress tensor P^K depends on the history variables. The partial variation of the first Piola–Kirchhoff stress tensor w.r.t. the history variables reads

$$ \delta_{h_n}\boldsymbol{\mathsf{P}}^{\mathsf{K}} = \frac{\partial\boldsymbol{\mathsf{P}}^\mathsf{K}}{\partial\boldsymbol{h}_n} \mathbf{Z}_n \delta\boldsymbol{s}. $$

(57)

The partial derivatives of the stress tensor w.r.t. the history variables $\boldsymbol {h}_{n} = \{\widetilde {\boldsymbol {\mathsf {C}}}_{\mathsf {p},n}^{\mathsf {-1}}, \alpha _{n}\}$ are given by

$$ \frac{\partial\boldsymbol{\mathsf{P}}^{\mathsf{K}}}{\partial\widetilde{\boldsymbol{\mathsf{C}}}_{\mathsf{p},n}^{\mathsf{-1}}} = \left( \mathbb{H} : \mathbb{B}_C\right) {~}^{\overset{21}{*}} \boldsymbol{\mathsf{F}}^{\mathsf{-T}}, $$

(58)

with

$$ \mathbb{B}_{C} = \boldsymbol{\widetilde{\mathsf{F}}} \mathbb{I}_s {~}^{\overset{21}{*}} \boldsymbol{\widetilde{\mathsf{F}}}^\mathsf{T} $$

(59)

and the tensor $\mathbb {H}$ from (32) for the internal variable $\widetilde {\boldsymbol {\mathsf {C}}}_{\mathsf {p}}^{\mathsf {-1}}$ and

$$ \frac{\partial\boldsymbol{\mathsf{P}}^\mathsf{K}}{\partial\alpha_n} = -2 \sqrt{\frac{2}{3}} \frac{\bar{\mu}}{f^{\prime}} k^{\prime} \boldsymbol{\mathsf{n}} \boldsymbol{\mathsf{F}}^\mathsf{-T} $$

(60)

for the internal variable α, with the first derivatives of (22) and (21)

$$ k^{\prime} = \sigma_\infty d \exp(-\beta \alpha)+H $$

(61)

and

$$ f^{\prime} = -2 \bar{\mu} \left[1+\frac{k^{\prime}}{3 \bar{\mu}}\right]. $$

(62)

4.4 Update of history variations

As (39) indicates, it is essential to update the total variations of the internal history variables at the end of each pseudo-time step that has caused plastic deformation. The values of the total derivatives are saved into the matrix Z, denoted as total design derivative of the history variables, cf. (40). The general update formula of the variations of the internal history variables reads

$$ \delta\boldsymbol{h} = \frac{\partial\boldsymbol{h}}{\partial\boldsymbol{u}} \delta\boldsymbol{u} + \frac{\partial\boldsymbol{h}}{\partial\boldsymbol{s}} \delta\boldsymbol{s} + \frac{\partial\boldsymbol{h}}{\partial\boldsymbol{h}_n} \delta\boldsymbol{h}_n. $$

(63)

Again, as for the pseudo load matrices, we have to distinguish between the choice of design variables.

Geometric design

For geometric design, (63) takes the form

$$ \delta\boldsymbol{h} = \frac{\partial\boldsymbol{h}}{\partial\boldsymbol{\mathsf{F}}} : \left( \delta_u\boldsymbol{\mathsf{F}} + \delta_X\boldsymbol{\mathsf{F}}\right) + \frac{\partial\boldsymbol{h}}{\partial\boldsymbol{h}_n} \mathbf{Z}_n \delta\boldsymbol{X} = \mathbf{Z} \delta\boldsymbol{X}, $$

(64)

with δ_uF = Gradδu. The yet unknown partial variations $\frac {\partial \boldsymbol {h}}{\partial \boldsymbol {\mathsf {F}}}$ and $\frac {\partial \boldsymbol {h}}{\partial \boldsymbol {h}_{n}}$ can be resolved to

$$ \begin{array}{@{}rcl@{}} \frac{\partial\widetilde{\boldsymbol{\mathsf{C}}}_{\mathsf{p}}^{\mathsf{-1}}}{\partial\boldsymbol{\mathsf{F}}} &=& -\left[\left( \widetilde{\boldsymbol{\mathsf{F}}}^\mathsf{-1} \otimes \widetilde{\boldsymbol{\mathsf{F}}}^\mathsf{-T}\right)^{\underset{\mathsf{T}}{23}} : \mathbb{F}\right] {~}^{\overset{21}{*}} \left( \widetilde{\boldsymbol{\mathsf{b}}}_{\mathsf{e}} \widetilde{\boldsymbol{\mathsf{F}}}^\mathsf{-T}\right) \\ &&+ \widetilde{\boldsymbol{\mathsf{F}}}^\mathsf{-1}\left[\left( \frac{1}{G} \mathbb{H} + \frac{1}{3} \boldsymbol{\mathsf{I}} \otimes \boldsymbol{\mathsf{I}}\right) : \mathbb{B}\right]{~}^{\overset{21}{*}} \widetilde{\boldsymbol{\mathsf{F}}}^\mathsf{-T} \\ &&- \widetilde{\boldsymbol{\mathsf{F}}}^\mathsf{-1} \widetilde{\boldsymbol{\mathsf{b}}}_{\mathsf{e}} \left[\left( \widetilde{\boldsymbol{\mathsf{F}}}^\mathsf{-T} \otimes \widetilde{\boldsymbol{\mathsf{F}}}^\mathsf{-T}\right)^{\overset{24}{\mathsf{T}}} : \mathbb{F}\right] \end{array} $$

(65)

for the internal variable $\widetilde {\boldsymbol {\mathsf {C}}}_{\mathsf {p}}^{\mathsf {-1}}$ w.r.t. the deformation tensor F and

$$ \frac{\partial\alpha}{\partial\boldsymbol{\mathsf{F}}} = \sqrt{\frac{2}{3}} \frac{G}{f^{\prime}} \left( \frac{2}{3} {\varDelta}\gamma \boldsymbol{\mathsf{I}} - \boldsymbol{\mathsf{n}}\right) : \mathbb{B}, $$

(66)

for the internal variable α.

Constitutive parameters

For material parameters as design variables, (63) takes the form

$$ \begin{array}{@{}rcl@{}} \delta\boldsymbol{h} &=& \frac{\partial\boldsymbol{h}}{\partial\boldsymbol{\mathsf{F}}} : \delta_u\boldsymbol{\mathsf{F}} + \left( \frac{\partial\boldsymbol{h}}{\partial\boldsymbol{m}} + \frac{\partial\boldsymbol{h}}{\partial\boldsymbol{h}_n}\mathbf{Z}_n\right) \delta\boldsymbol{m} \\ &=& \mathbf{Z} \delta\boldsymbol{m}. \end{array} $$

(67)

The partial derivatives of the internal history variables w.r.t. the material parameters $\frac {\partial \boldsymbol {h}}{\partial \boldsymbol {m}}$ can be identified to

$$ \frac{\partial\widetilde{\boldsymbol{\mathsf{C}}}_{\mathsf{p}}^{\mathsf{-1}}}{\partial\boldsymbol{m}} = \frac{1}{G}\mathbb{B}_{C}^{-1} : \left[\left( \boldsymbol{\mathsf{H}}_{\mathsf{m}}-\frac{\boldsymbol{\mathsf{s}}}{G}\right) \otimes\frac{\partial G}{\partial\boldsymbol{m}} -2 \sqrt{\frac{2}{3}} \frac{\bar{\mu}}{f^{\prime}} \boldsymbol{\mathsf{n}} \otimes \frac{\partial k}{\partial\boldsymbol{m}}\right] $$

(68)

for the internal variable $\widetilde {\boldsymbol {\mathsf {C}}}_{\mathsf {p}}^{\mathsf {-1}}$, where the second-order tensor H_m has been used, cf. (52). For the internal variable α, we find

$$ \begin{array}{@{}rcl@{}} \frac{\partial\alpha}{\partial\boldsymbol{m}} &=& \sqrt{\frac{2}{3}} \frac{1}{f^{\prime}} \left( \frac{2}{G} {\varDelta}\gamma \bar{\mu}-\boldsymbol{\mathsf{n}} : \text{dev}\boldsymbol{\widetilde{\mathsf{b}}}_{\mathsf{e}}^{\mathsf{tr}}\right) \frac{\partial G}{\partial\boldsymbol{m}} \\ &&+ \frac{2}{3 f^{\prime}} \frac{\partial k}{\partial\boldsymbol{m}}. \end{array} $$

(69)

The partial derivatives of the shear modulus and the nonlinear hardening function are already known, see (53)–(55).

In both cases, the partial derivatives of the internal variables w.r.t. their counterparts from the previous pseudo-time step are needed and read

$$ \frac{\partial\widetilde{\boldsymbol{\mathsf{C}}}_{\mathsf{p}}^{\mathsf{-1}}}{\partial\widetilde{\boldsymbol{\mathsf{C}}}_{\mathsf{p},n}^{\mathsf{-1}}} = \mathbb{B}_C^\mathsf{-1} : \left[\frac{1}{G} \mathbb{H} + \frac{1}{3} \boldsymbol{\mathsf{I}}\otimes\boldsymbol{\mathsf{I}} \right] : \mathbb{B}_C $$

(70)

and

$$ \frac{\partial\widetilde{\boldsymbol{\mathsf{C}}}_{\mathsf{p}}^{\mathsf{-1}}}{\partial\alpha_n} = -2 \sqrt{\frac{2}{3}} \frac{k^{\prime}}{f^{\prime}} \mathbb{B}_C^\mathsf{-1} : \boldsymbol{\mathsf{n}} $$

(71)

for the internal variable $\widetilde {\boldsymbol {\mathsf {C}}}_{\mathsf {p}}^{\mathsf {-1}}$. The fourth-order tensor $\mathbb {B}_{C}^{\mathsf {-1}}$ is given by

$$ \mathbb{B}_C^\mathsf{-1} = \widetilde{\boldsymbol{\mathsf{F}}}^\mathsf{-1}\mathbb{I}_s {~}^{\overset{21}{*}} \widetilde{\boldsymbol{\mathsf{F}}}^\mathsf{-T}. $$

(72)

For the internal variable α we find

$$ \frac{\partial\alpha}{\partial\widetilde{\boldsymbol{\mathsf{C}}}_{\mathsf{p},n}^{\mathsf{-1}}} = \sqrt{\frac{2}{3}} \frac{G}{f^{\prime}} \left[\frac{2}{3} {\Delta}\gamma \boldsymbol{\mathsf{I}} - \boldsymbol{\mathsf{n}}\right] : \mathbb{B}_C, $$

(73)

with the fourth-order tensor $\mathbb {B}_{C}$ from (59), and

$$ \frac{\partial\alpha}{\partial\alpha_n} = 1 + \frac{2}{3} \frac{k^{\prime}}{f^{\prime}}, $$

(74)

with the derivatives given in (61) and (62).

4.5 Stress sensitivity

The gradient of the stress triaxiality as defined in (1) can be computed following the same concepts, which leads to

$$ \delta\eta = \frac{\partial\eta}{\partial\sigma_{\mathsf{m}}} \delta\sigma_{\mathsf{m}} + \frac{\partial\eta}{\partial\sigma_{\mathsf{eq}}} \delta\sigma_{\mathsf{eq}}, $$

(75)

where the total variations of the mean and von Mises equivalent stress read

$$ \delta\sigma_\mathsf{m} = \frac{\partial\sigma_\mathsf{m}}{\partial\boldsymbol{\mathsf{\sigma}}} : \delta\boldsymbol{\mathsf{\sigma}}, \delta\sigma_{\mathsf{eq}} = \frac{\partial\sigma_{\mathsf{eq}}}{\partial\boldsymbol{\mathsf{\sigma}}} : \delta\boldsymbol{\mathsf{\sigma}}. $$

(76)

Thus, it is essential to calculate the total variation of the Cauchy stress tensor given by

$$ \delta\boldsymbol{\mathsf{\sigma}} = \delta_{u}\boldsymbol{\mathsf{\sigma}} + \delta_{s}\boldsymbol{\mathsf{\sigma}} + \delta_{h_{n}}\boldsymbol{\mathsf{\sigma}}. $$

(77)

This variation follows the same principles as the variation of the first Piola–Kirchhoff stress tensor, cf. Section 4.2, and is straight forward

$$ \begin{array}{@{}rcl@{}} \delta\boldsymbol{\mathsf{\sigma}} & =& \frac{\partial\boldsymbol{\mathsf{\sigma}}}{\partial\boldsymbol{\mathsf{F}}} : \delta_u\boldsymbol{\mathsf{F}} + \frac{\partial\boldsymbol{\mathsf{\sigma}}}{\partial\boldsymbol{s}} \delta\boldsymbol{\mathsf{s}} + \frac{\partial\boldsymbol{\mathsf{\sigma}}}{\partial\boldsymbol{h}_n} \delta\boldsymbol{h}_n \\ &=& \frac{\partial\boldsymbol{\mathsf{\sigma}}}{\partial\boldsymbol{\mathsf{F}}} : \frac{\partial\boldsymbol{\mathsf{F}}}{\partial\boldsymbol{u}} \delta\boldsymbol{u} + \frac{\partial\boldsymbol{\mathsf{\sigma}}}{\partial\boldsymbol{s}} \delta\boldsymbol{\mathsf{s}} + \frac{\partial\boldsymbol{\mathsf{\sigma}}}{\partial\boldsymbol{h}_n} \mathbf{Z}_n \delta\boldsymbol{s}\\ &=& \left[\frac{\partial\boldsymbol{\mathsf{\sigma}}}{\partial\boldsymbol{\mathsf{F}}} : \frac{\partial\boldsymbol{\mathsf{F}}}{\partial\boldsymbol{u}} \mathbf{S} + \frac{\partial\boldsymbol{\mathsf{\sigma}}}{\partial\boldsymbol{s}} + \frac{\partial\boldsymbol{\mathsf{\sigma}}}{\partial\boldsymbol{h}_n} \mathbf{Z}_n\right] \delta\boldsymbol{s}. \end{array} $$

(78)

5 Numerical treatment

5.1 $\bar {F}$ finite elements

Due to the assumption of incompressible plastic flow, volumetric locking occurs in the solution using pure displacement finite elements. Therefore, as proposed in de Souza Neto et al. (2009) and Elguedj and Hughes (2014), an $\bar {F}$ formulation is used. This formulation has advantages on implementation side as any constitutive law can be used. The $\bar {F}$ deformation gradient reads

$$ \bar{\boldsymbol{\mathsf{F}}} = \left( \frac{J_{{0}}}{J}\right)^{\frac{1}{3}} \boldsymbol{\mathsf{F}}=\kappa\boldsymbol{\mathsf{F}}, $$

(79)

where $J_{{0}} = \det \boldsymbol {\mathsf {F}}_{{0}}$ is the determinant of the deformation gradient evaluated at the element centroid. It is important to mention that all equations in Section 4 are evaluated at $\boldsymbol {\mathsf {F}} = \bar {\boldsymbol {\mathsf {F}}}$. Therefore, for exact linearization, the derivatives of $\bar {\boldsymbol {\mathsf {F}}}$ w.r.t. F and F₀ are needed. The total variation of $\bar {\boldsymbol {\mathsf {F}}}$ reads

$$ \begin{array}{@{}rcl@{}} \delta\bar{\boldsymbol{\mathsf{F}}} &=& \frac{\partial\bar{\boldsymbol{\mathsf{F}}}}{\partial\boldsymbol{\mathsf{F}}} : \delta\boldsymbol{\mathsf{F}} + \frac{\partial\bar{\boldsymbol{\mathsf{F}}}}{\partial\boldsymbol{\mathsf{F}}_{{0}}} : \delta\boldsymbol{\mathsf{F}}_{{0}} \\ &=&\frac{\kappa}{3} \left[3 \left( \boldsymbol{\mathsf{I}}\otimes\boldsymbol{\mathsf{I}}\right)^{\overset{23}{\mathsf{T}}}-\boldsymbol{\mathsf{F}}\otimes\boldsymbol{\mathsf{F}}^\mathsf{-T}\right] : \delta\boldsymbol{\mathsf{F}} \\ &&+\frac{\kappa}{3} \left[\boldsymbol{\mathsf{F}}\otimes\boldsymbol{\mathsf{F}}_{{0}}^\mathsf{-T}\right] : \delta\boldsymbol{\mathsf{F}}_{{0}} \\ &=& \mathbb{D} : \delta\boldsymbol{\mathsf{F}} + \mathbb{D}_{{0}} :\delta\boldsymbol{\mathsf{F}}_{{0}}. \end{array} $$

(80)

5.2 Geometry modeling and derivatives

In the case of shape optimization within a finite element framework, the sensitivity relations of the quantities of interest have to be computed regarding the nodal coordinates of the mesh. To ensure smooth boundaries and to reduce the number of design variables, it is often convenient to choose mesh controlling geometry parameters as design variables. The so-called design velocity matrix, defined as

$$ \mathbf{Q} := \frac{\partial\mathbf{X}}{\partial\mathbf{p}} \quad \in \mathbb{R}^{\mathsf{dof \times ndv}}, $$

(81)

connects the mesh coordinates X and chosen geometry controlling parameters p and can be computed analytically as far as the geometrical dependency in known, see e.g. Gerzen (2014) and Kijanski (2018). However, for the case that third-party programs are used for mesh generation, the design velocity matrix can be computed numerically using a finite difference scheme, e.g.

$$ \mathbf{Q} = \frac{\mathbf{X}(\mathbf{p}+\epsilon) - \mathbf{X}(\mathbf{p})}{\epsilon}, $$

(82)

using forward finite differences. Here, 𝜖 denotes a small perturbation number or step size.

Note, this approach can be numerically expensive but has the advantage of simplicity, especially for complex geometries. Using numerical differentiation of the design velocity matrix, the whole method of sensitivity analysis can not be called pure analytical. However, the numerical cost is still tiny compared to pure numerical approaches. Nevertheless, the authors recommend analytical derivation of the design velocity matrix whenever possible.

Listing 1

Matlab Code Ex. 1: numerical design velocity matrix

In Section 6.2, the geometry of the X0-specimen and the mesh have been built up using gmsh4.4.0. Thus, the design velocity matrix has been computed numerically. Matlab Code Ex. (1) summarizes the numerical treatment using MatlabR2018a running in a bash shell on Linux. Here, the file X0.geo stores the geometry of the X0-specimen that is parameterized by three radii (R1,R2,R3) and a penetration depth of the notch in thickness direction (d). The variable e represents the step size 𝜖 of (82). It is chosen here to e = sqrt(eps), which is the default forward finite difference step size in the Matlab solver fmincon. This reduces the number of error sources when analytical gradients are checked using fmincon.

Note that in Matlab eps =5⁻⁵² represents the distance of 1.0 to the next larger double-precision number.

Listing 2

Matlab Code Ex. 2: global history field

Listing 3

Matlab Code Ex. 3: call to lsqnonlin

5.3 Computational effort

For each integration point in each element, the vector of history variables $\mathbf {h} \in \mathbb {R}^{\mathsf {1 \times nhv}}$ has to be saved for the subsequent pseudo-time step for structural analysis. Additionally, for sensitivity analysis, the partial variations of the history variables $\mathbf {Z} \in \mathbb {R}^{\mathsf {nhv \times ndv}}$ w.r.t. design $\mathbf {s} \in \mathbb {R}^{\mathsf {ndv \times 1}}$ have to be saved. The structure of the global history field is realized in Matlab using a heterogeneous cell array of the form

$$ \mathcal{H} := \left\{\begin{array}{ccccccccc} \mathbf{h}_1^1 & {\dots} & \mathbf{h}^{\mathsf{ngp}}_1 && \vline && \mathbf{Z}_1^1 & {\dots} & \boldsymbol{Z}^{\mathsf{ngp}}_1 \\ {\vdots} & & {\vdots} && \vline && {\vdots} & & {\vdots} \\ \mathbf{h}_{\mathsf{nel}}^1 & {\dots} & \mathbf{h}^{\mathsf{ngp}}_{\mathsf{nel}} && \vline && \mathbf{Z}_{\mathsf{nel}}^1 & {\dots} & \boldsymbol{Z}^{\mathsf{ngp}}_{\mathsf{nel}} \end{array},\right\} $$

(83)

which can directly be used for parallel element evaluation. Matlab Code Ex. (2) shows how the history field is constructed. In total, in each pseudo-time step, we need to save $\mathsf {nh} = \mathsf {nel \times ngp \times nhv \times \left (1 + ndv\right )}$ scalar values. If the number of design variables ndv increases, the total amount of values that have to be stored nh explodes. Thus, the authors recommend parametrization of the mesh and utilizing the design velocity matrix on element level as explained in Appendix (A).

6 Optimization results

6.1 Parameter identification

In a first step, the method mentioned above is used to fit the material parameters of the constitutive model to the material behavior by means of a standard tensile test, cf. Fig. 4. The discrete model is pictured in Fig. 5. In this case, the design variables are chosen as the material parameters $\mathbf {s} := \mathbf {m} = [E \nu \sigma _0 \sigma _{\infty } H \beta ]$. The initial values of the material parameters are taken from Gerke et al. (2017) and summarized below:

$$ \begin{array}{@{}rcl@{}} E &=& {69}~\text{GPa}, \quad \nu = {0.3}, \quad \sigma_0 = {320}~\text{MPa},\\ \sigma_\infty &=& {100}~\text{MPa}, \quad H = {650}~\text{MPa}, \quad \beta = {30}. \end{array} $$

The geometry is discretized with 1539 hexahedral $\bar {F}$ elements. In total, the mechanical problem counts 4872 degrees of freedom. Due to symmetry, only one eighth of the geometry is modeled. Thus, the measured reaction force is divided by 4 for the comparison. Symmetric boundary conditions are applied and a maximum displacement of 10 mm is prescribed within 200 linear steps, that is, the point where the maximum reaction force has been measured. The time discretization implies that the number of data points that are compared with the simulation is n_dp = 200. The optimization problem is to minimize the error between the load-displacement curves of the experiment and the simulation, where F^R denotes the measured reaction forces and f^R(m) the simulated reaction forces depending on the set of material parameters. We can define the optimization problem as

$$ \begin{array}{@{}rcl@{}} &&\min_{m} \quad J(\mathbf{m}) = || \boldsymbol{f}^\mathsf{R}(\mathbf{m}) - \boldsymbol{F}^\mathsf{R} || \\ &&\text{s.t. } \quad m^i_l \leq m^i \leq m^i_u. \end{array} $$

(84)

The solution of the problem is obtained utilizing the Matlab solver lsqnonlin with the options shown in Matlab Code Ex. (3). Here, also the chosen values of upper and lower bounds of the design variables are declared. Note that the function is a separate Matlab routine that computes the value of the objective function J in (84). In Fig. 6 the history of the objective function during the curve fitting process in plotted. 22 iterations are needed to finally reach the stopping criterion, that is, the relative change of the objective compared to the prior iteration is less than tol = 1e − 6. Load-displacement curves are shown for different iterations in Fig. 7. The parameters are identified to

$$ \begin{array}{@{}rcl@{}} E &=& 61.92~\text{GPa},\quad~~ \nu = 0.31,\qquad\quad ~ \sigma_{0} = 313.3~\text{MPa},\\ \sigma_{\infty} &=& 119.2~\text{MPa},\quad H = 548.7~\text{MPa},\quad \beta =10.72. \end{array} $$

With this at hand, the next optimization task, that is the shape optimization of a chosen specimen, can be conducted.

Fig. 4

Tensile test

Fig. 5

FE model

Fig. 6

Objective history

Fig. 7

Load-Displacement curves (AlCuMgSi)

6.2 Specimen shape optimization

The chosen, so-called X0-specimen, is shown in Fig. 1. Due to symmetry, only one eighth of the geometry has to be discretized as pictured in Fig. 8. 5928 hexahedral $\bar {F}$ elements are used and the mechanical problem counts 23,535 degrees of freedom. Within 200 linear steps, the maximum displacement of 0.75 mm is applied to all sides of the specimen.

Fig. 8

FE model

Listing 4

Matlab Code Ex. 4: call to fmincon

Note, this displacement is much more than in the experiment. For more realistic results, the maximum displacement has to be adapted. The results in this paper can bee seen as a proof of concept.

As the shape of the specimen is to be optimized, geometric parameters are chosen as design variables s := p = [R₁R₂R₃d], namely

$$ \begin{array}{@{}rcl@{}} &&R_1 : \text{inner radius}, R_3 : \text{radius in thickness direction},\\ &&R_2 : \text{outer radius}, d : \text{penetration depth}. \end{array} $$

Aim of the shape optimization is to achieve a preferably homogeneous distribution of stress triaxiality in the cross-section. For this, the length of the vector of stress triaxiality in the cross-section is maximized. Additionally, during the optimization the cross-section has to stay constant, which is defined by the equality constraint

$$ c^{\mathsf{eq}} = t_{\mathsf{cs}} \cdot l_{\mathsf{cs}} - {12}~\text{mm}^{2} = 0, $$

(85)

cf. Fig. 2. To prevent the appearance of sharp corners at the edges of the notch in thickness direction, the penetration depth should not be larger than the radius in thickness direction. This is captured by the inequality constraint

$$ c^{\mathsf{in}} = d - R_{3} \leq 0. $$

(86)

Consequently, the optimization problem is defined as

$$ \begin{array}{@{}rcl@{}} \max_{X} \quad J(\mathbf{p}) &=& ||\boldsymbol{\eta}(\mathbf{p})||\\ \text{s.t. } \quad c^\mathsf{eq} &=& 0,\\ c^\mathsf{in} &\leq& 0, \\ {p}_{l}^{i} &\leq& p^i \leq {p}_{u}^{i}. \end{array} $$

(87)

The problem is solved with the Matlab solver fmincon. The options used as well as the bounds of the design variables are shown in Matlab Code Ex. (4). Here, the functions and are separate Matlab functions that compute the values of the objective function J and the geometry constraints in (85)–(87), respectively.

After 28 iterations, the solver found a solution. In Fig. 9, the stress triaxiality at the cross-section of the specimen is illustrated. Due to the equality constraint, the maximum stress triaxiality does not change significantly, but the distribution over the cross-section is much more homogeneous in the optimized case than it was for the initial shape. The values of the objective function during the optimization process are displayed in Fig. 10. Initial and optimized values of the design parameters are summarized in Table 1. The values of the constraints at converged state are

$$ c^{\mathsf{eq}} = 3.2925\mathrm{e}-10 \text{and} c^\mathsf{in} = -~{8.9378\mathrm{e}-4}. $$

Fig. 9

Objective history

Fig. 10

Triaxiality distribution at specimen cross section

Table 1 Initial and optimized design parameters

7 Summary and conclusions

In this work, an efficient method for the computation of design sensitivities for structures with elastoplastic material behavior has been presented. The method is based on a variational approach that allows simultaneous computation of structural response and response sensitivities. The obtained gradient information can be used to feed gradient-based optimization algorithms. The method has been applied to two different optimization problems. First, a parameter identification has been performed to fit the model parameters to the real material (AlCuMgSi). The curve fitting procedure resulted in an excellent agreement with the experimentally measured load-displacement curve. In a second example, the shape of the chosen X0-specimen has been optimized with the aim to achieve a preferably homogeneous stress triaxiality distribution. The stress triaxiality characterizes the stress state of the parts of the structure and can be used to classify different damage mechanisms. In the presented example, the length of the vector of stress triaxiality in the relevant part of the specimen has been maximized. The optimization resulted in a slightly different shape, which shows a much more homogeneous stress triaxiality distribution in the relevant part of the specimen. Future investigations will address the validation of the results in terms of manufacturing the optimized shapes and test them. Further, different specimen types and also different loading scenarios will be analyzed.

Appendix

1.1 A.1 Discrete gradient and divergence operators

For an 8-node hexahedral finite element, the discrete gradient operator at each FE node i reads

$$ \begin{array}{@{}rcl@{}} \mathbf{G}_i &= \left[\begin{array}{ccc} N_{i,1} & 0 & 0 \\ 0 &N_{i,1} & 0 \\ 0 & 0 & N_{i,1} \\ N_{i,2} & 0 & 0 \\ 0 &N_{i,2} & 0 \\ 0 & 0 & N_{i,2} \\ N_{i,3} & 0 & 0 \\ 0 &N_{i,3} & 0 \\ 0 & 0 & N_{i,3} \end{array}\right]\in \mathbb{R}^{\mathsf{9 \times 3}} \end{array} $$

(88)

with the interpolation functions N_i. Discrete matrix forms of non-symmetric fourth-order tensors follow the convention

$$ \mathbf{A} = \left[\begin{array}{cccccc} A_{1111} & {\dots} & A_{1131} & A_{1112} & {\dots} & A_{1133} \\ {\vdots} &{\ddots} & & & & \vdots\\ A_{3111} & & & & & A_{3133} \\ A_{1211} & & & & & A_{1233} \\ {\vdots} & & & &{\ddots} & {\vdots} \\ A_{3311} & {\dots} & A_{3331} & A_{3312} & {\dots} & A_{3333} \end{array}\right] \in \mathbb{R}^{\mathsf{9 \times 9}}, $$

(89)

as well as discrete vector forms of non-symmetric second-order tensors

$$ \boldsymbol{P} = [P_{11} P_{21} P_{31} P_{12} {\dots} P_{33}]^\mathsf{T} \in \mathbb{R}^{\mathsf{9 \times 1}}. $$

(90)

In the case of the shape sensitivity analysis, the partial variation of the deformation gradient w.r.t. geometry points is needed, cf. (46). Its discrete version takes the form

$$ {\mathbf{G}}_{i}^{X} = \mathbf{G}_i \text{Grad}\mathbf{u}^{h}, $$

(91)

with

$$ \begin{array}{@{}rcl@{}} \text{Grad}\mathbf{u}^h = \sum\limits_{i}^{\mathsf{n}}\mathbf{u}_i \mathbf{L}_i^\mathsf{T}, \mathbf{u}_i &= \left[\begin{array}{c} u^i_1 \\ u^i_2 \\ u^i_3 \end{array}\right] \mathbf{L}_i = \left[\begin{array}{c} N_{i,1} \\ N_{i,2} \\ N_{i,3} \end{array}\right] . \end{array} $$

The discrete divergence operator reads

$$ \text{Div}\boldsymbol{u}^h = \sum\limits_{i}^{\mathsf{n}} \mathbf{L}^T_i \mathbf{u}_i. $$

With the matrix forms D and D₀ of $\mathbb {D}$ and $\mathbb {D}_{{0}}$ from (80) using the convention stated in (89), we can define the two gradient operators for the $\bar {F}$ element

$$ \begin{array}{@{}rcl@{}} \bar{\mathbf{G}}_i &=&\mathbf{D} \mathbf{G}_i + \mathbf{D}_{{0}} \mathbf{G}_{{0},i}, \\ \bar{\mathbf{G}}_{i}^{X} &=&\mathbf{D} \mathbf{G}_{i}^{X} + \mathbf{D}_{{0}} \mathbf{G}_{{0},i}^{X}. \end{array} $$

(92)

1.2 A.2 Discrete equations

Using the discrete quantities and conventions stated in Appendix A.1, this section summarizes the discrete element vectors and matrices.

Element residual vector

$$ \mathbf{R}^e_{i} = \int\limits_{\mathcal{K}} \mathbf{G}_i^\mathsf{T} \mathbf{P}^\mathsf{K} \mathrm{d}V $$

(93)

Element stiffness matrix

$$ \mathbf{K}^e_{ij} = \int\limits_{\mathcal{K}} \mathbf{G}_i^\mathsf{T} \mathbf{A} \bar{\mathbf{G}}_j \mathrm{d}V $$

(94)

Element geometric pseudo load matrix

$$ \begin{array}{@{}rcl@{}} \mathbf{P}^e_{X,ij}& = & -\int\limits_{\mathcal{K}} \mathbf{G}_i^\mathsf{T} \mathbf{A} \bar{\mathbf{G}}_{j}^{X} \mathrm{d}V \\ &&-\int\limits_{\mathcal{K}} \boldsymbol{\mathsf{P}}^{\mathsf{K}} \left( \mathbf{L}_j \mathbf{L}_i^\mathsf{T} - \mathbf{L}_i \mathbf{L}_j^\mathsf{T}\right) \mathrm{d}V \end{array} $$

(95)

Note, by the choice of mesh controlling geometric parameters p as design variables, cf. Section 5.2, the design velocity matrix can be used on element level. This results in

$$ \mathbf{P}^e_{p,ik} = \mathbf{P}^e_{X,ij} \mathbf{Q}^e_{jk}, $$

(96)

which drastically decreases the sizes of the pseudo load, global sensitivity, and history sensitivity matrix, and thus saves a huge amount of memory. Additionally, in (64), we have to consider

$$ \delta\boldsymbol{\mathsf{F}} = \delta_u\boldsymbol{\mathsf{F}} + \delta_p\boldsymbol{\mathsf{F}} = \left[\frac{\partial\boldsymbol{\mathsf{F}}}{\partial\boldsymbol{\mathsf{u}}} \mathbf{S}^e + \frac{\partial\boldsymbol{\mathsf{F}}}{\partial\boldsymbol{X}} \mathbf{Q}^e\right] \delta\mathbf{p}. $$

(97)

Element material pseudo load matrix

$$ \mathbf{P}^e_{m,ij} = \int\limits_{\mathcal{K}} \mathbf{G}_i^\mathsf{T} \mathbf{M} \mathrm{d}V $$

(98)

Element history sensitivity matrix

$$ {\mathbf{H}}_{ij}^{e} = \int\limits_{\mathcal{K}} \mathbf{G}_i^\mathsf{T} \frac{\partial\boldsymbol{\mathsf{P}}^{\mathsf{K}}}{\partial\boldsymbol{h}_n} \mathbf{Z}^e_{n,j} \mathrm{d}V $$

(99)