I thank several people -the ones that were almost nished with their work when I was staring, the ones that were already there, the ones that started with me, the ones that started a er me.

In the rst category, I thank omas Corre for various discussions that we had at the beginning of my thesis and for a nice transition.

In the next category comes Baptiste Reyne, with whom I had the chance to have some very nice discussions about several topics in Mechanics. Tauno Tiirats introduced me to Abaqus and helped in e ciently linking the subroutines wri en in C++ with Abaqus. Not to miss in here are George Jagite and Erwan Grelier.

Next, the people who started their PhDs with me -Elie Eid, who happened to work on dynamic fracture as well and his thesis defended itself successfully the day a er mine. Raphaël Langlois (aka hungry) and Adrien Vinel, with whom I did not have several overlapping areas of interest at work, but had several common interests in comedy. Raphaël also introduced me to cheese and brought in several varieties of cheese for my a er defense. Not to forget Paris Mulye for interesting conversations about his and my works during the rst quarantine (for those people in future, yes the covid was during our thesis) and later on as well.

e people who started a er me and shared an o ce with me -Pauline Bourda, a fellow fan of e O ce, for various discussions and for nice presents on my farewell. Shaoqi Wu, for some interesting discussions on several topics.

On the other front, I'd thank Erwan Verron for our discussions regarding viscoelastic models (the details of which I couldn't unfortunately develop further) and other aspects, Julien Réthoré for his invaluable comments during my CSI meetings. Also, my supervisors Michel Coret and Nicolas Moës for giving me this opportunity and pu ing up with me for the three years.

e readers of this thesis, should they wish, can read the citations as that's what (s)he said, as Micheal Sco , the Regional Manager of Dunder Mi in, Scranton, would do. e template used in this thesis was developed by Baptiste Reyne and can be found at h ps://gitlab.com/breyne/latex-two-thirds.

General introduction

Elastomers are the materials that are capable of undergoing large deformations reversibly. Compared to bri le solids, high speed crack propagation (high speed when compared to the corresponding wave speeds) in these materials happen at large strain levels. When a membrane made of elastomer is stretched and a crack is introduced, the crack propagates with a speed that depends on the stretch level. ese types of experiments were widely performed in the literature on the specimen in pure shear geometry. e output of these experiments are the crack speeds. e crack speeds are usually correlated with the energy release rate, which is measured using the specimen geometry and the applied stretch, using a constitutive model. In some experiments performed like above, it was observed that the crack speeds exceeded the shear wave speed (computed from the hyperelastic material parameters) of the material. Such cracks are called Transonic cracks.

Linear Elastic Fracture Mechanics (LEFM) bounds the speed of the cracks from the above by the corresponding wave speeds depending on the mode of propagation. For the cracks that propagate in opening mode, the limiting speed is the Rayleigh wave speed. However, for the propagation under in-plane shear loading, the cracks were observed to propagate at speeds that exceed the shear wave speed. e theory of LEFM was later extended to include the description of such cracks. Also, during the cracks propagation in this regime, shock fronts that travel along with the crack tip were observed in the experiments. In the case of aforementioned transonic crack propagation in elastomers, no reports of such shock-fronts in the material exist. Instead, it is remarked in the literature that for the cracks that propagate in this regime, the crack faces (crack opening) are wedge-shaped as opposed to being parabolic-shaped in the subsonic regime. No speci c mention about the shock front in the material was made.

Some studies a ribute the existence of Transonic cracks in elastomers to the hyperelastic sti ening of the material.

ey hypothesize that the hyperelastic sti ening increases the wave speed locally near the tip allowing the cracks to propagate at Transonic speeds. Other studies a ribute the observation of Transonic cracks to the viscoelastic sti ening of the material. To the author's knowledge, no study exists that evaluates the two hypotheses to conclusively determine if either of these phenomena are responsible for the crack propagation in transonic regime.

In the experiments performed in literature on natural rubber, it was observed that the crack speeds became independent of the specimen geometry in the transonic regime. More precisely, it was observed that when the specimen of di erent geometries are subjected to the same stretch level, the crack speeds do not change from one specimen to the other when propagating in the transonic regime. When the crack speeds are smaller, the crack speed dependence on the specimen size was reported to be seen. Whether this observation holds for all the materials is not clear. Also, to the author's knowledge, no studies exist that examine the necessary ingredients to investigate this phenomenon using the Finite element method.

e energy partition during the fracture of elastomers is also of interest. Some part of the energy stored in the body is dissipated by the viscoelastic e ects of the material and some at the tip, in the fracture process zone. However, to the author's knowledge, how these energies compare with each other is not known. Also, the viscoelastic dissipation in the bulk leads to a raise in the temperature of the material, which in turn a ects the viscoelastic properties of the material. is leads to a fully coupled thermomechanical problem.

is thesis is divided into three parts. e rst part begins by introducing the basic notions of Linear Elastic Fracture Mechanics, followed by a brief description of the experimental and theoretical studies on elastomer fracture. e next part of this thesis describes the experiments performed on the pure shear samples of polyurethane. Chapter describes the experiments performed in literature on polyurethane samples one one geometry. Chapter describes the additional experiments performed by the current author on samples of di erent geometries. Chapter presents additional results and observations from the experiments performed in the literature. e other part involves performing FE simulations to understand and predict the results obtained from the experiments. Chapter examines the hypotheses of hyperelastic and viscoelastic sti ening using the data obtained from experiments. Chapter describes a rate dependent cohesive zone to predict the crack speeds. Chapter describes the plane stress version of the Finite viscoelastic model and the analysis of the energy expenditure during the dynamic fracture.

L I T E R AT URE REVIEW

is chapter introduces the notion of Linear Elastic Fracture Mechanics and the asymptotic crack tip elds in the case of a stationary and a moving crack. Di erent loading modes and the notion of limiting speeds in di erent modes will be discussed. Di erent criterion to predict crack propagation will be described. e extension of LEFM to moderate to large deformations will be described as well. It ends with a brief discussion on the notion of Transonic cracks. . Early analysis of cracks can be a ributed to the studies of Inglis (Inglis, ), Gri th (Gri th,

) and Irwin (Irwin, ). Inglis studied the problem of an elliptic hole in a plate and determined the stress concentration factor arising from the elliptical corners of the plate.

e stress concentration factor was seen to depend inversely on the radius of curvature of the elliptical edge. As the aspect ratio of the ellipse is increased, the ellipse tends to a crack with the radius of the corners tending to zero. e stresses hence tend to become singular as the limit of a crack is reached. Westergaard used complex analysis to develop the variation of stresses in a plate under plane stress and plane strain assumptions. Irwin later used this method to develop the expression of stresses in the vicinity of the crack tip.

is lead to the de nition of Stress Intensity Factor as the quantity the describes the amplitude of the stress singularity at the crack tip. Gri th arrived at a criterion for the propagation of crack in a medium. He stated that the crack will propagate in a medium when the rate of increase in the surface energy because of crack propagation is o set by the rate of decrease in the potential energy of the body under the external loading as a consequence of crack growth. Irwin developed a relation between the Stress Intensity Factor (SIF) and the energy criterion of Gri th.

e study of propagation of cracks under small strains in a linear elastic material comes under the umbrella of Linear Elastic Fracture Mechanics (LEFM). In this section, a brief review will be made about the theory and some results of LEFM.

( )

e study of simple crack in a linear elastic material has been performed by a number of researchers. An exhaustive discussion on the theory can be found in (Anderson, ) and when the dynamic e ects are included in (Freund, ; Broberg,

). e theory of LEFM states that the state of stresses at the crack tip are singular and the near tip variation of displacement, stress and strain elds are universal. e e ect of external loading and geometry are introduced through a quantity called Stress Intensity Factor (SIF), K = lim r →0 σ √ 2πr . e hypothesis of LEFM are as follows.

• e material behaves linearly and is elastic. Homogeneity and isotropy are also assumed.

• e failure and inelastic processes occur in a very small zone near the tip. is is called e hypothesis of Small Scale Yielding.

• e analysis is restricted to simple cracks. Events such as branching are discounted.

For a simple crack in a material, three modes of crack propagation can be de ned based on the external loading. ey can be seen in the gure . . In mode-I, also called the opening mode, the crack faces open perpendicular to the crack path. e shear stresses are zero along the line of symmetry. Mode-II is characterized by the in-plane shear loading along the crack propagation direction. In mode-III, also called the tearing mode, the crack faces open perpendicular to the plane of crack propagation. Using the notion of SIF developed by Irwin, and the small scale yielding hypothesis, the stress state in the vicinity of crack tip can be expressed in terms of the SIFs in three modes as (Anderson, ) σ i j = K I √ 2πr Σ I i j (θ ) + K I I √ 2πr Σ I I i j (θ ) + K I I I √ 2πr Σ I I I i j (θ )δ j3 (no sum on j) ( . )

where K denotes the SIF in the corresponding mode as indicated by the subscript.

Σ i s determine the angular dependence of the stress distribution. Owing to the linearity of the material and small strains, principle of superposition has been used.

δ denotes the Kronecker delta and equals to when j = 3 or to otherwise. e hypothesis of small scale yielding implies that the fracture processes occur at regions that are very close to the tip. e region where the stress distribution is dominated by the singularity lies outside the the region where the fracture processes occur.

Dynamic case

When inertial e ects are included, the asymptotic crack tip elds have the same singularity as in the static case, but the angular dependence is now a function of crack speed as well. Hence, the equation . now becomes (Freund, ) Figure . : Polar coordinates centered at the crack tip.

σ i j = K I √ 2πr Σ I i j (θ, ) + K I I √ 2πr Σ I I i j (θ, ) + K I I I √ 2πr
Σ I I I i j (θ, )δ j3 (no sum on j) ( . )

At a distance r from the tip, the displacement eld varies as √ r and hence, the particle velocity is proportional to K I /E √ r . So, the kinetic energy density varies as

K ∼ ρ 2 2 K 2 I r E 2 , ( .
)
where ρ and E are the density and the Young's modulus, respectively. e strain energy density can be seen to vary with r as

U ∼ 1 2 K 2 I r E . ( . )
e ratio of these two energy densities can be seen to be

K U ∼ 2 E/ρ . ( . )
e ratio is independent of r and so, it can be seen that the inertial e ects may not be important in studying the problem as long as the crack speeds remain less than about 1 3 rd of the elastic wave speed E/ρ . e introduction of inertial e ects results in the notion of the limiting speed for the crack. It will be shown in the following section that the Energy release rate, G changes its sign as the crack speed passes the Rayleigh wave speed (c R ). Since a positive value of G is necessary for the crack to propagate, c R is the upper bound on the crack speed (Freund, ).

.

To predict the propagation of a crack in a body, an additional relation is needed along with the constitutive model that describes the behavior of the bulk material. e additional relation is in terms of a propagation criterion. As mentioned earlier, (Gri th, ) established a crack propagation criterion based on the energy balance of the system. Formally, the Gri th's criterion can be stated as follows. Denoting the total mechanical energy of the structure as Π , the Energy Release Rate, G can be de ned to be G = -∂Π ∂A ( . )

where A denotes the area of the crack surface. De ning the fracture energy γ to be the energy required to create a unit surface area of the crack, Gri th's criterion can be expressed as G <γ , Crack does not propagate.

G =γ , Crack may propagate.

In other words, the crack propagates in the body if the loss of potential energy of the bulk material equals to the increase in surface energy resulting from the crack growth. Gri th's criterion can also be seen to imply that the equilibrium crack length is the one that minimizes the total energy of the system, E Π + γ A.

In quasi static case, the relation between the Energy release rate and the SIF has been obtained in Irwin ( ) to be

G = K 2 I E ( . )
where E is the Young's modulus of the material. e energy release rate, G, can also be related to the J-integral proposed by (J R Rice,

) using the Eshelby stress tensor (Eshelby, ) as = e 1 .

∫ Γ (ρ 0 ψ Iσ .∇u) .n 0 dS, ( .

)
where ψ is the energy density, σ is the stress tensor and ∇u denotes the displacement gradient. For a crack in linear elastic material and under quasi-static loading conditions, it has been proved that = G. e path independence of the contour integral can be used to obtain the G from far eld quantities. e relation in the equation . can then be used to obtain the SIF from far eld measurements. e crack propagation condition expressed above in terms of the energy release rate can also be expressed in terms of the SIF. Using the relation between the G and SIF, the propagation condition can be expressed to be (in the case of a pure mode-I loading)

K I = K crit , ( . )
where K crit can be expressed in terms of γ .

. . Energy ux integral

e expression for the net energy ux through a contour Γ moving at a speed along with the tip has been arrived at by (Freund, ). By taking a dot product of the momentum equation with velocity and using the divergence theorem in space and time, the expression for the energy ux through Γ has been obtained as

F (Γ ) = ∫ Γ [(U + K) n 1 + v.σn] dS, ( . )
where

U = ∫ t -∞ σ i j
∂ 2 u i ∂t ∂x j dt is the stress work density. In the case of elastic material, stress work density can be seen to be equal to the strain energy density, since σ i j = ∂W ∂ϵ i j . It can be seen (in the later sections) that the ux integral in the above equation is contour independent as long as steady state conditions prevail regardless of the material behavior. e dynamic energy release rate can be de ned using the de nition of ux integral to be

G = lim Γ →0 F (Γ ) ( . )
e value of the energy ux integral is independent of the shape of the contour as the contour is shrunk to the tip. Hence, a contour that looks like in gure . has been used where to shrink to the tip by rst le ing δ 2 → 0 and then δ 1 → 0. Using the asymptotic crack tip elds, the energy release rate has been obtained to be (Freund, )

G = 1 -ν 2 E A I ( )K 2 I + A I I ( )K 2 I I + 1 2µ A I I I ( )K 2 I I I ( . )
where

A I = 2 α d (1 -ν )c 2 s D , A I I = 2 α s (1 -ν )c 2 s D , A I I I = 1 α s . ( . ) D = 4α d α s -(1 + α s ) 2 . ( .
)
where

α d = 1 - 2 c 2 d , α s = 1 - 2 c 2 s
and D is the Rayleigh function the root of which is the Rayleigh wave speed.

In nite medium

For the case of crack propagation in an in nite medium, the SIF for a dynamic crack can be expressed in terms of a static crack under the same loading conditions and a factor that depends on crack speed and material parameters. In an in nite medium or in a nite medium before the waves re ected o the boundaries reach the crack tip, the SIF can be expressed to be (page of (Freund, ))

K I (t, l, ) = k( )K I (t, l, 0) ( . )
where k( ) is a function of crack speed and material properties. e energy release rate can be expressed similarly in terms of its value at equilibrium as

G(t, l, ) = A I ( )k( ) 2 G(t, l, 0) ( . )
where A I ( ) is another function of crack speed. e function A I ( )k( ) 2 can be approximated by 1c R , where c R is the Rayleigh wave speed. Chapter . Introduction to fracture mechanics

Hence, for the case of a mode-I crack moving in an in nite medium, a crack tip equation of motion can be wri en based on the expression developed for the energy release rate.

1 -ν 2 E A I ( )K I (t, l, ) 2 = γ . ( . )
Since, in an in nite medium, the SIF can be expressed in terms of its value for a stationary crack under the same loading, the above expression can be simpli ed further to Eγ (1ν 2 )K I (t, l, 0) 2 = A I ( )k( ) 2 ≈ 1 -/c R .

( . )

For the case where γ is a constant, it can be seen that since K I (t, l, 0) is an increasing function of crack length, the le side of the above equation tends to zero as the crack becomes longer. Hence, the right side tends to zero as well in the limit → c R . Hence, the theoretical upper bound on the crack speed in mode-I is the Rayleigh wave speed (c R ).

. . Propagation direction

By expressing the external loading in terms of modes I, II and III, the condition for determining the crack path can be uncoupled from the condition for crack propagation in terms of surface energy. en, the criteria for the selection of propagation direction can then be expressed in terms of the SIF in the corresponding modes (Co erell & Rice, ). One of the rst conditions to determine the direction of crack propagation was established by (Goldstein & Salganik, ) as 'Principle of local symmetry'. It suggests that the path taken by a crack in bri le homogeneous isotropic material is the one where the local stress eld is of a mode-I type. In other words, the crack chooses the direction along which the SIF in mode-II (K I I ) is zero.

Other criteria such as the maximum hoop stress criterion (Erdogan & Sih, ), maximum energy release rate criterion (Hussain, Pu, & Underwood, ), stationary Sih energy density factor (Sih, ) also exist. ese criteria are known to be consistent with each other (Co erell & Rice, ) in that all these criteria state that if K I I 0, the crack extends with an abrupt, non-zero, change in the tangent direction to the path.

In static conditions, these criteria are applied for the case of an inclined crack subjected to a mixed-mode loading. Experiments were performed accordingly and the results of the experiments were compared with the theoretical predictions in (Sih, ). As a special case, when the crack is subjected to a pure mode-II loading, the branch angle was found to be 70deg, which is consistent with the experimental ndings. Some other criteria were established to predict the crack turning in numerical simulations. For instance, (Belytschko, Chen, Xu, & Zi,

) uses the loss of ellipticity criterion to determine the direction of propagation of a crack. is criterion has been mentioned to be similar to the loss of material stability (Belytschko, Liu, & Moran, ). Other criteria based on con gurational forces have been used in Finite Element simulations ( Özenc ¸, Kaliske, Lin, & Bhashyam, ). e con gurational force vector is used to determine the propagation direction of the crack. e crack propagation takes place between the elements where a double node is introduced when the magnitude of the con gurational forces satis es some criterion (like that of the energy release rate criterion). e mesh is then adjusted so that the element edge lies along the direction of the con gurational force.

Crack propagation can also be modeled by techniques such as as the phase-eld method (Bourdin, Francfort, & Marigo, ; C. Miehe, Welschinger, & Hofacker, ). In this method, the propagation of crack and the direction of propagation are taken care of by the evolution of the phase-eld itself. Similarly, in other methods such as the gradient damage models (Marigo, Maurini, & Pham, ), integral type non-local damage models (Peerlings, De Borst, Brekelmans, & Geers, ), ick Level Set method (Moës, Stolz, Bernard, & Chevaugeon, ), a separate criterion for crack propagation direction is not required.

. e asymptotic crack-tip elds and the other quantities like SIF reported in the previous sections are based on the linearity of the material and small deformations. To understand the fracture of polymer gels where the strains are no longer small and the material is no longer linear at those strain levels, Weakly Nonlinear eory of fracture (WNLT) has been introduced in (Bouchbinder, Livne, & Fineberg, ). Fracture experiments were performed on bri le polyacrylamide gels in (Livne, Bouchbinder, & Fineberg, ). LEFM was used to predict the crack opening pro le and strains ahead of the tip. It was, however, observed that the results from the experiments did not match the predictions by LEFM in regions close to the tip.

e discrepancy was a ributed to the non-linear behavior of the material at the tip since the material undergoes large deformations in that region. To address this discrepancy, higher-order displacement gradient contributions were included to develop the asymptotic crack tip elds while also considering the material nonlinearity at the observed strain levels. A neo-Hookean model was used for this purpose. It was observed that the inclusion of the higher-order terms resulted in more singular strain terms (r -1 ) than the traditional square root singularity observed in LEFM. e inclusion of nonlinear corrections was seen to result in Figure . : Comparison of experimental results with that of LEFM and WNLT at di erent crack speeds. e red circles indicate the experimentally measured values, black do ed line the predictions of LEFM, blue line the predictions of WNLT. e crack speeds are 0.20c s in (a), 0.53c s in (b), and 0.78c s in (c). (Bouchbinder, Livne, & Fineberg, ) strain elds that match the observations from the experiments (see gure . ). Figure of (Bouchbinder, Goldman, & Fineberg, ), a comparison has been made between the LEFM, WNLT and fully non-linear theory (Tarantino,

) has been made. It was observed that the fully nonlinear asymptotic elds match the crack pro le closely even near the tip.

. Elastomers can undergo large deformations before a crack can propagate through. Also, the behavior of elastomers at such large strains cannot be described by the linear elastic material model and through small strain measures. Hence, the singularity at the crack tip may be di erent from what has been described in the earlier sections. Some studies establish the form of stresses at the crack tip under large strains using non-linear strain measures. Under quasi-static conditions, assuming the strain energy density to depend just on the rst invariant of Green strain and the material to be incompressible, (Knowles, ) establishes the asymptotic crack tip elds under mode-III loading conditions. (Geubelle & Knauss, ; Krishnan, Hui, & Long,

) obtained the asymptotic crack tip stress for a crack in a Generalized Neo-Hookean model for a crack in mixed mode loading. e model condenses to a Neo-Hookean model for n = 1. For that case, the singular crack tip Piola stresses are seen to vary as P 21 ∼ R -1 2 sin(Θ/2) ( . )

P 22 ∼ R -1 2 cos(Θ/2) ( . )
where R and Θ are measured from the crack tip to the particle in the undeformed con guration. Similarly, singular Cauchy stresses are seen to vary as

σ 21 ∼ r -1 2 -sign(θ )(2 -3/2 ) √ 2r sin 2 θ + cos θ 1 - cos θ 2r sin 2 θ + cos θ 1/2 ( . ) σ 22 ∼ r -1 1 2r sin 2 θ + cos θ ( . )
where r and θ are measured from the crack tip to the particle in the deformed con guration (see gure .

). e di erence in the singularity of the stresses can be noticed. Also, the singularity of the Cauchy stresses varies with the angle. For instance, σ 22 is singular as 1/r along θ = 0 while it is singular as 1/r 2 among θ = π /2. is is di erent from LEFM where there is no di erence between the deformed and the undeformed con gurations and as a result all the stress measures coincide.

e singularity is hence same in both the con gurations. e asymptotic crack tip elds in dynamic case was obtained by (Tarantino, ). Using a similar variable separable representation to that in (Geubelle & Knauss, ; Knowles, ), the leading singular terms were obtained. e deformation maps were observed to vary as

φ 1 ∼B 1 (t; ˜ )w(Θ)R + o(R) ( . ) φ 2 ∼A 2 (t; ˜ )ˆ (Θ; ˜ )R 1/2 + o(R) ( . )
where ˜ is a function of crack speed and material parameters (µ). w determines the angular variation of φ 1 . e function ˆ can be seen below.

ˆ (Θ; ˜ ) = sgn(Θ) (1 -˜ 2 sin 2 Θ) 1/2 -cos Θ 2 1/2 .
e PK stress components were seen to vary as

P 21 =2µR -1/2 A 2 1 2 ˆ cos Θ -ˆ ,Θ sin Θ + o(1) ( .
)

P 22 =2µR -1/2 A 2 1 2 ˆ sin Θ + ˆ ,Θ cos Θ + o(1) ( . )
It shall be noted that the order of singularity is similar to that of statics (R 1/2 ). e stress components that are singular are same as that in the static case.

.

As mentioned earlier, in LEFM, the crack speeds are restricted to c R from below in mode-I and c S in mode-III. e Energy release rate can be observed to change its sign to negative once the crack speed exceeds c R and c S in mode-I and mode-III, respectively. Since a positive energy supply is needed for the crack to propagate, the crack speeds are restricted to below the corresponding limiting speeds. e sign of energy release rate is positive even in mode-II when the crack speeds are smaller than c R . is translates to the condition that the stresses near the tip are singular as 1/ √ r . For crack speeds between c R and c S , the sign of G becomes negative as earlier, and hence speeds in this regime are forbidden. However, for crack speeds that exceed c S , it was observed that stresses at the tip are singular as 1/r m , where m ≤ 0.5 with the equality holding at a crack speed of √ 2c S . e shear stress and horizontal velocity component for such a case can be seen to be (Freund, )

σ x = K * 2 (t) (2π ) 1/2 cos(mθ l ) r m l + H (-ξ -β s |η|) sin(πm) (-ξ -β s |η|) m tan m , ( .
)

u x = K * 2 (t) 2µα l (2π ) 1/2 sin(mθ l ) r m l -2α l β s H (-ξ -β s |η|) sin(πm) (-ξ -β s |η|) m tan m , ( . )
where H denotes the Heaviside function, ξ , η, β, m are functions of crack speed and material parameters.

πm = arctan[4α l β s /(1 + α 2 s ) 2 ].
As mentioned earlier, ).

m = 1/2 when = √ 2c s . It can be seen that the stresses and particle speeds exhibit jumps along the surfaces -ξβ s |η| = 0 and are singular along the surface as well. Since the singularity of the stress is m(≤ 1/2), the energy release rate is 0 for all crack speeds above c S except at √ 2c S . is was seen to be the consequence of the sharp crack assumption and assuming that the failure happens over a nite region instead of just at the tip was seen to overcome this limitation. Overall, mode-II cracks are permi ed to travel at speeds that exceed c S . Such cracks have been observed experimentally, for instance, in (Rosakis et al., ; Gori, Rubino, Rosakis, & Lapusta, ). In such experiments, shock waves that travel along with the crack were observed (see gure . ) in the photoelastic fringes. A sharp rise in particle velocities a er the shock passes through were observed (Mello, Bhat, Rosakis, & Kanamori, ). ese cracks were modeled using Finite Elements in (Miller, Freund, & Needleman, ). Simulations performed using molecular dynamics such as in (Abraham et al.,

) also con rm the existence of such cracks. Transonic cracks have been observed in opening mode in fracture experiments on rubber. Such cracks have been rst reported to be observed by (Petersan, Deegan, Marder, & Swinney, ). More experiments performed by (Chen, Zhang, Niemczura, Ravi-Chandar, & Marder, ; Mai, Okuno, Tsunoda, & Urayama, ; Corre, Coret, Verron, Leblé, & Le Lay, ) con rm this. More details on this will be discussed later in this thesis.

In this chapter, the theory of LEFM was brie y described. e asymptotic crack tip elds were presented for static and dynamic cases. Di erent propagation conditions for crack were brie y presented. e crack tip elds for the case of large deformations and Neo-Hookean model have been presented in static and dynamic cases. A brief discussion about Transonic cracks has been given. is chapter introduces the available studies on the dynamic fracture. Di erent experimental specimen geometries used for elastomer fracture will be presented. Some available experimental and theoretical studies regarding the quasi-static and high speed fracture of elastomers will be presented. A brief summary of the notion of energy release rate and the energy ux integral will be presented. ; Ravi-Chandar and Knauss, a and the references cited therein. e fracture of Homalite-was studied in Ravi-Chandar and Knauss, a. e propagation of crack was studied under applied loading and the dimensions of the specimen were chosen so that considerable crack propagation occurred before the stress waves re ected from the boundaries interacted with the crack. e loading was applied on the crack faces through an electromagnetic device. e crack propagation was studied under di erent applied loads and loading rates.

One of the observations of that study is the change in the appearance of the fracture surfaces at di erent SIFs. It was observed that the crack speed remained almost the same during propagation but the size of the damage zone varied as the SIF varied. A change in the appearance of the fracture surfaces was noted in accordance with the change in the SIF (from photo-elastic measurements). e surface pa ern was reported to change from 'mirror' to 'mist' to 'hackle', where the roughness progressively increased. Crack branching was observed to occur at crack speeds of about % of c R . It was observed that the interaction with re ected waves is not necessary for the crack to undergo branching. A mechanism for branching has been proposed in Ravi-Chandar and Knauss, a by taking into account the interaction of micro branches and voids in the crack path. Branching was thought to occur when the small cracks deviate from the main crack plane. e results of the study were compared with that of Dally, (see gure . ). It can be seen that while Dally, reported a unique relationship between the SIF and crack speed, in the studies of Ravi-Chandar and Knauss, b, the crack speed remained constant in the experiment while the SIF was seen to vary (with a corresponding change in crack surface appearance).

Further experiments were performed on PMMA by Fineberg, Gross, Marder, and Swinney, . Some di erences can be noted between the experimental setups in Fineberg et al.,

and Ravi-Chandar and Knauss, a. In Ravi-Chandar and Knauss, a, a gap was machined into the material and the loading was directly applied onto the crack faces using an electromagnetic device. In Fineberg et al.,

, the specimen was rst stretched quasi-statically to a certain strain level and a crack is then introduced in the specimen. Also, in Ravi-Chandar and Knauss, a, the crack tip SIF is measured by using the method of caustics and the crack speed by a high speed camera. In Fineberg et al.,

, the crack speed is measured by a change in the electrical resistance of the thin layer of Aluminum deposited onto the specimen faces. Starting from about 0.34c R , the crack speeds were observed to oscillate and this was seen to coincide with the development of roughness on the crack surfaces. e authors term this as 'microbranching instability'.

Many other fracture experiments were performed in the context of dynamic fracture, for instance, to investigate the e ect of the fracture process zone on the terminal crack speed. In the references mentioned above, the crack was observed to undergo branching at a speed of about % to % of c R . is was a ributed to the development and propagation of micro branches from the main crack tip. To investigate if the cracks can propagate at higher speeds if the micro branches are suppressed, experiments were conducted in Washabaugh and Knauss, . Micro branching is suppressed using two methods. In the rst, the specimen is cut in half and the cut halves are machined to make the surfaces smooth. ey are then molded in a vacuum at high temperature till the interface is strong enough, but not as strong as the base material. In the other method, an array of holes were drilled into the material along the prospective crack path at various spacing. It was reported that using the above techniques, the crack speed was able to approach c R in the cases when the interface is 'weak enough' (see gure . ). A detailed review of the dynamic fracture of nominally bri le materials can be found in Ravi-Chandar and Yang, . e results of various models to describe the evolution of defects and propagation of crack have been discussed.

. ' '

Even though the upper bound on the crack speed is c R , in reality, cracks hardly reach that speed unless special measures are taken as discussed earlier. At crack speeds of about % to % of c R (Ravi-Chandar & Knauss, b), the tip is known to undergo branching where the main crack branches into multiple branches which later propagate. It was reported that the crack speeds remained at about % of c R even when the loads applied were higher and the crack did not branch as well. e authors a ribute the lower terminal speed observed and the events of branching to the growth and interaction of microcracks along with the main crack.

Microbranching instability

In the experiments performed on PMMA by Sharon, Gross, and Fineberg, ; Sharon and Fineberg, , it was observed that starting from a certain crack speed (called c ), the crack surfaces became rougher. Oscillations were also observed in the crack speed once it exceeded c (≈ 0.4c R ). e authors termed the phenomenon as 'Microbranching instability' (Fineberg & Marder, ; Fineberg et al., ). It was shown in those works that there was a corresponding increase in the area ). Oscillations in crack speed a er a certain limit can be seen.

of the fracture surface created once the crack speed exceeded c . e increase in fracture energy during fracture was regarded as a consequence of an increase in the mean surface area of the crack. e same behavior was reported for glass and polyacrylamide gels as well in (Bouchbinder et al.,

). It shall be noted that the cited references do not make any special comments about 'macro-branching' observed during the experiments as opposed to Ravi-Chandar and Knauss, b.

Oscillatory instability

It was also reported that the c mentioned earlier depended on the crack acceleration as well (page of Bouchbinder et al.,

). c was seen to be an increasing function of crack acceleration, increasing from its value (≈ 0.4c R ) at lowest accelerations. It was also seen that c can be increased by suppressing the 'instabilities', for instance, by decreasing the specimen thickness. To this e ect, fracture tests were performed on polyacrylamide gels in Livne, Ben-David, and Fineberg, . Specimen were stretched to 'very large strains' and a crack was introduced. Also, a thin specimen was used to suppress any instabilities. Using both these techniques, it was reported that the cracks traveling as fast as 0.9c s were achieved. It was observed, however, that instead of the 'micro-branching' instability mentioned above, an oscillatory the crack path was observed. Similar oscillatory cracks were also observed in the works of Deegan, Petersan, Marder, and Swinney, , although in a di erent context. . e investigation of fracture of elastomers can be traced back to Greensmith and omas,

. In that study, the fracture behavior and properties of elastomers were studied by using rubber specimen in 'tearing' and 'pure shear' con gurations. A brief description of the specimen con gurations and the fracture energy in those con gurations can be seen below.

. . Tearing

A 'trouser' specimen is used for testing in this con guration. e specimen contains two 'legs' which are held on to a tensile machine that are pulled apart. e fracture energy in this con guration is computed using the stretch level in each arm and the expression . . e force measurement can be obtained from the tensile machine and the stretch level can be obtained by measuring the nal lengths of arms. e strain energy density can then be calculated.

G = wh 0 - 2 e 0 λF . ( . )
e values of h 0 and e 0 can be obtained from the geometry as shown in the gure . . w is the strain energy density corresponding to the strain state at the loading condition.

. . SENT

Other specimen con guration that can be used is the 'Single Edge Notched Tension'. e geometry can be seen in the gure . . In this case, the energy release rate can be computed as

G = 2βwl, ( . )
where l denotes the crack length and β depends on the stretch level and can be approximated as 3λ -1/2 .

. . Pure Shear

In this case, a specimen that is wide enough (about times as tall as the specimen height) is clamped along its longer edges onto a tensile machine. e specimen is stretched to the target level and a crack is introduced. Once the crack reaches the center of the specimen, the energy release rate can be calculated by

G = wh 0 , ( . )
where h 0 is the initial specimen height (see gure . ).

.

Early experimental studies on the fracture of elastomers can be traced back to the works of R. S. Rivlin and omas, . In that study, specimens of natural rubber were torn in the con gurations described earlier to obtain a relation between the fracture energy vs speed. Further experiments were carried out by Greensmith and omas, ; Greensmith, and the references cited therein. Most of the works were carried out in the tearing con guration. Based on the nature of the traction separation curves observed and the fracture surfaces, the overall fracture behavior was classi ed into 'slip-stick' and 'steady' tearing. In 'slip-stick' tearing, large oscillations were observed in the crack speeds as well as the traction curves. It shall be observed that in all the studies mentioned above, a Gri th-like criterion has been used to interpret the experimental data in terms of fracture energy vs the crack speed.

On the other hand, the fracture of polymers was studied by Knauss and coworkers in W. Knauss and Mueller, . In those studies, the pure shear con guration was used. e material tested is Solithane . e results of the experiments were presented in terms of the macroscopic strain imposed vs the crack speed observed (see gure . ). ).

To interpret the data, an analytical investigation was carried out. Using a cohesive zone type model to describe the behavior of the material failing at the tip, a relation between the crack speed, the applied load, and the viscoelastic properties of the bulk material was obtained in W. Knauss, ; W. Knauss and Mueller, ; W. Knauss, . Two failure criteria were considered for crack propagationfailure strain criterion and an energy criterion.

In failure strain criterion, the cohesive zone was assumed to fail completely when the strain value reaches a certain limit. e tractions in the cohesive zone were taken to completely drop to zero at this value. is criterion can be wri en as u (0, 0, a) = u 0 .

( . )

In the above equation, the coordinates (0, 0) denote the physical crack tip, and a denotes the crack speed. u 0 is the strain value at which the failure occurs. In the energy criterion, a cohesive element was assumed to fail when the energy expended on the element by the surrounding material reached a certain value. is can be wri en as

- ∫ A c T n (X ). ∂u ∂X (X , 0, a) dX = γ , ( .
)
where A c denotes the area spanned by the failing material, T n are the tractions, u is the y-component of the displacement of the crack surface. γ denotes the fracture energy. e opening displacement of the crack faces under external loading for viscoelastic material was obtained by using the elastic solution under the same loading and correspondence principal. e constitutive behavior of the failing material is speci ed as a relation between the traction and separation. In W. Knauss, , a constant traction value has been assumed over a span of α. Using both the traction and separation, equations . and . can be evaluated. e displacement criterion hence becomes

