1 Introduction

Ductile sheet metals are of outstanding relevance in many engineering applications and consequently, the need to improve the lightweight design, to decrease the energy consumption and to increase the cost efficiency is evident and thus, the utilization factor of these materials has to be augmented. This leads to the demand to characterize the material behavior within the inelastic domain properly and to avoid early localization of irreversible strains as well as damage and fracture within structural components. The connected material failure processes strongly depend on the stress state: In tension dominated regions nucleation, growth and coalescence of micro-voids prevail, whereas shear stress states facilitate micro-shear-cracks. To analyze this behavior, experiments with carefully designed specimens covering a wide range of relevant stress states have to be performed.

Standard specimens are frequently designed to generate one pre-defined stress state within the region of interest whereas only one loading condition can be applied. Recent developments in this field indicate the request to develop new geometries for in-plane sheet metal testing in a representative range of loading conditions. In this context, the stress is frequently characterized by the stress intensity, the stress triaxiality and the Lode parameter which can be related to the material degradation processes. For instance, differently notched tensile specimens [1, 10, 16, 19] generate increasing stress triaxialities with decreasing notch radius. A similar effect can be induced by introducing a central hole in flat specimens [39]. Rather shear dominated stress states can be achieved with single [10, 19, 25] and double [39] connection one-dimensional specimens. The butterfly [35] specimen can be used in a special testing device producing a big variety of stress states from tension to shear domination. In this connection it has to be mentioned that the previously discussed specimens produce significant stress gradients within the region of interest, which might be unfavorable for material characterization.

Cruciform specimens facilitate the possibility to study different stress states with one geometry. Standard geometries have been applied successfully to determine the yield surface of ductile metals [28]. Unfortunately, these standard geometries tend to localize with ongoing deformation and unexpected failure of the specimen occurs. In this context, a new experimental program with biaxially loaded specimens has been proposed [9, 21], where stress states have been analyzed by corresponding numerical simulations to predict damage and fracture modes. In addition, [20] proposed several new geometries and first experimental results with the X0-specimen have been presented [13, 22, 23]. This specimen (Fig. 1) is characterized by a central opening and crosswise arranged notches producing four connecting parts. Consequently, the localization is pre-defined by the notched regions and the corresponding stress state can be determined by corresponding numerical simulations. However, to investigate the material behavior and the damage mechanisms in a more satisfactory way in a first cut it is preferable to reach on one hand high stress triaxialities (micro-voids) and on the other hand low stress triaxialities (micro-shear-cracks).

For this purpose, mathematical optimization techniques can be utilized to modify the shape of specimens in order to gain these preferred stress states in an area of interest. In this context, gradient based methods have proven to be efficient and stable, cf. [43]. Gradient information is determined variationally at continuous level as described in [2,3,4,5]. An enhanced viewpoint of kinematics that offers a rigorous separation of geometrical and physical quantities is convenient. Here, the classical continuum mechanical deformation mapping is split into a geometry mapping and a motion mapping. Therefore, the approach differs from other variational approaches. However, it can be linked to common methods like the Domain Parametrization Approach, cf. e.g. [38], or the Material Derivative Approach, cf. e.g. [46]. With any of the mentioned variational approaches, exact gradient information can be provided efficiently with moderate effort and additionally, it allows the computation of stress states and stress sensitivities simultaneously within a finite element framework.

In this work, the aforementioned X0-specimen is subject of investigation for shape optimization. The aim is to gain a distinct (high or low) stress triaxiality with a sufficiently homogeneous distribution in the notched regions of the specimen where damage and failure are expected to occur. Obviously, during the biaxial experiment, the specimen undergoes large plastic deformations before significant damage evolution takes place. The onset of damage evolution can be defined by a stress state damage criterion taking into account the different material degradation processes. Within this paper, special focus is given to the stress state at this point, i.e. in this first approach damage is not considered within the material model. Thus, the chosen material model, proposed in e.g. [42], captures large deformations and strains. Here, the deformation history deserves special attention in structural analysis as well as in sensitivity analysis. In the past, many researchers investigated into theoretical and computational details in the context of sensitivity formulations in elastoplasticity, cf. e.g. one of the pioneering contributions [34]. However, formulations in this work are based on the preparatory contribution [44].

Specimen shape optimization has already been tackled in literature. For instance in [27] the reduction of the stress concentration in the transition zone between straight-sided and wider-end zones has been subject of discussion. The relation between the geometric shape of specimens and parameter identification problems have been discussed in [7]. As outlined in e.g. [18] or [32], shape optimization of biaxially loaded specimens is more complex and challenging. In this paper the research on design sensitivity analysis applied to elastoplasticity as outlined in [29] is amplified.

Fig. 1
figure 1

X0-specimen

Full size image

The paper is structured as follows. Hints on the notation and the aforementioned underlying viewpoint of kinematics are given in Sect. 2. In continuation, Sect. 3 indicates the constitutive equations and the material used for the experiments is characterized. Section 4 summarizes the numerical model for structural and sensitivity analysis. The optimization problem is stated in Sect. 5 and the solution is depicted and discussed. In Sect. 6 the experimental setup is shown and results of the biaxial tests are illustrated. Furthermore, the results of the initial and shape optimized specimens are censoriously compared and discussed. Section 7 summarizes and draws a conclusion on the findings. Finally, an outlook to further investigations is given.

2 Preliminaries

2.1 Notation and operators

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

$$\begin{aligned} \varvec{v}= & {} v_i\,\varvec{e}_i, \quad {{\mathbf {\mathsf{{A}}}}} =A_{ij}\,\varvec{e}_i\otimes \varvec{e}_j, \nonumber \\ \mathbb {T}= & {} T_{ijkl}\,\varvec{e}_i\otimes \varvec{e}_j\otimes \varvec{e}_k\otimes \varvec{e}_l. \end{aligned}$$
(1)

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{aligned} {{\mathbf {\mathsf{{A}}}}}\,\varvec{v} = {{\mathbf {\mathsf{{A}}}}}\cdot \varvec{v} = A_{ij}v_j\,\varvec{e}_i,\ \mathbb {T} : {{\mathbf {\mathsf{{A}}}}} = T_{ijkl}A_{kl}\,\varvec{e}_i\otimes \varvec{e}_j. \end{aligned}$$
(2)

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

$$\begin{aligned} \mathbb {T}^{{\mathop {\textsf {T}}\limits ^{23}}} = T_{ikjl}\,\varvec{e}_i\otimes \varvec{e}_j\otimes \varvec{e}_k\otimes \varvec{e}_l. \end{aligned}$$
(3)

In some equations, the special product \(\,{{\,\mathrm{{\mathop {*}\limits ^{21}}}\,}}\,\) denoting a single contraction of the second basis of a fourth order tensor with the first basis of a second order tensor is used, i.e.

$$\begin{aligned} \mathbb {T} {{\,\mathrm{{\mathop {*}\limits ^{21}}}\,}}{{\mathbf {\mathsf{{A}}}}} = T_{ijlm}A_{jk}\,\varvec{e}_i\otimes \varvec{e}_k\otimes \varvec{e}_l\otimes \varvec{e}_m. \end{aligned}$$
(4)

The fourth order tensors \(\mathbb {I}_{\textsf {s}}\) and \(\mathbb {I}_{\textsf {d}}\) denote the symmetric and deviatoric projection tensors with Cartesian coefficients

$$\begin{aligned} I^\textsf {s}_{ijkl}= & {} \frac{1}{2}\,(\delta _{ik}\,\delta _{jl}+\delta _{il}\,\delta _{jk}),\nonumber \\ I^\textsf {d}_{ijkl}= & {} I^\textsf {s}_{ijkl}-\frac{1}{3}\,\delta _{ij}\,\delta _{kl}. \end{aligned}$$
(5)

Discrete matrix forms are denoted by upright bold face characters, e.g.

$$\begin{aligned} \mathbf{\mathrm{\mathbf{A}} } = [A_{ij}] \in \mathbb {R}^{\textsf {n} \times \textsf {m}}. \end{aligned}$$
(6)

2.2 Enhanced kinematics

In shape optimization, the material body is not considered with a fixed reference configuration. Thus, it is convenient to work with an enhanced viewpoint of kinematics within a general continuum mechanical framework. Motivated by arguments from differential geometry, cf. [2, 8, 37], the main idea is to rigorously separate physical and geometrical quantities. Introducing a fixed parameter space \({\mathcal{B}}\) with Cartesian basis \(\varvec{Z}_i\) and local coordinates \(\varvec{\Theta }\), the classical continuum mechanical deformation mapping \(\varvec{\varphi }\) that maps the referential coordinates \(\varvec{X}\) to the actual coordinates \(\varvec{x}\) at time t, defined as

$$\begin{aligned} \varvec{\varphi } : (\varvec{X},t) \mapsto \varvec{x}(\varvec{X},t), \end{aligned}$$
(7)

can be decomposed into two independent mappings, i.e. the design dependent geometry mapping \(\varvec{\kappa }(\varvec{\Theta },s)\) and the time dependent motion mapping \(\varvec{\mu }(\varvec{\Theta },t)\)

