Investigation of the stress distribution around a mode 1 crack with a novel strain gradient theory

Stress concentrations at the tip of a sharp crack have extensively been investigated in the past century. According to the calculations of Inglis, the stress ahead of a mode 1 crack shows the characteristics of a singularity. This solution is exact in the framework of linear elastic fracture mechanics (LEFM). From the viewpoint of multiscale modelling, however, it is evident that the stress at the tip of a stable crack cannot be infinite, because the strengths of atomic bonds are finite. In order to prevent the problem of this singularity, a new version of strain gradient elasticity is employed here. This theory is implemented in the commercial FEM code ABAQUS through user subroutine UEL. Convergence of the model is proved through consecutive mesh refinement. In consequence, the stresses ahead of a mode 1 crack become finite. Furthermore, the model predicts a size effect in the sense “smaller is stronger”.


Introduction
It is now more than 100 years ago that Inglis [1] presented his celebrated calculation of the stress field of an elliptical hole which takes on the shape of a sharp crack when the minor axis decreases relative to the major axis. At the tip of a sharp crack, the stress concentration factor led to a stress singularity. The solution of Inglis is exact within linear elastic fracture mechanics (LEFM). Or, to put it in other words, the solution is restricted to the case where Hooke's law is generally valid.
However, even in the elastic deformation regime deviations from Hooke's law are nowadays accepted. In the 1960s, theories of couple stresses and strain gradient elasticity were introduced [2][3][4][5][6]. In fact, the appearance of strain gradients is related to a bending of the microstructure. But since the microstructure tries to preserve its integrity, the material reacts with stiffness against bending. In consequence, couple-stresses are generated which appear in addition to the stresses following from Hooke's law. Thereby, the couple-stresses are connected to their corresponding components of the strain gradient tensor. The versions of strain gradient elasticity mentioned all have in common that the magnitude of stress concentrations is in general reduced. For the case of a circular hole in a field of simple tension Mindlin [5] obtained stress concentration factors depending on Poisson ratio and size of the hole. In the macroscopic limit for large holes, the result of strain gradient theory coincides with the classical result of Inglis. But for small holes, the stress concentration is reduced.
When the strain gradient theory of Toupin and Mindlin is applied to a sharp crack in a stress field, however, the stress singularity at the crack tip is still present. In spite of the general tendency to reduce stress concentrations, the strain gradient effects seem to be incapable of reducing the stress concentration at the tip of a sharp crack to a finite value [8]. Nevertheless, the aim of removing stress singularities seems to be an important issue of elasticity theory, and several researchers are still convinced that singularities of strain should somehow be removed at least for cases of static equilibrium [9,10].
In the present study, we therefore employ a new version of strain gradient elasticity which consequently leads to a removal of the stress singularity at the tip of a mode 1 crack. The main difference of the present approach compared to the older versions of strain gradient elasticity is that we consequently use a symmetric stress tensor throughout the model calculations, while Toupin and Mindlin postulated an asymmetric stress tensor.

Constitutive equation of isotropic strain gradient elasticity
The fundamental equation of our theory is an equation for the internal energy, whereby the deformation regime is restricted to small deformations. The model is reversible elastic. Thus, the internal energy writes as [11,12] (1) where summation is carried out over repeated indices. σ iJ is the Cauchy stress and ε iJ is the linear strain. η ijk are the strain gradients δ 2 u i /δx J δx k , where u i are the displacements and x i are the coordinates. μ ijk are higher order stresses related to the strain gradients. In stable equilibrium, the energy of the system takes a minimum. Per definition, our theory is linear. The relation between higher order stresses μ ijk and strain gradients η ijk reads as (2) where B is called the bending modulus. It should be noticed here that the energy density (3) related to the strain gradients is invariant with respect to rotations. Hence, the requirements for an isotropic material model are fulfilled. An even more general version of this theory was suggested by the present authors in reference [13], where B was consequently treated as sixth rank tensor leading to five independent bending moduli. The present version is a special case of the general theory. This means, we will here demonstrate that the consideration of one bending modulus is already sufficient to suppress the stress singularities at the tip of a mode 1 crack.

