## Abstract

We present a novel approach to analyse the fracture of fibre-reinforced composites. Experimental results on mode I fracture of glass fibre and carbon fibre unidirectional laminates presented here and published by others in the open literature formed the basis for the analytical and numerical results presented. When details of the external loading rate are explicitly accounted for, a new picture of fracture emerges, which encompasses the possibility for non-smooth crack growth and the necessity to relax the use of a critical strain energy release rate as a criterion for crack advancement. Results predicted by adopting the analytical model presented here are seen to capture a wide variety of fracture responses that have been observed previously.

## 1. Introduction

Layered materials are ubiquitous in nature, spanning a wide range of length scales as observed in human muscle to sedimentary rocks. Delamination is a commonly observed mode of failure in synthetic layered materials. Characterizing the delamination resistance of continuous fibre polymer matrix composite laminates (PMCLs) is of importance in assessing the structural integrity of components made of PMCLs. A popular contemporary approach in achieving this is through a combination of standardized coupon level tests in conjunction with finite-element (FE)-based numerical analyses of the test specimen/s that use cohesive zone modelling strategies [1–4]. With this approach, the modes I, II and III fracture toughness, *G*_{Ic},*G*_{IIc} and *G*_{IIIc}, and associated cohesive strengths, *σ*_{Ic},*τ*_{IIc},*τ*_{IIIc}, respectively, are obtained and used in subsequent structural analyses, as input data, in characterizing a composite structure and its ability to resist delamination growth. In those instances where the loading at the delamination tip involves the presence of more than one mode of crack growth, a suitable mixed-mode crack growth law is also needed.

While a substantial number of delamination crack-growth studies have focused on smooth quasi-static crack growth, regardless of the externally applied loading rate, there are many instances in which crack growth can occur in a non-smooth manner, where periods of stable smooth growth can occur bounded by regions of unstable fast crack growth even though the external loading rate is quasi-static. In these instances, there is a need to examine how rate dependence affects fracture toughness, cohesive strengths and the associated mixed-mode crack growth laws. Early investigations pertaining to examining the rate dependence of *G*_{Ic} were conducted by Aliyu & Daniel [5], while a conflicting result was later reported by Yaniv & Daniel [6]. Maikuma *et al.* [7] examined rate dependence of *G*_{IIc} and reported a decrease in value with an increase in loading rate. A comprehensive study, using double cantilever beam (DCB) specimens, on the effects of loading rate on mode I fracture of unidirectional carbon fibre composites was reported by Kusaka *et al.* [8]. They found that *G*_{Ic} decreased in a step-wise fashion with a transition region in which *G*_{Ic} showed a strong dependence on rate; however, above and below this region *G*_{Ic} was found to be fairly constant. Below the transition region, fracture was found to be unstable and non-smooth exhibiting stick-slip response, while above the transition region stable fracture was observed. The experimental findings were explained by incorporating rate dependence in *G*_{Ic} and kinetic energy contributions of the DCB loading arms. This finding was in contrast to that observed in Blackman *et al.* [9], who observed stable crack growth at low loading rates in a carbon fibre/PEEK composite that transitioned to unstable stick-slip behaviour as the loading rate increased. Stick-slip behaviour, but with a stochastic ductile to brittle transition, is also reported in the studies by Sun *et al.* [10], who examined fracture of adhesively bonded metallic specimens.

In the present work, the mode I fracture of glass fibre unidirectional composite (GC) samples and carbon fibre unidirectional composite (CC) samples is studied under quasi-static loading rates, using the ASTM standard [11] DCB test specimens. It is found that crack growth, under the examined quasi-static loading rates, exhibits both smooth and non-smooth responses, rendering the interpretation of an initiation toughness (*G*^{i}_{Ic}) problematic. Moreover, while an *R*-curve-like response is observed in the CC laminates, a plateau toughness, *G*^{p}_{Ic}, which is called for in the ASTM standard is not attained. In GC laminates, the observed test results are difficult to interpret with respect to identifying *G*^{i}_{Ic} and *G*^{p}_{Ic}. Non-smooth response has also been observed by Kusaka *et al.* [8], who explained the cause by appealing to rate dependence of *G*^{i}_{Ic}. Blackman *et al.* [9] also found non-smooth, stick-slip crack growth in a carbon fibre/PEEK composite as the loading rate was increased, suggesting that a definition of a critical load to back out a value of *G*^{i}_{Ic} would be problematic.

