The Green tensor of Mindlin's anisotropic first strain gradient elasticity

TL;DR

This paper derives a non-singular Green tensor for Mindlin's anisotropic strain gradient elasticity, applicable to arbitrary anisotropic materials, and demonstrates its numerical implementation and validation against molecular statics results.

Contribution

It provides the first general derivation of the Green tensor for anisotropic first strain gradient elasticity, including numerical implementation and validation.

Findings

Green tensor is non-singular at the origin.

The tensor converges to classical form away from the origin.

Numerical results agree with molecular statics calculations.

Abstract

We derive the Green tensor of Mindlin's anisotropic first strain gradient elasticity. The Green tensor is valid for arbitrary anisotropic materials, with up to 21 elastic constants and 171 gradient elastic constants in the general case of triclinic media. In contrast to its classical counterpart, the Green tensor is non-singular at the origin, and it converges to the classical tensor a few characteristic lengths away from the origin. Therefore, the Green tensor of Mindlin's first strain gradient elasticity can be regarded as a physical regularization of the classical anisotropic Green tensor. The isotropic Green tensor and other special cases are recovered as particular instances of the general anisotropic result. The Green tensor is implemented numerically and applied to the Kelvin problem with elastic constants determined from interatomic potentials. Results are compared to molecular statics calculations carried out with the same potentials.

Tables

Table 1. Table 1: Elastic and gradient-elastic constants obtained from the interatomic potentials Lee ( 2001 ) and Mendelev ( 2008 ) .

	Cu EAM	Cu MEAM	Al MEAM
$C_{1,1}$ $\left[\text{eV}/{\text{Å}}^3\right]$	1.0868	1.0994	7.1366$\cdot {10}^{-1}$
$C_{1,2}$ $\left[\text{eV}/{\text{Å}}^3\right]$	7.9386$\cdot {10}^{-1}$	7.7973$\cdot {10}^{-1}$	3.8649$\cdot {10}^{-1}$
$C_{4,4}$ $\left[\text{eV}/{\text{Å}}^3\right]$	5.2252$\cdot {10}^{-1}$	5.1043$\cdot {10}^{-1}$	1.9704$\cdot {10}^{-1}$
$D_{1,1}$ $\left[\text{eV}/\text{Å}\right]$	1.1182	6.5018$\cdot {10}^{-1}$	1.0855
$D_{1,2}$ $\left[\text{eV}/\text{Å}\right]$	3.5814$\cdot {10}^{-1}$	3.6659$\cdot {10}^{-1}$	1.4572$\cdot {10}^{-1}$
$D_{1,3}$ $\left[\text{eV}/\text{Å}\right]$	3.7951$\cdot {10}^{-1}$	2.4150$\cdot {10}^{-1}$	1.5934$\cdot {10}^{-1}$
$D_{2,2}$ $\left[\text{eV}/\text{Å}\right]$	4.7935$\cdot {10}^{-1}$	7.3885$\cdot {10}^{-1}$	8.4221$\cdot {10}^{-1}$
$D_{2,3}$ $\left[\text{eV}/\text{Å}\right]$	3.0103$\cdot {10}^{-1}$	2.0651$\cdot {10}^{-1}$	1.5671$\cdot {10}^{-1}$
$D_{2,4}$ $\left[\text{eV}/\text{Å}\right]$	1.2789$\cdot {10}^{-1}$	4.7496$\cdot {10}^{-1}$	7.1708$\cdot {10}^{-1}$
$D_{2,5}$ $\left[\text{eV}/\text{Å}\right]$	1.0652$\cdot {10}^{-1}$	-4.2545$\cdot {10}^{-2}$	-1.1434$\cdot {10}^{-2}$
$D_{3,3}$ $\left[\text{eV}/\text{Å}\right]$	4.3662$\cdot {10}^{-1}$	2.9055$\cdot {10}^{-1}$	2.7613$\cdot {10}^{-1}$
$D_{3,5}$ $\left[\text{eV}/\text{Å}\right]$	1.2789$\cdot {10}^{-1}$	-1.8275$\cdot {10}^{-2}$	-1.2408$\cdot {10}^{-1}$
$D_{16,16}$ $\left[\text{eV}/\text{Å}\right]$	1.4925$\cdot {10}^{-1}$	3.7419$\cdot {10}^{-2}$	1.6786$\cdot {10}^{-1}$
$D_{16,17}$ $\left[\text{eV}/\text{Å}\right]$	1.0652$\cdot {10}^{-1}$	3.7394$\cdot {10}^{-2}$	1.5006$\cdot {10}^{-1}$

Equations

$${\partial }_j{u}_i={\beta }_{ij}+{\beta }_{ij}^{\ast }.$$

$$\mathcal{W}(\mathbf{e},\nabla \mathbf{e})=\frac{1}{2}{\mathbb{C}}_{ijkl}{e}_{ij}{e}_{kl}+\frac{1}{2}{\mathbb{D}}_{ijmkln}{\partial }_m{e}_{ij}{\partial }_n{e}_{kl}.$$

$$e_{ij}=\frac{1}{2}\left({\beta }_{ij}+{\beta }_{ji}\right),$$

$${\mathbb{C}}_{ijkl}={\mathbb{C}}_{jikl}={\mathbb{C}}_{ijlk}={\mathbb{C}}_{klij},$$

