Modelling crack propagation during relaxation of viscoelastic material

A viscoelastic material which is stretched and is then held at constant elongation, normally results in decreasing stresses till the equilibrium has been reached. With the decreasing stresses a crack propagation is not expected as the energy of the system is decreasing. However, an initial damage could lead to an increase in the mechanical load on the undamaged chains during relaxation, leading to material degradation and crack propagation. While experimental investigations have been presented in the literature, modelling such an effect has not been thoroughly investigated. In this work, an initial framework for modelling the damage evolution during relaxation is presented. A mechanical model is coupled with a phase field to model the crack propagation. For simplicity, a linear viscoelastic model is implemented for the mechanical part. A mobility constant is employed to model the evolution of the phase field with the changing mechanical energy during relaxation. The evolution of phase field can be interpreted as the evolution with which the polymer chains get damaged. Different load conditions and geometries are simulated, which shows that the proposed framework is able to model the damage evolution during viscoelastic relaxation. Thus, with the help of the numerical model a physical explanation for the failure during relaxation is presented.


GRAPHICAL ABSTRACT Introduction
Polymers such as polyurethane and rubber are widely used as adhesives or sealants. These components are normally subjected to a constant deformation during their lifetime. Such loading conditions lead to relaxation of these viscoelastic materials. Therefore, a rupture of these components is ruled out from the theoretical point of view if the maximum stress applied is less than the ultimate strength of the material. However, in the literature, various experiments have shown that material rupture can occur even at a constant elongation and decreasing stress [1][2][3][4].
The phenomenon was observed in the early 40's by Tobolsky et al. [1] during the relaxation of Hevea rubber (natural rubber) where the stress decayed to 0 MPa during relaxation experiments at higher temperatures. It was theorised that the high temperature assists oxidation of the rubber which leads to bond breakage and the drop in the stresses. However, no change in the behaviour was noticed between the experiments conducted in presence of ordinary air atmosphere and in commercial nitrogen. Smith [5] measured the ultimate tensile strength for unfilled Styrene-Butadiene Rubber (SBR) at different temperatures and strain rates under tensile loads. With the help of the experimental data, a failure envelope was presented ( Fig. 1) that represented the stress-strain coordinate at which failure occurs for different strain rates and temperatures. According to the failure envelope if a specimen is elongated up to a certain strain before the failure point and then held at constant elongation (point b in Fig. 1), the stresses will relax to the failure point of a slower strain rate and lead to rupture of the specimen. For smaller strains the relaxation may lead to the equilibrium position (point a in Fig. 1). The concept was validated with several relaxation experiments conducted at various temperatures and strain rates, where the rupture of the specimens at constant elongation was noticed [2]. With the help of the time-temperature superposition principle, various failure envelopes were superimposed to obtain an envelop independent of the experimental conditions [3].
In much recent investigations, Neuhaus et al. [4] investigated the rupture of cross-linked polyurethane at constant elongations. The effect was not only reproducible in relaxation experiments, but also in cyclic experiments, where the specimens were unloaded after certain time of relaxation and then again subjected to relaxation. Although a statistical fluctuation with respect to the time of rupture can be seen in the experiments, but the stress level for the rupture was just outside the failure envelope for polyurethane. Friedrich [6] presented similar results for ethylene-propylene-diene monomer (EPDM) rubber where the time to rupture was dependent on the temperature, the strain rate and the maximum applied strain before start of relaxation. Failure during relaxation was limited in the range between two extreme points. At strains below the minimum strain, there was no failure observed in the material and at strains above the maximum point, the failure was observed in the loading phase instead of the relaxation phase. The concept of delayed fracture in polymer gels has also been discussed in the literature [7][8][9][10][11]. At a constant load the polymer gels don't J Mater Sci (2023) 58:6254-6266 show any signs of failure at a macroscopic level until crack initiation and then it leads to a sudden failure.
The classical fracture mechanics models presented by Griffith [12] and Irwin [13] predict crack propagation if the energy release rate reaches a critical value. To model the fracture with the finite element method (FEM), the path of the crack results in discontinuities in the mesh. These are modelled with the help of methods such as the extended finite element method (XFEM) [14] or stable generalised finite element [15]. The use of a phase field model to describe damage evolution has gained popularity in the last few years [16][17][18][19] as the problem of the discontinuity in the mesh is circumvented by introducing a field variable / that defines a diffused crack zone around the fracture. The variational methods for energy minimisation as introduced by Francfort et al. [20] and Bourdin et al. [21] are used to calculate the evolution of the phase field variable and in turn the crack propagation. Therefore, a fixed crack path does not need to be defined. Alternatively a Ginzburg-Landau type evolution equation [22] which is alternatively called the Allen-Cahn equation [23,24] has also been used to model the crack propagation for rate dependent dynamic fractures using the phase field model [25][26][27].
In this work, an evolution equation for the phase field is implemented to model the failure of the material during relaxation tests as noticed in the experiments. It is assumed that some polymer chains of shorter lengths break during the loading of the specimen on a microscopic level, and due to this other polymer chains carry more stress. During relaxation, the remaining short chains relax faster which leads to an increased loading of the longer chains. This eventually leads to reaching the critical load for the longer chains and a macroscopic failure during relaxation. The energy of the load carrying parts of the system continues to increase during the relaxation phase till it reaches the critical value sufficient for crack initiation and propagation. The increase in energy is modelled in the model with the increasing crack energy density with the evolution of the phase field variable. The mechanical part is modelled as a linear visco-elastic material extended by a degradation function and the evolution of the fracture energy during the relaxation phase is driven by the Ginzburg-Landau evolution equation.

