2.1

Eringen’s Stress Gradient Theory

Eringen’s Stress Gradient Theory (ESGT) of elasticity assumes an equivalent effect due to nearest neighbor interaction and beyond the single lattice in the sense of lattice average stress and strain. In this theory, it is assumed that the stress state at any reference point given by coordinates x = (x₁,x₂,x₃) in the body is assumed to be dependent not only on the strain state at this point, but also on the strain states at all other neighboring points given by coordinate vector x’ of the body. This is in accordance with atomic theory of lattice dynamics and experimental observations on phonon dispersion. The most general form of the constitutive relation in the nonlocal elasticity type representation involves an integral over the entire region of interest. The integral contains a nonlocal kernel function, which describes the relative influences of the strains at various locations on the stress at a given location. The constitutive equations of linear, homogeneous, isotropic, non-local elastic solid with zero body forces are given by ¹²

where

Eqn. (1) is the equation of equilibrium, where σ_kl,ρ,f^l and u_l are the stress tensor, mass density, body force density and the displacement vector at a reference point x in the body, respectively, at time t. Eqn. (2) is the classical constitutive relation where is the classical stress tensor at any point x′ in the body, which is related to the linear strain tensor at the same point through the Lame constants λ and μ. Eqn. (2) is the classical strain-displacement relationship. The only difference between the above three equations in Eqns. (2) and the corresponding equations of classical elasticity is the introduction of first of the Eqn.(2), which relates the global (or nonlocal) stress tensor σ_kl to the classical stress tensor using the modulus of nonlocal ness. The nonlocal modulus α(∣x – x^′∣, ζ) is the kernel of the integral in the first of Eqn. (2) and contains parameters which correspond to the nonlocal parameter ¹³. A dimensional analysis of first of the Eqn. (2) clearly shows that the nonlocal modulus has dimensions of (length)^-3 and hence it depends on a characteristic length ratio a/l, where a is an internal characteristic length (lattice parameter, size of grain, granular distance) and l is an external characteristic length of the system (wavelength, crack length, size or dimensions of sample)¹⁴, Ω is the region occupied by the body. Therefore the nonlocal modulus can be written in the following form:

where e₀ is a constant appropriate to the material and has to be determined for each material independently¹⁴.

After making certain assumptions¹⁴, the integro-partial differential equations of the stress gradient theory can be simplified to partial differential equations. For example, first of Eqn. (2) takes the following simple form:

where C_ijkl is the elastic modulus tensor of classical isotropic elasticity and ε_ij is the strain tensor. where ∇ denotes the second order spatial gradient applied on the stress tensor σ_kl,k and ξ = e₀1/a. The validity of Eqn. (4) is established by comparing the expressions for frequency of waves from the above ESGT model with those of the Born-Karman model of lattice dynamics¹⁴. Eringen reports a maximum difference of 6% and a perfect match for nonlocal constant value of e₀ = 0.39¹⁴. Eqn.(4) on simplification will yield the following constitutive equations for 1-D and 2-D problems is given by

where all the σ s and τ have usual meaning.

2.2

Strain Gradient Theory

As mentioned earlier there are many strain gradient theories reported in the literature. The most notable of them is the Mindlin theory¹⁵ formulated in 1964. This theory’s constitutive model gives a total of 1764 coefficients in the constitutive matrix out of which 903 are independent. Even for isotropic material assumption, the number of independent constants to be determined is equal to 18. It is indeed a daunting task to determine them and hence, this theory, although is rich in completeness and theoretical basis, has very little practical use. Hence, in order to increase its utility and reduce the number of constants to be determined, many researchers including Mindlin have made series of assumptions on deformation field of the microstructure and many such theories have been reported in the literature. For wave propagation analysis, we need a simple theory that not only has less number of constants in the constitutive model, but also is stable. That is, one of the problems associated with the strain gradient theories is that in some theories, the governing differential equations are such that they give unstable or divergent solutions (see¹⁶] for more details on stability issues associated with strain gradient theories). Hence, in this paper, we will outline a very simple strain gradient theory that can provide good insight into propagation of waves in nanostructures. One such theory is the Laplacian based Strain Gradient Theory.

Laplacian based strain gradient theories are extensively used in static analysis, especially in those structures involving cracks, primarily to overcome the effects of stress singularities near the crack tips. However its use in dynamic analysis is quite different, where it is primarily used to describe the dispersive wave propagation in a heterogeneous media. Laplacian based strain gradient theory can be derived using simple lattice model consisting of springs and masses shown in Fig 3.