$${\mathbb{D}}_{ijmkln}={\mathbb{D}}_{jimkln}={\mathbb{D}}_{ijmlkn}={\mathbb{D}}_{klnijm}.$$

$${\sigma }_{ij}$$

$${\tau }_{ijm}$$

$$\mathcal{L}=-\mathcal{W}-\mathcal{V}=-\frac{1}{2}\left({\mathbb{C}}_{ijkl}{\beta }_{ij}{\beta }_{kl}+{\mathbb{D}}_{ijmkln}{\beta }_{ij,m}{\beta }_{kl,n}\right)+u_i{b}_i,$$

$$\mathcal{V}=-u_i{b}_i$$

$$\frac{\delta \mathcal{L}}{\delta u_i}=\frac{\partial \mathcal{L}}{\partial u_i}-{\partial }_j\frac{\partial \mathcal{L}}{\partial ({\partial }_j{u}_i)}+{\partial }_k{\partial }_j\frac{\partial \mathcal{L}}{\partial ({\partial }_k{\partial }_j{u}_i)}=0.$$

$$\displaystyle\partial_{j}\big{(}\sigma_{ij}-\partial_{m}\tau_{ijm}\big{)}+b_{i}=0\,.$$

$$L_{ik}^{\text{M}}{u}_k+f_i=0.$$

$$L_{ik}^{\text{M}}={\mathbb{C}}_{ijkl}{\partial }_j{\partial }_l-{\mathbb{D}}_{ijmkln}{\partial }_j{\partial }_l{\partial }_m{\partial }_n,$$

$$f_i=b_i-\left[{\mathbb{C}}_{ijkl}{\partial }_j-{\mathbb{D}}_{ijmkln}{\partial }_j{\partial }_m{\partial }_n\right]{\beta }_{kl}^{\ast }$$

$$u_k=G_{kj}\ast f_j,$$

$$\left[{\mathbb{C}}_{ijkl}{\partial }_j{\partial }_l-{\mathbb{D}}_{ijmkln}{\partial }_j{\partial }_l{\partial }_m{\partial }_n\right]G_{km}+{\delta }_{im}\delta =0.$$

$$\hat{f}(\mathbf{k})$$

$$f(\mathbf{x})$$

$$f(\mathbf{x})=\frac{1}{(2\pi {)}^3}{\int }_{{\mathbb{R}}^3}\hat{f}(\mathbf{k})\cos \left(\mathbf{k}\cdot \mathbf{x}\right)\text{d}\hat{V}.$$

$$\left[{\mathcal{C}}_{ik}(\mathbf{k})+{\mathcal{D}}_{ik}(\mathbf{k})\right]{\hat{G}}_{kj}(\mathbf{k})={\delta }_{ij},$$

$${\mathcal{C}}_{ik}(\mathbf{k})$$

$${\mathcal{D}}_{ik}(\mathbf{k})$$

$$\mathbf{\kappa }=\frac{\mathbf{k}}{k},k=\sqrt{k_i{k}_i},{\mathbf{\kappa }}^2=1,$$

$$k^2\left[{\mathcal{C}}_{ik}(\mathbf{\kappa })+k^2{\mathcal{D}}_{ik}(\mathbf{\kappa })\right]{\hat{G}}_{kj}(\mathbf{k})={\delta }_{ij},$$

$$k^2\left[\mathcal{C}(\mathbf{\kappa })+k^2\mathcal{D}(\mathbf{\kappa })\right]\hat{\mathbf{G}}(\mathbf{k})=\mathbf{I}.$$

$$\hat{\mathbf{G}}(\mathbf{k})=\frac{{\left[\mathcal{C}(\mathbf{\kappa })+k^2\mathcal{D}(\mathbf{\kappa })\right]}^{-1}}{k^2}.$$

$$\mathbf{G}(\mathbf{x})$$

$$=\frac{1}{8{\pi }^3}{\int }_{\mathcal{S}}{\int }_0^{\infty }{\left[\mathcal{C}(\mathbf{\kappa })+k^2\mathcal{D}(\mathbf{\kappa })\right]}^{-1}\cos \left(k\mathbf{\kappa }\cdot \mathbf{x}\right)\text{d}k\text{d}\omega .$$

$$\mathcal{C}(\mathbf{\kappa })=\mathbf{R}(\mathbf{\kappa }){\mathbf{V}}^2(\mathbf{\kappa }){\mathbf{R}}^T(\mathbf{\kappa }),$$

$${\mathcal{C}}^{\frac{1}{2}}=\mathbf{R}(\mathbf{\kappa })\mathbf{V}(\mathbf{\kappa }){\mathbf{R}}^T(\mathbf{\kappa })$$

$$\mathcal{C}+k^2\mathcal{D}(\mathbf{\kappa })={\mathcal{C}}^{\frac{1}{2}}\left[\mathbf{I}+k^2{\mathbf{\Lambda }}^2(\mathbf{\kappa })\right]{\mathcal{C}}^{\frac{1}{2}},$$

$${\mathbf{\Lambda }}^2(\mathbf{\kappa })={\mathcal{C}}^{-\frac{1}{2}}(\mathbf{\kappa })\mathcal{D}(\mathbf{\kappa }){\mathcal{C}}^{-\frac{1}{2}}(\mathbf{\kappa })$$

