Gradient-Enhanced Modelling of Damage for Rate-Dependent Material Behaviour—A Parameter Identification Framework

The simulation of complex engineering components and structures under loads requires the formulation and adequate calibration of appropriate material models. This work introduces an optimisation-based scheme for the calibration of viscoelastic material models that are coupled to gradient-enhanced damage in a finite strain setting. The parameter identification scheme is applied to a self-diagnostic poly(dimethylsiloxane) (PDMS) elastomer, where so-called mechanophore units are incorporated within the polymeric microstructure. The present contribution, however, focuses on the purely mechanical response of the material, combining experiments with homogeneous and inhomogeneous states of deformation. In effect, the results provided lay the groundwork for a future extension of the proposed parameter identification framework, where additional field-data provided by the self-diagnostic capabilities can be incorporated into the optimisation scheme.


Introduction
The adequate numerical prediction of material behaviour under complex loads and boundary conditions is of high significance in the product development context, particularly in view of today's trends towards more efficient and cost-effective production, but also in view of component safety requirements and appropriate component lifetime predictions. To this end, not only the material model formulation itself needs to be physically sound and able to capture all relevant material phenomena, but the corresponding model parameters also need to be identified properly. This process of model calibration is typically carried out with the help of optimisation-based parameter identification schemes.
In this work, we present an optimisation-based framework that enables the identification of parameters for rate-dependent materials under large strains and subject to gradient-enhanced damage by using an efficient two-step approach. In the first step, basic constitutive parameters, in particular elastic parameters, are calibrated based on experiments reflecting homogeneous states of deformation. These parameters are then kept fixed in view of the second step of model calibration, which uses an optimisation scheme based on inhomogeneous states of deformation. The inhomogeneous states of deformation activate gradient terms within the gradient-enhanced damage formulation considered in this work, so that mesh-objective damage effects can be captured within finite element-based simulations. With damage gradient terms activated due to inhomogeneous states of deformation, we proceed to the identification of damage-related model parameters in the second step of the model calibration procedure.
The particular material considered in this work is an elastomer with self-diagnostic capabilities. Self-diagnostic poly(dimethylsiloxane) (PDMS) elastomers are produced by the addition of so-called mechanophore units, which generate a chemical response towards mechanical loads, as discussed by, e.g., Brighenti and Artoni [1]. The specific PDMS elastomer containing a supramolecular detection probe was presented and developed by Früh et al. [2]. In this case, the elastomer reacts with strain-dependent levels of fluorescence when illuminated with UV light, thus enabling an in situ quantification of the material's strain and, potentially, damage state. In this work, however, we restrict our elaborations to the purely mechanical properties of the elastomer, where solely the inhomogeneous displacement field and the load-deflection curve are considered within the model calibration process. A possible future extension of the framework proposed in this work could include an additional consideration of the fluorescence field emitted by the self-diagnostic specimen within the optimisation-based parameter identification procedure.
Different material models allow for the description of the mechanical behaviour of rate-dependent-or rather viscoelastic-materials such as the elastomer considered, mostly using either a convolution integral approach or a differential operator form of the constitutive relation. Reese and Govindjee [3], for example, discussed a non-linear viscoelastic model that utilises a non-linear evolution law that is valid for arbitrary deviations from thermodynamic equilibrium. Bergström and Boyce [4] proposed a constitutive model for the time-dependent behaviour of elastomers, motivated micromechanically by the relaxation of a single entangled chain in a polymer gel. The constitutive viscoelasticity model used in this work is based on a convolution integral approach. The model is based on the well-established finite hyperelasticity framework with uncoupled volumetric-isochoric response that was extended towards viscous effects by Simo and Hughes [5]. The coupling of finite viscoelasticity with gradient-enhanced damage is carried out within the damage regularisation framework provided by Ostwald et al. [6]. This framework is in line with the works of Forest [7] and Miehe [8]. Kiefer et al. [9] presented another approach for a coupling of gradient damage with inelasticity, and Sprave et al. [10] computed complex boundary value problems based on a coupled gradient-enhanced damage formulation. For further, advanced concepts in the context of gradient damage, see, e.g., the work by Liu and Jeffers [11] and the references cited therein.
A corresponding parameter identification scheme for finite viscoelastic material models was introduced by Kleuter et al. [12], where an error square functional-comparing experimentally-obtained inhomogeneous displacement fields with related finite element-based simulations-is minimised by using a gradient scheme. This gradient-based calibration process for finite viscoelasticity was studied in more detail and applied to different materials in [13]. The identification of material parameters for inelastic materials, using experimental data representing non-uniform stress and displacement distributions from the surface of the three-dimensional specimen, was established in the work of Mahnken and Stein [14]. An application of a parameter identification scheme for different constitutive models was presented in [15] where different identification methods, instabilities in least squares problems, and identification for small and non-uniform finite deformations were discussed. Hartmann [16] estimated the constitutive constants of hyperelasticity relations of the generalised polynomial-type and discussed the aspects related to the gradient and convexity behaviour in certain deformations. The particular optimisation methods that are generally applicable were illustrated in, e.g., [17][18][19][20][21].
The present paper is organised as follows: Section 2 introduces the experimental data used in this work, both for homogeneous deformations and for inhomogeneous deformations, as provided in Sections 2.1 and 2.2, respectively. In Section 3, the constitutive framework is introduced, starting with a brief review of the uncoupled finite hyperelasticity framework in Section 3.1. The extension towards viscous effects is provided in Section 3.2 before the coupling to gradient-enhanced continuum damage is established in Section 3.3. The parameter identification framework used, as well as the results obtained are described in detail in Section 4 before the paper is concluded in Section 5.

