Crack growth in viscoelastic media with large strains: further results and validation of nonlinear theory for rubber

This paper is a continuation of two recent publications on crack growth in viscoelastic media. It provides further theoretical results for large strains that enable prediction of crack opening displacement for comparison with experimental data in the region of the singularity. In order to achieve good agreement with experiment it was necessary to account for far-field viscoelasticity. Additionally, it is found that with large deformation throughout the singularity, the deformation consists of simple shearing and stretching normal to the crack plane. Thus, there is no significant displacement parallel to the crack plane; such simplicity exists for materials that stiffen or soften at high strains if the stress obeys a power law in strain at high strains. This finding means that, despite the frame-dependence of the theory, there is no local rotation in the singularity to affect the stress.


Introduction
This paper is the third in a series of three on crack growth in viscoelastic media. The first (Schapery 2022a) is concerned primarily with the effect of various models of the failure zone (FZ); i.e. the very small crack-tip process zone where the material comes apart. Using linear theory for the continuum, it is shown that the primary effect of the different models is to impart a small shift along the log-speed axis.
In the second paper (Schapery 2022b), hereafter cited as (RSb), it is shown that the author's nonlinear theory for viscoelastic crack growth with large deformation (Schapery 1984), later cited as (RSc), is in very good agreement with experimental data on carbon black (CB) filled rubber (Morishita et al. 2016), hereafter cited as (MTU); these data were collected using so-called pure shear specimens. This nonlinear theory is not much more involved than linear theory, and is easily used in predictions and in comparing theory with experiment. Citation of other related literature is in (Schapery 2022a) and will not be repeated here.
The primary simplification in the nonlinear theory is that all viscoelastic effects are characterized by one linear viscoelastic creep compliance. This assumption is quite well-satisfied by some highly-deformable materials, such as rubber. A second simplification is that, in the singularity, the strain energy obeys a power law in strain.
One objective of this paper is to predict deformation in the singularity and the associated J v integral. As a second objective, we compare experiment and theory for crack opening displacement (COD).
The experimental data (MTU) use millimeters for length and seconds for time. Our notation is simplified by interpreting all lengths and times in this paper as nondimensionalized with these dimensions. Additionally, far-field refers to points far outside of the singularity. Mathcad was used for all calculations.
Although the present analysis is concerned with a study of the opening-mode singularity, and thus is two-dimensional, the far-field may involve threedimensional deformation. It is J v , as found from the far-field, that determines the stresses and strains in the two-dimensional singularity. In other words, J v in the nonlinear viscoelasticity theory plays the same role as stress intensity factor in linear viscoelasticity theory.

Theory
The basic equations used in the crack growth theory are given in this section. Figure 1 shows the shape of the elastic COD for the opening mode within the singularity, which is the only mode considered here. The solid line is for a power law nonlinearity with exponent N = 1.8 based on the experimental data in (MTU), while the dotted line is the COD shape for a linear material, N = 1. The CODs are arbitrarily matched at the left end.

Geometry
Cartesian and polar coordinate systems are shown in Figs. 1 and 2. The X-Y system is fixed in the continuum, while the x-y and r-h systems are attached to the crack tip. All three systems define material points in the undeformed body. Note that the direction of positive x is opposite to that in (RSb) for consistency with the polar coordinate system.

Constitutive equation
Details of the theory are in (RSc). The Piola stress r(whose transpose is called the Lagrangian stress tensor by Fung (1965)), satisfies where W is a strain energy density for unit initial volume,ru R is the pseudo displacement gradient tensor, and u R is the pseudo displacement vector, as defined by the convolution integral, Gðt À sÞ ou os ds ð2Þ in terms of the displacement vector u. Although other relaxation moduli could be used, we use here the shear relaxation modulus, G(t), and its (long-time) equilibrium value, G e . The stresses and displacements are defined in terms of coordinates in the undeformed body for each current crack tip location. The inverse of Eq. (2) is It is pointed out in (RSb) that Eqs. (1) and (2) have been shown to be a special case of constitutive theory based on nonequilibrium thermodynamics (Schapery 1999, Sect. 3) if the stress therein is interpreted as the Piola stress and strain as the displacement gradient.
Analysis of the viscoelastic continuum is reduced essentially to elastic analysis when Eq. (1) can be used. One may check the history aspect of this equation experimentally by expressing displacement gradient (''nominal strain'') in terms of the pseudo variable. The resulting stress-pseudo strain behavior will be (essentially) independent of the strain history if this model is valid. Undoubtedly, it will not be valid for all histories, such as those that involve repeated loading and unloading, producing in rubber the socalled Mullin's effect (e.g. Mullins 1969;Mai et al. 2021) unless microstructural effects (such as breakage of entanglement points and microcracking) are accounted for (Schapery 1982).

