Comparison of discontinuous damage models of Mullins-type

The Mullins effect is a characteristic property of filled rubber materials whose accurate and efficient modelling is still a challenging task. Innumerable constitutive models for elastomers are described in the literature. Therefore, this contribution gives a review on some widely used approaches, presents a classification, proves their thermodynamic consistency, and discusses reasonable modifications. To reduce the wide range of models, the choice is restricted to those which reproduce the idealised, discontinuous Mullins effect. Apart from the theoretical considerations, two compounds were produced and tested under cyclic uniaxial and equibiaxial tension as well as pure shear. Based on this experimental data, a benchmark that compares the fitting quality of the discussed models is compiled and favourable approaches are identified. The results are a sound basis for establishing novel or improving existing rubber models.

Instead, these models should be considered as a starting point of more complex approaches. For instance, [4] proposed a model based on such an idealised softening and added a rate-independent, plastic contribution which generates equilibrium hysteresis and permanent set. Also the derivation of the elastoplastic model by [27] which was extended to viscoelastoplasticity by [36] starts from an idealised softening. Further, [19] presented a viscoelastic model whose equilibrium stress is affected by idealised softening. In spite of the limitations, Mullins-type damage models on their own are a reliable alternative for approximating the behaviour of filled rubber in certain applications. Since the choice of models is always a trade-off between numerical effort and the realism of their behaviour, the simulation purpose must be kept in mind. Here, a good prediction of the structural stiffness depending on the maximum load the material has experienced so far can be obtained. For instance, the structural stiffness is of crucial importance for damping components that either experience softening in use or are prestretched within the production process. Another example are seals whose closing pressure may be affected by the material softening. Moreover, the static behaviour after cyclic loading, e.g., interrupted fatigue tests can be estimated by damage models.
The present paper conducts a benchmark of the aforementioned damage models. Lists or rankings of rubber models have indeed been compiled in the recent years for example by [10], [28] or [39]. However, these publications deal with purely elastic models only. Looking for more sophisticated models, comparative reviews typically deal with two or three particular models only. For example, [17] compared only two damage approaches each combined with the polynomial strain energy function by [43]. In contrast, this paper pursues a more general comparison and takes fifteen approaches into account. In addition, four modifications are presented and included in the comparison.
Aiming for the best fit quality, each approach is tested with several hyperelastic models, see Sect. 3.1. Therefore, discontinuous damage models that can be added to an arbitrary, isotropic hyperelastic strain energy function Ψ 0 suitable for large deformations are taken into account. Moreover, models with increased numerical effort, e.g., numerical integration, sums over material/chain directions, numerical methods for differential equations or split of the deformation gradient which separates the damage kinematics are excluded. Finally, we restrict ourselves to isochoric damage since experimental findings indicate that the volumetric part can be treated as nearly perfectly elastic, see for example [16] or [37]. Summing up, this contribution aims for a comparison and discussion of numerical efficient, generalisable, isotropic, discontinuous modelling approaches for isochoric Mullins-type damage.
The paper is organised as follows: Sect. 2 provides the required basics and thermodynamics of damage models as well as an introduction to the fitting procedure. The theories and constitutive equations of the damage models and the basic hyperelastic models are briefly summed up in Sect. 3. After describing the experimental procedure in Sect. 4, the results are discussed in detail in Sect. 5.

Basic kinematics
Considering a material point with the position vector X in the reference configuration and position vector x in the current configuration, the according deformation gradient and the right Cauchy-Green tensor are given by respectively, where F T is the transposed deformation gradient. The square-root of the eigenvalues of C are called principal stretches and are denoted by λ x , λ y and λ z . However, the deformation can be multiplicatively decomposed into a volumetric and an isochoric part respectively, where I and det( · ) denote the unit tensor and the determinant. The isochoric right Cauchy-Green tensor and its principal invariants are defined as with the inverse isochoric right Cauchy-Green tensorC −1 and the trace tr ( · ). Accordingly, the square-root of the eigenvalues ofC are called isochoric stretchesλ x ,λ y andλ z .