The figure shows the deformation of the 1-D lattice at discrete points n+2, n+1, n,n-1,n-2, ........ Let us consider the 3 particles in this lattice at points n,n-1, and n+1. If we isolate these points, draw the free body diagram and apply Newton’s second law, we get

Next, we will convert this discritized equation motion into a continuum equation. For this, we will expand the deformation of the n + 1 and n – 1 particles having a spacing d between them in Taylor’s series as

We assume that homogenized lattice has a Young’s modulus E and density ρ and they can be expressed in terms of the lattice spacing constant K and lattice mass M as E = Kd/A and M = ρAd, where A is the area of the equivalent continuum and d is the spacing of particles in the lattice. Using Eqn. (8) in Eqn. (7) and ignoring the higher order terms in the Taylor series and assuming u_n = u(x, t), we get

The equation above can be re-written as

From the equation above, one can easily recognize that the left hand represents the constitutive model of the strain gradient theory and can be written as

The equation above represents the constitutive model for a 1-D case and this equation can be generalized to the 3-D case as

where g is the length scale parameter, which is normally expressed in terms of lattice parameter d and C_ijkl is the fourth order tensor. The constitutive model given in Eqn. (12) has been proposed by a number of researchers^17,18. The main problem with the above constitutive models is that the final solutions provided by certain order of this constitutive model are neither unique nor stable. For example, the constructive model for the second order strain gradient mode of 1-D nanorod, which is given by

is shown to be unstable (see¹⁰), while the fourth order strain gradient model of the same nanorod, whose constitutive model is given by

is found to be highly stable. However, Aifantis and his coworkers^19,20 derived the second order strain gradient elasticity constitutive model as

Note that the main difference between Eqn. (13) and Eqn. (15) is that the negative sign before the higher order strain terms. The above models said to be highly stable. The sign of the gradient term and the issues of uniqueness and stability as opposed to their ability to describe dispersive wave propagation have opened up serious debates across the elasticity community; for example, see for instance the early study of Mindlin and Tiersten²¹ and Yang and Gao²² where this dilemma was aptly named as sign paradox

3.1

Eringen’s Stress Gradient Nano Rod:

The first step in performing wave propagation analysis is to obtain the governing differential equation. Rods can sustain only the longitudinal motion along the axis of the nanorod system and hence the deformation field can be written as u(x,y,z,t)=(x,t). Using the constitutive model given by Eqn. (5) and theory of elasticity principle, one can write the governing differential equation for a stress gradient nanorod as

The first step in determining the wavenumber and group speeds is to first transform Eqn. (16) to frequency domain using Discrete Fourier Transform (DFT) using the transformation given by

where N is the number of FFT points with as circular frequency. Here . Substituting Eqn. (17) in Eqn. (16), the original PDE will be reduced to a set of ODEs, which are of constant coefficient type, whose solution will be of the form û(x,ω)= Ae^ikx, where k is the wavenumber, which requires to be computed. Substituting the assumed constant coefficient solution in the reduced set of ODEs, we get the characteristic equation for computation of wavenumber, which is given by

where Solution of Eqn. (18) will yield

The wave frequency is a function of wavenumber k, the nonlocal scaling parameter e₀a and the material properties (E & ρ) of the nanorod. If e₀a = 0, the wavenumber is directly proportional to wave frequency, which will give a nondispersive wave behavior. The cut-off frequency of this nanorod is obtained by setting k = 0 in the dispersion relation (Eqn. (18)). For the present case, the cut-off frequency is zero, that is, the axial wave starts propagating from zero frequency. The wavenumber as a function of nonlocal scale parameter e₀a is shown in Fig.4. From this figure, one can observe that at certain frequencies, the wavenumber is tending to infinity and value of this frequency, which is called the escape frequency, decreases with increase in the scale parameter. Its value can be analytically determined by looking at the wavenumber

Figure 4.

Dispersion relations for Eringen’s Nano Rods (s) Wavenumber relations (b) Group Speeds

The escape frequency is inversely proportional to the nonlocal scaling param-eter and is independent of the diameter of the nanorod. The variation of the escape frequency with nonlocal scaling parameter is shown in Fig. 4(a). Next, we will compute the wave speeds, namely the phase speeds and the group speeds, whose expressions are given by