The boundary value problem
Assume the solution for viscoelastic stress (which is equal to the elastic stress, with or without growing cracks(RSc)) has been found, which we can express as where T is the specified boundary traction and U R is the specified pseudo boundary displacement; the latter is obtained from the given physical displacement using Eq. (2). Here, we formulate this elasticity problem by expressing the (pseudo) strain energy density W ¼ WðE R Þ as a function of the pseudo displacement gradient tensor through the pseudo Green's strain tensor, where One way rigid body rotation can affect the stresses is through global rotation (which affects U R ); when pseudo displacements are derived from displacements due to (time-dependent) rotation they may produce stresses because the convolution integral distorts the rotation transformation and produces deformation. (An example of the effect of rotation is in the Appendix.) Suppose, for example, a structural element undergoes a large amount of uniform global rotation and translation. Then, by using coordinates attached to the element to define the ''initial'' geometry in computing pseudo displacements, the effect of this motion is removed. This is, of course, standard practice in the analysis of beam, plate and shell elements subject to large motions.
However, if there is local, nonuniform timedependent rotation it affects E R , and thus affects the stress, except for short and long time elastic behavior, and for geometrically linear theory. Fortuitously, we find in Sect. 4.0, for a wide range of material behavior, elements in the singularity initially aligned with the (x-y) axes do not rotate.

Equilibrium equations in polar coordinates
The standard form of the equilibrium equations must be modified to account for the unsymmetrical Piola stresses. This is easily done using the principle of virtual work, Upon using Eq. (1) we obtain, where the displacement gradient is Integration-by-parts of the left side of Eq. (8), for arbitrary virtual displacements du R , yields the equilibrium equations, which reduce to the standard form when the stress tensor is symmetrical.

Power law nonlinearity
As before (RSb), we assume that the elastic strain energy density within the singularity is a homogeneous function of degree N ? 1 in the pseudo displacement gradient, where c is a scalar and N is a constant. It readily follow that the 1-D stress-pseudo strain ðr À e R Þ behavior has the form (RSb) where r 0 is a constant with stress dimensions. The homogeneity assumption excludes dependence on the displacement gradient tensor through Green's strains unless the nonlinear component in the latter is negligible or dominant. The latter is assumed here.

Elastic displacements and stresses
in the singularity using polar coordinates.
There is no need to distinguish between elastic and viscoelastic displacements in this section; therefore, the R superscript will be omitted to simplify the notation. The singularity is defined by the r-dependence of displacements (RSb), which corresponds to the stress variation