Modelling
An introduction to the mechanical model and the phase field model is presented. The free energy density for both the parts is added to get the total free energy density and to derive the model equations.   [4] showing the damage of a cross-linked polyurethane system during relaxation.

Mechanical model
A rheological model is used to represent the linear viscoelastic behaviour of the material. A spring representing the basic elasticity is connected in parallel with a Maxwell element (Fig. 3a). In case of a small displacement u, the strain can be split additively into elastic and inelastic parts in the Maxwell element, with the elastic part e e acting on the spring and the inelastic strain e i stretching the dashpot. The inelastic strain evolves according to where s represents the relaxation time of the Maxwell element. With the evolution of the inelastic strain, the stress in the Maxwell element changes with time, thus representing the non-equilibrium part w neq of the mechanical energy density. At a constant strain, the stress in the Maxwell element decays to zero and the system reaches the equilibrium position. The equilibrium part of the mechanical energy density w eq is thus represented by the spring element connected in parallel to the Maxwell element. The total mechanical energy density is given by These energy densities are given by Hooke's law as and where k, l are the Lamé parameters for the equilibrium spring and k 1 , l 1 represent the Lamé parameters for the Maxwell element.

Crack propagation model
A regularised field given by the phase field variable / is used to model the diffused crack instead of a strong discontinuity. A value of / ¼ 1 defines the fully intact material and / ¼ 0 defines the crack in the material 1 .
The crack surface C (Fig. 3b) is regularised to a diffused crack volume which is achieved by the integration of a regularisation function cð/; grad /Þ The regularisation function which can also be defined as the crack surface density function as proposed by Miehe et al. [18,19] is given as The transition of the phase field from zero to one over the crack width l defines the surface of the crack. 1 Depending on the model it is also possible to model the crack with / ¼ 1 and the fully intact material with / ¼ 0, which would fit the ranges of the canonical damage parameters introduced by Kachanov oder Lemaitre The energy needed to create the crack surface is given by where G c is the Griffith-type critical energy release rate and gives the crack energy density.

Coupling
The effect of the phase field variable on the mechanical energy density is modelled with the help of a degradation function where j ! 0 is a stability parameter to prevent zero mechanical energy at the crack. A quadratic degradation function has been used in a general way. However a comprehensive study on the choice of the degradation function can be found in [28]. With decreasing value of /, the degradation function decreases and accordingly, the mechanical energy density is degraded, modelling the crack in the material. The total energy density is given by the contribution of the mechanical energy density degraded with the phase field and the crack surface energy density. This leads to the stress as a function of the strain and the degradation function. The evolution of the phase field is modelled according to the Ginzburg-Landau type evolution equation where M ! 0 is the mobility constant driving the rate of evolution and d u ðw t Þ denotes the variational of the total energy. Since the phase field evolves from / ¼ 1 to / ¼ 0, the evolution of the variable is negative. In this way, crack healing is prevented. Moreover, if / reaches negative values, then it is artificially corrected to the value of zero to maintain / 2 ½0; 1. By expanding equation (14) _ it can be seen that the evolution of / is coupled with the mechanical energy density and hence any mechanical deformation or change in stress can cause an evolution of the phase field.