$$\begin{aligned} \varvec{\kappa } : (\varvec{\Theta },s) \mapsto \varvec{X}(\varvec{\Theta },s) \text { and } \varvec{\mu } : (\varvec{\Theta },t) \mapsto \varvec{x}(\varvec{\Theta },t), \end{aligned}$$
(8)

with the corresponding tangent mappings

$$\begin{aligned} {{\mathbf {\mathsf{{K}}}}} = {\nabla _{\Theta }}\varvec{\kappa } = \varvec{G}_i \otimes \varvec{Z}^i \text { and } {{\mathbf {\mathsf{{M}}}}} = {\nabla _{\Theta }}\varvec{\mu } = \varvec{g}_i \otimes \varvec{Z}^i. \end{aligned}$$
(9)

Thus, the deformation mapping \(\mathbf {\varphi }\) and also the deformation gradient \({{\mathbf {\mathsf{{F}}}}} = {\nabla _{{X}}}\varvec{\varphi }\) can be decomposed and written as

$$\begin{aligned} \varvec{\varphi } = \varvec{\mu } \circ \varvec{\kappa }^\textsf {-1}\text { and } {{\mathbf {\mathsf{{F}}}}} = {{\mathbf {\mathsf{{M}}}}}\,{{\mathbf {\mathsf{{K}}}}}^\textsf {-1}= {{\mathbf {\mathsf{{g}}}}}_i \otimes {{\mathbf {\mathsf{{G}}}}}^i, \end{aligned}$$
(10)

cf. Fig. 2. Consequently, for volume mappings we obtain

$$\begin{aligned} \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_{{\Theta }}= \int \limits _{{\mathcal{B}}}J_M\,\mathrm{d}V_{{\Theta }}, \end{aligned}$$
(11)

with the determinants \(J = \det {{{\mathbf {\mathsf{{F}}}}}}\), \(J_K = \det {{\mathbf {\mathsf{{K}}}}}\) and \(J_M = \det {{\mathbf {\mathsf{{M}}}}}\). This viewpoint of kinematics has the main advantage that implicit dependencies do not arise until the definition of global equilibrium. Here, the ‘time-like’ variable s is used as a scalar design variable that parameterizes the material body in the reference configuration \({\mathcal{K}} = {\mathcal{K}}(s)\), as well as the referential points \(\varvec{X} = \varvec{X}(s)\).

Fig. 2
figure 2

Enhanced kinematics

Full size image

2.3 Variations and derivatives

The total variation of any quantity \((\bullet )(\varvec{u};\varvec{s};\varvec{h}_n)\) is defined as the sum of the partial variations

$$\begin{aligned} \delta (\bullet )(\varvec{u};\varvec{s};\varvec{h}_n)= & {} \delta _u (\bullet )(\varvec{u};\hat{\varvec{s}};\hat{\varvec{h}}_n)\nonumber \nonumber \\&+\, \delta _s (\bullet )(\hat{\varvec{u}};\varvec{s};\hat{\varvec{h}}_n)\nonumber \\&+\, \delta _{h_n} (\bullet )(\hat{\varvec{u}};\hat{\varvec{s}};\varvec{h}_n). \end{aligned}$$
(12)

A partial variation, e.g. w.r.t. \(\varvec{u}\) at fixed \(\hat{\varvec{s}}\) and \(\hat{\varvec{h}}_n\) is defined as the Gâteaux derivative

$$\begin{aligned} \delta _u(\bullet ) = \dfrac{\mathrm{d}}{\mathrm{d}\epsilon }(\bullet )(\varvec{u}+\epsilon \,\delta \varvec{u};\hat{\varvec{s}};\hat{\varvec{h}}_n)\Big |_{\epsilon =0}. \end{aligned}$$
(13)

This notation is used throughout the whole paper.

3 Constitutive equations and material

The elastoplastic material model can be found in the relevant literature, see e.g. [17, 40,41,42] and is briefly described here for better understanding of the sensitivity relations in Sect. 4.2. It is based on a multiplicative split of the deformation gradient \({{\mathbf {\mathsf{{F}}}}} = {{\mathbf {\mathsf{{F}}}}}_\textsf {e}\,{{\mathbf {\mathsf{{F}}}}}_\textsf {p}\) into an elastic, \({{\mathbf {\mathsf{{F}}}}}_\textsf {e}\), and plastic part, \({{\mathbf {\mathsf{{F}}}}}_\textsf {p}\), and captures finite deformations and strains.

The Kirchhoff stress tensor can be determined by differentiation of the strain energy function

$$\begin{aligned} \varvec{\tau } = \frac{\partial {W}}{\partial {{{\mathbf {\mathsf{{F}}}}}}}\,{{\mathbf {\mathsf{{F}}}}}^\textsf {T}. \end{aligned}$$
(14)

The strain energy function is given by

$$\begin{aligned} W= & {} U(J) + \widetilde{W}(\bar{{{\mathbf {\mathsf{{C}}}}}},\bar{{{\mathbf {\mathsf{{C}}}}}}_{\textsf {p}}^{\textsf {-1}})\nonumber \\= & {} \frac{K}{2}\,\left[ \frac{1}{2}\,(J^2-1)-\ln J\right] + \frac{G}{2}\,\left[ \bar{{{\mathbf {\mathsf{{C}}}}}} : \bar{{{\mathbf {\mathsf{{C}}}}}}_{\textsf {p}}^{\textsf {-1}}-3\right] , \end{aligned}$$
(15)

where \(\bar{{{\mathbf {\mathsf{{C}}}}}}_{\textsf {p}}^{\textsf {-1}}\) denotes the inverse of the isochoric, plastic contribution of the right Cauchy–Green tensor. K and G denote the compression and shear moduli, respectively. The von Mises yield condition

$$\begin{aligned} f^{\textsf {p}} = ||\,\varvec{\tau }_{\mathrm{dev}}\,|| - \sqrt{\frac{2}{3}}\,k(\alpha ) \end{aligned}$$
(16)

characterizes the plastic behavior, where \(\varvec{\tau }_{\mathrm{dev}}\) is the deviatoric part of the Kirchhoff stress tensor (14) and \(\alpha \) denotes the scalar isochoric hardening variable. The nonlinear function

$$\begin{aligned} k(\alpha ) = \sigma _0 + \sigma _\infty \,\left( 1-e^{-d\alpha }\right) + H\,\alpha \end{aligned}$$
(17)

describes hardening during plastic flow. The material parameters \(\sigma _0\) and \(\sigma _\infty \) represent the initial and saturation yield stress, respectively, H is the linear hardening slope and d is a dimensionless parameter. The deformation process is assumed to be quasi-static. However, a pseudo-time is introduced so as to solve the evolution of internal variables, namely \(\varvec{h} = \{\bar{{{\mathbf {\mathsf{{C}}}}}}_{\textsf {p}}^{\textsf {-1}}, \alpha \}\). The flow rules

$$\begin{aligned} \dot{\bar{{{\mathbf {\mathsf{{C}}}}}}}_{\textsf {p}}^{\textsf {-1}}= & {} -\frac{2}{3}\,\dot{\gamma }\,\left[ \bar{{{\mathbf {\mathsf{{C}}}}}} : \bar{{{\mathbf {\mathsf{{C}}}}}}_{\textsf {p}}^{\textsf {-1}}\right] \,{{\mathbf {\mathsf{{F}}}}}^\textsf {-1}\,{{\mathbf {\mathsf{{n}}}}}\,{{\mathbf {\mathsf{{F}}}}}^\textsf {-T}, \end{aligned}$$
(18)
$$\begin{aligned} \dot{\alpha }= & {} \sqrt{\frac{2}{3}}\,\dot{\gamma }, \end{aligned}$$
(19)

with the plastic multiplier \(\gamma \) and the flow direction

$$\begin{aligned} {{\mathbf {\mathsf{{n}}}}} = \nicefrac {\varvec{\tau }_{\mathrm{dev}}}{|| \varvec{\tau }_{\mathrm{dev}}||}, \end{aligned}$$
(20)

describe the evolution of the internal variables. A local Newton–Raphson scheme is used to solve the increment of the plastic multiplier within an implicit backward Euler time integration due to the nonlinear hardening function, cf. Eq. (17). The Kuhn–Tucker conditions

$$\begin{aligned} \Delta \gamma \ge 0, \qquad f^\textsf {p}(\varvec{\tau },\alpha ) \le 0,\qquad \Delta \gamma \,f^\textsf {p} = 0 \end{aligned}$$
(21)

define the solution point.

Onset and evolution of damage is characterized by the damage criterion

$$\begin{aligned} f^{\textsf {d}} = \alpha _{\textsf {d}}\,I_1 + \beta _{\textsf {d}}\,\sqrt{J_2} - \sigma _{\textsf {d}}=0 \end{aligned}$$
(22)

expressed in terms of the stress invariants \(I_1 = \mathrm{tr}\varvec{\tau }\) and \(J_2 = \frac{1}{2}\,\varvec{\tau }_{\mathrm{dev}}: \varvec{\tau }_{\mathrm{dev}}\) of the Kirchhoff stress tensor and the damage threshold \(\sigma _{\textsf {d}}\) , see [21,22,23] for further details. In Eq. (22) \(\alpha _{\textsf {d}}\) and \(\beta _{\textsf {d}}\) are the damage mode parameters corresponding to different stress state dependent damage processes acting on the micro-level. These parameters have been determined by numerical analysis of the deformation behavior of three-dimensionally loaded micro-defect-containing unit cells [11, 14]. The dependence of the functions on the stress state is expressed in terms of the stress triaxiality