Method of analysis
An approximate elastic solution will be obtained using the virtual work equation for a continuum of unit radius and thickness, accounting for symmetry with respect to the x-axis and using polar coordinates (r, h), where W T is the total strain energy in the upper halfplane, while T and u are the traction and displacement vectors on the boundary r ¼ 1.
A power law generalization of the neo-Hookean model for the strain energy density is used, as expressed in terms of the first invariant of the deformation tensor, in which g N is a constant modulus that may be selected such that Eq. (16) agrees with stress-strain data at large tensile strains; it is one-half the shear modulus for a neo-Hookean material. The principal extension ratios, k 1 and k 2 , are simply related to the principle Green's strains,E 1 and E 2 , (Fung 1964), where we define E 1 to be the largest principal strain. The third extension ratio is small, and thus can be neglected with incompressibility and the assumption that k 2 1 ) 1; assuming plane stress, there is no constraint connecting k 1 and k 2 . Thus, Note that (-1) is to be neglected to preserve homogeneity. The total strain energy in the upper half-plane is where G N = g N x unit thickness. Although motivated by the neo-Hookean model, Eq. (19) is the most general form of a power law that can be used, given the requirements of plane stress, homogeneity, incompressibility and initial isotropy. (The second deformation invariant I 2 is not needed because it is homogenous in the principal strains only when it is equal to I 1 .) Only terms that are quadratic in the displacement gradient ru are retained in Green's (pseudo) strains to preserve the homogeneity of the strain energy density, as noted earlier. These strains are where r is the gradient operator in Eq. (9). Thus, 3.2 The exact large strain solution for N = 1.
As a guide to the development of an approximate solution, it is possible to obtain an exact solution for N = 1. This is done most easily starting with Cartesian coordinates, as shown later in Sect. 4.1. In polar coordinates the displacements for N = 1 are found as and the Piola stresses are r r ¼ G N u 1 2 ffiffi r p cosð0:5hÞ À cosð1:5hÞ ½ ; where is the lower surface h-displacement at r = 1.
3.3 Approximate solution for N 6 ¼ 1 Displacements are assumed using free coefficients, which are then found from Eq. (15). Guided by the linear elastic solution's dependence on h (Williams 1957) as well as the exact, large strain solution for N = 1, Eq. (22), we assume u r ¼ r N Nþ1 0:5a 1 cosð0:5hÞ À 0:5a 2 cosð1:5hÞ þ a 3 cosð2:5hÞ þ a 4 cosð3:5hÞ f g where the dependence on r is exact, as found from the equilibrium equations. The terms following 1:5h are not in the linear solution, and are shown later to have very small coefficients. We should add that Eq. (15) does not converge when N = 1.
The surface displacement at r = 1, Eq. (24), is given as unity in solving for the coefficients.
All resulting displacement coefficients are proportional to u 1 when it is not unity as a consequence of the power law nonlinearity. Also, because G N appears as a factor on both sides of Eq. (15), its value has no effect on the value of the coefficients.
The radial displacement vanishes along the crack surface in the linear solution. This was found to be true for the nonlinear solution as well by adding terms with integer values of theta in the cosine and finding that this displacement was essentially zero. However, the radial strain at the crack surface is not zero because of the large deformation.
Using Eq. (25), the displacement gradient tensor as well as Green's strain tensor and its principle values can be found, after which the strain energy, Eq. (19), follows. The tractions in Eq. (15) may be expressed using the Piola stresses r, Eq. (1), and unit normal n (for an area element defined in the undeformed state), With index notation (i, j = 1,2) in Cartesian coordinates, Recalling that r is a transposed tensor, the first index on stress is the traction direction and the second is the surface normal direction, while recognizing that this stress tensor is not generally symmetric with large deformations.
With these considerations, Eq. (15) provides the coefficients shown in Table 1 when u 1 ¼ 1. The column for N = 1 is not, of course, for linear theory because of the large strains. Coefficients with magnitudes less than 0.001 are shown as zero. Table 1 The approximate solution is quite accurate in that the equilibrium equations are found to have a very small maximum error, relative to the maximum value of stress; e.g., on the order of 1x10 À3 for N = 1.8 and 5x10 À3 for N = 4. The traction and displacement boundary conditions at h ¼ AEp are satisfied exactly. Figure 3a shows the displacements, while Fig. 3b gives the largest principal extension ratio k 1 .The Piola stresses for G N ¼ 1 are in Fig. 4. The exponents N = 1 and N = 1.8 are used for these figures.

Some results based on
These graphs of elastic solutions are for a specified value of u 1 = 10, which is approximately the pseudo value for both materials (discussed later) used in the experimental work. Also, all graphs in polar coordinates use r = 1 because the solutions are proportional to in-common functions of r, which can be easily introduced.
We find 0 E 2 \10 À5 which means k 3 ' 1=k 1 due to incompressibility. Note that the largest value of k 1 is at the crack face. Indeed, one can easily confirm that it is the value one finds by simply calculating the extension ratio due to the stretching of the surface h ¼ AEp imposed by the crack opening with u r ðAEp; rÞ ¼ 0.
Also shown in Fig. 3b is the local rotation (negative for clockwise rotation) of elements with pure deformation, as rotation is normally defined. Biot (1965) develops it for large strains and gives it in Cartesian coordinates as  Note the weak dependence of the coefficients on N, even with a considerable amount of strain softening and hardening. The primary effect of N is in the r-and J v -exponents Note that the rotation vanishes if the shear gradients (i.e. gradient tensor components) are equal. For a general planar orthogonal system, where the subscript z denotes the out-of-plane component, Eq. (28) becomes In polar coordinates, Eq. (29) is (after converting it to degrees), When N = 1, it reduces to It is plotted in Fig. 3b, and seen to be quite large over a wide hrange.
The Piola stresses in Fig. 4 use the same given quantities as in Fig. 3. As expected, there is not much difference for the two values of N when r = 1.