Thermodynamics
The local, isothermal Clausius-Planck inequality per unit reference volume in terms of the Helmholtz free energy density 1 is given by with S :Ċ = S abĊab (sum over a and b), see for instance [18]. Basically, this inequality demands that for thermodynamic consistency the dissipation rate D m must be non-negative at any deformation state. D m is equal to the mechanical stress power 1 2 S :Ċ less the change of the free energyΨ . Herein, S denotes the 2nd Piola-Kirchhoff stress andĊ the material time derivative of C.
Assuming the free energy to be a function of C and a scalar history variable Γ , the derivative of the free energy with respect to time readsΨ = ∂Ψ ∂ C :Ċ + ∂Ψ ∂ΓΓ leading to Following the standard argumentation of [9], this inequality can be fulfilled for arbitrary deformations by defining and ensuring Equation (6) defines the 2nd Piola-Kirchhoff stress, whereas Eq. (7) is a restriction on the damage model. 1 For the sake of readability, the Helmholtz free energy density is briefly referred to as free energy.
The idealised Mullins effect is commonly modelled by introducing a scalar measure of the maximum load Γ as a history variable. This measure is typically chosen to be an invariant of a deformation tensor or the basic strain energy itself. Moreover, all considered models introduce an internal state variable whose evolution is a function of this measure. Depending on the type of the internal state variable, the models may be classified into three groups: models with 1. a virgin state variable η, 2. a damage variable d, 3. an amplification variable X .
However, as we assume that the damage affects the isochoric behaviour only, the free energy is additively split into an isochoric and a volumetric part Ψ = Ψ iso + Ψ vol such that Eqs. (6) and (7) can be rewritten as and Here, m(Γ ) can be replaced formally by η(Γ ), d(Γ ) or X (Γ ).

Parameter fitting
To find the best sets of parameters for the tested models in terms of the least mean squared error, a standard Trust-Region algorithm is used for the model calibration. The error is defined as the normalized deviation between the model and the experimental 1st Piola-Kirchhoff stress P = F · S in loading direction x, viz.
between P x x,mod = λ x S x x,mod and P x x,exp . That leads to a minimisation problem with the following cost function m ux , m ps , m bx are the number of experimental observations, viz. the load increments, considered for fitting. i − 1 and i refer to the beginning and the end, respectively, of the i-th increment. The abbreviations ux, ps and bx stand for the tested loading modes: uniaxial tension, pure shear and equibiaxial tension, cf. Sect. 4. The normalisation is carried out with respect to the maximum stress in the fitting range of each experiment. The division of the squared residuals r 2 i by the number of experimental data points ensures an equal weight of all loading modes. n is the number of parameters. The model stress S x x,mod,i is a function of the parameters p j ( j = 1, . . . , n), the history variable at the beginning of the current load step Γ i−1 (which must be updated to Γ i ) and the measured stretches. The stretches are turned into a corresponding deformation gradient F exp,i under the assumption of perfect incompressibility. This assumption is a good and pragmatic approximation for filled, technical elastomers up to moderate strains, see for example [26]. Due to the incompressibility, the volumetric stress contribution in Eq. (8) is not anymore computed from the volumetric free energy Ψ vol as but is given by where the hydrostatic pressurep stems from the stress-free condition in lateral direction.
To speed up the fitting procedure, an exact Jacobian is used within the least squares method. For material models with a history variable, it is obtained as Hence, at each time increment i, the history variable Γ i as well as its derivatives with respect to the material parameters must be updated and stored.
As the Trust-Region algorithm is only locally convergent, the whole fitting procedure is repeated ten times with different, randomly generated initial guesses according to the Latin hypercube sampling. More precisely, for each parameter of a model, an initial guess range is prescribed which is equally divided into ten intervals. In each interval a random number is drawn from a uniform distribution. These values are then randomly combined to ten different sets of parameters which serve as the initial guesses for the fitting procedure.