The present paper is organized as follows. Details of the experimental investigation including experimental results and observations are presented first. This is followed by an analysis of the test data using a cohesive zone FE model, while at the same time attempting to extract data according to the ASTM D 5528 standard. The numerical results, obtained using the cohesive zone FE model, are also used to illustrate the difficulty in defining critical fracture initiation, *G*^{i}_{Ic}, and propagation, *G*^{p}_{Ic}, parameters. Subsequently, the non-smooth crack growth response is explained by adopting a new analytical formulation of the DCB fracture problem that uses an alternative set of hypotheses associated with examining rate of crack growth energetics. A discussion of the results obtained through the different approaches is followed by concluding remarks to complete the paper.

## 2. Experiment

### (a) Quasi-static double cantilever beam test

Displacement-controlled DCB experiments were conducted on GC and CC unidirectional laminate coupons. The geometry of the sample is shown in figure 1, and table 1 lists the dimensions that follow the ASTM D 5528 standard. The initial crack was introduced using a Teflon film, inserted during manufacturing of the laminates. Steel blocks with a transverse through hole at the centre were bonded to the free end of the laminates for pin joint load transfer. The pin joint was lubricated using Teflon lubricant to reduce frictional loss. Clear markings on the specimen side surfaces were used to track crack growth. The modulus, *E*_{11}, in the fibre direction is 11.5 GPa for the GC and 117 GPa for the CC, while the major Poisson's ratio is 0.3 for the GC and 0.29 for the CC.

The tests were conducted on an Instron 4201 universal test machine with crosshead displacement rate of 5 mm min^{−1} for GC and 1 mm min^{−1} for CC. Load was measured continuously using a 1000 N high accuracy tension/compression load cell linked up to a data acquisition system. A high-resolution SLR camera, time synchronized with loading, was used to capture crack zone images with a framing rate of 1 frame per second. The images were then analysed manually, using a linear pixel measuring software calibrated against a reference grid pattern marked on a typical specimen. This method gives crack length in the time domain, which can be converted to plot crack length against load and load point displacement. After the test, each sample was checked for through the width variation and was found to conform to ASTM D 5528 standard [11]. Each test was terminated when the crack extended for approximately 30 mm from the initial position. Seven samples of each material were tested to ensure repeatability in test data.

### (b) Observations and results

Both GC and CC showed repeatable, linear behaviour up to the point of crack initiation with continued loading. GC specimens showed a higher compliance when compared with the CC specimens. Upon unloading, both types of specimens exhibited negligible permanent deformation. The mild permanent set is caused by some matrix micro-cracking and fibre bridging, the latter being the case for GC laminates. The majority of the energy dissipation can be attributed to the mode I crack growth, and thus effects of micro-cracks in the DCB arms are not considered in this study. Significant fibre bridging was observed in the GC laminates, whereas negligible fibre bridging was seen in the CC laminates, as shown in figure 2.

Figure 3 shows the normalized load–displacement response and figure 4 shows the crack growth versus the crosshead displacement for GC and CC laminates. Results are normalized to a value of *P*_{0}=120 N for ease of comparison. As seen in the crack extension plots (in figure 4), the failure history is quite different for the two different sets of coupons. In GC laminates, non-smooth type crack growth is observed. This non-smooth crack growth is also reflected in the force displacement plot, where each jump is associated with a sharp drop in force. This was an unexpected phenomenon because DCB tests by design are expected to result in smooth crack (no crack jumps) growth under quasi-static loading (0.5 mm s^{−1} to 5 mm s^{−1} as defined in ASTM standard [11]). Furthermore, others have reported smooth crack growth with *R*-curve response for GC laminates [12]. The fact that there is significant fibre bridging should facilitate smooth crack growth as reported by Morais & Pereira [13].

Kusaka *et al*. [8] have reported non-smooth crack growth in brittle carbon epoxy unidirectional laminates at loading rates of 0.01 mm s^{−1} which is below the rate used for defining quasi-static loading rate. However, they also demonstrated smooth crack growth with the same material system at higher loading rates of 5 mm s^{−1}. The reverse trend is reported in Blackman *et al.* [9], who show a transition from smooth to stick-slip behaviour as the loading rate increases. In CC laminates, even though the crack growth is smooth, the force continued to increase *even after initial crack growth* up to a certain point after which a gradual decrease in force was observed in the present experimental results. This result is similar to the results reported by Morais & Pereira [13], where an increase in fracture energy is seen with smooth crack progression. The key differentiator is the significant fibre bridging observed by Morais & Pereira, which is not present in the CC laminates studied here.

### (c) Energy release rate