Experiments
In the following, all accomplished types of experiments for the self-diagnostic PDMS elastomer are briefly presented in order to provide an overview of the experimental data required for the subsequent parameter identification of the implemented material model. For most of the experiments, a micro-testing machine by Kammrath & Weiss GmbH (Dortmund, Germany) was used; cf. Figure 1a. The machine can be combined with a CCD camera system from Chemnitzer Werkstoffmechanik GmbH (Chemnitz, Germany), which is of importance considering the need to track the displacement field for the inhomogeneous deformation tests. Figure 1b presents a screenshot of the software VEDDAC strain. The software was linked with the CCD camera and provided additional information on the displacement field during the experiment. The CCD camera additionally took a photo during each time step and delivered the changes in the displacement field directly to the DDS3 software of the testing machine. In the post-processing, the software VEDDAC could be used to extract all required displacement information from the photos taken during the experiment. Since the material was a transparent soft polymer, a stochastic pattern was sprayed onto the specimens with graphite spray. Consequently, the software was able to track all specific positions of the specimen up to the end of the experiment. The following experiments were conducted with virgin material samples, thereby avoiding the influences of Mullin's effect or the strain history on the material response. The environmental conditions were room conditions, including room temperature.  Experimental setting-the micro-testing machine by Kammrath & Weiss GmbH in combination with the CCD camera system by Chemnitzer Werkstoffmechanik GmbH providing load-displacement data, as well as the displacement field on the surface.

Homogeneous Deformation Tests
As a first step, experiments with homogeneous deformation states were conducted. In order to receive sufficient experimental data for the calibration of the viscoelastic material parameters, three different types of experiments were carried out. Apart from uni-axial tensile tests (cf. Figure  2), relaxation tests (see Figure 3), as well as creep tests (cf. Figure 4) provided the experimental data basis for the following parameter identification; see Section 4. Tensile tests with different strain rates were not considered here, as we instead focused on capturing the viscosity-related parameters of the material via the explicit incorporation of creep and relaxation characteristics within the error square functional. In the case of all these experiments, a homogeneous rectangular specimen with an initial length of 25 mm, depth of 2.2 mm, and width of 10 mm was used. An example of the specimen is shown in Figure 1b at the beginning of the corresponding loading path.

Inhomogeneous Deformation Tests
Apart from the experiments with homogeneous deformation states, an additional experiment with an inhomogeneous deformation state was carried out in order to identify the damage-related material parameters. Considering the regularised gradient-enhanced damage framework, specific information on localisation effects was needed for the identification. Hence, with regard to the possible specimen dimensions, a notched rectangular shape was chosen; see Figure 5a. A constant displacement rate of 20.0 µm/s was applied until the sample with an initial length of 10 mm, width of 10 mm, radius of 2 mm, and thickness of 2.2 mm tore apart. In the figure, prominent states of the experiment are visualised. Moreover, Figure 5e shows the corresponding load-displacement curve.   (e) Load-displacement curve.