Numerical implementation
The material model is implemented using the finite element method (FEM) to simulate different load conditions. The equilibrium condition div r ¼ 0 ð16Þ is solved for the mechanical equation where r is given by the equation (13). For the phase field model equation (15) is solved. The Newton method is used to solve both these coupled equations simultaneously. The weak form of the mechanical equation (16) gives the residual Here t is the traction acting on the surface of the body which is handled as a Neumann boundary condition. Similarly the residual related to the phase field is given by The boundary integral which represents the flux of the phase field is omitted as the flux is assumed ot be zero on the boundary. In equation (17), u u represents the test function for the displacement field u and in equation (18), u / represents the test function for the phase field / in the finite element implementation.
The time integration of _ / is achieved with the help of the Crank-Nicolson method resulting in The value of h ¼ 0:5 is taken to achieve high accuracy and stability. The quantity ½h t is the known quantity at the time step t and the quantity ½h tþ1 represents the unknown quantity at the time step of t þ Dt. Similarly equation (3) for the evolution of the inelastic strain is also calculated using the Crank-Nicolson method. By rearranging the terms, an equation for the unknown e tþ1 i can be formulated in terms of the current strain e tþ1 and substituted in equation (13). After solving the coupled equations, the value of the inelastic strain is stored locally as an internal variable which provides the value for e t i for the next time step. The global residual equations are transformed into the local form using the shape function interpolations. The linearised system of equation Du is solved with the Newton method for the coupled problem. The derivatives of the two residuals with respect to the two solution fields are calculated analytically. The coupled differential equations are implemented using the finite element libraries of deal.II [29,30].

Simulation results
Simulations are conducted on two different geometries to simulate crack propagation with and without an initial crack in the specimen. The results and their physical interpretation are presented in the following subsections.
Specimen with pre-existing cracks As a first example, a square of sides 1 mm with a slit half-way through its width is fixed at one face and loaded on the opposite face in a direction perpendicular to the direction of the face till failure (Fig. 4a).
The dimension of the specimen is irrelevant and is chosen randomly to be a unit cell. Similarly, the parameters for the material are chosen arbitrarily and are given in Tables 1 and 2. Five different strain rates were used to simulate the tensile test. The stressstrain curves for all the rates show an increase in the stress up to a particular strain, following which the stresses decrease with increasing strain, depicting crack initiation and propagation (Fig. 4b). The point of maximum stress is marked as it denotes the failure point for the material for that strain rate. These points can be joined to form the failure envelope as suggested by Smith et al. [5]. The equilibrium stressstrain curve is calculated by neglecting the Maxwell element and measuring the stress-strain response for just the equilibrium spring. After the failure point the crack propagates parallel to the slit during loading (Fig. 5) until complete fracture of the specimen. The parameters for the mechanical model that are given in Table 1 are selected so that the basic elasticity element has higher l than k, whereas the Maxwell element has a higher k. This produces a different volumetric and shear stress components in the specimen depending on the strain rate. At higher strain rates the influence of the Maxwell element is activated and hence along with the volumetric part the shear part of the stress also contributes to crack growth. In contrast, for slower strain rates the Maxwell element is not activated and only the volumetric part contributes to the crack growth. This results in different crack shapes. The phase field distribution at the time of complete fracture shows that the crack is wider for a higher strain rate of 10 s À1 (Fig. 5a) whereas it is narrower for the lower strain rate of 10 À4 s À1 (Fig. 5a).

