Material parametrization of natural rubber based compound and characterization of crack propagation under consideration of dissipative effects

Most concepts to characterize crack propagation were developed for elastic materials. When applying these methods to elastomers, the question is how the inherent energy dissipation of the material affects the cracking behavior. This contribution presents a numerical analysis of crack growth in natural rubber taking energy dissipation due to the visco-elastic material behavior into account. For this purpose, experimental tests were first carried out under different load conditions to parameterize a Prony series as well as a Bergström–Boyce model with the results. The parameterized Prony series was then used to perform numerical investigations with respect to the cracking behavior. Using the FE-software system ANSYS and the concept of material forces, the influence and proportion of the dissipative components were discussed.


Introduction
Due to their mechanical properties, especially their high deformability, elastomers have a wide range of applications. The design of components requires an understanding of failure and wear, which in turn requires an understanding of crack initiation and propagation. Nowadays there are several well-established concepts to characterize cracks in the context of brittle fracture like the J-integral [5,23]. However, when applied to elastomers, the question of the influence of dissipative material behavior on this concept arises. In the context of the present work, material forces are used to describe the influence of the dissipative material behavior on the crack. The material behavior of elastomers is presented by through constitutive models, which includes a suitable visco-elastic characterization. This material description is finally used to describe crack propagation experiments by material forces and to verify the applicability of the extended approach.
The foundation of material or configurational forces is the work of Eshelby [10], in which the elastic energy-momentum tensor is introduced to describe the force on defects and inhomogeneities in an elastic continuum. In a later work [9], he recognizes that the material forces on a crack in a purely elastic body under quasi-static load correspond to the well-known J-integral. The concept of material forces offers a more general interpretation of the J-integral and thus a unified approach to the description of cracks in complex materials such as hyperelastic [29], viscoelastic/elastoplastic [19,20] ,or magnetoelasticity [12].
From a phenomenological point of view, polymers exhibit both non-linear finite elastic and inelastic behavior including deformation induced damage effects. Describing this behavior with a constitutive model has been an active research area for many years. Though, most of the approaches can only account for some of these effects. Experimental observations, in which elastomers have been loaded by a step in strain, have shown that the stress relaxes towards a state of equilibrium (e.g. [11]). Based on such tests, the mechanical behavior of elastomers under infinitely slow loading is described using hyperelastic material models. There are Dr Michael Johlitz, Guest Editor of this Special Issue [New Challenges and Methods in Experimental Investigations and Modelling of Elastomers] confirms that where he was co-author of a research article in this Special Issue, he was not involved in either the peer review or the decision-making process of that particular article.

3
numerous different approaches for the formulation of such models, which are compared for example in [17].
Polymers under cyclic loading show that the non-equilibrium behavior is dependent on both rate and temperature [16,22]. A widely used way to model the time dependence of relaxation is to use dimensionless relaxation modules, which are usually represented by a Prony series (e.g. [13]). This approach shows good results, especially for small strains. Another approach is to decompose the deformation gradient multiplicatively into an elastic and a viscous part, which is used for example by the Bergström-Boyce model (BB model) [4]. Compared to the Prony series, this model is not based on a single hereditary integral with stress-like internal variables, but the inelastic part of the deformation represents strain-like tensorial state variables. This approach uses a non-linear evolution equation for the strain-like variable to better model the relaxation in the finite strain regime. This paper is organized as follows: In the second section, the mechanical behavior of carbon black-filled rubber is investigated. For this purpose, the BB model and the Prony are introduced and the experimental methodology is presented. Subsequently, the equilibrium and viscous behavior are analyzed separately. Then, in the third section, the cracking behavior of rubber is investigated, first introducing the methodology for determining the material forces and then presenting a numerical example for crack characterization. Finally, to sum up, a conclusion is drawn in section "Conclusions".