Modelling Framework-Finite-Strain Viscoelasticity Coupled with Gradient-Enhanced Continuum Damage
Viscoelastic material models can describe a large class of history-and rate-dependent materials. These materials differ from ideally elastic materials by time-dependent effects such as creep or relaxation. Polymeric rubber materials such as the elastomer considered in this work are typical representatives of this class of materials. Constitutive models for viscoelasticity either use a differential operator or convolution integral representation, where the constitutive model considered in this work is based on the latter approach. The convolution integral approach was described in detail in the works by Simo and Hughes [5], Kaliske and Rothert [22], amongst others. The model formulation considered in this work is introduced as follows: Section 3.1 briefly summarises the classic finite strain hyperelastic constitutive relations with volumetric-isochoric decomposition. Section 3.2 proceeds with the generalisation of these relations to a finite strain viscoelastic model. An extension towards gradient-enhanced continuum damage is then briefly covered in Section 3.3.

Finite Hyperelasticity with Volumetric-Isochoric Decomposition
This section briefly outlines the hyperelastic local constitutive setting considered in the work at hand, which will then be extended towards finite viscoelasticity. The constitutive framework at hand facilitates the accommodation of any arbitrary hyperelastic response based on a local Helmholtz free energy function ψ loc , enabling the use of higher order constitutive relations such as the Yeoh model.
Let F = ∇ X ϕ(X, t) denote the deformation gradient, where x = ϕ(X, t) ∈ B t reflects the current placement of material points with position X ∈ B 0 in the referential placement. The isochoric, i.e., volume-preserving, partF of the deformation gradient F is defined as In this context, J 1 3 I is referred to as the volumetric part of the deformation gradient, with I denoting the second-order identity tensor. The standard and isochoric right Cauchy-Green deformation tensors are introduced via respectively.
In view of the hyperelastic contribution to the constitutive response, a local Helmholtz free energy function ψ hyp of the form is considered, where an additive decomposition into volumetric and isochoric energy contributions, ψ vol and ψ ich , respectively, is taken into account. Within this modelling framework, the hyperelastic part of the Piola-Kirchhoff stress tensor follows as where the isochoric contribution to the hyperelastic part of the Piola-Kirchhoff stress tensor is abbreviated asS for notational simplicity. It is noted that the well-established relations were employed here. The corresponding hyperelastic part of the spatial Kirchhoff stresses, τ hyp = F · S hyp · F t , then follows as : I]I being the deviator operator in spatial representation.

Extension of Finite Hyperelasticity to Finite Viscoelasticity
In view of a generalisation of the volumetrically-isochorically decoupled finite strain hyperelastic modelling framework outlined in Section 3.1 to finite viscoelasticity, C(t) is now assumed to be a function of time t. As a relation for the viscoelastic stress response in terms of Piola-Kirchhoff stresses S at time t, the form is considered and is in line with classic works such as [5,23], where DEV(•) : Here, S hyp (t) denotes the hyperelastic stress contribution defined in (4) and Q i (t); i = 1, 2, . . . , N represent stress-type tensor-valued internal variables. This form of the stress extension leaves the volumetric part of the constitutive response unaffected by viscous effects. Relation (7) reflects a finite strain generalisation of a classic generalised relaxation model with N being the number of Maxwell elements, each consisting of a (linear) spring and a (linear) dashpot. Conceptually speaking, in a simplified and linear setting, each Maxwell element would be characterised by a given spring stiffness, say E i , and corresponding dashpot viscosity, say η i ; cf., e.g., [5].
For the evolution of the stress-type internal variables, rate equations of the typė are considered (cf. [5,23]), where τ i are referred to as relaxation times and where γ i represent relative moduli subject to the restriction with γ ∞ ∈ [0, 1) determining the quasi-static (hyperelastic) limiting case of the material response.
In a simplified and linear setting, as thoroughly introduced in [5], would be the parameters defining the properties of each Maxwell branch in the viscoelastic model. Here, E 0 = E ∞ + ∑ i E i is the stiffness obtained in the (hyperelastic) limiting case reflected by infinitely high strain rates. In the context of the three-dimensional non-linear formulation at hand, this contribution is represented by the general Helmholtz free energy potential ψ hyp that facilitates the incorporation of higher-order and even anisotropic constitutive relations. The linearity of evolution Equation (8) facilitates the standard closed-form convolution representation Combining (10) with (7) and integrating by parts yield an explicit expression for the stress response in terms of the Piola-Kirchhoff stress tensor of the form where the relaxation function is incorporated. The spatial counterpart of (11) is obtained by using , so that the convolution representation of the Kirchhoff stresses reads For more detailed information on the derivation of well-established relations in the context of (undamaged) finite strain viscoelasticity based on a volumetric-isochoric split of the underlying Helmholtz free energy function, the reader is referred to classic monographs such as [5,24], amongst other works.
The numerical solution of the convolution integrals introduced above is based on a transformation towards a recurrence formula, enabling the computation of the constitutive response based on the standard internal variable procedure, thereby bypassing the need to store the entire strain history that is individual to each material point in view of the later finite element application; see Appendix A for details on the numerical update procedure.