Full text

††journal: Materials Theory

1 Introduction

Green functions are objects of fundamental importance in field theories, since they represent the fundamental solution of linear inhomogeneous partial differential equations (PDEs) from which any particular solution can be obtained via convolution with the source term Green (1828). Moreover, Green functions are the basis of important numerical methods for boundary value problems, such as the boundary element method Becker (1992), and they provide “flexible” boundary conditions for atomistic simulations Trinkle (2008). In the context of linear elasticity, the Green function is a tensor-valued function of rank two, also known as the Green tensor. When contracted with a concentrated force acting at the origin, the Green tensor yields the displacement field in an infinite elastic medium. Lord Kelvin (1882) first derived the closed-form expression of the classical Green tensor for isotropic materials. For anisotropic materials, Lifshitz and Rosenzweig (1947) and Synge (1957) were able to derive the Green tensor in terms of an integral expression over the equatorial circle of the unit sphere in Fourier space. Barnett (1972) extended this result to the first two derivatives of the Green tensor, and showed that the line-integral representation is well suited for numerical integration (see also Bacon et al. (1979); Teodosiu (1982)).

The Green tensor and its derivatives are singular at the origin, ultimately as a consequence of the lack of intrinsic length scales in the classical theory of elasticity. The unphysical singularities in the elastic fields derived from the Green tensor hinder their applicability in nano-mechanics, including the elastic theory of defects such as cracks, dislocations and inclusions Mura (1987); Askes and Aifantis (2011). Generalized elastic field theories with intrinsic length scales have been proposed in the context of micro-continuum theories Eringen (1999), non-local theories Eringen (2002), and gradient theories Kröner (1963); Mindlin (1964, 1968, 1972); Mindlin and Eshel (1968). In particular, Mindlin’s anisotropic strain gradient elasticity has received renewed attention as a tool to solve engineering problems at the micro- and nano-scales for realistic materials Polizzotto (2018). Only recently, the structure of the gradient-elastic tensor has been rationalized for different material symmetry classes Auffray et al. (2013), and its atomistic representation and ensuing determination from interatomic potentials has become available Admal et al (2016).

The number of independent strain gradient elastic moduli ranges from 5 for isotropic materials, to 171 in the general case of triclinic materials. While simple expressions of the Green tensor exist for the isotropic case Rogula (1973); Lazar and Po (2018), and for simplified anisotropic theories Lazar and Po (2015, 2015), the Green tensor of the fully anisotropic theory of Mindlin’s strain gradient elasticity has remained so far a rather elusive object. Rogula (1973) provided an expression for the Green tensor in gradient elasticity of arbitrary order, which involves a sum of terms associated with the roots of a certain characteristic polynomial. However, such representation renders its numerical implementation rather impractical, and it conceals the mathematical structure of the Green tensor in relationship to its classical counterpart.

The objective of this paper is to derive a simple representation of the Green tensor of Mindlin’s anisotropic first strain gradient elasticity, whose integral kernel involves only matrix operations suitable for efficient numerical implementation. Following a brief summary of Mindlin’s anisotropic first strain gradient elasticity in section 2, we derive the matrix representation of the Green tensor in section 3. It is shown that the Green tensor is non-singular at the origin, while its first gradient is finite but discontinuous at the origin. The classical tail of the Green tensor, as well as its classical limit for vanishing gradient parameters are easily recovered from the non-singular expression. In section 4 we demonstrate that the Green tensor generalizes other expressions found in the literature. In section 5 we consider the Kelvin problem and compare the prediction of the Green tensor to atomistic calculations.

2 Mindlin’s anisotropic gradient elasticity

Let us consider an infinite elastic body in three-dimensional space and assume that the gradient of the displacement field $\bm{u}$ is additively decomposed into an elastic distortion tensor $\bm{\beta}$ and an inelastic111 The inelastic distortion comprises plastic effects, and is typically an incompatible field. When the inelastic distortion is absent the elastic distortion is compatible. eigen-distortion tensor $\bm{\beta}^{*}$:

[TABLE]

In the linearized theory of Mindlin’s form-II first strain gradient elasticity Mindlin (1964, 1968); Mindlin and Eshel (1968); Mindlin (1972), the strain energy density of an homogeneous and centrosymmetric222Due to the centrosymmetry, there is no coupling between $e_{ij}$ and $\partial_{m}e_{kl}$. material is given by

[TABLE]

The strain energy density (2) is a function of the infinitesimal elastic strain tensor

[TABLE]

and of its gradient $e_{ij,m}$. The tensor $\mathbb{C}$ is the standard rank-four tensor of elastic constants. By virtue of the symmetries

[TABLE]

it possesses up to 21 independent constants with units of $\text{eV}/\textup{\AA}^{3}$. The tensor $\mathbb{D}$ is the rank-six tensor of strain gradient elastic constants, with symmetries

[TABLE]

It has units of $\text{eV}/\textup{\AA}$. In the general case of triclinic materials the number of independent constants in the tensor $\mathbb{D}$ is equal to 171 Auffray et al. (2013).

The quantities conjugate to the elastic strain tensor and its gradient are the Cauchy stress tensor $\bm{\sigma}$ and the double stress tensor $\bm{\tau}$, respectively. These are defined as:

[TABLE]

