A theory of viscoelastic crack growth-Revisited (Revised)

A theory of viscoelastic crack growth developed nearly five decades ago is generalized to express traction in the so-called fracture process zone or failure zone as a function of the crack opening displacement (COD). In earlier work, except for minor exceptions, traction was specified as a function of location. The new model leads to a nonlinear double integral that has to be solved for the COD before crack growth can be predicted. First, a closed-form, accurate approximation is found for a linear elastic body. We then show that this COD may be easily and accurately extended to linear viscoelasticity using a realistic, broad spectrum creep compliance. An analytical relationship connecting the stress intensity factor to crack speed then follows. Consistent with earlier work, it is defined almost entirely by the creep compliance. Five different failure zone tractions are employed; their differences are shown to have little effect on crack growth other than through a speed shift factor. The Appendix discusses initiation of growth.


Introduction
The subject of crack growth in linear viscoelastic media has received considerable attention, both theoretically and experimentally. This may be seen in some recent and earlier reviews and analyses of the subject (Greenwood 2004;Kaminsky 2014;Knauss 2015;Rodriguez et al. 2020). A widely used idealization of the material that separates as the crack grows is that of a thin layer next to the crack tip called the fracture process zone or failure zone (FZ). This zone is where the continuum material begins to come apart and eventually separates entirely. In all cases known to the author, except for two, in which the FZ is modeled as a continuum, the traction is specified as a function of location. One exception is in Greenwood (2007), but it is limited to a 3-element model and uses an iteration method of solution; the other is in Schapery (1975b, Appendix) but only a highly approximate analytical solution is found. In this paper accurate analytical solutions are obtained using a realistic, broad spectrum compliance for five traction models.. It has been argued (Ciavarella 2020) that a realistic model should express FZ traction as a function of the opening displacement, as in Greenwood (2007). The author agrees, unless, of course, the material fails by yielding at a constant stress (the so-called Dugdale (1960) model, ).
As a result, the author re-examined his early work (Schapery 1973(Schapery , 1975a, as well as that of (Greenwood and Johnson 1981) and was able to construct an explicit and simple solution to the problem using realistic models. We should add that the author is aware of one publication (Kaminsky 2014) in which the problem of predicting the failure-zone traction was solved by dividing the zone into several individual filaments, thus requiring the solution of a large complex set of equations.
The basic mathematical problem for a continuum FZ involves solving a nonlinear single integral equation for an elastic body and then a nonlinear double integral equation for the viscoelastic body. An approximate analytical solution is found, but it is shown to be quite accurate.
First, we consider only elastic behavior in Section 2. The results are used to construct the viscoelastic solution in Section 3 for continuous crack growth. A so-called modified power law is employed to characterize the viscoelasticity. The models used for FZ traction and viscoelasticity are believed to be realistic, but the analysis method is not limited to them. The analysis leads to explicit relationships connecting the crack speed to the stress intensity factor and to the crack opening displacement. . The initiation of growth is briefly discussed in the Appendix. Mathcad is used for all calculations.

Crack tip model
The idealized crack tip region is shown in Fig.1, in which the crack tip is defined to be point P. The FZ is idealized as a thin layer of length α where all nonlinearity and bond breakage exist and may be elastic or a general rate-dependent material. Outside of this zone, the material is assumed to be linearly elastic or viscoelastic and isotropic. Extension to an orthotropic body (Brockway and Schapery 1978) is readily accomplished if its material axes are aligned with those in Fig. 1.
For a linear elastic material, after removing the stress singularity using Barenblatt's (1962) condition, the crack opening displacement (COD) given in Schapery (1975a) is and g C is the compliance for a locally plane strain state, defined in terms of Young's modulus g E and Poisson's ratio g ν . For later work on viscoelasticity, it is helpful to arbitrarily define the elastic properties as those in the short time glassy state, which is denoted with subscript g.
The singularity was removed by selecting the cohesive (or adhesive) traction f σ to satisfy a condition in terms of the stress intensity factor K and the length of the FZ, α: The traction f σ in this paper will be expressed in terms of the opening displacement v, which will then be found as a function of the distance x between a generic location in the FZ and the crack tip, P, Fig.1. Eq.(2.4) is valid for both elasticity and viscoelasticity.
The specific form of the traction used throughout this paper is 0 (1 / ) q fm vv σσ − = + (2.5) Greenwood and Johnson (1981)  It should be noted that the COD and fracture energy in this paper refer to only those above the crack plane, and thus are one-half of what is sometimes used by others in cohesive crack growth.
To find v we must solve the integral equation, Eq. (2.1), except for q=0. The writer is not aware of any existing analytical method for this. While an iterative or other numerical method could be used, as in Greenwood (2007) Expressing Eq. (2.1) using these normalized quantities, which was originally defined in (Schapery (1975a). It has its maximum value of 2 when q=0.
The energy release rate, ERR, and K are related through the equation (e.g.Anderson 2017) At the start of crack growth ERR is equal to the fracture energy. Thus, equating ERR to Γ , and using Eq.(2.9), we find at the critical point, 2 11 2 xI = Γ (2.12) In addition to the requirement that the critical solution satisfies Eq.(2.12), it must also equal the critical displacement at 1 x , We are now able to construct approximate solutions using the following form for all cases studied: 2 1.5 00 14) The exponent 1.5 is required for small x , according to the theory in Barenblatt (1962 For q=1 there is only one undetermined parameter 1 c ( with 0 c =0, 2 1.5 c = and 0 0 x = ) which was found by making 1 at c vv xx = = ; condition Eq.(2.12) was then found to be closely satisfied. The "Find" program in Mathcad was used in the other cases, which took less than one second per case on a laptop computer for each guessed 2 c ; the trial and error method to select 2 c took very little time because the acceptable range of 2 c was small.
The trial solution Eq.(2.14) is used in Eq.(2.8) to predict the solution. If they agree, then the result is the exact solution. If they are very close we consider Eq.(2.14) an acceptable representation. Fig.3 for 1 c v = shows that Eq.(2.14) is indeed acceptable; not shown are the , in which the input and output agree like that for q=3, c v =1.  Greenwood and Johnson (1981).
The largest c v for which Eq.(2.14) remains valid is determined primarily by the length of the traction tail; i.e. the value of ( )/ fc m v σσ . Viscoelasticity solutions, using the method in Section 3, start to degrade significantly when this ratio is less than 6% (which is the value for q=2 The elastic solutions do not begin to degrade until this ratio is somewhat smaller. However, a long tail with a limiting value of 6% may be unrealistic because of microscopic nonhomogeneity with cohesive cracks or surface roughness with adhesive cracks. In view of this, in the next section we omit predictions for q=3 ( 1.5) c v = to illustrate the effect of a long tail.
The viscoelastic opening displacement is from (Schapery 1975a), where 0 tt ≥ and 0 t is the time when the crack tip reaches the generic point X. Additionally,

Continuous crack growth
In this section we address only the problem of continuous growth, not the transient initiation of growth. The latter is briefly covered in the Appendix. Our prediction is simplified by the realistic assumption that the crack speed a  is essentially constant during the time it takes for growth equal to the FZ lengthα . Considering how small α is assumed to be (i.e. inside the singularity zone), this is not a serious restriction. Over the entire amount of growth, the speed a  may vary widely. Without loss in generality for a locally homogeneous body, there is no dependence on X, and we can set 0 0 t = for notational simplicity. Thus, in Eq.(3.1) Ctwill be used in the form of a modified power law, in which normalized quantities have been introduced, where and g e CC are the so-called glassy and equilibrium (elastic) compliances, respectively.
The time e t determines when the equilibrium compliance is approached and accounts for environmental effects, such as temperature and moisture, if the material is rheologically simple; these effects are assumed to be constant during the time taken for the crack to grow the current length of the FZ,α .
Elastic (n=0), viscous (n=1) and power law behavior ( 1) t appear as special cases. The creep compliance for polyurethane rubber in (Mueller and Knauss (1971)  After using normalizations, Eq.(3.1) becomes, where we have defined a normalized crack speed, a  , The problem is to find the v v that satisfies the double integral Eq.(3.9) when 1 xx ≤ .
This will be done by first defining the auxiliary displacement, and introducing an effective compliance, Rewrite Eq.(3.9) using the effective compliance, To motivate the development of a solution (, ) v v xa  to Eq.(3.12), change variables, this matching is done after eff C is derived. Eq.(3.15) is found to be an accurate solution for 1 0 xx ≤≤, given the FZ models in Table 1.
We next observe that the displacement for a propagating crack is needed only at 1 xx = . Let the viscoelastic solution at 1 x be written as, The end of the cohesive zone is designated here by the speed dependent parameter 1v x to distinguish it from that for elastic behavior, 1 x Note that () 1 fa≡  for the Dugdale model because Eq.(3.14) must reduce to 1 11 3.2 Characterization of eff C and f .
The problem has been reduced to finding (, ) eff C xa  and () fa  . We will show that eff C is wellapproximated by the creep compliance itself but shifted by a constant along the log( ) a  axis; this simplicity greatly aids the determination of () fa  . The process is further simplified by neglecting the low speed portion of the compliance for now, using x . Fig.5 shows that it is well-approximated by a simple power law, with p=2 (regardless of speed), in which the speed-dependent coefficient obviously has no effect on eff C . Upon using Eq.(3.20) in Eq.(3.11) we find the effective compliance, where f s is a shift factor on the log( a  axis). It can be expressed in terms of gamma functions, This is the same shift factor reported in (Schapery 1975, Eq.28). Fig. 6 covers the range of n values, as well as p values found later, used in this paper. Fig.5 is plotted using n=0.5 and a  =1; but all three n and all a  give the same agreement. The difference between the two curves in The complete effective compliance will be used to obtain f and make all subsequent predictions so that the long-time equilibrium value is approached smoothly, For the other q-values we find that the auxiliary function is also well-approximated by x , respectively; all resulting p-values are in Table 1.
After selecting the p-value using f=1, () fa  is found by matching the associated input, Eq.(3.15), and output, Eq.(3.12), at 1v x , given a set of a  spaced one-decade apart, and then fitting it to a sixth order polynomial in ( log( ) a  +constant).
Figs.7 and 8 show f and input-output pairs for two selected cases. The results for the other cases provide even closer agreement between input and output, and therefore are not shown. For q=1 and q=2 ( 1 c v = ) the deviation of () fa  and R I (Eq.(3.31)) from unity is less than that in Fig.   7. The speed chosen for Fig.8 is close to the maximum of the f-curves in Fig.7, which produces the largest input-output difference, although barely graphically distinguishable. which appears as a light straight line in Fig.11. The slope of this line for n=0.5, 1/6, is in agreement with the experimental data in Schapery (1975c) (neglecting the small effect of the speed dependence of 1 I and f), which was obtained by Mueller and Knauss (1971). As Fig. 11 shows, the effect of this speed dependence is quite small on logarithmic axes in a realistic range of n for polymeric solids.
Except for q=2, c v =3, there is very little difference in the a  position for each q because all are shifted on the log scale by approximately the same amount relative to the creep compliance.
In the former case the relative shift is noticeable because of the small value of q s in Table 1. The associated long tail in Fig.2 produced this small q s .

Direct comparison of traction models.
It is important to recall that α , x and a  have been defined in terms of physical variables using m σ , Eq.(2.7). Thus, to directly compare the physical response of different models where the correction factor, cf, relative to the Dugdale model is, It's values are listed in Table 1. Thus, for direct comparison of the models, they should all be plotted using common variables, say: The compliance is unaffected by this change because its argument is a ratio of these corrections.
This correction will change the relative horizontal position of each curve in Fig.11, but not the shape. Unless 0 v and c v are known, the specific model cannot be established from experimental data on K vs. a 

Comparison of CODs.
Although there may not be enough difference in () Ka  for the different models to identify the applicable model without precise experimental data, we will see if displacement measurements from α out through the far-field (singularity) zone will help. This zone is defined roughly by 5 x α < normalized distance from the crack tip to the nearest geometric feature.
where 1 I is in Eq. showing that the displacement shape can deviate greatly from elastic behavior. Fig.12 for n=0.5 illustrates the difference between the displacements for the Dugdale (q=0) and q=3 models. For both (a) elastic and (b) strong viscoelastic effects, the Dugdale COD joins the far-field solution at approximately five times the length of the FZ (0.7 decade). The q=3 model blends more quickly at approximately three times (0.5 decade) its FZ length. Note in Fig.12 that the slope of the line for the far field singular zone is 0.5 for elastic behavior, as expected, and approximately unity when 1 D a =  , as predicted from Eq.(3.50). Thus, in this latter case, the COD shape is essentially that of a flat-faced wedge. We have also found that unless measurements of the COD are made within a short distance from the FZ they cannot be used to identify the model-even when there is a large difference in the FZ traction shape. (In Schapery, 1975c, the FZ for the rubber is submicroscopic.) It should be added that the finite element method may be needed for realistic COD comparisons of theory and experiment due to the large local deformations near the crack tip.
As a final observation, we mention that the stress from the crack tip out through the singularity zone (x<0) is affected by the creep compliance due to its dependence on the traction in the FZ (Schapery 1975a, Eq.(17)).

Discussion and Conclusions
We have shown that crack growth in a linear viscoelastic body is defined almost entirely by the creep compliance when the FZ traction is constant or a monotonically decreasing function of the crack opening displacement. The different FZ tractions have only a small effect on log( ) . log( ) K vs a  , except through speed shift factors; generally, the combined shift, fq ss( Fig.6 and We have also examined unreported cases, using a partial sine function for () f v σ , in which the traction vs. COD gradually decreases and may increase near the critical value, with both cases causing the traction to increase as the critical value is reached (cf. Fig.1); the methods used in this paper were found to work well.
All examples utilized a creep compliance that was characterized in an intermediate time period by a single power law. For some materials the intermediate period is very broad, covering many decades, with a gradually changing log-log slope. As shown in Schapery (1975b), the local log-log slope of the creep compliance can be used in place of a single, constant exponent, n to predict the effective compliance because it has a narrow-band weight function.
Crack growth Eq.(3.30) is identical to that in Schapery (1975b) Schapery (1984), instabilities at these lower speeds are also predicted by embedding an FZ inside a thin process zone with viscoelastic properties (different from the continuum) and using the so-called pseudo J-integral as the crack driving force; the craze zone in some plastics could be a candidate fitting this description. This J-integral and crack growth theory for the shear-mode has been used to predict rubber friction with Schallamach waves (Schapery 2020a,b Referring to Fig.6, after growth begins it takes 34 times the initiation time for the crack to grow the length of the FZ for the Dugdale model, and slightly less for the other models. Thus, the time for the start of continuous growth may be negligible in many cases.    Fig.9. Ratio of approximate-to-complete effective compliances in the transition to long-time elastic behavior.