Incorporation of Damage
We considered a gradient-regularised damage formulation in view of a mesh-objective solution of boundary value problems by using the finite element method, where the constitutive response can include damage-induced softening. To this end, an overall Helmholtz free energy potential ψ is introduced, composed of local and non-local energy contributions, ψ loc and ψ nloc , in terms of ψ = ψ loc + ψ nloc .
The local Helmholtz free energy contribution, in extension of (3), takes the form where the volumetric and isochoric damage functions, f d ∈ (0, 1] and f n iso d ∈ (0, 1], are evaluated based on a common internal damage variable κ, but allow for an adjustment of the volumetric-isochoric damage distribution via the exponent n iso ∈ R. The exponent n iso is a material parameter that facilitates the modelling of a damage process that affects volumetric and isochoric contributions with different intensities. The non-local energy contribution ψ nloc enables the damage regularisation and consists of two contributions, namely a non-local gradient term ψ grad nloc (∇ X φ; F) and a non-local penalty term ψ Here, φ is the non-local damage variable that is coupled to the local damage variable κ via the penalty parameter β d . Moreover, c d represents a length scale-type regularisation parameter. With this particular choice of the non-local energy contributions, the Euler-Lagrange equations for the non-local damage field variable in spatial format take the form see, e.g., [6] for more details on the associated relations, as well as on the numerical treatment within appropriately coupled finite element formulations. The incorporation of the damage functions in the local free energy potential (14) induces Piola-Kirchhoff stresses of the form corresponding to spatial Kirchhoff stresses given by both of which are direct extensions of their purely hyperelastic counterparts introduced in (4) and (6). The damage function affecting the isochoric part of the free energy function furthermore affects the evolution of viscous stress-type internal variables viȧ which basically reflects the damage-extended version of (8). This induces a convolution representation in terms of see (10) for the corresponding expression that is associated with an undamaged setting. Finally, the Piola-Kirchhoff stresses that fully incorporate viscous and damage effects take the form which allows for the accommodation of arbitrary free energy functions, in particular energy functions of higher order such as the Yeoh model that appropriately captures the Treloar-type behaviour of materials such as the particular elastomer considered in this work. The numerical solution of the above-introduced convolution integrals is based on a transformation towards a recurrence formula, enabling the computation of the constitutive response based on the standard internal variable procedure, thereby bypassing the need to store the entire strain history that is individual to each material point in view of the later finite element application; see Appendix A for details on the numerical update procedure.

Parameter Identification
This section discusses the optimisation-based procedures employed for the identification of the constitutive parameters of the finite strain viscoelastic material model coupled with gradient-enhanced continuum damage as introduced in Section 3. Since the general framework can be applied to all kinds of energy functions, at first, considering the type of material and its behaviour, the chosen hyperelastic material model is briefly presented in Section 4.1. Next, the applied damage function f d is specified in Section 4.2. The main idea is to first carry out a parameter identification based on homogeneous states of deformation as described in Section 4.3. The obtained set of material parameters is then used as the initial value for the parameter identification based on inhomogeneous states of deformation, as elaborated in Section 4.4. The gradient terms of the gradient-enhanced damage framework are activated only in the case of inhomogeneous deformations.

