Computational shape optimisation for a gradient-enhanced continuum damage model

An isotropic gradient-enhanced damage model is applied to shape optimisation in order to establish a computational optimal design framework in view of optimal damage distributions. The model is derived from a free Helmholtz energy density enriched by the damage gradient contribution. The Karush–Kuhn–Tucker conditions are solved on a global finite element level by means of a Fischer–Burmeister function. This approach eliminates the necessity of introducing a local variable, leaving only the global set of equations to be iteratively solved. The necessary steps for the numerical implementation in the sense of the finite element method are established. The underlying theory as well as the algorithmic treatment of shape optimisation are derived in the context of a variational framework. Based on a particular finite deformation constitutive model, representative numerical examples are discussed with a focus on and application to damage optimised designs.


Introduction
The application of finite element simulations has become a well established framework to circumvent costly experiments in the industrial environment and to computationally predict material and structural behaviour. The accurate prediction of the evolution of deformation and material properties for different fields of application, ranging from e.g. simple forming processes to complex crash tests, requires advanced models capable of simulating inelastic effects of the particular material considered. In view of industrial applications, damage evolution and accumulation play a crucial role. To this end, the framework of continuum damage mechanics is adopted, where damage is defined as the degradation of effective material properties which evolve after the initiation of defects, such as microcracks, see e.g. [30]. Different damage models have been proposed in the literature and the interested reader is referred to, e.g., [30] for a detailed overview of different approaches. This work adopts the well-established [1 − d] approach in the context of a Lemaitre-type damage model which captures the coupling of elasticity and damage, see [27].
A common problem of all damage models, if embedded into the finite element method, is the mesh dependency of the simulation response together with a possible loss of ellipticity of the underlying governing equations. Various non-local damage formulations have been established to obtain regularised, mesh independent results. In [5,6] a non-local damage variable is defined as an integral value of a pointwise defined damage quantity, while in [10,40] regularisation is obtained by introducing additional gradient terms. Alternatively, regularisation can be achieved by means of rate dependent damage evolution equations, see e.g. [14] or the more recent works in [2,26]. Moreover, different gradient-enhancement strategies have been established in the literature. A popular approach introduces a non-local quantity, e.g. non-local damage, coupled to a local quantity and only takes the gradient of the non-local quantity into account as proposed in [33], later adopted in [11,38] and denoted micromorphic approach in [17]. The approach has been successfully extended to anisotropic damage formulations, see e.g. [41], and coupled to plasticity [8,12,24] to name but a few references. Implementation of this concept into commercial FEM software, e.g. Abaqus, by exploiting the similarity to the heat equation is presented in [32]. An alternative approach directly incorporates the gradient of the local quantity, see e.g. [34] introducing a quasi non-local energy release rate. Thereby, in the context of a finite element implementation, the underlying equations have to be solved at the global level, e.g. with an active-set strategy, see [28]. The drawback of using an active-set strategy is a non-constant size of the global system of equations. Alternatively, the underlying set of equations can be solved by means of a Fischer-Burmeister approach, see [16]. The Fischer-Burmeister approach has been proposed in [37] as a practical method of solving the Karush-Kuhn-Tucker (KKT) conditions in material modelling and was later applied to different modelling problems, see e.g. [3,9,13,23,24].
The application of optimisation techniques to damage mechanics can be traced in the literature in different forms of application. In [1] a numerical implementation of the Francfort-Marigo model of damage evolution in brittle materials is proposed. At each time step, the sum of an elastic energy and a Griffith-type dissipated energy is minimised. The interface between healthy and damaged phases is modelled by a level set function that is transported according to the shape derivative of the minimised total energy. Initially, the damaged zone is nucleated by using the so-called topological derivative. The overall algorithm is able to predict crack propagation, including kinking and branching. Furthermore, optimisation algorithms are applied to identify damage in an existing structure. The research focuses for example on the algorithmic aspects of nonlinear programming approaches and the formulation of the optimisation problem, cf. [35]. Another investigation highlights the benefits of sensitivity-based structural damage detection by using the matrix spectral decomposition theory, cf. [44].
While shape optimisation is generally broadly applied, its application to optimise damage accumulation and evolution is largely unexplored. Most construction parts contain stress concentrations at some changes in cross section or boreholes. One possible reason for failure are fatigue phenomena, especially in the case of dynamic loading. Fatigue is observed at load levels below the macroscopic yield stress. On a macroscopic scale the structural behaviour is elastic, but on a microscale plasticity is found at inclusions and other types of defects within the material. Based on a continuum damage mechanics model, shape optimisation minimising a cost function modelling fatigue has been successfully applied by [19]. An engineering definition of lifetime based on load-cycle counts has been maximised and verified by experiments. In [4] a local damage model for finite strains was used to optimise geometries in order to reduce damage. Thereby, the focus is set on the treatment of internal variables within the context of a variational approach of the sensitivities.
Similarly to shape optimisation, advancements of topology optimisation regarding damage are limited, while other nonlinear material effects, such as finite strain plasticity or viscoplasticity, c.f. [21,43], see increased research in the field of topology optmisation. A continuum micromechanics based approach for assessing changes in elastic properties and in the strength of the structure attributed to chemical reactions is proposed in [36]. It is found that the damage progression rate is mostly dependent on relative rates of sulfate ingress and calcium leaching. The relations between the factors and the damage progression rate are observed to be significantly nonlinear. A level set based topology optimisation approach to design structures exhibiting resistance to damage is outlined in [31]. As the damage process is irreversible, the structural responses are path-dependent and this dependency is accounted for in the sensitivity analysis. In [39] topology optimisation is applied to high-cycle fatigue. By constraining the maximum damage at each element via the p-norm of the damage values, improved topologies are generated by utilising sensitivities derived using an adjoint method.
The paper is structured as follows: A brief overview of the nonlinear kinematics is addressed in Sect. 2 with regard to the extension concerning the intrinsic formulation. Thereafter, the material model is discussed. The implementation of the model in the context of the finite element method is elaborated in Sect. 3 and the focus is set on the implementation of the Fischer-Burmeister approach. Sect. 4 deals with the optimisation approach and gives a brief overview of shape optimisation, the necessary derivations and the implementation for the numerical framework. To show the possibilities of the model in combination with the optimisation framework, Sect. 5 discusses representative numerical examples as well as optimised geometries with regards to reduced damage accumulation. Finally, Sect. 6 gives a summary and an outlook on possible future research work.

