A theory of viscoelastic crack growth: revisited

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 equation 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 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. The author prefers to use the term failure zone because in other more general models there may exist another zone surrounding the FZ where fracture processes start, called the process zone in Schapery (1984). 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 very approximate analytical solution is found. In this paper accurate analytical solutions are obtained using a realistic, broad spectrum compliance for five different 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 failurezone 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 Sect. 2. The results are used to construct the viscoelastic solution in Sect. 3 for continuous crack growth. A socalled 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 then to crack opening displacement (COD) and to length of the FZ. (Viscoelasticity is not limited to the neighborhood of the crack tip, and may exist throughout the body.) 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 a 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 where and C g is the compliance for a locally plane strain state, defined in terms of Young's modulus E g and Poisson's ratio m 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. For an incompressible material, C g is the compliance in simple shear. The singularity was removed by selecting the cohesive (or adhesive) traction r f to satisfy a condition in terms of the stress intensity factor K and the length of the FZ, a: The traction r 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. Equation (4) is valid for both elasticity and viscoelasticity.
The specific form of the traction used throughout this paper is Greenwood and Johnson (1981) used this form with q ¼ 3 to model the effect of intermolecular adhesion. Here, we shall use q = 0, 1, 2, and 3 in example predictions of the theory; the case q = 0 gives the socalled Dugdale traction. The quantity v 0 serves to nondimensionalise the displacement; it is a material constant. The nondimensional displacement v v=v 0 must increase to satisfy condition Eq. (4) as K increases up to the critical value, K c , for the start of crack growth. Its limiting value at the onset of crack growth is v c , which is assumed to be constant along with r m , for the purpose of simplifying the discussion; however, for both elastic and viscoelastic global behavior, they may vary with crack speed without changing the method of analysis. These tractions for q [ 0 are shown in Fig. 2 for v c = 1, and for one case, v c = 3.
The corresponding work input to a given column of material (or to the intermolecular force of adhesion) in the FZ, above the crack plane and with unit cross-section, when it breaks (i.e. when traction at x = a vanishes) is the fracture energy, which is obviously independent of crack speed if q, r m and v c are constant. For simplicity in later calculations the logarithm is avoided by using q ¼ 1 þ 10 À10 in Mathcad.
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.

Elastic analysis
To find v we must solve an integral equation, Eq. (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), this method would greatly complicate the viscoelastic solution. Here, we use an approximate, but accurate, method that is nearly as simple for viscoelastic media as for elastic media.
It is helpful to introduce additional dimensionless quantities. They, along with v v=v 0 , will be called normalized variables, For consistency with the use of x in the development of the theory, the normalized FZ length will be designated as x 1 instead of a.
Expressing Eq. (1) using these normalized quantities, Note that if the variable of integration in Eq. (8) is changed to y ¼ x 0 =x 1 , the result has the form v ¼ gðx=x 1 Þx 1 , which, due to the x 1 factor, does not lend itself to FZ scaling.
Partial normalization of Eq. (4) yields, where which was originally defined in Schapery (1975a). It has its maximum value of 2 when q = 0. The energy release rate in the glassy state, ERR g , and K are related through the equation (e.g. Anderson 2017), At the start of crack growth ERR g is equal to the fracture energy. Thus, equating ERR g to C, and using Eq. (9), we find at the critical point, In addition to the requirement that the critical solution satisfies Eq. (12), it must also equal the critical displacement at x 1 , We are now able to construct approximate solutions using the following form for all cases studied: The exponent 1.5 is required for small x, according to the theory in Barenblatt (1962). The solution for x x 1 enables the prediction for x [ x 1 using Eq. (8). The value of x 0 is found by invoking displacement continuity at this point, while Eq. (14) provides the critical normalized FZ length, Thus, there are three free constants, c 0 ; c 1 ; c 2 , while there are only two constraints, Eqs. (12) and (13). Although one could use a minimum square error between input, Eq. (14), and output, Eq. (8), as the third condition, it would take considerable computational time with the available search method. Instead, we evaluated input and output at 20 evenly spaced xpoints and manually selected c 2 to produce the smallest maximum absolute error over the range 0 x x 1 . For q = 1 there is only one undetermined parameter c 1 (with c 0 = 0, c 2 ¼ 1:5 and x 0 ¼ 0), which was found by making v ¼ v c at x ¼ x 1 ; condition Eq. (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 c 2 ; the trial and error method to select c 2 took very little time because the acceptable range of c 2 was small.
The trial solution Eq. (14) is used in Eq. (8) to predict the solution. If they agree, then the result is the exact solution. If they are very close we consider Eq. (14) an acceptable representation. Figure 3 for v c ¼ 1 shows that Eq. (14) is indeed acceptable; not shown are the results for q ¼ 2ðv c ¼ 3Þ and q = 3 ðv c ¼ 1:5Þ, in which the input and output agree like that for q = 3, v c = 1. The displacement for q = 3 (v c = 1) agrees graphically with that in Greenwood and Johnson (1981). Table 1 lists the constants; for completeness, fitting constants for q = 0 are included as well as constants defined in Sect. 3. Note that the exponent c 2 decreases as q and v c increase; increasing v c beyond those in the table shows that this exponent approaches that for the far-field, 0.5. The values listed for x 1 , and I are for the critical state.
The largest v c for which Eq. (14) remains valid is determined primarily by the length of the traction tail; i.e. the value of r f ðv c Þ=r m . Accuracy of the viscoelasticity solutions, using the method in Sect. 3, starts to degrade significantly when this ratio is less than 6% (which is the value for q = 2 ðv c ¼ 3Þ and q = 3 ðv c ¼ 1:5Þ. 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 ðv c ¼ 1:5Þ but include those for q = 2 ðv c ¼ 3Þ to illustrate the effect of a long tail.

Viscoelasticity
The viscoelastic opening displacement is from (Schapery 1975a where t ! t 0 and t 0 is the time when the crack tip reaches the generic point X. Additionally, C v is the creep compliance and

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 a. Considering how small a 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 t 0 ¼ 0 for notational simplicity. Thus, in Eq. (17) A realistic creep compliance C v ðtÞ will be used in the form of a modified power law, in which normalized quantities have been introduced, where C g and C e are the so-called glassy and equilibrium (rubbery) compliances, respectively. The time t e 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, a.
Elastic (n = 0), viscous (n = 1) and power law behavior ðt ( 1Þ appear as special cases. The creep compliance for polyurethane rubber obtained by Mueller and Knauss (1971), and reproduced in Schapery 1975c), is well-characterized by Eq. (20) with the values n ¼ 0:5; C n ¼ 145:7; t e ¼ 0:84 s ð0 CÞ. It is drawn in Fig. 4 along with the other two compliances used in the examples in this paper. However, the theory will be developed using the general form of Eq. (20).
After using normalizations in Eq. (17), compliance Eq. (20), and normalized crack speed, Eq. (17) where the length of the FZ is now designated with a v subscript, x 1v , to distinguish this speed dependent quantity from that for elastic behavior. Normalization changes Eq. (19) to, Equation (23) becomes, The problem is to find the v v that satisfies the double integral Eq. (25) when x x 1v . This will be done by first defining the auxiliary displacement, v w ðxÞ and then introducing an effective compliance, Rewrite Eq. (25) using the effective compliance, To motivate the development of a solution v v ðx; _ aÞ to Eq. (28), change variables, Thus, The solution to Eq. (30) would be the elastic solution with argument x replaced by C eff ðx; _ aÞx if the upper limit were y 1 ¼ C eff ðx 1v ; _ aÞx 1v . Nevertheless, it motivates us to propose a solution for 0 where f ð0Þ ¼ f ð1Þ ¼ 1, and v A is the elastic displacement, Eq. (14). The coefficient f ð _ aÞ is found by using Eq. (28) to match v v inside the integral (the input) to v v outside (the output) at x 1v ; this matching is done after C eff is derived. Equation (31) is later shown to be an accurate solution for 0 x x 1v , given the FZ tractions in Fig. 2.
We show in Sect. 3.3 that the displacement for a propagating crack is needed only at x ¼ x 1v . Let the viscoelastic solution at x 1v be written as, Using approximation Eq. (14), the solution is simply, where which selects x 1v to be at the end of the FZ as the crack

Characterization of C eff and f
The problem has been reduced to finding C eff ðx; _ aÞ and f ð _ aÞ. We will show that C eff is well-approximated by the creep compliance itself, but shifted by a constant along the logð _ aÞ axis; this simplicity greatly aids the determination of f ð _ aÞ. The process is further simplified by neglecting the low speed portion of the compliance for now, using As the first case, consider the Dugdale traction, for which q = 0. The auxiliary displacement Eq. (26) is the elastic displacement, Eq. (8), although it is dependent on the crack speed through x 1v . Figure 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 C eff . Upon using Eq. (36) in Eq. (27), with the approximate compliance from Eq. (35) utilized, we find the effective compliance, where s f 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 1975b, Eq. (28)). Figure 6 covers the range of n values and p values (including most found later) in this paper. Figure 5 is plotted using n = 0.5 and _ a = 1 to provide coordinate values; but all three n and all _ a give essentially the same agreement. The difference between the two curves in Fig. 5 produces only a 2% difference inĈ eff at x 1v ; the difference was found by evaluating Eq. (27) at x 1v using Eq. (36), and then comparing it with that predicted using the exact auxiliary displacement in Eq. (26) with q = 0.
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, which may be expressed more concisely as where C is the normalized creep compliance, Eq. (20). For the other q-values we find that the auxiliary function is also well-approximated by Eq. (36) in which p is essentially independent of crack speed. The value of p also has negligible dependence on f for its actual range of variation (1 f \1:35) and vice-versa. Thus, Eq. (31), with f ¼ 1, is substituted into Eq. (26), using Eq. (40) for compliance; p and c are selected such that Eqs. (36) and (26) agree closely in shape and match at x 1v , respectively; all resulting pvalues are in Table 1.
After selecting the p-value using f = 1, f ð _ aÞ is found by matching the associated input, Eq. (31), and output, Eq. (28), at x 1v , given a set of _ a spaced onedecade apart, and then fitting it to a sixth order polynomial in (logð _ aÞ ? constant). Figures 7 and 8 show f and input-output pairs for two selected cases; Table 2 lists the polynomial coefficients used for f ð _ aÞ, with n = 0.5, in terms of the variable (logð _ aÞ ? 5). Results for the other cases provide even closer agreement between input and output, and therefore are not shown. For q = 1 and q = 2 (v c ¼ 1) the deviation of f ð _ aÞ and I R (Eq. (48)) 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, but is barely graphically distinguishable.
Although s f was found without the low speed term in the compliance, we shall now show that the approximate effective compliance Eq. (40) is quite good for all _ a, compared to the use of the complete creep compliance, Eq. (20), in deriving C eff . It is sufficient here to express the latter using the single variablet x= _ a:  for auxiliary displacement and COD using n = 1 and _ a ¼ 0:1 Figure 9 shows, in the approach to equilibrium, the ratio R c , for the minimum and maximum values of p and three values of n. In this near-equilibrium period the difference in compliances is quite small, especially for the common range n 0:5. The case n = 1, realistically, characterizes only a viscous fluid, which has no equilibrium compliance; but it is used here to complete the upper end of the thermodynamically allowable range, 0 n 1. If the local logarithmic slope, d logðCÞ=d logðtÞ, is used in the shift factor, the difference is approximately one-half of that using a constant n, as illustrated for n = 0.5 in Fig. 9.
When Eq. (41) is used instead of Eq. (40) to find f ð _ aÞ there is negligible change because of a nearcancelling change in x 1v .

Prediction of crack growth
Equations (9) and (34) enable the prediction of the crack speed vs. stress intensity factor. First, we note that Eqs. (9) and (10) are valid for elastic and viscoelastic materials. In the former case one uses Eq. (16) for x 1 while in the latter case x 1 ¼ x 1v , where x 1v is found from Eq. (34). From now on a will be used instead of x 1v , The expression for stress intensity factor in Eq. (9), along with Eq. (12) for elastic behavior, yields, where I 1 and I 1g is for the glassy state, whose values are those of I 1 in Table 1. This integral is also equal to that in the rubbery (equilibrium) state as it is independent of viscoelastic properties. Next substitute Eq. (44) into Eq. (34) and use Eq. (40) to find, Table 2 Polynomial coefficients for f (n = 0.5) in Fig. 7 n = 0.5, q = 3, v c ¼ 1 n = 0.5, q = 2, v c ¼ 3 C0

C5
-5.8408 9 10 -5 -2.3736 9 10 -6 C6 1.2721 9 10 -6 6.0185 9 10 -6 (a) (b) Fig. 9 Ratio of approximate-to-complete effective compliances in the transition to long-time elastic behavior where Observe that is the normalized time for the crack to propagate the length of the FZ. The shift factor s f reduces this time to an effective time for creep softening. For an elastic body I R ¼ 1, which reduces the right side of Eq. (47) to the energy release rate. Introduce a normalized stress intensity factor, so that Eq. (47) becomes, Next, use Eqs. (48) and (50) in Eq. (44) to express a in terms of quantities in Eq. (51), a ¼ 1 where Values of this factor are in Table 1. Equation (16) shows that where a g is the critical glassy (elastic) FZ length; its values are listed in Table 1 which is drawn in Fig. 10 after using Eq. (56) below to find Kð _ aÞ. The dotted lines are for I R ¼ f ¼ 1, which shows that these functions have only a small effect on a. As expected, a ! 8=s q as _ a ! 1 in view of Eqs. (55) and (61) below.
Finally, the crack growth equation can be written in the simplified form, For the Dugdale model I R ¼ 1(because I 1 ¼ 2) and f ¼ 1 regardless of crack speed.
It is observed that s q appears as a shift factor for crack speed in Eq. (56). Equation (54) shows that it decreases with increasing length of the glassy FZ. With a long traction tail, the portion of the tail with a very small traction is not effective in its relationship to K. This is accounted for in Eq. (55), which can be written in the form, Thus, a long tail with viscoelasticity is reduced by the long tail for the glassy state. Equation (56) provides a relationship between K R and s q s f A that is independent of the FZ properties and depends only on the normalized creep compliance C. It may be solved explicitly for A, given the normalized compliance in Eq. (20). Thus, where Let glassy K g and equilibrium (rubbery) K e define the normalized stress intensity factors corresponding to the high (S = 0) and low(S = 1) speed limits, respectively. Since I R ¼ 1 at these limits, we find, Except for these limits, I R and A are rather involved functions of crack speed. However, this behavior is easily taken into account in predictions by using _ a as the independent variable.
Equation (51), together with Eq. (52), is plotted in Fig. 11 with solid lines while Eq. (56) is plotted using dotted lines; the dotted lines can be interpreted, of (d) (c) Fig. 11 Log(stress intensity factor) vs. log(crack speed); dotted lines are for f and I R ¼ 1 course, as plots of log K À Á vs: log _ a À Á when I R ¼ f ¼ 1: Because I R ¼ f ¼ 1 and s q ¼ 16 for the Dugdale model, the dotted lines in Fig. 11a are also the solution for this model except for a small difference in s f and s q .
In an intermediate speed range Eq. (59) obeys the power law, 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 I 1 and f). 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, v c = 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 s q in Table 1. The associated long tail in Fig. 2 produced this small s q .
When all significant viscoelastic effects are confined to the singularity zone, the behavior is said to be globally elastic. In this case the work/area required to drive the crack is simply the energy release rate based on the equilibrium compliance, which is an increasing function of speed due to speed dependence of K, as drawn in Fig. 11. Its range of variation is given by Eq. (61), However, it is doubtful that the behavior remains globally elastic as the upper limit is approached.

Direct comparison of traction models
It is important to recall that a, x and _ a have been defined in terms of physical variables using r m , Eq. (7). Thus, to directly compare the physical response of different models (assuming C; C v and v 0 are the same for all models) it is necessary to renormalize for the different values of r m using Eq. (6), where the correction factor, cf, relative to the Dugdale model, is Its 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 the corrected values. This correction will change the relative horizontal position of each curve in Fig. 11; but not the shape.
Unless v 0 and v c are known, the specific model cannot be established from experimental data on K vs. _ a alone, unless the data provide an unrealistically precise shape.

Comparison of CODs
Although there may not be enough difference in Kð _ aÞ for the different models to identify the applicable model without precise experimental data, we will see if displacement measurements from a out through the far-field (singularity) zone help. This zone is defined roughly by 5a\x ( normalized distance from the crack tip to the nearest geometric feature, such as a remote free surface or another crack tip. We find from Eqs. (2) and (28) that far behind the crack tip moving at constant speed, where I 1 is in Eq. (45); the values for s f ð0:5; nÞ and I 1 do not change in the conversion to a common length measure for model comparisons. This displacement for large x may not behave as the elastic singular displacement, ffiffi ffi x p . For example, when the compliance obeys a power law, showing that the displacement shape can deviate greatly from elastic behavior. Figure 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 farfield 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 nearly unity (0.9) when _ a D ¼ 1, as expected from Eq. (69). Thus, in this latter case, the COD shape is almost that of a flat-faced wedge; this shape agrees with the experimental COD for rubber reported in (Morishita et al. 2016, Fig. 6, with 0.09 and 0.17% volume fraction of carbon black). 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 andin Morishita et al. 2016, the FZ for the rubber studied is submicroscopic.) In Fig. 12(b) the far field displacements for the two models do not quite coincide primarily because the normalized FZ lengths are somewhat different.
One may expect with rubber that the finite element method would be needed for realistic COD comparisons of theory and experiment due to material nonlinearity and large local deformations near the crack tip. However, Schapery (2021) shows that a modified linear theory accurately predicts the measured CODs.
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)).

Concluding remarks
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ðKÞ vs: log( _ aÞ, except through speed shift factors; generally, the combined shift, s f s q (Fig. 6 and Table 1), decreases with increasing concavity of the FZ traction. Because the shift factors are combined with FZ properties, they cannot uniquely establish the form of the underlying function r f ðvÞ from only experimental data on crack growth. However, very precise measurements of K vs: _ a and COD close to the crack tip may be helpful in providing information on the underlying traction model.
We have also examined cases, using a partial sine function for r f ðvÞ, in which the traction vs. COD gradually decreases and then increases near the critical value (cf. Figure 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 (1975bSchapery ( , 2021, the local log-log slope of the creep compliance can be used in place of a single, (b) (a) Fig. 12 COD comparisons for q s= 3 and Dugdale (q = 0) models with n = 0 constant exponent, n to predict the effective compliance because its evaluation uses a narrow-band weight function; this point is illustrated in Fig. 9.
Crack growth Eq. (47) is identical to that in Schapery (1975b) if I 1 is constant and f = 1. Here, the model enables prediction of the crack speed dependence of I 1 and f.
It should be emphasized that the stress intensity factor itself will depend on global viscoelasticity if external boundary conditions are all or partly specified in terms of displacement history. A simple way of accounting for such conditions and some types of nonlinearity with the models in this paper is through use of pseudo displacements and J-integral (Schapery 1984(Schapery , 1990(Schapery , 2021, as done to predict rubber friction with Schallamach waves (Schapery 2020a, b).
The present model predicts unstable crack growth only when the glassy compliance is reached, or if one or more of the FZ properties decrease with speed at some point. Instabilities at speeds lower than glassy speed in filled elastomers have been reported and then explained by continuum nonlinearities (e.g. Morishita et al. 2016). We have found (with our linear theory; but with the energy release rate replacing K) that a drop in fracture energy C with increasing speed is sufficient to account for the observed discontinuity in speed (Schapery 2021). Alternatively, 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 J-integral as the crack driving force; the craze zone in some plastics could be a candidate fitting this description.

Appendix-Initiation of crack growth
We continue to use Eq. (17), but now have a transient problem. Referring to Fig. 1, the crack tip in the initial unloaded state is located at the coordinate system's origin, X = 0. Thus, the left end of the FZ is at the point, Other points in the FZ are located at, xðX; tÞ ¼ aðtÞ À X [ 0 With this notation Eq. (17) becomes, vðX; tÞ ¼ ð1=2pÞ where v w is defined as an auxiliary displacement, v w ðX; sÞ Z aðsÞ 0 r f ½vðX 0 ; sÞF aðsÞ À X 0 aðtÞ À X The problem of predicting the time at which crack propagation starts is somewhat analogous to that of continuous growth. However, it is greatly simplified if the so-called quasi-elastic approximation (Schapery 1965(Schapery , 1967 or if the Dugdale model is applicable. The former requires that both the creep compliance and v w ðX; tÞ (or their logarithms) exhibit small curvature on the logarithmic time scale. This approximation for the present problem reduces Eq. (72) to vðX; tÞ ¼ C v ðtÞ 2p v w ðX; tÞ ð 74Þ which is the elastic solution, but with creep compliance in place of the elastic compliance. Thus, the fracture initiation time t c is the root of where C t is not necessarily the elastic fracture energy. The Dugdale model, under the assumption that K(t) is nondecreasing, provides a simple, exact solution for initiation time t c (Schapery 1975b), which is the root of where C tef ðt c Þ ¼ 1 is an effective (or secant) compliance. When KðtÞ is a step function of time applied at t = 0, this result obviously reduces to Eq. (75), which does not require Dugdale specialization. We can now estimate the time t c relative to the time taken for the crack to propagate the length of the FZ, a _ a , immediately after growth initiates. Comparing Eqs. (47) and (75), assuming K is constant after it is applied stepwise, and that C t % C and I R % 1, Referring to Fig. 6, after growth begins it takes three to four 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.