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

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.


Introduction
In biaxial tests the stress state can be characterized by the stress triaxiality η defined by with the mean stress σ m = 1 3 I 1 (σ ) and the equivalent von Mises stress σ eq = √ 3 J 2 (devσ ). 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 Responsible Editor: Emilio Carlos Nelli Silva Jan Liedmann jan.liedmann@tu-dortmund.de 1 TU Dortmund University, August-Schmidt-Straße 8, Dortmund, 44227, Germany 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 2 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 Fig. 1 Initial geometry of X0-specimen is concerned with numerics. The method ofF 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 X0specimen in Section 6.2. Section 7 summarizes the paper and draws a conclusion on the findings. Discretization usinḡ F finite elements and the corresponding explicit discrete operators for structural as well as sensitivity analysis are given in the Appendix.

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.

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 p = [ R 1 R 2 R 3 d ] contains the design variables indicated in Fig. 2.
Here, l cs and t cs denote the length and thickness of the central cross-section, respectively. The objective function is chosen to be 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 which forces a constant cross-section areaĀ 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 3 to prevent sharp edges in the geometry and the finite element mesh. The results of the shape optimization are presented in Section 6.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 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 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.

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 Single contraction is either denoted as dot product for vectors or omitted for higher order tensors, double contraction is denoted as colon, e.g.
Special transpositions of fourth-order tensors are denoted as superscripts, e.g.
which represents the transposition of the second and third bases of a fourth-order tensor.
In some equations, the special product 21 * is used, which denotes a single contraction of the second basis of a fourthorder tensor with the first basis of a second-order tensor, i.e., The special fourth-order tensors I sym and I dev denote the symmetric and deviatoric projection tensors with the coefficients

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 B with Cartesian basis Z i and local coordinates Θ, the classical deformation mapping ϕ, defined as can be decomposed into two independent mappings, i.e., the design-dependent geometry mapping κ(Θ, s) and the time-dependent motion mapping μ(Θ, t) κ : (Θ, s) → X(Θ, s) and μ : (Θ, t) → x(Θ, t).
With the corresponding tangent mappings not only the deformation mapping but also the deformation gradient can be decomposed and written as Consequently, for volume mappings, we obtain with the determinants J = det F, J K = det K, and J M = det 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 K = K(s) and the referential points X = X(s). Note, due to the choice of the elastoplastic material model, the plastic intermediate configuration 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.

Variations and derivatives
The total variation of any quantity (•)(u; s; h n ) is given by the sum of the partial variations where, e.g., the partial variation with respect to u with fixed s andĥ n is defined as This notation is used throughout the whole paper.

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 e F p can be found in the relevant literature, e.g. Elguedj and Hughes (2014); Simo (1988a, b) and captures large deformations and strains.

Mechanical response
The nonlinear weak form of equilibrium R : V → R, commonly called residual, in the reference configuration K reads where Gradv is the gradient of the test function v ∈ V and u ∈ V is the displacement field. Here, V denotes the usual Sobolev space of all admissible displacements. The chosen strain energy function given by 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, C −1 p 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 b e = F e F T e using the isochoric deformation gradient F The von Mises-type yielding condition only depends on the deviatoric part of the stresses, and thus is pressure independent 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 α where the material parameters σ 0 and σ ∞ 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 quasistatic, 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 Evolution of the internal history variables is described by the flow rulė and the hardening laẇ 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 With the increment of the plastic multiplier at hand, the update of the deviatoric stress tensor can be determined as tr e 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 The total Kirchhoff stress tensor reads and the first Piola-Kirchhoff stress tensor is given by 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 theF method and in view of sensitivity analysis, it is convenient to use the algorithmic material tangent A = ∂P K ∂F , which can be split into a volumetric and deviatoric part The explicit forms of the volumetric and deviatoric parts read and respectively. The fourth order tensor S dev is given by with the factors and the identities The whole computational procedure is summarized in Algorithm 1 in pseudo code format.

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 where the two bilinear forms k : V × V → R and p : V × S → R represent the partial variations of the weak equilibrium condition w.r.t. displacements and design, respectively. The third bilinear form h : V × G → 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 S and 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.