Material model
This section presents the underlying theories for the damage model. The nonlinear kinematics for large deformations are briefly described and enhanced for later use regarding the calculation of the sensitivities. This enhancement, together with the introduction of the intrinsic formulation allows the decoupling of the quantities from the design and, as such, easier derivation. The material model for the non-local, gradient enhanced damage will then be derived.

Kinematics
The kinematics are presented here for later application, and to give a brief overview to the reader regarding the notation. Since only a broad overview is given, the reader is referred to [7] for the general kinematics, as well as [4,29] regarding the intrinsic formulation.
We consider a body in an undeformed, referential state with configuration B 0 ⊂ R m with referential position vectors of material points X ∈ B 0 . Any of these referential position vectors X can be mapped to position vectors x ∈ B t in the deformed, current configuration B t ⊂ R m at time t ∈ [0, T ] by the nonlinear mapping ϕ, see Fig. 1a, i.e. ϕ : For any differentiable function (•) the gradient and divergence operations with respect to the reference configuration read whereas the operations with respect to the current configuration read Based on this, the deformation gradient is then expressed by with u = x − X as the displacement vector field. The deformation gradient transforms an infinitesimal line element dX of the referential configuration, to the line element dx in the current configuration. The corresponding Jacobian J maps an infinitesimal volume element dV of the undeformed configuration to the deformed volume element dv, i.e.
Furthermore, we introduce the right Cauchy-Green tensor as For the derivation of the quantities necessary for the structural optimisation the introduction of the intrinsic formulation is helpful, which enhances the above setting by the domain B ζ ⊂ R n as illustrated in Fig. 1b. Therein, a new coordinate Θ ∈ B ζ is introduced, which allows the (sufficiently smooth) mapping to the parameter space in order to calculate the derivatives without any dependencies of the referential configuration. This leads to the additional operations GRAD(•) = ∇ Θ (•) and DIV(•) = ∇ Θ · (•).
The additional mappings read κ : for the nonlinear mapping from the parameter space to the reference configuration, as well as for the mapping from the parameter space to the current configuration. Following the derivation of the deformation gradient, the counterparts for the parameter space, as well as their Jacobians J • , are defined as respectively M := ∇ Θ μ and J M = det(M).
Combining the two mappings into a multiplicative decomposition of the deformation gradient simplifies the derivation of the total variation for the deformation gradient. For further detailed background information the reader is referred to the references cited above.

