The role of plastic strain gradients in the crack growth resistance of metals

TL;DR

This paper investigates how strain gradients influence crack growth resistance in metals, revealing the formation of elastic zones at crack tips that affect crack propagation and ductility.

Contribution

It introduces a model incorporating strain gradients that predicts elastic zones at crack tips, contrasting with previous analyses that neglected elastic strains.

Findings

Elastic zones form at crack tips due to strain gradients.

High stresses activate quasi-cleavage at the crack tip.

Elastic singularity reduces tensile ductility.

Abstract

Crack advance from short or long pre-cracks is predicted by the progressive failure of a cohesive zone in a strain gradient, elasto-plastic solid. The presence of strain gradients leads to the existence of an elastic zone at the tip of a stationary crack, for both the long crack and the short crack cases. This is in sharp contrast with previous asymptotic analyses of gradient solids, where elastic strains were neglected. The presence of an elastic singularity at the crack tip generates stresses which are sufficiently high to activate quasi-cleavage. For the long crack case, crack growth resistance curves are predicted for a wide range of ratios of cohesive zone strength to yield strength. Remarkably, this feature of an elastic singularity is preserved for short cracks, leading to a severe reduction in tensile ductility. In qualitative terms, these predictions resemble those of discrete dislocation calculations, including the concept of a dislocation-free zone at the crack tip.

Equations

$$\int_{V}\Big{(}\sigma_{ij}\delta\varepsilon_{ij}^{e}+q_{ij}\delta\varepsilon_{ij}^{p}+\tau_{ijk}\delta\varepsilon_{ij,k}^{p}\Big{)}dV=\int_{S}\left(T_{i}\delta u_{i}+t_{ij}\delta\varepsilon_{ij}^{p}\right)dS$$

$${\sigma }_{ij,j}=0$$

$${\tau }_{ijk,k}+{\sigma }_{ij}^{\prime }-q_{ij}=0$$

$$T_i={\sigma }_{ij}{n}_j$$

$$t_{ij}={\tau }_{ijk}{n}_k$$

$$\Psi \left({\epsilon }_{ij}^e,{\epsilon }_{ij,k}^p\right)=\frac{1}{2}{\epsilon }_{ij}^e{C}_{ijkl}{\epsilon }_{kl}^e+\frac{1}{2}\mu L_E^2{\epsilon }_{ij,k}^p{\epsilon }_{ij,k}^p$$

$$\delta \Psi ={\sigma }_{ij}\delta {\epsilon }_{ij}^e+{\tau }_{ijk}^E\delta {\epsilon }_{ij,k}^p$$

$${\sigma }_{ij}=\frac{\partial \Psi }{\partial {\epsilon }_{ij}^e}=C_{ijkl}\left({\epsilon }_{kl}-{\epsilon }_{kl}^p\right)$$

$${\tau }_{ijk}^E=\frac{\partial \Psi }{\partial {\epsilon }_{ij,k}^p}=\mu L_E^2{\epsilon }_{ij,k}^p$$

$${\dot{W}}^p=\Sigma {\dot{E}}^p$$

$${\dot{E}}^p={\left(\frac{2}{3}{\dot{\epsilon }}_{ij}^p{\dot{\epsilon }}_{ij}^p+L_D^2{\dot{\epsilon }}_{ij,k}^p{\dot{\epsilon }}_{ij,k}^p\right)}^{1/2}$$

$$\delta {\dot{W}}^p=\Sigma \delta {\dot{E}}^p=q_{ij}^D\delta {\dot{\epsilon }}_{ij}^p+{\tau }_{ijk}^D\delta {\dot{\epsilon }}_{ij,k}^p$$

$$q_{ij}^D=\frac{2}{3}\frac{\Sigma }{{\dot{E}}^p}{\dot{\epsilon }}_{ij}^p\text{and}{\tau }_{ijk}^D=\frac{\Sigma }{{\dot{E}}^p}{L}_D^2{\dot{\epsilon }}_{ij,k}^p$$

$$\Sigma ={\left(\frac{3}{2}{q}_{ij}^D{q}_{ij}^D+\frac{1}{L_D^2}{\tau }_{ijk}^D{\tau }_{ijk}^D\right)}^{1/2}$$

$$\Sigma ={\sigma }_F\left(E^p\right)V({\dot{E}}^p)$$

$${\sigma }_F={\sigma }_Y{\left(1+\frac{E{E}^p}{{\sigma }_Y}\right)}^N$$