Fracture energy was calculated using Modified Beam Theory (MBT), Compliance Calibration Method (CCM) and Modified Compliance Calibration (MCC) as specified by ASTM D 5528 [11]. The expressions for the strain energy release rate in the DCB test (*G*_{I}), based on MBT, CCM and MCC are
2.1
2.2
and
2.3where *P* is the measured load, *w* is the prescribed displacement, *b* is the width of the specimen, *a* is the measured crack length, *C* is the compliance, *κ* is the slope of compliance calibration curve [11], *A*_{1} is the slope of modified compliance calibration curve [11] and *h* is the height of the specimen. Figure 5 shows a typical strain energy release rate with respect to crack increment. Strain energy release rate values obtained from MBT are used for simulation and comparison purposes because MBT gives the upper bounds for energy release rates.

Figure 6 shows the *G*_{I} distribution with crack growth for GC and CC based on MBT calculation method. A cubic curve fit is used to obtain the *R*-curve in each case. In the case of GC laminates, an oscillatory pattern is seen for *G*_{I} with crack extension, which reflects non-smooth crack growth. The fitted *R*-curve of CC, on the other hand, conforms to a typical *R*-curve described in ASTM D 5528 [11]. The pattern observed shows values ranging from a minimum of 1.47 N mm^{−1} to a maximum of 2.2 N mm^{−1} and a mean value of 1.8 N mm^{−1}. In CC laminates, smooth crack growth is observed; hence, it is expected that the *R*-curve in the case of CC should also follow the one described in the ASTM standard [11]. However, the *R*-curve shows a monotonic increase with increase in crack extension. Similar results are also seen in experiments conducted by Morais & Pereira [13] for carbon-epoxy laminates with high fibre bridging. The increase in *G*_{I}, in the present study for CC laminates, is almost two times the initiation energy release rate. The *G*_{I} value ranged from 0.25 N mm^{−1} to 0.75 N mm^{−1} with a mean of 0.5 N mm^{−1}.

### (d) Experimental results: summary and discussion

In the previous section, two different crack growth histories are described for DCB tests conducted on GC and CC laminates corresponding to external loading regimes classified as quasi-static. The results show a mixture of smooth and non-smooth crack growth responses. Related experimental findings reported by others in the literature [8,9,13] also show a mixture of smooth and non-smooth crack growth histories. This poses a problem when critical energy release rates from coupon DCB tests are needed for predicting the response and failure of larger structures made of the same material or for comparisons between different material systems. For a typical DCB simulation, the critical value of *G*_{I} for initiation and propagation is identified from the experimental *R*-curve. As specified in ASTM D 5528, the point of nonlinearity in the *R*-curve is recommended for failure initiation, *G*^{i}_{Ic}, and the stabilized value is recommended for propagation, *G*^{p}_{Ic}. The pre-condition to obtaining the critical parameters from a DCB test is that crack propagation should be smooth; i.e. the crack growth rate with loading is constant. This definition of critical parameters can pose a significant constraint on different material systems because there is no clear initiation or propagation value that can be inferred from the experimental results, even with smooth crack growth. In the case of GC, *G*_{I} values have a large scatter and do not follow the typical *R*-curve behaviour given in ASTM standard [11]. In the case of CC laminates, there is a steady increase in *G*_{I}. Even though the initial point of nonlinearity can be used for *G*^{i}_{Ic}, the crack propagation *G*^{p}_{Ic} cannot be determined because there is no stabilization of the *R*-curve. Similar difficulties were also faced by Morais & Pereira [13] and Blackman *et al.* [9].

## 3. Numerical analysis

From the experimental results, the question as to the critical value of *G*_{I} that should be used for simulation and prediction of mode I fracture in the given material system has now to be posed. To answer this question and investigate the variance of *G*_{I} and its effect on the predictive capability, two approaches were used for the numerical simulation of the DCB test results. The first was to assume a fixed value of *G*_{Ic} along the length of the crack path (figure 7*a*), and the second was to vary the *G*_{Ic} value along the length of the crack path based on the observed experimental values (figure 7*b*). The latter was done to illustrate the limitations of using critical toughness values.

The numerical simulation of the DCB fracture experiment was conducted using the FE method in conjunction with the discrete cohesive zone method (DCZM). Several papers in the literature have employed nonlinear traction separation curves to explain the *R*-curve behaviour. The main reasoning behind these methods has been to incorporate fibre bridging into cohesive formulation. The ability of these methods to capture *R*-curve behaviour is highly dependent on the fitting parameters used for both *G*_{Icr} and *σ*_{Ic} which may not be unique, hence limiting the application of these methods. To limit the design parameter used for numerical simulation, in this paper we employ only a triangular traction separation law.