In the presence of a body forces density $\bm{b}$, the static Lagrangian density of the system becomes:

[TABLE]

where

[TABLE]

is the potential of the body force. The condition of static equilibrium is expressed by the Euler-Lagrange equation

[TABLE]

In terms of the Cauchy and double stress tensors, Eq. (10) takes the following form Mindlin (1964):

[TABLE]

Using Eqs. (1) (6) (7), Eq. (11) can be cast in the following equation for displacements:

[TABLE]

In Eq. (12), $L^{\text{M}}_{ik}$ denotes the differential operator of Mindlin’s anisotropic first strain gradient elasticity

[TABLE]

while

[TABLE]

is the forcing term. Note that the second term on the right hand side of Eq. (14) is an “effective” internal force due to the inelastic eigen-distortion, and arises in the presence of material defects, such as inclusions, cracks, and dislocations. This term is the gradient version of the internal force in Mura’s eigen-strain theory Mura (1987).

3 The Green tensor of Mindlin’s first strain gradient elasticity

In this section, we derive the three-dimensional Green tensor of the operator (13). To this end, we seek the solution to Eq. (12) in the form

[TABLE]

where the symbol $*$ indicates convolution over the three-dimensional space, and $\bm{G}$ is the Green tensor of Mindlin’s anisotropic differential operator $\bm{L}^{M}$. Substituting Eq. (15) into Eq. (12), one finds that $\bm{G}$ satisfies the following inhomogeneous PDE:

[TABLE]

In Eq. (16), $\delta_{ij}$ is the Kronecker symbol, while $\delta$ is the three-dimensional Dirac $\delta$-distribution.

Taking the Fourier transform333 The Fourier transform and its inverse are defined as, respectively Vladimirov (1971):

$\displaystyle\hat{f}(\bm{k})$ $\displaystyle=\int_{\mathbb{R}^{3}}f(\bm{x})\,\text{e}^{-\text{i}\bm{k}\cdot\bm{x}}\,\text{d}V\,,$

(17)

$\displaystyle f(\bm{x})$ $\displaystyle=\frac{1}{(2\pi)^{3}}\int_{\mathbb{R}^{3}}\hat{f}(\bm{k})\,\text{e}^{\text{i}\bm{k}\cdot\bm{x}}\,\text{d}\hat{V}\,.$

(18)

For a real-valued function, the inverse Fourier transform is

$\displaystyle f(\bm{x})=\frac{1}{(2\pi)^{3}}\int_{\mathbb{R}^{3}}\hat{f}(\bm{k})\,\cos\left(\bm{k}\cdot\bm{x}\right)\,\text{d}\hat{V}\,.$

(19)

of Eq. (16), we obtain the following algebraic equation for the Green tensor $\hat{G}_{kj}(\bm{k})$ in Fourier space

[TABLE]

where

[TABLE]

are symmetric matrices. If we further define the unit vector in Fourier space

[TABLE]

then (20) becomes:

[TABLE]

or equivalently, in matrix notation,

[TABLE]

Stability of the differential operator $\bm{L}^{M}$ requires that the matrix $\bm{\mathcal{C}}(\bm{\kappa})+k^{2}\bm{\mathcal{D}}(\bm{\kappa})$ be positive definite. Since this requirement must hold for all $k$ and $\bm{\kappa}$, then the matrices $\bm{\mathcal{C}}(\bm{\kappa})$ and $\bm{\mathcal{D}}(\bm{\kappa})$ must be individually positive definite. Under the assumption that $\bm{\mathcal{C}}(\bm{\kappa})$ and $\bm{\mathcal{D}}(\bm{\kappa})$ are symmetric positive definite (SPD) matrices, the solution of (25) in Fourier space clearly reads:

[TABLE]

The three-dimensional Green tensor in real space is obtained by inverse Fourier transform of Eq. (26). It reads:

[TABLE]

In Eq. (27), $\text{d}\hat{V}=k^{2}\,\text{d}k\,\text{d}\omega$ indicates the volume element in Fourier space, and $\text{d}\omega$ is an elementary solid angle on the unit sphere $\mathcal{S}$. Our objective now is to obtain an alternative expression of the matrix inverse $[\bm{\mathcal{C}}(\bm{\kappa})+k^{2}\bm{\mathcal{D}}(\bm{\kappa})]^{-1}$ which allows us to carry out the the $k$-integral analytically. By doing so, the non-singular nature of the Green tensor at the origin is revealed. We start by observing that, by virtue of its SPD character, the matrix $\bm{\mathcal{C}}(\bm{\kappa})$ admits the following eigen-decomposition

[TABLE]

where $\bm{R}(\bm{\kappa})$ is the orthogonal matrix of the eigenvectors of $\bm{\mathcal{C}}(\bm{\kappa})$, while $\bm{V}^{2}(\bm{\kappa})$ is the diagonal matrix of positive eigenvalues of $\bm{\mathcal{C}}(\bm{\kappa})$. Moreover, the matrix

[TABLE]

is also SPD. Using (29), let us consider the following identity:

[TABLE]

where

[TABLE]

is a SPD matrix with units of length squared. With this decomposition, the Green tensor in Fourier space becomes

[TABLE]

while in real space we obtain

[TABLE]

In order to carry out the $k$-integral, we make use of the following matrix identity:444 The proof of (35) descends from the fact that $\bm{\Lambda}^{2}(\bm{\kappa})$ is a real SPD matrix, and therefore it admits the eigen-decomposition