$$\begin{aligned} \eta = \frac{\sigma _\textsf {m}}{\sigma _{\textsf {eq}}} = \frac{I_1}{3\,\sqrt{3\,J_2}}, \end{aligned}$$
(23)

with the mean stress \(\sigma _\textsf {m} = \frac{1}{3}\,I_1\), the von Mises equivalent stress \(\sigma _{\textsf {eq}} = \sqrt{3\,J_2}\), and the Lode parameter

$$\begin{aligned} \omega = \frac{2\,\tau _2 - \tau _1 - \tau _3}{\tau _1 - \tau _3}, \qquad \text {with }\qquad \tau _1 \ge \tau _2 \ge \tau _3, \end{aligned}$$
(24)

where \(\tau _1\), \(\tau _2\) and \(\tau _3\) denote the principal Kirchhoff stress components. For the experimental validation of the stress state dependent parameters \(\alpha _{\textsf {d}}\) and \(\beta _{\textsf {d}}\), tests with biaxially loaded specimens undergoing a wide range of stress states have been performed [9, 12, 13]. However, in some cases non-homogeneous stress distributions appear in critical parts where damage and final fracture take place leading to difficulties in analysis of failure modes. Therefore, besides the need for distinct stress states, a nearly homogeneous distribution of stress triaxiality (23) and the Lode parameter (24) are required in notched parts, which will be realized by optimization of the geometry of the biaxially loaded specimens.

Experiments and the corresponding optimization have been performed with specimens cutted from aluminium alloy (AlSi1MgMn, EN AW 6082-T6) sheets with a thickness of 4.0 mm and Fig. 3 indicates a corresponding load–displacement curve. The experimental data is taken from [24] and the material parameters have been determined by a curve fitting procedure based on gradient information obtained following the same variational principles as described in this paper for shape sensitivities. The procedure is described in [31] in detail. The parameters found are indicated in Table 1.

Table 1 Material parameters for AlSi1MgMn
Full size table
Fig. 3
figure 3

Load–displacement curves

Full size image

4 Structural analysis and sensitivity relations

4.1 Elastoplastic stress response

In view of the sensitivity analysis presented in Sect. 4.2, here, the most important equations for the computation of the structural response based on an implicit Return-Mapping scheme, see e.g. [15, 17, 41, 42], are briefly summarized. Note that for the computation of the evolution of the internal variables, cf. Eqs. (18) and (19), a pseudo-time is introduced and discretized using an implicit backward Euler method. In the following, all quantities are evaluated at the current pseudo-time step \(t_{n+1}\), except quantities with subscript n, which indicates that these quantities are saved from the prior pseudo-time step \(t_n\).

The nonlinear weak equilibrium condition in the reference configuration reads

$$\begin{aligned} R(\varvec{u},\varvec{v})= & {} \int \limits _{{\mathcal{K}}}{{\mathbf {\mathsf{{P}}}}}^{\textsf {K}}(\varvec{u}) : {\nabla _{{X}}}\varvec{v}\,\mathrm{d}V\nonumber \\&- \int \limits _{{\mathcal{K}}} \varvec{b} \cdot \varvec{v}\,\mathrm{d}V- \int \limits _{\partial {\mathcal{K}}} \varvec{t} \cdot \varvec{v}\,\mathrm{d}A\nonumber \\= & {} R^{\textsf {int}} + R^{\textsf {ext}} = 0, \end{aligned}$$
(25)

with the internal part \(R^{\textsf {int}}\) and the external part \(R^{\textsf {ext}}\). Here, \({{\mathbf {\mathsf{{P}}}}}^{\textsf {K}} = \varvec{\tau }\,{{\mathbf {\mathsf{{F}}}}}^\textsf {-T}\) denotes the first Piola–Kirchhoff stress tensor, and \(\varvec{u},\varvec{v} \in {\mathcal{V}}\) are the state and test function, respectively, where \({\mathcal{V}}\) is the space of admissible states. The Kirchhoff stress tensor can be calculated using a Return-Mapping-Algorithm, where first an elastic trial state is assumed, indicated in the following with the superscript \(\textsf {tr}\). If the trial state violates the yield condition, cf. Eq. (16), the plastic corrector step is performed, in which the plastic multiplier is calculated within a local Newton–Raphson procedure to update the internal variables and the stress tensor. Once the increment of the plastic multiplier is determined, the stress tensor can be adapted

$$\begin{aligned} \varvec{\tau }_{\mathrm{dev}}= \varvec{\tau }_{\mathrm{dev}}^{\mathrm{tr}} - 2\,\bar{\mu }\,\Delta \gamma \,{{\mathbf {\mathsf{{n}}}}}, \end{aligned}$$
(26)

with

$$\begin{aligned} \bar{\mu } = \frac{G}{3}\,\bar{{{\mathbf {\mathsf{{C}}}}}} : \bar{{{\mathbf {\mathsf{{C}}}}}}_{\textsf {p},n}^{\textsf {-1}}= \frac{G}{3}\,\mathrm{tr}{{{\mathbf {\mathsf{{b}}}}}_\textsf {e}^\textsf {tr}}. \end{aligned}$$
(27)

As the plastic flow is assumed to be incompressible, the hydrostatic stress can be easily computed to

$$\begin{aligned} \tau _{\textsf {vol}} = J\,U^\prime (J) = \frac{K}{2}\,\left( J^2 - 1\right) . \end{aligned}$$
(28)

Finally, the first Piola–Kirchhoff stress tensor reads

$$\begin{aligned} {{\mathbf {\mathsf{{P}}}}}^\textsf {K} = \left( \tau _{\textsf {vol}}\,{{\mathbf {\mathsf{{I}}}}} + \varvec{\tau }_{\mathrm{dev}}\right) \,{{\mathbf {\mathsf{{F}}}}}^\textsf {-T}= \varvec{\tau }\,{{\mathbf {\mathsf{{F}}}}}^\textsf {-T}. \end{aligned}$$
(29)

Consistent linearization leads to the explicit form of the material tangent

$$\begin{aligned} \mathbb {A} = \frac{\partial {{\mathbf {\mathsf{{P}}}}}^{\textsf {K}}}{\partial {{\mathbf {\mathsf{{F}}}}}}= & {} K\,J^2\,{{\mathbf {\mathsf{{F}}}}}^\textsf {-T}\otimes {{\mathbf {\mathsf{{F}}}}}^\textsf {-T}\nonumber \\&- \varvec{\tau }\,\left( {{\mathbf {\mathsf{{F}}}}}^\textsf {-T}\otimes {{\mathbf {\mathsf{{F}}}}}^\textsf {-T}\right) ^{{\mathop {\textsf {T}}\limits ^{24}}} + \mathbb {S}_{\mathrm{dev}}{{\,\mathrm{{\mathop {*}\limits ^{21}}}\,}}{{\mathbf {\mathsf{{F}}}}}^\textsf {-T}, \end{aligned}$$
(30)

see [31] for details. Here, the partial derivative of the Kirchhoff stress deviator is given by

$$\begin{aligned} \mathbb {S}_{\mathrm{dev}} = \frac{\partial {\varvec{\tau }_{\mathrm{dev}}}}{\partial {{{\mathbf {\mathsf{{F}}}}}}} = \mathbb {H} : \mathbb {B}, \end{aligned}$$
(31)

where

$$\begin{aligned} \mathbb {H} = G\,\left[ \beta _0\,\mathbb {I}_{d} - \frac{2}{3}\,\Delta \gamma \,\beta _1\,{{\mathbf {\mathsf{{n}}}}} \otimes {{\mathbf {\mathsf{{I}}}}} + \beta _2\,{{\mathbf {\mathsf{{n}}}}} \otimes {{\mathbf {\mathsf{{n}}}}}\right] , \end{aligned}$$
(32)

with the factors

$$\begin{aligned} \beta _0 = 1-\frac{2\,\Delta \gamma \,\overline{\mu }}{||\,\varvec{\tau }_{\mathrm{dev}}^\textsf {tr}\,||},\qquad \beta _1 = 1+\frac{2\,\overline{\mu }}{f^\prime }, \qquad \beta _2 = \beta _1-\beta _0 \end{aligned}$$
(33)

and the identities