Constitutive relations
In this section the material model, combining isotropic elasticity with isotropic damage, is derived. For the gradientenhanced model a Helmholtz free energy is proposed as where d is the scalar variable associated to isotropic damage. The elastic energy is additively split into a volumetric ψ vol and an isochoric ψ iso contribution. Gradient enhancement is captured in the term ψ grd . The damage variable d deteriorates the elastic properties of the material model via so-called damage functions f • , i.e. The such that no additional constraints to limit the value of d are necessary. Following the framework of generalised standard dissipative materials together with the framework of the introduction of a non-locality residual, the Clausius-Duhem inequality reads with the Piola stresses P. Thereby, P denotes the nonlocality residual, see [15,22,28,34]. Moreover, the notatioṅ • represents the material time derivative of the quantity •.
The non-locality residual accounts for the energy exchange between the particles in the damaged area B d and has to satisfy the insulation condition which states that no energy is exchanged between particles inside and outside the damaged area. In the elastic domain B e := B\B d the non-locality residual has to be pointwise zero, i.e.
With the definition of the stresses and driving forces the dissipation inequality (16) reduces to Since d and ∇ X d cannot evolve independently from each other, the reduced dissipation inequality is assumed to take the bilinear form Comparing both forms of the reduced dissipation, (20) and (21), yields an expression for the non-locality residual in the form of Combining the insulation condition (17) with the vanishing non-locality residual in the elastic domain (18) renders the integral of P over the whole domain B to vanish. Additionally applying integration by parts and the Gauss theorem leads to with N being the referential outward unit normal vector. Since (23) has to be fulfilled for arbitrary damage ratesḋ one obtains the following conditions The quantityȲ can be interpreted as a non-local driving force. Consequently, damage potential φ d -governing onset and evolution of the damage variable d in this associated framework-is formulated in terms ofȲ such that the elastic domain is introduced as Homogeneous Neumann boundary conditions are applied to the total boundary and no Dirichlet boundary conditions are considered, see (25). The constrained minimisation problem ensuing from the postulate of maximum dissipation is solved with the Lagrange functional, i.e.
with Lagrange parameter λ ≥ 0. This leads to the Karush-Kuhn-Tucker (KKT) conditions. The evolution equation for the damage variable is obtained aṡ and the loading/unloading conditions are In this work, the elastic free Helmholtz energy is chosen to be of Neo-Hookean type, whereas the gradient part is a simple quadratic form, i.e.
where K is the bulk modulus, μ the shear modulus, c d shall be denoted as regularisation parameter and • = √ • · •.
Moreover, C iso = J − 2 3 C is the isochoric part of the right Cauchy-Green tensor C. The conditions for the damage function f • are fulfilled by the exponential function. We choose where η vol and η iso are material parameters controlling the speed of deterioration of the elastic properties. The stresses and driving forces are then given by Finally, damage potential φ d is chosen as so that damage driving forceȲ is compared to a constant threshold value y 0 . This results in the evolution equatioṅ cf. (28). The quantity d is introduced as a field variable in the algorithmic setting and the Lagrange multiplier λ is solved for by the implicit Backward-Euler scheme, which leads to the discrete update at time step n + 1 with the incremental time step Δt = t n+1 − t n . This means that the loading/unloading conditions cannot be fulfilled by an evolution of the damage variable d at integration point level within the finite element formulation discussed as this work proceeds-instead, the loading/unloading conditions need to be solved globally.

Finite element implementation
The previously described material model shall now be incorporated into a boundary value problem. Therefore, the strong and weak forms of the mechanical problem are briefly summarised and the governing equation for the damage field is derived. Afterwards, the coupled problem is discretised and linearised in the context of the finite element method.

Weak formulation
In order to derive the weak form of the quasi-static mechanical problem, the strong form is used as the starting point, i.e.
Here, B denotes the body forces, ρ 0 the referential mass density, and T the traction forces applied on the surface. Multiplication of (36) with test function η ϕ and integration over the body B 0 , as well as the application of Gauss's theorem and integration by parts, leads to the weak form of the mechanical problem which is used later on for the finite element implementation. The damage and field variable d is governed by the KKT-conditions (29). In order to formulate a single governing equation for the whole domain, without splitting the domain into a damage domain and an elastic domain [28], the problem is restated with the help of Fischer-Burmeister functions [3,9,16,37], i.e.
so that the inequality constraints are transformed to an equation. However, solving (39) in weak form-in analogy to the mechanical problem-requires the evaluation of the second gradient of the damage field d in the non-local driving forceȲ (24). Instead, damage potential φ d , cf. (33,24), and Lagrange multiplier λ are multiplied with test functions η φ and η λ respectively and integrated over the domain B 0 , identical to the approach in [28]. Inserting (33) and using integration by parts and Gauss's theorem results in with the last integral of (40) vanishing due to the boundary condition (25). The second derivative of damage field d, as present within the expression ∇ X · Y , is no longer explicitly required. This framework is only applicable if the damage potential φ d is linear in the non-local driving forceȲ . The global form of the loading/unloading conditions (40) and (41) are now enforced by means of the Fischer-Burmeister complementarity function which results in Finally, this leads to the coupled problem, the residual form of which takes the representation