Root mean squared error and correlation matrix
Since the normalized residual of the minimisation problem (10) is not very descriptive, the root mean squared error (RMSE) with respect to the absolute error is used for the discussion of the results. As explained in Sect. 5, not all experimental data points are considered for the fitting procedure. Instead, the data outside of the fitting range (i.e., inside the predicting range) is used to evaluate the predictivity of the models. For this purpose, Eq. (15) is applied analogously to the data in the predicting range yielding the root mean square prediction error (RMSPE). One desirable property of material models is a low correlation between the parameters. The Pearson correlation matrix ρ i j is obtained from the covariance matrix D i j as In case of least square methods, the covariance matrix can be computed from the approximated Hessian H i j , i.e., the 2nd derivative of the cost function F p j and the scaled mean square error see for example [14]. The correlation matrix component ρ i j ∈ [0, 1] describes the correlation between the i-th and j-th parameter. The matrix ρ i j is symmetric and ρ ii = 1 ∀i. Values close to 0 represent a low correlation. Whereas values close to 1 indicate a high correlation stemming from either an overparameterised model or from an improper experiment which provides insufficient information about the material behaviour.

Damage models
3.1 Choice of basic hyperelastic models As described in Sect. 1, the considered damage models must be combined with a basic hyperelastic model defined by its strain energy function Ψ 0 . To obtain a purely isochoric damage, Ψ 0 must be defined in terms of the isochoric right Cauchy-Green tensorC. Due to the endless list of proposed strain energy functions, we restrict ourselves to three promising models: (a) The polynomial approach with terms up to 6th order of deformation, i.e., up toλ 6 whereĪ 1 ∝λ 2 andĪ 2 ∝λ 4 [20] reads with five parameters c 10 , c 20 , c 30 , c 01 , c 11 ∈ [0, ∞). The initial shear modulus is given by G 0 = 2 (c 10 + c 01 ). All material parameters units are shown in Table 6. (b) Alternatively, the invariants themselves can be exponentiated, see [40]. This approach leads to [40], the exponents are fixed to certain values rather than fitted because the Trust Region algorithm mostly tends to push them to very high or low values. This may lead to numerical problems and high parameter correlations. Here, the exponents are fixed according to the approach of [7]: (c) A simplified version of the extended non-affine tube model from [21] was proposed by [35]: The criteria for the choice of these three models are as follows: The strain energy function must be formulated in a closed-form and in terms of elementary functions since some damage models use the basic strain energy as measure of maximum load. For this reason, the models by [1] as well as [25], among others, are not included here. The stress computation must not depend on numerical integration methods to minimize the computational effort. Therefore, the micro-sphere model by [30] is not considered here.
The following Sects. 3.2 to 3.4 explain in detail the three classes of damage models given in Sect. 2.2. In each of the sections, suitable model formulations from the literature are introduced which will be used for the comparison thereafter. To point out the different motivations of the model classes, Fig. 2 depicts qualitatively their behaviour.