The triangular traction separation law (figure 8) is implemented through a user-defined interface element (UEL) subroutine in Abaqus [3]. UELs are placed at the interface between the two arms of the DCB, which is the intended crack path. The critical traction corresponding to the mean *G*_{Icr} was found in each case by matching the simulation's first crack initiation point with the experimental crack initiation point. This resulted in, for GC *σ*_{Ic}=13 MPa and for CC *σ*_{Ic}=47 MPa.

### (a) Double cantilever beam numerical simulation: constant *G*_{Ic}

As mentioned earlier, the first approach to analyse the variations seen in the experimental *G*_{Ic} was to use a fixed value of critical strain energy release rate along the crack path (figure 7*a*). The upper bound, mean and lower bound of the experimental *G*_{Ic} were used. Figure 9*a*,*b* shows the numerical simulation results for GC and CC specimens, respectively. Simulations of both material systems showed similar behaviour. When the minimum experimental *G*_{Ic} value is used for simulation, the crack initiation point is captured reasonably well; however, the post-initiation curve is severely under-predicted. If the mean experimental *G*_{Ic} value is used, the results show a higher load for crack initiation and the response follows the minimum load line in the propagation phase. In the case of the maximum experimental *G*_{Ic} value, the results over-predict both the crack initiation point and the propagation phase. For the GC simulation, the use of fixed *G*_{Ic} is not able to predict the *non-smooth* crack growth as seen in experiments. In the case of CC laminates, the simulation is not able to replicate the increase in load with crack progression. It is clear that the general assumption of a fixed critical strain energy release rate cannot capture the increase in load during crack propagation, and while this can be addressed by appealing to *R*-curve response, the question of *non-smooth* crack growth still begs to be resolved. Even though non-smooth crack growth has been observed here and also reported by others in the literature, a satisfactory model that can replicate non-smooth crack growth histories is presently unavailable.

### (b) Double cantilever beam numerical simulation: variable *G*_{Ic}

Here, we demonstrate that oscillatory *G*_{Ic} values along the intended crack path are necessary to capture non-smooth crack growth using a DCZM modelling approach. Simulations were conducted by varying the *G*_{Ic} value along the length of the intended crack path according to the values obtained from the experiment. For simulation of GC laminates, a saw-tooth function was used for the *G*_{Ic} distribution over every 5 mm, with an amplitude of 0.73 N mm^{−1}, and a baseline value of 1.47 N mm^{−1}, as shown in figure 7*b*. In case of CC laminates, the experimental *R*-curve (figure 7*b*) was used to vary the critical energy release rate along the crack path. Figure 10*a*,*b* shows the simulation results for GC and CC laminates, respectively. By varying the *G*_{Ic} values, both crack initiation and propagation are matched well. What is not clear is whether such an approach, when used to *predict* the fracture response of another structure made of the same material system, can replicate the experimental findings. Furthermore, the usefulness of the classical strain energy release rate as a suitable metric to drive cracks, at least for the two material systems studied here, is now brought into question.

### (c) Discrete cohesive zone method simulations: summary and discussion

From the simulation results presented, it can be inferred that the notion of a critical strain energy release rate (constant value) may not exist for a given material system even when the loading rates are quasi-static. The classical assumption that the critical strain energy release rate is a constant value is not true for the material systems studied here. More troubling is the case corresponding to *non-smooth* crack growth. These issues show that ‘critical’ parameters obtained from coupon level tests (following as closely as possible the suggested ASTM standard) may not be reliable for failure prediction in other structures. Many authors who were faced with similar issues have attributed discrepancies between experiment and analysis to rate dependence of *G*_{Ic} [8,13]. It is possible that rate dependence of *G*_{Ic} can explain the overall fracture response in the case of *smooth crack growth*. Non-smooth crack growth, spanning a large range of loading rates, including those that have been hitherto specified as ‘quasi-static’, as stated in ASTM D 5528, cannot be explained by appealing to rate dependence of *G*_{Ic}.

## 4. Revisiting classical fracture theory

From the experimental and numerical analysis presented, it is evident that the critical strain energy release rate approach might be sufficient to predict the onset of fracture but it is not sufficient for predicting the subsequent growth of a crack because *G*_{Ic} is not constant in all cases studied here. Furthermore, a model that can capture both smooth and non-smooth crack growth is required to distinguish between the variety of responses that have been observed in the results presented here and elsewhere, in the literature.

### (a) Double cantilever beam energy release rate