Discretisation
For the finite element implementation, domain B 0 is discretised by n el finite elements and n np nodes, i.e.
where each finite element is characterised by n ϕ en nodes with three displacement degrees of freedom and n d en nodes with one non-local damage degree of freedom. Applying the isoparametric concept, the geometry description X, the placements ϕ and the damage field d are approximated by using the nodal values of the field and the shape functions N • , so that one obtains and the values for their respective gradients as According to the Bubnov-Galerkin approach, the above approximations are also applied to the variations, respectively test functions, as well as their gradients, i.e. and Inserting the above approximations into (38) yields the mechanical residual in element e connected to element node I , where the first integral can be identified as the internal forces f int eI and the sum of the latter integrals as the external forces f ext eI . Since the (49) has to be fulfilled for any variation η ϕ I , the residual can be written as The same discretisation is applied to the global forms of the loading/unloading conditions (40) and (41), i.e.
where the test functions are already removed. The element residuals r ϕ eI , r φ eK and r λ eK are now assembled into their respective global finite element residuals, i.e.
After assembling, residuals r φ and r λ are used to compute the Fischer-Burmeister residual for every global node point to

Linearisation
To solve the nonlinear system of equations an iterative Newton-Raphson scheme is employed. A Taylor series expansion at iteration step m-terms of order higher than linear are neglected-yields where the nodal increments of the residual are given by with Δϕ j = ϕ j m+1 − ϕ j m and Δd l = d l m+1 − d l m being the increments of the nodal degrees of freedom. Defining the global stiffness matrix as enables the entries connected to the total derivative of the mechanical residual, K ϕϕ and K ϕd , to be computed via the assembly of the element stiffness matrices, i.e.
with the element stiffness matrices given by with I being the identity tensor of second order. The entries in the stiffness matrix connected to the damage residual need to be treated in a different manner due to the modification of the assembled residual in (54). Using the chain rule, the incremental form of the damage residual takes the representation Setting up diagonal matrices for the partial derivatives of the Fischer-Burmeister function in the form of and computing the total derivatives from the assembly of the element stiffness matrices, i.e.
allows the specification of the global stiffness entries K dϕ and K dd as The element stiffness matrices K φϕ eK J , K φd eK L and K λd eK L can be computed per element by The total derivatives of the stresses and driving forces are given in Appendix A. Finally, this leads to the system of linear equations to be solved in every iteration of the Newton-Raphson scheme Due to the application of the Fischer-Burmeister function for the nodal residual of the damage evolution, the above system of equations remains of constant size-in contrast to other solution methods, like the active-set method, where this is not always the case. This is due to the fact, that the Fischer-Burmeister framework does not explicitly distinguish between "active" and "inactive" constraints of the KKT-conditions (29) but rather transforms the inequality constraints to equations.

Remark 1
Due to the Fischer-Burmeister approach solving the underlying set of equations on global finite element level, the nodal residuals have to be calculated before the final stiffness matrix can be correctly assembled, see (69). This has to be taken into account especially when integrating this model into established FEM-programs, since modifications affecting program parts beside element and material model formulation are necessary.

Optimisation
In this section the focus is set on the shape optimisation including the previously described material model. The sensitivities will be derived by making use of the enhanced kinematics presented in Sect. 2.1. Finally, the sensitivities are used for the objective functions and constraints to allow shape optimisation to reduce and control damage for optimal designs with reduced damage accumulation.