Matrix form of sensitivity relations
From the matrix form of (37) where R ∈ R ndof , the sensitivity matrix S ∈ R ndof×ndv can be derived, that is, the total derivative of the structural response u w.r.t. design changes s Here, the matrix Z n ∈ R nhv×ndv is the total design derivative of the history variables δh n = Z n δ s, see Section 4.4 for details. K ∈ R ndof×ndof denotes the tangent stiffness matrix, P ∈ R ndof×ndv denotes the so-called pseudo load matrix and H ∈ R nhv×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 = 0 and δh 0 = 0 at time t = t 0 , from (39) we obtain 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.

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 where J K = det 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 the determinant of the geometry gradient the first Piola-Kirchhoff stress tensor and the deformation gradient cf. Barthold (2002), Materna (2009), and Kijanski (2018), we can write where we have used the identity dV = J K dV Θ .

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 but where M is the third-order tensor connecting the first Piola-Kirchhoff stress tensor with the material parameters. Thus, we can write For reasons of brevity, without detailed derivation, for the vector of material parameters m = [E ν σ 0 σ ∞ H β], the explicit form of the tensor M reads Here, the second-order tensor H m reads The partial derivatives of the bulk modulus K, the shear modulus G, and the nonlinear hardening function k are given by with with c 1 = 1 − exp(−β α), c 2 = σ ∞ α exp(−β α), respectively.

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 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 δ h n P K = ∂P K ∂h n Z n δs.
The partial derivatives of the stress tensor w.r.t. the history variables h n = { C −1 p,n , α n } are given by with and the tensor H from (32) for the internal variable C −1 p and for the internal variable α, with the first derivatives of (22) and (21) and

Update of history variations
As (39) 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
with δ u F = Gradδu. The yet unknown partial variations ∂h ∂F and ∂h ∂h n can be resolved to for the internal variable C −1 p w.r.t. the deformation tensor F and for the internal variable α.
Constitutive parameters For material parameters as design variables, (63) takes the form The partial derivatives of the internal history variables w.r.t. the material parameters ∂h ∂m can be identified to for the internal variable C −1 p , where the second-order tensor H m has been used, cf. (52). For the internal variable α, we find 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 pseudotime step are needed and read and for the internal variable C −1 p . The fourth-order tensor B −1 C is given by For the internal variable α we find with the fourth-order tensor B C from (59), and with the derivatives given in (61) and (62).

Stress sensitivity
The gradient of the stress triaxiality as defined in (1) can be computed following the same concepts, which leads to where the total variations of the mean and von Mises equivalent stress read Thus, it is essential to calculate the total variation of the Cauchy stress tensor given by

5.1F 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), anF formulation is used. This formulation has advantages on implementation side as any constitutive law can be used. TheF deformation gradient reads where J 0 = det 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 F =F. Therefore, for exact linearization, the derivatives of F w.r.t. F and F 0 are needed. The total variation ofF reads

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 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.
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.
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 X0specimen 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 −52 represents the distance of 1.0 to the next larger double-precision number.

Computational effort
For each integration point in each element, the vector of history variables h ∈ R 1×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 Z ∈ R nhv×ndv w.r.t. design s ∈ R ndv×1 have to be saved. The structure of the global history field 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 nh = nel × ngp × nhv × (1 + ndv) 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).

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. The geometry is discretized with 1539 hexahedralF 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 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 @Obj 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 With this at hand, the next optimization task, that is the shape optimization of a chosen specimen, can be conducted.

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 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 1 R 2 R 3 d], namely R 1 : inner radius, R 3 : radius in thickness direction, R 2 : outer radius, d : penetration depth.
Aim of the shape optimization is to achieve a preferably homogeneous distribution of stress triaxiality in the crosssection. 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 eq = t cs · l cs − 12 mm 2 = 0, cf. Fig. 2. To prevent the appearance of sharp corners at the edges of the notch in thickness direction, the penetration depth (a) Initial geometry (b) Optimized geometry  should not be larger than the radius in thickness direction. This is captured by the inequality constraint Consequently, the optimization problem is defined as 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 @Obj and @Geomconstr 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 eq = 3.2925e − 10 and c in = − 8.9378e − 4.

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 loaddisplacement 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.
Funding Information Open Access funding provided by Projekt DEAL.

Compliance with Ethical Standards
Conflict of interest The authors declare that they have no conflict of interest.

Replication of results
The described method is outlined in all significant steps naming the continuous formulation as well as resulting discrete equations. Furthermore, most important matrices are stated and important algorithmic features are described showing the corresponding Matlab code.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creativecommonshorg/licenses/by/4.0/.

A.1 Discrete gradient and divergence operators
For an 8-node hexahedral finite element, the discrete gradient operator at each FE node i reads as well as discrete vector forms of non-symmetric secondorder tensors P = [P 11 P 21 P 31 P 12 . . . P 33 ] T ∈ R 9×1 .
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 with The discrete divergence operator reads With the matrix forms D and D 0 of D and D 0 from (80) using the convention stated in (89), we can define the two gradient operators for theF element

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 stiffness matrix
Element geometric pseudo load matrix 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 P e p,ik = P e X,ij Q e jk , 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)