Consider a DCB of thickness 2*h*, with a starter crack of length *a*_{0} at the centre line and subjected to displacement control loading at the ends as shown in figure 11. The elastic energy of the system, based on simple beam theory, can be written as
4.1Differentiating (4.1) with respect to time gives
4.2Now, by definition, the strain energy release rate is *G*_{I}=−(1/*b*)(*dU*/*da*); hence, (4.2) can be rearranged and written as
4.3

A similar analysis, using higher order (including shear deformation of the DCB arms) beam theory, leads to somewhat different expressions; however, here we wish to use Euler–Bernoulli beam theory to illustrate the main points of this work. Equation (4.3) shows that the energy release rate is dependent on both crack velocity and loading rate. For quasi-static loading, it is assumed that because, either for fixed grip conditions or for dead loading, the crack is assumed to advance in infinitesimal steps, with each step leading to an adjacent state of *equilibrium*. The crack propagation history is thus assumed to be a series of equilibrium configurations where Griffith's theory [14] is assumed valid. However, it is important to note that Griffith's theory is valid only for a crack initiating from an equilibrium position. Once the crack has initiated and starts to grow, the system is no longer in equilibrium. In fact, crack advancement is due to a lack of *equilibrium*. Hence, during propagation, the application of Griffith's theory is an assumption; when this assumption is relaxed, must be taken into account in calculations. Note that these arguments remain unchanged also for cracks that propagate at large velocities. In those instances, dissipation due to kinetic energy must also be accounted for.

Let us look at the different scenarios that emerge from (4.3), for a displacement-controlled loading test . Figure 12 refers to each scenario explained as follows.

If the crack velocity is significantly larger than the loading rate, i.e. , then and .

If the crack velocity is significantly smaller than the loading rate, i.e. , the second term can no longer be ignored and the overall energy release rate will be lowered.

If the crack velocity is of the same order as the loading rate, i.e. ., then the energy release rate will be lower than in case (1), i.e. .

If but is of the same order of magnitude and is increasing with time, then the energy release rate starts at a lower value and increases at the same rate as that of the crack velocity.

If but is of the same order of magnitude and increases to reach an asymptotic value, then the energy release rate will also increase and reach an asymptote.

From the various scenarios given above, the first refers to the case when the rate of work done by the external load is negligible compared with the rate of change of elastic strain energy due to crack growth, which is the case in unstable crack growth. Case (4) provides the situation that corresponds to the results seen in the experimental *R*-curve of CC, where an asymptotic value is not reached. Case (5) corresponds to the case where the classical description of an *R*-curve is observed (with smooth crack growth), where *G*_{I} is seen to increase and reach an asymptote.

### (b) Fracture theory: summary and discussion

Equation (4.3) not only shows that *G*_{I} is dependent on loading rate and crack velocity but also provides the phenomenological reasoning behind the existence of an *R*-curve behaviour seen in mode I DCB tests. The dependence on crack velocity is shown to be inherent even at loading rates that are considered to be quasi-static. The traditional characterization of loading rate as being quasi-static or dynamic, based solely on the nature of crack growth (stable or unstable and driven by inertial effects), is not sufficient. Instead, the rate of loading with respect to crack velocity should be used as a discriminator even in cases that have hitherto been characterized as quasi-static.

## 5. Alternative hypotheses and a novel approach to fracture analysis

Griffith's [14] original motivation to formulate the problem of determining the maximum load required to advance an existing crack in a solid stemmed from an inability to use an available local solution to the stress field [15], around an elliptical cutout in a planar body of infinite extent, owing to its singular nature. This led to a formulation that invoked a macroscopic energy balance, with the result being a criterion for crack growth *initiation*. Noting the development related to the DCB in the previous section, we now pose four basic questions that one faces in a combined continuum mechanics/fracture analysis.

When (in terms of external loads) will an existing crack start to grow?

How much will the growing crack propagate by?

In which direction will the crack propagate?

Is the crack growth stable, in the sense of continued growth or arrest?

Because we are examining a symmetrical DCB specimen in this study, for which crack propagation is along the mid-plane, we will examine answers to questions (a), (b) and (d).

### (a) Crack initiation criteria