Sensitivities
In the sense of the structural optimisation, not only are the field variables ϕ and d allowed to change, but the designthe initial geometry X-is also allowed to change over the course of the optimisation. This expansion can be pictured mathematically, i.e. the residual term in (43) has to always be fulfilled, even if a change in geometry is determined. This leads to the requirement because a change in any of the three variables ϕ, d and X has to yield a physical admissible solution R = 0. The additional variations δ d R and δ X R require the enhancements made in Sect. 2.1. The specific expressions for the additional variations are given for the mechanical residual by where the body forces and traction forces are omitted. The additional variations of the damage residual are given by where the variation of the KKT-condition (40) is obtained as and the variation of condition (41) as Applying the finite element discretisation introduced in Sect. 3.2, leads to the following nodal quantities of the pseudo-load matrix. Like the stiffness matrix, the pseudoload matrix can be split into two parts, i.e the mechanical part where H = ∇ X u, and the damage parts The special dyadic product ⊗ is defined as with A as any arbitrary tensor of second order and a and b as arbitrary vectors.
In the same fashion as the stiffness matrix, the pseudo-load matrices can be assembled, i.e.
where the damage quantities are then multiplied with the Fischer-Burmeister matrices, see (64), to obtain the complete damage part of the pseudo-load matrix, i.e.
The the global pseudo-load matrix is then defined as where n = [n ϕ dim +1] n np is equal to the dimension of the stiffness matrix K and m = n ϕ dim = n X dim . Since the derivatives of the Fischer-Burmeister function are required, the respective residuals have to be calculated and enter the calculation of the pseudo-load matrix as well. In addition, Remark 1 regarding the numerical implementation of the stiffness matrix K also applies here for the calculation of the pseudo-load matrix P.
Variation of (74) in matrix representation can then be summarised as Rearranging the above equation allows the introduction of the sensitivity matrix S δϕ δd which describes the response of the field variables ϕ and d at each node, when a change in the referential geometry is made.

Remark 2 For the optimisation, Computer Aided Geometric
Design (CAGD) is used to define the shapes for the optimisation. Henceforth, the control points of those shapes are used as the design variables s, which can be incorporated into the sensitivity matrix by means of the chain rule, i.e.
With this, Bézier-surfaces are defined to generate the shape for the optimisation problem in order to allow shape optimisation.
Here, the direct differentiation sensitivity analysis is conducted. One could also choose the adjoint method which is generally used when problems involve many design variables and few constraints are applied, cf. [20,25]. In this paper however, the problem is of opposite type, i.e. few design variables in comparison to a large number of constraints, e.g. (100).

Objective functions and constraints
For the optimisation, objective functions and constraints are necessary to generate geometries with specific properties. Since the central aim of this work is to generate structures with reduced damage, damage variable d is of main importance with regard to the minimisation of damage. The reduction of damage can be incorporated either by direct minimisation of the overall accumulated damage or by utilising the nodal damage values as a side constraint.
Starting with the direct minimisation, an objective function of scalar value is chosen. To be specific, the 2-norm of the damage vector d-d is the vector with all nodal damage values, see (73) for comparison-is chosen in order to generate the scalar valued function The derivative with respect to the design variables can then be calculated by utilising the chain rule as well as the sensitivity matrix S, to be specific For the indirect damage minimisation, damage variable d is restricted at every FEM-node by the constraint The derivatives for these constraints are stored in the damage part S d of the sensitivity matrix.
To incorporate the indirect damage minimisation, any suitable objective function may be chosen, e.g. the compliance, or compliance-like quantity since typically the factor 2 is multiplied in the linear elastic case, cf. [18]. It is equivalent to the negative energy, such that where the external part is neglected, since no external forces are applied to the example problems. Other definitions of compliance exist, see e.g. [42], which gives a brief overview of additional compliance definitions and possible alternatives. Minimisation of the compliance is equivalent to maximising the stiffness of a given design and as such leads to a stiffer structural response. The derivative only requires known quantities and the discretised form results in which is assembled in the standard manner within the FEM.
To additionally control the design during the optimisation, the volume is constrained to allow comparison between the different designs, which is of special interest regarding the compliance minimisation, since the initial volume of the reference design shall not be exceeded.
The volume is calculated by the sum of the element volumes leading to the discretised derivative With such objective functions at hand, the general minimisation problem can be represented as

Examples
First, the capabilities of the gradient enhanced model regarding the mesh-independence are presented. Afterwards, the given example of a plate with a hole is optimised with the aim of damage reduction. By choosing the damage as the objective function, or as the constraint in the optimisation of another objective function, the damage accumulation in the optimised structure can be reduced. Finally, the generated geometries are tested by 3d-printing the structures to generate experimental results of the optimised geometries so as to (qualitatively) validate the results of the simulations.