3K √ α √ 2π E ∞ ϑ α a = u 0 , ( .
)
where K is the SIF corresponding to the applied load, E ∞ is the long term elastic modulus and ϑ is a parameter that depends on the viscoelastic properties of the material. e energy criterion becomes

3 4 K 2 Θ R α a = E ∞ γ , ( . )
where Θ R is a function that depends on the viscoelastic properties of the material. e predictions of the above criterion were compared with the experiments whereby the values of parameters like γ and α were obtained. It shall be noted that a similar development was carried out by Schapery, b, a, c. e above-mentioned studies use a failure zone at the tip to describe the fracture processes. Some other studies exist that follow Gri th's failure criterion that the energy to failure comes from the reduction of the energy in the bulk Christensen, . In such a case, the global conservation of energy gives

aγ + ∫ V Λ d + U = 0. ( . )
In the above equation, V is the volume of the material, Λ is the energy dissipation rate in the material, and U is the rate of change in strain energy. Inertial e ects have been neglected. An earlier study in this direction was done by Kostrov and Nikitin, . However, that study leads to a conclusion that the sharp crack modeling approach or Gri th's failure criterion leads to a conclusion that the viscoelastic properties do not a ect the crack propagation (James R. Rice,

). is can be seen by examining equations . or . in the limit of α → 0. In that case, the equation . , for instance, reduces to

3 4 K 2 Θ R (0) = E ∞ γ . ( . )
It can be seen that only the initial (glassy) viscoelastic properties enter the equation as opposed to the case when α is not zero. However, Christensen, a ributes this conclusion to the omission of heat source term which indirectly leads to elimination of dissipation in the formulation of Kostrov and Nikitin, . e comparison between the theory and experiments has only been done for a limited range in Christensen, . Other studies where the sharp crack modeling technique has been used are that of B. N. Persson and Brener, ; B. N. J. Persson, Albohr, Heinrich, and Ueba, . In this approach, the point of departure is the conservation of energy as earlier. In this case, however, the singularity at the crack tip has been considered explicitly as opposed to Christensen,

, where multiple terms of a power series were considered to approximate the crack tip singularity. e point of departure of B. N. Persson and Brener, is

G a = G 0 a + P, ( . )
where the rst term on the le indicates the global energy expenditure, the rst term on the right indicates the energy consumed by the crack tip during the fracture process and P denotes the energy dissipation in the bulk material. e dissipation in the bulk material was calculated explicitly by using a power-law model for the relaxation function. e formulation was compared with the experimental results from Gent, b. e slope of the energy release rate vs the crack speed curve was found to be . , which was found to be in good agreement with the theory.

More experiments were conducted by Lake, omas, and Lawrence, ; Lake, Lawrence, and omas, ; Tsunoda, Bus eld, Davies, and omas, on the fracture of elastomers. In Lake et al.,

, the phenomenon of cavitation of elastomers has been investigated. More discussions on this will follow in later chapters. Other works such as that of Lake et al.,

; Tsunoda et al., investigate the effect of various parameters such as the specimen geometry (thickness in particular), the chemical composition, the e ect of llers, etc. on the fracture properties of the specimen. e crack speeds in Tsunoda et al., are considerably smaller than that in Lake et al.,

. Hence, inertial e ects were not considered in that study.

High speed fracture of elastomers

Studies in which the crack speeds are much smaller than the wave speed in the material do not need the consideration of inertial e ects. When the crack speed is higher, typically more than about % of the shear wave speed, inertial e ects become important. Such studies were performed by Stevenson and omas, ;

Kadir and omas,

; Gent and Marteny, a. Since the crack speeds in the studies are higher, a comparison to wave speeds was made. In Stevenson and omas, , fracture tests were performed by using balloons made of natural rubber. A needle was used to initiate the crack once a target strain level was reached. e wave speeds were also measured by performing retraction experiments on rubber sheets. However, the crack speeds and wave speeds reported in that study are on the deformed con guration. It was observed and remarked that the crack speeds obtained in that study are considerably higher than the wave speeds measured and even the 'strain dependence of the wave speed' does not account for the observed higher crack speeds. e high rate of extensions at the crack tip was thought to modify the material sti ness and hence the wave speeds.

is was regarded as a cause of the observed higher crack speeds. No speci c mention of 'Transonic cracks' was made. e crack surface was found to have wavelike series of undulations, -µm apart and -µm deep over most of the fracture surface. Bi-axial strain state was stated to be the probable reason for the undulations. Livne et al., report such oscillatory crack path, but does not mention the study of Stevenson and omas, . In Gent and Marteny, a, more experiments were made on un lled and carbon-black-lled natural rubber specimens in pure shear geometry. Wave speeds were measured by performing retraction tests in a di erent study (Gent & Marteny, b). e crack speeds are compared to 0.3 s (Mo , ), where s is the velocity of stress pulse in the material. e crack speeds and stress pulse speeds were reported in the deformed con guration as well. Also, specimens of di erent geometry (speci cally heights) were subjected to the same stretch level and the crack speeds were measured. It was observed that the crack speeds become independent of the specimen height starting from a certain height. No mention was made about Transonic cracks as well. . ese experiments were perhaps the rst to report the existence of 'Transonic cracks' in the opening mode in elastomers. It was reported that in natural rubbers, starting from a certain stretch level, the crack speeds exceeded the Shear wave speed. It was also observed that the crack speeds became independent of the specimen height in this regime. More on this will be discussed later in this thesis.

. e Energy release rate is related to the crack speeds observed in the experiments throughout the literature. It is measured by using the strain energy density in the material far ahead of the tip and far behind it -that the energy lost by the body as the crack propagates through it is consumed by the fracture processes.

e expression for the energy release rate can also be obtained though contour integrals established in the literature.

is section brie y reviews two of such integrals to bring out their similarities and di erences.

. . Energy ux integral

In Freund, , an expression was obtained to compute the instantaneous energy ux through a contour that moves along with the crack with a speed along x-axis.

e expression can be seen to be

F (Γ ) = ∫ Γ σ ji n j ∂u i ∂t + (U + T ) n 1 dS, ( .
)
where U is the stress work density and T is the kinetic energy density. σ i j denotes the ij-th component of Cauchy stress in the Cartesian basis and n j is the j-th component of the outward normal to the contour Γ . e expressions for U and T are given by ) where the summation on repeated indices is implied.

U = ∫ t -∞ σ ji ∂ 2 u i ∂t ∂x j dt and T = 1 2 ρ ∂u i ∂t ∂u i ∂t , ( .
Path (in)dependence e path dependence of the above contour integral was examined by considering the di erence in the value of ux integral evaluated along two di erent contours Γ 1 and Γ 2 . the di erence can be seen to be

F (Γ 2 ) -F (Γ 1 ) = ∫ A 12 ∂U ∂t + ∂U ∂x 1 + ρ ∂u i ∂t ∂ 2 u i ∂t 2 + ∂ 2 u i ∂t ∂x 1 dA, ( . )
where Γ 2 is outside Γ 1 and A 12 is the area enclosed between the two contours.

In the special case of a steady crack growth, the elds remain invariant in a frame traveling with the crack tip. In such a case, the elds depends on x 1 and t only through the combination, ζ = x 1t. us, the terms in the parenthesis of the equation . go to zero and the contour integral is hence independent of Γ chosen. Also, in that case, using

∂ ∂t = - ∂ ∂x 1
, the contour integral reduces to

F (Γ ) = ∫ Γ (U + T )n 1 -σ ji n j ∂u i ∂x 1 dS. ( . )
U and T reduce to

U = ∫ ∞ x 1 -t σ ji ∂ 2 u i ∂x j ∂ζ dξ and T = 1 2 ρ 2 ∂u i ∂x 1 ∂u i ∂x 1 . ( . )
It shall be noted that the path independence for a steady state crack propagation has been established without using a constitutive model and hence is valid even for plastic, viscoplastic and viscoelastic materials.

Energy release rate e energy release rate has been de ned using the expression in equation . as (Freund, )

G = lim Γ →tip F (Γ ) ( . ) = lim Γ →tip 1 ∫ Γ σ ji n j ∂u i ∂t + (U + T ) n 1 dS ( . )
Under steady state, the above equation reduces to

G = ∮ tip (U + T )n 1 -σ ji n j ∂u i ∂x 1 dS. ( . )
. .

Eshelby's integral

In addition to the energy ux integral described above, other contour integrals exist in fracture mechanics that have been developed from the perspective of the mechanics in material manifold (Gurtin, ; Steinmann, ). e development in this direction relies on the change in the energy of the system as a result of the change in the material con guration (denoted B X (t)) of the system (Runesson, Larsson, & Steinmann,

). e notions of a time invariant reference con guration (denoted B ξ ) and a time varying spatial con guration (denoted B x (t)) have been de ned for this purpose. e three con gurations are related by maps as de ned.

e maps φ :

(B ξ , R + ) → B X (t) such that (ξ , t) φ -→ X φ(ξ , t)
between the reference con guration and the material con guration, φ :

(B ξ , R + ) → B x (t) such that (ξ , t) φ -→ x φ(ξ , t)
between the reference con guration and the spatial con guration and φ :

(B X (t), R + ) → B x (t) such that (X , t) φ -→ x
φ(X , t) between the material con guration and the spatial con guration are de ned. e three maps can be seen to be related by the composition φ(ξ , t) = (φ • φ)(ξ , t). As mentioned earlier, the intention is to study the dissipation in the system as a consequence of the change in material con guration, that is a variation of the map φ. Assuming the material to be homogeneous, free of body forces and elastic, the dissipation can be seen, under isothermal conditions to be

D(w 0 ) = ∫ ∂Ω X C k n 0 .w 0 dS, ( . )
where w 0 ∂X tip ∂t is, in this case, the crack speed in the material con guration, n 0 is the outward normal to the contour ∂Ω X and C k is the Eshelby stress tensor in dynamics, de ned as

C k (ψ + k)I -F T P . ( . )
Here, ψ is the strain energy density per unit volume in material con guration, k is the kinetic energy density per unit volume in material con guration, F ∂x ∂X is the deformation gradient and P is the rst Poila Kirchho stress tensor. e contour ∂Ω X is taken as shown in the gure . . From the gure, it can be seen that ∂Ω X = ∂Ω X,ext ∪ Γ ± c ∪ Γ tip . Γ ± c are the top and bo om of the portion of crack face that are not enclosed within Γ tip . Γ tip is a contour that goes around the tip from the bo om face of the crack to the top face. Realizing that w 0 = o on ∂Ω X,ext and Γ ± c , the expression for dissipation reduces to just over the contour enclosing the tip, shrunk to the tip.

D(w 0 ) = ∮ tip C k n 0 .w 0 dS. ( . )
is allows to de ne the 'con gurational forces' acting at the tip that is conjugate to the crack speed. e dissipation in the above equation can hence be expressed as

D(w 0 ) = w 0 . f tip , ( . )
where

f tip ∮ tip C k n 0 dS, ( . )
is the con gurational force acting at the tip.

If the material is inelastic, the dissipation in the bulk also contributes along with the dissipation at the tip. In that case, the dissipation expression becomes ) where D bulk = ∫ Ω X \Ω tip ∂ψ ∂k k dV . k indicates the array of internal variables and , a suitable contraction operation. Ω tip is the volume of material enclosed within Γ tip .

D(w 0 ) = ∮ tip C k n 0 .w 0 dS + D bulk = w 0 . f tip + D bulk , ( .

Steady state

For the case of a straight crack propagating in the X 1 -direction under steady state conditions, as stated previously, the elds depend on X 1 and t through the combination ζ = X 1w 0 t. Hence, the spatial and temporal derivatives are related as

∂ ∂t = -w 0 ∂ ∂X 1
. Hence, the bulk dissipation in the equation . can be wri en as

D bulk = -w 0 ∫ Ω X \Ω tip ∂ψ ∂k ∂k ∂X 1 dV ( . )
e total dissipation becomes

D(w 0 ) = w 0 e 1 . ∮ tip C k n 0 dS - ∫ Ω X \Ω tip ∂ψ ∂k ∂k ∂X 1 dV = w 0 e 1 . ∫ ∂Ω X C k n 0 S dS (
. ) A similar expression can be found in Simha, Fischer, Kolednik, and Chen, ; Simha, Fischer, Shan, Chen, and Kolednik, in the context of plasticity.

Path (in)dependence

Considering another contour Γ 2 that encloses Γ tip , a relation can be found between the integrals on the two contours during the steady state.

w 0 e 1 . ∫

Γ 2 C k n 0 dS - ∮ tip C k n 0 dS = w 0 e 1 . ∫ Γ tip ∪Γ 2 ∪Γ ± c C k n 0 dS = w 0 e 1 .
∫ since on crack faces e 1 .n 0 = 0. e divergence can be evaluated as below.

Ω 2 \Ω tip div 0 C k dV , ( . )
(C k ) i j, j = (ψ + k)δ i j -F ki P k j , j = ∂ψ ∂k ∂k ∂X i . ( . )
Steady state has been assumed. Hence, e 1 . div 0 C k = ∂ψ ∂k ∂k ∂X 1 .

Hence, the di erence in contour integrals can be seen to be

w 0 e 1 . ∮ tip C k n 0 dS - ∫ Γ 2 C k n 0 dS = -w 0 e 1 . ∫ Ω 2 \Ω tip ∂ψ ∂k ∂k ∂X 1 dV . ( . )
e RHS term can be identi ed as the dissipation in the volume between the two contours (see equation . ). Hence

w 0 e 1 . ∮ tip C k n 0 dS - ∫ Γ 2 C k n 0 dS = D Ω 2 \Ω tip bulk . ( . )
Hence, unless the material is elastic, which would mean zero bulk dissipation, the contour integral is not path independent.

. .

Small strains

Under small strain conditions, the integral in . can be simpli ed. Under small strain conditions, F ≈ I + ∇u, P ≈ σ .

∮ tip C k n 0 dS ≈ ∮ tip (ψ + k)I -(∇u) T σ n 0 dS - ∮ tip σn 0 dS. ( . )
If, at the tip, the stress components are singular as r p where p < 1, the second term on the right side goes to o (Gurtin, ). Hence, the contour integral can be expressed as

∮ tip C k n 0 dS = ∮ tip (ψ + k)I -(∇u) T σ n 0 dS, ( . )
and the component along e 1 can be found to be

e 1 . ∮ tip C k n 0 dS = ∮ tip (ψ + k)n 1 -σ i j n j ∂u i ∂X 1 dS. ( . )
. .

Comparison between the two integrals

At this point, a comparison can be made between the two integrals in equations . and . . e former integral contains

U = ∫ t -∞ σ ji ∂ 2 u i ∂t ∂x j
dt , which is the stress work density, while the le er integral contains ψ , which is the strain energy density. In the elastic case, the two integrals, hence, coincide. Hence, the integral in . can be seen as the energy ux to the tip and the con gurational force at the tip as well.

Also, the integral in . can be seen to be path dependent except in the case of an elastic material in steady state. e integral in . , however, can be seen to be path independent in the steady state scenario regardless of the constitutive behavior of the material. e la er also has the signi cance of representing the energy ux through the contour Γ , while the former doesn't, even in the steady state (except if the material is elastic). For instance, in the case of crack propagation under steady state conditions in a strip (see gure . ), the rst integral (in equation . ) when evaluated on a contour far from the tip results in

e 1 . ∫ Γ ∞ C k n 0 dS = [ψ (+) -ψ (-)] h, ( . )
which is o en stated as the energy release rate in literature (Greensmith & omas, ). e integral in . on the same contour gives a di erent value since it involves evaluation of U along the part of the curve labeled (-). Along this curve, U (-) ψ (-), since the evaluation of U requires integration over the entire history of the strain states. But along the part of the curve (+), U (+) = ψ (+). Hence, as stated earlier, unless the material is elastic,

∫ Γ ∞ (C k ) 1j (n 0 ) j dS ∫ Γ ∞ (U + T )n 1 -σ ji n j ∂u i ∂x 1 dS ( . )
An illustration of the scenario where the inelastic e ects can a ect the evaluation of the contour integral can be seen in Freund and Hutchinson, . In that study, the propagation of a crack in rate-dependent plastic solid was studied.

e plastic deformation was assumed to remain within the region dominated by the SIF. Assuming a steady state, the relation between the near eld and far eld Energy release rates de ned according to the equation . was obtained to be

G tip = G - ∫ h -h U * dx 2 , ( .
)
where h is the height of the plastic wake le behind by the propagating crack. In the above equation, G tip has been de ned as

G tip = f (m) 1 -ν 2 E K 2 tip . ( . )
K tip denotes the amplitude of the near-tip stress singularities. Because of the constitutive assumption of the material behavior, the elastic strain rates dominate the plastic strain rates and hence, the solution near the tip is similar to that in an elastic material. G is de ned as the far eld energy release rate de ned as

f (m) 1-ν 2 E K 2 .
Note the di erence in the SIF used (K for far eld vs K tip for near eld). U * (x 2 ) = lim x 1 →-∞ U (x 1 , x 2 ) is the stress work density locked in the wake behind the tip. is chapter presents an overview of the experiments of dynamic fracture performed on materials like Homalite-, PMMA, and glass.

is is followed by a description of the experimental setups used to study polymer fracture. Some theories on polymer fracture in the literature were presented. It was followed by an overview of the experimental studies available on high speed fracture of elastomers. It ends with a discussion about the integral formulations of energy release rate in the context of viscoelastic fracture.

is chapter provides a brief introduction to the polymer materials, dependence of their behavior on the strain rate. A brief description of the hyperelastic and viscoelastic models has been given as well.

. ). e chains are linked to each other through cross-linking. Elastomers are o en reinforced by llers such as Carbon black or Silica. is increases the sti ness of elastomers. Elastomers are characterized by their ability to undergo large deformations elastically.

When an elastomer is subjected to uniaxial tension the response typically looks like in the gure . . e slope of the strtess strain response can be seen to initially decrease from its value at zero strain level. A er a certain strain level, however, the slope of the curve starts to raise again to higher values. As a result of sliding between the molecular chains, the behavior of the elastomers are a function of temperature and the loading rate. Based on the temperature and loading rate, the behavior of the elastomers can be classi ed into three categories.

• Bri le and rigid, which is called a glassy state.

• So and highly deformable, which is called the rubbery state.

• Transition behavior, where the polymer chains are highly mobile which makes the behavior rate dependent.

When subjected to high loading rates and/or smaller temperatures, the material has a modulus E 0 , which is called glassy modulus. e modulus at small loading rates and higher temperatures is called rubbery modulus, E ∞ .

. A nondefective material body is a simply connected region B of a three-dimensional Euclidean manifold M, called the material manifold. e elements of this manifold are so-called material points X (Marsden & Hughes, ; Gérard A Maugin, ). A con guration of B is a mapping ϕ : B → S ⊂ R 3 . e set of all con gurations of B is denoted C . A motion of a body is a curve in C , that is, a mapping t → ϕ t ∈ C .

Le ing ϕ : B → S to be smooth, the deformation gradient is de ned as the tangent map of ϕ and is denoted F = T ϕ (Yavari,

). At each point X ∈ B, it is a linear map, F (X ) : T X B → T ϕ(X ) S, where T X B and T ϕ(X ) S denote the tangent spaces to B and S at X and ϕ(X ), respectively.

Assuming {X A } and {x a } to be the local charts on B and S, respectively, the components of F are

F a A (X ) = ∂ϕ a (X ) ∂X A . ( . )
e transpose of F is de ned by

F T : T x S → T X B, Fv 0 , v = v 0 , F T v G ∀v 0 ∈ T X B, v ∈ T x S, ( . )
for X ∈ B, x ∈ S. It has been assumed that B and S are Riemannian manifolds with inner products , G and , , respectively. e right Cauchy-Green deformation tensor is de ned by

C(X ) : T X B → T X B, C(X ) = F (X ) T F . ( . )
In component form,

C A B = (F T ) A a F a B .
Isotropic hyperelasticity e theory of hyperelasticity assumes the existence of a strain energy density functional, ψ , that depends on the deformation gradient ψ (F ). e nominal stress,

P : T * X B → T x S, N → T = P N . ( . )
is related to the strain energy density function as

P = ∂ψ (F ) ∂F . ( . )
Requirements of objectivity require the strain energy density to depend on F through the combination F T F = C. e Piola-Kirchho stress tensor (S) is de ned as

S : T * X B → T X B, N → T = SN . ( . )
such that P = FS. It is related to the strain energy density function as

S = 2 ∂ψ (C) ∂C . ( . )
Requirements of isotropy further dictate the strain energy functional to depend on C through its invariants. e three invariants of C are de ned as

I 1 = tr(C), ( .
)

I 2 = 1 2 (tr(C)) 2 -tr(C 2 ) , ( .
)

I 3 = det C. ( . ) Hence S = 2 3 i=1 ∂ψ ∂I i ∂I i ∂C . ( . )
Elastomers are treated to be incompressible materials and hence, det F = 1. As a consequence of incompressibility constraint, the PK stress can now be expressed as

S = 2 ∂ψ (C) ∂C + γ C -1 , ( . )
where γ is a Lagrange multiplier that enforces incompressibility. It coincides with the pressure, p, if the strain energy functional is homogeneous of order 0, ∂ψ ∂C : C = 0 (Bonet, ). is can be achieved by making it a function of the

deviatoric part of C, C -2 3 C. De ning ψ (C) = ψ ( C), S = 2 ∂ ψ (C) ∂C + p C -1 . ( . )
is can be further simpli ed as

S = 2 -2 3 DEV ∂ψ ( C) ∂ C + p C -1 , ( .
)
where DEV • = • -•:C 3 C -1 . For isotropic materials, ψ is now a function of the rst two invariants of C.

Hyperelastic models

An overview on the available hyperleastic models and the physics behind them can be found in the review articles of Marckmann and Verron, ; Verron, . A brief summary will be presented below.

In the Neo-Hookean hyperelastic model, the strain energy is a function of the rst invariant, I 1 , only (Treloar, ). is can be seen as an extension of Hooke's law to large deformations.

ψ N H = µ(I 1 -3).

( . )

is model has been developed by using statistical mechanics for polymer networks.

A model that involves both the invariants has been developed by Mooney (Mooney, ). For this model, the strain energy functional is wri en as

ψ M R = C 10 (I 1 -3) + C 01 (I 2 -3) ( . )
where C 10 > 0 and C 01 ≥ 0. e polynomial model has been developed by Rivlin (R S Rivlin, ), where the strain energy functional is expressed as a polynomial function of the invariants.

ψ P ol = ∞,∞ i=0, j=0 C i j (I 1 -3) i (I 2 -3) j , ( .
)
where C 00 = 0.

In Ogden model (Ogden, ), the strain energy density is expressed directly in terms of the eigen values of C instead of its invariants.

ψ O = n i=1 µ i α i λ α i 1 + λ α i 2 + λ α i 3 -3 ( . )
where µ i α i > 0 for all i.

More recent models include the -chain model by Arruda and Boyce, , limiting chain extensibility model by Gent, a, the extended tube-model (Kaliske & Heinrich, ) and many others.

. e theory of viscoelasticity can be easily described under small strain se ing. Under small strains, all the stress measures described in the previous section coincide. e strain measures are replaced by the small strain tensor, ϵ, which is de ned as the symmetric part of the gradient of the displacement.

ϵ = ∇u + (∇u) T 2 . ( . )
In the case of a linearly elastic material under small strains, the stress becomes a unique function of the strain. Roughly, in a viscoelastic material, the stress is seen to depend on the stain as well the strain rate, σ (ϵ, ϵ).

A viscoelastic material is typically represented by a spring-dash pot mechanical system in various con gurations such as by Kelvin or a Maxwell branch or a combination of both. eir behavior is characterized by time dependent material functions. More generally, the stress is made to depend on the current strain and its history. In di erential representation of linear viscoelasticity, the general equation relating the stress and strain can be wri en as

∞ n=0 u n ∂ n σ (t) ∂t n = ∞ m=0 q m ∂ m ϵ(t) ∂t m , ( . )
where u n and q m are constant coe cients. e above equation, wri en in an operator form, can be seen to be

U[σ (t)] = Q[ϵ(t)], ( . )
where U and Q are di erential operators. e di erence between the elastic and viscoelastic behaviors of the materials can be seen through Dynamic Mechanical Analysis (DMA) tests, where the stresses (reaction forces) are monitored on a test sample on which a cyclic strain is imposed. In the case of the elastic material, the stress is seen to vary in phase with the strain cycles. In a viscoelastic material, however, as a consequence of the time dependence of the material, the stress cycles are seen to lag the strain cycles, the lag depending on the frequency. is can be seen by taking the Laplace transform of equation . . e transformed equation can be seen to be ū(s)σ (s) = q(s)ε(s).

( . )

For suddenly imposed strain, ϵ(t) = ϵ 0 H (t), where H denotes the Heaviside function, the relation between the transformed stress and strain can be seen to be

σ (s) = Q(s) s ϵ 0 , ( .
)
where Q(s) = q(s) ū(s) . e quantity Q (s) s can be de ned to be Ḡ(s), the Laplace transformed relaxation function. is function can be de ned in the time domain to be

G(t) = L -1 Q(s) s . ( . )
Now, for a harmonic excitation, ϵ(t) = ϵ 0 exp(iωt), the transformed stress can be seen to be

σ = ϵ 0 Q(s) s -iω . ( . )
Steady state response can be obtained by se ing s = iω and inverting the transform.

σ ss (ω) = Q(iω)ϵ(ω). ( . )
e complex modulus can be de ned as G * (ω) = Q(iω). e real part of the complex modulus is termed as the storage modulus and the imaginary part as loss modulus.

For a standard linear solid, a er some algebraic manipulations, the storage and loss modulus can be seen to be

G s (ω) = G 0 1 - N i=1 γ i + G 0 N i=1 γ i τ 2 i ω 2 1 + τ 2 i ω 2 , ( . ) G l (ω) = G 0 N i=1 γ i τ i ω 1 + τ 2 i ω 2 .
( . )

An elastic material is characterized by τ i → ∞ in which case, the storage and loss modulus can be seen to be G 0 and 0 respectively. ey can be observed to be frequency independent.

e dynamic modulus is then de ned as

G d (ω) G 2 s + G 2
l . e loss tangent, δ , is de ned as δ arctan(G l /G s ). e typical variation of dynamic modulus and loss tangent for an elastomer can be seen in the gure . . e DMA experiments are thus used to characterize the viscoelastic behavior of the elastomer and calibrate the relaxation times and sti ness ratios. e dynamic modulus at small frequencies is called rubbery modulus to indicate the material behaves in a rubbery fashion at these frequencies. At larger frequencies, the modulus is called glassy modulus to indicate that the material behaves in a glassy fashion.

e loss tangent can be seen to increase from a small value at small frequencies to a maximum at the intermediate frequencies and decrease to a smaller value again at large frequencies. In the gure . , it can be seen that the loss tangent a ains a maximum at about ω = 10 9 Hz.

When in large strains, the e ects of viscoelasticity can be considered by either introducing an internal variable (Coleman & Gurtin, ) or by considering the current state of stress to be a function of current deformation state and the deformation history (Green & Rivlin, ; Coleman & Noll, ). A detailed discussion on the di erent approaches can be found in Gerard A Maugin and Muschik, ; Germain, Nguyen, and Suquet, and the references therein. In the framework of Standard Generalized Materials, the internal variables are evolved using a convex dissipation functional and its (Legendre-Fenchel) dual (Halphen & Nguyen, ).

Convolution integral formulation

A multiple integral constitutive equation has been proposed in Green and Rivlin, to consider the e ect of deformation history on the current state of stress. e state of stress at any time has been expressed as

σ (t) = F (t) G E(t -s)| ∞ s=0 F (t) T ( . )
Assuming that the functional G is continuous in the sense of Drapaca, Sivaloganathan, and Tenti, , a form of stress has been obtained in the form of multiple integral series (Wineman, ).

G E(t -s)| ∞ s=0 = ∫ t -∞ K 1 (t -s 1 ) dE(s 1 ) + ∫ t -∞ ∫ t -∞ K 2 (t -s 1 , t -s 2 ) dE(s 1 ) dE(s 2 ) + ∫ t -∞ ∫ t -∞ ∫ t -∞ K 3 (t -s 1 , t -s 2 , t -s 3 ) dE(s 1 ) dE(s 2 ) dE(s 3 ) + . . . ( . )

Finite Linear Viscoelasticity

Using the hypothesis of fading memory in Coleman and Noll, , the stress was made to depend on recent deformations. To account for history e ects, the relative deformation gradient, F t (τ ), has been introduced. F t (τ ) can be seen to describe the con guration of system at a time τ < t with respect to the current con guration at time t (Coleman & Noll,

). e stress at the current time is then expressed as a functional of the current strain and of the history.

σ (t) = F (t) k 2 [C(t)] + ∫ t -∞ K 2 [C(t), t -s] F (t) T (C t (s) -I )F (t) ds F (t) T ( . )
e rst term in the above expression can be taken to represent the 'long term' behaviour of the material. e second term characterizes the history dependence. It shall be noted that the use of the above from of stress is restricted to in nitesimal deformations if applied for a 'fast' process. Using it for large deformations restricts its application to 'slow' process.

e notions of 'fast' and 'slow' can be found in (Coleman & Noll, ). Many other formulations of linear and non-linear viscoelastic models exist. A review together with the notions of 'weak' and 'strong' fading memory can be found in Wineman, and Drapaca et al., .

Internal variable formulation

In the internal variable formulation of viscoelasticity, the history dependence accounted for by an internal variable. e usual approach in these models is to decompose the deformation gradient into an elastic and a viscoelastic part (F = F e F ). e internal variable is taken as the strain in the dash pot (F ) in a Maxwell element.

e evolution of the internal variable is represented by a di erential equation (Halphen & Nguyen, ). Many viscoelastic models have been formulated using this approach as in (Reese & Govindjee, ; Bergström & Boyce,

).

In this chapter, a brief introduction to the polymer materials was provided. Some hyperelastic and viscoelastic models have been presented brie y. Di erent approaches to model the inelastic behavior of polymers have been presented brie y as well.

In the current thesis, the fracture experiments performed in Polyurethane elastomers in Corre, will be analyzed rst. It will be determined if any shock-frontlike features are present in the material. e velocity and strain elds obtained from the experiments will be examined.

More experiments performed on the elastomer samples of di erent geometries will be presented next. ese experiments have been conducted to con rm the observations of Chen et al., for polyurethane elastomers. e importance of including viscous e ects while performing the analysis of dynamic fracture of elastomers will be discussed next. In that part of the thesis, the crack speeds will be imposed implicitly imposed on the model using the displacement elds extracted from the experiments.

A rate dependent cohesive zone will be presented next and will be used to predict the crack speeds.

e crack speeds will be obtained and compared for various stretch levels and geometries.

e non linear viscoelastic model by Reese and Govindjee, will be implemented next under plane stress conditions and used to compute the viscoelastic dissipation during the crack growth. e strain energy loss in the body will be compared with the viscous dissipation and the energy expended in the cohesive zone. . is chapter presents an overview of the experimental protocol followed and the results obtained in that thesis. e same setup will be used later on to conduct more tests on di erent geometries.

E XPERIMEN TS

. e setup can be found in Corre, , but is repeated here for convenience. e specimen has dimensions such that its length (l 0 ) is greater than its height (h 0 ) which, in turn, is greater than the thickness (e 0 ). It is held on its longest side in the jaws of a tensile machine and stretched to the target level using the supports molded onto it. See gures . and . . e strain state is homogeneous at the center, but the homogeneity is lost towards the right and le edges.

e deformation gradient at the center of the specimen a er being stretched (with respect to the coordinate system in gure . ) is

F = 1 0 0 0 λ 0 0 0 1 λ e x ,e ,e z ( . )
where λ = h f h 0 , h f is the nal height of the specimen a er being stretched. e specimens tested in Corre, were mm long, mm tall and mm thick. A razor blade is used to initiate a seed crack once the specimen is stretched. A stand is used to mount the razor blade, the height of which is adjusted to the required level depending on the specimen stretch.

. e specimen is initially painted in speckles that allow tracking the displacement elds of speckle particles using a DIC so ware. Once painted, the specimen is mounted on the tensile machine using the supports on its ends. Once secured, the top end is pulled at a speed of mm/min. Once, the target stretch level is reached, a seed crack is initiated using a razor blade. It shall be noted that the tensile machine continues to pull the specimen when the crack is propagating through it. However, as will be seen later, as the crack propagation process typically lasts for about ms, the displacement of the top end during this time is negligible.

e entire process is monitored by using two cameras -a High Resolution (HR) and a High Speed (HS) camera. e HR camera monitors the initial stretching part of the process, while the HS camera monitors the fracture part. e HR camera

Test sample Razor blade

captures the images at about frames/s and the HS camera at to frames/s. A trigger device is used to stop the HS camera once the crack completely passes through the specimen. e images from the two cameras are then used to obtain the displacement elds all through the crack propagation phase using the Digital Image Correlation (DIC) technique using VIC-D so ware ("VIC-D, Correlated Solutions", ). is procedure can be found in Corre, .

. e primary outputs of the experiments are the crack speeds. e use of the DIC technique enables to access the displacement elds all through the duration of the experiment. Other quantities such as velocities and strains can be derived from the displacement elds through their di erentiation in space and time.

e crack speeds obtained from the experiments in Corre, can be found below. Distinguishing the cracks into either subsonic or transonic requires the de nition of the shear wave speed. For linear elastic material subjected to small strains, this can be computed trivially as G ρ , G and ρ are the shear modulus and density, respectively. However, elastomers are so materials that can be subjected to large deformations and a non-linear model is necessary to describe their behavior at these strain levels. Hence, the expression for shear wave speed in a linear elastic material cannot be used for the current scenario.

Mooney-Rivlin model (Ronald S Rivlin & Saunders, ) can be used to describe the behavior of elastomers until a strain level of about %. e expression for shear wave speed for this material model has been obtained by Boulanger and Hayes,

. e cited reference obtains the speed of a nite perturbation that travels in a material described by the Mooney-Rivlin model. No assumption has been made in regards to the amplitude of perturbation. Hence, it can be used to obtain the speed of even large perturbations that are superposed on a base strain state. e strain energy density of the material in Mooney-Rivlin model is

ψ = C 10 (I 1 -3) + C 01 (I 2 -3), ( . )
where I 1 and I 2 are the invariants of the right Cauchy green tensor, C = F T F . C 10 and C 01 are the material parameters. e shear wave speed in this case can be obtained to be

c 2 s = 1 2ρ 0 C 10 n.bn + C 01 a.b -1 a , ( .
)
where b = F F T is the le Cauchy green tensor corresponding to the background deformation on which the perturbation is superposed, ρ 0 is the density of the material. n and a are the propagation and polarization directions respectively. For a shear wave, hence, n.a = 0. It shall be noted that the wave speed in the above expression is on the deformed state. e speed in the undeformed state can be obtained by using the expression wv = Fw 0 , where w and w 0 are the speeds in the deformed and undeformed con guration respectively. v is the particle velocity.

A ention should be drawn towards a couple of assumptions made during the development of the expression for shear wave speed in equation . . First, the background deformation state on which the perturbation has been imposed is assumed to be homogeneous.