We see from Eqn. (21) that both phase and group speeds depend on the non-local parameter and the speeds depend on the frequency indicating the dispersiveness of the waves. When e₀a →0, we obtain the local rod solution, that is . For the present analysis, a SWCNT is assumed as a nanorod. The values of the radius, thickness, Young’s modulus and density are assumed as 3.5 nm, 0.35 nm, 1.03 TPa, and 2300 kg/m³, respectively. Fig. 4(a) shows the variation of the axial wavenumber of a nanorod. The thick lines represent the real part and the thin lines show the imaginary part of the wavenumbers. From the figure, for a nanorod, it can be seen that there is only one mode of wave propagation that is, the axial or the longitudinal. For local or classical model, the wavenumbers for the axial mode has a linear variation with the frequency, which is in the THz range. The linear variation of the wavenumbers denotes that the waves will propagate non-dispersively. On the other hand, the wavenumbers obtained from nonlocal elasticity have a non-linear variation with the frequency, which indicates that the waves are dispersive in nature. However, the wavenumbers of this wave mode have a substantial real part starting from the zero frequency. This implies that the mode starts propagating at any excitation frequency and does not have a cut-off frequency. At the escape frequency, which was defined earlier, the wavenumbers tends to infinity as shown in Fig. 4(a). Hence, the nonlocal elasticity model shows that the wave will propagate only up to certain frequencies and after that the wave will not propagate.

From Eqn. (21), we see that escape frequency are purely the function of non-local scaling parameter. From Eqn. (20), we see that as the value of e₀a increases, the numerical value of escape frequency decreases. Fig. 4(b) shows the variation of phase speeds with frequency. We can clearly see the shifting of the dispersion curves with the increase in the non-local parameter. In addition, one can see the extent of non-dispersiveness of waves in Eringen’s nano rods.

3.2

Strain Gradient Nano Rod:

Using the constitutive equations given by Eqn. (13) and (14), we can write the governing equations for the second and fourth order strain gradient rods as

From the equations above, the spatial order of the governing equations has increased from two in local model to four in second order strain gradient model and six in the fourth order model. Let us first consider the second order strain gradient model. Following the procedure outlined in the previous section, we can write the characteristic equation for wavenumber computation as

Fourth order equation in spatial direction will yield four wave modes (2 incident waves and the rest reflected wave components). The corresponding phase speed and group speeds can be written as

As before in the case of Eringen’s rod case, one can see that we can recover local rod solution by substituting e₀ = 0.We will next compute one of the important feature of wave propagation of Second order SGT nanorods, which is its Critical Wavenumber, which is the wavenumber at zero frequency, which can be obtained from Eqn. (24) and its expression can be written as

which is inversely proportional to the non-local parameter. As mentioned earlier, this rod will yield no-unique unstable solution (see ¹⁰ for details). Next we will derive the dispersion relation for the fourth order strain gradient model, which is given by

Corresponding phase and group speed expressions can be written be as

Next, we will next plot the wavenumber and speeds with respect to the frequency. Such plots are generated using the diameter (d), Young’s modulus (E) and density (ρ) are assumed as 5 nm, 1.06 TPa, and 2270 kg/m³, respectively. The wavenumber relations are generated for e₀a = 0.2 and it is shown in Fig 5.

Figure 5.

Wavenumber relations for Strain gradient Nano Rods

The waves in both second and fourth order strain gradient are non-dispersive unlike local rod. It can also be seen that the fourth-order strain gradient model give improved approximation over the second order strain gradient model, as compared to the classical continuum model. They are also compared with the results are of the Born-Karman model²³ (atomistic model) and the Eringen’s nano rod model presented previously. The instability of the Second order SGT model can be seen developing in Fig. 5 for wavenumbers larger than , where the angular frequency and the phase velocity become imaginary. This means that waves with larger wavenumbers (or, equivalently, with smaller wavelengths) cannot propagate through this medium. Instead, the imaginary frequency and velocity imply that the response occurs everywhere in the medium instantaneously. This is physically unrealistic. Therefore, these smaller wavelengths should not be considered. Filtering shorter waves occurs automatically in a discrete medium, where wavelengths smaller than two times the particle size cannot be monitored. However, in a continuous medium, all wavelengths can in principle be present. Especially when shock waves are investigated, waves with all wavelengths are triggered by the loading. The imaginary angular frequency (or phase velocity) of these high-frequency waves prohibits a proper wave propagation simulation with this model. The cut-off value for the wavenumber occurs at (see Fig. 5). In the response of the Fourth order SGT model, the effect of these high frequency waves are of minor importance. In summary, the wave behavior in Second order FGT and Fourth order SGT models are quite different and in addition, instability in the second order model can be clearly seen.