Regularisation behaviour
For all results presented in the following, a plate with a hole is used as the representative example, with its design leading to an inhomogeneous deformation and stress state. The plate is depicted in Fig. 2. Width w and height h are set to 10 mm, while the thickness t is set to 1 mm. The radius of the inner hole is equal to r = 2.5 mm. For the load, a prescribed displacement of u pre = 0.5 mm is applied at the top and bottom of the specimen. Due to its symmetric properties, it is sufficient to reduce the simulation to the meshed part in Fig. 2, by applying appropriate boundary conditions. The parameters used for the model in the mesh independence study are presented in Table 1.
To analyse the mesh-independent features of the model, three different mesh sizes with 1350, 3750, and 7350 hexahedral elements with linear shape functions (displacements and damage) are used. For the initial comparison, gradient parameter c d is kept at the value denoted in Table 1. The results are presented in Figs. 3 and 4a. The different meshes show a nearly identical response behaviour, demonstrating the regularising effect on the damage evolution. Figures 4b and 5 demonstrate the influence of the gradient parameter c d on the material behaviour, as it governs the damage regularisation. As expected, lowering the parameter value leads to a more localised evolution, while increasing the parameter value leads to a wider band of damage evolution.
For this set of material parameters and for the design of the underlying problem, small values of the gradient parameter (c d < 1) lead to localising models, which cannot be solved by means of a standard Newton-Raphson scheme but which need the application of solvers such as the arc-length method to follow the solution path. Increasing the c d value leads to a stable solution process and displays the impact of the parameter regarding the regularisation behaviour but also leads to wider damage bands, which may be inappropriate for the desired material to be modelled.

Damage minimisation
So far, only the structural response and the capabilities of the damage model have been analysed. Now, the same exemplary plate problem is used for optimisation purposes. For the optimisation, the initial design has to be enhanced by CAGD in order to automatically generate new geometries during the optimisation. The initial design, along with its structural response and the CAGD, are presented in Fig. 6. A 3d model is used for the optimisation as well, however the thickness is omitted for optimisation purposes and remains constant. Since the results are compared to physical 3d-prints in Sect. 5.2.3, the material parameters are adapted for a better match of the behaviour of PLA (polylactide), see Table 2. Since the model can at this stage only capture damage effects, rather than an additional combination with plasticity better capturing the behaviour of the particular material considered, the response of the reference sample could not be matched perfectly. As such, the parameters were chosen such that the final reaction force behaviour is similar to the experiments within the region of steady damage evolution. The initial elastic stiffness does therefore not perfectly match with the experimental results, see Fig. 13a in the validation section. The optimisation is conducted with a sequential quadratic programming (SQP) algorithm.
Two different minimisation problems are presented and compared in this section. For the first optimisation, the overall damage is minimised with a constraint regarding the volume, which shall not exceed the initial volume of  In the second optimisation, the problem is enhanced by additionally constraining the maximum damage at each node, i.e. the damage value shall not exceed a critical damage value, here chosen as d crit = 0.8, which is lower than the maximum damage value of d max = 1.140, see Fig. 7, of the first optimisation. This allows the additional indirect minimisation for each node, while simultaneously reducing the overall damage evolution in the given structure. In addition, this value can also be interpreted as a critical damage value which may lead to failure in a physical specimen. The related optimisa-  Table 2 Material parameters for the optimisation The optimised designs for the two optimisations are presented in Fig. 7. The overall reduction in the damage accumulation, compared to the damage state in the reference design, is directly distinguishable. The general shape of the two optimised designs are very similar at the top and the right region of the plate. The hole in the centre is generally enlarged and a curvature on the right border is generated. The enlargement of the centre hole by thinning the upper part of the plate is simply restricted by the initial CAGD and the allowed design space for the control points. The main difference between the two designs lies in the lower part of the hole. While the first design leads to a kink in the lower left area, the second design generates a smoother shape which leads to the desired fulfilment of the damage constraint. The first design reduces the maximum damage from the initial value of d ref max = 1.741 to d max = 1.140 which can be further reduced to the chosen critical value of d max = d crit = 0.8 in the second design.
The force-displacement curves of the new designs are presented in Fig. 8a. While the two new designs show a softer initial response, the damage evolution is delayed and leads to an ultimately higher stiffness due to the reduced overall damage after complete loading. The additional constraint for the second design, as based on (100), has no visible impact on the final stiffness which can be seen in Fig. 8b, where the value of the objective function with respect to the number of iterations is shown. The objective functions of both designs are more than halved with respect to the initial value of J init = 27.37 to J min = 10.75 in the first design and J min = 10.99 in the second design. Due to the small difference in the values both objective functions take in the respective optimised states, the structural responses of both optimised structures are very similar. In addition, due to the additional damage constraint, the solution of the second optimisation problem is found after only seven steps, whereas the first problem requires 23 steps for the optimal design.