Fracture is a process involving separation and severing of atomic bonds. The rate of fracture is proportional to the number of bond ruptures per unit time. Thus, atomistically informed calculations (for instance, those that consider interactions among atoms at atomistic length scales) are necessary to properly capture precise physical processes that take place at such length scales. An approach to connect the two scales (atomistic and continuum) is to apply the concept of an equivalent continuum [16,17]. The equivalence between the continuum and the discrete atomic level, in a dynamically deforming particle system, includes conservation of internal, external and inertial work. This equivalence can be provided by introducing the concept of a virial stress, which, for atomic systems, is a tensor quantity that measures the time rate of change of momentum for spatial regions. For the discussion in this paper, we assume that there exists a critical internal stress beyond which the atomic separation is large enough to break the attractive bond. Furthermore, we assume that this internal stress is proportional to the continuum based notion of Irwin's stress intensity factor [18]. This gives us the first hypothesis that we introduce.
Hypothesis 1: crack initiation takes place when the internal stress at the crack tip

*σ*_{t} reaches or exceeds a critical value *σ*_{cr}. Alternatively, the crack initiates when the internal strain at the crack tip *ϵ*_{t} reaches a critical value *ϵ*_{cr}. The two criteria are related through an appropriate constitutive relation of the solid.

Parallels can be drawn to yield criteria used in plasticity (Von Mises or Tresca) or Irwin's stress intensity factors (critical *K*_{I/II/III}) in linear elastic fracture mechanics [19]. For simplicity and ease of analysis, we assume the system to be linear elastic, and energy loss in the system is only due to crack growth. For the condition of irreversibility to be satisfied, crack healing is not allowed.

Since the proposed hypothesis requires a measure of internal stress at the crack tip which can be obtained using atomistic or molecular dynamic simulations, and is beyond the scope of this study, we employ a more workable form of hypothesis 1, where we assume that there is a linear relationship between the internal atomistically informed stress at the crack tip and the stress state around the crack tip at a finite distance away, thereby allowing the use of Irwin's stress intensity factor as an equivalent measure of stress state around the crack tip.

We now cast a necessary condition for crack initiation: 5.1

### (b) Crack progression criterion

Many researchers have developed crack-velocity formulations especially when pertaining to dynamic fracture studies. Since crack growth releases stored energy, stress waves are emitted that transmit this energy as kinetic energy. These waves influence the further motion of the crack within the body. Hellan [20] provides a relation for particle velocities in a body based on stress waves in a one- and two-dimensional setting. Taking the simplest case of one dimension, the particle velocity due to stress wave loading is given by
5.2where *C*_{0} is a material constant and *E* is the elastic modulus. Hellan [20] also provides the relation between crack velocity and particle velocity given in (5.3). From hypothesis 1, we may assume the stress corresponding to crack formation as the critical stress *σ*_{cr}:
5.3

Combining the two equations, we get a relation for crack velocity in terms of the ratio of instantaneous stress at the crack tip, and the critical fracture stress at the crack tip. Since crack propagation cannot be reversed, the negative sign is removed. The resulting equation forms the second hypothesis.

Hypothesis 2: crack progresses with a velocity proportional to the ratio of instantaneous crack tip stress to the critical crack tip stress. 5.4