$$\begin{aligned} \mathbb {B}:= & {} \frac{\partial {\bar{{{\mathbf {\mathsf{{b}}}}}}_{\textsf {e}}^{\textsf {tr}}}}{\partial {{{\mathbf {\mathsf{{F}}}}}}} = \left[ \mathbb {F}\,{{\,\mathrm{{\mathop {*}\limits ^{21}}}\,}}\,\bar{{{\mathbf {\mathsf{{C}}}}}}_{\textsf {p},n}^{\textsf {-1}}\,\bar{{{\mathbf {\mathsf{{F}}}}}}\quad {{{\mathbf {\mathsf{{F}}}}}}^\textsf {T}+ \bar{{{\mathbf {\mathsf{{F}}}}}}\,\bar{{{\mathbf {\mathsf{{C}}}}}}_{\textsf {p},n}^{\textsf {-1}}\,\mathbb {F}^{{\mathop {\textsf {T}}\limits ^{12}}}\right] , \end{aligned}$$
(34)
$$\begin{aligned} \mathbb {F}:= & {} \frac{\partial {\bar{{{\mathbf {\mathsf{{F}}}}}}}}{\partial {{{\mathbf {\mathsf{{F}}}}}}} J^{-\frac{1}{3}}\,\left[ \left( {{\mathbf {\mathsf{{I}}}}} \otimes {{\mathbf {\mathsf{{I}}}}}\right) ^{{\mathop {\textsf {T}}\limits ^{23}}} - \frac{1}{3}\,{{\mathbf {\mathsf{{F}}}}}\otimes {{\mathbf {\mathsf{{F}}}}}^\textsf {-T}\right] . \end{aligned}$$
(35)

The tensor \(\bar{{{\mathbf {\mathsf{{F}}}}}} = J^{-\frac{1}{3}}\,{{\mathbf {\mathsf{{F}}}}}\) denotes the isochoric contribution of the deformation gradient.

4.2 Sensitivity analysis

Utilizing gradient based optimization algorithms, gradient information of objective and constraint functions have to be provided to the optimization solver. This section outlines the computation of sensitivities using a variational approach at continuous level, i.e. continuous in space.

4.2.1 Sensitivity of the stress state

As aforementioned, the quantity of interest is the stress triaxiality, cf. Eq. (23). Thus, the computation of the gradient of the stress triaxiality is essential. Its total variation reads

$$\begin{aligned} \delta \eta= & {} \frac{\partial {\eta }}{\partial {\sigma _{\textsf {m}}}}\,\delta \sigma _{\textsf {m}} + \frac{\partial {\eta }}{\partial {\sigma _{\textsf {eq}}}}\,\delta \sigma _{\textsf {eq}}\nonumber \\= & {} \frac{1}{\sigma _{\textsf {eq}}}\,\left( \delta \sigma _{\textsf {m}} - \eta \,\delta \sigma _{\textsf {eq}}\right) , \end{aligned}$$
(36)

where the total variations \(\delta \sigma _{\textsf {m}}\) and \(\delta \sigma _{\textsf {eq}}\) are given by

$$\begin{aligned} \delta \sigma _{\textsf {m}}= & {} \frac{\partial {\sigma _{\textsf {m}}}}{\partial {\varvec{\tau }}} : \delta \varvec{\tau } = \frac{1}{3}\,{{\mathbf {\mathsf{{I}}}}} : \delta \varvec{\tau }, \end{aligned}$$
(37)
$$\begin{aligned} \delta \sigma _{\textsf {eq}}= & {} \frac{\partial {\sigma _{\textsf {eq}}}}{\partial {\varvec{\tau }}} : \delta \varvec{\tau } = \frac{3}{2\,\sigma _{\textsf {eq}}}\,\varvec{\tau }_{\mathrm{dev}}: \delta \varvec{\tau }. \end{aligned}$$
(38)

By means of Eq. (31), the total variation of the Kirchhoff stress tensor can be revealed as

$$\begin{aligned} \delta \varvec{\tau } = \left[ K\,J^2\,{{\mathbf {\mathsf{{I}}}}} \otimes {{\mathbf {\mathsf{{F}}}}}^\textsf {-T}+ \mathbb {S}_{\mathrm{dev}}\right] : \delta {{\mathbf {\mathsf{{F}}}}} + \frac{\partial {\varvec{\tau }}}{\partial {\varvec{h}_n}}\,\delta \varvec{h}_n. \end{aligned}$$
(39)

As geometric design changes are considered, the design vector is chosen as the referential coordinates \(\varvec{X}\). From literature, e.g. [2, 33], we find

$$\begin{aligned} \delta {{\mathbf {\mathsf{{F}}}}} = \delta _u{{\mathbf {\mathsf{{F}}}}} + \delta _X{{\mathbf {\mathsf{{F}}}}} = {\nabla _{{X}}}\delta \varvec{u} - {\nabla _{{X}}}\varvec{u}\,{\nabla _{{X}}}\delta \varvec{X}. \end{aligned}$$
(40)

The second addend in Eq. (39) corresponds to the deformation history that is captured by the internal variables of the prior load step \(\varvec{h}_n = \{\bar{{{\mathbf {\mathsf{{C}}}}}}_{\textsf {p},n}^{\textsf {-1}}, \alpha _n\}\). Note that this notation is an abbreviation for

$$\begin{aligned} \frac{\partial {\varvec{\tau }}}{\partial {\varvec{h}_n}}\,\delta \varvec{h}_n := \frac{\partial {\varvec{\tau }}}{\partial {\bar{{{\mathbf {\mathsf{{C}}}}}}_{\textsf {p},n}^{\textsf {-1}}}} : \delta \bar{{{\mathbf {\mathsf{{C}}}}}}_{\textsf {p},n}^{\textsf {-1}}+ \frac{\partial {\varvec{\tau }}}{\partial {\alpha _n}}\,\delta \alpha _n. \end{aligned}$$
(41)

The partial derivatives of the Kirchhoff stress tensor w.r.t. the internal variables read

$$\begin{aligned} \frac{\partial {\varvec{\tau }}}{\partial {\bar{{{\mathbf {\mathsf{{C}}}}}}_{\textsf {p},n}^{\textsf {-1}}}}= \mathbb {H} : \mathbb {B}_C \qquad \text {and }\qquad \frac{\partial {\varvec{\tau }}}{\partial {\alpha _n}} = -2\,\sqrt{\frac{2}{3}}\,\bar{\mu }\,\frac{k^\prime }{f^\prime }\,{{\mathbf {\mathsf{{n}}}}}, \end{aligned}$$
(42)

with the tensor

$$\begin{aligned} \mathbb {B}_C = \bar{{{\mathbf {\mathsf{{F}}}}}}\,\mathbb {I}_\textsf {s} {{\,\mathrm{{\mathop {*}\limits ^{21}}}\,}}\bar{{{\mathbf {\mathsf{{F}}}}}}^\textsf {T}\end{aligned}$$
(43)

and the first derivatives of Eqs. (16) and (17)

$$\begin{aligned} k^\prime= & {} \sigma _\infty \,d\,e^{(-d\,\alpha )}+H, \end{aligned}$$
(44)
$$\begin{aligned} f^\prime= & {} -2\,\bar{\mu }\,\left[ 1+\frac{k^\prime }{3\,\bar{\mu }}\right] . \end{aligned}$$
(45)

In Eq. (40), the total variation of the structural response \(\delta \varvec{u}\) occurs. This implies that the total variation of the structural response has to be computed before the gradient information of the stress state can be obtained.

4.2.2 Structural response sensitivity

The weak equilibrium condition, Eq. (25), has to hold for any change in design. Thus, its total variation has to vanish, cf. [2, 30, 31, 44]. We assume that external forces are design independent (\(\delta R^{\textsf {ext}} = 0 \rightarrow \delta R = \delta R^{\textsf {int}}\)). Consequently, the total variation of the weak equilibrium reads

$$\begin{aligned}&\delta R(\varvec{v},\delta \varvec{u},\delta \varvec{X},\delta \varvec{h}_n) = \delta _u R^\textsf {int} + \delta _X R^\textsf {int} + \delta _{h_n} R^\textsf {int}\nonumber \\&\quad = k(\varvec{v},\delta \varvec{u}) + p(\varvec{v},\delta \varvec{X}) + h(\varvec{v},\delta \varvec{h}_n) = 0. \end{aligned}$$
(46)

The bilinear forms \(k : {\mathcal{V}} \times {\mathcal{V}} \rightarrow \mathbb {R}\), \(p : {\mathcal{V}} \times {\mathcal{X}} \rightarrow \mathbb {R}\) and \(h : {\mathcal{V}} \times {\mathcal{H}} \rightarrow \mathbb {R}\) denote the tangent stiffness, pseudo load and history sensitivity, respectively.

Tangent stiffness. The tangential stiffness operator represents the partial variation of the internal residual \(R^{\textsf {int}}\) w.r.t. the displacements \(\varvec{u}\)

$$\begin{aligned} k(\varvec{v},\delta \varvec{u}) = \int \limits _{{\mathcal{K}}} {\nabla _{{X}}}\delta \varvec{u} : \mathbb {A} : {\nabla _{{X}}}\varvec{v}\,\mathrm{d}V. \end{aligned}$$
(47)

Pseudo load. The pseudoload operator is the partial variation of the internal residual \(R^{\textsf {int}}\) w.r.t. design \(\varvec{X}\). Based on the enhanced viewpoint of kinematics described in Sect. 2, a pull-back operation to parameter space \({\mathcal{B}}\) can be performed and all required variations can easily be obtained, see e.g. [2, 26, 31, 33] for details. Finally, a push-forward to the reference configuration leads to