Hyperelastic Material Model
Considering the nearly incompressible, non-linear elastic material behaviour (see for example the stress-strain-response in Figure 2b), the Yeoh hyperelastic material model was chosen [25,26]. In [27], Yeoh and Flemming argued that Rivlin's theory [28], where the strain energy function depends on the first two invariants of the right Cauchy-Green tensor C, to be specific runs into difficulties with respect to the identification of the material parameters since the contributions cannot be perfectly determined separately in the experiments. Here, C ij are the corresponding material parameters, and I 1 and I 2 are the first and second principal invariants related to the isochoric right Cauchy-Green tensor, respectively. Yeoh neglected the contribution of the second invariant assuming that no serious error would arise since the contribution would be sufficiently small with respect to the impact of the first invariant. In this work, however, the Yeoh energy potential serves as the isochoric contribution to overall energy potential, where C i denote the underlying material parameters. This phenomenological material model was developed for rubber elasticity and is applied in this work as a first choice with respect to the mentioned material behaviour. The volumetric energy contribution ψ vol is chosen in a standard manner as where K denotes the bulk modulus of the material. The related energy contributions ψ ich and ψ vol can, however, be straightforwardly replaced within the general constitutive framework elaborated in this work in order to account for different particular material characteristics.

Damage Function
The damage function mentioned in Section 3 was chosen to be of an exponential type, following [6,29]. Thus, a damage initiation threshold, as well as a damage saturation rate can be used in the function fulfilling the requirements for the damage function to be restricted to f d : R + → (0, 1]. Here, d denotes the classic damage variable with d = 0 indicating zero damage and d = 1 for 100% damage. The variable κ represents a related internal damage variable, while κ d > 0 is the damage threshold parameter, and η d introduces the damage saturation parameter. In addition, • := max{0, •} represents the Macaulay brackets. Consequently, damage is obtained in the material if the local damage variable exceeds the damage threshold. The damage function for the isochoric energy contribution matches this chosen function except for the exponent n iso . Thus, a relation between both functions is still given and only adjusted via this material parameter. The evolution of the internal damage variable κ is based on the associated formκ where λ denotes a proper Lagrange multiplier, Φ d (F, φ, ∇ X φ, κ) = q(F, φ, ∇ X φ, κ) − κ ≤ 0 represents the damage condition, and q = − ∂[ψ loc + ψ nloc ]/∂d is the energy release rate. Hence, apart from the balance of linear momentum, the Euler-Lagrange equation governing the non-local damage variable needs to be solved simultaneously. For further information regarding this regularised damage model and its implementation, see Ostwald et al. [6]. The implementation was performed in a user material subroutine (UMAT) in Abaqus following [6].

Parameter Identification Based on Homogeneous States of Deformation
In the context of parameter identification procedures, it is advantageous to compute the Piola stress tensor as a quantity, facilitating the comparison of experimental and simulated material response; cf. Section 2. The Piola stress tensor P is related to the Piola-Kirchhoff stress tensor S and the Kirchhoff stress tensor τ via where the sensitivity of the Piola stress with respect to the deformation gradient is denoted as A, in particular, or, expressed in index notation, Moreover, A relates to the spatial elasticity tensor e via It is noted that the tangent operator, as typically required within implicit finite element formulations, respectively so-called constitutive drivers, is based on A = dP/dF. The Yeoh-type hyperelastic energy function was combined with the constitutive viscoelasticity model, presented in Section 3, including two Maxwell elements. Apart from the Yeoh material parameters C 1 , C 2 and C 3 , the viscoelastic parameters, γ 1 and γ 2 , denoting the relation of the stiffness of each Maxwell element with respect to the pure elastic Young's modulus, and τ 1 , τ 2 , representing the ratio of the viscosity to the stiffness in each Maxwell element, need to be identified. Regarding the parameters, it has to be mentioned that the material was assumed to be nearly incompressible, and thus, Poisson's ratio was fixed to ν = 0.49. Thus, seven material parameters needed to be identified.
Since these seven material parameters could be identified via experiments of homogeneous deformation states, strain-and stress-driven constitutive drivers were implemented in MATLAB instead of full finite element (FE) simulations in order to reduce the computational cost within each iteration of the parameter identification. The pseudo-codes of both types of constitutive drivers are depicted in Algorithms 1 and 2. The fminsearch-algorithm in MATLAB was used for the parameter identification, minimising the goal function with respect to the difference in the simulated and experimental reaction force over all load steps, i.e., where w T , w R , and w C denote the weighting factors for the tensile, creep, and relaxation test, respectively. The additional weighting factors w T t , w R t , and w C t were introduced to emphasise specific time steps t of each experiment.
In the following, the viscoelastic material parameters are identified via the three tests based on homogeneous deformation states, presented in Section 2.1.
If the parameters of the Yeoh-model were solely fitted with respect to the tensile test, the simulated material response of the identified parameters matched the experimental data perfectly. Next, the relaxation parameters were obtained with respect to the relaxation and creep test, though the material behaviour highly depended on the previously identified Yeoh parameters. Thus, the experimental relaxation and creep response could not be sufficiently matched via optimisation of the relaxation parameters only. Consequently, all of the material parameters were identified simultaneously for all three different experiments in order to obtain the best parameter set reflecting the complete material behaviour. The comparison of the simulated and experimental material behaviour is shown in Figure 6 for the optimised parameter set. The corresponding parameter set is presented in Table 1.
In the graphs shown in Figure 6, the experimental results of the relaxation and creep tests are nearly perfectly matched, and only slight differences in the diagrams can be seen at certain points. The simulated material response of the tensile test overestimated the experimental curve, but captured the general trend. The tensile test could be captured better by increasing the associated weight within the objective function (33) at the cost of the accuracy with which creep and relaxation tests were captured. The focus here was set on the rate-dependent constitutive characteristics that were captured with high accuracy. Though beyond the scope of this work, a precise capturing of tensile, creep, and relaxation test at the same time could be achieved by introducing higher order energy functions and a larger number of independent Maxwell branches within the constitutive framework provided. The accompanying increase in parameter identification complexity due to a significantly increased number of constitutive parameters could then be dealt with using parameter correlation matrices.