Relating the COD to J v
The nonlinear elastic displacement distributions are found by substituting the coefficients in Table 1 into Eq. (25) and multiplying by the given COD/2 at r = 1. However, this COD is not usually a known quantity. Rather, it is the viscoelastic J-integral, J v , that is known or specified.
The relationship between J v and COD R is found by starting with the definition of the two-dimensional version of J v for a body of unit thickness, where T is the traction vector and S is the path of counterclockwise integration surrounding the crack tip, which starts on one crack face and ends on the opposite face, for a crack along the x axis. As above, the superscript R is omitted in the following elasticity analysis to simplify the notation. The integral is independent of path; we select the path to be circular with unit radius, The tractions and displacements are written in the Cartesian coordinate system (x,y) with displacements ðu x ; u y Þ. Thus, with the definitions, c cosðhÞ; s sinðhÞ ð 35Þ the relationships are, The Piola stresses in the (x,y) system are found from those in the ðr; hÞ system using Eq. (1) and a standard tensor transformation but with subscripts on the shear stresses interchanged, r x ¼ r r c 2 þ r h s 2 À ðr rh þ r hr Þcs r y ¼ r r s 2 þ r h c 2 þ ðr rh þ r hr Þcs r xy ¼ ðr r À r h Þcs þ r rh c 2 À r hr s 2 r yx ¼ ðr r À r h Þcs À r rh s 2 þ r hr c 2 together with the relationship between coordinate systems, Given the displacements Eq. (25) and strain energy density Eq. (18), the Piola stresses in the ðr; hÞ system may be found in terms of u 1 . In turn, Eqs. (34)-(39) provide Jv.
We next define a coefficient whose value is the COD at r = 1 when J v =G N ¼ 1.This coefficient depends only on N. Its values are given in Table 1, where it is seen to be a very weak function of N. As a comparison, when N = 1 and the deformation is small, we find C J = 4= ffiffiffiffiffi ffi 3p p ' 1:30. Thus, with large deformation for all values of J v =G N and r, and reintroducing the superscript R, 4 Elastic displacements and stresses in the singularity using Cartesian coordinates.
To predict crack growth in the x-direction it is necessary to use a Cartesian system. Equations (37)-(39) convert the displacements and stresses to those in the Cartesian system. Figure 5a shows graphs of representative displacements for N = 1.8. It is seen that u x is three orders of magnitude smaller than u y . The same relationship is found for the gradients; thus only the gradients of u y are drawn in Fig. 5b and c. The u x and its gradients were found to be negligible for the entire N-range in Table 1.
The principal direction is in Fig. 5d, where it is seen that it is closely aligned with the crack face when y = 0.1, except for small x; this is consistent with Fig. 3b in which the principal extension ratio is greatest at h ¼ p. (In a study of the linear elastic solution (unpublished) it is found that u x and its gradients are comparable in magnitude to u y and its gradients, although mostly smaller. In this linear case locally small rotation does not affect the stresses.)

An exact solution
Let us look for an exact solution when u x 0: Using this condition in Eq. (20) we find These strains produce the principal strains, The strain energy density is found from Eq. (18), The Piola stresses in the x-direction obviously vanish, For N = 1 the y-stresses greatly simplify, They must satisfy the y-equilibrium equation, while the x-equation is automatically satisfied. Upon use of Eq. (46) in Eq. (47) we obtain Laplace's equation for u y , which is solved most easily in polar coordinates by assuming Upon substitution in Eq. (48) we find that the relevant solution is where A is an arbitrary constant at this point. Finally, from Eq. (37) and trigonometric identities, we obtain Eq. (22) after selecting A to satisfy the displacement condition at h ¼ p; r ¼ 1. The stresses in Eq. (23) then follow from Eqs.(45), (46) and the inverse of Eq. (38).