3.3

Stress Gradient Nano Beam:

In this section we will derive the governing differential equation for a stress gradient nano rods. The beam deformation field is based on Euler-Bernoulii beam theory and is given by

where w and u⁰ are the transverse and axial displacements of the point (x,0) on the middle plane (that is, z = 0) of the beam, where z is the thickness coordinate. The only nonzero strain of the Euler-Bernoulli beam theory, accounting for the strain is the axial strain given by

Using the constitutive model given in Eqn. (5) and the theory of elasticity procedure, we get the following governing equation do stress gradient rod, which is given by

we can clearly see that when e₀a=0, we recover the governing equation of the local Euler Bernoulli beam.

Next, we will transform Eqn. (31) to frequency domain using DFT and this procedure will transform the Governing PDE to a set of ODEs, which will again be constant coefficient type. Adopting the same procedure as was done for nano rods, we get the following characteristic equation for wavenumber computation, which is given by

Fourth order equation will give four modes, two of which are forward moving wave modes, while the other two correspond to backward moving or reflected wave components. Solving Eqn. (32), we get

From the above equation, we can see that unlike Eringen’s nano rod, the Eringen’s nano beam does not exhibit any escape frequency and in addition, cut-off frequency also does not exist for this kind of beam. The phase speed can be calculated from C_p =ω/k. The expression for group speed is given by

where k_α =1,2,3,4 corresponding to four wave modes. From this expression we can plot the dispersion relation. This plot will give full description of the wave propagation in beams. Both the speeds are also a function of nonlocal scaling parameter and wave circular frequency.

The spectrum curves (wavenumber vs. frequency) and dispersion curves (group speed vs. frequency) are shown for both local and nonlocal continuum models in Fig. 6 (a) and Fig. 6 (b).

Figure 6.

Wavenumber and dispersion relations for Stress gradient Nano beams

As mentioned earlier, the plot does not show the existence of escape or cut-off frequency and the behavior is pretty much the same as local beam solution. The main difference is that, the scale parameter changes the slope of the wavenumber plot, which means substantial change in the group speeds as shown in Fig. 6(b). That is, introduction of scale parameter, reduces the group speeds in comparison with local Euler-Bernoulli beam.

4.1

Governing Equations for Flexural Wave Propagation in Mono- layer Graphene Sheets

Graphene sheet is idealized as 2-D plate that is subjected to both in-plane and out-of-plane loading. As a result, the Graphene sheet can undergo both axial, bending and shear deformation. Fig. 7(a) shows a rectangular Graphene sheet and Fig. 7(b) shows its equivalent continuum model.

Figure 7.

Single-layered Graphene sheet: (a) Discrete model (a monolayer Graphene of 40 ×40, consists of 680 carbon atoms arranged in hexagonal array) (b) Equivalent continuum model.

Here, the Graphene is treated as an isotropic plate. If w(x,y,t) is the transverse displacement, the governing differential equation for such a structure considering Eringen’s non-local model constitutive relations (Eqn.(6)) is given by (see ¹⁰ for more details)

where

with E,p,v and z being the Young’s Modulus, density, Poisson’s ratio and thickness coordinate, respectively.

4.2

Wave Dispersion Analysis:

The deformation field is transformed to frequency domain and it is assumed that the model is unbounded in Y –direction (although bounded in X–direction). Thus the assumed form is a combination of Discrete Fourier transforms in Y–direction and Fourier transform in time, which is written as

The ω_n and the η_m are the circular frequency at n^th sampling point and the wavenumber in Y–direction at the m^th sampling point, respectively. The N is the index corresponding to the Nyquist frequency in FFT, and . Substituting Eqn.(36) in Eqn.(35), we get the following ODE

where

Since this is an ODE is having constant coefficients, its solution can be written as , where k is the wavenumber in the X–direction, yet to be determined and is an unknown wave coefficient. Substituting this assumed form of ŵ(x) in the ODE gives for ), we get the characteristic equation for computation of wavenumber as

Solving the above equation, we get

Next, we compute the group speed C_g=dω/dk, which is given by

where

From the above expressions, we see that both the group speeds and wavenumber are functions of Y-directional wavenumber and the non-local scale parameter.

Looking at Eqn. (37), the term H₀ indicates the possibility of a waveguide having cut-off frequencies. This is obtained by setting k = 0 in the dispersion relation (Eqn. (40)) that is, for the present case one can set H₀= 0, which gives the cut-off frequency expression as