Relaxation test on specimen with preexisting crack
Numerical simulations are conducted to find out the effectiveness of the model to reproduce the damage during relaxation tests. The simulations are conducted with the geometry given in Fig. 4a. Strain is applied at the rate of 10 s À1 up to a certain value and is then held constant. With the constant strain condition the change in the stress during relaxation decreases the mechanical free energy. However, the already intiated crack drives the evolution of the phase field. This results in the desired effect of crack propagation during the relaxation phase. However, a certain initial strain needs to be applied to observe the crack propagation. Three different levels of strains before the start of relaxation are simulated (Fig. 6a). When the strain is higher than that of the failure point of the equilibrium curve, then the crack propagates and the stress decays to 0 MPa. However, when the strain is below the failure point for the equilibrium curve, then the stress relaxes to the equilibrium stress (Fig. 6b).  (a) (b) Figure 5 The crack propagates parallel to the slit in the specimen. The crack is much wider for the strain rate of a) 10 s À1 , where as it is narrower for b) a strain rate of 10 À4 s À1 .
The point of intersection of the stress-strain relaxation curve for the highest initial strain and the failure envelope is marked with the symbol 'Â' in Fig. 6b. For the corresponding point in Fig. 6a, a change in the rate of stress decay can be noticed. This indicates that the failure began at this point of time, which validates the failure envelope. If the crack length is taken as a measure of failure, then it can be seen in Fig. 7 that the point of crack initiation lies outside the failure envelope. With further increase in time the crack propagates and the stress reduces to zero.