Compliance minimisation
The second type of optimisation is the damage controlled minimisation. Here, the compliance is chosen as the objective function in order to generate stiffer structures. For comparison, the first optimisation is carried out without any damage constraints and only the volume is constrained in order to allow the comparison of the different designs generated with compliance optimisation, i.e. minimise s l ≤s≤s u C(s) The second problem is then enhanced to take the damage reduction into account, i.e. the damage is constrained at each node point, as it was done for the second design in the previous section. In order to compare the results, the volume constraint is applied as well, which results in minimise s l ≤s≤s u C(s) (102) Figure 9 shows the generated optimal designs. Due to the different constraints for compliance optimisation, i.e. without and with constraining damage, the changes in design are dras- tically different from one another. The first problem results in a shape which allows a higher load bearing capacity in the direction of the applied load, due to an increase of the cross section in loading direction and a reduction of the volume perpendicularly to it. This leads to a stretched hole in the centre with the right side of the plate remaining almost flat. Even though the damage is not directly controlled in this problem, the maximum damage is reduced to a value of d max = 1.542 compared to the initial value of d ref max = 1.741. Furthermore, the norm of damage reaches a value of d = 25.491 which is a little lower than the initial value of d = 27.37. The second design, however, leads to a solution that is comparable to the second one of the direct damage minimisation in Fig. 7. The right side shows a similar curvature, but the hole in the centre is smaller than the one in the direct minimisation. In addition, while the right side of the plate is similar to the damage minimisation, the complete section is shifted to the left and shows a bulkier design and a larger cross section in loading direction. The norm of the damage is almost half of the initial value with d = 14.395.
A look at the graphs in Fig. 10a shows the desired increase in stiffness. Both problems generate an overall stiffer response than the initial value. Interestingly enough, the second design yields a lower stiffness during elastic loading compared to the first design, but after initiation of damage a higher stiffness remains. The green curve also displays the behaviour of the direct damage minimisation curves, where the damage evolution after initiation is reduced, yielding the stiffer final response. Since the mathematical optimisation is gradient based, the solver may only find a local solution. Due to the additional constraint in the nodal damage values, the second problem is artificially pushed to a more suitable local minimum which leads to the stiffer response after complete loading. This is also visible in the evolution of compliance over the iterations in Fig. 10b, where the second problem finally results in a lower value of C = −19.5016 J than C = −18.6253 J for the first compliance optimisation problem.

Remark 3
As mentioned above, the result of the optimisation according to (101) unexpectedly results in a (local) minimum with higher compliance than the additionally constrained problem (102). To further investigate this, the optimisation problem (101) is recalculated with an initial design chosen identical to the result of optimisation problem (102). One would expect the new optimal compliance to be lower than the compliance of either of the two previously achieved designs. However, using the same algorithm, i.e. SQP, unexpectedly results in the same solution as the original design problem, namely optimisation problem (101) with reference geometry as initial design. Other algorithms can be applied, e.g. based on the Matlab fmincon toolbox. To be specific, the active-set method converges into a result with a lower compliance, as expected, see Figs. 15, 16 and 17 in the Appendix. The validation discussed in Sect. 5.2.3, however, considers the original result of optimisation problem (101). Further- Fig. 11 3d-printed specimens. The geometries from left to right: reference, damage optimised with volume and damage constraint, compliance optimised with volume constraint, and compliance optimised with volume and damage constraint more, the detailed investigation of the effect of different optimisation algorithms on the solution exceeds the scope of this paper but constitutes future research in the context of continuum damage mechanics.

Validation
With the optimised geometries at hand, physical specimens are created by means of 3d printing with PLA. Thereby, the simulated response shall be validated with experimental data. A photo of the specimens is displayed in Fig. 11. The physically printed specimens match the numerically determined geometry, but have a longer section attached to the discretised domain in order to clamp the specimens in the micro-tension device, see Fig. 12. The specimens are loaded by displacement controlled monotonous tension. Due to the limitations of the material model considered at this stage, the accurate behaviour of PLA cannot be reproduced exactly in the simulations. The material parameters used during the optimisation are manually fitted by using the experimental data of the reference geometry in such a way that the point of damage initiation and evolution of damage can be predicted, see Fig. 13a, see the previous Table 2. The experimental specimens are generated from the optimised CAGD by exporting stl-files-a file-format often used in 3dprinting-which are then used to generate the data for the 3d printer. During optimisation, the maximum prescribed displacement in the simulation is u pre = 0.5 mm. However, the samples in the experiments are loaded until failure. In order to enable a proper comparison, the simulations are continued until a prescribed displacement of u pre = 1.5 mm is reached.
The results for the comparison of the damage optimised sample is presented in Fig. 13b. The validation of the geometries obtained from the damage minimisation focuses solely on the geometry resulting from the optimisation with constrained damage, since both shapes of the optimised geometries from damage minimisation only differ slightly. Qualitatively, the response of the simulation shows a good agreement with the experimentally obtained result. Especially in comparison to the reference geometry, the same trend in the force-displacement response can be observed. The stiffness is reduced slightly, while the point of peak force is reached at higher (displacement) loading level. In addition, the force peak becomes less sharp, albeit not as much as in the experiments. The small force plateau observed in the experimental data can most likely be attributed to plastic effects which, however, are not considered in the simulation.
The compliance optimised results in Fig. 14 show the desired increased stiffness. The experimental data of the test with the geometry corresponding to the compliance optimisation without the damage constraint in Fig. 14a displays an even higher stiffness than predicted by the simulation. However, while the maximum reaction forces are drastically increased, damage seemingly starts to evolve very quickly leading to fracture shortly after damage initiation. This cannot be directly predicted by the model, where only a steady loss of stiffness can be modelled. Concerning the tendency of the force-displacement response in comparison with the reference geometry, a good agreement between experiments and simulation for the geometry with optimised compliance subject to a damage constraint is obtained, see Fig. 14b. The initial stiffness of reference and optimised geometry nearly coincide and the maximum reaction forces also reach a similar value. However, the plateau observed in the experimental response cannot be predicted with this model. Overall, the comparison between reference and optimised geometries reveals-with limited applicability of the material model towards PLA-a good validation of the optimisation response. The analysed optimised geometries all confirm the tendencies of the simulations.