e spatial derivatives of quantities like b have hence been ignored (see equations -of Boulanger and Hayes,

). Another assumption is that of the existence of a perturbation of the form f = h(x .nc s t) + k(x .n + c s t). Physically, this means that the form of perturbation should remain unchanged with respect to an observer traveling with the wave. It has not been discussed whether such a waveform can indeed exist. Saccomandi, mentions that in non-linear materials, solutions that retain their form may not exist. Perhaps a perturbation of the form f = h(x .nc s t) + k(x .n + c s t) can exist as long as it is in nitesimal.

e shear wave speed can be computed using the equation . . It can then be used to distinguish the crack speeds to belong to either subsonic or transonic regime. together with the shear wave speed can be seen in gure . . As will be discussed later, the shear wave speed seen in the gure . is slightly di erent from what has been obtained in Corre et al.,

. As can be seen, the general trend is for the crack speed to increase with the applied strain. At larger strains, the increase in crack speed is not as big as that at smaller strains for the same strain increment. e crack speeds remain smaller than c s for stretches smaller than . . Starting from a stretch value of about . , the crack speeds exceeded the shear wave speed. e former cracks are said to be Subsonic, while the la er is Transonic. e crack opening pro les for crack propagating in di erent regimes can be found in gures . and . of Corre, . It is generally known (Chen et al., ; Petersan et al., ) that the crack opening changes from parabolic to a wedge shape once the crack speed exceeds c s .

Some more data that can be obtained from that study is the variation of crack speed along the length of the specimen during the fracture process. e speeds in the deformed and undeformed con gurations can be seen in the gure . . e speeds in the two con gurations are related by we x = w 0 Fe X . It can be seen that the speeds in both the con gurations are almost the same at the center of the specimen where F takes the form

λ x 0 0 0 λ 0 0 0 1 λ x λ
, where λ x ≈ 1. At the edges, . Some more experiments were performed where the specimen was made nonhomogeneous by introducing a pro le into the specimen (see gure . ). e crack, initiated at the le end of the specimen, propagates at a higher speed till it reaches the geometric non-homogeneity. It slows down, propagates at a smaller speed in the region where the stretch level is smaller because of the inserted pro le and accelerates back to a di erent speed towards the right end.

In another case, the crack was initiated closer to the top of the specimen instead of at the middle between the two supports to study the path it takes (see gure . ).

e crack, a er a brief instance of traveling straight, turned towards the center of the specimen. e vertical speed of the crack gradually goes to zero once the crack reaches the center of the specimen. Some snapshots from the experiment can be found in gure . and the variation of horizontal and vertical components of crack speed with time can be seen in gure . . . In that article, the shear modulus was measured through the experiments and was then used to compute the shear wave speed. It was reported that the cracks traveled about -% faster than the shear wave speeds ). e experiments were performed on natural rubber specimen. e specimens were stretched to the target level and the crack propagation was studied. An additional observation was reported in that study that the crack speeds become independent of the specimen height in the Transonic regime. Meaning, the crack speeds were seen to be a function of specimen height for a given stretch level when the crack speeds are smaller than the shear wave speed. Once the crack enters the Transonic regime, its speed was seen to not depend on the specimen height. Meaning, in this regime, crack speeds in specimens of di erent geometry subjected to the same stretch level will remain same. Some similarity can be noted between this observation and that in Gent and Marteny, b (see gure . ). In Gent and Marteny, b, specimens of di erent heights are all subjected to the same stretch level and crack speeds were measured. It was observed that a er a certain specimen height, the crack speeds did not increase. . In the gure a , crack speeds were plo ed against applied stretch while in the gure b , crack speeds were plo ed against the Energy density, measured as wh 0 .

Whether the two observations are similar is unclear and will be discussed in the later chapters.

Further experiments in this direction were also carried out by Corre et al., ; Mai et al.,

. Transonic cracks were reported in these studies as well. However, experiments on di erent geometries have only been conducted in Chen et al., . Hence, more experiments have been performed on the specimens of Polyurethane elastomer to see if the phenomenon of height independence is observed here as well. ese will be reported in the next chapter.

In this chapter, the experiments performed in Corre et al., have been presented. Samples that are 40 mm tall were tested and Transonic cracks were observed starting from a stretch level of about . . More experiments were performed where the sample geometry was made non-homogeneous and where the crack initiation was made nearer to top edge of the sample.

e height independence of crack speed in Transonic regime wasn't investigated and it will be the point of departure of the next chapter.

To test the e ect of specimen geometry on the crack speeds, we tested Polyurethane samples of di erent geometries to obtain the crack speeds at di erent loads. is chapter will present the additional experiments performed on di erent geometries of the polyurethane elastomers. It also describes the tests performed using di erent test protocols. Temperature measurements performed on some test specimens will be presented as well.

. . ose will be reported below.

.

mm 60 mm

Batch e other specimens tested were mm and mm tall and belonged to a di erent batch (to be called batch ). e length ( mm) and the thickness ( mm) were the same as that of the mm samples. e test protocol followed was the same as that followed for the mm samples. e specimens had bulges molded into them which have been used to hold them in a tensile machine. ey were then stretched to di erent stretch levels and a crack is introduced. e initial stretching and the crack propagation phases were monitored with an HR and an HS camera as was done for the mm specimen tests. e specimens were pulled with the same speed while testing the and mm tall samples as the mm tall sample, that is mm/min.

As opposed to the tests in Corre, , in the tests performed on the new test samples, the test setup has been changed to include strain gauges at the le and the right ends of the support. See gure . . is enables to obtain the evolution of reaction forces (from the strains from the strain gauges) at the le and right ends separately as the crack propagates through the specimen.

e crack speeds obtained for these two specimen sizes can be seen in gure . along with the mm tall specimen results from Corre, . It can be seen that the overall trend is an increase in crack speed with an increase in the stretch level.

However, the trend is not monotonic. In some experiments, it was observed that the crack did not accelerate to a higher speed a er initiation. Instead, it continued to propagate at a speed that is smaller than that at smaller stretches. In some experiments, the crack speed varied drastically as it progressed along the specimen. It was also observed that the crack surfaces were considerably rougher in the places where the crack slowed down. e variation of crack speed along the length of the specimen in such cases can be seen in the gure . . e cases where the crack speed varied all through the propagation phase are indicated by black dots and circles in gure . . e crack speeds reported for such cases are for when the crack is at the center of the specimen, ie at about mm from the le end. It shall be noted that these experiments were performed only on the and mm specimens and the uctuations in the crack speeds were observed in both the geometries. e variation of reaction forces measured from the strains in the strain gauges as the crack propagated through can be seen in the gure . . e results reported correspond that presented in gure . ( mm tall specimen and λ = 3.0). As the crack begins to propagate and breaks a part of the specimen, the reaction force computed from the le strain gauge begins to drop. A er a li le delay, the force from the right strain gauge begins to drop as well. Corresponding to the large variations observed in the crack speed, oscillations can be seen in the force evolution as well.

e intent in performing the above experiments is to obtain the crack speeds in . However, it can be seen that because of the signi cant variations of crack speeds observed in and mm geometries, the speeds in those geometries could not be compared with the mm geometry as well as with each other. Hence, it was decided to procure another batch of samples, this time to test all the sizes.

. mm Since the mm tall specimens were already tested in Corre, , other kinds of experiments were performed on that geometry (still belonging to batch ). e test protocol is di erent from the previous section and is as follows. An initial seed crack of about mm is introduced in the specimen before it is stretched. e specimen is then held in the jaws of the tensile machine and is pulled at di erent speeds from the previous section. Tests were conducted at cross-head speeds of mm s -1 , mm s -1 , mm s -1 , mm s -1 and mm s -1 . It was intended to see if the crack initiation and speeds obtained depend on the initial stretching rate of the specimen.

It was observed that the crack did not propagate in a straight line like the other tests and instead, branching occurred at di erent locations. In some cases, one of the branches stopped propagating and the other propagated alone while in the other cases, both the branches propagated. Result for one such instance where the specimen is pulled at 30 mm s -1 and the crack is near the right end of the specimen can be seen in gure . . Similar observations have been made in the other tests as well, when the crack was at di erent locations. e locations where the crack a empted to branch in di erent experiments can be seen in gure . . However, the results could not be reproduced when the experiments were repeated multiple times. e experiments which resulted in crack branching earlier did not do so the second time under similar loading conditions. Some mm geometry samples were then tested under the same protocol as the rest of the other geometries ( and mm) to check if the properties of the specimen in batch-are same as that in Corre, . e specimen was stretched to a stretch level of . and a seed crack was introduced. An a empted branching was observed when the crack was at the center of the specimen! See gure . . It shall be noted that no such observations were made prior to this except for the studies of Moulinet and Adda-Bedia, . However, in Moulinet and Adda-Bedia, , experiments were conducted on balloons where the stress state is bi-axial and hence, there is a considerable stress along the crack path in addition to perpendicular to it.

. ( )

Since the results from batch-were not as expected, experiments were conducted on specimen of a di erent batch (to be noted as batch-). e test protocol followed is exactly the same as that in section the previous sections (as in Corre, ), where the specimens were stretched and a crack is introduced. e crack propagation was again monitored by an HR and an HS camera. e crack speeds obtained from this batch of specimen can be seen in gure . .

It can be seen that the crack speeds in this new batch are entirely di erent from those obtained from the batch-as well as in Corre,

. For instance, at a stretch level of . , the crack surface was considerably rougher and the crack speeds were considerably lower for all the geometries than in the previous batch. It might be possible that the specimen properties varied considerably between e arrow points to the location where the branching was a empted. Some branching events were successful while the others were not. An initial seed crack was present in all the cases. e cross-head of the tensile machine was pulled at di erent rates as indicated in the captions. di erent batches. No visible defects were present in the specimen.

e uniaxial response of the material from the new batch is determined by extracting the test sample from the 60mm tall specimen. e response can be seen in the gure . . e response of the sample in pure shear con guration can also be seen. It can be observed that starting from a strech level of about 3.0, the behavior of the material from the new batch considerably di ers from that of the previous batch in the uniaxial case. Similarly, in pure shear case, the behaviors of the two batches can be seen to be similar until about a stretch of 3.0, beyond which they were observed to diverge. is observed di erence in the material behavior may be one of the reasons for the observed di erences in the crack speeds between the two batches.

. In addition to obtaining the displacement elds in the specimen all throughout the experiment, a empts were made to measure the temperature change in the specimen (of batch ) as the crack propagated. Some other instances of measurements of temperature raise during the crack propagation can be seen in Fuller, Fox, and Field, ; D'Amico, Carbone, Foglia, and Galie i, . Two techniques were used in Fuller et al., for this. e rst method is called Liquid crystal lm technique, where a thin, even lm of the liquid crystal is spread on the surface of the specimen. e color of the liquid crystal is sensitive to temperature and hence, a change in the color of the layer was used to determine the temperature raise. is technique was used to measure the temperature raise of the PMMA specimen as a crack passed through it. e second method is using thermocouple junctions, where one junction was placed about . mm above the crack path while the other junction is placed remotely. is setup is used to measure the temperature raise during crack propagation. An Infrared detector was used in the same study as well. In D'Amico et al.,

, the temperature raise due to the crack propagation in a viscoelastic solid has been measured. e material used is SBR. An IR camera was used. A temperature raise of about 1 • C was reported. However, the crack speeds in that study are of the order of mm s -1 .

e temperature measurements were made on the samples in batch-. It was not possible to measure the temperature for many experiments because of limited number of samples. e measurement was done for the mm sample for an applied stretch level of 3. For this batch, the crack speed for this stretch is about m s -1 , which is much smaller than observed in Corre, . Hence, these temperature measurements should only be considered qualitatively. Due to the low acquisition rate of the IR camera, only one frame was obtained when the crack was propagating through the specimen. e results for that case can be in gures . and . . It can be seen from the gures that there is a large increase in temperature in the material at the vicinity of the crack tip. e temperature of the specimen during the stretching process changed from about to • C. is can be observed from the temperature di erence in the material far away from the tip in front of and behind it. e temperature of the material in the vicinity of the tip (and hence the crack faces) can be seen to rise to about • C. is can perhaps be a ributed to the viscoelastic dissipation in the vicinity of the tip as a consequence of high strain rates in that region. In the later chapters, an a empt will be made to identify and compare the region in which the material dissipates to the results just presented.

. e reason for performing additional experiments on samples of di erent sizes is to verify the phenomenon of the crack speed (in-)dependence on the specimen size in Transonic regime as reported in Chen et al.,

. However, as was presented, the crack speeds for the mm and mm (batch-) samples did not follow a monotonic trend as observed in Corre,

. A large variation of crack speed was observed in some experiments and there was a corresponding change in the appearance of the crack surface. e crack surface was observed to become rougher with a decrease in the crack speed. However, a quantitative relation between the two has not been obtained. Upon visual inspection, later on, it was observed that some specimens contained air bubbles throughout the specimen (and so, probably along the prospective crack path). Some specimens contained air bubbles that were initially invisible but became visible once they were stretched. is might be one of the reasons for the observed variation in the crack speeds. Di erent type of experiments were conducted on the mm specimen from batch-as data was already available for this geometry from Corre, . A seed crack was introduced into the specimen and the specimen was then pulled at di erent speeds. It was observed that the crack underwent branching in multiple cases. Since this was unexpected, the mm specimen was then subjected to the same protocol to check if similar results as that of Corre, can be obtained. However, a di erent behavior was observed. In Corre, , the crack, once initiated travels straight through the specimen. e crack speed remained constant in the middle of the specimen in most cases. In the current study, an a empted branching was observed when the crack is at the center. ree branches can be seen to originate in gure . . e top and bo om branches stopped propagating while the center crack moved on. As already mentioned, this is in contrast to the observations from Corre, . More tests, however, could not be conducted as there were no more samples available from batch-.

Another batch (batch-) of samples of all the geometries, , , and mm were procured to repeat the experiments, this time on the samples from the same batch. Experiments conducted on this new batch, however, resulted in crack speeds that are entirely di erent from the previous batches. At the same stretch levels, the crack speeds were considerably smaller than in the previous batches. e crack surface was considerably rougher as well. However, it can be seen from gure . that, in the experiments conducted on the new batch, the crack speeds in the three geometries se led down at a certain speed a er certain stretch level. Beyond a stretch of %, the crack speeds se led down at about m s -1 in all the three geometries. Perhaps, this can be compared with the results obtained by Chen et al.,

. e uniaxial and pure shear responses of the samples from batch di ered from that in Corre, (see gure . ). is may explain the observed di erences in the crack speeds.

e images taken from the HR and the HS camera during the initial stretching and the crack propagation for samples from batch-and batch-can be found at Kamasamudram and Coret, a, b. Temperature measurements have been taken on some of the samples from batch-. It has been observed that the temperature of the material in the vicinity of the crack tip rises to about 70 • C. e rest of the sample remains at about 26 • C. A change in temperature of the sample has been observed during the initial stretching process as well, which corresponds to the entropic elasticity of polymers. is e ect can also be observed when the crack is at the center of the specimen. In regions far from the tip, a change in temperature can be observed in front of and behind the crack tip. is temperature di erence is about 2 • C.

As a consequence of di erent results from di erent batches, the results from these two batches will not be considered for the analyses that will be performed in this thesis.

e additional experiments performed on di erent geometries and batches of polyurethane specimens have been described in this chapter. As a consequence of the variation of material properties between di erent batches and presence of defects in specimens, the results from the experiments could not be compared with those from Corre, . Additionally, temperature measurements have been taken on the specimens from batch-. A considerable increase in the temperature near and behind the crack top along the crack faces has been observed. e experiments performed in Corre et al., will be revisited in this chapter. e intent is to examine the velocity and the strain elds to look for any shock waves in the transonic regime. e variation of surface roughness with the crack speed in di erent experiments will be presented as well.

. Last, we note that when a moving object travels in an elastic medium faster than the wave speed, it creates a shock wave in the form of a Mach cone. e wedgelike crack opening typical of cracks in rubber (shown in Fig. ) is strikingly similar to the Mach cone, suggesting that rubber cracks are shocklike, exceeding some response speed of the medium.

is has been further studied in Marder, , where la ice models were used to simulate the crack propagation. A Neo-Hookean model was used to describe the hyperelastic part of the material behavior. Kelvin dissipation is added to describe the viscoelastic dissipation of the material. e author comments

Buehler et al. (

) have proposed that hyperelasticity plays a critical role in dynamic fracture, and that in particular an increase of sound speed near a crack tip can allow cracks to travel faster than the distant shear wave speed. In the analysis of this paper, there is an increase of elastic modulus near the tip of the rupture, but it is not in the form that Buehler et al. proposed. e theory here involves an increase in the imaginary modulus of rubber with frequency, rather than an increase in the real modulus of rubber with extension.

is is somewhat similar to what will be discussed in this thesis. However, a FE analysis will be used along with a continuum based viscoelastic material model.

Other studies of this nature were made in Chen et al., . Rubber sheets were used to perform fracture experiments in that study. However, as natural rubber tends to strain-crystallize at high strains, the crack could not propagate at those strains (Zhang, Niemczura, Dennis, Ravi-Chandar, & Marder,

). To circumvent this, the experiments in Chen et al., were conducted at 85 • C. High temperatures prevent the crystals from forming at higher strains. Numerical studies were also performed in that study using la ice models. Kelvin dissipation was introduced as in Marder, , but the model was calibrated based on retraction tests performed on rubber in Niemczura and Ravi-Chandar, . However, it shall be noted that the experiments in Niemczura and Ravi-Chandar, were carried out at 24 • C. It is not clear whether any adjustments were made to take the temperature di erence into account.

Other It shall be noted that in mode-II cracks propagating along weak planes in the Transonic regime, a shock front was observed that travels along with the crack (Rosakis et al.,

). e shock waves were clearly visible in the density of the photoelastic fringes observed. e crack speed was seen to exceed c s , reach c l and se le down at about √ 2c s . Particle velocities were seen to rise sharply when the Mach front passes through in the material (see gure of Mello et al., ). Photoelastic fringes con rm the passage of Mach front corresponding to the velocity jump observed.

In the studies performed on the fracture of elastomers, such shock fronts have not been reported. Instead, as already mentioned, the crack faces behind the tip were noted to be of wedge shape and resemble a shock front. Since no reports of a Mach front were found in the literature and since the displacement, strain and velocity pro les are available for the experiments in Corre, from DIC technique, they can be used to determine the presence or absence of the Mach fronts. It shall be noted that a jump in particle speeds should be observed across any such mach fronts. Hence, examining the velocity pro les for the presence of such jumps will help determine the presence of mach fronts.

. , e experiments in Corre, reveal that the crack speed at a stretch level of 1.7 is about 17m s -1 which is subsonic, while the crack speed at a stretch level of 3.5 is about 56m s -1 , which is transonic. Hence, examining the velocity pro les from these two experiments will help determine the presence, or absence thereof, shock fronts (Kamasamudram, Coret, & Moës, ).

. .

Horizontal displacement

Presence of a Mach front in experiments will result in discontinuity of velocity and strain elds. e displacement elds however, do not exhibit any jumps. However, for the sake of comparison later on with the simulations, the horizontal displacement elds for the two experimental cases of λ = 1.7, 3.5 will be presented. Only top half of the specimen is included.

e horizontal displacement elds can be seen for both the cases in gures . and . when the crack is about the center of the specimen. In both the cases, a distinct kidney-bean-shaped pro le can be observed in the material just behind the tip. A similarity in both the distributions can be seen even though the former crack is subsonic and the la er is transonic.

. .

Velocity elds e velocity elds for stretches . and . when the crack is at about the center of the specimen can be seen in the gures . a and . b. It can be seen that the velocity elds do not show a signi cant di erence in form between the two experiments. e peak value, however, is higher in the case of . as a consequence of a higher crack speed. However, the velocity is continuous.

To obtain a be er comparison, the velocity magnitudes are extracted at about mm above the crack path and plo ed in gure . . e velocity elds are continuous. A jump was not observed as was the case in the experiments in mode-II (Mello et al., ). And as mentioned earlier, the maximum value of particle velocity is di erent between the two cases because of di erent stretch levels and crack speeds. -6.69 -5.38 -4.07 -2.76 -1.45 -0.14 1.17 - Like in the case of velocity, a sharp jump could not be seen even for the deformation gradient. Similar to what has been done for velocities, the variation of yy-component of deformation gradient at about mm above the crack path in the undeformed con guration has been extracted and plo ed in gure . . It can be seen that F varies continuously along the specimen length for both λ = 1.7 and λ = 3.5.

An additional observation has been made from the fracture tests on the mm tall specimen. As mentioned earlier, a Pure Shear specimen is stretched to a required value and a seed crack is initiated using a razor blade on the le edge. e crack then runs from the le to the right edge. When the crack is just near the right edge (for example at about mm from the right edge), it has been observed that the specimen at the le does not go back to zero strain level -some residual strain is le behind (Figure . ), which eventually goes to zero a er about a few minutes to few days a er the experiment. e strain level at the le end at the instant described depends on the initial stretch and hence the crack speed. is indicates the presence of relaxation times that are larger than the duration of the experiment, which is typically about ms. e yy-component of the deformation gradient near the le edge are seen to vary between . and . for prestretch levels between . and . .

.

In the studies on the fracture of elastomers, another aspect of interest is the appearance of the fracture surface once the specimen is broken. roughout the literature, there are reports of changes in the appearance of the crack surfaces with crack speed and in some cases, the thickness of the specimen used for the study. One of the earlier studies that discuss this aspect is Greensmith and omas, . In that study, the propagation of a crack in an elastomer in 'tearing' con guration was studied and the process of crack propagation was classi ed into steady and stick-slip modes based on the reaction vs time behavior observed. A smooth surface was observed in the steady tearing mode, while a rougher surface was observed in the stick-slip mode. e stick-slip mode is a ributed to the presence of a toughening structure as a result of high strains at the tip, typically crystallization.

Another such observation can be found in Stevenson and omas, , where fracture of elastomers in balloon geometry was studied. e surfaces were reported to contain 'undulations' that are about -µm apart and -µm deep, but nothing more was discussed about them. Another such discussion can be found in Kadir and omas, , where again the fracture surfaces of specimen tested in the pure-shear con guration were reported to appear smooth but with small 'paraboloidal' marks pointing in the direction of propagation at higher crack speeds and a rougher surface at smaller crack speeds. However, Gent and Marteny, a reports smooth crack faces despite the smaller crack speeds observed in some load cases.

Cavitation was regarded as the reason for the change of surface roughness with the crack speed in Kadir and omas, . e smoother surface observed at higher crack speeds was a ributed to the absence of cavitation because of the 'viscoelastic sti ening' of the region in the vicinity of the tip as a consequence of high strain rates as the crack moves.

is theory was further investigated in Lake et al.,

, where an SBR sample was placed in a chamber lled with Ethylene Glycol which was used to apply external pressure to the sample. e crack propagation was studied under di erent applied pressures to study its in uence on the occurrence or the suppression of cavitation. e results, however, were reported to be unclear in that study for the lack of extensive data.

In Lake et al., , the phenomenon of cavitation was examined by studying the e ects of specimen thickness on the surface roughness. It was found that the specimen thickness a ected the crack speeds and the surface roughness at smaller crack speeds and that they remain independent of the thickness in the highspeed regime. e observed behavior in the low-speed regime was a ributed to the 'transverse stresses' that develop at the tip as a result of an increase in the specimen thickness which in turn promotes cavitation. Further studies in this direction were performed by Tsunoda et al.,

, where the e ects of specimen thickness and viscoelastic properties of SBR on the crack speed and surface roughness were examined. e SBR samples were swollen to various levels by Dibutyl Adipate to modify the viscoelastic properties of the material. In other cases, unswollen specimens with various thicknesses were used. In the case where the specimen is swollen, the surface was observed to become rougher at smaller crack speeds. In the case where specimens of di erent thicknesses were used, the crack speed was observed to change with thickness in the low-speed regime. High-speed regime was seen to be una ected by thickness. Both these results were a ributed to the presence (or absence) of cavitation.

When the fracture surfaces of the specimens broken in Corre, were examined using Scanning Electron Microscope, some variations in the surface roughness with crack speed were observed. e variations in the surface roughness observed can be seen in gures . and . . Before presenting the observations from these gures, it shall be noted that the thickness of the specimen that is used in this study is mm.

is value is higher than what has been used in the literature, like in Chen et al.,

; Gent and Marteny, a; Mai et al., ; Greensmith and omas, . Also to be noted is the fact that the surface pa erns presented in each gure are from the di erent regions of the same experiment as opposed to the cited literature where the roughness transition is observed between di erent experiments as the crack speed is varied. Other experiments that study the crack propagation in high speed regime do not report such an observation. Now, regarding the surface pa erns in the gures, the arrow in the gures . g and . g gives the direction in which the crack runs. Each sub gure is associated with a region in the specimen through an alphabet. In each of the sub gures, the direction of crack propagation is from right to le .

e results presented here are for the stretches of and where the crack speeds are about m s -1 and m s -1 respectively. In each case, the crack starts from the right end of the specimen, accelerates to some speed, which it maintains for most of the experiment. For the case of stretch level of , it was observed that the crack slightly decelerates as it approaches the other end, while for the case of stretch level of , it was seen to almost remain at the same speed (see gure . and Corre,

). e corresponding variation in the surface roughness can also be seen. In gure . , at the right end near the initiation, a very rough surface can be seen, probably from the razor blade used for the initiation. Once initiated, it can be seen from gures . b, . c and . d that the surface becomes smoother. During this phase, the crack speed can be seen to increase as well. However, as the crack approaches the le end, the surface can be seen to becomes rougher again ( gures . e and . f). e crack speed towards the le end can be seen to decrease for this stretch level.

For the case of λ = 3.0, a similar trend in the surface roughness can be observed in the initiation and acceleration phases ( gures . a, . b and . c). When the crack is at the center of the specimen, the surface is smooth at the magni cation level used. However, at the le end of the specimen, it can be seen that the surface remains almost smooth as opposed to the case when λ = 2.0. Correspondingly, the crack speed at the le edge for λ = 3.0 can be seen to remain at the same value as it was at the center. Between the two cases, it can be seen that the surface is considerably smoother for the case of λ = 3.0 than for λ = 2.0. A quantitative analysis has not been performed to link the roughness of the surface at a location (in terms of RMS roughness) with the crack speed at that location. One of the reasons for not having such observations from literature could be because of the smaller specimen thicknesses used in those studies. As stated at the beginning of this section, the specimen thickness is seen to in uence the crack speeds and the surface pa ern in the low-speed regime which would be the case when the crack starts from the rest at the right end of the specimen. us, an interesting question would be if the same pa ern will be obtained if a specimen of di erent (smaller) thickness is used for the current study or if the surface would be completely smooth in such a case.

. It can be seen from gures . and . that the velocity and strain vary continuously across the specimen even in the Transonic regime. A jump cannot be seen in either of these quantities. A di erence in the crack opening has been noted between the two regimes. In Petersan et al.,

; Corre et al.,

, it has been noted that the crack opening changed from a parabolic to a wedge shape as the crack propagation regime changes from subsonic to transonic. It has been remarked in Petersan et al., that the crack opening resembles a shock wave that propagates with the crack tip. is is, however, di erent from the shock front in the traditional sense where jumps in particle velocities and strains can be observed along a moving front within the material.

is is di erent from the Transonic cracks propagating in mode-II along weak planes. In studies where this crack propagation is studied (Rosakis et al., ; Rosakis,

), a mach cone can be clearly observed from the photoelastic fringes. In Mello et al., , the particle velocities are monitored just above the crack path, and a sharp rise in particle speeds was observed in the Transonic regime.

e experiments on elastomers do not report any such observations. Hence, the following questions might be asked (Kamasamudram et al., ) -Why don't we observe shocks in the bulk material even when a crack propagating in an elastomer travels faster than the shear wave (in the transonic regime)? Which material model more appropriately captures the behavior of the bulk material?

is chapter presents additional observations made regarding the velocity and strain elds together with the variation of crack surface roughness in the studies of Corre, . A variation of surface roughness has been observed with the crack speed. Jumps have not been observed in either the velocities or strains. is will be explored further in the following chapter.

S IMULATIONS

is chapter examines the phenomena of hyperelastic and viscoelastic sti ening to determine the phenomenon responsible for observing transonic cracks in polyurethane elastomers. Crack speeds are implicitly imposed on the model using the data from the experiments as described in this chapter. e predictions of the hyperelastic and the viscoelastic model will be compared with the experimental observations. e methodology presented in this chapter can also be viewed as a method of calibrating the viscoelastic model. In the previous chapter, it has been observed that the particle velocity and strain elds did not exhibit jumps in the Transonic regime. Traditionally, two main hypotheses are proposed for crack speeds in elastomers to enter the Transonic regime. e rst of these rely on the dependence of the wave speed on the strain state because of the material non-linearity (Gent & Marteny, a; Stevenson & omas,

). e stress-strain response of elastomers rst exhibits a drop in the slope at moderate strains followed by an increase at larger strains. e upturn in the stress-strain response then increases wave speed. Presence of a crack in the material results in a strain concentration near the tip. is implies that the material near the tip is subject to larger strains than the bulk material. Since the stress-strain response of the material exhibits an upturn at higher strains, this would in turn mean that the wave speeds in the vicinity of the crack tip are higher than that in the bulk material far from the tip. is hypothesis has been investigated in M. Buehler, Abraham, and Gao,

, where the propagation of a crack in a material whose modulus increases beyond a certain strain level (called the o set strain) is studied. Crack speeds that exceeded the Rayleigh wave speed (and also the shear wave speed) of the bulk material were observed depending on the values of the onset strain used. However, it was reported in that study that the crack tip is followed by mach cones corresponding to the shear waves.

In M. J. Buehler, , the problem of crack propagation along the interface of two solids with a large di erence (a factor of ) in sti ness was studied. It was observed that the crack propagated at speeds greater than the shear wave speed of the so er material. It was observed that the crack was able to propagate in the supersonic regime (with respect to the so er material) while remaining sub-Rayleigh with respect to the harder material. A mother-daughter propagation mechanism was observed similar to in mode-II. Shock fronts can be observed in the so er material ( gures . , . , and . of the reference). However, it was remarked that this hasn't been observed in the experiments conducted on the crack propagation through the interface in a bi-material system.

In Guo, Yang, and Huang, , the propagation of a mode-III crack in a solid that exhibits an upturn in the stress-strain response was investigated. In the vicinity of the crack tip, the material is sti er than the bulk. e distribution of shear stress in the bulk was presented in that study. e equations of motion lose their ellipticity in the region where the crack speed is higher than the shear wave speed. As a result, discontinuities in the stresses can be observed along with the characteristics.

e other hypothesis proposed for observing transonic cracks comes from Stevenson and omas, . In that reference, the author suggests that the higher crack speeds observed in that study are a consequence of the viscoelastic sti ening of the material as a result of the high strain rates present in the vicinity of the crack tip. A suggestion of this kind was also made in Marder, , where a la ice based method is used to simulate the propagation of crack in rubber. However, to the knowledge of the author, no continuum based numerical studies exist that examine this hypothesis.

e velocity and strain elds analyzed in the previous chapter did not show the presence of any jumps. e maximum value of velocity varied between the experiments but it varied continuously along the length of the specimen. Hence, in this section, an analysis will be performed to determine the appropriate material behavior responsible for transonic cracks. A part of this chapter is an extract from the article Kamasamudram et al., . .

In this section, the hypothesis of hyperelastic sti ening leading to transonic cracks will be examined. Since the displacement and velocity elds in the bulk material are obtained from the experiments, a hyperelastic model that exhibits an upturn in stress-strain response can be used to compute stresses and in turn momentum balance. A weak form of the momentum balance equation will be used for this purpose. It is intended to see if a hyperelastic model with an upturn in stress strain response, when applied for the crack propagation in transonic regime, will result in zero residual forces. e region where a non-zero residual is observed can be regarded as where the constitutive assumption is inappropriate. us, the hypothesis of hyperelastic sti ening can be examined. A weak form of the equilibrium will be developed below (Belytschko et al.,

). Beginning with the strong form of the momentum balance equation

div 0 P = ρ 0 a, ( .
)
where P is the rst Piola Kirchho stress, ρ 0 is the density in the undeformed con guration and a is the acceleration. De ning the space of admissible functions that are su ciently smooth to be V,

V = w that is su ciently smooth|w = o on ∂Ω D 0 , ( .
)

Ω D
0 is the part of the boundary on which the Dirichlet boundary conditions are described. e weak form of momentum balance equation can be wri en as

∫ Ω 0 v (div 0 P -ρ 0 a) dV = 0 ∀ v ∈ V.
( . )

Using the divergence theorem and that on Ω N 0 , Pn 0 = t, the above equation can be wri en as

∫ ∂Ω N 0 v.t dS - ∫ Ω 0 P : ∇ 0 v dV = ρ 0 ∫ Ω 0 v.a dV . ( . )
To discretize the above equation using nite elements, v is expressed in terms of its nodal values and shape functions. Hence

i = i I N I (X ), ( . )
where i = 1, 2, 3 indicates the spatial direction and I indicates the node number. Summation over repeated index is implied. Using this and also expressing the displacement eld in a similar way, the internal, external and inertial nodal forces can be wri en as

f int i I = ∫ Ω 0 ∂N I ∂X j P ji dΩ 0 , ( .
)

f ext i I = ∫ Ω N 0 N I t i dS, ( .
)

f kin i I = ∫ Ω 0 ρ 0 N I a i dΩ 0 , ( . )
and the momentum equation becomes

f ext i I -f int i I = f kin i I . ( . )
For the sake of discussion that follows, the residual is de ned as follows:

r i I = f ext i I -f int i I -f kin i I . ( . )
e computation of f ext i I and f kin i I require the external force applied and the time derivatives of displacement respectively which are readily available from the experiments and so, no constitutive assumption is needed regarding the material behavior. e computation of f int i I , however, requires the computation of P which depends on the constitutive behavior of the material. Hence, if an incorrect constitutive behavior of the material is assumed, non-zero residual forces will be observed.

Computation of stresses

e material is assumed to be perfectly incompressible. Hence, the strain energy depends only on the rst two invariants of the deviatoric part of the Cauchy-Green tensor. De ning C = -2/3 C, strain energy can be expressed as ψ = ψ ( Ī1 , Ī2 ). e PK stresses can be computed as (G. A. Holzapfel, ; Bonet & Wood, )

S = -pC -1 + 2 -2/3 DEV ∂ψ ∂ C , ( . ) DEV • • - 1 3 [• : C]C -1 , ( .
)
where p is the lagrange multiplier that enforces incompressibility condition. In plane stress conditions, p can be found by enforcing S 33 = 0. As a consequence of incompressibility and plane stress conditions, the third component of C can be de ned to be

C 33 = 1/det(C 2d ), where C 2d is the in-plane part of C. C = C 2d o o T C 33 ( . )
Instead of this, the strain energy function can be made to depend on the three invariants of C 2d where in the incompressibility condition is stricly imposed by taking

C 33 = 1/det(C 2d ) (Bonet & Wood, ).
is can be done by realiz-

ing that Ī1 = tr( C) = tr(C 2d ) + 1/det(C 2d ). Likewise, Ī2 = 1 2 Ī 2 1 -C : C = 1 2 tr(C 2d ) 2 + 2 tr(C 2d )C 33 -C 2d : C 2d .
In such a case, PK stress can be simply wri en as