Continuum mechanical description of the BB model
Each body B consists of a set of material particles and is described through its configurations i . Then let ⃗ X be a one-to-one correspondence between a Particle and the Point that B occupies within 0 at reference time t 0 . The function maps all Points ⃗ X of this reference configuration onto points ⃗ x = ( ⃗ X, t) of the current configuration. To describe the transition between the reference and the instantaneous configuration, the deformation gradient is introduced by = ⃗ x∕ ⃗ X and the volume change J, which is defined by J = det( ) . Furthermore, the eigenvalues of the deformation gradient i in the main strain space describe the elongation of the body B . The BB model is based on a multiplicative decomposition of the deformation gradient, which is based on experimental observations. As a consequence, the simplified material behavior of a polymer can be modeled as two networks, as shown in Fig. 1. The elastic equilibrium behavior is modeled by a single spring in A, whereas the spring and damper in B represent the time-dependent viscous component. Thus, an applied deformation affects both networks in equal measure, while the deformation in B can be decomposed by the relationship into the deformation of the damper v and the spring e . With the viscous part v , an intermediate configuration is introduced, in which virtual unloading of the elastic part in B takes place and thus a stress-free relaxed state prevails. If one compares the BB model with one-term Prony series, the difference becomes apparent in the description of the damper. In the Prony series, the strain rate in the damper is described as a linear function of the Stress in B. Whereas in the BB model the rate is based on a non-linear relationship, which is derived from a Reptition-type tube theory.
An important finding from this micromechanically inspired idea is that the viscous flow is energetically activated, as the tube-like conformation areas of the polymer chains act as an energy barrier. Consequently, the relaxation process is described through which defines the effective flow rate in network B. The first part of the Eq. (2) describes the dependence on the actual chain elongation via the material parameters c and ̇0 as well as the principle macroscopic stretch state chain . The principal macroscopic stretch state is based on the 8-chain assumption [1] and defined by chain = √ I 1 ∕3 . The second part of the Eq. (2) describes the connection of ̇B to the applied stress, which is assumed to be energetically activated. This relation is modeled through a generically expressed power-law, depending on the material parameters ̂ and m and an effective stress measure ‖ dev ( B )‖ F , where ‖‖ F is the Frobenius norm and dev ( ) is deviatoric part of the Cauchy stress tensor. N() is a function to provide the direction of the viscoelastic flow depending on the viscous stress in B. For a more detailed explanation of the BB model, One dimensional rheological representation of the BB model especially with regard to the numerical implementation, the reader is referred to [6]. For a constitutive description of the mechanical behavior, a function of the Helmholtz free-energy density is introduced which can be split into a volumetric and isochoric part. The isochoric part of the deformation is characterized by where the scalar can be interpreted as hydrostatic pressure. For the isochoric response of the material, i can be defined by a phenomenological approach like the Ogden model [21] where p and p are material parameters and n the number of terms. A micromechanically inspired approach is the Extended Tube Model [14] where and are material parameters, describing the finite extension of the polymer and global rearrangements of cross-links upon deformation, Ĩ 1 =̃2 1 +̃2 2 +̃2 3 is the first invariant of the isochoric part of the Cauchy-Green tensor, G c the shear modulus characterized by the cross-linking of the network, and G e the shear modulus defined by the confining tube constrains.
To reproduce the viscous material behavior it is useful to split the isochoric free energy further by into an elastic i,∞ and a viscous i,v part. The elastic part characterizes the equilibrium state. Where the viscous part defines a dissipative potential, describing a non-equilibrium state. By deriving Eq. (7) with respect to the deformation gradient a relationship can be found for the spherical v and deviatoric i,∞ part of the 1-Piola-Kirchhoff stress. Furthermore, the quantity i,v can be introduced, which can be interpreted as non-equilibrium stress.