Viscoelastic analysis
In this section we use the superscript R on the displacements and its gradients, recognizing that all of Sects. 3 and 4 apply here when this modification is made. The stress and J v integral notation is unchanged. A very important result from the elasticity analysis is that u R x ffi 0, and thus the corresponding physical displacement is negligible. This means that the local deformation consists primarily of only normal stretching and shearing in the y-direction, as depicted in Fig. 6.
Obviously, this element does not rotate. The rotation in Fig. 3b is not relevant. The deformation in the entire singularity can be essentially accounted for using only the two gradients ou R y =ox and ou R y =oy in the strain energy. This finding helps to show why the theory and experiment in (RSb) are in such good agreement with very large strains.
If ou R x =oy was not zero, the smallest (in magnitude) of the two shear gradients (after converting it to the physical gradients using Eq. (3)) could be used in estimating the effect of rotation on the stresses.
In one comparison of theory to experiment in (RSb), the shape of the COD in the singularity was predicted and found to agree with experiment for rubber filled with two different volumn fractions of carbon black (CB), 0.09 and 0.17.
However, the COD itself was not predicted because the relationship between J v and COD R was not available; but it is now in Eq. (41). Here we shall provide details on the prediction of the COD R for the 0.09 and 0.17 CB rubber.The CODs will then be compared to experimental values.
Additionally, the viscoelastic deformation field will be shown at selected distances from the crack plane.

Linear viscoelastic creep compliance
The normalized master curves of shear creep compliances are in Figs.7a and 8a using.logarithmic (base 10) coordinates. They are plotted for a reference temperature of 25 0 C, which is the temperature used in the pertinent experiments. The equilibrium compliances are Cs e ¼ 0:88 ðMPaÞ À1 for 0:09CB and Cs e ¼ 0:60 ðMPaÞ À1 for 0:17CB: The compliances were developed from the master curves of the real part of the dynamic shear modulus by using a cubic spline fit to the data. Each master  curve in (MTU) has a data spread of approximately 12%, while the fit was to its centerline. Figure 9 shows the nominal stress-strain data at 25 0 C for both materials from (RSb). These quantities are the Piola axial stress r and the strain e ou y =oy. In (RSb) it was assumed that the material was essentially in an equilibrium (elastic) state; this was believed (incorrectly) to be acceptable because of the low rate of stress relaxation during the test period. The stress-strain curves (from a so-called pure shear test) were generated using the same strain rate as in the crack growth study (Urayama 2022). Specifically, specimens were loaded at a displacement rate of 1 mm/s, with an initial specimen height of H = 20 mm. In crack growth tests, when the desired constant displacement was reached, it took only 1-2 s while an edge crack was cut, followed immediately by crack growth (Urayama 2022). The time scale of each test was not long enough to neglect viscoelastic effects in the stress-strain behavior, especially in predicting COD, as done here. This condition necessitates a modification of the data reduction used previously, this time using pseudo displacements instead of displacements. This correction is needed because, according to theory, the value of J v is equal to the energy release rate based on pseudo displacement, not displacement itself. This correction is done by first modifying the stress-strain curve using pseudo strain for the abscissa. According to Eq. (2),