16:
(k) F n+1 = (k) F 11 e 11 ⊗ e 11 + (k) F n+1 . 17: partition the tangent operator: 18: . 19: update the transverse deformation gradient: 20: 21: end while 22: assemble the deformation gradient and the stress tensor F n+1 , P n+1 Considering the curve of the tensile test, especially the values and the signs of the Yeoh parameters are important, since the first parameter C 1 weights the linear behaviour of the stress strain-response, whereas C 2 scales the quadratic term of the isochoric energy contribution resulting in the decreasing slope. The third parameter C 3 corresponds to the third-order polynomial contribution of the energy function and thus yields the increasing slope at the end of the loading path. Thus, different values for the Yeoh parameters would probably fit the tensile response way better; the dependence on the relaxation parameters, however, would result in a worse fit of the relaxation and creep behaviour. Higher order contributions to the Helmholtz free energy, in particular the isochoric contribution, would therefore provide a potential in order to capture the tensile test response better. n , ∀ i = 1, . . . , N. 2: get material parameters C 1 , C 2 , C 3 , γ 1 , τ 1 , . . . , γ i , τ i . 3: initialize F n+1 . 4: for every time step t n do 5: given: Piola-Kirchhoff stress tensor P exp . 6: while P − P exp < tol do 7: compute (k) F n+1 , (k) J n+1 , (k) C n+1 , (k) F n+1 , (k) C n+1 . 8: compute initial elastic Kirchhoff stress tensor τ hyp n+1 . 9: compute algorithmic internal variablesS compute (k) S n+1 according to (A6). 11: compute the spatial tangent moduli (k) e n+1 . 12: compute Piola-Kirchhoff stress tensor (k) P n+1 . 13: compute tangent operator (k) A n+1 . 14: compute (k) ∆P n+1 = (k) P n+1 − (k) P exp n+1 . 15: update the deformation gradient 16:  The identified viscoelastic material parameters were used as fixed values for the subsequent optimisation of the damage-related parameters.