S 2d = 2 ∂ψ ∂C 2d . ( . )
It shall be noted that nding p by enforcing S 33 = 0 and substituting back into the equation . will lead to the same result as in the equation . .

. . Calibration of hyperelastic model

A hyperelastic model that exhibits the upturn in stress-strain response is rst calibrated using the data from uniaxial tensile and pure shear tests. A polynomial model of third degree has been picked for this. e resulting model prediction vs the experimental data can be seen in the gure . . It can be seen that there is a slight deviation between the model prediction and the experimental data.

In the studies performed in Corre, , the hyperelastic model, Ogden model for instance, has been calibrated to match the uniaxial test data without any deviations ( gure . of Corre,

). Simulations have been performed using that model and parameters where a rectangular specimen of dimensions mm

times

mm is stretched to a level of . . e deformed shape of the specimen along with the deformed shape from the experiments can be seen in gure . a. A signi cant di erence can be seen. It was found to be necessary that the model coe cients be modi ed to bring the deformed shape closer to the experiments. e deformed shape predicted by the polynomial model with current coe cients can be seen in gure . b. ey can be seen to be quite close. Hence, these coe cients (in table . ) will be used in the current study.

. . Checking momentum balance

Using the hyperelastic model, the residuals in equation . can be computed for di erent experiments. As a recollection, if the hyperelastic model is correct and appropriate, the momentum balance equation should be satis ed exactly. is is rst carried out for the case of λ = 3.5. e residual forces are rst computed just before the crack is introduced into the specimen. Since the strain rates during this time are quite small, it is expected that the residual computed will be small throughout the sample. e residual forces thus computed can be seen in gure . a. Non zero and large forces can be seen on the boundaries of the specimen. e forces on the top and the bo om are the consequence of the applied loading and hence can be considered the reaction forces, acting along these boundaries, f ext i I = f int i I . Some forces can also be observed on the le and right edges of the specimen.

is is because of the lack of data along the boundary from DIC technique. ese forces can be seen to be the reaction forces applied by the missing material on the analyzed part of the material. Inside the material, small residual forces of about 0.02N (compared to that of about 17N along the top and bo om boundaries) can be observed. ese can be considered to be the consequence of the errors induced by the DIC technique and the constitutive model. e residual forces when the crack is at the center of the specimen can be seen in gure . . e scale has been truncated in gure . b. It can be seen that towards the right of the crack where the material is relatively stationary, the residual forces are comparatively smaller than in the material to the le of the crack. Similar observations can be seen in gure . , where the residual forces are computed similarly for the case of λ = 2.5. It shall be noted that higher residuals may also be observed in the material behind the tip even in the case of λ = 1.7, where the crack speed is subsonic. is is a consequence of the residual strains in the material a er the crack passed through (see gure . ).

It can be seen from gure . , where the vertical component of deformation gradient is presented for two experiments, that the material behind the crack tip does not completely relax to zero strain level even a er the crack tip has passed and the material is broken. is, along with a hyperelastic model, will result in non-zero stresses in this region since in a hyperelastic model, a nonzero strain will imply a non zero stress. is contributes to the observed residual forces in the region behind the tip. Also, the presence of crack tip in the material will result in strain concentration in the vicinity of the tip. Hence, the material sti ness and hence the wave speeds are higher in this region when compared to the bulk. is implies that the equations of motion remain their ellipticity in the region near the tip, but the ellipticity is lost in the far eld region (Guo et al., ). Hence, a shock-front-like feature along which the eld variables like velocity and strain exhibit discontinuities is expected. But since the elds observed in the experiments are continuous, this contributes to the residual forces as well. In short, the residual forces presented are a combination of contributions from non-zero residual strains in the region behind the tip and loss of ellipticity in the bulk away from the tip. .

In the analysis performed so far, the bulk material is assumed to be hyperelastic with an upturn in the stress strain response. e residual forces computed in this case contains a combined contribution of the residuals due to the strains behind the tip and due to the loss of ellipticity. To gain a be er understanding on the two hypotheses, another approach will be presented in this section. e e ect of including viscoelastic e ects in the bulk will be examined as well (Kamasamudram et al., ).

. . Methodology

e analysis methodology is as follows. Displacement elds are available from the experiments throughout the duration of the experiment. e intent is to use these elds to implicitly impose the crack speeds on the geometry and study the response of the bulk material by assuming di erent constitutive behaviors. e geometry of the specimen and the FE model used for this purpose can be seen in the gure . . e crack propagation path in the experiments is indicated by a dashed red line in the gure (in the undeformed con guration). Extracting the displacement elds from the experiments along this line and using them as (time varying) displacement boundary conditions in the FE model is equivalent to imposing the crack speed in the FE model. However, the DIC technique does not provide the displacement elds along the boundary. Hence, the extraction line is moved to just (about . mm) above the crack path. e new data extraction region can be seen as a green line in gure . . e top part of the specimen is constrained from moving horizontally but pulled vertically till the target stretch level is reached. Once the target stretch level is reached, the displacement is held constant on the top face.

At this point of time, it shall be noted that the crack propagation has not been explicitly modeled. e boundary conditions imposed on the model are similar to when the crack passes through the material (Kamasamudram et al., ). e same analysis can be made by using a cohesive zone model where the damage pro le is moved through the material. is approach, however, requires an additional assumption to be made regarding the behavior of the failing material. e approach discussed above bypasses the need to make a constitutive assumption regarding the behavior of the failing material.

e integration of the momentum equation in time has been performed using HHT-α time integration scheme (Hilber, Hughes, & Taylor, ). e simulations performed using value of 0 for α did not result in convergence. Hence, a value of -0.33 has been used. is value of α results in maximum numerical dissipation.

. . Simulations

Simulations have been performed using Finite Elements by the methodology prescribed above. Two cases have been examined. One for λ = 1.7, where the crack speed is about 17m s -1 and other for λ = 3.5 where the crack speed is about 56m s -1 . In the former case, the crack speed lies in subsonic regime while in the la er case, it is Transonic.

With hyperelastic bulk e polynomial hyperelastic model described and calibrated earlier has been used. As already mentioned, hyperelastic sti ening was deemed to be one of the reasons for observing Transonic cracks (M. Buehler et al.,

). e velocity pro le obtained from the simulation for the case of λ = 1.7 can be seen in gure . (to be compared with gure . a). A good agreement can be seen between both the cases.

e velocity pro le obtained from the simulation for the case of λ = 3.5 can be seen in gure . (to be compared with gure . b). A sharp rise in the particle speeds can be seen along a front, which corresponds to the existence of a Mach front.

To obtain a be er comparison with the experiments for the case of λ = 3.5, . , a sharp rise in particle speeds can be seen just behind the tip corresponding to the Mach front. e variation of particle speeds from the experiment can also be seen.

ere is a signi cant di erence between both the observations. Also, oscillations in particle speeds can be seen behind the mach front which is typical of the FE simulations containing a jump. Numerical Experimental e horizontal displacements in both the cases can be seen in gures . and . (to be compared with gures . and . respectively). Similar to what has been seen for velocity pro les, a di erence can be seen in the horizontal displacements between simulation and experiment for the case of λ = 3.5.

e results obtained are coherent with what has already been observed in the literature (Guo et al., ; M. Buehler et al., ). In both the cited studies, where the propagation of a transonic crack in a material with an upturn in stress strain response is studied, the crack was observed to be accompanied by Mach fronts along which the velocity and strain elds exhibit jumps. However, as already shown earlier, the elds obtained from the experiments do not exhibit such jumps. Hence, it can be concluded that at least for the current material, hyperelastic sti ening is not su cient for the cracks to enter the transonic regime.

With viscoelastic bulk

Since it has been shown that the upturn in the stress strain response is insu cient to reproduce the displacement and velocity elds in the bulk material, the simulations above will be rerun with viscoelastic bulk. A Finite Linear Viscoelastic (FLV)

*

-6.68 -5.37 -4.06 -2.75 -1.44 -0.14 1.17

•10 -3 -1.33 -1.03 -0.73 -0.43 -0.14 0.16 0.46

•10 -2 model (Simo, ) will be used for this purpose. Later on, the simulations will be performed again by using a Finite Viscoelastic (FV) model (Reese & Govindjee, ; Bergström & Boyce, ; Dal & Kaliske, ) in the bulk. In FLV model, the stresses can be expressed as a sum of a long term hyperelastic part and a viscous over-stress. e total stress in this case is

S = -pC -1 + -2/3 DEV t H , ( . )
H ∫ t 0 (t -s) ∂ ∂s DEV s 2 ∂ψ 0 ∂ C (s) ds, ( . ) DEV α • • - 1 3 [• : C α ]C -1 α , ( . ) (s) = ∞ + N i=1 i exp(-s/τ i ). ( . )
It shall be noted that the DEV t in equation . shall be evaluated using the quantities at the current time, while the DEV s in equation . shall be evaluated using the quantities at time s. In the equation . , ∞ and i s are the ratio of the modulus of the hyperelastic and viscous arms to the glassy modulus respectively.

ψ 0 is the strain energy corresponding to the instantaneous behavior of the material. Under small strain loading, the above constitutive model reduces to the linear viscoelastic model. In such a scenario, the model can be seen to represent the material by a spring that is in parallel to the Maxwell elements. An extension of such an interpretation of this model to nite strains is not possible. p in the equation is Lagrange multiplier that enforces incompressibility. In plane stress scenario, p can be found by enforcing that the third component of stress, S 33 is zero. See Appendix A. e model is initially calibrated by using the DMA tests performed on the material in Corre,

. When the simulations in the previous section are performed with that model parameters, the Mach fronts appearing in gure . disappear. However, the displacement and velocity elds were observed to be over-predicted. Hence, the model parameters were further adjusted to bring the results closer to the experiments. e adjusted model parameters can be seen in table . . e velocity eld predicted by the FLV model for the case of λ = 3.5, where the crack speed is transonic can be seen in gure . (to be compared with gures

τ i (s) E- E- E- E- E- E- E- E- E- E-
. and . b). e variation of particle speeds at about mm above the crack path can be seen in gure . . It can be seen that the velocity distribution obtained with a viscoelastic bulk more closely matches the experimental results than with a hyperelastic bulk. e displacement elds obtained from the FE simulations can be ). e results from FE simulations with a viscoelastic bulk can be seen to be similar to the experiments.

e particle velocity pro les plo ed on the deformed con guration can be seen Once the crack passes through the material, the faces behind the crack tip are (ideally) traction free. Since the displacements from the experiments are used as an input in the FE analysis, it is expected that the tractions predicted by the model along the crack faces will remain close to zero as well. e tractions resulting from hyperelastic and viscoelastic bulk can be seen in the gure . .

It can be seen that in the hyperelastic case, the stresses remain far from zero for a signi cant region behind the tip whereas in the viscoelastic case, the stresses a: Experimental result. Since the viscoelastic model is 'calibrated' using the experimental data for the case of λ = 3.5, it has been applied to the other experimental case (λ = 2.5) to verify its validity. e variation of particle speeds about mm above the crack path can be seen in gure . . When compared with gure . , a sharp rise in the particle speeds cannot be seen for the hyperelastic case. A smooth variation to the maximum value can be seen. However, the velocity variation observed in the It can be clearly seen that as in the case of λ = 3.5, the viscoelastic model predicts smaller stresses than the hyperelasic case. Hence, it can be concluded that the viscoelastic model 'calibrated' using the experiments of λ = 3.5 is able to predict the results of the λ = 2.5 with a good accuracy.

Sensitivity of results with respect to the data extraction line

In the analyses performed thus far, data are extracted from the experiments at about . mm above the crack path. Here, it will be determined if the results of the FE analysis obtained is sensitive to the position of the location along which the data are extracted from the experiments. For this purpose, another set of analyses have been performed where the data extraction from the experiments is done at about the midway (about mm above the crack path) of the half of the specimen.

is has been done for λ = 3.5. e results can be seen in gures . and . . In the results for horizontal displacements, a distinct kidney bean shaped pro le can be observed in the experimental results which can also be seen in the results from simulation with a viscoelastic bulk. Simulation with hyperelastic bulk does not exhibit this feature. Similarly, the velocity pro le obtained with a viscoelastic bulk can be seen to be closer to the experiments than with a hyperelastic bulk. is demonstrates the robustness of the viscoelastic model as well as the simulation methodology presented. . e previous sections contain the results from the experiments together with the results from FE simulations with hyperelastic and viscoelastic models. Experiments indicate that a 'shock front' was not observed even a er the crack speeds exceed the elastic shear wave speed -a jump has not been observed either in velocity or strain elds which would have been the case if a shock front was present.

e ). An interesting comparison can be made regarding the ratios of crack speed to the shear wave speed at a prestretch level of . in di erent studies. In Mai et al.,

, the speed ratio can be seen to be . , while in Chen et al., , the ratio is smaller than . . In the current study, the ratio is about . . is indicates perhaps an obvious conclusion that the material properties (mainly viscoelasticity) determine the extent to which the crack speed can exceed the shear wave speed.

e hyperelastic model calibrated using the uniaxial and pure shear cases is used to perform FE simulations with the displacements extracted from the experiments in Corre et al.,

. It shall be noted that the hyperelastic model does not exactly follow the Uniaxial stress-strain curve. e calibrated model, however, is seen to exhibit an upturn in the stress-strain curve which is usually deemed to be one of the reasons for crack speeds to go into Transonic regime (Gent & Marteny, a; Stevenson & omas, ; M. Buehler et al.,

). e FE result using this model, however, is seen to produce 'jumps' in both the velocity and strain elds for the cases when the crack speeds are Transonic -a result not seen in the experiments.

Including viscoelasticity is seen to bring the results closer to the experiments. e strain and velocity elds no longer exhibit the discontinuities that were observed in the hyperelastic case. e maximum horizontal displacements and the velocity magnitudes have been observed to be within % of the experiment results.

e crack faces behind the tip, being a free edge, are ideally traction-free. As presented in gure . , the residual strains behind the tip are non-zero. e computation of tractions from strains involves making a constitutive assumption about the material. Using a hyperelastic model in the simulations results in nonzero stresses along the edge as a consequence of the non-zero strain (despite it being a free edge in the experiments). See gure . , where the component of Cauchy stress is plo ed along the data extraction location ( . mm above the crack path). However, using a viscoelastic model with relaxation times that are comparable to the duration of the experiment, the stresses along the edge behind the tip quickly drop to small values within % of the far-eld stresses (and smaller than in the hyperelastic case). Re ning the model further can perhaps bring even those stresses to zero.

Hence, at least for the material in the current study, the cracks can go Transonic as a result of the viscoelastic 'sti ening' of the material. A suggestion of this kind can be found in Stevenson and omas, . However, no evidence has been presented in that study as to whether the Transonic cracks are a result of Viscoelastic sti ening or hyperelastic sti ening. A similar comment has been made in Marder, . However, in that work, a Kelvin type model has been assumed for the elastomer where the wave speeds increase unboundedly with the frequency. Hyperelastic sti ening has been deemed to be not necessary, but its su ciency has not been discussed. In studies where this su ciency has been studied M. Buehler et al., , 'Mach cones' corresponding to shear waves were reported that tail the crack tip (Figure of the reference).

e notion of limiting speed for a crack comes from the inclusion of inertial e ects in the study of the problem (Freund,

). e LEFM establishes the limiting speed as the Rayleigh wave speed in mode-I, while in mode-III, it is the Shear ; Willis, ), the limiting speed is seen to depend on the glassy modulus of the material rather than on the rubbery modulus (Figure . ). Hence, the crack is allowed to exceed the rubbery shear wave speed. Such a result can also be observed in Geubelle, Danyluk, and Hilton,

Rubbery Speeds | c R | c s | c d Glassy Speeds | A | c * R | c * s | B | c * d
, where spectral methods are used to investigate the problem of viscoelastodynamic mode-III crack. Using a rate independent cohesive zone type model, an analysis has been made on the e ect of relaxation times on crack speeds. From the Figure of the reference, it can clearly be seen that the crack speed exceeds the rubbery shear wave speed (c s ), but remains smaller than the glassy shear wave speed (c * s ). However, no speci c comments have been made in regards to the presence or absence of shock fronts in those references.

In Fisher and Gurtin, , the propagation of waves of order N in viscoelastic media has been studied. Waves of order N have the solution, u, which is N -1 times continuously di erentiable and exhibit discontinuities in the N t h derivative along a hypersurface. It has been determined in that study that such waves, should they exist, travel with a speed that is derived based on the glassy modulus of the material rather than the rubbery modulus. Shock fronts, by de nition, are waves of order and hence travel with that speed as well. Hence, it might be possible that shock fronts can be observed in viscoelastic material only when the crack speed exceeds even the glassy shear wave speed of the material. In short, a weak discontinuity, should it exist, would be made of high frequency waves at the sharp front. As a consequence of viscoelasticity, this region behaves with a high modulus which in turn roughly raises the shear wave speeds to a value that is (equal to or) greater than the crack speed. is, in turn, prevents a shock wave from developing.

Perhaps, the cracks whose speeds exceed c s while still remaining below c * s should be called r -Transonic (r stands for rubbery) to distinguish them from cracks that travel faster than c * s . One such instance may be found in Gori et al., , where a mode-II crack propagating in a PMMA specimen along a weak plane has been studied experimentally. It was observed that because of the viscoelastic behavior of the material, the crack speed exceeds c * s (hence c s ) and c d as well. Based on the inclination of the 'mach fronts' observed, the shear wave speed was computed and was seen to be in good agreement with the shear wave speed computed based on the modulus corresponding to the strain rates in that region.

As already mentioned, in the current chapter, the crack speeds are implicitly imposed on the material by using the displacement evolution from the experiments. In the next chapter, a cohesive zone model will be used to predict the crack speeds using the viscoelastic model calibrated in this chapter to represent the bulk material.

e crack speeds thus obtained will be compared with those from the experiments.

In this chapter, the hypotheses of hyperelastic vs viscoelastic sti ening for the crack speeds to enter the transonic regime were tested. e momentum balance check performed by using a hyperelastic sti ening showed the presence of large residual forces in the region behind the tip. is is expected since the material behind the tip was observed to not completely relax. e residual strains together with a hyperelastic model results in residual stresses as can be observed. However, the residual obtained from the analysis can be regarded as the combined e ect of the residual strains behind the tip and the 'loss of ellipticity' in the regions far from the tip. Hence, some more analysis were performed where the crack speed is used as an input in the model. e displacement eld data from the experiments obtained using DIC technique were used as boundary conditions in the FE model. is way, the crack speed and hence the crack propagation are implicitly modeled. Results from the FE analysis indicate that the hyperelastic model with an upturn in the stress-strain response resulted in a Mach front across which the particle speeds and strains exhibited 'jumps'. Results with a viscoelastic bulk closely match the experimental observations. e viscoelastic model was re ned further to match the results for one experimental case.

e same model parameters when used to simulate another load case predicted displacement and velocity elds close to the experiment demonstrating the robustness of the model and the 'calibration technique'.

In this chapter, a rate dependent cohesive zone model will be used to predict the crack speeds at di erent applied stretches. e Finite Linear Viscoelastic model will be used for the bulk material. e simulations will be performed at di erent geometries and di erent stretches to determine their e ect on the crack speed. . In other scenarios, the cohesive elements are inserted between all the elements of the FE model Miller et al., ; X. P. Xu and Needleman, to predict the crack path as well. In the la er scenario, some mesh dependency of the crack path was observed. e CZM technique involves specifying the behavior of the material undergoing failure. e output of the simulation is the crack speed or in some cases, the time to failure Schapery, c. CZM have been successfully used in the simulation of bri le fracture in the references cited earlier, in the case of plastic materials in Hutchinson and Tvergaard, and in the case of viscoelastic fracture in W. Knauss and Mueller, ; Schapery, b; Elmukash and Kroon, . Failure process in bri le elastic solids are usually rate independent -the fracture energy per unit area created remains constant. e fracture energy increases with crack speed as a consequence of an increase in the area of crack surface created due to micro branches (Sharon et al.,

). A correlation was found between the area of the crack surface and the crack speed once the crack speed exceeds about % of c R . It was a ributed to the development of micro branches that a empt to propagate along with the main branch, but are arrested. When the fracture surfaces of elastomers are examined, the opposite pa ern is observed. e crack surfaces tend to be rougher at smaller crack speeds and get smoother with an increase in crack speed. Crack surfaces of Polyurethane elastomer analyzed using SEM can be found in the earlier chapter con rm this observation. e variation in the surface roughness is usually a ributed to cavitation. At higher cracks speeds, the material in the vicinity of the crack tip sti ens as a consequence of viscoelastic e ects and hence the cavitation is suppressed. To the author's knowledge, no conclusive study however exists that con rms this e ect.

During the fracture of an elastomer, the total amount of energy dissipated is usually divided into dissipation in bulk (as a consequence of viscoelastic e ects) and in the fracture processes (in the creation of new surfaces). e fracture processes involve formation and coalescence of micro voids and chain scission. It is possible that the fracture processes themselves are rate dependent as a consequence of rate dependent e ects during fracture. Analyzing the peel tests of polymers, it was determined in Rahulkumar et al., ; Jagota, Bennison, and Smith, that the fracture processes are indeed rate dependent. A rate dependent cohesive model based on that of W. G. Knauss and Losi, has been used in those studies. Hereditary integral type model together with damage parameter has been used. Other studies that perform this type of analysis are Elmukash and Kroon, , . In those studies, it was found that the cohesive zone parameters are to be modi ed between di erent experiments to match the observed crack speeds with the experiment.

e rate dependent traction separation relation used in those works are based on a dash pot in parallel with a spring as opposed to the earlier mentioned hereditary integral type.

eoretical and experimental analysis of crack growth in viscoelastic materials have been performed by several researchers. A exhaustive review of such studies can be found in W. Knauss, . Knauss broadly classi es the studies into the ones that involve global energy balance in the srtructure (R. S. Rivlin & omas, ) and the ones that involve examining the crack tip elds and establishing criterion based on a cohesive zone like model (W. Knauss & Mueller, ; W. Knauss, ; Schapery, b, a, c) at the tip. In fact, there have been some discussions, such as in Kostrov and Nikitin, ; James R. Rice, , regarding the need to eliminate the singularity at the crack tip in viscoelastic fracture. e discussion goes as follows -a sharp crack model predicts singular stresses, strains and strain rates at the crack tip regardless of the speed of the crack. is would imply that the material near the tip behaves in a glassy fashion (as a result of the singular strain rates at the tip) even at low crack speeds. Also, using the failure criterion of Knauss and Schapery (Equation of W. Knauss, ), if the length of failure zone is taken to tend to zero, it can be shown that only the glassy modulus of the material enters the equation (page of Freund, ) and the overall viscoelastic properties do not seem to in uence the crack speed. ere have, of course, been some works that disagree with this analysis, for example Christensen, , that suggests that the discrepancy that arose is a consequence of the assumption that material behaves in a glassy fashion even far away from the tip. With a sharp crack model, even though the material near the tip behaves in a glassy fashion, a viscoelastic zone surrounds and follows the crack tip which controls the energy input and hence the speed of the crack. Other such instances of sharp crack model are B. N. Persson and Brener, ; B. N. J. Persson et al., , where the dissipation of the material has been explicitly taken into account. Hence, the whole spectrum of viscoelastic properties are seen to a ect the crack propagation behaviour. A point to be noted is that the cited references so far do not take the inertial e ects into consideration.

In dynamic fracture, LEFM bounds the speeds of the crack from below by Rayleigh Wave speed (c R ) for modes I and II, shear wave speed c S for mode-III (Freund,

). e theory is then extended to include the Transonic cracks observed during the laboratory experiments of mode-II cracks and during earthquakes (Rosakis,

). Some experiments performed on literature (Petersan et al ) to model the propagation of crack, an a empt will be made to reproduce the experimental results and understand the nature of fracture processes and dissipation mechanisms that accompany a moving crack. Similar studies have been performed by Kroon, ; Elmukash and Kroon, , ; Yin and Kaliske, ; Loew, Peters, and Beex, on fracture of elastomers using a Cohesive Zone or a Phase eld method. However, the emphasis of those studies are not to examine the transonic cracks or to investigate the e ect of geometry on crack dynamics in that regime. Infact, except for the analysis in Elmukash and Kroon, , , the crack speeds in the other studies are signi cantly smaller when compared to the elastic shear wave speed and hence, the inertial e ects do not play a signi cant role. e studies of Elmukash and Kroon, , use a rate dependent cohesive model to analyze the dynamic fracture of elastomers. But the parameters in the traction separation relations used in those studies are modi ed as the crack speed changes to match the experimentally observed crack speeds. In this work, an a empt will be made to reproduce the experimentally observed crack speeds using traction separation relation that relies on a set of parameters that do not change with crack speed. is will be done using the cohesive zone model proposed in W. G. Knauss and Losi, .

. e Finite Linear Viscoelastic model that was used in Kamasamudram et al., and in the earlier chapters has been used in this analysis as well. e hyperelastic part of the material is represented by a Polynomial model with N = 3 as

W = N i+j=1 C i j ( Ī1 -3) i ( Ī2 -3) j , ( . )
where Ī1 and Ī2 are the rst and second invariants of the deviatoric part of Green strain tensor, C = -2 3 C and is the determinant of the deformation gradient. C i j s are the model parameters. e model calibration has been done in Kamasamudram et al., and the parameters can be found in the previous chapters. Plane stress conditions are assumed to prevail. As a consequence, the stresses can be completely wri en in terms of the in-plane components of C. More details on this can be found in the next chapter (see equation . ).

e nite linear viscoelastic model of Simo, has been used for the bulk. e expression for the stress can be wri en as

S = -pC -1 + -2/3 DEV t H , ( . ) H ∫ t 0 (t -s) ∂ ∂s DEV s 2 ∂ψ 0 ∂ C (s) ds, ( . ) DEV α • • - 1 3 [• : C α ]C -1 α , ( . ) (s) = ∞ + N i=1 i exp(-s/τ i ). ( . )
e plane stress version of the above model can be obtained by obtaining the expression of p using the condition S 33 = 0. e expression of p when put back in the above equations gives the plane stress version of the constitutive model. is algebra has been performed in the appendix A. e nal expression for the stress can be seen to be

S 2d = ∫ t -∞ (t-s) ∂ ∂s DEV 2D s Ŝs + C t 33 3 (C t 2d ) -1 ∫ t -∞ (t-s) ∂ ∂s tr( Ŝs C s 2d )(C s 33 ) -1 .
(A. ) In the above, S 2d represents the restriction of the stresses to the in-plane components. e above equation can be integrated using the algorithm in Taylor, Pister, and Goudreau, .

. e experiments performed by Corre et al., have been analyzed in the previous chapter and in Kamasamudram et al., to study the importance of viscoelastic e ects in the bulk material in transonic regime. e crack speed has been used as an input in the analysis. is chapter aims at predicting the crack speeds by using a viscoelastic model in the bulk and a Cohesive Zone Model to predict the crack growth. Similar a empts were made earlier in Elmukash and Kroon, , where the simulations were performed to match the experiments performed in Gent and Marteny, a. A rate dependent cohesive zone model based on a Maxwell element was found to be necessary in those studies to be able to match the crack speeds in the experiment. e rate dependent Cohesive Zone was seen to be a consequence of the rate dependence of the failing material. However, in those studies, no particular emphasis was made on Transonic cracks or on their behavior in specimens of di erent geometries.

A rate independent cohesive model will be initially used here to con rm the observations from the earlier studies. It will be determined if the experimental results can be reproduced with this simple model using the same parameters for all the experiments. A model with initial sti ness and a linear unloading response will be used. e traction can be expressed in this model as

t = (1 -d)K∆, ( .
)
where d is the damage variable, K is the sti ness and ∆ is the separation vector. d is taken to depend on the current separation in the cohesive zone, ∆ [u], and is represented in the case of a linear cohesive law as

d =            0, if ∆ max < ∆ o ∆ f (∆ max -∆ o ) ∆ max (∆ f -∆ o ) , if ∆ max ∈ [∆ o , ∆ f ] 1, if ∆ max > ∆ f ( . )
where

∆ max (t) = max τ ∈[o,t ] ∆(τ )
is taken as the maximum value of separation over the entire time history to enforce the irreversibility constraint on the damage variable. ∆ o and ∆ f denote the value of separations at which the damage begins and the material completely fails, respectively. e energy dissipated in the cohesive zone can be found to be G = 1 2 K∆ o ∆ f . is cohesive zone model can be seen as a special case of a rate-dependent model, the implementation of which has been described later. Hence, only the results of this model will be presented in this section.

. . Specimen geometry and loading

e specimen that is 40 mm tall and 200 mm long has been analyzed. Plane stress conditions have been assumed to prevail. Cohesive elements have been inserted in the middle all along the length of the specimen. An element size of 50 µm has been used in the cohesive zone. Linear plane stress elements (CPS4) have been used to mesh the bulk region. e degrees of freedom of Cohesive elements are tied to be bulk elements.

e bo om part of the specimen is held xed while the top if stretched at a speed of 0.3 mm s -1 till the target stretch level is reached. A er reaching the stretch level, a seed crack that is about 3 mm long is imposed at the le end of the specimen (by reducing the ∆ 0 for those elements) which then propagates and breaks it into two pieces. e position of the crack tip is then extracted which is then used to compute the crack speed.

. .

Observed crack speeds with di erent cohesive zone parameters e parameters of the Cohesive model (∆ o and ∆ f ) are calibrated to match the crack speed at a stretch level of 3.5. e initial sti ness of cohesive elements, K is taken to be 1 × 10 7 MPa m -1 . It shall be noted that the crack speed in visoelastic fracture is a ected by not just the cohesive energy but also either the strain at failure or the fracture initiation stress (specifying one of them together with the fracture energy determines the other). For a given cohesive energy, di erent initiation stresses will result in di erent crack speeds. e e ect of cohesive parameters such as the strength, cohesive energy and length on steady-state viscoelastic crack growth has been studied in T. D. Nguyen and Govindjee,

. A similar study, in the context of peeling has been studied by Rahulkumar et al.,

. In both these studies the e ect of non-dimensional parameters based on the cohesive zone and bulk model has been studied. An analytical investigation has been made by D. B. Xu, Hui, and Kramer, . A cohesive strength of 27 MPa (or ∆ o = 2.7 • 10 -6 m) and the separation at failure of 2 mm was used for the analysis initially. is corresponds to an energy of 27 kJ/m 2 . is yielded a crack speed of about 58 m s -1 at a stretch level of 3.5, which is approximately the experimental value (56 m s -1 ). However, the same cohesive parameters gives a crack speed of about 5 × 10 -3 m s -1 at a stretch level of 1.7 while a crack speed of about 17 m s -1 was observed in the experiments.

Further a empts were made to decrease the cohesive strength and increase the separation at a stretch level of 3.5 to maintain the crack speeds close to the experiment. It was observed that the increase in the cohesive strength had a larger e ect on reducing the crack speeds than increasing the failure separation for the same cohesive energy. Hence, when the cohesive strength was reduced to 20 MPa (or ∆ o = 2 • 10 -6 m), it was found that the separation needed to obtain a crack speed comparable to the experiments was 3 mm. e crack speed was about 59 m s -1 . is, however, did not result in an increase in crack speed at a stretch of 1.7, which was still observed to be about 6 × 10 -3 m s -1 . Hence, the cohesive zone parameters were further modi ed, this time the cohesive strength to 15 MPa (or ∆ o = 1.5 • 10 -6 m) and critical separation to 5 mm resulting in crack speed at a stretch of 3.5 to be 59.5 m s -1 . No improvement was observed at a stretch level of 1.7 even at this point. It shall be noted that the cohesive energies varied for di erent cases as 27 kJ/m 2 , 30 kJ/m 2 and 37.5 kJ/m 2 while the crack speed was maintained in the range 58 ± 2 m s -1 . e parameter combinations tried can be seen in the table . . Table . : Crack speeds obtained for various cohesive zone parameter combinations.

∆ 0 (m) ∆ f (m) G coh (kJ/m 2 ) λ w 0 (m s -1 )
2.7 × 10 -7 2 × 10 -3 27 3.5 1.7 58 5 × 10 -3 2 × 10 -7 3 × 10 -3 30 3.5 1.7 59 6 × 10 -3

1.5 × 10 -7 5 × 10 -3 37.5 3.5 1.7 59.5 6 × 10 -3

It can be observed that using the same cohesive parameters across di erent loading does not result in the desired crack speed variation as observed in the experiments. Hence, rate dependence will be introduced in the subsequent sections to see if the same cohesive zone parameters can be used across di erent loading to obtain crack speeds closer to the experiments. e need for rate dependence in cohesive model has been arrived at by a similar analysis in Rahulkumar et al.,

; Elmukash and Kroon, as mentioned earlier.

. e rate dependence of the energy dissipated in polymer fracture is a ributed to the dissipation in bulk or the failure zone or both. e form of energy dissipated is usually taken as a product of an intrinsic energy (G 0 ) and a term that depends on crack speed, G = G 0 (1 + f ( ,T )) (Gent, b; B. N. Persson & Brener, ). Other forms for this dependence are available like in Schapery, b; W. Knauss, where the rate dependence factor is obtained from the rheological properties of the material. e total energy in these studies is wri en as a product of intrinsic energy and a term that depends on the crack speed (instead of the sum mentioned earlier). However, the intrinsic fracture energy or the energy that is dissipated at the tip during the fracture process is taken as a constant.

In the case of failure of adhesively bonded joints, the rate dependence is solely a ributed to the rate dependence arising due to the viscoelastic nature of the failing material alone. In such instances, the failure processes are taken into account by rate dependent cohesive zones (Musto & Alfano, ; Rahulkumar et al., ; Geißler & Kaliske, ). In some studies as in T. D. Nguyen and Govindjee, ; Rahulkumar et al.,

, the rate dependence is included in both the bulk as well as in the failure zone. e studies are however, restricted to the crack propagation in thin strips where the interaction of the bulk dissipation with the boundaries of the specimen is studied. In the context of dynamic fracture of elastomers, rate dependence in the cohesive zone is deemed to be necessary in the studies made by Elmukash and Kroon, , . In a similar manner presented above, the authors of the cited reference note that using a rate independent cohesive zone is insu cient to predict the crack speeds under all the loading conditions. Rate dependence in the cohesive zone was introduced through a Kelvin-Voigt element. No speci c link has been made between the rate dependence property of the cohesive zone and the bulk material. Also, no speci c a ention was paid to the transonic regime. In the context of phase eld modeling of polymer fracture, rate dependence of fracture energy was included in Loew et al., by the inclusion of damage rate in the damage evolution equation. More studies on polymer fracture using phase eld method have been made by Christian Miehe and Schänzel, ; Christian Miehe, Schänzel, and Ulmer, . However, the crack speeds in these studies are not high enough for dynamic e ects to become important. e rate dependence in the cohesive zone can be introduced in a number of ways. e cohesive strength can be made rate dependent while keeping the separation constant or the la er can be made rate dependent while keeping the former the same or both of them can be varied with the crack speed. Di erent criteria are known to result in di erent crack speeds as a consequence of viscoelastc bulk (W. Knauss,

). Some models used to prevent localization problems like delay-damage model (Su s, Lubrecht, & Combescure, ) can also be used to obtain rate dependent cohesive laws. In this study, a strain (crack opening in this case) based criterion will be used to predict damage initiation in the cohesive zone while keeping the failure to separation to be constant. is results in the maximum traction in the cohesive zone to change with the crack speed. As mentioned earlier, increasing the cohesive strength has been observed to result in a larger decrease in the crack speed while a larger variation in separation was found to be needed. e traction can be made a function of separations as well as their rates which is equivalent to introducing a dashpot in parallel with the spring that indicates a Kelvin type model (Geißler & Kaliske, ; Elmukash & Kroon, ). Or a convolution integral type law can be chosen where the traction can be made to depend on the (entire) separation history like in materials with memory (Rahulkumar et al., ; Zreid, Fleischhauer, & Kaliske, ). In W. G. Knauss and Losi, , the failure material is modeled by multiplying the stress strain law that is used in the bulk with a degradation function that depends on strain that results in the so ening of the material. In the current study, the la er model has been chosen and the kernels in the integrals (resembling W. G. Knauss and Losi, ; Rahulkumar et al., ) have been chosen to be same as that for the bulk material thereby establishing a link between the cohesive and the bulk properties. As will be seen later, this type of cohesive model is able to predict the crack speeds for the entire range of applied stretches and also di erent geometries.

