On the relation between the energy of a distorted crystal lattice and the bending modulus of strain gradient elasticity

It is a well known fact that linear elastic fracture mechanics (LEFM) predicts stress singularities at the tips of sharp cracks, at sharp edges, at corners and at the surface of material transitions. However, from the viewpoint of the strengths of atomic bonds it is clear that only finite forces may be present at the tip of a stable crack. Therefore, theories of strain gradient elasticity were developed which reduce the values of stress concentrations. Within these theories a bending modulus is postulated which introduces an increased stiffness of the crystal lattice against bending. In the present study, the value of this bending modulus is evaluated on the basis of the electrostatic energy of a bent crystal lattice. This is done for the face centred cubic structure of NaCl. In fact, results for the bending modulus could be obtained although they depend on the crystal size.


Introduction
Nowadays, stress singularities are commonplace in linear elastic fracture mechanics. Inglis [1] was the first to evaluate the stress ahead of a sharp crack in a transverse stress field. This calculation of Inglis is said to be exact within linear elastic fracture mechanics (LEFM). Since then singularities of stress were found by other authors on numerous occasions. On the other hand, there are others who believe that singularities of stress or strain should be avoided [2][3][4][5] within continuum mechanics. But this can only be achieved through a change of the physical assumptions made in Hooke's law. One approach of this kind was proposed by Mindlin [6] who assumed an elastic bending stiffness in addition to the elastic response described in Hooke's law. Thereby, the bending stiffness was related to a bending modulus which has the dimension of a force. In fact, the theory of Mindlin could reduce the values of stress concentrations in some cases. But unfortunately, the theory could not reduce the stresses around sharp cracks to finite values. Nevertheless, the stress singularities can successfully be removed after a modification of Mindlin's theory, which was recently proposed by the authors of the present article [7,8]. The main idea behind this approach is that the asymmetric stress tensor used by Mindlin is replaced by a symmetric stress tensor. This version of the theory, which was implemented in the commercial FEM code ABAQUS through user subroutine UEL, has successfully been applied to several cases of stress concentrations [7][8][9]. Thereby, stresses singularities were consequently removed.
However, one disadvantage of strain gradient elasticity seems to be that reliable values for the bending modulus are still not known. In this context it should be said that nowadays there exist several versions of strain gradient elasticity which can include different numbers of bending moduli. In the case of isotropic materials, the most general version includes 5 independent bending moduli, while simplified versions may be formulated using only one bending modulus. In the anisotropic case the calculations become very complicated. Mindlin [10] reported that in his version of the theory there are all together 903 independent material constants including the usual 21 elastic constant of the triclinic system. The large number of constants makes it difficult to measure their values. Until today, precise values for the bending moduli are not known. Therefore, the present article is devoted to a numerical computation of the bending modulus. The computation is performed for the NaCl crystal structure, because the electrostatic energy of this crystal lattice may be computed in straightforward manner. The computations reveal that the value of the bending modulus is governed by long range electrostatic interactions.

Assumptions at the level of continuum theory
The constitutive equation for the energy density w of the strain gradient theory used in the present investigation writes as [7] (1) where summation is carried out over repeated indices. c iJ is the Cauchy stress and s 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. The relation between higher order stresses μ ijk and strain gradients η ijk reads as (2) where B is called the bending modulus. Further, the contribution of higher order stresses to the internal energy may be rewritten as (3) It should be mentioned here that also more general types of strain gradient elasticity have been developed. In reference [8] an energy density was used which may be considered as the most general isotropic invariant expression. Therefrom, a sixth rank bending tensor was derived which has 5 independent components. The model used in the present investigation is a special case of this general theory. Aside from the isotropic models also anisotropic versions of strain gradient elasticity have been developed [10]. But for sake of simplicity, we will here concentrate on a model with only one bending modulus.

Energy of the deformed crystal lattice
There are several contributions to the internal energy of a bent crystal lattice. First, one has to consider the cohesive energy of the material, second there are the ordinary strains and third we get the contribution of the additional bending stiffness. All these contributions are a result of the electrostatic interactions between the atoms together with repulsive forces originating from quantum mechanical principles or from the kinetic motion of atoms.

Inter-atomic potential of Na+ and Clions and the electrostatic energy
For alkali halides the inter-atomic potential of two atoms may be written in the form [11] (4) where we can distinguish between short range and long range interactions. When the coefficients of this equation are known for some material, then one may determine the value of the elastic constants through molecular dynamic simulations. However, in the present study we are more interested in the value for the bending modulus. In this respect long range interactions play an important role. We will therefore try to extract the value of the bending modulus from the electrostatic energy of a bent crystal lattice. Thereby, the cohesive energy of the undeformed crystal and the energy of the ordinary strains are to be subtracted from the overall energy in order to obtain the energy contribution of bending. But in contrast to molecular dynamic simulations, we will here neglect the influence of the atomic motion. We claim that the repulsive forces caused by the motion of atoms are very similar in homogeneously strained and in bent crystals. Therefore, we here concentrate on the potential energy instead of simulating the atomic motion. Thereby, the amount of computation time is drastically reduced.