Parameter Identification Based on Inhomogeneous States of Deformation
In the following, the previously identified Yeoh and relaxation parameters are used to identify the damage-related material parameters η d , κ d , n iso , and c d via the tensile test with inhomogeneous deformation states presented in Section 2.2. In contrast to the identification of the viscoelastic parameters, the constitutive driver is not sufficient enough to capture the material behaviour of the gradient-enhanced damage model. Thus, a finite element (FE) formulation is required. Therefore, the material model was implemented into a UMAT in Abaqus, as already mentioned in Section 4.2.
In addition, a parameter identification tool was implemented in Python, using the Nelder-Mead-Simplex algorithm of the scipy-package. In contrast to the calibration of the Yeoh and relaxation parameters, it did not suffice to solely consider the difference in the Piola stress P 11 or the (clamping) displacement u 1 , depending on the experiment, for the goal function f obj of the optimisation. Since the gradient terms-capturing the mesh-objectivity of damage-of the regularised damage framework are activated by inhomogeneous deformation states, the difference in the displacement field between the simulated and experimental material response needs to be added to the goal function apart from the difference in the reaction force where n np denotes the number of node points considered, w u the weighting factor for the displacement contribution, and w F the weighting factor for the impact of the reaction force. The number of node points was the element nodes on the surface of the specimen in the FE simulations. The experimental data was interpolated onto those node points via a 2D-interpolation scheme prior to the parameter identification following Kleuter [13]. For the purpose of reducing the computational cost of the FE simulation within each iteration of the parameter identification, the symmetry properties of the sample were used; cf. Figure 7a. Furthermore, considering the material properties of the soft polymer, only the segment shown in the figure was used to exclude boundary effects of the clamping jaws. In order to still use the experimental boundary conditions in the simulations, the experimentally measured displacements, taken via the CCD camera system, were applied to the right boundary of the specified segment. As can be seen from the photos of the specimen during the experiment (cf. Figure 5) and from the sketch of the segment, the displacements were not uniform over the boundary of the chosen segment of the sample.
The initial guess for the damage-related model parameters was η d = 0.002, κ d = 0.4, n iso = 1.0, and c d = 0.04 and resulted in no damage initiation in the material; cf. Figure 8a. In Figure 7b, the comparison of the load-displacement curves is presented. Since the stress-strain path of the homogeneous tensile test already overestimated the experimental curve, the simulated response for the tensile test with inhomogeneous deformation states lied above the experimental result as well. Furthermore, considering the large total stretch values present in combination with the identified Yeoh material parameters, the increasing response for the reaction force was triggered by the parameters C 2 and C 3 weighting the hyperelastic energy contributions of second and third order.
Since parameter n iso strongly influences the simulated deformation behaviour of the specimen in each iteration, n iso was fixed to 1.0 in a first step, thereby neglecting a different damage contribution of the volumetric and isochoric part. The corresponding optimised parameter set, denoted as Final #1, η d = 0.202531, κ d = 0.15185, and c d = 0.52593, was rounded to five decimal places. Considering the load-displacement curves, the response of the final set lied closer to the experimental behaviour than the initial guess; see Figure 7b. The stress distribution of the initial guess and the final set was comparable; the magnitude of the final set, however, was less than half of the stress of the initial set; cf. Figure 9. The von Mises stress distribution within the specimen-with the maximum stress value obtained at the upper surface, corresponding to the region dominated by the initially circular notch-was in line with observations made by, e.g., Kleuter [13].
The contour plot of the damage value d in Figure 8 compares the damage value obtained for the initial parameter set with the damage value obtained for the Final #1 set. Initially, no damage was initiated, and with the Final #1 set, a damage value of nearly 75% occurred in the symmetry plane of the loading direction.
In Figures 7b and 10, the influence of the different damage functions for the volumetric and isochoric contributions is visible. The set Final #2 included n iso = 0.9, and Final #3 considered n iso = 1.1. A higher value for n iso flattens the load-displacement curve, while a lower value provides an increased path. Another important impact of the parameter is indicated in Figure 10. Apart from the difference in the stress magnitude, the deformation of the sample was already different for slight changes in the parameter with regard to the necking of the sample at the symmetry plane of the loading direction.
The advantage of the regularised damage framework is the mesh objectivity of the simulated material response. To demonstrate this feature, the boundary value problem was calculated with three different meshes for the optimised parameter set. In addition to the mesh containing 1548 elements, used for the results in Figures 7-10, a coarser mesh with 1092 and a finer mesh with 2613 elements were used in order to analyse the mesh sensitivity of the results. As can be seen in Figure 11, only marginal differences in the contour plots of the damage function values and the von Mises stress, as well as the load-displacement curves are visible. Thus, the gradient-enhanced damage model is working properly. In contrast, it is noted that the local damage model diverged at different load steps, depending on the mesh discretisation. In the case of the fine mesh, the solution diverged at a displacement of 3.2 mm, for the basis mesh at 3.375 mm and for the coarse mesh at 3.415 mm. The displacement was taken from the bottom node of the right boundary in the symmetry plane.