This equation is similar to an empirical power law relationship given by Evans [21] for slow crack growth in brittle materials. According to Evans, it is found that the velocity of the crack *V* is proportional to the stress intensity factor, raised to some power, *n*, i.e. *K*^{n}. This equation is also similar to the dislocation velocity relation proposed by Stein & Low [22] and Gilman & Johnston [23], where dislocation velocity *V* _{d}=(*τ*/*D*)^{m}, where *τ* is the applied resolved shear stress, and *D* and *m* are material constants. Using activation energy theory, as shown in Krausz & Krausz [24] and as discussed in Kotousov [25], it is theoretically more sound to express the crack velocity with a hyperbolic sine dependence instead of a power-law relationship. However, the power-law form is favoured for reducing test data and extending the analysis to treat fatigue crack growth.

We can see from (5.4) that the material constant *C*_{0} should have units of velocity. Hence, there exists a characteristic crack velocity for a given material configuration around the crack tip. This implies that when the stress at the crack tip just reaches the critical stress value, i.e. *σ*_{t}/*σ*_{cr}=1, the crack will progress with a velocity *C*_{0}. A similar conclusion is also provided by Kotousov [26].

Combining the two hypotheses (equations (5.1) and (5.4)), we get the necessary condition for continued crack growth, 5.5

### (c) Rate equation of energy

From the previous section, we observe that the crack velocity is an important factor in strain energy release rate calculations. A more general explanation of smooth and non-smooth crack growth can be made if we examine the rate form of the energy balance equation. Consider the DCB described in figure 11. The total potential energy (*Π*) in the system, prior to crack initiation due to a load *F*, can be written as
5.6where *U* is the elastic strain energy in the body and *W*_{ext} is work done by external load. Now, if the body develops a crack, the energy released can be attributed to dissipation associated with the crack formation. We can re-write the equation for potential energy of the system by deducting the energy dissipated *W*_{dis} to create additional area attributed to the crack:
5.7where *γ* is the fracture energy^{1} per unit area. Griffith [14], using the notion of surface energy, stated that the above equation is only valid for a system where fracture has occurred. Also, according to the theory of minimum potential energy, during the fracture process, the system will reduce its potential energy by crack growth. However, there is a subtle but key observation that needs to be made. Griffith [14] stated that the total reduction in potential energy due to fracture is equal to the increase in strain energy less the increase in surface energy. The point to note here is that the increase in strain energy has two components: first, the increase in strain energy due to a change in system stiffness and second, the increase in energy due to continued loading of the system. Let us look at the rate form of (5.7). Differentiating both sides with respect to time gives
5.8The first term on the right-hand side is the strain energy stored due to continued loading and is dependent on the loading rate. The second term corresponds to the strain energy stored owing to change in stiffness, which is dependent on instantaneous crack velocity. The third term accounts for energy dissipation due to creation of fracture surfaces, which is also a function of instantaneous crack velocity. From equation (5.8), three scenarios emerge.

If

This scenario reflects the condition where instantaneous dissipative work is larger than the cumulative energy due to external load and strain energy stored. The overall potential energy in the body will decrease, which will reflect as a drop in the force displacement plot. Indicated as curve section I in figure 13.

If

Here, the instantaneous dissipated energy equalizes the sum of energies due to external load and strain energy stored, causing the overall potential energy in the body to remain constant, which will reflect as a horizontal line in the force displacement plot. Indicated as curve section II in figure 13.

If

In this scenario, the instantaneous dissipated energy is smaller than the sum of energies due to external load and strain energy stored. The overall potential energy in the body will increase, which will reflect as an increase in the force displacement plot. Indicated as curve section III in figure 13.

Now, the continued growth of the crack in time can be visualized as a sequence of these scenarios in an order decided by the crack velocity. By combining the two hypotheses and the rate equation (5.8), we get the discrete form of instantaneous potential energy: 5.9

### (d) Crack stability and crack arrest

From equation (5.9), it is evident that non-smooth crack growth is afforded by the first scenario with very high instantaneous dissipated energy. For a given increment in time, if the stress at the crack tip equals or exceeds the critical value, then, based on the instantaneous crack velocity, one of the above three scenarios will occur. If in the next time increment, the stress condition (hypothesis 1) is still met, then again one of the three scenarios will be met. However, if the condition is not met, then the crack will arrest and the system will follow Hooke's Law.

### (e) Algorithm for computing crack growth

To summarize the information in the two previous sections, it is seen that the crack will only initiate when the stress at the crack tip (*σ*_{t}) reaches a critical value (*σ*_{cr}), as given in (5.1). After attaining the critical stress, the increment in crack growth is determined by the resistance provided by the material to crack growth, which is in turn determined by the crack propagation velocity owing to the stresses at the crack tip (5.4). We can now develop a general algorithm for fracture analysis of the DCB as shown in the flow chart in figure 14. From this flow chart, we see that there is no energy calculation required to determine crack initiation and crack progression. Since the algorithm uses incremental time steps (*N*) to calculate both the change in loading and the increment of crack growth, a sample DCB geometry was used to check for convergence and is shown in figure 15.

## 6. Implementation

### (a) Double cantilever beam formulation

Consider the case of the DCB (figure 11) with dimensions symmetrical about the neutral axis, which is also the fracture path. The relation between force and displacement is given by
6.1The stress intensity at the crack tip is given by
6.2where *b* denotes the width of the specimen and *I* the moment of inertia of arms. The critical value of the stress intensity (*K*_{cr}) can be calculated from the average critical strain energy release rate, at initiation, obtained from the experiment:
6.3Therefore, for implementation of the algorithm, the crack initiation criterion will now be *K*_{t}≥*K*_{cr}. The second parameter, the characteristic velocity *C*_{0}, was found by trial and error by matching the load–displacement curves of the experiment corresponding to each material system.

### (b) Double cantilever beam analysis

Figure 16 shows the analysis results for the two laminates, when computed using the algorithm specified in the prior section. Table 2 gives the parameters used for the analytical model.

In the case of the GC laminates, non-smooth crack growth behaviour is captured quite well. As the crack velocity is high in GC, there is a rapid loss of energy from the system after initiation. The sudden release of energy also relaxes the system (reduction in stiffness) and the stress at the crack tip reduces below the critical threshold, causing the crack to arrest (figure 17*a*). Since the body is under constant displacement loading and the crack has been arrested, the energy of the system will increase, thereby gradually increasing the stress at the crack tip. This cycle repeats to provide the non-smooth crack growth that is observed in the GC laminate.

In the CC laminate, the phenomenon is quite unique. After the crack initiates, the system does not have a net loss of energy. This is because the rate of energy input to the system due to loading is more than the rate of energy loss from the system due to crack propagation, thus giving a smooth crack growth (figure 17*b*). Because the stiffness of the system is also reducing with crack propagation, there is a point after which the input to the system due to loading drops below that lost due to crack propagation.

Because the new crack growth algorithm provides the load–deflection relation as an output, it is a simple task to *calculate* *G*_{I} as the crack progresses. Figure 18 shows the calculated *G*_{I} values for GC and CC laminates using the present analytical formulation. We now see that the release of strain energy is a *consequence* of crack growth and not the *cause* of crack growth. Thus, its use as a critical parameter to advance cracks is brought into question.

## 7. Discussion

In the cohesive zone modelling framework, implemented using the FE method, two parameters are needed: the cohesive strength and the energy release rate. It has been shown that the energy release rate may not be a constant value during crack propagation; hence, an *R*-curve is also required. Using the stress intensity factor alone, only crack initiation can be predicted. The new approach proposed here, on the other hand, required only two material constants, the critical stress at the crack tip and the characteristic velocity.

As mentioned earlier, most of current theory, with the exception of molecular dynamics or atomistic simulations, uses a macroscale approach to fracture. Current macroscopic approaches use either energy-based or strength-based techniques, both of which are very accurate at predicting the onset of cracking (crack initiation). This is because the continuum continuity definition holds until the final equilibrium point (crack initiation) is reached; hence, the energies are quantifiable and easily calculated. However, these techniques fail to provide accurate predictions for crack progression, in general. The reasons stem from the fact that during crack propagation the system is not in equilibrium. By assuming a constant strain energy release rate or stress intensity factor during crack propagation, equilibrium is enforced, which is only possible if every crack in every material propagates from one equilibrium point to another in negligible time, i.e. at extremely high crack velocities, which is not the case in actual materials.

The theory proposed in this paper overcomes this apparent disconnect by acknowledging that the process of crack growth is not an equilibrium process and introduces the crack velocity as an internal variable in the calculation of crack growth. The validity of the crack velocity as an internal variable is shown to be sound in the analytical solutions; however, the question of existence of a unique characteristic velocity *C*_{0} can be argued. The characteristic crack velocity term, *C*_{0}, is in essence a parameter that captures the effects of the process zone (dislocations, plasticity, micro-damage due to micro-cracking, etc.) and its ability to dissipate energy, which is specific to a given material system. Phenomenologically, one can argue that *C*_{0} should be a fraction of the Rayleigh wave speed and hence be a constant for a given material system, but a rigorous way to establish this—for example, using atomistic or micro-mechanical simulations—is needed.

## 8. Conclusions

We have introduced a novel approach to predict non-smooth crack growth observed in mode I fracture of a composite material DCB specimen. By postulating the existence of a critical fracture internal stress or strain, and a critical crack velocity, non-smooth crack growth (and smooth crack growth, as a special case) observed in DCB experiments was predicted. The method is able to predict both smooth and non-smooth crack growth. Using the critical fracture internal stress or strain, and a critical, material-specific crack velocity, a simple algorithm was formulated to explain and replicate the experimental results well. In our approach, the strain energy released is an *outcome* of the crack growth process instead of a value that is postulated to drive and control crack advancement. Numerical simulation, using cohesive zone models that use constant *G*_{Ic} values, was not able to capture the non-smooth crack growth seen in GC laminates nor was it able to capture the steady increase in load *after* crack initiation, observed in CC laminates. By including the rate terms in the *G*_{I} formulation, it was shown that strain energy release rate is a function of both loading rate and crack velocity. This formulation explains the reason attributed to *R*-curve behaviour and the linear decrease in *G*_{Ic} with increasing loading rate. Both phenomena are explained by the use of the crack velocity as a parameter in the energy calculation.

## Footnotes

One contribution of 15 to a Theme Issue ‘Geometry and mechanics of layered structures and materials’.

↵1 Fracture energy is more general in that it includes all mechanisms of dissipation associated with the creation of a crack, whereas Griffith's surface energy accounts only for that portion associated with the creation of free surfaces.

- This journal is © 2012 The Royal Society