e traction separation relation for the cohesive law has been chosen to be

t = (1 -d) ∫ t -∞ (t -s)K ∂∆ ∂s ds ( . )
Here, the memory kernel (t) = ∞ + N i=1 i e -t τ i . i s are the sti ness ratios while τ i s are the relaxation times. It is intended to see if the same sti ness ratios as that of the bulk material can be used in the cohesive zone as well. is way, cohesive zone can be regarded as the bulk material that is undergoing failure.

e evolution equation for damage is same as earlier as in equation . . But it can be seen that the traction in equation . depends not just on the current separation but its entire history. Since, storing all the separation data until the current time step demands large storage space, the integration is performed by taking advantage of the exponential dependence of the memory kernel using the technique presented in Taylor et al.,

or Simo and Hughes, . e integration technique is as follows.

From the equation . , the traction can be wri en as

t = (1 -d) K∆ + K N i=1 ∫ t -∞ i e -t -s τ i ∂∆ ∂s ds ( . )
It shall be noted that the i s in the above equation are to be scaled with respect to the rubbery modulus instead of the glassy modulus in the previous section. Each of the integrals in the parenthesis can be evaluated by using the Exponential mapping technique presented in Simo and Hughes, . Denoting the i th term to be h i (t),

h i (t) = ∫ t -∞ i e -t -s τ i ∂∆ ∂s ds ( . a) = ∫ t -∆t -∞ i e -t -s τ i ∂∆ ∂s ds + ∫ t t -∆t i e -t -s τ i ∂∆ ∂s ds ( . b) = e -∆t τ i h i (t -∆t) + ∫ t t -∆t i e -t -s τ i ∂∆ ∂s ds ( . c)
Computation of the second term e second term in the above expression can be computed in a number of ways.

Simo and Hughes (Simo & Hughes,

) propose integrating the term using a midpoint rule.

∫ t t -∆t e -t -s τ i ∂∆ ∂s ds ≈ e -t -s τ i ∂∆ ∂s s=t -∆t 2 ∆t = e -t 2τ i [∆(t) -∆(t -∆t)] ( . ) Hence ∫ t t -∆t e -t -s τ i ∂∆ ∂s ds ≈ e -∆t 2τ i [∆(t) -∆(t -∆t)] ( . )
In Taylor et al., , the second term is integrated assuming that strain rate, ∂∆ ∂s , is constant during ∆t.

∫ t t -∆t e -t -s τ i ∂∆ ∂s ds ≈ ∆(t) -∆(t -∆t) ∆t e -t τ i ∫ t t -∆t e s τ i ds ( . )
Evaluating the integral in the above equation,

∫ t t -∆t e -t -s τ i ∂∆ ∂s ds ≈ 1 -e -∆t /τ i ∆t/τ i [∆(t) -∆(t -∆t)] ( . )
In the case when the ratio ∆t/τ i in the equation . falls below 10 -7 , the expression

1-e -∆t /τ i ∆t /τ i is approximated by 1 -∆t 2τ i . e coe cients of the displacement increment in the equations . and . can be noted. e coe cients will be generally denoted as m i (t) for convenience from hereon. Hence

∫ t t -∆t e -t -s τ i ∂∆ ∂s ds ≈ m i (t) [∆(t) -∆(t -∆t)] ( . )
Using the above integration techniques in . c, the expression for h i (t) can be expressed as (no sum on i)

h i (t) ≈ e -∆t τ i h i (t -∆t) + i m i (t) [∆(t) -∆(t -∆t)] ( . )
e expression for traction can hence be found to be

t = (1 -d) K∆ + K N i=1 h i (t) ( . )
.

In this section, the Finite Element implementation of Cohesive model will be discussed. Beginning with the strong form of momentum balance equation

div 0 P = ρ 0 a ( . )
De ning the space of admissible functions that are su ciently smooth to be V, V = w that is su ciently smooth|w = o on ∂Ω D 0 , ( . )

Converting the momentum equation into weak form

∫ Ω 0 (div 0 P -ρ 0 a) δv dΩ 0 = 0 ∀ δv ∈ V ( . )
which a er simpli cations gives

∫ ∂Ω 0 t .δv ds - ∫ Ω 0 (P : ∇ 0 δv + ρ 0 a.δv) dΩ 0 = 0 ( . )
e boundary of the body can be divided into external boundary and the crack surfaces, ∂Ω 0 = ∂Ω ex t ∪ Γ ± 0 . Using this in the above equation

∫ ∂Ω e x t t .δv dS + ∫ Γ ± 0 t .δv dS - ∫ Ω 0 (P : ∇ 0 δv + ρ 0 a.δv) dΩ 0 = 0 ( . )
Since they are initially coincident, the integral on both the crack faces can be replaced by an integral on one crack face as

∫ Γ ± 0 t .δv dS = ∫ Γ 0 [δv].t dS ( . )
where [•] denotes the jump of •.

Using their values at nodes and the shape functions, δv can be expressed to be

δv = N i δv i ( . )
In one element, the jump in displacement can be expressed similarly as

∆ = [u] = [B]{u} ( . )
{u} is a vector containing the nodal values of quantities, expressed as

{u} = [u x 1 u 1 u x 2 u 2 u x 3 u 3 u x 4 u 4 ] T ( . )
and taken into consideration that the nodes 1 and 4, and the nodes 2 and 3 coincide with each other in the undeformed con guration. e term ∫ Γ 0

[B] = N 1 0 N 2 0 -N 2 0 -N 1 0 0 N 1 0 N 2 0 -N 2 0 -N 1 . (
[δv].t dS can be expressed using the nodal values of actual and virtual displacements as

∫ Γ 0 [δv].t dS = {δv} T ∫ Γ 0 [B] T t dS = {δv} T { f coh } ( . )
e integration in . is performed numerically as

{ f coh } = n p i=1 w i [B] T t | i ( . )
where w i s are weights and [B] T t | i indicates that the quantity [B] T t is computed at the corresponding integration point.

Using . , the jump in displacement can be expressed in terms of nodal displacements. e above equations can hence be expressed as

t = (1 -d) [K][B]{u} + [K] N i=1 h i (t) ( . )
h i (t) is computed recursively from the value at the previous time step as . ) where d{u} = {u(t)} -{u(t -∆t)}. e expression for traction can be further simpli ed as

h i (t) = e -∆t τ i h i (t -∆t) + i m i (t)[B] d{u} (
t = (1 -d)[K] [B] {u} + N i=1 i m i (t) d{u} + N i=1 e -∆t τ i h i (t -∆t) ( . )
Hence, tractions at the integration point can be evaluated by computing [B] and h i at those points. e m i (t) obtained from Simo and Hughes, is known to be sensitive to the time step size and hence the m i (t) from Taylor et al., will be used for this study. e tangent computation will be presented next.

. . Computation of residual

e weak form of momentum equation is integrated in time using the HHT-α (Hilber et al., ) time integration technique. Since a zero thickness cohesive zone has been considered for the analysis, there is no contribution to the mass matrix from the cohesive zone. e only computation from the cohesive comes from the term f coh . e residual on the weak form of momentum equation can hence be wri en as

r n+1 = -(1 + α)f n+1 coh + α f n coh . ( . )
e above equation together with the contribution from the bulk is solved using Newton's method in Abaqus. e above equation is hence linearized about the current displacement.

r n+1 k+1 = r n+1 k + ∂r k+1 ∂u ∆u = 0. ( . )
e tangent matrix can hence be seen to be

(1 + α) ∂ f coh ∂u .
. .

Computation of tangent

Tangent computation involves di erentiating . with respect to displacement.

∂ f coh ∂u = ∫ Γ 0 [B] T ∂t ∂u dS = ∫ Γ 0 [B] T ∂t ∂∆ ∂∆ ∂u dS = ∫ Γ 0 [B] T ∂t ∂∆ [B] dS ( . )
e traction derivative with respect to separation is computed as

∂t ∂∆ = (1 -d) [K] + [K] N i=1 ∂h i ∂∆ -t ∂d ∂∆ T ( . a) = (1 -d) (1 + N i=1 i m i (t))[K] -t ∂d ∂∆ T ( . b)
From . , the derivative of damage with respect to separation can be seen to be

∂d ∂∆ = ∆ 0 ∆ f -∆ 0 ∆ f ∆ 3 ∆ ( . )
e nal expression for the tangent can be seen to be

(1 + α) ∫ Γ 0 [B] T (1 -d) (1 + N i=1 i m i (t))[K] - ∆ 0 ∆ f -∆ 0 ∆ f ∆ 3 t∆ T [B] dS ( . )
Since only mode-I fracture is considered in this study, ∆ can be approximated by ∆ , the vertical component of the displacement jump.

. . Implementation into user subroutine

e above formulation of the cohesive zone has been implemented in to the user subroutine UELMAT of Abaqus (Dassault,

). e implementation details can be seen below.

It can be seen from the equation . that the computation of traction vector requires the computation of the jump vector, ∆ and the history variables, h i . e computation of h i at the current time step requires its value at the previous time step as well as the displacement increment between the two time steps. Hence, the history variables h i (of size 2 × 1) are stored in the state variable array (SVARS) that is supplied by Abaqus. e displacement increment between the two steps is supplied by Abaqus as the variable DU. In addition to this, the half step residual check from Abaqus at the end of each time step and the computation of residual for time integration require the force vector (of size 8 × 1) computed at the previous iteration. is is hence stored in the SVARS array as well. e damage variable is also stored to take the irreversibility constraint into consideration.

is computation is done at each integration point to get the right hand side of the force vector in the equation . . Since the current formulation uses 10 viscous arms, the total size of the state variable array at each integration point is 29 (= 1 (damage) + 8 (force) + 10 × 2 (internal variables)). Since two integration points are used for each element, the size of the state variable array is 58. e structure of the SVARS array can be seen in the gure . . e quantities to be returned to Abaqus is determined from the LFLAGS array of the subroutine. Per the documentation, the internal variables are updated when LFLAGS(3) = 1. In that case, the sti ness matrix and the force vector are computed at each integration point and are combined with the weights from the integration rule to get the overall sti ness matrix and the force vector, which is then supplied to Abaqus through the arrays AMATRX and RHS, respectively. e internal variables are updated using the equation . . . e implementation of the CZM has been rst checked by making it rate independent and comparing the results against the results from the COH2D4 element of Abaqus. In the rate independent case, the results from the two simulations are expected to match. is was tested for the case of λ = 2.5 and with ∆ 0 = 1 × 10 -6 m and ∆ f = 1.5 × 10 -3 m. Both the analyses resulted in the crack speed of about 59 m s -1 and similar displacement and velocity pro les.

mm specimen e cohesive model presented has been implemented into the UELMAT subroutine of ABAQUS. e model parameters that are to be calibrated are the separation at initiation (∆ 0 ) and the separation at failure (∆ f ). e material in the bulk has been taken to be the same as that in the previous section. As mentioned earlier, the i s and τ i s of the cohesive zone are taken to be the same as that of the bulk material.

e values of 1 × 10 -6 m and 1.5 × 10 -3 m for ∆ 0 and ∆ f , respectively, were seen to give a crack speed of m s -1 at a stretch level of . and a crack speed of about m s -1 at a stretch level of . , which is much larger than what was observed with a rate independent cohesive zone. e crack speed observed in the experiments varied between and m s -1 . Hence, the rate dependent model with the aforementioned choice of parameters was used to obtain crack speeds for the whole range of loading conditions to compare them with the experiments. e results can be seen in the gure . . It can be seen that the crack speeds from the simulation follows the experimental trend at least up to a stretch level of . . A slight upturn in the trend has been observed in the simulations a er the stretch level of . . e data point available from the experiments does not seem to follow this trend a er the stretch of . .

Specimen of other geometry

In crack propagation experiments performed by Chen et al., , it was observed that the crack speed becomes independent of the specimen height in the transonic regime. La ice model was used as well in that study to simulate the propagation of crack. In this section, it will be tested if the rate dependent cohesive model presented earlier will be able to predict the same result for polyurethane as well. For this purpose, the crack propagation in specimen of other heights such as mm, mm, mm, mm and mm will be studied. e lengths of the mm and mm specimen have been increased to mm to make sure that the steady state can be achieved when the crack is at the center. e results of the simulations can be seen in the gure . . It can be seen that at a given stretch level, the rate of increase in crack speed decreases as the specimen height is increased. In fact, it can be seen that the crack speeds saturate at some height and does not increase further with an increase in height. is can be observed if the crack speeds are plo ed against the specimen height instead of the stretch level (see gure . ).

As mentioned earlier, the crack speeds seem to saturate at a particular height and does not increase any further. a. In that study experiments were conducted on specimen of di erent heights all stretched to the same level. e crack speeds were observed to not increase once a speci c height level is reached as is the case in the simulations performed in the current study. Traction separation curves Since the traction separation law was made rate dependent, a variation in crack speed results in a variation in the traction separation law. It is expected that the maximum stress in the cohesive zone is larger than the cohesive strength in the rate independent case. e traction separation behavior of the cohesive zone is extracted from the simulations and plo ed vs the crack speed in gure . .

It can be seen that the maximum traction value increases rapidly to a value of about 21 MPa for a crack speed of 5.6 m s -1 . In the other cases where the crack speeds are about 41 m s -1 and 56 m s -1 , the critical stress can be similarly seen to increase to about 40 MPa and 46 MPa respectively. is, of course, leads to an increase in fracture energy as well. e fracture energies can be computed as the area under the traction separation curves. For the results presented in gure . , they can be seen to be 12.89 kJ/m 2 , 18.5 kJ/m 2 , and 21 kJ/m 2 , for λ = 1.7, 2.5, and 3.5, respectively.

Traction inside the cohesive zone

e variation of cohesive traction and the damage variable inside the cohesive zone for the three cases presented in the gure . can be seen in gures . and . . As can be seen from the gure . , the stress at the entry of the cohesive zone increases slightly with the crack speed. From its value at the entry, the cohesive traction increases to a maximum value and starts to decrease again. As the material in the cohesive zone starts to damage, the crack faces start to separate faster resulting in an increase in traction as a consequence of viscous e ects. However, as the crack faces separate, the damage variable increases and reaches 1 at the end of cohesive zone. Hence, it can be seen that the traction increases at the beginning of the cohesive zone as a consequence of viscous e ects and decreases once the damage variable increases. A similar trend can be observed for all the crack speeds 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6

•10 -3 In the gure . , the cohesive tractions are presented with the distance along the zone normalized with the length of cohesive zone. e length of the cohesive zone is measured in the undeformed con guration from the location where the separation reaches ∆ 0 till the location where it reaches ∆ f . e length of the cohesive zone can be seen to be 4.5 × 10 -4 m, 7 × 10 -4 m, and 7.5 × 10 -4 m for λ = 1.7, 2.5 and, 3.5, respectively. From the gure . , it can be observed that the normalized distance at which the cohesive traction a ains the maximum moves to the right with an increase in the crack speed. e damage pro le can also be seen to move towards the right with an increase in the crack speed.

Comparison of pro le of crack faces, displacement and velocity elds

e pro le of crack faces behind the tip obtained from the experiments and the simulations for the 40 mm sample for the case of λ = 3.5 can be seen in the gure . . e crack speed observed in the experiment is about 56 m s -1 and the speed obtained from the simulation is about 57 m s -1 .

e pro le obtained from the simulation is steeper than what was observed in the experiments. is could be a consequence of the viscoelastic model being 'less sti ' than expected at the strain rates experienced. However, as can be seen, the crack speed obtained from the simulation is about the same as in the experiments. e displacement and velocity elds when the crack is at about the center of the specimen can be seen in the gures . and . (to be compared with gures . a and . ), respectively. It can be seen that in the displacement and velocity elds over predicted in the simulations when compared with the experiments.

.

In this chapter, the crack propagation through the elastomer has been simulated using a cohesive zone model. Finite Linear Viscoelastic model has been assumed for the bulk material. Initial a empts were made to model the crack propagation using a rate independent cohesive zone. It was found that the cohesive zone parameters which result in crack speeds that match the experiment at λ = 3.5 either does not result in the crack propagation at λ = 1.7 or result in a very low crack speed -11.75 -8.99 -6.23 -3.47 -0.71 2.05 4.81

(orders of magnitude smaller than the experiments) at that stretch level. Hence, the entire range of crack speeds could not be predicted using the same cohesive zone parameters with a rate independent cohesive zone.

Studies performed in the literature using a cohesive zone model to predict the crack growth in polymers or adhesives (Rahulkumar et al., ; Geißler and Kaliske, ; Elmukash and Kroon, and the references cited therein) underline the need for a rate-dependent cohesive zone. e rate dependence in the cohesive zone can be a ributed to the rate dependence in the fracture processes as a consequence of viscoelastic e ects. A rate dependent cohesive zone has hence been employed in the current study as well. e rate dependence has been introduced by using convolution integrals as in the bulk. e relaxation times and the sti ness ratios of the viscous arms have been picked to be the same as that of the bulk. is just leaves two parameters in the cohesive zone to be calibrated, ∆ 0 and ∆ f . e parameters were chosen to reproduce the experimental crack speed at a stretch level of 3.5. en, using the same parameters, the cohesive zone model was able the entire range of the crack speeds observed in the experiments. See gure . . When used to predict the crack speeds for di erent geometries, it was observed that the crack speeds became independent of the specimen height at a given stretch level as was reported in Chen et al.,

. e height independence was, however, observed for specimen heights starting from about 60 mm. For specimens of smaller heights, the crack speeds still changed with a change in the specimen height for the same stretch. See gure . . e specimen heights tested in Chen et al., were 2 inch (50.8 mm), 4 inch (101.6 mm), and 7 inch (177.8 mm). Hence, it is possible that the phenomenon of height independence observed in that study was the consequence of the geometries used for the experiments. For instance, in this thesis, the geometries that are 20 mm, 40 mm, and 60 mm were tested by the authors. From the gure . , it can be seen that for this height range, the crack speeds have yet not become independent of the specimen height in the transonic regime.

e traction-separation curves in the rate dependent cohesive zone for di erent crack speeds can be seen in the gure . along with the rate independent case. It can be seen that the maximum traction in the cohesive zone is a function of crack speed and it increases to about 47 MPa for the case of λ = 3.5, where the crack speed is about 57 m s -1 . e variation of traction in the cohesive zone for di erent crack speeds can be seen in the gures . and . . In all the cases, the traction at the entry into cohesive zone can be seen to slightly increase with the crack speed. It reaches a value of 16 MPa for a crack speed of about 57 m s -1 .

e location of maximum traction can be seen to be behind the mathematical crack tip (where ∆ = ∆ 0 ). It can also be observed that the location moves away from the mathematical tip with an increase in the crack speed.

e maximum occurs at about x/l c < 0.2, ≈ 0.2, and 0.37 for three crack speeds presented. e results from the current study in gure . can be compared with the results from W. G. Knauss and Losi, in gure . . Only one case from W. G. Knauss and Losi, is presented here. All the cases can be found in the gure of that publication. In that study, the stresses at the entry of cohesive zone was found to vary non-monotonically with the crack speed. However, for range of crack speeds in the gure . , it can be seen that the stress at the entry of the cohesive zone is a weakly increasing function of crack speed as is the case in the current study. e maximum traction can also seen to increase with the crack speed till a certain crack speed. In the current study, the maximum traction increases with crack speed as well, and may also probably se le down for higher crack speeds (at higher stretches).

e cohesive model chosen in this study is intrinsic (Kubair & Geubelle, ) since it has a nite initial sti ness. e damage initiation criterion in this case has been described in terms of a critical opening displacement (d = 0 and d > 0 when ∆ = ∆ 0 and ∆ > 0). As the crack approaches a material point, the strain levels (and hence the separation) increases as a result of the strain concentration near the tip. is results in the stress at the entry of the cohesive zone being a function of crack speed as a consequence of the rate dependence. e cohesive element is taken to fail completely when d = 1, which happens at ∆ = ∆ f . However, if an extrinsic cohesive zone is used for the analysis, the separation remains at 0 till a certain critical stress is reached. In such a case, the damage initiation criterion used in this study may not be applicable. One alternative is to make the initiation stress as a (unknown) function of crack speed. However, as can be seen from the results of this study, even with a critical opening criterion, the stresses at the entry of the cohesive zone do not vary drastically with the crack speed (as compared to the maximum traction). Hence, the damage initiation stress in the case of an extrinsic cohesive model can be taken to be a constant while still retaining the rate dependent behavior from this study.

Another feature of interest is the length of the cohesive zone. e lengths are 4.5 × 10 -4 m, 7 × 10 -4 m, and 7.5 × 10 -4 m for λ = 1.7, 2.5 and, 3.5, respectively. In W. G. Knauss and Losi, , the cohesive zone length was reported to be a non-monotonic function of crack speed. However, at large crack speeds, the length was observed to increase with the crack speed (see gure . ). As mentioned already, in the current study, for the crack speeds considered, the cohesive zone length increases with the crack speed. If the failure of the material is taken to occur within the process zone, the roughness of the crack surface presented in the earlier chapter in the gures . and . can be taken to represent the variation of cohesive zone with the crack speed. Hence, this observation is in contrast with the result from the simulations.

e opening pro le of the cracked surface behind the tip from the experiments and simulation can be seen in the gure . . It can be seen that the pro le obtained numerically is steeper than the experimentally observed pro le. Also, the magnitude of velocity and the horizontal displacement can be seen in the gures . and . , respectively. e FE simulations can be seen to over predict these quantities when compared with the experiments. e reason for this can be described as follows. From the gures . and . , the reaction forces behind the tip do not go to 0 directly behind the tip. Instead, they go to 0 gradually over some length behind the tip. e magnitude of these stresses are positive as well, indicating that they are 'holding back' the crack faces from moving upward. In the simulations, however, since these faces are traction free, they are free to move until the material in front resists their motion. is can be seen to be the reason behind obtaining a steeper pro le in the simulations. As mentioned in the earlier chapter, if the viscoelastic model is further re ned to make the reaction in the gures . and . go to 0 over a smaller span, the pro le from the simulations may match closely with the ones from the experiments. e total energy expenditure in the bulk and the cohesive zone during the dynamic fracture, however, cannot be computed with the Finite Linear Viscoelastic model. Hence, the Finite Viscoelastic model has been implemented and this will be the focus of the next chapter. e energy expenditure during the viscoelastodynamic fracture will be examined. e portion of energy lost due to viscoelastic e ects of the bulk and the fracture processes will be examined.

A rate dependent cohesive model has been used to predict the crack speeds in Polyurethane elastomer. A convolution integral type model has been used for the cohesive zone, the coe cients of which are taken to be the same as that of the bulk material. e results from the simulations match the experimental results from Corre et al., for the 40 mm geometry. When applied for the case of other geometries, it was observed that for a given stretch level, the crack speeds did not depend on the specimen geometry starting from the specimen size of 60 mm.

is coincides with the observations from Chen et al., where the specimens of the natural rubber were tested. e evolution of traction within the cohesive zone, the maximum traction, the traction at the entry of the cohesive zone are similar to the observations from W. G. Knauss and Losi, . e area under the traction separation curves was observed to increase with the crack speed which is analogous to an increase in fracture energy with crack speed. To conclude, the hyperelastic model along with viscous contributions modeled in terms of Finite Linear Viscoelasticity and a convolution type cohesive zone was observed to predict crack speeds similar to the experimentally observed speeds.

In this chapter, the Finite Viscoelastic model (Reese & Govindjee, ; Bergström & Boyce,

) will be presented and its implementation will be described in a plane stress se ing. At rst, the model will be presented in a D se ing and the di erences arising in the implementation in plane stress will be discussed. A er this, some comparisons will be made between the current model and the Finite Linear Viscoelastic model. Finally, energetic analysis will be performed on the experiments and the simulations using the current model. e Finite Linear Viscoelastic (FLV) model has been used in the previous chapters to describe the viscoelastic behavior of the bulk material. However, the energy dissipated in the bulk material during the propagation of a crack cannot be computed with that model as it does not involve an explicit expression for energy or dissipation. To the author's knowledge, a plane stress implementation of the model by Reese and Govindjee, does not exist. Hence, the current chapter discusses its plane stress implementation and its use to compute the energy dissipated in the material during the propagation of a dynamic crack.

. In order to describe the processes that are far from equilibrium, a viscoelastic model has been proposed in Bergström and Boyce, ; Reese and Govindjee, in a thermodynamically consistent way. In this model, the strain energy density is taken to be a function of the deformation as well as some internal variables that characterize the rate dependence of the material. e strain energy density can be wri en as

ψ = ψ (C, Q 1 , Q 2 , ..., Q n ), ( .
)
where C is the right Cauchy Green deformation tensor and Q i are the internal variables. e evolution of internal variables are described by n equations of the form

Q k = f k (C, Q 1 , Q 2 , ..., Q n ). ( . )
e evolution equations and the expression for the internal energy shall satisfy the dissipation inequality

D 1 2 S : C -ψ ≥ 0. ( . )
Using a Maxwell type spring-dash pot model to represent the material, the strain energy function can be split into equilibrium and non-equilibrium parts as

ψ ψ EQ ( C) + ψ N EQ (C e ) = ψ EQ ( C) + ψ N EQ (F -T i CF -1 i ), ( . )
where C e is the elastic part of the deviatoric Cauchy Green tensor, C FT F . It has been assumed that F = -1 3 F admits a decomposition into an elastic and an inelastic part, F = F e F i (Dal & Kaliske,

). Using the above expressions in the internal dissipation inequality gives the expressions for stress as

S = -p C -1 + -2/3 DEV S EQ + S N EQ .
( . )

S EQ and S N EQ denote the elastic and viscous contributions to the total stress. DEV • represents the deviatoric projection and the rst term is a consequence of the incompressibility constraint. ese stresses are given by

S EQ 2 ∂ψ EQ ∂ C , ( . )
S N EQ 2F -1 i ∂ψ N EQ ∂C e F -T i , ( . )
and the deviatoric projector by

DEV • • - • : C 3 C -1 . ( . )
e Kirchho stress is then obtained as

τ = FSF T = -p I + τ iso = -p I + P : τ , ( . )
where τ = τ EQ + τ N EQ are de ned as

τ EQ 2 ∂ψ EQ ( b) ∂ b b, ( . ) τ N EQ 2 ∂ψ N EQ (b e ) ∂b e b e . ( . )
P is the deviatoric projector in the deformed con guration de ned as

P I - I ⊗ I 3 . ( . )
e dissipation inequality becomes

D = -τ N EQ : 1 2 L b e .b -1 e ≥ 0, ( . )
where L b e denotes the Lie derivative of b e de ned as

L b e = F C -1 i FT . ( . )
e expression . can be satis ed by specifying the evolution equation as

- 1 2 L b e .b -1 e = γ 0 V -1 : τ N EQ . ( .
)
where V is a fourth order isotropic positive de nite tensor possibly a function of b e and γ 0 > 0. A slightly di erent equation has been proposed by Bergström and Boyce, . In fact, the model of Reese and Govindjee, can be seen to be a special case of the model of Bergström and Boyce, .

. e integration of equation . is carried out by a predictor-corrector type algorithm. In the elastic predictor step, the inelastic strains are taken to be xed and so,

(C -1 i ) t r = (C -1 i ) t n-1 =⇒ b t r e = F (C -1 i ) t n-1 FT . ( . )
In the inelastic corrector step, the total deformation is assumed to be held xed and so, L b e = b e . Using this in equation . gives

b e b -1 e = -2γ 0 V -1 : τ N EQ . ( . )
e above equation can be integrated using the exponential mapping technique (Weber & Anand,

). e resulting expression is

b e = exp -2γ 0 ∫ t n t n-1 V -1 : τ N EQ dt b t r e , ( . )
(b e ) t n ≈ exp -2γ 0 ∆t(V -1 : τ N EQ ) t =t n b t r e . ( . )
e above equation is rst order accurate.

Since the material is assumed to be isotropic, b e , b t r e and hence τ N EQ share the Eigen space. Since V -1 1 2η D I -I ⊗I 3 is isotropic, equation . can be wri en in Eigen basis as

λ 2 Ae = exp - γ 0 ∆t η D dev(τ A ) (λ 2 Ae ) t r . ( . )
Taking logarithm of both sides,

ϵ Ae = - γ 0 ∆t 2η D dev(τ A ) + (ϵ Ae ) t r , ( . )
where ϵ Ae = ln λ Ae , (ϵ Ae ) t r = ln(λ Ae ) t r . e above equation is non-linear if τ A is a non-linear function of ϵ e . Hence, Newton iterations are used to solve it as below. De ning

r A ϵ Ae + γ 0 ∆t 2η D dev(τ A ) -(ϵ Ae ) t r = 0, ( . )
it can be solved by linearizing around ϵ Ae = (ϵ Ae ) k as ) where