Experimental procedure
To model complex material behavior correctly tensile tests under uniaxial (UA, 1 ; 2 = 3 = 1∕ √ 1 ), plane strain (PS, 1 ; 2 = 1; 3 = 1∕ 1 ) and equibiaxial (EB, 1 = 2 ; 3 = 1∕ 2 1 ) loads were carried out to parametrize the material models. The PS and EB tests were carried out on the biaxial test machine, built by Coesfeld GmbH & Co, Dortmund, Germany, which is described in more detail by [25]. Besides, uniaxial tensile tests were carried out on an universal tension-testing machine Z010 of ZwickRoell, Ulm, Germany. A NR20 (natural rubber, 20 phr carbon black) material mixture from Weber & Schaer GmbH & Co. KG. was used. The compounding recipe is shown in Table 1. The samples were vulcanized with sulfur as a crosslinker for 10 minutes at 150 • C . For the PS and EB test, a square specimen of 77 x 77 x 1.8 mm 3 with a bulge (diameter 5.5 mm ) was used. For the UA tests, a dumbbell specimen with a length of 75 mm was stamped out of the square specimen after preconditioning.
In preliminary tests, the strain fields of the test specimens were determined through an optical method of grey-value correlation, using the ARAMIS software from GOM GmbH, Braunschweig. Thus, a second-order polynomial, describing the mean strain as a function of displacement, was used to evaluate the experimental strains. The material parameters are determined in two stages: For the purely elastic characterization, the test specimen is loaded step by step with relaxation times in between to a maximum strain and then also unloaded again, as shown in Fig. 2. In a second test, the rate-dependent characterization is done by series of tests, each with one load cycle of different maximum strains and strain rates, which are listed in Table 2. The determination of the material parameters represents a non-linear optimization problem. To solve this problem, a suitable measure for the quality of the fit and a non-linear optimizer is required. All experiments were strain-controlled, thus the difference between simulated and experimental stresses was used to determine the error. A common measure of error is the normalized mean absolute deviation ("NMAD"), which is estimated according to [15] by This error measure is considered robust and can be used for the comparison of several tests with different strain and stress amplitudes due to its standardization. Within this work, the "Differential-Evolution" and "Basin-Hopping" algorithms from the Python-based software SciPy are used as optimizers. The stress for the simulated material behavior for the purely elastic behavior is determined in an own Python routine. The cyclic experiments were simulated in a FE analysis with a single second-order hexahedral element and corresponding boundary conditions.

Equilibrium behavior
The results of the step test in Fig. 2 show that the stress on the loading path relaxes faster than on the unloading path. In [16] this phenomenon is interpreted as a proof that the equilibrium behavior can only be represented by a hysteresis. Another explanation is given in [2]. The relaxing parts of the polymer on the load path are subjected to tensile forces, whereas those on the unloading path are under compressive load. This causes a fundamentally different behavior. Another aspect is the different distance to the equilibrium state. Although the total elongation at each stage is the same due to the same external displacement, this does not imply the same loading of the viscous parts, since the material cannot completely relax at any stage. Experimental confirmation of one explanation seems possible only by infinitely long relaxations test, which is not possible. However, the following investigation aims to determine the most suitable test parameters to determine the equilibrium curve. The equilibrium stress is in the range between the stresses after the steps with the same strain. Determining the value from the arithmetic mean leads to an overestimation of the equilibrium, because the stress on the loading path relaxes faster. Suitable test parameters lead to a small difference in stress P and therefore to the smaller error in determining the equilibrium points. The studied parameters are t D and t r , which accordingly describe the finite time to apply the deformation and the relaxation time on each load step. The Fig. 3 shows the influence of the loading time on the different stresses from loading and unloading with the same strain magnitude of a PS tensile test.
The Fig. 3 also shows, that the difference decreases with increasing strain. This is due to the fact that the relaxation process is energetically activated. Because at higher strains there are also significantly higher stresses, the material relaxes much faster towards equilibrium. Furthermore, the figure also shows that a decreasing loading time leads to a strong decrease in the difference. This becomes particularly clear comparing test with t D = 1 s and t D = 10 s . An analogous comparison of relaxation times of 1 min and 10 min was done. The difference in stresses decreases with increasing relaxation time. However, the minimization of the difference and thus of the error does not seem to be justified by the significantly higher experimental effort. For the determination of the equilibrium stress, it is, therefore, appropriate to carry out step tests with t r > 60 s and t D > 10 s , which is equivalent to stretch rates of ̇< 0.01 s −1 .