$$V({\dot{E}}^p)=\{\begin{array}{ll}{\dot{E}}^p/\left(2{\dot{\epsilon }}_0\right)& \text{if }{\dot{E}}^p/{\dot{\epsilon }}_0\le 1\\ 1-{\dot{\epsilon }}_0/\left(2{\dot{E}}^p\right)& \text{if }{\dot{E}}^p/{\dot{\epsilon }}_0>1\end{array}$$

$$u_i=\sum _{n=1}^{D_u}{N}_i^n{U}^n$$

$${\epsilon }_{ij}^p=\sum _{n=1}^{D_{{\epsilon }_p}}{M}_{ij}^n{\epsilon }_p^n$$

$${\left[\begin{array}{c}\mathbf{u}\\ {\mathbf{\epsilon }}_{\mathbf{p}}\end{array}\right]}_{t+\Delta t}={\left[\begin{array}{c}\mathbf{u}\\ {\mathbf{\epsilon }}_{\mathbf{p}}\end{array}\right]}_t-{\left[\begin{array}{cc}{\mathbf{K}}_{u,u}& {\mathbf{K}}_{u,{\epsilon }^p}\\ {\mathbf{K}}_{{\epsilon }^p,u}& {\mathbf{K}}_{{\epsilon }^p,{\epsilon }^p}\end{array}\right]}_t^{-1}{\left[\begin{array}{c}{\mathbf{R}}_u\\ {\mathbf{R}}_{{\epsilon }^p}\end{array}\right]}_t$$

$${\mathbf{R}}_u^n={\int }_V{\sigma }_{ij}{B}_{ij}^{n}dV-{\int }_S{T}_i{N}_i^{n}dS$$

$${\mathbf{R}}_{{\epsilon }^p}^n={\int }_V\left[(q_{ij}-{\sigma }_{ij}^{\prime })M_{ij}^n+{\tau }_{ijk}{M}_{ij,k}^n\right]dV-{\int }_S{t}_{ij}{M}_{ij,k}^{n}dS$$

$$u_i=\frac{K}{E}{r}^{1/2}{f}_i\left(\theta ,\nu \right)$$

$$f_x=\frac{1+\nu }{\sqrt{2\pi }}\left(3-4\nu -\cos \theta \right)\cos \left(\frac{\theta }{2}\right)$$

$$f_y=\frac{1+\nu }{\sqrt{2\pi }}\left(3-4\nu -\cos \theta \right)\sin \left(\frac{\theta }{2}\right)$$

$$R_p=\frac{1}{3\pi }{\left(\frac{K}{{\sigma }_Y}\right)}^2$$

$${\epsilon }_{yy}^p=A+B{r}^{3/2}$$

$${\Gamma }_0=\frac{1}{2}\hat{\sigma }\left({\delta }_c+{\delta }_2-{\delta }_1\right)$$

$$R_0=\frac{1}{3\pi (1-{\nu }^2)}\frac{E{\Gamma }_0}{{\sigma }_Y^2}$$

$$K=\frac{1.15E{u}^{\infty }\sqrt{\pi a}}{(1-{\nu }^2)H}$$

$$a>2.5\frac{K^2}{{\sigma }_Y^2}$$

$$\frac{u^{\infty }}{H}<\frac{{\sigma }_Y(1-{\nu }^2)}{1.8\sqrt{\pi }E}$$

Full text

††journal: Journal of the Mechanics and Physics of Solids

1 Introduction

Conventional plasticity theories, such as $J_{2}$ flow theory, predict that the tensile stress state ahead of a mode I crack in an elastic-perfectly plastic solid is on the order of three times the yield stress $\sigma_{Y}$. The factor of 3 arises from plastic constraint effects, and is explained in terms of the Prandtl stress field for a flat punch, see for example Rice (1968). Ductile fracture by void growth is promoted by this high hydrostatic stress (Rice and Tracey, 1969; Hancock and Mackenzie, 1976; McClintock, 1968). This level of tensile stress is far below the cleavage strength (typically on the order of $10\sigma_{Y}$) yet cleavage fracture in the presence of significant plastic flow has been observed, see for example Elssner et al. (1994), Bagchi and Evans (1996), Korn et al. (2002). Additional physics is needed to explain the occurrence of cleavage failure in the presence of plasticity. As argued by Wei and Hutchinson (1997) and Jiang et al. (2010), cleavage can occur if the stress elevation due to plastic strain gradients is sufficient to attain the ideal strength. This mechanism is reassessed in the present study.

The micromechanical basis for strain gradient effects is the elevation of dislocation-based hardening, and thereby of stress levels, as first appreciated by Nye (1953) and Cottrell (1964). Additional, dislocation-based arguments were advanced by Ashby (1970) and Brown and Stobbs (1976). In broad terms, plastic strain gradients demand the existence of geometrically necessary dislocations (GNDs), and this elevation in dislocation density increases the flow strength by mechanisms such as forest hardening (Fleck et al., 1994). Such strain gradient strengthening can explain a wide range of phenomena such as the Hall-Petch size effect, see for example Shu and Fleck (1999). The Hall-Petch effect is based on grain-to-grain plastic anisotropy such that strain gradients are present at the grain-to-grain level when the macroscopic strain field is uniform (or non-uniform). Additional strain gradient effects arise when the macroscopic strain field is non-uniform, as near a crack tip, or in simple test geometries such as a wire in torsion (Fleck et al., 1994), a beam in bending (Stölken and Evans, 1998) or at the tip of an indenter (Stelmashenko et al., 1993; Poole et al., 1996; Nix and Gao, 1998). A large literature has emerged on strain gradient plasticity (SGP) formulations (Aifantis, 1984; Fleck and Hutchinson, 1993, 2001; Gao et al., 1999). The pivotal step in constructing these phenomenological theories is to express the plastic work in terms of both plastic strain and plastic strain gradient, thereby introducing a length scale into the material description. Recent SGP models incorporate both dissipative (that is, unrecoverable) and energetic (that is, recoverable) gradient contributions (Gudmundson, 2004; Gurtin and Anand, 2005; Fleck and Willis, 2009a, b).

Recently, the effect of strain gradients in elevating crack tip stress levels has been emphasized in a number of numerical investigations on stationary cracks (Jiang et al., 2001; Komaragiri et al., 2008; Martínez-Pañeda and Betegón, 2015; Martínez-Pañeda and Niordson, 2016). It has been suggested that this elevation in stress influences fatigue damage (Sevillano, 2001; Brinckmann and Siegmund, 2008; Pribe et al., 2019), microvoid cracking (Tvergaard and Niordson, 2008), and hydrogen embrittlement (Martínez-Pañeda et al., 2016b, a). Stress elevation due to plastic strain gradients is also relevant to propagating cracks. Wei and Hutchinson (1997), and Wei et al. (2004), quantified the dependence of steady state fracture toughness $K_{SS}$ upon material length scale $\ell$ for the Fleck and Hutchinson (1997) gradient theory and Gao et al. (1999) gradient theory, respectively. Recently, Seiler et al. (2016) computed the initial stages of the crack growth resistance curve for a viscoplastic, strain gradient plasticity theory (Huang et al., 2004), and investigated the dependence of R-curves on viscoplastic constitutive parameters and on the intrinsic material length scale.555We note in passing that the Huang et al. (2004) theory is a lower order theory that neglects higher order stresses. The present study assumes the presence of higher order stresses that are work conjugate to plastic strain gradients. The recent strain gradient theory of Gudmundson (2004) (see also Gurtin and Anand, 2005; Fleck and Willis, 2009a) has additional features that can significantly influence crack growth resistance: this motivates the present paper. First, the recent asymptotic analysis of Martínez-Pañeda and Fleck (2019) for a stationary crack in a dissipative strain gradient solid reveals the existence of an elastic crack tip zone, reminiscent of a dislocation-free zone (Suo et al., 1993). Second, both energetic and dissipative length scales enter the constitutive relations; their influence on fracture problems has not yet been assessed. Both features are explored here in the context of both stationary and propagating cracks. In addition, we explore the effect of crack length in relation to the material length scales and to the fracture length scale of the crack tip process zone. Thereby, insight is gained into the role of strain gradients on the behaviour of short cracks.

The present study is structured as follows. The constitutive model is presented in Section 2, including the phenomenological formulation of strain gradient plasticity, and the implicit finite element implementation. The asymptotic response at the tip of a stationary crack in a strain gradient solid is investigated in Section 3. Crack growth is explored in Section 4 in two steps. First, R-curves are computed for a long crack by means of a cohesive zone, and the relative role of energetic versus dissipative strain gradient terms is quantified. Second, the short crack case is examined and we compute the sensitivity of the macroscopic stress versus strain response to crack length $a$ and to the material length scale $\ell$. The sensitivity of the tensile ductility to the presence of a short crack is emphasized: it is shown that strain gradients play a major role. Finally, concluding remarks are given in Section 5.

2 Strain Gradient Plasticity

2.1 Flow theory

2.1.1 Variational principles and balance equations

We adopt a small strain formulation. The total strain rate $\dot{\varepsilon}_{ij}$ is the symmetric part of the spatial gradient of the displacement rate $\dot{u}_{i}$, such that $\dot{\varepsilon}_{ij}=\left(\dot{u}_{i,j}+\dot{u}_{j,i}\right)/2$; $\dot{\varepsilon}_{ij}$ decomposes additively into an elastic part, $\dot{\varepsilon}^{e}_{ij}$, and a plastic part, $\dot{\varepsilon}^{p}_{ij}$. Write $\sigma_{ij}$ as the Cauchy stress, $q_{ij}$ as the so-called micro-stress tensor (work-conjugate to the plastic strain $\varepsilon_{ij}^{p}$) and $\tau_{ijk}$ as the higher order stress tensor (work-conjugate to the plastic strain gradient $\varepsilon_{ij,k}^{p}$). For a solid of volume $V$ and surface $S$, the principle of virtual work reads,

[TABLE]

The right-hand side of Eq. (1) includes both conventional tractions $T_{i}$ and higher order tractions $t_{ij}$. Write $\sigma_{ij}^{\prime}$ as the deviatoric part of $\sigma_{ij}$, and write $n_{k}$ as the unit outward normal to the surface $S$. Then, upon making use of the Gauss divergence theorem, equilibrium within $V$ reads

[TABLE]

and on $S$ reads,

[TABLE]

2.1.2 Constitutive description

Gudmundson (2004) and Fleck and Willis (2009a) explain that both $q_{ij}$ and $\tau_{ijk}$ can have dissipative and energetic contributions: $q_{ij}=q_{ij}^{D}+q_{ij}^{E}$ and $\tau_{ijk}=\tau_{ijk}^{D}+\tau_{ijk}^{E}$, where the superscripts $D$ and $E$ denote dissipative and energetic, respectively. In general, the Cauchy stress $\sigma_{ij}$, along with $q_{ij}^{E}$ and $\tau_{ijk}^{E}$, are derived from the bulk free energy of the solid $\Psi$. In the present study, we shall assume that $q_{ij}^{E}$ vanishes and thus limit attention to a solid that displays isotropic hardening in the absence of a strain gradient. The significance of a finite value of $q_{ij}^{E}$ (with $\tau^{D}_{ijk}\equiv 0$) has been explored in the recent study of Martínez-Pañeda and Fleck (2018); here, we limit our focus to the role of kinematic hardening associated with the gradient of plastic strain. Accordingly, the bulk free energy $\Psi$ of the solid depends upon the elastic strain $\varepsilon_{ij}^{e}$ and the plastic strain gradient $\varepsilon_{ij,k}^{p}$ but not upon the plastic strain $\varepsilon_{ij}^{p}$, such that

[TABLE]

Here, $C_{ijkl}$ is the isotropic elastic stiffness tensor, $\mu$ is the shear modulus and $L_{E}$ is the energetic constitutive length parameter. Upon noting that

[TABLE]

the energetic stress quantities follow as

[TABLE]

Now consider plastic dissipation. For both the rate dependent case, and the rate independent limit, we define the plastic work rate as

[TABLE]

where $\Sigma$ is an effective stress, work-conjugate to a gradient-enhanced effective plastic strain rate $\dot{E}^{p}$. The latter is defined phenomenologically as

[TABLE]

where $L_{D}$ is the dissipative length scale. Upon noting that

[TABLE]

the constitutive relations for the dissipative stress quantities read

[TABLE]

The effective stress $\Sigma$ is readily obtained by substitution of (11) into (9) to give

[TABLE]

2.2 Numerical implementation

A robust and efficient finite element framework is now presented in order to model crack propagation in a rate independent gradient plasticity solid. An implicit time integration scheme is developed for both energetic and dissipative higher order contributions.

Gradient plasticity theories are commonly implemented within a rate-dependent setting, thereby taking advantage of computational advantages and circumventing complications in the corresponding time independent model associated with identifying active plastic zones (see, for example, Nielsen and Niordson, 2014). The mathematical foundations and associated variational structure for both the rate dependent and rate independent cases are given by Fleck and Willis (2009a, b). Here, we make use of the viscoplastic law by Panteghini and Bardella (2016), and exploit the fact that it adequately approximates the rate-independent solution in a computationally efficient manner. The effective flow resistance is related to the gradient-enhanced effective plastic flow rate through a viscoplastic function,

[TABLE]

where the current flow stress $\sigma_{F}$ depends upon an initial yield stress $\sigma_{Y}$ and on $E^{p}$ via a hardening law. Here, we adopt the following isotropic hardening law,

[TABLE]

in terms of the Young’s modulus $E$ and strain hardening exponent $N$ $(0\leq N\leq 1)$. Following Panteghini and Bardella (2016) the viscoplatic function $V(\dot{E}^{p})$ is defined as666This choice has the advantage that the consistent stiffness matrix, as defined in the Supplementary Material, remains finite as $\dot{E}^{p}\to 0$.

[TABLE]

in terms of a reference strain rate $\dot{\varepsilon}_{0}$. A sensitivity study for a sufficiently small choice of $\dot{\varepsilon}_{0}$ is conducted to ensure that the rate independent limit is attained in all the results presented subsequently. The reader is referred to Panteghini and Bardella (2016) for a more detailed interpretation of $\dot{\varepsilon}_{0}$.

The finite element scheme takes displacements and plastic strains as the primary kinematic variables, in accordance with the theoretical framework. $C_{0}$-continuous finite elements are adopted since the differential equations are of second order. The displacement field $u_{i}$ at position $\bm{x}$ is written in terms of the shape functions $N_{i}^{n}$ and associated nodal displacements $U^{n}$, where $n$ denotes the degree of freedom, such that

[TABLE]

Here, $D_{u}$ is the total number of degrees of freedom for the nodal displacements. Likewise, the plastic strain field $\varepsilon_{ij}^{p}$ is expressed in terms of the shape functions $M_{ij}^{n}$ and associated nodal quantities $\varepsilon_{p}^{n}$ as

[TABLE]

where $D_{\varepsilon_{p}}$ denotes the total number of degrees of freedom for the nodal plastic strain components. Quadratic shape functions are employed for interpolation of both displacements and plastic strains. Accordingly, the plastic strain gradient $\varepsilon_{ij,k}^{p}$ and the total strain $\varepsilon_{ij}$ are related to the nodal plastic strains and displacements through $M_{ij,k}^{n}$ and the strain-displacement matrix $B_{ij}^{n}$, respectively; see the Supplementary Material for further details.

The non-linear system of equations is solved iteratively by the Newton-Raphson method from time step $t$ to $\left(t+\Delta t\right)$

[TABLE]

where the residuals comprise the out-of-balance forces,

[TABLE]

[TABLE]

and the components of the consistent stiffness matrix $\bm{K}$ are obtained by differentiating the residuals with respect to the incremental nodal variables. Details are given in the Supplementary Material.

The numerical scheme is implemented in the commercial finite element package ABAQUS by means of a user element subroutine. To the best of the authors’ knowledge, it constitutes the first Backward Euler implementation of the Gudmundson (2004) class of strain gradient plasticity theories, including energetic and dissipative higher order contributions.777The code is made freely available at www.empaneda.com, hoping to facilitate research and enabling readers to reproduce the results. The reader is referred to Danas et al. (2012) and Dahlberg and Faleskog (2013) for implicit implementations for the case of dissipative higher order stresses (with $\tau_{ijk}^{E}=0$).

3 Stationary crack analysis

We assume that small scale yielding conditions prevail and we make use of a boundary layer formulation to prescribe a remote $K$ field. Consider a crack with its tip at the origin and with the crack plane along the negative axis of the Cartesian reference frame $(x,y)$. The elastic response of the solid is characterised by the Young’s modulus $E$ and Poisson’s ratio $\nu$. Then, an outer $K$ field is imposed by prescribing nodal displacements on the outer periphery of the mesh as

[TABLE]

where the subscript index $i$ equals $x$ or $y$, and the functions $f_{i}\left(\theta,\nu\right)$ are given by

[TABLE]

and

[TABLE]

in terms of polar coordinates $(r,\theta)$ centred at the crack tip. A representative value for the plastic zone size $R_{p}$ is given by the Irwin expression

[TABLE]

for a stationary crack in an elastic, ideally plastic solid. Upon exploiting reflective symmetry about the crack plane, only half of the finite element model is analysed. A mesh sensitivity study reveals that the domain is adequately discretised by means of 5200 plane strain, quadratic, quadrilateral elements. The characteristic element size is $R_{p}/7500$ and the outer radius of the boundary layer is $5000R_{p}$, ensuring small scale yielding conditions.

A representative small scale yielding solution is presented in Fig. 1 for the choice $N=0.1$, $\sigma_{Y}/E=0.003$, and $\nu=0.3$. Insight is gained into the relative role of $L_{E}$ and $L_{D}$ by considering the three cases $L_{E}=L_{D}=0.05R_{p}$, $L_{E}=0.05R_{p}$ ($L_{D}=0$), and $L_{D}=0.05R_{p}$ ($L_{E}=0$). An elastic zone exists directly ahead of the crack tip if $L_{E}$ and/or $L_{D}$ is finite. The plastic strain $\varepsilon_{yy}^{p}$ reaches a plateau value over $0<r<\ell$, see Fig. 1b. Consequently, the stress state within this crack tip zone is elastic in nature. This finding is supported by a plot of tensile stress $\sigma_{yy}$ as a function of $r$ directly ahead of the crack tip ($y=0$), see Fig. 1a. The stress component $\sigma_{yy}$ scales as $r^{-1/2}$ for sufficiently small $r$. Likewise, the elastic strain component $\varepsilon_{yy}^{e}$ scales as $r^{-1/2}$ for $r/\ell<1$; from Hooke’s law and Fig. 1 it is clear that the elastic strain dominates the plastic strain $\varepsilon_{yy}^{e}>>\varepsilon_{yy}^{p}$. Beyond the plastic zone ($r/R_{p}>1$) the stress state again converges to the elastic $K$-field and $\sigma_{yy}$ scales as $r^{-1/2}$. Thus, both an outer and an inner $K$ field exist. The plastic strain distribution $\varepsilon^{p}_{yy}(r)$ is relatively insensitive to the choice of values of $L_{E}$ and $L_{D}$ in Fig. 1b. In all three cases, the plastic strain is almost constant over $0<r<\ell$. In their recent asymptotic analysis, Martínez-Pañeda and Fleck (2019) find that the leading order terms of the plastic strain $\varepsilon_{yy}^{p}$ along $\theta=0^{\circ}$ are

[TABLE]

for the case $L_{D}\neq 0$, $L_{E}=0$, where $(A,B)$ are functions of $R_{p}$. In the present finite element study, it is also found that the plastic strain is finite at the crack tip when energetic higher order terms are present. We note in passing that the plastic strain is not sufficiently singular to contribute to the $J$-integral as the crack tip is approached. Instead, the $J$-integral is determined solely by the elastic strain state near the crack tip.

The relative insensitivity of the stationary crack response in Fig. 1 to the ratio $L_{E}/L_{D}$ leads us to focus on a single reference length scale $L_{E}=L_{D}=\ell$. The tensile stress $\sigma_{yy}$ directly ahead of the crack tip is shown in Fig. 2a for selected values of $\ell/R_{p}$, with $\ell/R_{p}=0$ corresponding to the conventional plasticity limit. In all cases, except for $\ell/R_{p}=0$, the asymptotic stress state is elastic in nature, with the tensile stress exhibiting an $r^{-1/2}$ singularity. Now place a cohesive zone at the crack tip; then, a cohesive zone strength $\hat{\sigma}$ on the order of $4\sigma_{Y}$ is sufficient to prevent crack advance in the conventional solid but not in the strain gradient case. In broad terms, the presence of strain gradients elevate stress and diminish the degree of plastic straining near the crack tip. To illustrate this, the crack opening profile for the strain gradient solid is compared to that of the conventional elasto-plastic solid ($\ell/R_{p}=0$) and to that of an elastic solid in Fig. 2b. The opening profile in the strain gradient plasticity solid ($\ell>0$) is close to the elastic case as $r\to 0$, and is close to the conventional elasto-plastic solid as $r/R_{p}\to 1$.

The above results for the asymptotic crack tip fields are in marked contrast to those obtained by Chen et al. (1999). They considered the asymptotic crack tip singular field for a mode I crack in a rigid power-law hardening strain gradient solid, as introduced by Fleck and Hutchinson (1997). Chen et al. (1999) neglected elastic strains by assuming, a priori, that the crack tip plastic strain field dominates the elastic strains. They find that the crack tip plastic strain field scales as $r^{N/(N+1)}$ in order for the strain energy density to scale as $r^{-1}$ as the crack tip is approached (thereby giving a finite value of the $J$-integral at the crack tip). Consequently, the plastic strain vanishes as $r\to 0$. The asymptotic analysis of Chen et al. (1999) carries over directly to our case if we assume that elastic strains are negligible in comparison with plastic strains at the crack tip. But in so doing, we find that the plastic strain vanishes at the crack tip and consequently the elastic strain vanishes at the crack tip also. This result is unphysical: the crack tip is sharp and will give rise to a strain concentration. We conclude that the elastic strains must dominate the plastic strains as the crack tip is approached. Consequently, an elastic $K$-field exists at the crack tip, such that the elastic strains and Cauchy stresses scale as $r^{-1/2}$. The finite element results fully support this finding, and reveal that the crack tip plastic strain is finite.

Chen et al. (1999) argued that the asymptotic field is not a physical representation over a small region (a small fraction of $\ell$) from the crack tip on the basis that the traction is negative in that zone. We draw an alternative conclusion: within a zone of order $\ell$, the crack tip field is elastic in nature. The asymptotic field of Chen et al. (1999) has no zone of validity as it neglects elastic straining.

4 Analysis of a growing crack

In the current study, we investigate crack growth from either a short or a long crack by making use of strain gradient plasticity theory. In the long crack case, R-curves are obtained and the present study thereby extends the results of Tvergaard and Hutchinson (1992) by incorporating the role of plastic strain gradients. Failure by cleavage, by void growth or by other mechanisms is idealised by an assumed traction $T$ versus separation $\delta$ law along a cohesive strip directly ahead of the crack tip, see Fig. 3a. Following Tvergaard and Hutchinson (1992), a trapezoidal shape is assumed for the $T(\delta)$ relation, as characterised by three salient values of opening ($\delta_{1}$, $\delta_{2}$, $\delta_{c}$) and a strength $\hat{\sigma}$ (see Fig. 3b). We hold fixed the ratios $\delta_{1}/\delta_{c}=0.15$ and $\delta_{2}/\delta_{c}=0.5$, and thereby treat $(\delta_{c},\hat{\sigma})$ as the two primary parameters that define the cohesive zone law. The work of fracture $\Gamma_{0}$ is the area under the $T(\delta)$ curve, as given by

[TABLE]

It follows directly from the surface work terms on the right hand side of (1) that, in general, a cohesive zone can support both tractions $T_{i}$ and higher order tractions $t_{ij}$. We assume that the tensile traction $T$ on the cohesive zone depends only upon the crack opening displacement $\delta$. Further, we assume that the higher order traction $t_{ij}$ vanishes on the surface of the cohesive zone; this is a natural boundary condition in the finite element formulation. The use of a cohesive zone model embedded within an elasto-plastic solid gives insight into both ductile fracture and cleavage by suitable choices of the cohesive zone parameters $\hat{\sigma}$ and $\Gamma_{0}$.

We proceed to evaluate the influence of crack length, material length scale $\ell$ of the strain gradient solid, and a representative fracture process zone size888$R_{0}$ corresponds to the plastic zone size in a conventional solid at the onset of crack growth.

[TABLE]

on the fracture response. Regimes of behaviour are sketched in non-dimensional space $(\ell/R_{0},R_{0}/a)$ in Fig. 4a. Our analysis spans the regimes of small scale yielding (for which an outer $K$-field exists), $J$-controlled fracture and large scale plasticity. A representative crack tip plastic zone, computed at crack initiation, is shown in Fig. 4b for the case of small scale yielding. The plastic zone size is defined by the contour along which the von Mises effective stress equals the initial yield stress. Crack growth resistance is assessed for three distinct regimes in $(\ell/R_{0},R_{0}/a)$ space, as shown by the ellipses in Fig. 4a. Section 4.2 deals with the fracture response of a strain gradient plasticity solid with a long crack while the mechanics of short flaws and the influence of crack length on the fracture response are addressed in Section 4.3.

4.1 Boundary value problem

We investigate crack initiation and subsequent growth in an edge-cracked plate loaded in uniaxial tension under plane strain conditions, see Fig. 5. The same geometry is used for the study of long and short pre-cracks. The specimen has a height-to-width ratio of $H/W=4$ and an initial crack length of $a/W=0.1$. The cohesive zone model outlined above is employed to model crack initiation and growth. Following Wei and Hutchinson (1997), micro-free boundary conditions $t_{ij}=0$ are adopted on the symmetry plane. Cohesive elements with 6 nodes and 12 integration points are implemented by means of a user element (UEL) subroutine, as described elsewhere (del Busto et al., 2017). The finite element mesh is refined ahead of the initial crack tip to ensure that the element size is able to resolve the fracture process zone. Specifically, the model consists of approximately $10^{6}$ degrees of freedom and the characteristic element length equals $R_{0}/100$. Post-processing of the results is performed with Abaqus2Matlab (Papazafeiropoulos et al., 2017).

It is widely appreciated that elastic snap-back instabilities can arise when cohesive elements experience stiffness degradation, complicating the modelling of the post-instability behaviour. The simultaneous reduction of the remote load and the prescribed displacement inevitably triggers convergence problems in quasi-static finite element computations. A numerical strategy to overcome these instabilities lies in prescribing a quantity that increases monotonically throughout the loading history while making the remote load an output of the model (Tvergaard, 1976; Segurado and LLorca, 2004). In the present study, a control algorithm is used to prescribe the crack tip opening and obtain the displacement at the remote boundary by ensuring global force equilibrium. Details are given in the Supplementary Material.

4.2 Small scale yielding response

The $K$ calibration for the specimen geometry of Fig. 5 was determined as follows. Consider the elastic solid, absent a cohesive zone and apply a uniform remote displacement $u_{y}=u^{\infty}$ on the top edge, with $T_{x}\equiv 0$. A linear elastic finite element calculation reveals that the average traction $\bar{T}$ on the top edge is $\bar{T}\approx E^{\prime}u^{\infty}/H$ where $E^{\prime}=E/(1-\nu^{2})$, assuming plane strain conditions. A contour integral evaluation of the stress intensity factor $K$ at the crack tip gives $K=1.15\bar{T}\sqrt{\pi a}$. Thus, for the small scale yielding case of limited crack tip plasticity, the remote $K$ value for the geometry of Fig. 5 is given by,

[TABLE]

Small scale yielding prevails when,

[TABLE]

in accordance with ASTM E1820. This places an upper limit on the value of $u^{\infty}/H$ for small scale yielding; rearrangement of (28) and (29) implies,

[TABLE]

This condition was satisfied in the following determination of the R-curve under small scale yielding conditions. Consider a long crack subjected to a remote load $K$. Crack initiation occurs within the cohesive zone at a value of $K$ equal to

[TABLE]

Dimensional analysis implies that the crack growth resistance for a long crack depends on the following dimensionless groups

[TABLE]

where $(N,\,\sigma_{Y}/E,\,\nu)$ are held fixed in the present study, along with the values of $\delta_{1}/\delta_{c}$ and $\delta_{2}/\delta_{c}$ in (26). The computed crack growth resistance curves for $\hat{\sigma}/\sigma_{Y}=3.8$ and for selected values of the constitutive length scales $L_{E}=L_{D}=\ell$, relative to $R_{0}$, are shown in Fig. 6a. The influence of plastic strain gradients in lowering the fracture resistance is evident: the steepness of the R-curve and the steady state value $K_{SS}/K_{0}$ diminish with increasing $\ell/R_{0}$. Seiler et al. (2016) considered the initial stage of the R-curve for a visco-plastic solid whereby the viscoplastic strain rate $\dot{\varepsilon}^{VP}$ scales as $\sigma^{m}$ where $1<m<\infty$. They showed that the sensitivity of the R-curve to the material length scale $\ell$ increases with increasing $m$. In the present study we consider the rate independent limit, $m\to\infty$, and a high sensitivity of the R-curve to length $\ell$ is, indeed, observed.

The crack tip opening angle has been used as a criterion for crack growth resistance in metallic alloys (Kanninen and Popelar, 1985). The dependence of the crack tip opening angle upon crack extension is shown in Fig. 6b. Here, the crack opening angle $\alpha$, as defined in the inset of Fig. 6b, is almost independent of $\Delta a$ after an initial transient phase. The steady state value of $\alpha$ decreases with increasing $\ell/R_{0}$, consistent with the crack opening profile for a stationary crack, as shown in Fig. 2b. It is clear that the crack tip opening angle is sensitive to strain gradient effects. In turn, this is due to the sensitivity of the plastic strain field to strain gradient effects. This is now explored in detail.

The plastic field surrounding the tip of a crack propagating at steady state is examined in Fig. 7. A von Mises measure of plastic strain is defined as,

[TABLE]

and its contours are plotted in Fig. 7 for strain gradient plasticity, with $\ell/R_{0}=0.05$, and also for the conventional plasticity case $\ell=0$. For the choice $\ell/R_{0}=0.05$, plastic strains attain a plateau value of $\varepsilon_{p}/\varepsilon_{Y}=3$ at a distance on the order of $\ell$ from the crack tip. Furthermore, the maximum level of plastic strain is not attained at the crack tip, a feature which also observed in discrete dislocation plasticity (Chakravarthy and Curtin, 2010). This contrasts with the conventional plasticity case, see Fig. 7b. In addition, plastic strains are approximately one order of magnitude larger than for the strain gradient plasticity case.

The dependence of $K_{SS}/K_{0}$ upon $\hat{\sigma}/\sigma_{Y}$ is given in Fig. 8 for selected values of $\ell/R_{0}$. There is a qualitative change when $\ell/R_{0}$ is increased from zero to a finite value. For $\ell/R_{0}=0$, continued crack advance (at $K=K_{SS}$) is precluded for $\hat{\sigma}/\sigma_{Y}>4$; the level of crack tip stress is unable to overcome the cohesive strength when $\hat{\sigma}/\sigma_{Y}\geq 4$. In contrast, when strain gradients are taken into account, the crack tip stresses can attain any value of cohesive strength, and $K_{SS}/K_{0}$ increases monotonically with increasing $\hat{\sigma}/\sigma_{Y}$. However, the degree of elevation of the R-curve, $K_{SS}/K_{0}$, decreases with increasing $\ell/R_{0}$ for any given $\hat{\sigma}/\sigma_{Y}$; this is consistent with the results shown in Fig. 6a for the choice $\hat{\sigma}/\sigma_{Y}=3.8$. Recall that the choice of $\hat{\sigma}/\sigma_{Y}\approx 10$ is representative of the mechanism of quasi-cleavage in metallic alloys: the crack tip advances by cleavage, but surrounded by a plastic zone. The predictions of Fig. 8 show that a shallow R-curve can exist for such a case: $K_{SS}/K_{0}$ equals 4 for $\hat{\sigma}/\sigma_{Y}=10$ and $\ell/R_{0}=0.06$. The qualitative response is similar to that obtained by Wei and Hutchinson (1997) for the case of Fleck and Hutchinson (1997) strain gradient theory. However, significant quantitative differences arise. If we consider a cohesive strength of $\hat{\sigma}/\sigma_{Y}\approx 10$ in both studies, then a value of $K_{SS}/K_{0}$ on the order of 4 is achieved for $\ell/R_{0}$ an order of magnitude smaller than that found by Wei and Hutchinson (1997).

Finally, we investigate the relative influence of energetic and dissipative gradient contributions to the R-curve. Crack growth resistance curves are shown in Fig. 9a for three cases: (i) $L_{E}=10L_{D}=\ell$, (ii) $L_{D}=10L_{E}=\ell$, and (iii) $L_{D}=L_{E}=\ell$ (i.e., the reference case). All of the R-curves are for $\hat{\sigma}/\sigma_{Y}=5$, and results are given for the two choices $\ell/R_{0}=0.03$ or $\ell/R_{0}=0.05$. The R-curve is steepest for $\ell/R_{0}=0.03$ and $L_{D}=10L_{E}=\ell$, for which dissipative hardening dominates. Combined energetic and dissipative hardening with $L_{E}=L_{D}=\ell$ emphasizes the role of strain gradients and leads to a less steep R-curve; the choice $L_{E}=10L_{D}=\ell$ (energetic hardening dominant) is the intermediate case. Consistent with the results shown in Fig. 8, for which $L_{E}=L_{D}=\ell$, the R-curve is less steep and $K_{SS}/K_{0}$ drops with increasing $\ell/R_{0}$ for all 3 choices of $L_{D}/L_{E}$.

4.3 Short crack limit

We now turn our attention to crack advance from a short pre-crack, for which $a/R_{0}<1$, recall Fig. 4a. Such cracks commonly arise at grain boundaries, at cracked carbide particles or as machining damage in structural alloys. The cracks are sufficiently short for no $K$-field (or $J$-field) to exist and are accompanied by plastic collapse at the structural level. Thus, failure occurs at a stress level somewhat above the yield strength, and the question of interest becomes: what is the dependence of macroscopic failure strain (below the necking strain) on $a/R_{0}$ and $\ell/R_{0}$?

First consider the case of a short crack of length $a/R_{0}=0.38$. Dimensional analysis implies,

[TABLE]

where $\sigma^{\infty}$ is the macroscopic remote stress on a tensile specimen (recall Fig. 5) and $\varepsilon^{\infty}$ is the work-conjugate remote tensile strain. A series of finite element simulations have been performed for $\hat{\sigma}/\sigma_{Y}=5$ and $N=0.1$, for illustrative purposes. The $\sigma^{\infty}$ versus $\varepsilon^{\infty}$ response is given in Fig. 10a for selected values of $\ell/R_{0}$ in the range 0 to 0.03. The tensile response is very sensitive to the choice of $\ell/R_{0}$, as follows. For $\ell/R_{0}=0$, the tensile response is almost identical to the material stress versus strain curve, and no failure is predicted. In contrast, the failure strain drops to about 1% when plastic strain gradients are accounted for. This is emphasized by the plot of failure strain $\varepsilon_{f}$ versus $\ell/R_{0}$ in Fig. 10b: $\varepsilon_{f}$ drops steeply from $\varepsilon_{f}/\varepsilon_{y}=3.7$ at $\ell/R_{0}=0.007$ to $\varepsilon_{f}/\varepsilon_{y}=1.45$ at $\ell/R_{0}=0.09$. Thus, strain gradient plasticity theory, along with a cohesive zone model, gives mechanistic insight into the drop in ductility when the fracture length scale $R_{0}$ drops (e.g., due to embrittlement) in relation to the plasticity length scale $\ell$.

In order to interrogate the source of the dramatic drop in ductility with increasing $\ell$ we examine the stress field ahead of a stationary short crack, absent the cohesive zone. The tensile stress $\sigma_{yy}$ is plotted as a function of $r/a$ in Fig. 11 for a fixed value of $\varepsilon^{\infty}=0.005$ such that material remote from the crack tip has fully yielded. Results are shown for $\ell/a=0,0.1$ and it is clear that the asymptotic stress field is similar to that for the long crack, as shown in Fig. 1. For the strain gradient solid, an elastic zone of extent on the order of $\ell$ exists at the crack tip. A crack tip $K$-field is evident, as for the long crack case, and it is this feature that results in the drop in ductility for the growing crack case of Fig. 10.

It remains to explore the dependence of failure strength $\sigma_{f}/\sigma_{Y}$ upon crack length $a/R_{0}$. It is anticipated that, for sufficiently large $a/R_{0}$, small scale yielding applies and failure occurs at $K=K_{SS}$ for a long pre-crack, such that $\sigma_{f}\approx K_{SS}/\sqrt{\pi a}$. With diminishing crack length, $\sigma_{f}/\sigma_{Y}$ rises until, for sufficiently small pre-cracks ($a/R_{0}<25$) the $K$-field ceases to exist and a $J$-analysis is necessary for a fracture mechanics assessment. A further reduction in $a/R_{0}$ leads to the short crack regime, and the full trajectory of $a/R_{0}$ is labelled as transition in Fig. 4a. The above qualitative discussion is now made precise by a series of calculations for selected values of $a/R_{0}$.

The predicted failure strength $\sigma_{f}/\sigma_{Y}$ is plotted as a function of $a/R_{0}$ in Fig. 12a for the case of fixed $\ell/R_{0}=0.02$. As expected, $\sigma_{f}/\sigma_{Y}$ increases from the small scale yielding value to the plastic collapse value $\sigma_{f}/\sigma_{Y}\approx 1$ with diminishing $a/R_{0}$. The regimes of validity of $K$ and $J$ are shown for completeness. A transition crack length can be identified by equating the fracture strength from plastic collapse theory $\sigma_{f}=\sigma_{Y}$ to the fracture strength from $K=K_{SS}$; such that $\sigma_{Y}\sqrt{\pi a_{T}}=K_{SS}$. Thus,

[TABLE]

The dependence of $a_{T}/R_{0}$ upon $\ell/R_{0}$ is plotted in Fig. 12b upon making use of (35). Upon recalling (27), the relation (35) reduces to

[TABLE]

Thus, the sensitivity of $a_{T}/R_{0}$ to $\ell/R_{0}$ arises directly from the dependence of $K_{SS}/K_{0}$ upon $\ell/R_{0}$.

5 Conclusions

The current study highlights the role of plastic strain gradients in influencing the R-curve for a long crack under small scale yielding and the tensile response in the presence of a short crack. An asymptotic analysis of the elastic-plastic stress state at the tip of a stationary crack in a strain gradient solid reveals that an elastic zone is present in the immediate vicinity of the crack tip. Consequently, the tensile stress immediately ahead of the crack tip displays an inverse square root singularity, in contrast to the HRR field of a conventional solid. This has immediate implications for a cohesive zone analysis of a growing crack: crack advance is predicted for cohesive strengths much greater than the yield strength. These predictions are consistent with observations of quasi-cleavage fracture with limited plasticity (Elssner et al., 1994; Bagchi and Evans, 1996; Korn et al., 2002).

Our study also reveals that the elastic crack tip singularity persists for the short crack case. Consequently, the tensile stress ahead of the short crack can far exceed the yield strength and overcome the cohesive strength of a cohesive zone placed at the crack tip. In turn, this leads to a significant drop in tensile ductility.

6 Acknowledgments

The authors would like to acknowledge the funding and technical support from BP (ICAM02ex) through the BP International Centre for Advanced Materials (BP-ICAM). E. Martínez-Pañeda also acknowledges financial support from the Ministry of Economy and Competitiveness of Spain through grant MAT2014-58738-C3 and the People Programme (Marie Curie Actions) of the European Union’s Seventh Framework Programme (FP7/2007-2013) under REA grant agreement n*∘* 609405 (COFUNDPostdocDTU). C.F. Niordson acknowledges support from the Danish Council for Independent Research through the research project “Advanced Damage Models with InTrinsic Size Effects" (Grant no: DFF-7017-00121). N.A. Fleck and V.S. Deshpande are grateful for additional financial support from the European Research Council in the form of an Advance Grant (MULTILAT, 669764).