r A + ∂r A ∂ϵ Be ∆ϵ Be = 0 =⇒ K AB ∆ϵ Be = -r A . ( .
K AB = ∂r A
∂ϵ Be . e above equation is solved to obtain ∆ϵ e , which is then used to update the elastic strain as (ϵ e ) k+1 = (ϵ e ) k + ∆ϵ e .

. As already mentioned, to the author's knowledge, a plane stress implementation of the FV model does not exist. is and the further sections discuss this implementation. It shall be noted that no changes to the model will be made. Rather, all the expressions for the stresses and the tangents will be rewritten so that they can be computed only using the in-plane components of the deformation gradient (F ) and its elastic part (F e ).

. . Stress and tangent computation

Stress computation

In plane stress scenario and for an incompressible material, the computation of stress can be simpli ed. e value of p in the equation . can be found by using the condition that τ 33 = 0. is condition can be imposed separately for the elastic and the viscous branches and the results can be combined. Beginning with the elastic branch, the term ∂ψ

∂ b b can be computed as ∂ψ ∂ b b = ∂ψ ∂I 1 I + ∂ψ ∂I 2 (I 1 I -b) b, ( . )
where the incompressibility of the material has been taken into account. e deviatoric projection of the above term is

P : ∂ψ ∂ b b = ∂ψ ∂I 1 I + ∂ψ ∂I 2 (I 1 I -b) b - 1 3 ∂ψ ∂I 1 I 1 + 2 ∂ψ ∂I 2 I 2 I . ( . )
e contribution of the elastic branch to p, denoted p e , can be wri en as (using

τ 33 = 0) p e = 2 ∂ψ ∂I 1 b 33 + ∂ψ ∂I 2 (I 1 b 33 -b 2 33 ) - 2 3 ∂ψ ∂I 1 I 1 + 2 ∂ψ ∂I 2 I 2 . ( . )
e total elastic part of the stress can then be found as

-p e I + 2P : ∂ψ ∂ b b = 2 ∂ψ ∂I 1 (b -b 33 I ) + 2 ∂ψ ∂I 2 I 1 (b -b 33 I ) -(b 2 -b 2 33 I ) . ( . )
e above equation, wri en with its components restricted to within the plane, can be seen to be

τ e = 2 ∂ψ ∂I 1 (b 2d -b 33 I 2d ) + 2 ∂ψ ∂I 2 I 1 (b 2d -b 33 I 2d ) -((b 2d ) 2 -b 2 33 I 2d ) . (
I 1 = tr(b 2d ) + 1/det(b 2d ), ( .
)

I 2 = 1 2 (I 1 (b 2d )) 2 -(b 2d : b 2d + 1/det(b 2d ) 2 ) , ( . ) the
S = F -1 τ F -T . ( . )
e above equations are restricted to in-plane components and the superscript 2d has been eliminated for convenience. e total stress in the equation . can be wri en for the case of multiple (say, N ) viscous branches simply as τ = τ e + N i=1 (τ ) (i) , where each of the (τ ) (i) s now denote the viscous stress in the corresponding viscous arm, de ned as (τ ) (i) 2 ) is the le cauchy green tensor in the i th viscous arm.

∂ψ (i) ∂(b 2d e ) (i) (b 2d e ) (i) . (b 2d e ) (i
e two invariants wri en in terms of principal stretches become

I 1 = λ 2 A + λ 2 B + λ 2 C = λ 2 A + λ 2 B + 1/λ 2 A λ 2 B , ( . )
I 2 = λ 2 A λ 2 B + λ 2 B λ 2 C + λ 2 C λ 2 A = λ 2 A λ 2 B + 1/λ 2 A + 1/λ 2 B , ( . )
where λ 2 C = 1/λ 2 A λ 2 B has been used. Using the above, the Kirchho stress can be wri en in principal basis to be

τ A = 2λ 2 A ∂ψ ∂λ 2 A = 2λ 2 A ∂ψ ∂I 1 ∂I 1 ∂λ 2 A + ∂ψ ∂I 2 ∂I 2 ∂λ 2 A , ( .
)

τ B = 2λ 2 B ∂ψ ∂λ 2 B = 2λ 2 B ∂ψ ∂I 1 ∂I 1 ∂λ 2 B + ∂ψ ∂I 2 ∂I 2 ∂λ 2 B , ( . )
and τ C = 0 as a consequence of the plane stress assumption. e partial derivatives of the invariants can be evaluated as

∂I 1 ∂λ 2 A = 1 -1/λ 4 A λ 2 B , ∂I 1 ∂λ 2 B = 1 -1/λ 2 A λ 4 B , ( .
)

∂I 2 ∂λ 2 A = λ 2 B -1/λ 4 A , ∂I 2 ∂λ 2 B = λ 2 A -1/λ 4 B . ( . )
e stresses can be expressed in global Cartesian basis by using

τ = τ 1 n 1 ⊗ n 1 + τ 2 n 2 ⊗ n 2 , ( . )
where n 1 and n 2 are the eigen vectors of b = F F T . As a recollection, the stresses obtained this way are the total stresses, τ = -p I + τ iso = -p I + P : τ .

Stresses can be similarly computed in viscous branches where λ is replaced by λ e , which are the eigen values of b e = F e F T e and the ψ replaced by the strain energy of the corresponding viscous arm.

Integration of evolution equation

For viscous branches, the evolution equation remains same even in the plane stress scenario. e residual can be wri en, similar to the equation . as

r A = ϵ Ae + γ 0 ∆t 2η D dev(τ A ) -(ϵ Ae ) t r = 0, ( . )
Since plane stress condition is assumed to prevail, only the in-plane components of the above equation are considered. Also, the deviatoric part of the Kirchho stress can be expressed as dev(τ ) = τ + pI , where p is the Lagrange multiplier that enforces incompressibility, which is found by using the condition that τ 3 = 0. τ can be evaluated by using the procedure in the previous section. e expression for p can be seen to be

p = - τ 1 + τ 2 3 . ( . )
e above equation can be obtained by taking the trace of the equation . and realizing that τ 3 = 0 as a consequence of plane stress assumption and that trace of the deviatoric projector is 0. e in-plane evolution equations then become

r 1 = ϵ 1e + γ 0 ∆t 2η D (τ 1 + p) -(ϵ 1e ) t r = 0, ( .
)

r 2 = ϵ 2e + γ 0 ∆t 2η D (τ 2 + p) -(ϵ 2e ) t r = 0. ( . )
In the third direction, the evolution equation becomes ϵ 3e + γ 0 ∆t 2η D p -(ϵ 3e ) t r = 0. It can be shown that solving the rst two equations exactly will result in the third equation being satis ed automatically. Adding the equations . and . will result in ϵ 1e +ϵ 2e +

γ 0 ∆t 2η D (τ 1 + τ 2 + 2p)-[(ϵ 1e ) t r + (ϵ 2e ) t r ] = 0.
is, in conjunction with the assumption of incompressibility and plane stress condition, results in

ϵ 3e + γ 0 ∆t 2η D p -(ϵ 3e ) t r = 0,
which is the third equation. e equations . and . are solved iteratively using Newton method.

r (k+1) A = r (k) A + ∂r A ∂ϵ Be ∆ϵ Be = 0. ( . )
is in turn requires the evaluation of K AB = ∂r A ∂ϵ Be . is can be evaluated as

K AB = ∂r A ∂ϵ Be = δ AB + γ 0 ∆t 2η D ∂τ A ∂ϵ Be + ∂p ∂ϵ Be . ( . )
e pressure derivative can be computed from equation . as

∂p ∂ϵ Ae = - 1 3 ∂τ 1 ∂ϵ Ae + ∂τ 2 ∂ϵ Ae . ( . )
e computation of derivative ∂τ A ∂ϵ Be can be carried out as shown in the following sections.

. . Tangent computation

For elastic branch e computation of tangent rst involves the computation of C = 2 ∂S ∂C . e derivative can be computed by noting that (Bonet, )

S = ∂S ∂C : C. ( . ) Since C = 2 i=1 λ 2 i N i ⊗ N i , C = 2 i=1 ∂λ 2 i ∂t N i ⊗ N i + λ 2 i N i ⊗ N i + λ 2 i N i ⊗ N i . N i = 2 j=1 W i j N j
, where W i j = -W ji are the components of a skew symmetric tensor. Hence

C = 2 i=1 ∂λ 2 i ∂t N i ⊗ N i + 2 i, j=1,i j W i j (λ 2 i -λ 2 j )N i ⊗ N j . ( . )
As a consequence of isotropy, S and C share the eigen vectors. Hence, following the same procedure,

S = 2 i, j=1 2 ∂ 2 ψ ∂λ 2 i ∂λ 2 j ∂λ 2 j ∂t N i ⊗ N i + 2 i, j=1,i j W i j (S i -S j )N i ⊗ N j . ( . )
e tangent can hence be wri en as

C = 2 i, j=1 4 ∂ 2 ψ ∂λ 2 i ∂λ 2 j N i ⊗N i ⊗N j ⊗N j + 2 i, j=1,i j S i -S j λ 2 i -λ 2 j (N i ⊗N j ⊗N i ⊗N j +N i ⊗N j ⊗N j ⊗N i ).
( . )

Its push forward to the spatial con guration can be seen to be

c = 2 i, j=1 (C i j -2σ i δ i j )n i ⊗ n i ⊗ n j ⊗ n j + 2 i, j=1,i j σ i λ 2 j -σ j λ 2 i λ 2 i -λ 2 j (n i ⊗ n j ⊗ n i ⊗ n j + n i ⊗ n j ⊗ n j ⊗ n i ), ( . )
where

C i j = ∂ 2 ψ ∂ ln λ i ∂ ln λ j = ∂τ i ∂ϵ j and σ i = τ i , since the material is incompressible.
e components of the above fourth order tensor can be stored in a matrix as

[c] =       1111 1122 1112 2211 2222 2212 1211 1222 1212       n 1 ,n 2 . ( . )
e components of the tangent can be converted to Cartesian basis by using the transformation (Reese & Wriggers, )

[c] (e 1 ,e 2 ) = [P][c] (n 1 ,n 2 ) [P] T , ( . )
where

[P] =       Q 2 11 Q 2 12 2Q 11 Q 12 Q 2 21 Q 2 22 2Q 21 Q 22 Q 11 Q 21 Q 12 Q 22 Q 11 Q 22 + Q 12 Q 21       . ( . )
Here, Q i j s are the elements of

[Q] matrix which is the transpose of [ Q], [Q] = [ Q] T .
e columns of [ Q] matrix are the components of eigen vectors of b in cartesian basis.

e tangent to be supplied to abaqus (C ( K ) ) corresponds to the Jaumann rate of the Kirchho stress (N. Nguyen & Waas,

) wri en as

∇ τ ( K ) = τ + τW -W τ = C ( K ) : d. ( . ) C ( K ) is related to c as C ( K ) i jkl = c i jkl + 1 2 (σ i j δ kl + σ kl δ i j + σ il δ jk + σ jk δ il ). ( . )
e computation of c requires the computation of ∂τ i ∂ϵ j , which can be carried out as follows

∂τ i ∂ϵ j = 2λ 2 j ∂τ i ∂λ 2 j , i,j= , . ( .
) ) e partial derivatives can be further evaluated as

∂τ i ∂λ 2 j = 2 ∂ 2 ψ ∂λ 2 j ∂I 1 λ 2 i ∂I 1 ∂λ 2 i + ∂ψ ∂I 1 ∂ ∂λ 2 j λ 2 i ∂I 1 ∂λ 2 i + ∂ 2 ψ ∂λ 2 j ∂I 2 λ 2 i ∂I 2 ∂λ 2 i + ∂ψ ∂I 2 ∂ ∂λ 2 j λ 2 i ∂I 2 ∂λ 2 i . ( .
∂ 2 ψ ∂λ 2 i ∂I j = ∂ 2 ψ ∂I 1 ∂I j ∂I 1 ∂λ 2 i + ∂ 2 ψ ∂I 2 ∂I j ∂I 2 ∂λ 2 i ( . ) ∂ ∂λ 2 A λ 2 A ∂I 1 ∂λ 2 A = 1 + 1/λ 4 A λ 2 B , ∂ ∂λ 2 B λ 2 A ∂I 1 ∂λ 2 A = 1/λ 2 A λ 4 B , ( . ) ∂ ∂λ 2 A λ 2 A ∂I 2 ∂λ 2 A = λ 2 B + 1/λ 4 A , ∂ ∂λ 2 B λ 2 A ∂I 2 ∂λ 2 A = λ 2 A . ( . )
Remark. e total stress and tangent computation in the case of plane stress condition for incompressible materials involves the computation of derivatives a er enforcing all the material (incompressibility) and geometric (plane stress) constraints.

In case of equal eigen values λ 1 = λ 2 , the second term of equation . takes a 0 0 form and so, L'Hospital's rule is used to compute it.

lim λ 2 →λ 1 σ 1 λ 2 2 -σ 2 λ 2 1 λ 2 1 -λ 2 2 = 1 2 ∂ 2 ψ ∂ϵ 2 ∂ϵ 2 - ∂ 2 ψ ∂ϵ 1 ∂ϵ 2 -σ 2 = 1 2 ∂τ 2 ∂ϵ 2 - ∂τ 2 ∂ϵ 1 -σ 2 . ( . )
For viscous branches For the viscous branches, the following procedure will be used to compute the tangent similar to as in Reese and Govindjee, . All the stress and strain components are now restricted to within the plane.

In the elastic trial state, since the inelastic strain is held xed,

F n = F t r e F n-1 i =⇒ C n = (F n-1 i ) T C t r e F n-1 i . Hence, ∂S I ∂C K L = ∂S I ∂(C t r e ) α β ∂(C t r e ) α β ∂C K L = (F n-1 i ) -1 K α (F n-1 i ) -1 Lβ ∂S I ∂(C t r e ) α β . ( . )
where the symmetry of C t r e has been used. Since in viscous branches,

S = F -1 τ F -T = (F n-1 i ) -1 . (F t r e ) -1 τ (F t r e ) -T S .(F n-1 i ) -T , ( . )
the stress derivative can be further re ned as

∂S I ∂(C t r e ) α β = (F n-1 i ) -1 Iγ (F n-1 i ) -1 δ ∂ Sγ δ ∂(C t r e ) α β . ( . )
Hence,

2 ∂S I ∂C K L = 2 (F n-1 i ) -1 Iγ (F n-1 i ) -1 δ (F n-1 i ) -1 K α (F n-1 i ) -1 Lβ ∂ Sγ δ ∂(C t r e ) α β
. ( . ) e push-forward of above by F results in

c i jkl = 2F i I F j F k K F l L ∂S I ∂C K L = 2(F t r e ) iγ (F t r e ) jδ (F t r e ) kα (F t r e ) l β ∂ Sγ δ ∂(C t r e ) α β . ( . )
S, wri en in Eigen basis is

S = 2 A=1 τ A (λ Ae ) 2 t r Ñ A ⊗ Ñ A .
( . )

It is to be noted that the τ A in the above equation is a function of ϵ e . e development from here is similar to that used to arrive at the equation . except that C i j in that equation will be replaced by C al i j , which will be de ned below. A crucial factor is that b e and (b e ) t r share the same eigen space as a consequence of isotropy.

c = 2 i, j=1 (C al i j -2σ i δ i j )n i ⊗ n i ⊗ n j ⊗ n j + 2 i, j=1,i j σ i (λ j ) 2 t r -σ j (λ i ) 2 t r (λ i ) 2 t r -(λ j ) 2 t r (n i ⊗ n j ⊗ n i ⊗ n j + n i ⊗ n j ⊗ n j ⊗ n i ). ( . )
In the above,

C al AC = ∂τ A ∂(ϵ Ce ) t r
. Since τ A s are a function of ϵ e s, the derivative is computed using chain rule.

∂τ A ∂(ϵ Ce ) t r = ∂τ A ∂ϵ Be ∂ϵ Be ∂(ϵ Ce ) t r . ( . )
e derivative ∂ϵ Be ∂(ϵ Ce ) t r can be computed by realizing that the equations r B = 0 are satis ed at the end of local Newton iterations and hence are valid at all the times during the global Newton iterations. Hence, during the global Newton iterations, ∂r B ∂(ϵ Ce ) t r = 0 as well. Hence,

∂ϵ Be ∂(ϵ Ce ) t r = K -1 BC , ( . )
where K BC is de ned in equation . . e expression for dissipation becomes

D = 1 η D dev{τ N EQ } : dev{τ N EQ } = (τ + pI ) : (τ + pI ) + p 2 , ( .
)
which can be seen to be positive since η D > 0. In the above equation, τ is as evaluated in the equation . .

. . Implementation details

e above has been implemented into a UMAT subroutine of Abaqus. e subroutine computes the updated stress and tangent for a given time step and also updates the internal variables. e implementation details will be discussed in this section.

e deformation gradient is obtained as an input to the subroutine from Abaqus. e internal variables (C i ) at the beginning of the step are stored in an array. e working of the subroutine can be seen below.

As can be seen, the internal variables are updated and the strain energy density and dissipation are computed with the updated values of b e and τ . e value of TOL has been chosen to be 10 -5 .

e C++ implementation of the model can be found at Kamasamudram, and in the appendix C.

. . . Hyperelastic case

e correctness of the implementation of the material model in the UMAT subroutine will be checked rst. e polynomial strain energy functional calibrated in the

Data: F n , C n-1 i , ∆t Result: τ n , C n i , C (J K) , ψ n , D n (b e ) t r = F n (C n-1 i ) -1 (F n ) T ; Compute (λ e ) t r s and n A s so that (b e ) t r = 2 A=1 (λ Ae ) 2 t r n A ⊗ n A ; De ne (ϵ Ae ) t r = ln((λ Ae ) t r ); k = 0, ϵ Ae ← (ϵ Ae ) t r ;
do Compute τ and p from ϵ Ae ;

r A ϵ Ae + ∆t 2η (τ A + p) -(ϵ Ae ) t r = 0; Compute ∂τ A ∂ϵ Be , ∂p ∂ϵ Be ; K AB = δ AB + ∆t 2η ∂τ A ∂ϵ Be + ∂p ∂ϵ Be ; ∆ϵ k Ae = -K -1 AB r A ; ϵ k +1 Ae ← ϵ k Ae + ∆ϵ k Ae ; k ← k + 1 while ||r A || > TOL; Update τ n , p, λ Ae = exp(ϵ Ae ), and K AB ; b e = 2 A=1 λ 2 Ae n A ⊗ n A , b -1 e = 2 A=1 λ -2 Ae n A ⊗ n A ; C n i = (F n ) T b -1 e F n ; Compute C (J K) using C al , τ A ,
and λ Ae ; Compute ψ n and D n ; Algorithm : Steps followed in UMAT for Finite Viscoelastic model chapter will be used. e results of the predictions from UMAT will be compared with that from Abaqus to check the sanity of the implementation.

e geometry along with the mesh can be seen in the gure . . e top end of the sample is subjected to displacement u x = 40 mm, and u = 50 mm, while the bo om end is held xed. e 22 component of Cauchy stress, σ predicted by the UMAT and the internal Abaqus implementation can be seen in gure . . It can be seen that they are identical. e other components of the Cauchy stress have been found identical as well and will be presented in the appendix B. . . Viscoelastic case -small strains and small perturbations Under small strain conditions, the deformation gradient can be approximated by F ≈ I +ϵ and C ≈ I +2ϵ. e principal stretches can be approximated as λ 2 1 ≈ 1+2ϵ 1 and λ 2 2 ≈ 1 + 2ϵ 2 , where the ϵ i s are the eigen values of ϵ, the small strain tensor. e small perturbation assumption allows similar approximations for the viscous branches as well. F e ≈ I + ϵ e and C e ≈ I + 2ϵ e . Similarly, λ 2 1e ≈ 1 + 2ϵ 1e and λ 2 2e ≈ 1 + 2ϵ 2e . Plane stress conditions are assumed to prevail. Material is taken to be incompressible. In that case, the stresses in the principal directions in the equations . and . become

τ A ≈ 2C 10 (4ϵ A + 2ϵ B ) , (
r 1 = ϵ 1e + ∆t 2η D 4C 10 ϵ 1e -(ϵ 1e ) t r = 0 =⇒ ϵ 1e = (ϵ 1e ) t r 1 + ∆t τ , ( .
) ) where η D 2C 10 τ has been used. e above equations can be seen to be the backward Euler discretization in time of the PDE

r 2 = ϵ 2e + ∆t 2η D 4C 10 ϵ 2e -(ϵ 2e ) t r = 0 =⇒ ϵ 2e = (ϵ 2e ) t r 1 + ∆t τ , ( .
ϵ Ae + ϵ Ae τ = ϵ A . ( . )
is can be proved by using ϵ Ae =

ϵ n Ae -ϵ n-1 Ae ∆t and ϵ A = ϵ n A -ϵ n-1 A ∆t
. e equation . then becomes

ϵ n Ae -ϵ n-1 Ae ∆t + ϵ n Ae τ = ϵ n A -ϵ n-1 A ∆t =⇒ ϵ n Ae = (ϵ Ae ) t r 1 + ∆t τ , ( . )
which is the equation . . In the above equation, (ϵ

Ae ) t r = ϵ n A -ϵ n-1 A = ϵ n A - ϵ n-1 A -ϵ n-1
Ae . e FLV model in equation A. can be wri en in the case of small strains and small perturbations by realizing that

C s 33 = 1 det C s 2d ≈ 1, ∀ 0 ≤ s ≤ t, ( . ) (C s 2d ) -1 ≈ I -2ϵ s , ∀ 0 ≤ s ≤ t, ( . ) tr( Ŝs C s 2d ) = tr(τ s ), ∀ 0 ≤ s ≤ t, ( .
)

Ŝs ≈ τ s ∀ 0 ≤ s ≤ t . ( . )
Using all the above, equation A. can be wri en as

τ A = ∫ t -∞ (t -s) ∂τ s A ∂s ds = τ ∞ A + τ Ae , ( . )
where τ ∞ A is the stress in the elastic branch and τ Ae = 2C 10 (4ϵ Ae + 2ϵ Be ) is the stress in the viscous branch, ϵ Ae = ∫ t -∞ exp(-t -s τ ) ∂ϵ A ∂s ds. is, in turn, can be seen to be the analytical solution of equation . . Hence, both FLV and FV models can be seen to give the exactly same result in the limit ∆t τ → 0. However, when ∆t τ does not tend to 0, the two models can result in slightly di erent result as the FLV model uses the semi-group property of the exponent to integrate the equations while FV model uses backward Euler. e predictions of FLV and FV models are examined for the geometry in gure . . e top portion of the model is xed in the x-direction, u x = 0 and is displaced in y-direction as follows. It is subjected to a displacement of 5 mm in 5 × 10 6 sec to make sure that the viscoelastic e ects do not come into the picture. It is then subjected to a cyclic displacement A + B sin(2πωt), where A = 5 mm, B = 2 mm, and ω = 1000 s -1 . e constants 1 and 2 for FLV have been taken to be 0.3 and 0.4, respectively, while τ 1 and τ 2 to be 1 × 10 3 s and 1 × 10 4 s, respectively. e sti ness ratios of FV model are adjusted accordingly using the equation . and the η D s are computed as 2C 10 τ . Time step size has been set to 1 × 10 -4 s. Hence, ∆t τ is about 10 -7 . e evolution of σ 22 and ϵ 22 with time for the two models along with the di erence between them at the centroid of the element shaded in brown in the gure . can be seen in the gure . . As can be seen, the two models predict the same result within plo ing accuracy. e di erence between the stresses predicted by them lies within 3% of the maximum value and the di erence between the strains lies within 0.2% of the maximum value.

is demonstrates the accuracy of implementation of the FV model into the UMAT subroutine.

. . Viscoelastic case -large strains and small perturbations

In this case, the approximations of F e and C e in the previous section still hold. However, the background strains are taken to be nite and so, they cannot be approximated using small strain measures. In this case, the evolution equations for FV model for viscous branches remain similar to the equation . , except that the total strains, ϵ A , are now logarithmic, ϵ = 1 2 ln C. e evolution equations in this case can be shown to be exactly that of Lubliner, . e model in the gure . is again used to test this scenario. e top edge of the model is xed in the x-direction, u x = 0 and is displaced rst vertically by 40 mm in a span of 5 × 10 7 sec. It is then subjected to a cyclic displacement A + B sin(2πωt), where A = 40 mm, B = 2 mm, and ω = 1000 s -1 . e model coe cients are as in the previous section and the time step has been set to 1 × 10 -4 s as well. It can be seen from the gure . that the two models predict slightly di erent values of stresses at the same strain level. e FV model predicts a slightly smaller stress level than the FLV model.

e di erence between the FLV and the FV models in this case can be seen by writing the FLV model for the case of small perturbations superposed on large background deformations. For this case, FLV model becomes

S 2d ≈ Ŝ∞ 2d + Ŝ 2d . ( . )
In the above, Ŝ 2d ≈ ∫ t -∞ (ts)

∂ Ŝs ∂s ds. It has been assumed for times t > t 1 , where t 1 is the time when the specimen is stretched to the target pre-stretch level (in this case 5 × 10 7 s), F (t) can be wri en as F (t) = F (t 1 )F r el , where F r el is the relative deformation gradient. Since the perturbations imposed are in nitesimal, F r el has been assumed to be of the form I + ϵ, where the magnitude of ϵ is much smaller than 1. Hence, the higher order terms in ϵ have been neglected.

Ŝs is computed according to the equation A. at time s. e evolution equations for the FV model remains same as that in the equation . . Hence, it can be seen that the evolution equations di er for the FV and the FLV models in this case.

. To compare the predictions of FLV and FV models when applied to large strains and large perturbations, the methodology employed in the chapter will be used, where the crack speed is imposed on the model using the boundary conditions from the experiment. e FLV model has been used in that chapter and was 'calibrated' to match the observations from the experiments. Here, the same sti ness ratios and the relaxation times will be used for the FV model and the results will be compared to the predictions from FLV model. e sti ness ratios in the chapter are with respect to the glassy modulus of the material, while the sti ness ratios in the FV model are with respect to the rubbery modulus. Hence, they are scaled using the expression

( i ) rubbery = ( i ) glassy 1 -N i=1 ( i ) glassy . ( . )
e parameter γ 0 in the equation . has been taken to be 1. e above model parameters were then used to perform some simulations to understand and compare the behaviors of FLV and FV models.

Table . : FV model parameters with γ 0 = 1.

Branch(i) Sti ness ratio, ( i ) rubbery . . . . . . . . Relaxation time, τ i (s) E- E- E- E- E- E- E- E-
e scaled parameters and the relaxation times can be seen in the table . . e FV model with these parameters are then used to determine the displacement and velocity elds for the case of λ = 3.5, where the crack speed is about 56 m s -1 . e predictions from the model can be seen in the gures . and . (to be compared with gures . and . a, respectively).

e velocity pro le can be seen to be sharper than that from the FLV model but blunter than in the hyperelastic case. Also, the horizontal displacements can be seen to be under predicted when compared to the experiments and the FLV model. However, the pro le can be seen to have some resemblance with the hyperelastic case (see gure . ). -14.6 -11.71 -8.83 -5.95 -3.07 -0.18 2.7

As in the previous chapters, the variation of particle speeds at about 5 mm above the crack path has been extracted from the FE simulations and the result can be seen together with the results from experiments, hyperelastic FE simulations, and FLV simulations in the gure . . It can be observed from the gures that the predictions from the FV model do not match the experiments or the results from the simulations with FLV model. In fact, the velocity magnitude can be seen to closely follow the hyperelastic case without any jumps.

e elds can be seen to be continuous as opposed to the hyperelastic case.

Further discussions

It can be seen from the above that the results from the FV model with γ 0 = 1 are closer to the hyperelastic case, but without a shock front. In order to determine if the viscous e ects are active in these simulations and the extent to which they are active, a viscous indicator parameter β is de ned as follows.

β ψ ∞ (b) + N i=1 ψ i (b (i) e ) ψ ∞ (b) . ( . )
N denotes the number of viscous arms. e ψ i s in the above equation indicate the strain energy densities in each of the viscous arms. It has been indicated that the ψ i s depend on the corresponding b (i) e s, the le Cauchy-Green tensor based on the corresponding elastic parts of the deformation gradient, F (i) e = F F (i) -1 s. Hence, the further the value of β from 1, the further is the material behavior from the hyperelastic case. A value of 1 for β indicates a purely hyperelastic behavior. For the result presented above for the case of λ = 3.5, the value of β in various regions of the body can be seen in the gure . . e scale has been truncated at 1.3 to bring out the features of the distribution. It can be seen that the value of β remains very close to 1 (about 1.0008) in front of the front-like region to the right of the crack. Its value then increases to about 1.2 and then behind another front-like feature, it again increases beyond 1.3. It shall be noted that the two front-like features in the gure . coincide with the shock front-like features observed for the particle velocities in the gure . b. Hence, in the case of the simulations with FV model with the current parameters, it can be said that the viscous e ects di use or blunten the shock front-like feature observed in the hyperelastic simulations (the velocity variations in gure . are closer to the hyperelastic simulations, but without the jumps) by increasing the sti ness of the material as a consequence of higher strain rates (like discussed in the chapter and in Kamasamudram et al., ). e FV model is expected to displace the relaxation times to smaller values at higher strains. Hence, it can be is is consistent with the predictions from Govindjee and Reese, , where the relaxation functions of FLV and FV models are compared for di erent strain levels by subjecting a square object with a side of 1 mm to a certain displacement under plane strain conditions and observing the evolution of stress with time.

e same numerical experiments will be performed for the models in plane stress conditions and the results will be presented in the appendix B. For the case of FLV model, the authors of Govindjee and Reese, comment

As the strain of the loading is increased the curve shi s upward as an indication of the nonlinear sti ening of the material in the current (spatial) con guration. Note that the curves merely shi upwards and do not translate simultaneously along the time axis.

While for the case of FV model, e curves simultaneously map le and upwards as the amount of strain increases. is mapping to the le is a result of the nonlinear structure of the viscoelasticity model itself and not due to the explicit incorporation of any shi -factors.

e shi towards le of the relaxation curves at higher strain levels indicate the shi of relaxation times to smaller values with an increase in strain. is implies that the FLV model behaves in a sti er (or more viscous) way when compared to the FV model at larger strains for a given strain rate.

In addition to the above, the activity of i th viscous arm in di erent regions of the body can be examined using the ratio ψ i ψ ∞ , the ratio of strain energy density in the viscous arm to that of the energy in the hyperelastic arm. e ratios for the case of λ = 3.5 can be seen in the gure . . Activity can be seen along the shock-front-like features in multiple branches.

. . Model with adjusted parameters

As was observed in the simulations with the imposed crack speed, the predictions of the FLV and the FV model di er when applied for the stretch ratio of 3.5. is is the consequence of the shi in the relaxation times with the applied strain in the la er model. To counter this, the parameter γ 0 in the equation . can be made less than 1. It shall be noted that the reducing γ 0 from 1 has the same e ect of increasing the relaxation times (moving the relaxation curves to the right in the sense of Govindjee and Reese, ). Hence, the model becomes sti er than when γ 0 = 1 at a given strain rate. A er performing some simulations with various values of γ 0 < 1, it was found that choosing γ 0 as 1 60 , while modifying one of the a:

ψ 1 ψ ∞ 0 1.
ψ i
∞ for all the viscous branches (i = 1 to 8) for λ = 3.5 plo ed on deformed con guration, γ 0 = 1. e position of crack tip is indicated by a circle. sti ness ratios (corresponding to the relaxation time of 1 × 10 -7 s) results in the predictions of the FV model matching the results from the experiment. e new model parameters can be found in the table . . Only 7 viscous arms were found ψ ∞ , any branches with values smaller than 10 -5 have been discarded) for the simulations with the new model parameters.

e simulations from the previous section are repeated with the new parameters and the predictions of the model can be seen in the gures . and . . -14.62 -11.48 -8.34 -5.2 -2.07 1.07 4.21 e velocity pro le in the gure . can be found to be blunter than in the case of γ 0 = 1 ( gure . ) and closer to the FLV ( gure . ). e horizontal displacement can also be seen to be closer to the experimental and the FLV results. However, FV model with γ 0 = 1/60 slightly under predicts the displacements. e variation of particle velocity magnitude at about 5 mm above the crack path can be seen in the gure . together with the results from the experiments, FLV model, and FV model with γ 0 = 1. It can be seen that the predictions of FV model with the new model parameters are close to that of the FLV model and the experiments.

e variation of parameter β that determines how much the viscoelastic e ects are active can be seen in the gure . . e scale has been truncated so that it can be compared to the results from the old model parameters ( gure . ). e variation of β at about 5 mm above the crack path for both the cases can be seen in the gure . . It can be seen that the viscous e ects start to act further ahead in the material with γ 0 = 1/60 than with γ 0 = 1. is is expected, since a value of less than 1 for γ 0 moves the relaxation times to larger values which would result in the viscous e ects starting to act at smaller strain rates. e increase in β in the gure . is also more gradual. Also, the value of β exceeds 1.3 in the region just behind the tip than when compared to the previous case indicating more prevalent viscous e ects in that region. e resulting horizontal displacements can be seen in the gure . . e displacements from the experiment for λ = 2.5 can be seen in gure . . Likewise, the velocity evolution at about 5 mm above the crack path, together with the results from the experiments and the FLV, can be seen in the gure . . -10.56 -8.28 -6 -3.71 -1.42 0.86 3.14 -10.56 -8.28 -6 -3.71 -1.42 0.86 3.14

Further discussions e FV model with the parameter γ 0 = 1/60 and one adjusted sti ness ratio (to be called shi ed FV model from hereon) can be seen to predict the velocity and displacement elds close to the experiments when compared to γ 0 = 1. Choosing the parameter γ 0 less than 1 has the same e ect of moving the relaxation times to higher values. In short, choosing γ 0 = 1/60 has the same e ect as retaining γ 0 = 1 and multiplying each of the relaxation times by 60. One of the consequences of γ 0 1 is that the FV model does not coincide with that of ( nite) linear viscoelastic model at small strains and small perturbations. e relaxation times of the la er model are to be modi ed to make their predictions match at these conditions.

e shi ed FV model can be seen to slightly under predict the displacements at λ = 3.5 and slightly over predict them at λ = 2.5.

is indicates that the viscoelastic model tends to get sti er as the strains are reduced, as expected. However, the shi ed FV model can be seen to predict the experimental results with a 'good enough' accuracy. Hence, energetic analysis of those experiments will be performed using this model in the next section.

. Since the shi ed FV model was able to predict the displacements and velocities from the experiments of λ = 2.5 and 3.5 with a good enough accuracy, it will be used to perform an energetic analysis of those experiments. e FLV model from the previous chapter could not be used to do so as a consequence of lack of an expression for the strain energy function and the energy dissipation in that model.

. . From the experiments

First, the energy balance will be performed for the case where the crack speed is imposed onto the FE model using the displacements from the experiments. It shall be noted that the analysis can also be performed directly on the displacements and velocities as done in Corre, . e methodology is exactly the same as in the section . , but with the FLV model in the bulk replaced by the shi ed FV model.

For the cases of λ = 2.5 and 3.5, the evolution of strain energy, kinetic energy, and the dissipation in the material with time have been extracted for the simulations presented in the previous section with the shi ed FV model. eir evolution with time can be seen in the gures . and . once the initial transients have passed. It shall be noted that only half of the specimen was used to obtain these results taking symmetry into account. Using this data, energy balance can be wri en as

P ex t = dSE dt + dD dt + dKE dt , ( . )
where P ex t is the power of the external forces, SE, D, and KE denote the strain energy, visoelastic energy dissipation, and kinetic energy, respectively. In the current case, since only the part of material that is about 1 mm away from the crack path has been modeled, P ex t includes the power of external forces, the viscoelastic energy dissipation in the material that is not included, and the power dissipated during the fracture process. e time derivative of the quantities for the case of λ = 2.5 and 3.5 can be seen in the table . . It shall be noted that the sign of P ex t is negative in both the cases. is indicates that the work has been extracted from the system or work has been done by the system during the crack propagation. As already mentioned, this work includes the work in the fracture process zone and the dissipation in the dt by the crack speed gives Energy release rate, which is experimentally computed (usually) as ψh 0 , where ψ is the strain energy density is the material far ahead of the crack tip. h 0 is the undeformed specimen height. In the current case, the strain energy behind the tip may not be zero as a consequence of viscous e ects. Hence, the quantity ψ h 0 = (ψ (∞)ψ (-∞)) h 0 (see gure . ) will have a di erent value from ψh 0 . In the current case, ψ (∞) and ψ (-∞) have been computed at about 50 mm in front of and behind the tip, respectively. All the computed quantities can be seen in the table . . e factor 2 above has been used since only half of the specimen geometry has been used in the analysis. From the tables . and . , and the gure . , it can be observed that the majority (about half) of the strain energy is taken up by the dissipating material. e kinetic energy takes up very li le energy when compared with the viscoelastic dissipation. e rest (P ex t ) is taken up by the material along and a li le away from the crack path in the form of viscoelastic dissipation and the fracture processes at the crack tip. e quantities -2 w 0 dSE dt and ψ h 0 are quite close to each other for both the cases presented. Again, it shall be noted that the strain energy behind the tip used in the computation of ψ is the total strain energy contained in the elastic as well as the viscous arms. When the above results are compared with that of the analysis performed in Corre, (see gures . and . of the cited reference), some di erences can be noticed. e analysis in the cited reference uses a hyperelastic material to compute the evolution of the total energy. e energy release rate (obtained by evaluating d(SE + KE) dt , where the strain energy has only the elastic contribution) from that analysis can be seen to match the quantity ψh 0 in the same analysis as well as from the table . . However, when viscous e ects are included, some strain energy is locked in the viscous arms, which a ects the computation of d(SE + KE) dt . is quantity is hence smaller in the current study and it can be seen to match ψ h 0 , rather than ψh 0 .

. .

From the simulations with cohesive model e previous sections analyze the experiments performed in Corre et al., . However, as can be noted, only the material starting from about 1.5 mm above the crack path has been modeled. Hence, the term P ex t included the viscous dissipation in the excluded material. In the simulations performed using the cohesive model using the FLV model presented in the previous chapter, this problem does not exist. Data can be extracted right along the crack path from those simulations which then can be imposed onto the current model, where the bulk material has been modeled by the FV model. P ex t computed this way will solely comprise of the work done in the cohesive zone and the dissipation, D, will include the total viscoelastic dissipation in the material. It shall be noted that the crack propagation has not been modeled explicitly in this section. Only half of the bulk material has been used with the crack speeds imposed implicitly . is procedure has been followed for the case of λ = 2.5 and 3.5. e results of energy evolution for both the cases can be seen in the gures . and . . When comparing the results from the experiments, it can be seen that the strain energy release rate is about the same in both the cases. Viscoelastic dissipation is higher in the results from the simulations.

is is expected, as the data in the case of FE simulations have been extracted right along the crack path and hence includes the material where higher dissipation takes place. e P ex t simply becomes P coh , the energy spent during the fracture processes. Even in this case, the two parameters, -2 w 0 dSE dt and ψ h 0 can be seen to be close. e computed value of G coh can be compared with the area under the traction separation curves ( gure . ) in the previous chapter. e areas under the traction separation curves can be approximately seen to be 18.5 kJ/m 2 and 21 kJ/m 2 , for λ = 2.5 and 3.5, respectively. e results from the shi ed FV model can be seen to over predict these quantities (as 24 kJ/m 2 and 30 kJ/m 2 , respectively). Performing a FE simulation with the FV model in bulk together with the CZM may result in a more accurate estimate of this energy. However, that analysis has not been done as a result of higher computational cost involved.

e region of the body in which the viscoelastic dissipation takes place can be seen in the gure . . Signi cant dissipation can be seen in the material that is within about 2 mm from the crack path.

is can explain the di erences in the viscoelastic dissipation observed from the results in tables . and . . Also, the temperature measurements made on the test samples from batch-revealed a signi cant increase in temperatures near the tip. us, the viscous dissipation region in the gure . can be seen to be a qualitative match to the experiment.

.

FV vs FLV models

In this chapter, the plane stress version of the Finite Viscoelastic model of Reese and Govindjee, has been used to understand the energy budget during the viscoelastodynamic fracture. e model has been implemented into UMAT subroutine of Abaqus (Dassault,

). Initially, the analyses were performed using the same sti ness ratios and the relaxation times as that of the FLV model. It was observed that the velocity and displacement predictions from FV model did not match the experimental observations. e velocity raise was seen to be sharper than in the FLV case and the displacements were seen to be under predicted and closer to the hyperelastic case. ere was, however, no shock-front-like feature in the results from FV model unlike for the hyperelastic case. Using the viscous indicator, β, which is de ned as the ratio of the total strain energy to that in the hyperelastic branch, the regions in the body where the viscoelastic e ects are active have been identi ed. It was observed that the value of β is very close to 1 in front of a front-like region. It value increased to about 1.3 in the region behind another front like region. e two front-like features observed in the gure . coincide with the shock-front-like features observed in the velocity pro les in the hyperelastic case. is, when combined with the velocity pro les, indicate that the e ect of viscoelasticity is to blunten the shock-front-like feature in the velocity distributions observed in the hyperelastic case. e di erences from FLV case can be a ributed to the relaxation times shi ing to the smaller values with an increase in the strain in the FV model.

e results from the FV model were then brought closer to the experiments by decreasing the value of the coe cient γ 0 (= 1 initially) to 1 60 . is can be viewed as multiplying the relaxation times by 60 to compensate for their decrease at higher strains. is, whoever, will lead to a di erence in the predictions between the FV and the FLV models at small strains (FV model with γ 0 = 1/60 will behave sti er than the FLV model at small strains with the same sti ness ratios and relaxation times). But decreasing γ 0 was seen to bring the results from simulations closer to the experiments. e velocity and displacement elds for λ = 3.5 and 2.5 were seen to be closer to the experiments than when γ 0 = 1. Even with the reduced γ 0 , the FV model can be seen to slightly under predict the displacements at λ = 3.5 and over predict them at λ = 2.5. By making the model agree with the experiments at λ = 3.5, it was made to slightly sti er than needed at 2.5.

One way to compensate this would be to make γ 0 a decreasing function of either b e = F e F T e , corresponding to the elastic part of F or C, the overall deformation tensor, such that γ 0 (b e or C = I ) = 1. is is equivalent to slightly compensating for the shi in the relaxation times as a consequence of the non-linearity in the FV model. is will result in the model predicting the same results as that of FLV model at small strains and agreeing with the experimental results at higher strains. However, when γ 0 is made a function of C (its invariants, more precisely, as a consequence of isotropy), the viscous evolution equations then become

ϵ 1e + γ 0 (C)∆t 2η D (τ 1 + p) -(ϵ 1e ) t r = 0, ( . )
is poses no problem when solving for the ϵ Ae s in the local Newton iterations as C does not change during those iterations and hence γ 0 can be treated as a constant. However, when the global tangent is to be computed as a part of global Newton iterations, γ 0 cannot be treated as a constant. e contribution of the viscous evolution equations to the global tangent comes in the form of the derivative ∂ϵ Ae ∂(ϵ Be ) t r

. Since the overall deformation changes during the global iterations, γ 0 (C) contributes to the tangent as well. An a empt can be made to evaluate the derivative by observing that . However, this requires the knowledge of F n-1 , which is not known. Only C n-1 is known and hence F n-1 can be evaluated up to rigid body motions. One other alternate is to make γ 0 a function of the invariants of b. In this case, b n = F n (F n ) T = (F e ) t r (b n-1 ) T ((F e ) t r ) T . Again, this requires the knowledge of (F e ) t r , which can be determined only up to a rigid body motion. Hence, this path was not explored further in this thesis. e other alternate is to make γ 0 a function of b e like in Bergström and Boyce, . However, the Bergstorm-Boyce model may lead to a faster decay of stresses in the relaxation tests than that of the Reese and Govindjee, model, which is equivalent to shi in the relaxation times further to the le .

C n = (F n ) T F n = (F n-1 ) T (C e ) t r F n-1 . Since F n-
e cohesive zone model used in the previous chapter is based on the use of convolution integrals, where the relaxation times and the sti ness have been chosen to be same as that of the bulk material. However, in the current chapter, it can be seen that the inclusion of non-linearity in the viscoelastic model shi s the relaxation times of the material with the applied strain. How this shi a ects the relaxation times in the cohesive zone is currently not known. However, this shi can be included in the cohesive zone by using 'strain clock' models (Wineman,

), where t -s in the equation . is replaced by ξ (t)ξ (t), where ξ is a function of strain and possibly its history.

Energy budget in viscoelastodynamic fracture

e shi ed FV model can be seen to predict the displacements and velocities close to the experiments for λ = 2.5 and 3.5. Hence, it has been used to analyze the energy evolution for those experiments. e strain energy release rate was seen to match the quantity ψ h 0 and smaller than ψh 0 . When compared with the computations from Corre, , the value of -2 w 0 dSE dt in the current case can be seen to be smaller than what can be seen in the gure . of that study. is can be a ributed to the viscoelastic e ects included in the current study. Using a similar method, the results from the cohesive zone simulations performed in the previous chapter have been analyzed as well. A signi cant portion of the strain energy can be seen to be dissipated by the viscoelastic e ects in the material. e rest of the energy can be seen to be taken up by the fracture processes. e shi ed FV model over predicts the energy dissipated by the fracture processes when compared with the results from the cohesive zone simulations from the previous chapter. Perhaps performing the simulations with the cohesive zone and the shi ed FV model in the bulk will make in both the results to agree.

e region in the body where the viscoelastic dissipation takes place can be seen in the gure . . e region near the crack path can be seen to dissipate the majority of the energy, while the dissipation in the material that is a bit far can be seen to be smaller. An element along the crack path and in the vicinity of the tip undergoes very high rates of deformation, which can be seen to result in higher dissipation. is can be seen as a reason for the di erence in the viscoelastic dissipation computed in tables . and . . Also, the temperature measurements from the experiments performed on the materials of di erent batch reveal a large increase in temperature along the crack faces (see gure . ). is is consistent with the regions where the viscous dissipation is concentrated in gure . . A quantitative comparison of the temperature raise has not been made as the samples are from di erent batches.

e plane stress version of the FV model has been presented and implemented in Abaqus. Using the relaxation times that are di erent from the FLV model, it was observed that the FV model was able to predict the results from the experiments reasonably well. e energy evolution during the crack propagation in the experiments from Corre, and the simulations with the cohesive zone from the previous chapter has been computed. e slope of the strain energy evolution vs time curve has been observed to match ψ h 0 , which stresses the importance of including the viscous e ects in the analysis. Examining the parameter β, the viscous e ects were seen to be active starting from the region just in front of the crack tip. e region along the crack faces behind the tip can be seen to dissipate considerable amount of total energy. is region is comparable to the region where there was a considerable temperature raise in the experiments.

Conclusions and perspectives

e analysis of dynamic fracture of elastomers has been made in the current thesis. First, regarding the origin of transonic cracks in elastomers, examining the velocity and strain elds from the experiments revealed no presence of shock fronts in the material. Simulations have been performed to explain this observation using the experimental data and implicitly imposing the cracks speeds. In one case, the bulk material has been assumed to be hyperelastic that exhibits an upturn in its stress strain response.

e velocity elds obtained from this case were observed to contain jumps along a shock-front-like feature, which is absent in the experiments. When the bulk has been made viscoelastic, the shock-front-like features disappeared and the displacement, velocity and the strain elds resembled those from the experiments. e hyperelastic branch of the material still included the upturn in its stress-strain response. When it is replaced by an model that does not exhibit the upturn, the elds were still observed to be continuous, devoid of any shock-front-like features. Hence, it has been concluded that the observation of transonic cracks in polyurethane elastomers can be a ributed to the viscoelastic sti ening of the material rather than the hyperelastic sti ening.

Experiments performed on natural rubber specimens of di erent geometries in literature revealed that the crack speeds were independent of the specimen geometry (speci cally the height) beyond a certain threshold. Whether such a conclusion is valid for all the materials is not known. To address this, the current author performed more experiments on Polyurethane samples of di erent geometries. However, as a consequence of either defective material or the material with di erent properties from di erent batches, no conclusions could be drawn from the experiments performed along with the experiments from the literature. But, the crack speeds from batch-for di erent geometries seemed to se le down at a certain level starting from some stretch level, indicating that they became height independent. e cohesive zone model has been used to predict the crack speeds observed in the experiments since the crack path is already known in advance if initiated at the middle of the specimen. Simulations have been performed with cohesive zone inserted all along the prospective crack path and the bulk material modeled by the Finite Linear Viscoelastic model. Rate dependence was found to be necessary in the cohesive zone to match the crack speeds observed in the simulations with that of the experiments.

e rate dependence has been introduced in terms of convolution integrals with the relaxation times and sti ness ratios chosen to be same as that of the bulk material. Using this model, the crack speeds observed in the experiments at di erent stretches could be reproduced. e tractions in the cohesive zone and the energy expended in the cohesive zone have been found to be an increasing function of crack speed for the cases studied. Also, when the specimens of di erent heights were simulated, it has been observed that the crack speeds became independent of the height from a certain threshold like in the experiments from the literature.

In the next part, a Finite Viscoelastic model has been implemented in the plane stress se ing. Some di erences have been observed in the predictions of the FV and the FLV models as a consequence of the non-linearity of the evolution equations in the former model. From the numerical relaxation tests performed, it has been observed that the relaxation times in the FV model are a function of strain as a consequence of the non-linearity of the model. For the FLV model, the stress vs time curves during the numerical relaxation tests at di erent strain levels were observed to just translate (upward along the stress axis) without any change in shape, indicating an independence of relaxation times on the strain level. e material parameters have been modi ed in the FV model so that the predictions of this model matches the observations from the experiments. Using this FV model, the strain energy evolution and the dissipation in the material during the propagation of crack has been studied. It has been found that the majority of the energy is consumed as the viscoelastic dissipation of the bulk material. A comparatively small portion is consumed in the cohesive zone. is has been the case for both the cases studied in this thesis. Whether the energy consumed in the cohesive zone can be computed directly from the experiments in any other way is not known at the moment.

From the temperature measurements made from the experiments, a large raise in temperature has been observed in the material in the vicinity of the crack tip. is region also corresponds with the regions where the dissipation takes place, as identi ed by the FV model.

e viscoelastic properties of the material are expected to be a function of temperature, which in turn a ects the dissipation and possibly the crack speeds as well. Solving a fully coupled thermoviscoelastic problem will help gain insights into how these phenomena interact with each other. Also, the FV model used in this thesis has been calibrated using the data from the crack propagation experiments. As expected, some di erence has been noticed between the predictions of the FLV and the FV models. However, whether the actual behavior of the material is closer to one model or the other is not known. Perhaps, performing simple relaxation tests like in the appendix B can help understand the behavior of the material. Once the model is calibrated to predict these experiments accurately, it can perhaps be further re ned using the data from the crack propagation experiments.

e other aspect would be to consider the shi in relaxation times while modeling the behavior of the failing material. In the convolution integral type models used in this study, this can perhaps be achieved by using strain clock models, where the times in the relaxation functions are replaced by shi ed times depending on the strain level. How this will a ect the crack speeds and other observations from the simulations is not known.

One other aspect of interest would be to relate the surface roughness variations observed with the crack speed and the viscoelastic properties of the material. e role played by cavitation and the e ects of viscoelasticity on cavitation has to be investigated. e e ect of specimen thickness on the crack speeds and roughness in the high speed regime needs to be investigated as well. ). Peel process simulation of sealed polymeric lm computational modelling of experimental results.

Engineering Computations (Swansea, Wales), ( ), -. doi: . / ). Crack propagation in a linearly viscoelastic strip.