Parameter identification of the material model
The step tests were carried out and the corresponding points on the equilibrium curve were determined. These points were used to parameterize both an extended Tube model and a second-order Ogden model, which lead to the material parameters in Table 3. The points are shown together with the result of the extended Tube model in Fig. 4. The extended tube model has the smallest normalized error of NMAD = 0.054 . Especially the behavior during uniaxial and equibiaxial tension can be well reproduced by the model. On average, the predicted stresses deviate less than 0.02 MPa from those in the experiments. On the other hand, the behavior in plane tension can be described much worse, resulting in a mean error of 0.131 MPa . The results of the Ogden model with two parameters show a similar quality with NMAD = 0.057. The cyclic tests were carried out and the results were used to parameterize a BB model using the Ogden model as hyperelastic part in the FE software ABAQUS. In comparison to this, a Prony series with 10 links based on the extended tube model was characterized in ANSYS. The error for both models is NMAD = 0.096 . A comparison of some experiments with the corresponding simulated results is shown in the Fig. 5. The corresponding material parameters are listed in Table 4. The illustration shows that viscous behavior can be represented quite well in some load cases like Fig. 5a, b whereas the other two cases are not described sufficiently.
One reason for the deviations between the calculated and the experimental stress curves is the path dependent relaxation. As described in the last section, the relaxation is faster on the loading path than during unloading. The models cannot represent this behavior properly. A further problem is the dependence of the relaxation on the current stretch. Again the Prony series cannot model this. The BB model describes

Material forces
The starting point for this approach is the Eshelby stress tensor or energy momentum tensor of elastostatics, which is defined as where is the Kronecker delta and ∇⃗ u is the gradient of the displacement. The sink in the field of at a crack tip can be interpreted as a driving force for the crack growth and is called material force. Material forces differ from the "normal" forces in that they are described by changes in material and not in the actual physical reference system, see [18] for an extensive explanation. In a two-dimensional elastic body with a crack, a shrinking path S around the tip is defined, like in Fig. 6. Applying the divergence theorem the material forces ⃗ R can be calculated by where r describes the radius of the shrinking path S and ⃗ n the normal vector to the path. For a general calculation of ⃗ R , a closed integration path via an arbitrary path , including the crack edges + and − , will be defined. This contour characterizes the defect-free material area A , which does not contain the crack tip. For this path, the Eq. (11) is extended by The surface integral in Eq. (12) is called dissipation force and defines the inelastic body force. With purely elastic material behavior and straight crack growth, the contour integral over the edges and the surface integral vanish, and the Eq. (12) is simplified to a path independent vectorial J-integral. A prerequisite for the use of this concept is that the constitutive material laws are derived from the Helmholtz free energy and that a change in non-elastic strains is accompanied by a dissipative potential [19].