Finite Element implementation in ABAQUS through User subroutine UEL
The Finite Element implementation of this theory has already been explained elsewhere [12]. Therefore, only the fundamental principles of this approach are explained here. We have defined square nine node elements consisting of isoparametric four node sub-elements as depicted in figure 1.
The values for the nodal forces are derived from a principle of virtual work which requires that the work done by external forces equals the energy stored in the volume element.
In strain gradient elasticity neighbouring elements must fulfill boundary conditions in order to achieve C 1 continuity. Otherwise strain gradients would occur mainly at the element boundaries, where their energy would be neglected. In order to avoid this problem, we have introduced an overlapping mesh technique, which is considered as the most precise approach to handle this issue.

Definition of an intrinsic material length scale
Owing to the fact that Young's modulus E and bending modulus B have different dimension, a size effect of strain gradient elasticity is expected. This behavior may be described by the definition of an intrinsic material length l 0 [ 6]. Samples or structures with size larger than l 0 will behave similar as bulk material, and Hooke's law will provide a reasonable good approximation of the material properties. On the other hand, structures with size of l0 or below are dominated by a size effect in the sense smaller is stronger. In order to define the material length l0, we compare the contributions to the elastic energy arising either from strains or from strain gradients. Let us therefore assume a quadratic volume element with constant strain gradient ^1 12 = ^m. The strain in the middle of this volume element is assumed to be 0. The elastic energy of this square with side length l may be calculated using equations (1) and (2). Thus, one gets for the elastic energy of the strains The stresses occurring in equations (4a, b) are obtained from the condition ε yy = 0 in combination with the plane stress approximation. Taking into account η 112 = η 21 = η, one obtains (5) where v is the Poisson ratio. For the same volume element the elastic energy related to the strain gradients is The total energy of this volume element is W trran + W g . Now we are in the position to define the intrinsic material length l 0 by the condition that the energies W strain and W sg shall be equal for the square with side length l0. This condition yields (7)

Results
Now, the theory is applied to a crack in a transverse field of uniform uniaxial tension. According to the definition of a mode 1 crack, one should consider an infinite sample. But since the simulation of an infinite sample leads to numerical difficulties, we will here present a series of simulations where length and width of the sample are twice the crack length. Afterwards, it is demonstrated how an increase of sample size relative to the crack length influences the results for the stress at the crack tip. We here assume a Young's modulus E of 68 GPa and a Poisson ratio v of 0.33. These are values for aluminum, which seems to be a good example for a solid with nearly isotropic elasticity. Unfortunately, an experimental value for the bending modulus B of aluminum is not known. But when the material length l 0 is used as unit of the length scale, the value of the bending modulus may be obtained from equation (7). Therefrom, we get B = 4.245 [GPa l 0 2 ]. In figure 2, one can see a plot of the von Mises stress for a sample exposed to a tensile stress of 1 MPa with a crack length equal to l 0 . For this simulation a rather coarse mesh consisting of 40500 nodes was used. The maximum of the von Mises stress appeared at the crack tip where a value of 1.2728 MPa was found.
In order to check the mesh convergence of the solution another two simulations were performed with refined meshes using 251001 and 1002001 nodes for a quarter of the sample, respectively. Due to the mesh refinement, the von Mises stress at the crack tip was increased to 1.2867 MPa in the second and to 1.2885 MPa in the third simulation. This is only a small increase in stress compared to a reduction of the element length to 1 /5 and 1 /10 of the first mesh, respectively. The convergence norm reflecting the largest imbalance of nodal forces divided by the average of nodal forces was better than 10 -6 for each of the 3 simulations. A plot of the stress distribution using the finest mesh is depicted in figure 3. For comparison, a simulation of the model without strain gradient effects is shown in figure 4. Here, the maximum of von Mises stress is 23.42 MPa. This plot shows the typical characteristics of a stress singularity, where the maximum value for the stress observed in the plot is governed by the element size.
Let us now investigate how the crack length influences the results for the stress distribution within strain gradient elasticity. Simulations were performed for crack lengths between l /10 and 100 l, whereby the proportions of crack length to sample length and width remained the same, respectively. The results of this comparison are shown in table 1. They indicate a size effect in the sense smaller is stronger. In the case of short micro-cracks stress concentrations are drastically reduced, and in the limit to infinitesimal cracks they disappear completely.    In general, an increase of the sample size relative to the crack length slightly reduces the stress value at the crack tip, since the additional material increases the stiffness of the structure.

Summary and conclusions
A new approach of strain gradient elasticity was employed to study the stress distribution around a mode 1 crack. Thereby, the stress singularity at the crack tip has been removed, and a pronounced size effect has been found. The stress values at the tips of short cracks are lower compared to stresses observed at longer cracks.