$$\begin{aligned} p(\varvec{v},\delta \varvec{X})= & {} - \int \limits _{\partial {\mathcal{K}}} \left[ \left( {\nabla _{{X}}}\varvec{u}\,{\nabla _{{X}}}\delta \varvec{X}\right) : \mathbb {A}\right] : {\nabla _{{X}}}\varvec{v}\,\mathrm{d}V\nonumber \\&- \int \limits _{\partial {\mathcal{K}}} {{\mathbf {\mathsf{{P}}}}}^{\textsf {K}} : ({\nabla _{{X}}}\varvec{v}\,{\nabla _{{X}}}\delta \varvec{X})\,\mathrm{d}V\nonumber \\&+ \int \limits _{\partial {\mathcal{K}}} {{\mathbf {\mathsf{{P}}}}}^{\textsf {K}} : {\nabla _{{X}}}\varvec{v}\,{\nabla _{{X}}}\cdot \,\delta \varvec{X}\,\mathrm{d}V. \end{aligned}$$
(48)

History sensitivity. The partial variation of the internal residual \(R^{\textsf {int}}\) w.r.t. the history variables of the prior pseudo-time step \(\varvec{h}_n\) is denoted as the history sensitivity operator. Here, only the first Piola–Kirchhoff stress tensor depends on the history variables. Thus,

$$\begin{aligned} h(\varvec{v},\delta \varvec{h}_n) =\int \limits _{\partial {\mathcal{K}}} \delta _{h_n}{{\mathbf {\mathsf{{P}}}}}^{\textsf {K}} : {\nabla _{{X}}}\varvec{\eta }\,\mathrm{d}V. \end{aligned}$$
(49)

We can identify the partial variation to

$$\begin{aligned} \delta _{h_n}{{\mathbf {\mathsf{{P}}}}}^{\textsf {K}} = \frac{\partial {{{\mathbf {\mathsf{{P}}}}}^\textsf {K}}}{\partial {\varvec{h}_n}} \,\delta \varvec{h}_n \end{aligned}$$
(50)

and determine the partial derivatives of the stress tensor w.r.t. the history variables \(\varvec{h}_n = \{\bar{{{\mathbf {\mathsf{{C}}}}}}_{\textsf {p},n}^{\textsf {-1}}, \alpha _n\}\) to

$$\begin{aligned} \frac{\partial {{{\mathbf {\mathsf{{P}}}}}^{\textsf {K}}}}{\partial {\bar{{{\mathbf {\mathsf{{C}}}}}}_{\textsf {p},n}^{\textsf {-1}}}}= & {} \left( \mathbb {H} : \mathbb {B}_C\right) {{\,\mathrm{{\mathop {*}\limits ^{21}}}\,}}{{\mathbf {\mathsf{{F}}}}}^\textsf {-T}, \end{aligned}$$
(51)
$$\begin{aligned} \frac{\partial {{{\mathbf {\mathsf{{P}}}}}^\textsf {K}}}{\partial {\alpha _n}}= & {} -2\,\sqrt{\frac{2}{3}}\,\bar{\mu }\,\frac{k^\prime }{f^\prime }\,{{\mathbf {\mathsf{{n}}}}}\,{{\mathbf {\mathsf{{F}}}}}^\textsf {-T}, \end{aligned}$$
(52)

respectively.

Update of history sensitivities. At the end of each plastic step the internal history variables evolve. Thus, its variations w.r.t. design changes have to be computed and saved for the subsequent step, cf. Eq. (50). Recalling Eq. (40), the algorithmic update formula reads

$$\begin{aligned} \delta \varvec{h} = \frac{\partial {\varvec{h}}}{\partial {{{\mathbf {\mathsf{{F}}}}}}} : \delta {{\mathbf {\mathsf{{F}}}}} + \frac{\partial {\varvec{h}}}{\partial {\varvec{h}_n}}\,{{\mathbf {\mathsf{{Z}}}}}_n\,\delta \varvec{X} = {{\mathbf {\mathsf{{Z}}}}}\,\delta \varvec{X}, \end{aligned}$$
(53)

where also the relation \(\delta \varvec{u} = \mathbf{\mathrm{S} }\,\delta \varvec{X}\) has been used, which is explained later. The only unknown quantities are the partial derivatives \(\displaystyle \frac{\partial {\varvec{h}}}{\partial {{{\mathbf {\mathsf{{F}}}}}}}\) and \(\displaystyle \frac{\partial {\varvec{h}}}{\partial {\varvec{h}_n}}\). For the former we find

$$\begin{aligned} \frac{\partial {\bar{{{\mathbf {\mathsf{{C}}}}}}_{\textsf {p}}^{\textsf {-1}}}}{\partial {{{\mathbf {\mathsf{{F}}}}}}}= & {} -\left[ \left( \bar{{{\mathbf {\mathsf{{F}}}}}}^\textsf {-1}\otimes \bar{{{\mathbf {\mathsf{{F}}}}}}^\textsf {-T}\right) ^{{\mathop {\textsf {T}}\limits ^{23}}} : \mathbb {F}\right] {{\,\mathrm{{\mathop {*}\limits ^{21}}}\,}}\left( \bar{{{\mathbf {\mathsf{{b}}}}}}_{\textsf {e}}\,\bar{{{\mathbf {\mathsf{{F}}}}}}^\textsf {-T}\right) \nonumber \\&+ \bar{{{\mathbf {\mathsf{{F}}}}}}^\textsf {-1}\left[ \left( \frac{1}{G}\,\mathbb {H} + \frac{1}{3}\,{{\mathbf {\mathsf{{I}}}}} \otimes {{\mathbf {\mathsf{{I}}}}}\right) : \mathbb {B}\right] {{\,\mathrm{{\mathop {*}\limits ^{21}}}\,}}\bar{{{\mathbf {\mathsf{{F}}}}}}^\textsf {-T}\nonumber \\&- \bar{{{\mathbf {\mathsf{{F}}}}}}^\textsf {-1}\,\bar{{{\mathbf {\mathsf{{b}}}}}}_{\textsf {e}}\,\left[ \left( \bar{{{\mathbf {\mathsf{{F}}}}}}^\textsf {-T}\otimes \bar{{{\mathbf {\mathsf{{F}}}}}}^\textsf {-T}\right) ^{{\mathop {\textsf {T}}\limits ^{24}}} : \mathbb {F}\right] , \end{aligned}$$
(54)
$$\begin{aligned} \frac{\partial {\alpha }}{\partial {{{\mathbf {\mathsf{{F}}}}}}}= & {} \sqrt{\frac{2}{3}}\,\frac{G}{f^\prime }\,\left( \frac{2}{3}\,\Delta \gamma \,{{\mathbf {\mathsf{{I}}}}} - {{\mathbf {\mathsf{{n}}}}}\right) : \mathbb {B}, \end{aligned}$$
(55)

for the internal variables \(\bar{{{\mathbf {\mathsf{{C}}}}}}_{\textsf {p}}^{\textsf {-1}}\) and \(\alpha \), respectively, w.r.t. the deformation gradient. The latter can be identified to

$$\begin{aligned} \frac{\partial {\bar{{{\mathbf {\mathsf{{C}}}}}}_{\textsf {p}}^{\textsf {-1}}}}{\partial {\bar{{{\mathbf {\mathsf{{C}}}}}}_{\textsf {p},n}^{\textsf {-1}}}}= & {} \mathbb {B}_C^\textsf {-1}: \left[ \frac{1}{G}\,\mathbb {H} + \frac{1}{3}\,{{\mathbf {\mathsf{{I}}}}}\otimes {{\mathbf {\mathsf{{I}}}}} \right] : \mathbb {B}_C, \end{aligned}$$
(56)
$$\begin{aligned} \frac{\partial {\bar{{{\mathbf {\mathsf{{C}}}}}}_{\textsf {p}}^{\textsf {-1}}}}{\partial {\alpha _n}}= & {} -2\,\sqrt{\frac{2}{3}}\,\frac{k^\prime }{f^\prime }\,\mathbb {B}_C^\textsf {-1}: {{\mathbf {\mathsf{{n}}}}}, \end{aligned}$$
(57)

with the fourth order tensor

$$\begin{aligned} \mathbb {B}_C^\textsf {-1}= \bar{{{\mathbf {\mathsf{{F}}}}}}^\textsf {-1}\mathbb {I}_s {{\,\mathrm{{\mathop {*}\limits ^{21}}}\,}}\bar{{{\mathbf {\mathsf{{F}}}}}}^\textsf {-T}, \end{aligned}$$
(58)

and

$$\begin{aligned} \frac{\partial {\alpha }}{\partial {\bar{{{\mathbf {\mathsf{{C}}}}}}_{\textsf {p},n}^{\textsf {-1}}}}= & {} \sqrt{\frac{2}{3}}\,\frac{G}{f^\prime }\,\left[ \frac{2}{3}\,\Delta \gamma \,{{\mathbf {\mathsf{{I}}}}} - {{\mathbf {\mathsf{{n}}}}}\right] : \mathbb {B}_C \end{aligned}$$
(59)
$$\begin{aligned} \frac{\partial {\alpha }}{\partial {\alpha _n}}= & {} 1 + \frac{2}{3}\,\frac{k^\prime }{f^\prime }. \end{aligned}$$
(60)

4.3 Discrete matrix forms