Numerical example
For the numerical investigation of cracks, the tests on the biaxial testing machine were simulated. The discretization for the FE model is shown in Fig. 7, which uses plane  quadrilateral elements with a second order shape functions. The crack tip was modeled by collapsed quadrilateral elements.
The material forces were calculation in an own usersubroutine. For this, the line integral in Eq. (11) was converted into a surface integral via a virtual crack extension as described in [26]. All required parameters were determined or interpolated at the integration points. The integration was done numerically via the Gaussian quadrature.
In the first step, the relation between material forces in the direction of crack growth and the radius of the line integrals is analyzed. A displacement was applied to the upper side of the model to stretch the sample outside the crack region up to y = 1.25 . The strain rate is initially set to ̇y = 0.5 s −1 . The material is described by the extended-tube model extended by a prony series, which was defined in section" Parameter identification of the material model". Figure 8 shows the results of this investigation.
The graph shows a local maximum point with 8.3 N∕mm at the first element row around the crack tip and decreases to a value of 7.5 N∕mm at r = 0.3 mm . This phenomenon has been described in the literature with elasto-plastic material, see [20,27]. The area with the same strain magnitude around the crack tip is kidney-shaped. However, the evaluation of the material forces are carried out by circular line integrals, which is going through areas with different strain rates and thus also different material stiffness. From a distance of r = 0.3 mm , R x begins to rise to a value of 9.2 N∕mm at the very edge of the specimen. This increase is due to the viscous behavior of the material, converting energy independent of crack propagation. In addition, the rate of change decreases with distance to the crack tip, since the strain rates and thus also the dissipation becomes smaller. This area acts as a material shield, since parts of the energy introduced from the outside do not flow further to the crack, but dissipate on the surface. The phenomenon of material shielding or reverse anti-shielding is described in [28].
In a second step, three quantities and their dependencies on stretch rate are compared: R x,T is calculated by a path of the first element row around the crack, R x,F is determined by a path along the edge of the sample and the vectorial J-integral J x . R x,T characterizes the driving force of the crack growth. R x,F consists of both the driving and the global dissipation forces. Therefore, this quantity includes all dissipative processes even those far away from the crack tip, analogous to the tearing energy [24]. These values are compared to the J x , which was calculated with a purely elastic material without Prony extension. J x characterizes the crack behavior under infinitely slow loading rates. Again the model was loaded by a displacement with the same magnitude of y = 1.25 , but at different rates. Figure 9 shows that J x underestimates the driving force R x,T . At very low strain rates such as ̇y = 0.005 s −1 , this error is only 17% . However, at higher strain rates like ̇y = 2.5 s −1 , the error increases up to 40% . Furthermore, this approach fails to examine the influence of the load rates. On the other hand, R x,F overestimates the value of R x,T . At ̇y = 2.5 s −1 this error is 7% but increases to 30% at ̇y = 0.005 s −1 .
The Fig. 9 also shows that the material forces acting directly on the crack tip are of the same order of magnitude as those acting in the region far from the crack tip. This is again an indication that under the investigated boundary conditions the dissipative effects for the crack propagation are of the same order of magnitude as the viscous effects far from the crack tip. Additionally, R x,F depends on the sample size. Therefore, for an exact description of the crack behavior of elastomers, these energetic parts must be considered separately.
Both R x,T and R x,F increase with the loading rates. For the investigated range the relation between the two quantities and the logarithm of the stretch rate is approximately linear. The increase of R x,F is due to the fact that the stiffness of the material increases with increasing rate. If the test specimen is loaded with the same strain amplitude but at different strain rates, the external forces do more work, which in turn leads to an increase in material forces. Compared to R x,F , the value of R x,T increases much faster. This finding suggests that at higher load rates a higher proportion of the work is done independently of crack growth compared to lower rates.

Conclusions
In this contribution, a BB model as well as a Prony series were parameterized based on experiments to describe the rate-dependent viscous behavior of natural rubber. To model even complex loading conditions not only UA and PS, but also EB tests were carried out. The models show some difficulties in correctly characterizing the dissipative behavior, for which reasons were discussed.
Common methods for crack characterization such as the J-integral or the tearing energy fail to include the influence of dissipative effects in inelastic materials or to isolate them from effects far away from the crack tip. Therefore, here the concept of material forces and dissipation rate was used to separate these energetic quantities within a framework of finite element analysis. In the next step, this separation has to be performed experimentally and compared with the numerical results through thermographic measurements like [7,8].
Funding Open Access funding enabled and organized by Projekt DEAL.

Conflict of interest
The authors have no conflicts of interest to declare that are relevant to the content of this article.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/. Fig. 9 Material forces at the crack tip R x,T , at the edge of the specimen R x,F and J x for a hyperelastic material depending on the rate of the load