The cut-off frequency is directly proportional to the Y–directional wavenumber (η_m) and is also dependent on the nonlocal scaling parameter. For η_m = 0, the wavenumbers of the flexural wave mode have a substantial real part starting from the zero frequency, which implies that the mode starts propagating at any excitation frequency and does not have a cut-off frequency. For η_m = 0, the flexural wave mode, however, has a certain frequency band within which the corresponding wavenumbers are purely imaginary. Thus, the wave mode does not propagate at frequencies lying within this band. These wavenumbers have a substantial imaginary part along with the real part, thus these waves attenuate as they propagate.

Next, the term H₄ in Eqn. (37) indicates the possibility of a waveguide having escape frequencies. Its value can be analytically determined by looking at the wavenumber expression and setting k → ∞. This accounts to setting the H₄ = 0, which gives

Next, we will investigate the relations derived above and plot the wave behavior in the monolayer Graphene sheet. For the present wave propagation analysis, the material properties of the Graphene are assumed as follows: Young’s modulus E = 1.06 TPa and density ρ = 2300 kg/m³. The choice of effective wall thickness t of nanostructures such as CNT, Graphene, etc., is a long-standing issue in nano mechanics. One of the best approaches to estimate the thickness of CNT (that is, rolled Graphene sheet) is to model single-wall CNTs as linear elastic thin shells ²⁷. The shell thickness t is determined by fitting the atomistic simulation results of tensile rigidity and bending rigidity of single-wall CNTs. Such an approach gives the CNT thickness t much smaller than the graphite inter-layer spacing 0.34 nm, ranging from 0.06 to 0.09 nm. The scattered CNT thickness 0.06 – 0.09 nm depends on the inter- atomic potential as well as simulation details. The thickness of the Graphene chosen here is equal to t = 0.089 nm, obtained by Kudin et. al.²⁷ via ab inito computations. They defined the effective thickness of Graphene or CNTs as .

The flexural wavenumber dispersion with wave frequency in the Graphene s shown in Figs. 8 (a) and 8 (b), respectively, obtained from both local and ESGT theories. For the present analysis, the nonlocal scaling parameter is assumed as e₀a = 0.5 nm.

Figure 8.

Wavenumber dispersion in monolayer graphene sheets (a) local elasticity (e₀a =0 nm) (b) nonlocal elasticity (e₀a =0.5 nm). Wavenumber variation at lower frequencies is shown separately for clear visibility

The spectrum curves shown in Fig. 8(a) is for η_m = 0 (1D wave propagation), 3, 5, 10 nm^-1. The local elasticity calculation shows that the flexural wavenumber follow a nonlinear variation at low values of wave frequency; and at higher frequencies, it varies linearly as shown in Fig. 8(a). This nonlinear variation indicates that the waves are dispersive in nature, that is, the waves will change their shape as they propagate. The linear variation indicates that the waves are in non-dispersive nature. For η_m = 0, the wavenumbers of the flexural wave mode does not have a cut-off frequency. As η_m increases, all the waves are still dispersive in nature as shown in Fig. 8(a). As the Y–directional wavenumber increases from 0 to 10 nm^-1, the wave modes are having a frequency band gap region, within which the corresponding wavenumbers are purely imaginary. Thus, the flexural mode does not propagate at frequencies lying within this band. Hence, these wavenumbers have a substantial imaginary part along with the real part, thus these waves attenuate as they propagate. It can also be seen that from Fig. 8(a), that the frequency band also increases with increase in η_m.

The wavenumber dispersion with frequency obtained from nonlocal elasticity (e₀a = 0.5 nm) is shown in Fig. 8(b). The observations made in local elasticity are still valid in nonlocal elasticity also. The only difference is that, because of nonlocal elasticity, the wavenumbers of the flexural wave become highly non-linear and tends to infinity at escape frequency. It can be seen that the wavenumbers before escape frequency are real and after that imaginary. The cut-off frequency of the flexural wave for η_m = 3, 5 and 10 nm^-1, respectively, occurs at 0.8087 THz, 2.2280 THz, and 8.6820 THz in local/classical elasticity, and at 0.4578 THz, 0.8392 THz, and 1.7090 THz in nonlocal elasticity (e₀a = 0.5 nm). The nonlocal scale highly affects the frequency band gap of the flexural waves in Graphene sheet. The escape frequency of this flexural wave is 6.9580 THz for e₀a = 0.5 nm. It has also been observed that the escape frequencies are independent of η_m from Fig. 8(b).