With the approximate functions \(\varvec{u}_h, \varvec{v}_h \in {\mathcal{V}}_h \subset {\mathcal{V}}\), \(\varvec{X}_h \in {\mathcal{X}}_h \subset {\mathcal{X}}\) and \(\varvec{h}_{n,h} \in {\mathcal{H}}_h \subset {\mathcal{H}}\) based on a conforming Galerkin method, cf. e.g. [6, 45], and the discrete approximations, i.e. the nodal vector of test functions, displacements and displacement variation \(\mathbf{\mathrm{\mathbf{v}} },\mathbf{\mathrm{\mathbf{u}} }, \delta \mathbf{\mathrm{\mathbf{u}} } \in \mathbb {R}^{\textsf {ndof}}\), the nodal design vector and its variation \(\mathbf{\mathrm{\mathbf{X}} }, \delta \mathbf{\mathrm{\mathbf{X}} } \in \mathbb {R}^{\textsf {ndv}}\), as well as the vector of variations of global internal variables \(\delta \mathbf{\mathrm{\mathbf{h}} }_{n,h} \in \mathbb {R}^{\textsf {nhv}}\), we find the discrete versions of Eqs. (25), (47), (48) and (49), respectively

$$\begin{aligned}&R(\varvec{u}_h,\varvec{v_h}) = \mathbf{\mathrm{\mathbf{v}} }^\textsf {T}\,\mathbf{\mathrm{\mathbf{R}} }, \end{aligned}$$
(61)
$$\begin{aligned}&k(\varvec{v_h},\delta \varvec{u}_h) = \mathbf{\mathrm{\mathbf{v}} }^\textsf {T}\,\mathbf{\mathrm{\mathbf{K}} }\,\delta \mathbf{\mathrm{\mathbf{u}} }, \end{aligned}$$
(62)
$$\begin{aligned}&p(\varvec{v_h},\delta \varvec{X}_h) = \mathbf{\mathrm{\mathbf{v}} }^\textsf {T}\,\mathbf{\mathrm{\mathbf{P}} }\,\delta \mathbf{\mathrm{\mathbf{X}} }, \end{aligned}$$
(63)
$$\begin{aligned}&h(\varvec{v_h},\delta \varvec{h}_{n,h}) = \mathbf{\mathrm{\mathbf{v}} }^\textsf {T}\,\mathbf{\mathrm{\mathbf{H}} }\,\delta \mathbf{\mathrm{\mathbf{h}} }_n. \end{aligned}$$
(64)

The global quantities \(\textsf {ndof}\), \(\textsf {nhv}\) and \(\textsf {ndv}\) correspond to the global number of degrees of freedom, history variables and design variables, respectively. With this, the discrete matrix form of Eq. (46) reads

$$\begin{aligned} \delta \mathbf{\mathrm{\mathbf{R}} }= \mathbf{\mathrm{\mathbf{K}} }\,\delta \mathbf{\mathrm{\mathbf{u}} } + \mathbf{\mathrm{\mathbf{P}} }\,\delta \mathbf{\mathrm{\mathbf{X}} } + \mathbf{\mathrm{\mathbf{H}} }\,\delta \mathbf{\mathrm{\mathbf{h}} }_n = \mathbf{\mathrm{0} }, \end{aligned}$$
(65)

with the tangent stiffness matrix \(\mathbf{\mathrm{\mathbf{K}} } \in \mathbb {R}^{\textsf {ndof} \times \textsf {ndof}}\), the pseudoload matrix \(\mathbf{\mathrm{\mathbf{P}} } \in \mathbb {R}^{\textsf {ndof} \times \textsf {ndv}}\) and the history sensitivity matrix \(\mathbf{\mathrm{\mathbf{H}} } \in \mathbb {R}^{\textsf {ndof} \times \textsf {nhv}}\). Rearranging Eq. (65), the total sensitivity matrix \(\mathbf{\mathrm{\mathbf{S}} } \in \mathbb {R}^{\textsf {ndof} \times \textsf {ndv}}\) can be derived

$$\begin{aligned} \delta \mathbf{\mathrm{\mathbf{u}} }= & {} -\mathbf{\mathrm{\mathbf{K}} }^\textsf {-1}\,\left[ \mathbf{\mathrm{\mathbf{P}} }\,\delta \mathbf{\mathrm{s} } + \mathbf{\mathrm{\mathbf{H}} }\,\delta \mathbf{\mathrm{\mathbf{h}} }_n\right] \nonumber \\= & {} -\mathbf{\mathrm{\mathbf{K}} }^\textsf {-1}\,\left[ \mathbf{\mathrm{\mathbf{P}} }+\mathbf{\mathrm{\mathbf{H}} }\,\mathbf{\mathrm{\mathbf{Z}} }_n\right] \,\delta \mathbf{\mathrm{\mathbf{X}} } = \mathbf{\mathrm{\mathbf{S}} }\,\delta \mathbf{\mathrm{\mathbf{X}} }. \end{aligned}$$
(66)

It connects design changes with changes in the structural response. Note that \(\mathbf{\mathrm{\mathbf{Z}} }_n\) is the matrix form of the quantity \({{\mathbf {\mathsf{{Z}}}}}_n\) in Eq. (53) and is updated algorithmically. Note that within a finite element procedure, the product \(\mathbf{\mathrm{\mathbf{Q}} } = \mathbf{\mathrm{\mathbf{H}} }\,\mathbf{\mathrm{\mathbf{Z}} }_n = \bigcup _e \mathbf{\mathrm{\mathbf{H}} }_e\,\mathbf{\mathrm{\mathbf{Z}} }_\mathbf{n,e} \) can be computed on element level, which makes the assembly of \(\mathbf{\mathrm{\mathbf{H}} }\) unnecessary and therefore saves memory, as for each integration point of each finite element, the history variables have to be stored and the total number of history variables \(\textsf {nhv}\) can be huge, depending on the finite element discretization. With the assumption \(\mathbf{\mathrm{\mathbf{h}} }_0 = \delta \mathbf{\mathrm{\mathbf{h}} }_0 = \mathbf{\mathrm{\mathbf{0}} }\) at time \(t = t_0\), from Eq. (65), we obtain

$$\begin{aligned} \mathbf{\mathrm{\mathbf{S}} }_0 = -\,\mathbf{\mathrm{\mathbf{K}} }^\textsf {-1}\,\mathbf{\mathrm{\mathbf{P}} }\qquad \text {and }\qquad \mathbf{\mathrm{\mathbf{Z}} }_0 = \mathbf{\mathrm{\mathbf{0}} } \end{aligned}$$
(67)

for the first pseudo-time increment. Hence, the sensitivity part corresponding to the deformation history is considered in each subsequent load step that has caused plastic yielding.

The overall computational procedure for a pseudo-time interval \(\left[ 0, N\right] \) is illustrated in the flow chart in Fig. 4. Note that the mathematical optimizer in step 4 can be chosen differently depending on the problem size and the form of the objective and constraint functions. See e.g. [31] for details on the implementation utilizing the Matlab function fmincon available in the Matlab Optimization Toolbox. Further details on the approach, discretization using \(\bar{F}\)-finite elements, cf. [15], and important implementation aspects can be found in [31].

Fig. 4
figure 4

Overall computational procedure

Full size image

5 Optimization problems

5.1 Initial geometry

The initial geometry of the X0-specimen with notches inclined by 45\(^{\circ }\) to the loading axis has outer dimensions of 240 mm by 240 mm, see Fig. 1. The depth of the notches is 1 mm, reducing the thickness here from 4 to 2 mm at its thinnest point, Fig. 5. First polished micrograph images indicate a homogeneous grain size distribution and consequently the remaining material can be seen as representative. The connectors between the specimen legs have a width of 6 mm and are centrally arranged, quite robust and thus, reduce the liability to fabrication inaccuracies .

Fig. 5
figure 5

a Initial geometry, b Optimized geometry V1 (LC 1/1), c Optimized geometry V2 (LC 1/-1)

Full size image

5.2 Forward problem

An FE model of the initial specimen is given in Fig. 6. Here, the symmetries in longitudinal, lateral and thickness direction are exploited. Consequently, only one eighth of the whole geometry is modeled. Symmetric boundary conditions are adopted at the symmetry planes and predefined displacements are applied at the top and the right edges of the mesh. In total, 5376 \(\bar{F}\)-finite elements are used and the mesh counts 21,375 degrees of freedom.

Two different load cases are considered, namely (LC 1/1) and (LC 1/-1). In the first case, both axes are equivalently loaded in positive direction (\(\bar{u}_2=\bar{u}_1\)), which produces a high tensile stress state in the notched specimen area. In the latter case, the axis 2 is loaded in negative direction (\(\bar{u}_2=-\bar{u}_1\)), which leads to a shear stress state in the notched area. The loads are applied in terms of prescribed displacements. The maximal values are \(\bar{u}^{\textsf {max}}_{\text {(1/1)}} = {0.15}\) mm and \(\bar{u}^{\textsf {max}}_{\text {(1/-1)}} = \pm \,{0.625}\) mm.

The forward boundary value problem is defined by Eq. (25) together with the boundary conditions indicated in Fig. 6. With the solution vector \(\mathbf{\mathrm{\mathbf{u}} }(\mathbf{\mathrm{\mathbf{p}} })\), the vector of stress triaxiality values at the cross section of the specimen \(\varvec{\eta }^{\textsf {cs}}(\varvec{\mathbf{\mathrm{\mathbf{u}} }}(\varvec{\mathbf{\mathrm{\mathbf{p}} }}))\) can be computed in a postprocessing step.