Comparison of homogeneously strained and bent crystals
When the energies of homogeneously strained and bent crystals are compared, a size effect must be taken into account. Therefore, only samples with equal numbers of atoms should be compared. Moreover, the strain of a bent crystal can only be represented by an ensemble of homogeneously strained samples. Consequently, the bent crystal is divided into fibers of different lengths, and for each length a homogeneously strained crystal is constructed which has the same length as the fiber. Thereby, the cross contraction is considered in either case: in the bent crystal as well as in the homogeneously strained samples. Thus, the electrostatic energy of the homogeneously strained crystal corresponding to the strain of the bent crystal is obtained from the ensemble average representing all fiber lengths. The atomic structure of the bent crystal is depicted in figure 1, schematically. A neutral fiber of the same length as the undeformed sample is assumed in the middle of the cross section. Cross contractions were considered according to a Poisson ratio of v = 0.252. The lattice constant of NaCl is a = 0.564 nm. Computations were performed for a two-dimensional and for a three-dimensional model. But due to a size effect, only the three-dimensional model seems to be relevant for the real material properties. One derives the contribution of strain gradients to the internal energy through the equations (6) where the strain gradients j ijk are inserted according to their geometrical values used in the computation. This means, we have used the equation points exactly along the direction of the y-axis. Since the derivatives with respect to δxδy and δyδx must be equal, one gets η 112 = η 121 . Furthermore, we have assumed that η 222 = v η 112 , because of the cross contraction. The remaining components of the strain gradient tensor were assumed to be zero. During the computation it has been verified that the values which we get from equation (6) for the bending modulus are independent from the radius of curvature chosen for the samples. But since our theory of strain gradient elasticity is restricted to small deformations, the radii of curvature in the computation were always much larger than the sample size. The computation was repeated for different sample sizes. The samples always had cubic shape, and the largest sample included 64000 atoms. It was found that the value obtained for the bending modulus strongly depends on the sample size. Figure 2 shows the results for the energy E sg per atom in units of e 2 /πε 0 a , where e is the elementary charge, s 0 is the vacuum permittivity and a is the lattice constant. The lattice constant of NaCl is twice the atomic distance between neighbouring ions of Na and Cl. In this computation the radius of curvature for the neutral fibre was chosen to be 1000 times the lattice constant. The charges of the Na and Cl ions were assumed to be e and -e, respectively. The maximum number of atoms considered in the model corresponds to a side length of 40 atoms. Obviously, the computation shows a pronounced size effect. Unfortunately, the results obtained from the numerical calculation are only available for microscopic sample sizes due to the large computation time of the algorithm. In this context it should be noticed that the electrostatic energy of the homogeneously strained crystal was derived through an ensemble average over up to 40 replicas of the model representing the different strain values of the bent crystal. However, we are here interested in the model results which may be expected for mesoscopic sample sizes. The results show a strong dependence on the size of the model due to long range interactions.

Results
We therefore extrapolate the results in order to estimate the bending modulus of a crystal with side length of 10 μm. Hence, we replot figure 2 to derive an estimate for larger samples: In fact, figure 3 2nd

Discussion
It must be admitted that the extrapolation of the calculated value for the bending modulus to mesoscopic scale is somewhat speculative, because we have evaluated only a fraction of the desired sample size. On the other hand, it must be said that the computation is very time consuming due to the large number of summands occurring in equation (5). Within molecular dynamic simulations such problems are usually handled by introducing a cutoff radius for long range interactions. But such a method cannot be used here, since the value of the bending modulus mainly depends on long range forces. Nevertheless, the computation led to the result that there exists some bending stiffness in excess of that stiffness what would be expected from Hooke's law alone. The main purpose of the inclusion of this bending stiffness in strain gradient elasticity is the removal of stress singularities at sharp edges, corners, material transitions and crack tips [7][8][9]. In this context it is sufficient that there exists some bending modulus B > 0. Insofar, the present investigation may be viewed as justification of the strain gradient approach even though a value for the bending modulus on mesoscopic level could only be estimated.
However, the size effect predicted by strain gradient elasticity needs to be reconsidered. So far, it has been assumed that a size independent bending modulus will introduce a strong size effect in the sense smaller is stronger. On the other hand, the present study shows that the bending modulus itself is size dependent. This may lead to competing interactions with various types of size effects in different materials.
The geometrical reason for the size effect observed in the computation seems to be that crystallographic layers which are further away from the middle show a larger missorientation of crystallographic directions due to bending of the structure. But for crystals with high dislocation density or high concentration of other lattice defects the present calculations might not be accurate.

Summary and conclusions
We have here elaborated a relationship between the theory of strain gradient elasticity and the electrostatic energy of a distorted crystal lattice. The results obtained here may be seen as a justification for strain gradient theory. Moreover, a quantitative value for the bending modulus could be estimated. The evaluations led to the conclusion that this value is size dependent or at least grain size dependent. This discovery came unexpected, since earlier publications assumed bending moduli which are size independent constants [12]. However, the result concerning size dependence is only valid if the concentration of crystallographic defects is not too high. If the periodicity of the lattice is disturbed, then long range interactions will occur randomly distributed. In consequence, the value for the bending modulus could saturate already at a microscopic length scale. Nevertheless, the existence of a bending modulus B > 0 leads to the consequence that the stress singularities of continuum mechanics are suppressed.
In this investigation we have used a simple version of strain gradient elasticity which uses only one bending modulus. For the future it will be interesting to extend such evaluations to more general versions of this theory including more independent bending moduli and also anisotropic behaviour. Furthermore, it seems eligible to perform such calculations for larger samples.