Conclusions
In this work, a parameter identification framework for gradient-enhanced damage in rate-dependent materials under finite strains was introduced. In order to identify the set of over 12 material parameters, an efficient scheme for the process was carried out. At first, the basic constitutive model parameters, namely elastic parameters, were calibrated with respect to experiments displaying homogeneous states of deformation. Considering the deformation states, the computational cost of the identification was further reduced by using constitutive drivers for the simulations of all three required types of experiments within each iteration of the calibration process. In the next step, the already identified parameters were fixed during the identification of the damage-related material parameters.
For the purpose of developing a general framework for the calibration of rate-dependent materials coupled to damage, the constitutive viscoelasticity model by Simo and Hughes [5] was coupled to the gradient-enhanced damage formulation. Thus, the advantages of both models were combined, enabling different types of hyperelastic energy formulations-considering the energy-independent viscoelastic derivations by Simo and Hughes [5]-as well as mesh-independent results with respect to the damage evolution due to the damage regularisation framework by Ostwald et al. [6].    To investigate whether the framework works even for very large deformations, a self-diagnostic elastomer incorporating the mechanophore units within the polymeric microstructure was chosen as an application. With regard to the material behaviour of the soft polymer, the hyperelastic model by Yeoh and Fleming [27] was applied to the viscoelastic damage model.
The results of the parameter identification in Section 4 showed that the general framework of the implemented material model, as well as the parameter identification tool works properly. Even though the chosen self-diagnostic elastomer was only a first attempt for the application of the framework, the experimental results were already basically matched; cf. Figures 6 and 7b. Since the simulations for the experiments with inhomogeneous deformation states were conducted before the sample started to tear in the necking region, the damage distribution with the optimised parameter set-showing damage evolution of nearly 75% in the symmetry plane of the loading direction, i.e., the necking region-reflected such a material behaviour.
Nevertheless, further improvements could be achieved for example by using different weighting factors for the contributions of the three different experiments with homogeneous deformation states to the goal function of the optimisation in order to improve the calibration of the parameters with respect to the tensile test. However, after all, the largest improvement could be achieved by a refinement of the experiments. The material properties led to difficulties in the sample manufacturing-such as uniformly cutting the shape-and in the clamping of the specimen, since the clamping jaws already generated deformations in the soft polymer. An improvement would be if the samples were directly moulded or if clamping were conducted with a very accurate torque spanner. Furthermore, experimental data based on different strain rates could be included in the identification framework.
In addition, in the next step, the load path for the inhomogeneous experiments could be changed to cyclic loading, improving the damage characterisation, though increasing the computational cost within each iteration of the parameter identification procedure. Nevertheless, further experiments including unloading should be included in the calibration process to enhance the accuracy. Until now, only damage activated prior rupture has been taken into account. As a long-term goal, the current framework should be combined with a rupture model to improve the overall prediction of the material behaviour.
After all, the general parameter identification framework provides an efficient scheme for the calibration of gradient-enhanced damage in rate-dependent material models under finite strains. In the future, an interesting extension of the parameter identification framework is the incorporation of the self-diagnostic properties of the elastomer in terms of additional field data that is considered within the objective function.  where H (i) (t), i = 1, 2, . . . , N, i.e., we have one stress-type tensor-valued internal variable H (i) for each Maxwell element considered in the generalised relaxation model; cf. Simo and Hughes [5].
Using the semigroup property exp([t + ∆t]/a) = exp(∆t/a) exp(t/a) and the property of the additivity of the integral over the interval of integration, we obtain the recurrence relation n iso d (κ)DEV (2 ∂Cψ ich (C(s))) ds . (A2) Finally, the integral (A2) is approximated using the midpoint rule, so that using the abbreviations S n+1 hyp := DEV n+1 (2 ∂Cψ ich (C n+1 )) andS n hyp := DEV n (2 ∂Cψ ich (C n )) , where DEV n and DEV n+1 are computed with respect to C n and C n+1 , respectively, consistent with The algorithmic approximation of the Piola-Kirchhoff stress tensor at time t n+1 then takes the form hyp representing relation (A4) evaluated at the current time t n+1 and with κ n+1 the (updated) internal variable representing the evolution of damage within the material. For details on the straightforward time-integration of κ via a standard Euler backward scheme, the reader may refer to, e.g., [6,29]. The spatial counterpart of (A6), i.e., the current approximation of the Kirchhoff stress tensor, is then computed using a standard push-forward operation in terms of