\displaystyle{\color[rgb]{1,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,1}\pgfsys@color@cmyk@stroke{0}{1}{0}{0}\pgfsys@color@cmyk@fill{0}{1}{0}{0}\bm{\Lambda}^{2}(\bm{\kappa})={\bm{Q}}(\bm{\kappa})\bm{D}^{2}(\bm{\kappa}){\bm{Q}}^{T}(\bm{\kappa})\,,}

(34)

where $\bm{D}^{2}(\bm{\kappa})=\text{diag}\left\{\lambda^{2}_{i}(\bm{\kappa})\right\}$ is the diagonal matrix of the positive eigenvalues of $\bm{\Lambda}^{2}(\bm{\kappa})$, and {\color[rgb]{1,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,1}\pgfsys@color@cmyk@stroke{0}{1}{0}{0}\pgfsys@color@cmyk@fill{0}{1}{0}{0}{\bm{Q}}(\bm{\kappa})} is the orthogonal matrix of its eigenvectors. With this observation, we immediately obtain

$\displaystyle\int_{0}^{\infty}\left[\bm{I}+k^{2}\bm{\Lambda}^{2}(\bm{\kappa})\right]^{-1}\,\cos(k\bm{\kappa}\cdot\bm{x})\,\text{d}k$ $\displaystyle=\int_{0}^{\infty}\left[\bm{Q}(\bm{\kappa})\left(\bm{I}+k^{2}\bm{D}^{2}(\bm{\kappa})\right)\bm{Q}^{T}(\bm{\kappa})\right]^{-1}\,\cos(k\bm{\kappa}\cdot\bm{x})\,\text{d}k$

$\displaystyle=\bm{Q}(\bm{\kappa})\int_{0}^{\infty}\text{diag}\left\{\frac{\cos(k\bm{\kappa}\cdot\bm{x})}{1+k^{2}\lambda^{2}_{i}(\bm{\kappa})}\right\}\,\text{d}k\,\bm{Q}^{T}(\bm{\kappa})\,.$

With the help of the definite integral 3.767 in Gradshteyn and Ryzhik (2007), we obtain

$\displaystyle\int_{0}^{\infty}\left[\bm{I}+k^{2}\bm{\Lambda}^{2}(\bm{\kappa})\right]^{-1}\,\cos(k\bm{\kappa}\cdot\bm{x})\,\text{d}k$ $\displaystyle=\frac{\pi}{2}\,\bm{Q}(\bm{\kappa})\,\text{diag}\left\{\frac{\text{e}^{-|\bm{\kappa}\cdot\bm{x}|/\lambda_{i}(\bm{\kappa})}}{\lambda_{i}(\bm{\kappa})}\right\}\bm{Q}^{T}(\bm{\kappa})$

$\displaystyle=\frac{\pi}{2}\,\bm{Q}(\bm{\kappa})\,\text{diag}\left\{\text{e}^{-|\bm{\kappa}\cdot\bm{x}|/\lambda_{i}(\bm{\kappa})}\right\}\bm{D}^{-1}(\bm{\kappa})\,\bm{Q}^{T}(\bm{\kappa})$

$\displaystyle=\frac{\pi}{2}\,\bm{Q}(\bm{\kappa})\,\exp\left\{-|\bm{\kappa}\cdot\bm{x}|\,\bm{D}^{-1}(\bm{\kappa})\right\}\,\bm{Q}^{T}(\bm{\kappa})\bm{\Lambda}^{-1}(\bm{\kappa})$

$\displaystyle=\frac{\pi}{2}\,\exp\left\{-|\bm{\kappa}\cdot\bm{x}|\,\bm{\Lambda}^{-1}(\bm{\kappa})\right\}\,\bm{\Lambda}^{-1}(\bm{\kappa}).$

In the last step we have used the property that the matrix exponential is an isotropic tensor-valued function of its argument.

[TABLE]

With this identity, the Green tensor takes the form

[TABLE]

Next, Eq. (36) is further simplified noting that the integration kernel is an even function of $\bm{\kappa}$. Therefore, the integral over the unit sphere $\mathcal{S}$ is twice the integral over a hemisphere.

At the origin, any arbitrary hemisphere $\mathcal{H}$ can be chosen, and the Green tensor assumes the value

[TABLE]

This noteworthy result shows that the Green tensor is non-singular at the origin, in contrast to classical elasticity.

Away from the origin, we can choose the hemisphere having the direction $\bm{x}$ as the zenith. This is a convenient choice because all points $\bm{\kappa}$ on such a hemisphere satisfy the condition $\bm{\kappa}\cdot\bm{x}\geq 0$. This hemisphere can be parameterized by the zenith angle $\theta$ and the azimuth angle $\phi$, as shown in Fig. 1. In this reference system, the unit vector $\bm{\kappa}$ can be expressed as

[TABLE]

where $\hat{\bm{e}}_{3}=\bm{x}/x$. Finally, letting $q=\cos\theta$, the elementary solid angle becomes

[TABLE]

and

[TABLE]

Therefore the Green tensor of the anisotropic Mindlin differential operator of first order finally becomes

[TABLE]

3.1 The first two gradients of the Green tensor