Stress-strain data and determination of J v
which requires the relaxation modulus; note that e R is proportional to the linear viscoelastic stress. This modulus may be easily constructed from the creep compliance using the well-known equation (Ferry J 1980), where is the local log-log slope. It is very accurate for functions spread smoothly over many decades (Schapery and Park 1999). The relaxation moduli are drawn in Figs. 7b and 8b. For the purpose of numerical integration they were fit using a cubic spline. The factor G N is also needed. The strain energy Eq. (16), for the test conditions k 2 1 ) 1;k 2 2 ¼ 1; k 2 3 ( 1; e¼ ou y oy becomes, in terms of pseudo strain, Thus, The numerical values of various quantities for CB = 0.09 and CB = 0.17 will have 9 and 7 subscripts, respectively. With time as the common parameter, the stressstrain curve is converted to stress-pseudo strain. When e ¼ 3, the pseudo strains are e R 9 ¼ 4:90 and e R 7 ¼ 6.47. The complete stress-pseudo strain curves are the solid lines in Fig. 10. Matching Eq. (57) to them at the e ¼ 3 state gives, in MPa-mm(= kJ/m^2) units, These values are used in Eq. (57) to draw the dotted lines in Figs.10 and 11.
The solid lines in Fig. 11 are the stress-pseudo strain curves having a maximum value at the pseudo strain used to photograph the COD. That the end point touches the power law means that the far-field strain is high enough to be in the power law regime, assuring that N = 1.8 for all strains in the singularity.
For the pure shear test, the value of J v is equal to the specimen height (20 mm) multiplied by the area under the solid lines in Fig. 11. This yields, in MPa-mm units, which are the J v values when the CODs were photographed.

Prediction of the elastic COD
The elastic COD in the singularity is given by Eq. (41), From Eqs. (58) -(60), which provides the values at x = -1, to be used in the elastic displacements, Eq. (25), which are now designated as u R r and u R h .

Prediction of the viscoelastic COD
Equation (3) provides the viscoelastic COD, According to Fig. 1, xðX; tÞ ¼ X À aðtÞ, with X fixed in location, then dt ¼ Àdx= _ a. Assuming crack speed _ a is constant, at least during the time it takes for the crack to propagate the length of the singularity (RSc), and selecting X = 0, without loss in actual generality, COD becomes As done in (Schapery 1975), we change the integration variable to and find where the lower limit is LL ¼ 1 and a weight function w(z, p) is defined, Also Now, numerical integration of Eq. (67) does not converge unless the lower limit is finite. We find the factor 10 z in the weight function enables the lower limit to be changed to -3 with negligible error. Additionally, the narrow width of the weight function leads to (RSb) for x ! 0.The factor s f is a quantity that depends on N and the local log-log slope of the creep compliance; for the materials used here s f % 2 for the relevant range of x= _ a. We can now predict the COD for the two materials. The applicable set of parameters is _ a 9 ¼ 2510 mm=s; J v9 ¼122 kJ=m 2 ; G N9 ¼ 0.0979 _ a 7 ¼ 531 mm=s; J v7 ¼183 kJ=m 2 ; G N7 ¼ 0.134 ð71Þ and, for both materials, With these numbers, we find at the end points x 9 ¼ À6 and x 7 ¼ À4, using Eqs. (61) and (67), while approximate Eq. (70) yields The experimental CODs are, Thus, the theoretical COD is 7.5% greater than the experimental value for CB = 0.09, while it is 4% less than the measured value for CB = 0.17. Recognizing that there are significant specimen-to-specimen differences in stress-strain behavior (Urayama 2022), such as the 12% mentioned earlier for the linear viscoelastic modulus, it may be concluded that the theory satisfactorily agrees with the data.
At the end points, showing that viscoelasticity greatly suppresses the COD compared to the pseudo value (which uses the equilibrium value, C ¼ 1) Moreover, this strong effect of viscoelasticity means that the neglect of first order terms in Green's strains (in terms of pseudo gradients) produces considerably less error than the actual strains would in the singularity.

Effect of viscoelasticity on J v
In (RSb) the J v integral was calculated directly from the data in Fig. 9, using whereê is the applied strain used for the COD pictures. In this case, which are considerably less than the values in Eq. (59). The data analysis of crack growth in (RSb) made use of only the ratio J ¼ J v =C, where C is the intrinsic fracture energy, as found by shifting data. Consequently, the fracture energy values reported in (RSb) should be increased by the factors These factors have been found to be essentially the same (within 5%) for all crack speed data points in Fig. 7b,c in (RSb); even if this 5% variation were accounted for, there would be no observable effect on the log-log plots.
Thus, the only significant corrections to (RSb) are the values of C and the associated change in value of the failure zone parameters, r m and v 0 , such that r m v 0 ¼ C. Also, the secant stiffness r m =v 0 doesn't change; thus, r m and v 0 $ ffiffiffi ffi C p . Table 2 gives the corrected values for low (L) and high (H) crack speeds.
Also listed in Table 2 are values for k D , the nondimensionalized version of C J , as introduced in (RSb). We find, where, for plane stress C e ¼ 4Cs e =3; thus C e9 ¼ 1:17 and C e7 ¼ 0.80: For a linear material k D ¼ 1.
Crack growth for a third material, CB = 0.05, was analyzed in (RSb), although the COD was not available. The correction factor for its fracture energy has been found as 1.40; corrected critical state values are in Table 2.
The increase in J v due to far-field viscoelasticity, which is rigorously established in (RSc), is clearly needed if J v is to be proportional to K 2 when the theory is reduced to linear theory, where K is the stress intensity factor. Additionally, if the stress-strain data are not corrected for far-field viscoelasticity, as in (RSb), the COD predictions are, which are much smaller than the experimental values, Eq. (75).

Prediction of the viscoelastic deformation field
Given a constant crack speed for all time we may easily predict the viscoelastic displacement and gradient fields. In contrast to the COD, it is necessary to consider all past time; thus, Use x ¼ À _ at, and make a further change of integration variable like Eq. (66), but now for positive and negative x, With these two changes Eq. (82) becomes In which lower LL and upper UL limits must be large, but finite, for convergence of the numerical integration. The weight function is for x \0; In the following examples u R 1 ¼ 10 is used, as before. The viscoelastic displacement and gradients are drawn in Figs. 12-14 for CB = 0.09 using values in Eq. 71. The weight functions w ux and w uy are predicted by simply replacing u R y with ou R y =ox and ou R y =oy, respectively, in Eq. (85). The weight functions are drawn in Figs. 12-14 for selected (x, y); but they are similar for all (x, y). We can use LL = 3 and UL = 6.6 because there is no significant contribution to the integrals beyond these limits. The numerical integration for ou y =oy does not converge for UL [ 6.6. The farther the location from the crack tip, the larger the limits must be, especially UL. These weight functions are much broader than that for the COD, Eq. (68), which is approximately one decade. As a result, approximations. such as Eq. (70), will not be nearly as accurate.

Comments on frame-dependence
If a 3-D (2-D) deformation field is defined essentially by normal gradients and, at most, three (one) shear gradients, each in a different plane, the issue of framedependence does not apply. In this case the strain energy may be expressed as a function of, at most, the six pseudo deformation gradients that completely define the deformation. It is clear that any effect of rotation should be based on rotation of an element's edge, not with the rotation of an element suffering pure deformation-the standard definition of rotation. Of course, if co-planar shear gradients are essentially equal there is no significant rotation in that plane.
Generally, it is the smallest shearing component of the time-dependent displacement gradient tensor in each plane, and at each point, which defines the amount of rotation that may affect the stress prediction (neglecting interaction between local elements). This fact provides a way to estimate the effect of rotation on the stress from the deformation solution; i.e. the pseudo rotation would be converted to the physical rotation using the transformation in Eq. (3). An example of this estimate is in the Appendix.

Comments on plane strain
In this section the notation in (RSb) is used.
In order to determine the length of the failure zone (FZ) a from experimental data on crack speed, the dimensionless constant k J defined in (RSb, Eq. (17)), k J 4a J v is needed; the overbar denotes normalized value and k J is unity in linear theory. Here k J will be estimated for N = 1.8.
Very close to the crack tip, on the scale of the failure zone (FZ) length a (assumed to be small compared to the specimen thickness), the local deformation due the crack is complex and threedimensional. In the far field where k 2 = 1 and plane stress exists, the incompressibility condition gives k 3 = 1/k 1 , which is a small quantity for the test conditions. We assume k 3 is not small close to the crack tip; but instead k 3 & 1 for the material affected by the failure zone. The COD produces this increase in the local specimen thickness. This approximate plane strain condition, incompressibility and k 1 ) 1; lead to Eq. (18) but with k 2 = 1/k 1 & 0. The latter result implies E 2 & -0.5. However, to preserve homogeneity (which may or may not be true close to the crack tip), we must set E 2 = 0 in the virtual work analysis.
(Although E 2 was found to be very small for plane stress, its presence had a significant effect on the virtual work.) With this one change in the virtual work, the constants for N = 1.8 are: a 1 = 1.202, a 2 = 1.107, a 3 = -0.046, a 4 = 0, b 1 = -.005, b 2 = 0.049, b 3 = 0.051, C J = 1.42. Note that the value of C J is reduced to 0.90 of that for plane stress. This reduction is only a little less than ffiffiffiffiffiffiffiffiffiffiffiffi ffi 1 À m 2 p ¼ 0:87 as predicted from linear theory. With these constants the length of the failure zone can be estimated. The method consists of equating the load supported by the failure zone, r m a, to the load produced by the singular stress acting over the length a. This does not account for the change in singular stress due to the FZ; the change is approximately accounted for by using a reduction factor for a from linear theory; i.e. the ratio of the exact a (e.g. Schapery 2022b) to that from the uncorrected singular stress, p 2 /16 = 0.62.
Theses considerations, combined with values in Table 1, produce the following results for the two materials in the low(L) and high speed (H) ranges: CB = 0.09: k JL = 3.44, k JH = 2.94 CB = 0.17 k JL = 1.56, k JH = 1.52 The normalized FZ length is somewhat smaller for plane stress with the same J v , CB = 0.09: k JL = 2.74, k JH = 2.34 CB = 0.17: k JL = 1.25, k JH = 1.22 These results should be considered tentative in view of the idealizations used. A nonlinear, 3-D finite element analysis would be helpful in checking the estimated values of k J.
We found u x , ou x =ox and ou x =oy with E 2 = 0 are not quite as small as those for plane stress; but they are small enough that the predicted mechanical state of the viscoelastic rubber is not significantly affected by local rotation.

Concluding remarks
We have shown that the theory (RSc) is in good agreement with experimental data on rubber filled with carbon black under very large deformation, previously for crack speed and shape of the COD (RSb), and now for the COD itself. It was found that viscoelasticity of the far-field has a significant effect on the COD. Without accounting for it, the predicted COD is in considerable error; in one case it is less than half of the measured value.
Local material rotation in the singularity is essentially zero for the entire range of exponents in Table 1, 0:1 N 10 (in unpublished studies).
The simplicity of the model facilitates its use with microstructural changes, such as microcracking. To the author's knowledge, this extension has been experimentally validated only in the globally geometrically linear range (e.g. Ha and Schapery 1998), although the microstructure may have large strains, as in (Schapery 1982).
All of the complexities of an application, if any, appear in only the elasticity analysis. Considering the finite element method, for example, it would be needed only for an elasticity analysis if the present theory is applicable; viscoelastic effects may be introduced through pre-and postprocessing. This fact greatly reduces processing time when the theory is applicable.

Appendix
An example is given to illustrate the effect of local element rotation on the pseudo deformation which, in turn, affects the stresses. A rigid, time-dependent, counter-clockwise rotation x of a square element produces the displacements, The associated displacement gradients are Green's strains are, of course, unaffected by the rotation for an elastic material. However, to predict stresses in the viscoelastic material it is necessary to replace the gradients by pseudo gradients, using Eq. (2). In order to illustrate how rotation affects pseudo strains using a realistic gradient history for crack growth, we specify a pure rotation (no physical strains) due to a fraction of the gradient ou y =ox in Fig. 13b for y = 0.2 and CB = 0.09. For example, using a given fraction a, the angle of rotation is, which then provides the other gradients through Eq. (87). The resulting pseudo gradients then follow from Eq. (2). For example, The Green's strains in Eq. (6), retaining both linear and quadratic terms, then follow; because the gradient tensor is anti-symmetric and the diagonal elements are equal, Also, the normal strains are positive and even functions of the rotation. Thus, the effect of rotation produces straining analogous to that from thermal expansion.
The engineering pseudo strain (extension ratio-1) of this uniform expansion is, Given the relaxation modulus in Fig. 7b, the predicted engineering pseudo strains are drawn in Fig. 15 for three values of a.The associated rotation angles are in Fig. 16.
Of particular interest is the (x,y) location where the maximum strain due to rotation alone occurs. Figure 3b shows the maximum extension ratio at the radius r = 1 without rotation. By reducing the radius to correspond to the (x,y) at which the maximum strain due to rotation occurs, we find k 1 À 1 ' 8.
We thus estimate that if the gradient ou R x =oy is not zero, but instead corresponds to the angle shown in Fig. 16, the effect of rotation probably can be neglected if the rotation is less than ten degrees.
Author contributions R.A. Schapery is the sole author.

Declarations
Competing interests The authors declare no competing interests.