Damage models with virgin state variable
The first group of damage models assumes that the basic strain energy function Ψ 0 describes the virgin state curve. For this purpose, a virgin state variable η ∈ [0, 1] is introduced, such that where η = 1 indicates a primary loading and η = 1 an un-or reloading. [33] coined the name pseudoelasticity for this model type. The idea was first described by [11] and generalised by [32]. Furthermore, the latter authors proved that the only feasible choice for the scalar measure of maximum load is the basic strain , the maximum value of the basic strain energy within the loading history. The total free energy reads such that and To define the function η Ψ 0 , Ψ 0,max , two general approaches are presented in the literature. On the one Here, the basic strain energy function must be restricted to Ψ 0 ≥ 0 to ensure Ψ 0 Ψ 0,max ∈ [0, 1]. Demanding η to be a monotonically increasing function (24). This approach includes the model by [11]: Unfortunately, the authors did not give an explicit formula for η but plots instead, see their Figs. 1 and 7. Thus, the function given above is a reconstructed approximation reproducing these plots. To ensure η(1) = 1 and to consider only the principal branch of the tangent, following three fitting parameters are defined: c ∈ 0, π 2 , To prove the capability of the presented approximation, Fig. 3a depicts the reconstructed plots from [11].
< 0 is sufficient to fulfill inequality (24), see [32]. A discrepancy between these demands for monotonicity and experimental findings are discussed by [22]. Based on uniaxial tension tests with η = (Sxx−Svol,xx) S 0,xx , they showed that η is in general non-monotonic. The second approach is applied for instance by The material parameters are r ∈ [0, 1] and m, q ∈ [0, ∞) and the units of all parameters are given in Table 6. The evolution equations are compared in Fig. 3b. Models that do not fulfil for instance [41]. 2 The parameter q = 1 2 is fixed to reduce the number of fitting parameters. Note that modified evolution equations will be indicated by an asterisk. 3 The parameter r is introduced in the present paper to improve the model capability. the amplification variables X and X max according to models 3.1a/3.1b and 3.2/3.2* with X 0 = 10, X ∞ = 1, and γ such that X (5) = 3 and X max (5) = 3, respectively

Damage models with damage variable
The second class of damage models, as given in Sect. 2.2, introduces a damage variable which depends on the scalar measure of maximum load d = d(Γ ). The desired stress response is Here, d = 0 denotes the virgin state and d = 1 a totally damaged material. The basic stress response S 0 acts as an upper limit (if d → 0). The corresponding damage evolution maps d : [0, ∞) → [0, 1] with d(0) = 0 and the free energy is given by with the dissipation rate SinceΓ ≥ 0, the restriction Ψ 0 ∂d ∂Γ ≥ 0 has to hold true for any deformation to fulfil inequality (28). This requirement cannot be satisfied without further restrictions on the basic strain energy because any constant term can be added to Ψ 0 without affecting the stress-deformation relation Eq. (9) but the extent and sign of dissipation. Thus, equal conditions are enforced by demanding Ψ 0 (I) to be a global minimum with Ψ 0 (I) = 0. The restriction reduces then to ∂d ∂Γ ≥ 0, i.e., d has to be a monotonically increasing function. In contrast to the models using a virgin state variable, the measure of maximum load is not restricted to the basic strain energy and, hence, can be chosen more arbitrarily. The following list compiles the considered models along with the measure of maximum load and the damage evolution (cf. Fig 4a): The material parameters are α ∈ [0, ∞) and β ∈ [0, 1]. See Table 6 for the parameter units.