Note that the vector \(\mathbf{\mathrm{\mathbf{p}} }\) stores the geometric design variables described in the next section. As the sensitivity relations have been derived w.r.t. the referential coordinates, a connection between these coordinates and the design variables has to be established, which leads to the so-called design velocity matrix \(\mathbf{\mathrm{\mathbf{Q}} }\)

$$\begin{aligned} \delta \mathbf{\mathrm{\mathbf{X}} } = \frac{\partial {\mathbf{\mathrm{\mathbf{X}} }}}{\partial {\mathbf{\mathrm{\mathbf{p}} }}}\,\delta \mathbf{\mathrm{\mathbf{p}} } = \mathbf{\mathrm{\mathbf{Q}} }\,\delta \mathbf{\mathrm{\mathbf{p}} }, \end{aligned}$$
(68)

described in detail in [31]. Consequently, in the following the number of design variables \(\textsf {ndv}\) is the length of the design vector \(\mathbf{\mathrm{\mathbf{p}} }\).

Fig. 6
figure 6

FE mesh and boundary conditions

Full size image

5.3 Inverse problem

As aforementioned, the shape of the X0-specimen is to be changed so as to gain a distinct (high or low) preferably homogeneous stress triaxiality in the cross section of the notched area. These geometry changes have to be fulfilled under consideration of the producibility at reasonable costs; i.e. the fabrication process constrains the design variables in a certain way. Consequently, as design variables the inner and outer radii \(R_i\), \(R_o\) as well as the radius in thickness direction and the penetration depth \(R_t\) and D are chosen, see Fig. 7. Thus, the design vector reads \(\mathbf{\mathrm{\mathbf{p}} } := \left[ {\begin{array}{llll} R_i&R_o&R_t&D \end{array}} \right] .\)

Fig. 7
figure 7

Design variables in the notched area

Full size image

For each of the mentioned load cases, one individual optimum is desired. Thus, we state the two optimization problems for the respective load case.

Load case (1/1)

In this load case, the stress state is desired to be tensile dominated. Thus, the norm of the vector \(\varvec{\eta }^{\textsf {cs}}\) that collects the stress triaxiality values in this area is to be maximized.

Problem 1

$$\begin{aligned} \max _{X}&J\left( \mathbf{\mathrm{\mathbf{u}} }(\mathbf{\mathrm{\mathbf{p}} })\right) = ||\,\varvec{\eta }^{\textsf {cs}}\left( \mathbf{\mathrm{\mathbf{u}} }(\mathbf{\mathrm{\mathbf{p}} })\right) \,||_{_2}\nonumber \\ \text {s.t. }&c^\textsf {eq}(\mathbf{\mathrm{\mathbf{p}} }) = 0,\nonumber \\&\mathbf{c}^\textsf {in}(\mathbf{\mathrm{\mathbf{p}} }) \le 0, \nonumber \\&p^i_\textsf {l} \le p^i \le p^i_\textsf {u}, \quad i=1, \dots \textsf {ndv}, \end{aligned}$$
(69)

with the box constraints

$$\begin{aligned} \mathbf{\mathrm{\mathbf{p}} }_\textsf {l} = \left[ {\begin{array}{llll} 1.0&1.0&1.0&1.0 \end{array}} \right] , \mathbf{\mathrm{\mathbf{p}} }_\textsf {u} = \left[ {\begin{array}{llll} 4.0&4.0&3.0&1.5 \end{array}} \right] . \end{aligned}$$

Load case (–1/1)

In contrast to the former load case, by inverting the loading direction of one axis, the stress state in the notched specimen area is shear dominated. Hence, a minimization of the stress triaxiality in that area is desired.

Problem 2

$$\begin{aligned} \min _{X}&J(\mathbf{\mathrm{\mathbf{u}} }(\mathbf{\mathrm{\mathbf{p}} })) = ||\,\varvec{\eta }^{\textsf {cs}}(\mathbf{\mathrm{\mathbf{u}} }(\mathbf{\mathrm{\mathbf{p}} }))\,||_{_2}\nonumber \\ \text {s.t. }&c^\textsf {eq}(\mathbf{\mathrm{\mathbf{p}} }) = 0,\nonumber \\&\mathbf{c}^\textsf {in}(\mathbf{\mathrm{\mathbf{p}} }) \le 0, \nonumber \\&p^i_\textsf {l} \le p^i \le p^i_\textsf {u}, \quad i=1, \dots \textsf {ndv}, \end{aligned}$$
(70)

with the box constraints

$$\begin{aligned} \mathbf{\mathrm{\mathbf{p}} }_\textsf {l} = \left[ {\begin{array}{llll} 1.0&1.0&1.0&0.5 \end{array}} \right] , \mathbf{\mathrm{\mathbf{p}} }_\textsf {u} = \left[ {\begin{array}{llll} 5.5&5.5&4.0&1.5 \end{array}} \right] . \end{aligned}$$

In both cases, the equality and inequality constraints are identical. The equality constraint forces the cross section area in the notched region to stay constant at 12 mm\(^{2}\), i.e.

$$\begin{aligned} c^{\textsf {eq}} = l^{\textsf {cs}} \cdot t^{\textsf {cs}} - {12}\,{\hbox {mm}^{2}} = 0. \end{aligned}$$
(71)

Three inequality constraints are chosen to prevent destruction of the geometry and FE mesh due to design changes. Clearly, the inner and outer radii should always be greater than the radius in thickness direction and the penetration depth should always be smaller than the radius in thickness direction. This leads to

$$\begin{aligned} \mathbf{\mathrm{\mathbf{c}} }^{\textsf {in}} = \left[ {\begin{array}{l} R_t - R_i \\ R_t - R_o \\ E - R_t \end{array}} \right] \le \varvec{0}. \end{aligned}$$
(72)

The optimization problems are solved utilizing the interior-point solver of the internal Matlab function fmincon, available in the Optimization Toolbox. Analytical gradients are provided to fmincon that are computed utilizing the method described in Sect. 4.2.

5.4 Optimization results

The solutions of Problem 1 and Problem 2 are given in Table 2 and illustrated in Fig. 5b, c. The optimized geometry for the load case 1/1 is denoted with V1 and for the load case \(-1/1\) with V2. V1 (Fig. 5b) is characterized by a smaller notch radius in thickness direction \(R_t={1.0}\) mm, whereas the penetration depth remains \(D={1.0}\) mm. Furthermore, the outer radius has been reduced to \(R_o={1.0}\) mm and the inner radius remains \(R_i={3.0}\) mm. Due to the reduced notch radii, a more elevated stress triaxiality under 1/1 loading and a reduced zone with increased strains can be expected. In contrast, V2 is characterized by a reduced penetration depth of \(D={0.5}\) mm which leads to a smaller width of the notched area of 4.0 mm based on the constant cross section of \({12.0}\,{\hbox {mm}}^2\), see Eq. (71). In addition, the notch radii \(R_t={3.0}\) mm and \(R_o={4.0}\) mm have been enlarged whereas \(R_i={3.0}\) mm remains constant. Consequently, it can be expected that V2 conducts ’softer’, i.e. zones with elevated strain increase and under 1/1 loading a reduced stress triaxiality can be expected.

Table 2 Initial and optimal values of design variables for the two different load cases
Full size table

The relative change in the objective function is given as

$$\begin{aligned} c_{\textsf {rel}} = \frac{|\,J_0 - J^*\,|}{|\,J_0\,|} = \frac{\Delta J}{|\,J_0\,|}, \end{aligned}$$
(73)

were \(J^*\) denotes the value of the objective function in the solution point. In addition, Fig. 8 summarizes the value of the objective function during the optimization process. Note that for Problem 1 the maximization problem has been transformed into a minimization problem. Thus, the value of the objective function in Fig. 8a is negative. For both load cases the optimization procedure stops if either the step tolerance is less than 1e–6 or if the first-order optimality criterion is less than 1e–6, see e.g. [36] for details on the stopping criteria.

Fig. 8
figure 8

Objective history

Full size image

Figure 9 displays the stress triaxiality \(\eta \) (see Eq. 23) and Fig. 10 the Lode parameter (see Eq. 24) in the notched area of the X0-specimen for the initial geometry ((a) and (d)) as well as for the optimized geometries V1 ((b) and (e)) and V2 ((c) and (f)) under 1/1 loading ((a–c)) and \(-1/1\) loading ((d–f)) at the end of the simulations. Consequently, Fig. 9b displays the objective quantity \(\eta \) on the corresponding optimized geometry for 1/1 loading and Fig. 9f for \(-1/1\) loading; further results are given as reference.

For 1/1 loading the optimized geometry (Fig. 9b) indicates at the center of the notched region stress triaxialities up to 1.12, whereas the initial geometry (Fig. 9a) indicates values up to 0.8, i.e. through the geometry changes, mainly the sharper notch in thickness direction, the V1 geometry indicates substantially higher tension dominated stress triaxialities. Whereas the V2 geometry (optimized for \(-1/1\) loading) reaches only values up to 0.68, see Fig. 9c. Hence the effect of the optimization process can be clearly seen. In contrast, the influence of the optimization process is less obvious for \(-1/1\) loading (Fig. 9d–f). All three geometries show values in the range from 0.04 to 0.1 in a similar distribution throughout the cross section characterizing a shear dominated stress state.