Journal of Applied Mechanics, , -. doi: . / . Knowles, J. K. (

). e nite anti-plane shear eld near the tip of a crack for a class of incompressible elastic solids. International Journal of Fracture, ( ), -. doi: . /BF Kostrov, B. V. & Nikitin, L. V. (

). Some General problems Of Mechanics Of Bri le Fracture. Archiwum Mechaniki Stosowaniki, ( ). Krishnan, V. R., Hui, C. Y., & Long, R. (

). Finite strain crack tip elds in so incompressible elastic solids. Langmuir, ( ), -. doi: . / la e Kroon, M. (

). Energy release rates in rubber during dynamic crack propagation. ). e elasticity of a network of long-chain molecules-II.

International

Transactions of the Faraday Society, ( ), -. doi: . /TF Tsunoda, K., Bus eld, J. J., Davies, C. K., & omas, A. G. (

). E ect of materials variables on the tear behaviour of a non-crystallizing elastomer. ). Finite deformation constitutive equations and a time integration procedure for isotropic, hyperelastic-viscoplastic solids. where the rst DEV t • projector is evaluated at time t and the projector inside the integral is evaluated at time s. In this appendix, the plane stress version of the above model will be derived. First, the above equation is simpli ed to move the time derivative from the DEV • to the rst term by using integration by parts. (A. )

As a consequence of the plane stress assumption, the deformation gradient can be wri en in a block form to be . ) e above form is assumed to be valid at all times. Now, the expression for S, wri en in terms of D-strain invariants can be seen to be . ) where ) Check e correctness of the above equation can be checked by se ing (t) = ∞ , which would be the case of a hyperelastic material. . ) which is the result expected for a hyperelastic material.