Summary
In this work, an isotropic non-local gradient enhanced damage model for large deformations was derived and used to generate optimised geometries in order to reduce damage accumulation. The damage formulation is regularised by adding an additional term to the free energy incorporating the gradient of the damage field variable. By introducing a nonlocality residual and applying the insulation condition, an enhanced damage driving force is generated and used to determine damage evolution. The corresponding Karush-Kuhn-Tucker conditions were solved globally by means of a Fischer-Burmeister approach. The overall computational framework avoids the necessity of additional internal damage variables and local iteration schemes to determine their evolution in time. Deducing the final FEM problem, however, leads to some additional assembly procedures after the standard FEM assembly, which has to be accounted for. By means of a variational approach, utilising the concept of an intrinsic formulation, the analytical gradients for the mathematical optimisation were derived, and the numerical aspects presented. The additional analytical work allows a faster computation time, bypassing numerical gradients.
The applicability of the model was presented for the inhomogeneous plate with a hole boundary value problem for which the desired regularisation abilities were shown yielding identical results for different mesh sizes. Application to shape optimisation allowed for new designs with optimised damage properties. Direct minimisation of damage, by choosing a scalar valued objective function that incorporates damage, allowed for direct reduction of damage in the structure. This led to new geometric designs with less damage and stiffer structural responses under maximum load. Using the nodal damage values as a side constraint in the direct minimisation of damage or by using another objective function such as the compliance, allowed further enhancement of, and more control within, the optimisation problem. To validate the optimised designs, specimens were 3d-printed and subjected to tensile loading. The experimental results were compared to simulations. Even though the exact behaviour of the 3d-printed material was not reproduced by the considered material model, the obtained simulation results showed similar behaviour compared to those obtained from the experiments considered.
For future work, the formulation allows for an enhancement regarding the modelling of material behaviour such as plasticity, viscosity or anisotropy. Other effects such as reduced damaging under pressure can also be of interest. Moreover, the general framework established in this work can be extended and combined with load optimisation in order to generate not only optimised structures but also optimised loads and loading paths to generate less damaged designs. Since the results of compliance optimisations discussed in this paper (may) differ for different solution strategies, the influence and choice of the particular algorithm should be further investigated in the context of continuum damage mechanics and constitutes future research.

A Derivatives
The calculation of the residual and stiffness terms in Sect. 2 require the total derivatives of the stresses and driving forces given in (19). These derivatives are split into a mechanical and a damage part in the following. The derivatives of the mechanical part read where the special dyadic products are defined via

B Fischer-Burmeister derivatives
In addition, the derivatives of the Fischer-Burmeister residual r d are required for the FEM-implementation, as well as in the derivation of the sensitivity matrix. The specification of these results in These (two) derivatives lead to numerical difficulties when the denominator becomes zero, e.g. in the first inelastic step.
One straightforward approach to circumvent this algorithmic singularity is to add a small constant 1 to the denominator, such that the denominators in the above equations are replaced by [r φ ] 2 + [r λ ] 2 + . Figures 15,16 and 17 depict the result of optimisation problem (101), where the initial design is chosen identical to the result of problem (102), as obtained by the active-set method from the Matlab fmincon toolbox. Fig. 17 Values of the objective functions as shown in Fig. 10 together with the result referred to the design shown in Fig. 15 (Volume New)