Specimen without pre-existing crack
Since no slit or initial crack was present in the experimental specimens, the model should be able to simulate crack propagation without a slit in the geometry. To this end, an hour-glass geometry with a reduced cross section in the middle is simulated. The geometry resembles the mid-section of a tensile test sample. The boundary conditions are similar to the previous geometry but for faster convergence of the simulation, the horizontal component of the displacement is not constrained (Fig. 8). The parameters for the material model are as given in Table 1 and  Table 2. However, to investigate the interaction between the dynamics of crack propagation with the dynamics of viscoelastic relaxation, the ratio of the relaxation time of the Maxwell element s and the mobility constant M is varied. There are two evolution equations that are solved simultaneously in time. The evolution equation for the phase field model is driven by the mobility constant, whereas the evolution of the inelastic strain is controlled by the (a) (b) Figure 6 Relaxation tests carried out with the strain rate 10 s À1 plotted on a) stresstime plot and on b) stressstrain plot along with the failure envelope. relaxation time of the Maxwell element. With the interaction between these two parameters, the time for crack propagation during relaxation can be controlled. Eight different ratios of s=M are tested and in order to compare their results, the initial strain is kept constant for all ratios at 0.2 before the start of relaxation. Along with the different s=M ratios, two different values of mobility constant are used.
The change in s leads naturally to change in the relaxation time and as can be seen in Fig. 9, the relaxation time is almost non-existent for the ratio of s=M ¼ 10 À7 for the orange curve and increases till the ratio of s=M ¼ 10 À4 . At higher relaxation times (for s=M [ 10 À4 ), the stress that is generated during the initial loading is high enough to initiate the crack. Therefore, the stress in the material starts to degrade before the start of relaxation and the crack propagates during the relaxation phase. In comparison to M ¼ 50 the crack initiates quite early for M ¼ 500. The higher mobility constant accelerates the phase field evolution and leads to a faster failure. Thus, to achieve the varying time to material failure during relaxation experiments, the parameters M and s and their ratio are instrumental.
Apart from that, the value of the critical energy release rate G c also influences the time to failure. With a s=M ratio of 10 À4 and with M ¼ 50 the relaxation simulation was conducted with a smaller G c of 0.002 Jmm À1 . It can be seen in Fig. 10 that the loading phases remains identical for the two values of G c . However, during the relaxation phase, the crack energy increases with the parameter G c and hence the lower value leads to a faster crack propagation.

Comparison with polyurethane specimen
The model is applied to simulate the experimental results given in Fig. 2 [4]. The experiments were conducted on a polyurethane (PU90/10) specimen with the geometry given in Fig. 11.
The geometry is a modified form of a general tensile test specimen. It has a curvature in the middle so that the crack can initiate. The mesh at this narrow area is finer in comparison to the mesh in the outer areas for accurate and efficient numerical results. The material parameter such as the stiffnesses and relaxation time was determined with the help of the longest available relaxation data and are listed in Table 3.
The relaxation test was conducted at a lower elongation of 10 mm at a speed of 0.1 mm/s. With the same boundary conditions, simulations were conducted and compared to the experimental results. The parameters were optimised using a Nelder-mead algorithm [31]. The results for the best fit can be seen in Fig. 12(a). It should be noted that the simulation has been carried out with just one Maxwell element. Hence the experiment for longer durations cannot be adequately captured. However, for modelling the failure at higher initial loads, such long relaxation time would not be needed.
To determine the fracture related variables more experiments are needed, specially in tensile directions. The fracture toughness is therefore assumed to be G c ¼ 0:02 J/mm as per the value available in literature for other cross-linked polymers [32]. To apply the suggested model for a real material, the effect of chain breakage with time needs to be captured in the model. With increasing damage at the microscopic level, the mobility of the chains increases, which results in faster relaxation and faster failure of the material. This is modelled numerically by the evolution of the mobility constant M with the increasing phase-field. Thus, is substituted in equation (14). A value of M o ¼ 2:0 Â 10 À5 was used to simulate the failure during relaxation. The power n ¼ 6 was used to match the time of failure during relaxation. These parameters were determined using the experimental data given in Fig. 2. Therefore, these are not a validated set of parameters. They are used to demonstrate the ability of the model to reproduce experimental results. The experiment in Fig. 2 were conducted with an initial deformation of 50 mm with 0.1 mm/s deformation rate. With the same boundary conditions as the experiment and with the identified parameters, the model was simulated. Considering the simplicity of the model, the simulation result shows an acceptable similarity with the experimental results as can be seen in Fig. 12 (b).

Conclusion
Experiments showing damage during the relaxation of a viscoelastic material have been presented in the literature by various authors. A probable explanation for such a behaviour is an initial breaking of some polymer chains during the loading of the specimen, which leads to higher stresses on remaining polymer chains causing them to break during the relaxation phase. The concept of failure envelope was proposed Figure 10 The influence of G c on the crack propagation.  by Smith et al. [2] to explain the failure of the material at constant elongation and relaxing stresses. However, a more detailed attempt to model this phenomenon is not to be found in the literature up to now.
In this work, a material model coupling the viscoelastic behaviour with the crack propagation is presented. A rheological model with a Maxwell element connected in parallel with a spring element is used to model the viscoelasticity. For simplicity, the strain is split additively in the elastic and the inelastic part of the strain, and the stress response is calculated using Hooke's law. An evolution equation for the inelastic strain is responsible for the time dependent response of the material. Along with the viscoelastic part, a phase field model is proposed to model the crack propagation in the material. The evolution of the phase field is based on the Ginzburg-Landau type equation, which is alternatively known as the Allen-Cahn equation. The evolution leads to a change in the phase field variable during loading as well as during the relaxation of the material. A mobility constant is used to couple the evolution with the energy available in the system. The coupled model is solved numerically using the finite element method and numerical examples are presented. A failure envelope is reproduced for a given set of parameters by simulation at different strain rates. Relaxation tests for the same parameters are conducted validating the concept of the failure envelope. Further numerical examples are presented to test the effect of the ratio of the relaxation time for the Maxwell element and the mobility constant for the phase field evolution. Thus, these two parameters can be varied for different materials to reproduce the different times for failure during relaxation. Further, the experimental results on PU90/10 specimen were reproduced numerically with the help of an evolving mobility constant.
The presented model can reproduce the crack propagation during relaxation as observed in the experiments, however further improvements are required to be able to identify the parameters from an experiment and apply it for predictions. The model is based on a linear approach for small deformation which is not the ideal choice for modelling crosslinked polymers. The mechanical model should therefore be extended to geometrically non-linear deformations.

Data and code availability
The data and the code used for this work are available and can be obtained by contacting the authors.

Declarations
Conflict of interest All authors agree with the content and give their consent to submit. There are no conflicts of interest that could potentially influence or bias the submitted work.

Ethical approval No experiments involving human tissue had been carried out.
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://creativecommons.org/licen ses/by/4.0/.