The first two gradients of the Green tensor are computed directly by differentiation of (36). The first gradient reads

[TABLE]

In components this is:

[TABLE]

Note that, because of the presence of the sign function, the gradient of the Green tensor is finite but discontinuous at the origin. From a computational perspective, it is more convenient to express this result in reference system of Fig. 1. Doing so we find the alternative representation

[TABLE]

The second gradient of the Green tensor reads

[TABLE]

Letting $\bm{n}(\phi)=\bm{\kappa}(\pi/2,\phi)$ be a unit vector on the equatorial plane $\bm{\kappa}\cdot\bm{x}=0$, we finally obtain

[TABLE]

Note that the second gradient of the Green tensor is singular at the origin.

3.2 The classical limit

It is now shown that Green tensor (36) converges to the classical Green tensor $\bm{G}^{0}$ Lifshitz and Rosenzweig (1947); Synge (1957) when the field point $\bm{x}$ is sufficiently far from the origin compared to the characteristic length scales, that is when

[TABLE]

where $\lambda_{i}$ is an eigenvalue of $\bm{\Lambda}$, and $i=1,2,3$. This important property guarantees that the non-singular Green tensor (41) regularizes the classical anisotropic Green tensor in the far field. Moreover, as a special case satisfying condition (47), the classical Green tensor $\bm{G}^{0}$ is also recovered in the limit of vanishing tensor of strain gradient coefficients $\mathbb{D}$. The classical Green tensor $\bm{G}^{0}$ is readily recovered if we consider the limit555 Using the eigen-decomposition (34):

$\displaystyle\lim_{\||\bm{\kappa}\cdot\bm{x}|\,\bm{\Lambda}^{-1}\|\rightarrow\infty}\,\exp\left\{-|\bm{\kappa}\cdot\bm{x}|\,\bm{\Lambda}^{-1}(\bm{\kappa})\right\}\,\bm{\Lambda}^{-1}(\bm{\kappa})=$

$\displaystyle\lim_{\||\bm{\kappa}\cdot\bm{x}|\,\bm{D}^{-1}\|\rightarrow\infty}\,\bm{Q}(\bm{\kappa})\exp\left\{-|\bm{\kappa}\cdot\bm{x}|\,\bm{D}^{-1}(\bm{\kappa})\right\}\,\bm{D}^{-1}(\bm{\kappa})\bm{Q}^{T}(\bm{\kappa})=$

$\displaystyle\lim_{|\bm{\kappa}\cdot\bm{x}|/\lambda_{i}\rightarrow\infty}\bm{Q}(\bm{\kappa})\,\text{diag}\left\{\frac{\exp\left\{-|\bm{\kappa}\cdot\bm{x}|/\lambda_{i}(\bm{\kappa})\right\}}{\lambda_{i}(\bm{\kappa})}\right\}\,\bm{Q}^{T}(\bm{\kappa})=$

$\displaystyle\bm{Q}(\bm{\kappa})\,\frac{2\bm{I}}{x}\,\delta(\bm{\kappa}\cdot\hat{\bm{x}})\,\bm{Q}^{T}(\bm{\kappa})=\frac{2\bm{I}}{x}\,\delta(\bm{\kappa}\cdot\hat{\bm{x}})\,.$

[TABLE]

where $\hat{\bm{x}}=\bm{x}/x$ and $\bm{I}$ is the identity tensor. In fact, the substitution of (48) into (36) yields

[TABLE]

Here we used again the notation $\bm{n}(\phi)=\bm{\kappa}(\pi/2,\phi)$ to indicate a unit vector on the equatorial plane $\bm{\kappa}\cdot\bm{x}=0$. Note that the span of integration can be reduced to the range $0\leq\phi\leq\pi$ using the symmetry $\bm{\mathcal{C}}^{-1}(\bm{n})=\bm{\mathcal{C}}^{-1}(-\bm{n})$.

4 Special cases

In this section we show that the Green tensor (36) generalizes other results obtained in the literature.

4.1 The weakly non-local Green tensor $\mathbf{G}^{\text{NL}}$

Lazar and Po Lazar and Po (2015) have considered a simplified strain gradient elasticity theory under the assumption

[TABLE]

a framework which was named Mindlin’s strain gradient elasticity with weak non-locality because of its relation to non-local theories Lazar et al (2018); Lazar and Agiasofitou (2011). The Green tensor (36) recovers our previous result as a special case. In fact, under the previous assumption, we have

[TABLE]

and

[TABLE]

Therefore the Green tensor becomes

[TABLE]

which is the expression given in Lazar and Po (2015).

4.2 The Green tensor of anisotropic gradient elasticity of Helmholtz type $\mathbf{G}^{\text{H}}$

An even simpler theory, named Mindlin’s gradient elasticity of Helmholtz type, has been proposed by Lazar and Po (2015). The theory is characterized by only one gradient length scale parameter $\ell$, which renders the tensor $\bm{L}$ diagonal:

[TABLE]

The non-singular Green tensor of this theory is obtained by substituting (53) in (52), thus yielding

[TABLE]

which coincides with the expression given in Lazar and Po (2015).

4.3 The isotropic Green tensor $\mathbf{G}^{\text{I}}$

The isotropic tensor $\mathbb{C}$ has components

[TABLE]