F = F 2d o o T F 33 . ( A
S = 2 ∂ψ ∂I 1 I + 2 ∂ψ ∂I 2 [I 1 I -C] , ( A
I
S 2d (t )=1 = ∞ DEV 2D t Ŝt + ∞ 3 tr( Ŝt C t 2d )(C t 2d ) -1 = ∞ Ŝt , ( A
A.

e expression in equation A. can be integrated in time by using the semi-analytic technique in Taylor et al.,

. Denoting the term inside the time derivative as h and denoting the integral as H , the time integral can be wri en as relaxation times to smaller values at larger strains in the case of FV model. Also, the shapes of the curves at various stretch levels can be seen to be similar in the case of FLV model. e relaxation times can be seen to be not a ected in this case. is observation is in coherence with what was presented in Govindjee and Reese, (see gure . of the reference).

H (t) = ∫ t -∞ (t -

B. . Analysis

e results from the relaxation test above can be analyzed as follows. In the case of FLV model, the stress can be wri en as e loading can be divided into two steps, where in the rst step, the body is rapidly loaded to the target strain level and in the second step, it is held at that strain level.

S 2d = ∫ t -∞ (t-
e strain (C tensor) is assumed to be increased from I to C t 1 in a span t 1 that is much smaller than the corresponding relaxation time, τ , in the rst step. Assuming only one viscous branch for convenience, this step corresponds to the glassy material behavior. is can be seen from the fact that t -s τ → 0 ∀ t, s ≤ t 1 τ and so, the exponential term will almost be equal to 1. e integrals can hence be evaluated directly. e total stress at the end of rst step can be wri en as

S 2d = Ŝ∞ 2d | t =t 1 + 1 DEV 2D t 1 Ŝt 1 + 1 3 tr( Ŝt 1 C t 1 2d )(C t 1 2d ) -1 = Ŝ2d | t =t 1 , (B. )
since ∞ + 1 = 1. Hence, it can be wri en that

∫ t 1 -∞ 1 exp(-(t -s)/τ ) ∂ ∂s DEV 2D s Ŝs = 1 DEV 2D t 1
Ŝt 1 , and (B. )

C t 1 33 3 (C t 1 2d ) -1 ∫ t 1 -∞ 1 exp(-(t -s)/τ ) ∂ ∂s tr( Ŝs C s 2d )(C s 33 ) -1 = 1 3 tr( Ŝt 1 C t 1 2d )(C t 1 2d ) -1 . (B. )
In the second step, the strain is held constant and the evolution of total stress is studied in the intervals of ∆t. For this step, the two integrals in the stress term . r e a l ( ) ( s e q ( , ) , ) ; v = e s . e i g e n v e c t o r s ( ) . c o l ( ) . r e a l ( ) ( s e q ( , ) , ) ; v = e s . e i g e n v e c t o r s ( ) . row ( ) . r e a l ( ) ( , s e q ( , ) ) ; v = e s . e i g e n v e c t o r s ( ) . row ( ) . r e a l ( ) ( , s e q ( , ) ) ; e i g s t r ( ) = e s . e i g e n v a l u e s ( ) ( , ) . r e a l ( ) ; e i g s t r ( ) = e s . e i g e n v a l u e s ( ) ( , ) . S t i f f p r i n c i p a l ( , ) ) ; p r e s s u r e t a n g e n t ( , ) = -. / . * ( S t i f f p r i n c i p a l ( , ) + S t i f f p r i n c i p a l ( , ) ) ; p r e s s u r e t a n g e n t ( , ) = p r e s s u r e t a n g e n t ( , ) ; p r e s s u r e t a n g e n t ( , ) = p r e s s u r e t a n g e n t ( , ) ; C a l g = S t i f f p r i n c i p a l * ( T a n g e n t p r i n c i p a l . i n v e r s e ( ) ) ; m a t t a n ( s e q ( , ) , s e q ( , ) ) . a r r a y ( ) = C a l g -. * M a t r i x d ( t a u p r i n c i p a l . a s D i a g o n a l ( ) ) ; i f ( e i g s t r ( ) == e i g s t r ( ) ) { m a t t a n ( , ) = . * ( C a l g ( , ) -C a l g ( , ) ) -t a u p r i n c i p a l ( ) ; } e l s e { m a t t a n ( , ) = ( t a u p r i n c i p a l ( ) * e i g s t r ( ) -t a u p r i n c i p a l ( ) *

.

  Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Crack Propagation condition . . . . . . . . . . . . . . . . . . . . . . . Energy ux integral . . . . . . . . . . . . . . . . . . . . . . . Propagation direction . . . . . . . . . . . . . . . . . . . . . Weakly Nonlinear eory of fracture . . . . . . . . . . . . . . . . . Crack tip elds in elastomers . . . . . . . . . . . . . . . . . . . . . . Transonic cracks . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure

  Figure . : A body with a crack.

Figure . :

 . Figure . : Polar coordinates centered at the crack tip.

  Figure . : Opening mode (I), Shear mode (II) and Anti-plane shear mode (III).

Figure . :

 . Figure . : Contour Γ surrounding the tip.

Figure . :

 . Figure . : Rectangular contour Γ surrounding the tip. It is shrunk to the tip rst by le ing δ 2 → 0 and then δ 1 → 0.

Figure . :

 . Figure . : Polar coordinates centered at the crack tip in undeformed and deformed con gurations.

Figure . :

 . Figure . : Comparison of experimental results (A) with the theory from (Freund, ) with a sharp crack assumption (B) and the theory with anite failure zone (C). Reproduced from (Rosakis et al.,).

Figure . :

 . Figure . : SIF vs crack speed. Comparison of results from Dally, and Ravi-Chandar and Knauss, b. Adapted from Ravi-Chandar and Knauss, b.

Figure . :

 . Figure . : Crack speeds obtained when using a weak interface and holes (Washabaugh & Knauss,).

Figure . :

 : Figure . : Microbranching instability in PMMA (Sharon et al.,). Oscillations in crack speed a er a certain limit can be seen.

Figure

  Figure . : Oscillatory instability observed in Livne, Ben-David, and Fineberg, . A sequence of photographs, shown at . ms intervals, of a propagating crack undergoing a transition from linear (top pictures) to oscillatory motion. Photographs of XY pro le (top) and (XZ) fracture surface (bo om) of (b) a . mm thick gel sample where oscillations developed and (c) a . mm thick gel where the crack retained its straight line trajectory.

Figure . :

 . Figure . : Tearing con guration.

Figure . :

 . Figure . : Single Edge Notched Tension specimen con guration.

Figure . :

 . Figure . : Crack propagation data for multiple production runs, prior to applying the time-temperature superposition process (W. Knauss,).

Figure . :

 . Figure . : Description of cohesive zone from W. Knauss, .

Figure . :

 . Figure . : Undulations on fracture surfaces reported in Stevenson and omas, . λ = 2.3, 3.7, 5.0 in top, middle and bo om gures.

Figure . :

 . Figure . : Crack propagation studied in pure shear con guration in Gent and Marteny,a. Also can be seen is the result of the crack speed vs the specimen height for the same stretch level.

Figure . :

 . Figure . : Reference con guration, B ξ ( xed), Material con guration, B X (t) and Spatial con guration, B x (t)

Figure . :

 . Figure . : Contour on which the integration is performed in the material con guration.

Figure . :

 . Figure . : Contour enclosing the contour at the tip.

Figure . :

 . Figure . : Crack propagation in an in nite strip.

Figure . :

 . Figure . : Typical stress-strain response of an elastomer, Polyurethane in particular. Extracted from Corre, .

Figure . :

 . Figure . : Evolution of dynamic modulus and loss tangent with frequency for polyurethane elastomer. Extracted from Corre, .

Figure . :

 . Figure . : Specimen geometry used in the experiments. Adopted from Corre, .

Figure . :

 . Figure . : Specimen held in tensile machine and razor blade mounted.

Figure

  Figure . : Crack speed vs stretch results obtained in Corre,. Shear wave speed can also be seen as a blue line.

Figure . :

 . Figure . : Crack speeds in deformed (c) and undeformed con gurations (c 0 ). Figure adopted from Corre, with axes relabeled.

Figure . :

 . Figure . : Specimen with a pro le drilled-in. Figure adopted from Corre, with labels in English.

Figure . :

 . Figure . : Results of the experiment with non-homogeneous geometry. (a) Evolution of distance travelled by the crack with time elapsed. (b) Crack speed vs the distance travelled by the crack. Figure adapted from Corre, with axes relabeled in English.

Figure . :

 . Figure . : Propagation of an oblique crack. Figure adopted from Corre, . a crack initiated at the le near the top end. b crack propagates straight before turning. c crack turns and d propagates straight.

Figure . :

 . Figure . : Horizontal and vertical crack speeds for the turning crack. Figure adapted from Corre, with labels renamed in English.

Figure . :

 . Figure . : Results from Chen, Zhang, Niemczura, Ravi-Chandar, and Marder,. In the gure a , crack speeds were plo ed against applied stretch while in the gure b , crack speeds were plo ed against the Energy

Figure

  Figure . : Experimental setup including strain gauges.

Figure . :

 . Figure . : Evolution of reaction forces vs time in the le and the right LVDT.is corresponds to the mm tall specimen and λ = 3.0.

  (a) Cross-head speed of mm min -1 . (b) Cross-head speed of mm min -1 . (c) Cross-head speed of mm min -1 . (d) Cross-head speed of mm min -1 .

Figure . :

 . Figure . : Branching a empts made by crack in the mm specimens (indicated on the broken samples).e arrow points to the location where the branching was a empted. Some branching events were successful while the others were not. An initial seed crack was present in all the cases. e cross-head of the tensile machine was pulled at di erent rates as indicated in the captions.

Figure . :

 . Figure . : Crack branching a empt in mm specimen (batch-) following the normal test protocol. e two side branches were arrested and the center branch propagated further.

Figure

  Figure . : Results of the experiments from batches and . Results from Corre,can also be seen.

Figure . :

 . Figure . : Uniaxial and pure shear response of the material from batch-.Corresponding results from the material from Corre, Coret, Verron, Leblé, and Le Lay, can also be seen.

Figure

  Figure . : Temperature distribution in the vicinity of the tip (in • C). Temperatures as high as • C can be seen (ambient at about • C).

Figure

  Figure . : Horizontal displacement (in m) for λ = 1.7 plo ed on undeformed con guration.

Figure . :

 : Figure . : Horizontal displacement (in m) for λ = 3.5 plo ed on undeformed con guration.

  (a) Particle velocity magnitude in m s -1 for a stretch of 1.7 and crack speed of 17 m s -1 . Particle velocity magnitude in m s -1 for a stretch of 3.5 and crack speed of 56 m s -1 .

Figure . :

 . Figure . : Particle velocity magnitudes in sub-and tran-sonic cases.

Figure . :

 . Figure . : Comparison of particle velocity magnitude at about mm above the crack path for stretches . and . .

Figure . :

 . Figure . : Particle velocity magnitudes in sub-and tran-sonic cases.

Figure . : 2 3

 .2 Figure . : Comparison of F at about mm above the crack path for stretches . and . .

  (a) Crack initiation at the right end. (b) Just a er initiation. (c) About mm from the right. (d) In the middle of the specimen. (e) About mm from the le end. (f) e le end. Locations where the roughness were measured.

Figure

  Figure . : Variation in surface roughness for the case λ = 2.0.

  (a) Crack initiation at the right end. (b) About mm from the right end. (c) About mm from the right end. (d) About mm from the right end. (e) In the middle of the specimen. (f) e le end. Locations where the roughness were measured.

Figure

  Figure . : Variation in surface roughness for the case λ = 3.0.

Figure

  Figure . : Experimental results vs Model prediction, Uniaxial and Pure Shear case.

Figure . :

 . Figure . : Deformed shape comparison, λ = 3.5. e black lines indicate the deformed shape from experiment.

  (a) Residual force (in N) a er initial stretching (λ = 3.5). (b) Residual force (in N) with changed scale a er initial stretching (λ = 3.5).

Figure . :

 . Figure . : Residual nodal forces plo ed on deformed con guration a er initial stretching for λ = 3.5.

  (a) Residual force (in N) a er initial stretching (λ = 3.5). (b) Residual force (in N) with changed scale a er initial stretching (λ = 3.5).

Figure . :

 . Figure . : Residual nodal forces plo ed on deformed con guration a er initial stretching for λ = 3.5.

Figure . :

 . Figure . : Residual nodal forces plo ed on deformed con guration a er initial stretching for λ = 2.5.

Figure . :

 . Figure . : Data extraction from experiments as input to FE model.

Figure . :

 . Figure . : Particle velocity magnitude in m s -1 for λ = 1.7 plo ed on undeformed con guration.

Figure

  Figure . : particle velocity magnitude in m s -1 for λ = 3.5 plo ed on undeformed con guration.

Figure . :

 : Figure . : Horizontal displacement (in m) for λ = 1.7 plo ed on undeformed con guration (hyperelastic bulk).

Figure . :

 : Figure . : Horizontal displacement (in m) for λ = 3.5 plo ed on undeformed con guration (hyperelastic bulk).

Figure . :

 . Figure . : Particle velocity magnitude in m s -1 for λ = 3.5 plo ed on undeformed con guration (viscoelastic bulk).

Figure . :

 . Figure . : Velocity magnitude in m s -1 for λ = 3.5 plo ed on undeformed con guration.

Figure

  Figure . : Horizontal displacement (in m) for λ = 3.5 plo ed on undeformed con guration (viscoelastic bulk).

  Particle velocity magnitude (m s -1 ) for λ = 3.5 plo ed on deformed con guration.

Figure . :

 . Figure . : Cauchy stress (σ ) in Pa for λ = 3.5 plo ed on undeformed con guration.

Figure . :

 . Figure . : Velocity magnitude in m s -1 for λ = 2.5 plo ed on undeformed con guration.

Figure . :

 . Figure . : Cauchy stress (σ ) in Pa for λ = 2.5 plo ed on undeformed con guration.

  Horizontal displacement (in mm) from viscoelastic simulation.

Figure . :

 . Figure . : Horizontal displacement elds for a di erent data extraction location, λ = 3.5.

  Velocity (in m s ) from viscoelastic simulation.

Figure . :

 . Figure . : Velocity magnitude elds for a di erent data extraction location, λ = 3.5.

Figure . :

 . Figure . : Rubbery and glassy wave speeds in linear viscoelasticity. R stands for Rayleigh, s for shear, d for dilataional. A and B indicate the possible range of the wave speeds that depend on the glassy to rubbery modulus ratio. wave speed based on the elastic properties of the material. In studies where the inertial e ects are included together with Viscous e ects (Graham & Walton, ; Atkinson & Popelar,; Willis, ), the limiting speed is seen to depend on the glassy modulus of the material rather than on the rubbery modulus (Figure.). Hence, the crack is allowed to exceed the rubbery shear wave speed. Such a result can also be observed in Geubelle, Danyluk, and Hilton,, where spectral methods are used to investigate the problem of viscoelastodynamic mode-III crack. Using a rate independent cohesive zone type model, an analysis has been made on the e ect of relaxation times on crack speeds. From the Figure of the reference, it can clearly be seen that the crack speed exceeds the rubbery shear wave speed (c s ), but remains smaller than the glassy shear wave speed (c * s ). However, no speci c comments have been made in regards to the presence or absence of shock fronts in those references.In Fisher and Gurtin, , the propagation of waves of order N in viscoelastic media has been studied. Waves of order N have the solution, u, which is N -1 times continuously di erentiable and exhibit discontinuities in the N t h derivative along a hypersurface. It has been determined in that study that such waves, should they exist, travel with a speed that is derived based on the glassy modulus of the material rather than the rubbery modulus. Shock fronts, by de nition, are waves of order and hence travel with that speed as well. Hence, it might be possible that shock fronts can be observed in viscoelastic material only when the crack speed exceeds even the glassy shear wave speed of the material. In short, a weak discontinuity, should it exist, would be made of high frequency waves at the sharp front. As a consequence of viscoelasticity, this region behaves with a high modulus which in turn roughly raises the shear wave speeds to a value that is (equal to or) greater than the crack speed. is, in turn, prevents a shock wave from developing.Perhaps, the cracks whose speeds exceed c s while still remaining below c *

Figure . :

 . Figure . : Body with a crack.

  Figure . : Cohesive element.

Figure . :

 . Figure . : Structure of SVARS array of the UELMAT subroutine.

  Figure . : Crack speeds vs Stretch. e symbols indicate the experimental results obtained in Corre, Coret, Verron, Leblé, and Le Lay, . e solid line indicates the crack speeds from the simulation obtained using a rate dependent CZM.

  Figure . : Crack speeds vs Stretch, all heights

  Figure . can be compared with Figure of Gent and Marteny,

Figure

  Figure . : Crack speeds vs Stretch, all heights

Figure . :

 . Figure . : Traction separation curve for various crack speeds

Figure . :

 . Figure . : Traction separation curve for various crack speeds

Figure . :

 . Figure . : Comparison of the crack opening pro le between the experiment and simulation for the case of λ = 3.5 when the crack is at about the center of the specimen.

Figure . :

 . Figure . : Particle velocity magnitude in m s -1 for λ = 3.5 plo ed on deformed con guration.

Figure . :

 . Figure . : Evolution of traction within the cohesive zone. Extracted from W. G. Knauss and Losi, . h denotes the thickness of the layer used in that study. α denotes the position of the mathematical tip .

  . ) b 2d is the restriction of b to within the plane. As a consequence of plane stress assumption, b has been assumed to be of the form b = b 2d o o T b 33 . Expressing the rst and the second invariants in terms of in-plane components as (realizing that b 33 = 1/det b 2d )

  term b 2db 33 I 2d in the equation . can now be simply wri en as ∂I 1 ∂b 2d b 2d and the term I 1 (b 2db 33 I 2d ) -((b 2d ) 2b 2 33 I 2d ) as ∂I 2 ∂b 2d b 2d . e total stress in the equation . then simply becomes τ exercise for the viscous branches, the viscous contribution to the total stress becomes . and . become τ = τ e + τ , and ( . )

Figure . :

 . Figure . : Geometry and mesh to test the implementation of UMAT.

  Figure . : σ 22 (in Pa) predicted by UMAT and the polynomial model from Abaqus. en τ A + p = 4C 10 ϵ A and τ B + p = 4C 10 ϵ B . e evolution equations . and . become, with γ 0 = 1,

Figure . :

 . Figure . : σ 22 and ϵ 22 vs time for FLV and FV models. e FLV model has already been implemented within Abaqus while the FV model in the UMAT.

  Figure . : σ 22 and ϵ 22 vs time for FLV and FV models. Small perturbations are imposed on large background strains. ϵ 22 is the 22 component of the logarithmic strain.

Figure . :

 . Figure . : Particle velocity magnitude in m s -1 for λ = 3.5 plo ed on deformed con guration, FV model, γ 0 = 1.

Figure . :

 : Figure . : Horizontal displacement magnitude in mm for λ = 3.5 plotted on undeformed con guration, FV model, γ 0 = 1.

Figure . :

 . Figure . : Velocity distribution along about 5 mm above the crack path for λ = 3.5, γ 0 = 1.

  Figure . : e viscous indicator parameter β for λ = 3.5 plo ed on deformed con guration, FV model, γ 0 = 1.

Figure . :

 . Figure . : Particle velocity magnitude in m s -1 for λ = 3.5 plo ed on deformed con guration, FV model, γ 0 = 1/60.

Figure . :

 . Figure . : Velocity distribution along about 5 mm above the crack path for λ = 3.5, various cases.

FigureFigure . :

 . Figure . : e viscous indicator parameter, β for λ = 3.5, plo ed on deformed con guration, FV model, γ 0 = 1/60.

Figure . :

 : Figure . : Horizontal displacement magnitude in mm for λ = 2.5 plo ed on undeformed con guration, γ 0 = 1/60.

Figure . :

 : Figure . : Horizontal displacement magnitude in mm for λ = 2.5 obtained in the experiments.

Figure . :Figure . :

 .. Figure . : Velocity magnitude in m s -1 for λ = 2.5 plo ed on undeformed con guration.

  Figure . : Energy evolution for the case of λ = 3.5 with the shi ed FV model. e crack is at about the center of the specimen at about 3.5 ms.

Figure

  Figure . : Strain (Str), Viscous (Vis), External work (Ext), and Kinetic (Kin) Energy release rates for the data from the experiments.

Figure . :

 . Figure . : Energy evolution for the case of λ = 2.5 with the shi ed FV model. e crack is near the center of the specimen at about 2.2 ms.

Figure

  Figure . : Strain (Str), Viscous (Vis), External work (Ext), and Kinetic (Kin) Energy release rates for the data from FLV simulations with cohesive zone.

Figure . :

 . Figure . : Viscoelastic energy dissipation density for λ = 3.5, crack speed imposed from CZM simulations.

Gent

  new constitutive relation for rubber. Rubber chemistry and technology, ( ), -. Gent, A. ( b). Adhesion and Strength of Viscoelastic Solids. Is ere a Relationship between Adhesion and Bulk Properties? †. Langmuir, e e ect of strain upon the velocity of sound and the velocity of free retraction for natural rubber. Journal of Applied Physics, ., Nguyen, Q. S., & Suquet, P. ( ). Continuum thermodynamics. Journal of Applied Mechanics, Transactions ASME, ( ), -. doi: . / . Geubelle, P., Danyluk, M. J., & Hilton, H. H. ( ). Dynamic mode III fracture in viscoelastic media. International Journal of Solids and Structures, strains at the tip of a crack in a sheet of hyperelastic material: I. Homogeneous case. Journal of Elasticity, ( -), -. doi: . /BF Goldstein, R. V. & Salganik, R. L. ( ). Bri le fracture of solids with arbitrary cracks. lnt. J. Fracture, , . Gori, M., Rubino, V., Rosakis, A. J., & Lapusta, N. ( ). Pressure shock fronts formed by ultra-fast shear cracks in viscoelastic materials. Nature communications, . & Reese, S. ( ). A presentation and comparison of two large deformation viscoelasticity models. Journal of Engineering Materials and Technology, Transactions of the ASME, ( ), -. doi: . / . Kamasamudram, V. & Coret, M. ( a). Investigation of Dynamic Fracture of Elastomers -IIa. doi:h ps://doi.org/ . /zenodo. Kamasamudram, V. & Coret, M. ( b). Investigation of Dynamic Fracture of Elastomers -IIb. doi:h ps://doi.org/ . /zenodo. Kamasamudram, V., Coret, M., & Moës, N. (). e role played by viscoelasticity in the bulk material during the propagation of a dynamic crack in elastomers. . G. & Losi, G. U. (). Crack propagation in a nonlinearly viscoelastic solid with relevance to adhesive bond failure. Journal of Applied Mechanics,

  Journal of Solids and Structures, analysis of extrinsic and intrinsic cohesive models of dynamic fracture. International Journal of Solids and Structures, . J., Lawrence, C. C., & omas, a. G. (). High-Speed Fracture of Elastomers: Part I. Rubber Chemistry and Technology, ( ), -. doi: . / . Lake, G. J., omas, A. G., & Lawrence, C. C. (). E ects of hydrostatic pressure on crack growth in elastomers. Polymer, ., Bouchbinder, E., & Fineberg, J. (). Breakdown of Linear Elastic Fracture Mechanics near the Tip of a Rapid Crack. Physical Review Le ers, . J.,Peters, B., & Beex, L. A. ( ). Rate-dependent phase-eld damage modeling of rubber and its experimental parameter identi cation. Journal of the Mechanics and Physics of Solids, model of rubber viscoelasticity. Mechanics Research Communications, ( ), -. Mai, T. T., Okuno, K., Tsunoda, K., & Urayama, K. ( ). Crack-Tip Strain Field in Supershear Crack of Elastomers. ACS Macro Le ers, ( ), -. doi: . /acsmacrole . c Nguyen, N. & Waas, A. M. ( ). Nonlinear, nite deformation, nite element analysis. Zeitschri fur Angewandte Mathematik und Physik, ( ). D. & Govindjee, S. [S.]. ( ). Numerical study of geometric constraint and cohesive parameters in steady-state viscoelastic crack growth. the response of rubbers at high strain rates-III. E ect of hysteresis. Journal of the Mechanics and Physics of Solids, deformation isotropic elasticity-on the correlation of theory and experiment for incompressible rubberlike solids. Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences, K., Kaliske, M., Lin, G., & Bhashyam, G. ( ). Evaluation of energy contributions in elasto-plastic fracture: A review of the con gurational force approach. . H., De Borst, R., Brekelmans, W. A., & Geers, M. G. ( ). Localisation issues in local and nonlocal continuum approaches to fracture. European Journal of Mechanics, A/Solids, . N. J., Albohr, O., Heinrich, G., & Ueba, H. ( ). Crack propagation in rubber-like materials. Journal of Physics Condensed Ma er, ( ). doi: . / -/ / /R Persson, B. N. & Brener, E. A. ( ). Crack propagation in viscoelastic solids. Physical Review E -Statistical, Nonlinear, and So Ma er Physics, ( ), -. doi: . /PhysRevE. . Petersan, P. J., Deegan, R. D., Marder, M., & Swinney, H. L. ( ). Cracks in rubber under tension exceed the shear wave speed. Physical Review Le ers, ( ), -. doi: . /PhysRevLe . . . arXiv: [cond-mat] Rahulkumar, P., Jagota, A., Bennison, S. J., & Saigal, S. ( ). Cohesive element modeling of viscoelastic fracture: Application to peel testing of polymers.International Journal of Solids and Structures, , K. & Knauss, W. ( a). An experimental investigation into dynamic fracture: I. Crack initiation and arrest. International Journal of Fracture, ( ), -. doi: . /BF Ravi-Chandar, K. & Knauss, W. ( b). An experimental investigation into dynamic fracture: III. On steady-state crack propagation and crack branching. International Journal of Fracture, ( ), -. doi: . /BF Ravi-Chandar, K. & Yang, B. ( ). On the role of microcracks in the dynamic fracture of bri le materials. Journal of the Mechanics and Physics of Solids, nite element method for stability problems in nite elasticity. International journal for numerical methods in engineering, Independent Integral and the Approximate Analysis of Strain Concentration by Notches and Cracks. Journal of Applied Mechanics, elastic deformations of isotropic materials IV. Further developments of the general theory. Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences, ( ), -. Rivlin, R. S. [R. S.] & omas, A. G. ( ). Rupture of Rubber. I. Characteristic Energy for Tearing. Journal of Polymer Science, . S. [R. S.] & omas, A. G. ( ). Rupture of Rubber. I. Characteristic Energy for Tearing. Collected Papers of R.S. Rivlin, . S. [Ronald S] & Saunders, D. W. ( ). Large elastic deformations of isotropic materials VII. Experiments on the deformation of rubber. Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences, ( ), -. Rosakis, A. J. ( ). Intersonic shear cracks and fault ruptures. Advances in Physics, ( ), -. doi: . / Rosakis, A. J., Samudrala, O., & Coker, D. ( ). Intersonic shear crack growth along weak planes. Materials Research Innovations, ( ), -. doi: . / s Runesson, K., Larsson, F., & Steinmann, P. ( ). On energetic changes due to con gurational motion of standard continua. International Journal of Solids and Structures, amplitude waves in nonlinear elastodynamics and related theories: A personal overview. In Cism international centre for mechanical sciences, courses and lectures (Vol. A theory of crack initiation and growth in viscoelastic media -I. eoretical development. International Journal of Fracture, ( ), -. doi: . /BF Schapery, R. A. ( b). A theory of crack initiation and growth in viscoelastic media -III. Analysis of continuous growth. International Journal of Fracture, theory of crack initiation and growth in viscoelastic media II. Approximate methods of analysis. International Journal of Fracture. doi: . /BF Sharon, E. & Fineberg, J. ( ). Microbranching instability and the dynamic fracture of bri le materials. Physical Review B, basic problems in fracture mechanics and new concepts. Engineering fracture mechanics, ( ), -. Simha, N. K., Fischer, F. D., Shan, G. X., Chen, C. R., & Kolednik, O. ( ). J-integral and crack driving force in elastic-plastic materials. Journal of the Mechanics and Physics of Solids, ( ), -. doi: . /j.jmps. . . Simha, N. K., Fischer, F., Kolednik, O., & Chen, C. ( ). Inhomogeneity e ects on the crack driving force in elastic and elastic-plastic materials. Journal of the Mechanics and Physics of Solids, a fully three-dimensional nite-strain viscoelastic damage model: Formulation and computational aspects. Computer Methods in Applied Mechanics and Engineering, of material forces to hyperelastostatic fracture mechanics. I. Continuum mechanical se ing. International Journal of Solids and Structures, A., Lubrecht, T. A., & Combescure, A. ( ). Damage model with delay e ect analytical and numerical studies of the evolution of the characteristic damage length. International Journal of Solids and Structures, the nite motions generated by a mode I propagating crack. Journal of Elasticity, ( ), -. doi: . /A: Taylor, R. L., Pister, K. S., & Goudreau, G. L. ( ). ermomechanical Analysis of Viscoelastic Solids. (February ), -. Treloar, L. R. G. (

  hyperélastiques pour le comportement mécanique des élastomères. Techniques de l'ingénieur, AM V ( ), . VIC-D, Correlated Solutions. ( ). Washabaugh, P. D. & Knauss, W. ( ). A reconciliation of dynamic crack velocity and Rayleigh wave speed in isotropic bri le solids. International Journal of Fracture, ( ), -. doi: . /BF Weber, G. & Anand, L. (

∫

  m a t t a n = MatrixXd : : Zero ( , ) ; m a t t a n v o l = MatrixXd : : Zero ( , ) ; R o t a t i o n m a t = MatrixXd : : Zero ( , ) ; } v o i d v i s c o u s b r a n c h : : s e t i n e l a s t i c s t r a i n ( MatrixXd C i ) { gamma = . ; / / U s e r i n p u t C i = C i ; e t a = . * r e l a x a t i o n t i m e * c [ ] * gamma ; } ; v o i d v i s c o u s b r a n c h : : c o m p u t e b e t r ( MatrixXd F ) { b e t r = F * ( C i . i n v e r s e ( ) * F . t r a n s p o s e ( ) ) ; E i g e n S o l v e r <M a t r i x d > e s ( b e t r ) ; v = e s . e i g e n v e c t o r s ( ) . c o l ( )

  r e a l ( ) ; R o t a t i o n m a t ( , ) = v ( ) * v ( ) ; R o t a t i o n m a t ( , ) = v ( ) * v ( ) ; R o t a t i o n m a t ( , ) = v ( ) * v ( ) ; R o t a t i o n m a t ( , ) = v ( ) * v ( ) ; R o t a t i o n m a t ( ,) = v ( ) * v ( ) ; R o t a t i o n m a t ( , ) = v ( ) * v ( ) ; R o t a t i o n m a t ( , ) = . * v ( ) * v ( ) ; R o t a t i o n m a t ( , ) = . * v ( ) * v ( ) ; R o t a t i o n m a t ( , ) = v ( ) * v ( ) + v ( ) * v ( ) ; R o t a t i o n m a t tr a n s p = R o t a t i o n m a t . t r a n s p o s e ( ) ; } ; v o i d v i s c o u s b r a n c h : : c o m p u t e i n v a r ( MatrixXd C , b o o l a ) { i f ( a ) { d o u b l e c = . / C . d e t e r m i n a n t ( ) ; I = i n v a r [ ] = C ( , ) + C ( , ) + c ; M a t r i x d CC = C * C ; I = i n v a r [ ] = . * ( i n v a r [ ] * i n v a r [ ] -(CC ( , ) +CC ( , ) + c * c ) ) ; } i f ( ! a ) { d o u b l e c = . / C ( ) / C ( ) ; I = i n v a r [ ] = C ( ) + C ( ) + c ; I = i n v a r [ ] = C ( ) * C ( ) +C ( ) * c + c * C ( d v i s c o u s b r a n c h : : c o m p u t e d e r i v a t i v e ( ) { I = i n v a r [ ] ; I = i n v a r [ ] ; d e r i v a t i v e [ ] = c [ ] + . * c [ ] * ( I -) + . * c [ ] * s q ( I -) + . * c [ ] * ( I -) * ( I -) ; d e r i v a t i v e [ ] = c [ ] * s q ( I -) ; d I d A = . -. / s q ( lambda A ) / lambda B ; d I d B = . -. / s q ( lambda B ) / lambda A ; d I d A = lambda B -. / s q ( lambda A ) ; d I d B = lambda A -. / s q ( lambda B ) ; } ; d W d I I = d W d I I ; S t i f f p r i n c i p a l ( , ) = . * dW dI * d AI dAdA + . * lambda A * d I d A * ( d W d I I * d I d A + d W d I I * d I d A ) ; S t i f f p r i n c i p a l ( , ) + = . * dW dI * d AI dAdA + . * lambda A * d I d A * ( d W d I I * d I d A + d W d I I * d I d A ) ; S t i f f p r i n c i p a l ( , ) * = * lambda A ; S t i f f p r i n c i p a l ( , ) = . * dW dI * d AI dAdB + . * lambda A * d I d A * ( d W d I I * d I d B + d W d I I * d I d B ) ; S t i f f p r i n c i p a l ( , ) + = . * dW dI * d AI dAdB + . * lambda A * d I d A * ( d W d I I * d I d B + d W d I I * d I d B ) ; S t i f f p r i n c i p a l ( , ) * = * lambda B ; S t i f f p r i n c i p a l ( , ) = . * dW dI * d B I d B d B + . * lambda B * d I d B * ( d W d I I * d I d B + d W d I I * d I d B ) ; S t i f f p r i n c i p a l ( , ) + = . * dW dI * d B I d B d B + . * lambda B * d I d B * ( d W d I I * d I d B + d W d I I * d I d B ) ; S t i f f p r i n c i p a l ( , ) * = * lambda B ; S t i f f p r i n c i p a l ( , ) = . * dW dI * d BI dBdA + . * lambda B * d I d B * ( d W d I I * d I d A + d W d I I * d I d A ) ; S t i f f p r i n c i p a l ( , ) + = . * dW dI * d BI dBdA + . * lambda B * d I d B * ( d W d I I * d I d A + d W d I I * d I d A ) ; S t i f f p r i n c i p a l ( , ) * = * lambda A ; } ; v o i d v i s c o u s b r a n c h : : c o m p u t e t a n g e n t p r e s s u r e ( ) { p r e s s u r e t a n g e n t ( , ) = -. / . * ( S t i f f p r i n c i p a l ( , ) +

  d v i s c o u s b r a n c h : : c o m p u t e t a n g e n t n r p r i n c i p a l ( ) { c o m p u t e t a n g e n t p r e s s u r e ( ) ; T a n g e n t p r i n c i p a l = d e l t a t / . / e t a * ( S t i f f p r i n c i p a l + p r e s s u r e t a n g e n t ) + MatrixXd : : I d e n t i t y ( , ) ;} ; v o i d v i s co u s b r a n c h : : c o m p u t e r e s i d u a l p r i n c i p a l ( ) { c o m p u t e i n v a r ( e i g s , ) ; c o m p u t e d e r i v a t i v e ( ) ; c o m p u t e s t r e s s t a u p r i n c i p a l ( ) ; I = i n v a r [ ] ; I = i n v a r [ ] ; dW dI = d e r i v a t i v e [ ] ; dW dI = d e r i v a t i v e [ ] ; t a u d e v = t a u p r i n c i p a l . a r r a y ( ) +p ; r e s p r i n c i p a l = e p s i l o n + d e l t a t / . / e t a * t a u d e ve p s i l o n t r ; } ; v o i d v i s c o u s b r a n c h : : m a t t a n p r i n c i p a l ( b o o l e l a s t i c ) {

  Contents I Literature review Introduction to fracture mechanics . Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Crack Propagation condition . . . . . . . . . . . . . . . . . . . . . . Weakly Nonlinear eory of fracture . . . . . . . . . . . . . . . . Additional experiments . Results of 20 mm and 60 mm samples . . . . . . . . . . . . . . . . . Other experiments on the mm specimen . . . . . . . . . . . . . . Samples of another batch (Batch ) . . . . . . . . . . . . . . . . . . Temperature measurements . . . . . . . . . . . . . . . . . . . . . . Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . role of viscoelsaticity in the bulk . Using momentum balance . . . . . . . . . . . . . . . . . . . . . . . Using the cracks speed as an input . . . . . . . . . . . . . . . . . . Discussions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . zone model . Bulk material . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Rate independent cohesive zone model . . . . . . . . . . . . . . . . Viscous Cohesive zone model . . . . . . . . . . . . . . . . . . . . . Finite Element Implementation . . . . . . . . . . . . . . . . . . . . . Crack speeds prediction and traction -separation relations . . . . . Discussions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	Cohesive Finite Viscoelastic model
	.

. Crack tip elds in elastomers . . . . . . . . . . . . . . . . . . . . . . Transonic cracks . . . . . . . . . . . . . . . . . . . . . . . . . . . . Experimental and theoretical studies . Early studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 'Instabilities' in dynamic fracture . . . . . . . . . . . . . . . . . . . Fracture of elastomers . . . . . . . . . . . . . . . . . . . . . . . . . . Experimental and theoretical studies on polymer fracture . . . . . . Energy release rate through integrals . . . . . . . . . . . . . . . . Polymers . Elastomers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Nonlinear elasticity . . . . . . . . . . . . . . . . . . . . . . . . . . . viscoelasticity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . II Experiments Fracture of polyurethane . Experimental setup . . . . . . . . . . . . . . . . . . . . . . . . . . . Procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Complex geometry and crack path . . . . . . . . . . . . . . . . . . . Comparison with literature . . . . . . . . . . . . . . . . . . . . . . Revisiting the fracture experiments . Displacement, Velocity and Strain elds . . . . . . . . . . . . . . . . Observations on the variation of surface roughness . . . . . . . . . Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . III Simulations e ermodynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . Integration of the evolution equation . . . . . . . . . . . . . . . . . Plane stress formulation . . . . . . . . . . . . . . . . . . . . . . . . Model checks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Application of FV model to the experiments . . . . . . . . . . . . . Energy budget in viscoelastodynamic fracture . . . . . . . . . . . . Discussions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A Plane stress version of Finite linear viscoelastic model A. Time integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . A. Computation of tangent . . . . . . . . . . . . . . . . . . . . . . . . B Additional results for Finite viscoelastic model B. Stress components for the hyperelastic case . . . . . . . . . . . . B. Comparison of relaxation tests between FLV and FV . . . . . . . . C C++ implementation of the UMAT

  Elastomers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Nonlinear elasticity . . . . . . . . . . . . . . . . . . . . . . . . . . . viscoelasticity . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	.
	Elastomers are polymers that are generally amorphous comprised of long chains of macro-molecules (A. G. Holzapfel, ; Verron,

  As already mentioned, in Chen et al.,, it was reported that the crack speeds remained independent of the specimen height in the Transonic regime. No other studies exist that report a similar observation. Hence, experiments have been conducted on the Polyurethane elastomers of di erent geometries in addition to what has been done in Corre,

Results of 20 mm and 60 mm samples . . . . . . . . . . . . . . . . . Other experiments on the mm specimen . . . . . . . . . . . . . . Samples of another batch (Batch ) . . . . . . . . . . . . . . . . . . Temperature measurements . . . . . . . . . . . . . . . . . . . . . . Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  Crack propagation in Polyurethane sheets has been studied in corre . e specimens that are mm tall were stretched and a crack was introduced in the specimen. Transonic cracks were observed in that study starting at a stretch value of . . Transonic cracks were rst reported in the studies of Petersan et al., . In that study, natural latex rubber sheets were stretched and a crack is introduced. It was observed that at a stretch level of about . , the crack speed was able to exceed the shear wave speed. e authors ofPetersan et al., remark

Displacement, Velocity and Strain elds . . . . . . . . . . . . . . . . . Horizontal displacement . . . . . . . . . . . . . . . . . . . . . Velocity elds . . . . . . . . . . . . . . . . . . . . . . . . . . . Strain elds . . . . . . . . . . . . . . . . . . . . . . . . . . . Observations on the variation of surface roughness . . . . . . . . . Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  notable studies are that of Corre et al., ; Mai et al., . In both the studies, Transonic cracks were observed. Mai et al., considers the crack opening to be of a wedge shape and to resemble a shock front as in Petersan et al., . e experiments of Corre et al., will be analyzed in this thesis to look for shock fronts in the Transonic regime.

Table . :

 . Model Co-e cients

	Parameter C 10	C 20	C 30	C 21
	Value (Pa)	. E	-.	E	.	.

Table . :

 . Viscoelastic model parameters

	Branch(i) i	.	.	.	.	.	.	.	.	.	.

  crack speeds have been observed to exceed the shear wave speed in other studies in the literature like Petersan et al., ; Chen et al., and Mai et al., . e material tested in Chen et al., is Latex rubber and the tests were performed at • C to prevent strain crystallization. SBR was tested in Mai et al., . e current study is performed on Polyurethane. e strain elds reported in Mai et al., do not exhibit any jumps as is the case in the current study even for crack speeds in Transonic regime. Chen et al., does not report any such shock fronts as well, but instead report that the crack faces in this regime are wedge shaped (which were noted to resemble a shock front in Petersan et al.,

  In Finite Element simulations, Cohesive Zone Models (CZM) are used extensively to predict the crack growth in scenarios where the crack path is known in advance. is might be in the case of a crack growth in mode-II along a weak

	interface Needleman, testing Rahulkumar, Jagota, Bennison, and Saigal, or delamination Corigliano and Ricci, ; Geißler, Kaliske, Nase, and or for peel Grellmann,

. Bulk material . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Rate independent cohesive zone model . . . . . . . . . . . . . . . . . Specimen geometry and loading . . . . . . . . . . . . . . . . . Observed crack speeds with di erent cohesive zone parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Viscous Cohesive zone model . . . . . . . . . . . . . . . . . . . . . Finite Element Implementation . . . . . . . . . . . . . . . . . . . . . . Computation of residual . . . . . . . . . . . . . . . . . . . . . Computation of tangent . . . . . . . . . . . . . . . . . . . . . Implementation into user subroutine . . . . . . . . . . . . . Crack speeds prediction and traction -separation relations . . . . . Discussions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Table . :

 . FV model parameters with γ 0 = 1 60 .

	Branch(i) Sti ness ratio, ( i ) rubbery Relaxation time, τ i (s)	. E-	. E-	. E-	. E-	. E-	. E-	. E-
	necessary (by examining the ratio	ψ i

Table . :

 . Energy rates with the shi ed FV model in MJ m -1 s -1 .

		λ	dSE dt	dD dt	dKE dt	P ex t
		2.5 -2.81 1.03 0.18 -1.6 3.5 -7.4 3.7 0.08 -3.62
	material in the vicinity of the crack path that has been excluded.
	Dividing -	d(SE + KE)

Table . :

 . Energy release rates for the experiments in kJ/m 2 .

	λ	-2 w 0	dSE dt	2 w 0	dD dt	2 w 0	dKE dt	-2 w 0 P ex t ψh 0	ψ h 0
	2.5	130	49	8.5	76	166	120
	3.5	260	129	2.0	127	348	268

Table . :

 . As earlier, the time derivatives of the quantities can be seen in the table . . e energy release rates can be computed and can be seen in the table . , plo ed in the gure . . Energy rates with the shi ed FV model for the FE simulations in MJ m -1 s -1 .

					λ	dSE dt	dD dt	dKE dt	P ex t = P coh
					2.5 -2.63 2.02 0.07 -0.54 3.5 -7.18 6.22 0.06 -0.9
	-1 ) Energy per thickness (J m	0 0.5 1 1.5	1.5 •10 4	1.6 Strain energy 1.7 Dissipation energy 1.8 Kinetic energy	1.9 time (s)	2	2.1	2.2	2.3 •10 -3
	-1 ) Energy per thickness (J m	0 1 2 3	1 •10 4	1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 time (s) Strain energy Dissipation energy Kinetic energy	2 •10 -3	Figure . : Energy evolution for the case of λ = 3.5 with the shi ed FV model. e crack is at about the center of the specimen at about 1.7 ms.

Table . :

 . Energy release rates for the FE simulations in kJ/m 2 .

	λ	-2 w 0	dSE dt	2 w 0	dD dt	2 w 0	dKE dt	G coh = -2 w 0 P ex t ψh 0	ψ h 0
	2.5	119	91	3.1	24	166	127
	3.5	250	220		2	30	348	252

  1 is constant during the global iterations, γ 0 can be viewed as a function of (C e ) t r to evaluate the derivative ∂ϵ Ae ∂(ϵ Be ) t r

  Nonlinear Viscoelastic Solids-A Review. doi: . / Xu, D. B., Hui, C. Y., & Kramer, E. J. ( ). Interface fracture and viscoelastic deformation in nite size specimens. Journal of Applied Physics, ( ), Kaliske, M. [Michael]. (). Fracture simulation of viscoelastic polymers by the phase-eld method. Computational Mechanics. doi: . Zhang, H. P., Niemczura, J., Dennis, G., Ravi-Chandar, K., & Marder, M. (). Toughening e ect of strain-induced crystallites in natural rubber. Physical Review Le ers, ( ), -. doi: . /PhysRevLe . . Zreid, I., Fleischhauer, R., & Kaliske, M. (). A thermomechanically coupled viscoelastic cohesive zone model at large deformation. International Journal

	e PK stress in Simo,	's model is given by
	S = -p C -1 + DEV t	∫ t -∞	(t -s)	∂ ∂s	DEV s 2	∂ψ ∂C	(s) ds ,	(A. )
	of Solids and Structures, ( -),	-		. doi: .		/j.ijsolstr.	. .
	Computer Methods in Applied Mechanics and Engineering, ( ), doi: . / -( ) -Willis, J. ( ). Crack propagation in viscoelastic media. Journal of the Mechanics -. and Physics of Solids, ( ), -. doi: . / -( ) -Wineman, A. ( ). -. doi: . / . Xu, X. P. & Needleman, A. ( ). Numerical simulations of fast crack growth in bri le solids. Journal of the Mechanics and Physics of Solids, ( ), -. doi: . / -( ) -Yavari, A. ( ). A geometric theory of growth mechanics. Journal of Nonlinear Science, ( ), -. doi: . /s ---y. arXiv: . Yin, B. & /s ---

  and DEV t DEV t • = DEV t • , since DEV is a projector. A superscript has been added to indicate the time at which the quantity is evaluated. e deviatoric projector inside the integral can be evaluated as It has been made clear that the quantities are to be evaluated at the time s. e deviator at time t can now be wri en as ) tr((C s ) -1 C t )(C t ) -1 ds.

	t -∞ Hence, where, DEV 2 (t-s) ∂ ∂s ∂ψ DEV s 2 ∂C (-∞) = O has been used. De ning S 2 ∂ψ ∂C ds = (t-s)DEV 2 ∂ψ ∂C t -∞ -∫ t -∞ (0)DEV 2 ∂ψ ∂C (t) -∫ t ∂ ∂s (t -s)DEV 2 ∂ ∂s (t-s)DEV s 2 ∂ψ ∂C , (A. ) ∂ψ ∂C . (A. ) -∞ ∂ψ ∂C , the rst equation becomes S = -pC -1 + (0)DEV t S(t) -DEV t ∫ t -∞ h(t -s)DEV s Ss ds, (A. ) p = (0) St 33 C t 33 -1 3 tr( St C t ) -∫ t -∞ h(t -s) Ss 33 C t 33 ds+ 1 3 ∫ t h(t -s) tr( Ss C t ) ds -∞ + 1 3 ∫ t -∞ h(t -s) tr( Ss C s )(C s 33 ) -1 C t 1 ∫ t 33 ds -9 -∞	=
	where h(t -s)	∂ ∂s	(t -s) DEV s	Ss =	Ss -	1 3	tr(	Ss C s )(C s ) -1 .	(A. )
		e integral hence becomes
	∫ t -∞	h(t -s)DEV s	Ss ds =	∫ t -∞	h(t -s)	Ss ds -	1 3	∫ t -∞
	-	DEV t 1 3 ∫ t -∞	∫ t -∞ h(t -s) tr( h(t -s)DEV s ∫ t -∞ h(t -s) Ss C s )(C s ) -1 ds+ Ss ds = Ss ds -1 3 ∫ t -∞ 1 9 ∫ t -∞	h(t -s) tr( h(t -s) tr( Ss C s (A. ) Ss C t )(C t ) -1 ds

h(ts) tr( Ss C s )(C s ) -1 ds. (A. ) h(t -s) tr( Ss C s ) tr((C s ) -1 C t ) ds.

  1 and I 2 are the invariants of C, since the material is incompressible. e inplane stress components can now be expressed in terms of the invariants of the in-plane strain tensor as C 33 )I 2d -C 2d , (A. )whereC 2d = F T 2d F 2d . I 2d 1 = tr(C 2d ), C 33 = 1/det(C 2d ).Hence, is the glassy strain energy density functional now expressed completely in terms of the invariants of C 2d . e expression can be seen to be ds. (A. ) Inserting all the above in equation A. , the total stress can be obtained to be

	Integrating the above by parts and realizing that h(t -s) = equation becomes	∂ ∂s	(t -s), the above
	-(0) 33 ) -1 S 2d = Ŝt + ∫ t (t -s) ∂ ∂s DEV 2D Ŝs s -∞ + C t 33 3 (C t 2d ) -1 ∫ t -∞ (t -s) ∂ ∂s tr( Ŝs C s 2d )(C s ∫ t -∞ (t-s) ∂ ∂s DEV 2D s Ŝs + C t 33 3 (C t 2d ) -1 ∫ t -∞ (t-s) ∂ ∂s tr( Ŝs C s 2d )(C s 33 ) -1 (A.
				S2d = 2	∂ψ ∂I 1	I 2d + 2 1 + S33 = 2 ∂ψ ∂I 2 (I 2d ∂ψ ∂I 1 + 2 ∂ψ ∂I 2	I 2d 1	(A. )
	De ning
	where ψ Ŝ = 2	∂ψ ∂I 1	I -	1 det C 2d	C -1 2d + 2	Ŝ = 2 ∂ψ ∂I 2 I 2d ∂ψ ∂C 2d 1 I -	, det C 2d I 2d 1	C -1 2d +	1 det C 2d	(A. ) I -C 2d . (A. )

  integration technique as in Taylor et al., will be used to integrate the above equation. Since (t) = ∞ + 1 exp(-t τ ), the integration technique will be demonstrated on the exponential term. De ning I (t) ∫ t -∞ e -t -s

					s)	∂ ∂s	h(s).			(A. )
										τ	∂ ∂s	h(s),
	I (t) = = e ∫ t -∆t t -∆t τ I (t -∆t) + e -t -s τ ∂ ∂s h(s) + ∫ t t -∆t e second term in the above equation will be integrated approximately. ∫ t e -t -s τ ∂ ∂s h(s) t -∆t e -t -s τ ∂ ∂s h(s).	(A. ) (A. )
	∫ t t -∆t	e -t -s τ	∂ ∂s	h(s) ≈	h(t) -h(t -∆t) ∆t	e -t τ	∫ t t -∆t	e	s τ	(A. )

e ≈ [h(t)h(t -∆t)] 1e -∆t /τ ∆t/τ . (

A

. )

  N i=1 i = 1. ψ is the glassy strain energy density functional. De ning ψ ∞ ∞ ψ to be the rubbery strain energy functional, the above expression for stress can be simpli ed as

	s) where (s) = ∞ + N ∂ ∂s DEV 2D s i=1 i exp(-s/τ i ) and Ŝs + C t 33 3 (C t 2d ) -1	∫ t -∞	(t-s)	∂ ∂s	tr(	Ŝs C s 2d )(C s 33 ) -1 . (A. )
			Ŝ = 2	∂ψ ∂C 2d	,		(A. )
	S 2d =	Ŝ∞ 2d + + C t 33 N i=1 3 (C t ∫ t -∞ 2d ) -1	i exp(-(t -s)/τ i ) N i=1 ∫ t -∞ i exp(-(t -s)/τ i ) ∂ ∂s DEV 2D s ∂ ∂s	Ŝs tr(	Ŝs C s 2d )(C s 33 ) -1 . (B. )

∞ +

Acknowledgements . ermodynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . Integration of the evolution equation . . . . . . . . . . . . . . . . . Plane stress formulation . . . . . . . . . . . . . . . . . . . . . . . . . Stress and tangent computation . . . . . . . . . . . . . . . . . Tangent computation . . . . . . . . . . . . . . . . . . . . . . . Implementation details . . . . . . . . . . . . . . . . . . . . . Model checks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Hyperelastic case . . . . . . . . . . . . . . . . . . . . . . . . . Viscoelastic case -small strains and small perturbations . . . Viscoelastic case -large strains and small perturbations . . Application of FV model to the experiments . . . . . . . . . . . . . . Model with adjusted parameters . . . . . . . . . . . . . . . Energy budget in viscoelastodynamic fracture . . . . . . . . . . . . . From the experiments . . . . . . . . . . . . . . . . . . . . . . From the simulations with cohesive model . . . . . . . . . . Discussions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

e 33 component of the above deviator is given by

h(t-s) tr( Ss C s ) tr((C s ) -1 C t )(C t 33 ) -1 ds.

(A. )

It has been assumed that the strain tensors can be expressed in the form of A. to write the above equation. e expression for p can now be developed by using All the above can be combined to obtain the total stress. e terms with (0) as coe cient can be evaluated as . ) e restriction of above equation to in-plane components will yield (from equation A. )

(A. ) e le over terms, a er canceling the like terms, can be seen to be

Restricting to in-plane components,

h(ts) tr( Ss C s )(C s 2d ) -1 ds (A. ) e term 2 of A. becomes

e term 3 of A. becomes

e term 4 of A. becomes

) -1 and combining the above equations, A. becomes

Hence, . ) e rst term on the right hand side of A. is evaluated as (denoting it as T 1 )

e second term on the right hand side of A. inside the integral is evaluated as (denoting it as T 2 )

e entire integral can hence be evaluated as

e terms T 1 and T 2 are updated as per equations A. and A. .

A.

e tangent is computed as

. It involves the computation of the derivative of the deviatoric projector.

, where

is the fourth order identity tensor. e derivative of the second term involves the computation of Here, additional results for the plane stress version of the FV model in chapter will be presented.

B.

In section . . , only the 22 component of stress has been compared between the UMAT and Abaqus for the hyperelastic case. Here, the other stress components will be presented. From the gures B. and B. , it can be seen that the stress distributions are exactly the same between the UMAT and the Abaqus's internal implementation.

B.

Similar to what has been done in Govindjee and Reese, , a comparison has been made between the FLV and FV models using a relaxation type test. A similar setup as in Govindjee and Reese, has been used, except that the plane stress conditions are assumed to prevail as opposed to that study. A square sample of dimensions 1 mm × 1 mm is held xed at the bo om and a displacement of δ has been applied in the -direction to the top face in a span of 0.01 s. It is then held there for 30 s. e evolution of -component of the cauchy stress tensor has been monitored. e relaxation function is de ned as E(t) σ (t)/ln(1 + δ /L), where L = 1 mm. e e ect of applied stretch on the relaxation function has been studied.

e relaxation test can be considered to come under the case of large strains and large perturbations, since the strains in the viscous arms are equal to the applied strain at the end of rst step.

e material has been taken to be incompressible. e hyperelastic and the viscoelastic branches have been assumed to be represented by a Neo-Hookean model. is is achieved by taking N = 1 and C 01 = 0 in the polynomial model.

C 01 has been taken to be 10 MPa. C 01 for the viscous branch has been taken to be time (s) B. that the modulus, E, decays faster for larger strains in the case of FV model when compared to the FLV model. As the strains are reduced, the gap between the two models decreases as well. In the case of in nitesimal strain, the models can be seen to predict almost same results. is is a consequence of the shi in (B. ) can be wri en as

where • denotes the terms in the rst and the second integral. It shall be noted that for for all times bigger than t 1 , ∂• ∂s = 0, since the total strain is held constant and the terms inside the derivative are functions of the total strain only. e rst integral on the RHS in the equation B. can be recognized as either of B. or B. . Hence, using the above recursively at time steps 2∆t, 3∆t,• • • ,n∆t,

e total stress for all times greater than t 1 can then be wri en as

)

∞ + 1 = 1 has been used. When n = 0, t = t 1 and the above expression becomes

, the fully relaxed state. Hence, for the FLV model, the relaxation times can be seen to be 'independent of the strain level'. e curves in the gure B. for the FLV model hence merely translate upward for di erent strain levels.

In the case of FV model, in the rst step, when ∆t < t 1 τ , the term γ 0 ∆t 2η D (τ 1 + p) in the equation . tends to 0. So, ϵ Ae = (ϵ Ae ) t r = ϵ A . Hence, the strain in the viscous arms will be equal to the total strain, which implies the glassy state. During the second step, the evolution of the strain in the viscous arms is governed by the evolution equations . and . . However, since the strains in the viscous arms are no longer small, the evolution equations do not condense to the ones in the equation . and hence are non-linear. Hence, it can be seen to be di erent from the equation B. and so the observed di erences in the gure B. . . * D * d e l t a t / rho / C v r e t u r n s t r e s s ( STRESS , t a u f i n a l ) ; r e t u r n t a n g e n t ( DDSDDE , t a n g e n t f i n a l ) ; r e t u r n i n t e r n a l v a r ( STATEV , i v a r , n u m b e r b r a n c h e s ) ; r e t u r n ; } ;

. ) { * SSE = . ; * S C D r a t e = . ; f o r ( i n t i = l e f t ; i < b r a n c h e s ; ++ i ) { * t a u f i n a l = * t a u f i n a l + * ( t a u + i ) ; * t a n g e n t f i n a l = * t a n g e n t f i n a l + * ( t a n g e n t + i ) ; * SSE += * ( S S E + i ) ; * S C D r a t e += * ( SCD+ i ) ;

. u s i n g namespace s t d ; v o i d M a t e r i a l : : s e t m a t e r i a l ( d o u b l e s t i f f r a t i o , d o u b l e r e l a x a t i o n t i m e ) { / / C o n s t r u c t o r s t i f f r a t i o = s t i f f r a t i o ; r e l a x a t i o n t i m e = r e l a x a t i o n t i m e ; f o r ( i n t i = ; i < , v , v , v , v , v ; d