Damage models with amplification variable
The previous two model classes are phenomenologically motivated. In contrast, the third class (cf. Sect. 2.2) is based on the physical concept of strain amplification. Due to the presence of rigid filler particles, an inhomogeneous deformation field and a locally amplified strain of the polymer chains is assumed. The introduction of an amplification variable X seems likely such that with the macroscopic stretchλ 0 , see [31]. Moreover, a deformation of the rubber will lead to a breakdown of filler clusters and to a debonding of polymer chains from the filler particles, see for instance [12] for physical interpretations of the Mullins effect. The ruptured filler particles have a less amplifying effect on the polymer network, i.e., X decreases. Note that the considered approaches do not model the filler stress response itself, that is, the filler network has no load-bearing capacity (for a counterexample see [35] who added a filler contribution to the basic strain energy). Hence, the basic hyperelastic response S 0 = 2 ∂Ψ 0 ∂ C can be interpreted as the response of an unfilled rubber which is equivalent to a totally damaged state containing broken filler clusters only. Thus, S 0 is a lower limit for X → 1.
From a numerical point of view, the strain amplification described above is unattractive since an eigenvalue computation is required for every model. Therefore, a more practicable amplification of the principle invariants is introduced. In contrast to [3] who presented an invariant amplification restricted to the Neo-Hooke model, a more general approach suitable for many invariant-based strain energy functions is proposed here. Noting thatĪ 1 ∝λ 2 andĪ 2 ∝λ 4 , a feasible, consistent amplification is given by following replacements x denotes the Tresca invariant of the isochoric right Cauchy-Green tensor. 5 · denotes the Frobenius norm.
The former replacement Eq. (30) is applicable to the basic strain energy function (a) from Sect. 3.1, whereas the latter replacement Eq. (31) is suitable for strain energy function (b). In case of (c), Eq. (30) is used for theĪ 1 -part and Eq. (31) for theĪ 2 -part. Amplified strain energy functions that are obtained by applying the replacements Eq. (30) or (31) to a basic strain energy functions Ψ 0 will be denoted by Ψ * 0 X,Ī 1 ,Ī 2 . Since the limitations pointed out in Sect. 1 lead to narrow confines of the paper's scope, a lot of highly specialised approaches are not considered here. Thus, besides the amplification concept, no further physically or micro-mechanically motivated models are discussed. For an overview of such sophisticated models, the reader is referred to [12].
The principle of strain amplification was considered for instance by [23]. Originally, they applied the strain amplification according to Eq. (29) to a generalized tube model. However, the more convenient invariant amplification (cf. Eqs. (30) and (31)) is used here such that ( 3 2 ) [23] proposed two evolution equations reading can be satisfied by ∂Ψ iso ∂ X ≥ 0 and ∂ X ∂Γ ≤ 0. The former inequality is fulfilled if Ψ 0 is a monotonic function w.r.t. the principal invariants and an amplification according to Eq. (30) or Eq. (31) is assumed. Whereas the latter inequality demands a monotonically decreasing evolution equation X (Γ ).
Another generalisable model is outlined by [34]. They assumed a non-uniform distribution of strain amplifying filler domains on the micro-mechanical scale and presented a mathematical homogenisation to obtain a representative, macroscopic free energy function. Transferring their approach to arbitrary strain energy functions in terms of principal invariants, one obtains where P denotes the distribution of the local amplification X which can be represented as The integration limits 1 and X max represent the range from non-amplified to the maximum amplified domains.
[34] suggested a power law distribution Their evolution equation and the measure of maximum load read Note that the parameters X max,0 = X max (0) = 1000 and X max,∞ = X max (Γ → ∞) = 1 are fixed to minimize the number of fitting parameters and their correlation. The remaining parameters are χ ∈ [1, ∞) and γ ∈ [0, ∞).
The mechanical dissipation rate is obtained by Since − ∂ X max ∂Γ ≥ 0 andΓ ≥ 0 by construction, ∂Ψ iso ∂ X max ≥ 0 has to be proven: Using Eq. (35) 3 , this is equivalent to The first mean value theorem states, ∃k ∈ [1, X max ] such that and, hence, Note that, if X max ] holds true and Eq. (41) is fulfilled for any k. Thus, similar to the approach by [23], an amplification according to Eq. (30) or Eq. (31) along with a monotonically increasing function Ψ * 0 w.r.t. the principal invariants are sufficient conditions for thermodynamic consistency.