where $\lambda$ and $\mu$ are the Lamé constants. On the other hand, the isotropic tensor $\mathbb{D}$ reads

[TABLE]

where $a_{1}$, $a_{2}$, $a_{3}$, $a_{4}$, $a_{5}$ are the gradient parameters in isotropic Mindlin’s first strain gradient elasticity theory Mindlin (1964) (see also Mindlin (1968); Lazar and Po (2018)). Therefore, the matrices $\bm{\mathcal{C}}(\bm{\kappa})$ and $\bm{\mathcal{D}}(\bm{\kappa})$ become, respectively

[TABLE]

The two characteristic lengths $\ell_{1}$ and $\ell_{2}$ introduced above are defined as

[TABLE]

Owing to the special structure666 Consider a matrix $\bm{A}$ with structure

$\displaystyle A_{ij}=a\kappa_{i}\kappa_{j}+b(\delta_{ij}-\kappa_{i}\kappa_{j})\,.$

(61)

If $a>b>0$, then the matrix is SPD, and a unique SPD square root of $A_{ij}$ exists with form

$\displaystyle A_{ij}^{\frac{1}{2}}$ $\displaystyle=\sqrt{a}\kappa_{i}\kappa_{j}+\sqrt{b}(\delta_{ij}+\kappa_{i}\kappa_{j})\,.$

(62)

Moreover, the inverse of $A_{ij}$ reads

$\displaystyle A_{ij}^{-1}$ $\displaystyle=\frac{1}{a}\kappa_{i}\kappa_{j}+\frac{1}{b}(\delta_{ij}-\kappa_{i}\kappa_{j})\,.$

(63)

of $\bm{\mathcal{C}}(\bm{\kappa})$ and $\bm{\mathcal{D}}(\bm{\kappa})$, the following results are easily obtained:

[TABLE]

The matrix $\bm{\Lambda}^{-1}$ admits the eigenvalue $1/\ell_{1}$, corresponding to the eigenvector $\hat{\bm{v}}_{1}=\bm{\kappa}$. The degenerate eigenvalue $1/\ell_{2}$ has multiplicity two, corresponding to two arbitrary eigenvectors $\hat{\bm{v}}_{2}$ and $\hat{\bm{v}}_{3}$ perpendicular to $\bm{\kappa}$. Choosing such eigenvectors to be mutually orthogonal, the matrix $\bm{\Lambda}^{-1}$ admits the eigen decomposition $\bm{\Lambda}^{-1}=\bm{Q}\bm{D}^{-1}\bm{Q}^{T}$. Here

[TABLE]

is an orthogonal matrix whose columns are the eigenvectors of $\bm{\Lambda}^{-1}$, and

[TABLE]

is the diagonal matrix of its eigenvalues. This special form of $\bm{Q}$ yields the identity

[TABLE]

Using these results in (36), we obtain

[TABLE]

Because they form an orthonormal basis, the three eigenvectors satisfy the identity $\hat{\bm{v}}_{1}\otimes\hat{\bm{v}}_{1}+\hat{\bm{v}}_{2}\otimes\hat{\bm{v}}_{2}+\hat{\bm{v}}_{3}\otimes\hat{\bm{v}}_{3}=\bm{I}$, hence we have

[TABLE]

The integral over the unit sphere is carried out using the relation

[TABLE]

where the scalar function $A(x,\ell)$ is

[TABLE]

The scalar function $A(x,\ell)$ can be regarded as a regularized distance function in the sense that $A(x,\ell)$ tends to $x$ when $x/\ell\gg 1$, while it smoothly approaches to $2\ell$ for small $x$, as shown in Fig. 2. By sake of (71), the Green tensor finally becomes:

[TABLE]

This result can also be obtained by direct inverse Fourier transform of (26), as shown in Appendix A. A more detailed analysis of the isotropic Green tensor (73) can be found in Lazar and Po (2018).

5 A comparison with Molecular Statics: The Kelvin problem

In this section, we compare the Green tensor obtained from Mindlin’s strain gradient elastic theory to that obtained from an atomistic system. This study was carried out using Minimol Tadmor and Miller (2011) which is a KIM-compliant molecular dynamics (MD) and molecular statics (MS) program. The Open Knowledgebase of Interatomic Models (KIM) is a project focused on creating standards for atomistic simulations including an application programming interface (API) for information exchange between atomistic simulation codes and interatomic potentials Tadmor et al. (2011, 2013).

We choose face-centered-cubic Aluminum and Copper for this comparison, and consider the following two interatomic potentials: the modified-embedded-atom-method (MEAM) by Lee (2001), and the embedded-atom-potential by Mendelev et al. (2008), which are archived in the OpenKIM repository. Elastic and gradient-elastic constants for these potentials were computed using the method described in Admal et al (2016), and they are available on the KIM repository Lee (2001); Mendelev (2008). For convenience, the values of the independent elastic and gradient-elastic constants are reported in table 1. These components are used to populate the elastic tensors $\mathbb{C}$ and $\mathbb{D}$ Admal et al (2016); Auffray et al. (2013). The Voigt structure of the resulting tensors $\mathbb{C}$ and $\mathbb{D}$ is shown in Fig. 3.