Note that Fig. 8 indicates no significant change of the overall homogeneity of the stress triaxiality distribution for all three geometries, i.e. gradients clearly occur. However, in the central part of the V1 geometry the high intense stress values are locally homogeneous within a limited area. For more significant results in terms of an overall homogeneous distribution over the whole cross section area, the optimization problem has to be reformulated in terms of a reconsidered objective function or additionally formulated constraints.

Fig. 9
figure 9

Stress triaxiality at cross section

Full size image

The Lode parameter \(\omega \) has not been involved into the optimization process as objective value. However, it also characterizes the stress state. Under 1/1 loading (Fig. 10a–c) the influence of the geometry can be clearly noted: The initial and the V1 geometry, both with a cross sectional area of 6.0 mm by 2.0 mm, indicate values between \(-0.1\) and \({-1.0}\) whereas the V2 geometry (Fig. 10c) with a cross-sectional area of 4.0 mm by 3.0 mm indicates only values close to \({-1}\). Again, under \(-1/1\) loading the influence is less significant and a homogeneous distribution with values close to \({-1}\) is shown.

Fig. 10
figure 10

Lode parameter at cross section

Full size image

6 Experimental investigations

6.1 Experimental setup

All experiments have been performed with the biaxial test machine LFM-BIAX 20 kN produced by Walter+Bai, Switzerland. It contains four electromechanically, individually driven cylinders with load extrema of ±20 kN. During the experiments the specimens are clamped in the four heads of the cylinders and the machine reports its displacements as well as the applied forces of each cylinder. To avoid non-symmetric behavior during the experiments a stable, mainly displacement driven procedure has been used, see [20] for details. In this context the nominal displacements \(u_{i.j}\) (indicated by the red dots in Fig. 11) allow the introduction of the relative displacements \(\Delta u_{\text {ref}.i} = u_{i.1} - u_{i.2}\) as adequate displacement measures and the averaged forces \(F_i=\nicefrac {(F_{i.1}+F_{i.2})}{2}\) are introduced as corresponding force measures.

Fig. 11
figure 11

Notation of experimental technique. (Color figure online)

Full size image

During the experiments the displacement fields of the specimen surfaces have been monitored with a Q-400 digital image correlation (DIC) system provided by Limess/Dantec. For the setup 6.0 Mpx cameras equipped with 75 mm lenses have been applied. The forces \(F_{i.j}\) and the machine displacements \(u^\text {M}_{i.j}\) have been transmitted to the DIC-system and stored with the corresponding DIC data sets. The nominal displacements \(u_{i.j}\) have been extracted by post processing the DIC data. Furthermore, macro photographs of the fractured specimens have been taken and the fracture surface has been analyzed by scanning electron microscopy (SEM). Further details on the experimental techniques can be found in [22].

6.2 Forces and displacements

Figure 12 displays the numerically (sim) and experimentally (exp) obtained force–displacement curves for (a) 1/1 loading and (b) \(-1/1\) loading. The relative displacements \(\Delta u_{\text {ref}.i}\) of the numerical simulations have been extracted at corresponding nodes to the experimentally obtained ones, see red dots in Fig. 11.

Under both loading conditions the influence of the specimen geometry can be clearly seen within the elastic region: The V1 geometry indicates the stiffest reaction due to the sharp notch radius in thickness direction with \(R_t={1.0}\) mm causing a relatively small region with elevated elastic deformations. Under 1/1 loading the initial geometry (\(R_t={2.0}\,{\hbox {mm}}\)) reacts less stiff and the V2 geometry with \(R_t={1.0}\) mm and a reduced penetration depth of \(D={0.5}\) mm reacts softest. Under \(-1/1\) loading this influence is present, but less pronounced. Overall the agreement between numerical simulation and experiment within the elastic region is good.

The inelastic behavior is predicted in a reasonable way, having in mind that the material parameters have been determined with respect to the tension test and that the numerical simulations considered elastic–plastic material behavior based on the von Mises yield criterion and the isochoric inelastic deformations, i.e. the material degradation due to damage has not been included numerically. For 1/1 loading it can be noted that the simulations of the initial and V2 geometry underestimate the load whereas the simulation of the V1 geometry overestimates the load. This can be seen in relation to the increased stress triaxiality of the V1 geometry under 1/1 loading, see Fig. 9. Under \(-1/1\) loading more elevated relative displacements \(\Delta u_{\text {ref}.i}\) occur, see Fig. 12b. For the tension as well as for the compression axis the elastic behavior is estimated in good accordance. For the compression axis also the inelastic behavior is estimated in a reasonable way, whereas for the tension axis starting from a relative displacement \(\Delta u_{\text {ref}.2}={0.4}\) mm the loads are overestimated, see Fig. 12b.

Fig. 12
figure 12

Global load–displacement diagrams

Full size image

6.3 Deformation and fracture behavior

The deformation behavior of the specimens surfaces has been monitored by digital image correlation (DIC) and the reported first principal strains are displayed in Fig. 13 for the different geometries and load cases shortly before fracture occurrence. Furthermore, the photos given in Fig. 14 display the resulting fractured specimens and give a good impression of the overall deformation behavior of the central part of the specimen. It can be noted that the geometry has a significant influence on the deformation and fracture behavior, specially the sharp notch in thickness direction of the V1 geometry leads to a reduced area of elevated strains and consequently less overall deformation before fracture.

Fig. 13
figure 13

First principal strains reported by digital image correlation (DIC): load case (1/1): a initial geometry, b V1 geometry, c V2 geometry; load case \((-1/1)\): d initial geometry, e V1 geometry, f V2 geometry

Full size image
Fig. 14
figure 14

Fractured specimens

Full size image

The fracture surfaces have been analyzed by scanning electron microscopy (SEM) and representative pictures are given in Fig. 15 for 1/1 loading and in Fig. 16 for \(-1/1\) loading. The upper row (a–c) reflects a representative area at the center of the fracture surface and the lower row (d–f) reflects a representative area towards the outer radius \(R_o\).

Under 1/1 loading, the fracture surfaces indicate failure due to void nucleation, growth and coalescence which corresponds to the indicated tension dominated stress state, see Figs. 9a–c and 10a–c. For the optimized V1 geometry significantly higher stress triaxialities have been predicted at the center of the notch and this is reflected by the material failure. Figure 15b indicates significantly remarkably larger voids and a more brittle behavior compared with the initial (Fig. 15a) and the V2 geometry (Fig. 15c). Towards the outer radius for all geometries a similar stress triaxiality of approximately \(\eta =0.4\) has been predicted and the SEM images (Fig. 15d–f) indicate the corresponding material failure characterized by smaller voids without significant differences between the different geometries.

Fig. 15
figure 15

Scanning electron microscopy images load case (1/1): a initial geometry, central; b V1 geometry, central; c V2 geometry, central; d initial geometry, boundary; e V1 geometry, boundary; f V2 geometry, boundary

Full size image

For \(-1/1\) loading stress triaxialities without significant gradients have been numerically predicted (Fig. 9d–f) with stress triaxialities of approximately \(\eta =0.05\) and a Lode parameter of approximately \(\omega =-1\) (Fig. 10d–f) which corresponds to a shear dominated stress state. The SEM images of the fracture surface in Fig. 16 indicate the corresponding shear behavior initiated by micro shear cracks. Overall the fracture surfaces are rather homogeneous and no significant differences can be noted neither with respect to the geometry nor with respect to the location on the fracture surface.

Fig. 16
figure 16

Scanning electron microscopy images load case \((-1/1)\): a initial geometry, central; b V1 geometry, central; c V2 geometry, central; d initial geometry, boundary; e V1 geometry, boundary; f V2 geometry, boundary

Full size image

7 Summary and outlook

The paper has discussed the shape optimization of the biaxially loaded X0-specimen with respect to a homogeneous stress state to study the damage and fracture behavior. Experiments with the optimized specimens have been conducted and the obtained results emphasize the efficiency of the optimization process.

The proposed variational sensitivity analysis has been integrated into a gradient based optimization algorithm ensuring reasonable computing time for all considered numerical investigations. The chosen optimization problems modeling the mechanical intention of distinct and preferably homogeneous stress states have been solved requiring only a small number of iterations. Both remarks emphasize the practicability of the chosen computational approach generating quantitative results unavailable by engineering intuition.

The experiments with the initial and the two load case dependent optimized X0-specimens confirmed the numerically predicted results. Especially, the increased stress triaxiality of the V1 geometry under 1/1 loading indicated bigger micro-pores and a less ductile behavior revealed by scanning electron microscopy.

The presented promising results mark the onset of a research program of shape optimization with respect to biaxial specimens to study the damage and fracture behavior of ductile metals. In future investigations, shape optimization problems with respect to other loading condition needs to be studied and if possible standardized geometries for certain loading conditions need to be proposed. Furthermore, the damage evolution and the constraint of limiting elevated deformations have to be considered within the optimization processes.