Experiments
For illustration of the applicability of damage models, two filled elastomers possessing different mechanical behaviour are experimentally tested: a highly silica filled, sulphur-linked solution-styrene butadiene rubber (S-SBR, referred to as stiff compound, see Fig. 5a) and a lowly carbon black filled, sulphur-linked polybutadien/polyisoprene blend which tends to strain-induced crystallisation (BR/IR, referred to as soft compound, see Fig. 5b). Herein and in the following, quantities without indices refer to the loading direction.
For each material, quasi-static (λ = 0.01 1 s ) multihysteresis experiments were conducted under uniaxial, planar (pure shear) and equibiaxial tension. For this purpose, the material was loaded in a cyclic manner with  three cycles per amplitude. Here, the amplitude levels for the uniaxial case were spaced logarithmically in the stretch domain because the extend of softening increases approximately logarithmically. The amplitudes are given in Table 1. To find equivalent strain amplitudes for the planar and equibiaxial deformation which generate a similar extend of softening, the concept of the invariant radiusĪ r is introduced. For this purpose, the principle invariants are plotted in a parametric representation for each loading mode, see Fig. 6. Then, circles with the center (3, 3) are drawn in the invariant plane. The circle's radii are defined by the uniaxial amplitude levels, viz.Ī 2 r = Ī 1,ux − 3 2 + Ī 2,ux − 3 2 . The intersection with the planar and equibiaxial curves provide the corresponding amplitudes of these deformation modes. The obtained data of the multihysteresis experiments are used for the parameter fitting whose results are discussed in the following section.

Stiff compound
For the stiff compound, the first three amplitude levels of all three deformation modes were considered for fitting. The fourth amplitude levels of the uniaxial and planar tension test were used to prove the predictivity of the models. Fig. 7a and Table 2 show the ranking of all damage models sorted by the cost function value (10).
Only the best combinations regarding the hyperelastic model are considered for the ranking.
Studying the effect of the basic model in Fig. 7a, two general things can be stated: The choice of the basic hyperelastic model is related to the model class of the damage model, e.g., all models using a damage  variable provide the best results if combined with the hyperelastic approach (a) by [20]. Generally, the basic hyperelastic model (a) is a reasonable overall choice for the stiff compound. However, the virgin state variable models are less sensitive to the employed basic model and tend slightly to prefer strain energy function (b). The models 3.1a/b are insensitive to the basic model, too, whereas the combination of models 3.2/3.2* with strain energy function (b) should be avoided. A second conclusion is that those models with a virgin state variable and models 3.2/3.2* are clearly the better choice for this material. However, note that the approach 2.4 by [4] is just one part of a more sophisticated model and, for sure, yields better results taking the full model into account.
Notably is also the difference between the performance of the damage models with an amplification variable. Apparently, the simpler amplification approach 3.1a/b is not sufficient to reproduce the complex softening behaviour of the tested materials. In contrast, model 3.2 provides a good trade-off between fitting quality and number of parameters. Though, the high number of model calls needed to find optimal parameter set and their quite high correlations are a slight drawback. The reader should keep in mind that the models 3.1a/b  Table 2. In case of models using a virgin state variable, it can be explained by the provided data. For these models, the basic strain energy describes the virgin curve, see Fig. 2. Since the virgin curve of the stiff compound does not include an upturn, see Fig. 5a, all parameters related to the upturn are indeterminate to some extent.
We like to point out another property of the virgin state variable models which does not appear in the ranking. In contrast to the remaining approaches, these models provide a pretty small correlation between the parameters stemming from the basic strain energy and from the damage model itself, see also Fig. 8d. That is, there is no undesired coupling which may lead to physically unlikely or unforeseen behaviour. Moreover, the η-functions of the best ranked virgin-state variable models 1.6 and 1.4 (as well as their modifications) are defined via a square root and an exponent q < 1, respectively, leading to an infinite slope for Ψ 0 → Ψ 0,max , cf. Fig. 3b. This behaviour seems to be essential for good fitting results since it represents the steep stress decrease after the turning point between a virgin loading and an unloading. Though, a slight drawback is to be expected considering a finite element implementation that requires a material tangent, i.e., the second derivative of the free energy with respect to the deformation. Since the derivative of η is indeterminate due to the infinite slope, a special numerical treatment is needed. All fitted parameters are given in Table 4. The results of the best ranked model 1.6* are depicted in Fig. 8. One should keep in mind that each amplitude was cycled three times. Therefore, a much higher weight is given to un-and reloading curves in comparison with the virgin load.