The atomistic system is constructed by stacking together $15\times 15\times 15$ unit cells resulting in $13500$ atoms. A force of $0.0116$ eV/Å in the $x_{1}$ direction is imposed on the central atom of the system, and displacement boundary conditions are imposed on five layers of atoms close to the boundary using the classical solution given in Eq. (49). The padding atoms thickness is 0.15 times the size of the box. A MS simulation is carried out using the above-mentioned boundary conditions resulting in a deformed crystal. The resulting displacement field normalized with respect to the force on the central atom yields the atomistic Green tensor component fields.

Simulation results are shown in Fig. 4, where we compare the Green tensor components $G_{11}(x_{1},0,0)$ and $G_{22}(x_{1},0,0)$. Despite the fact that these potentials were never fitted to gradient-elastic constants, it can be observed that the analytical predictions are in good agreement with MS calculations, with a maximum error at the origin in the order of 5-30%, depending on the potential used. It should be noted that, compared to the EAM potential, the MEAM potential better compares to the analytical results, possibly as a result of artifacts in gradient-elastic constants evaluated by EAM potentials Admal et al (2016).

6 Summary and Conclusions

In this paper we have derived an expression for the Green tensor of Mindlin’s anisotropic strain gradient elasticity, which possesses up to 21 elastic constants and 171 gradient elastic constants in the general case of triclinic media. The Green tensor is found in terms of a matrix kernel integrated over the unit sphere in Fourier space. Such representation is similar to that of the classical anisotropic Green tensor, which requires integration over the equatorial plane of the unit sphere. In contrast to its classical counterpart, however, the Green tensor of Mindlin’s anisotropic strain gradient elasticity is non-singular at the origin, while its gradient is finite but discontinuous at the origin. It is shown that the non-singular Green tensor converges to the classical tensor a few characteristic lengths away from the origin. Therefore, the Green tensor of Mindlin’s first strain gradient elasticity can be regarded as a physical regularization of the classical anisotropic Green tensor. Moreover, existing expressions of the Green tensor found in the literature are recovered as special cases. Because the Green tensor regularizes its classical counterpart without unphysical singularities, it offers a more realistic description of near-core elastic fields of defects in micro-mechanics, and it provides more accurate boundary conditions for atomistic and ab-initio energy-minimization calculations. As an illustrative example, we have computed the displacement field induced by a concentrated force acting at the origin (Kelvin problem), and compared the analytical predictions to atomistic calculations when the elastic and gradient-elastic moduli are consistently derived from the interatomic potentials. Despite the fact that these potentials were not fitted to gradient-elastic constants, it is shown that the analytical predictions are in good agreement with MS calculations, with a maximum error at the origin in the order of 5-30%, depending on the potential used.

List of Abbreviations

PDE: partial differential equation. SPD: symmetric positive definite. KIM: Open Knowledgebase of Interatomic Models. API: application programming interface. EAM: embedded atom method. MEAM: modified embedded atom method.

Declarations

Availability of data and materials. Elastic and gradient-elastic material constants used to obtain the results in section 5 are freely available as part of the Open Knowledgebase of Interatomic Models (KIM).

Competing Interest. The authors declare that they have no competing interests.

Funding. G.P. acknowledges the support of the U.S. Department of Energy, Office of Fusion Energy, through the DOE award number DE-SC0018410, the Air Force Office of Scientific Research (AFOSR), through award number FA9550-16-1-0444, and the National Science Foundation, Division of Civil, Mechanical and Manufacturing Innovation (CMMI), through award number 1563427 with UCLA. N.A. acknowledges the support of the US Department of Energy’s Office of Fusion Energy Sciences, Grant No. DE-SC0012774:0001. M.L. gratefully acknowledges a grant from the Deutsche Forschungsgemeinschaft (Grant No. La1974/4-1).

Authors Contribution. G.P. and M.L. obtained the expression of the Green Tensor. N.A. and G.P. carried out the numerical analysis. All authors read and approved the final manuscript.

Appendix A Direct derivation of Mindlin’s isotropic strain gradient elasticity of form II

Plugging (57) and (58) into (26) we have

[TABLE]

Owing to its special structure (see footnote 6), the matrix in the numerator can be easily inverted. In index notation the result is

[TABLE]

Using the the general form of the Fourier transform of the derivative, the Green tensor in real space is obtained as

[TABLE]

Now consider the identity

[TABLE]

Using (77) in (76), the Green tensor (73) is readily recovered.

List of Figures

1. 1 The unit sphere in Fourier space. The unit vector $\bm{\kappa}(\theta,\phi)$ is defined by the azimuth angle $\phi$, and the zenith angle $\theta$ measured from the axis $\hat{\bm{e}}_{3}=\bm{x}/x$.
2. 2 Plot of the regularized distance function $A(x,\ell)$.
3. 3 Voigt representation of the elastic tensors $\mathbb{C}$ and gradient-elastic tensor $\mathbb{D}$ for fcc Al and Cu, computed from the interatomic potentials Lee (2001) and Mendelev (2008). and Cu for EAM potential Mendelev (2008). and Cu for MEAM potential Lee (2001). and Al for MEAM potential Lee (2001).
4. [4 Components of the Green tensor for fcc Al and Cu, and comparison to atomistic calculations obtained from the interatomic potentials Lee (2001) and Mendelev (2008).

• Cu for EAM potential Mendelev (2008).
• Cu for MEAM potential Lee (2001).
• Al for MEAM potential Lee (2001).](#S5.F4)