Soft compound
Seven amplitude levels of each of the equibiaxial, planar, and uniaxial tests were considered for the fitting with the soft compound. The data of the eighth amplitude of the uniaxial and planar tension test provide an insight into the predictivity. Comparing the ranking regarding the soft compound in Fig. 7b and Table 3 with the results from the stiff compound, the much smaller RMSE of all models is conspicuous. This is related to the fact, that the soft material is dominated by strain-induced crystallization rather than the Mullins effect. As a consequence, the permanent set which cannot be reproduced by the idealized damage models, cf. Fig. 1, is less pronounced, cf. Fig. 5b. Hence, better fits will be obtained for this material. Like for the stiff compound, the virgin state variable models 1.4/1.4* [42] and 1.6/1.6* [15] combined with the strain energy function (b) by [40] achieved top rankings. They also show a quite low RMSPE indicating a good predictivity. The other virgin state variable models, except model 1.1, provide sound fitting quality, too.
In contrast to the stiff compound, strain energy function (c) is a good overall choice of the basic hyperelastic model. However, the virgin state variable models are again largely insensitive to the basic model and tend to better results if combined with strain energy function (b).
The model 3.2 by [34] in combination with basic model (c) shows also remarkable performance and is placed closely behind model 1.6. The results are depicted in Fig. 9. It is a promising alternative to the models with a virgin state variable. Note that in general models 3.2/3.2* should be preferably combined with strain energy function (c), cf. Fig. 7.
Models with a damage variable are less satisfactory again. However, in contrast to the stiff compound results, they are not far behind the models with a virgin state variable and may also be combined with basic model (c). Comparing models 2.2, 2.3, 2.5, 2.6 which use the same exponential evolution equation but different measures of maximum load, the approach 2.3 by [8] with Γ = Ī 1,max 3 − 1 seems to be the better choice. Further, the comparison between model 2.1 and 2.2 which both define Γ = Ψ 0,max but employ different evolution equations indicates that an exponential evolution as used by model 2.2 is probably not the best choice.

Conclusion
Fifteen common approaches to model the characteristic Mullins-effect of rubber materials in a discontinuous way were analysed. After discussing the possible range of applications, these models were classified by the type of their internal variable, proven to be thermodynamically consistent, and reasonable modification are proposed. Then, their ability for reproducing multihysteresis experiments of two different filled rubber compounds, referred to as stiff and soft compound, were tested. For this purpose, each damage model was combined with three different basic hyperelastic models and fitted to the experimental data. Thus, 19 × 3 × 2 = 114 fittings were conducted. Rankings for both the soft and stiff compound sorted by the fitting quality were compiled. Moreover, the root mean square prediction error (RMSPE), parameter correlations, number of parameters as well as the number of model calls needed for the fitting procedure were discussed.
For the presented materials, the damage models with a virgin state variable in combination with the hyperelastic function by [40] referred to as basic model (b) seem to be a good overall choice. In particular, the models by [15] and [42] (as well as their modifications) ranked in the top five for both materials. An alternative to the models with a virgin state variable is the approach by [34] based on the concept of strain amplification. Its combination with the simplified tube model by [35] referred to as basic model (c) provides sound fitting results. Models based on a damage variable lead to less satisfactory fits. In summary, even if several aspects of the complex mechanical material behaviour of elastomers are neglected, the considered damage models offer quick insights and a reasonable predictability for certain applications at low numerical costs. Moreover, the presented investigations may be used as a starting point for improving more complex models for rubber materials. Table 4 Fitted parameters of the damage models for the stiff compound (see Table 6 for the parameter units)
Parameters Mod. Parameters   Table 5 Fitted parameters of the damage models for the soft compound (see Table 6 for the parameter units)