A semianalytical approach for determining the nonclassical mechanical properties of materials

1 Introduction

Phonons, which describe vibrational motions in condensed matter, are collective excitations in the periodic elastic arrangement of atoms [1]. Phonon propagation in materials, either single crystals or polycrystalline materials, has been used widely for characterizing materials because it illustrates the vibrational modes as well as the atomic interactions in materials. Studying phonon dispersions has become more and more important and favored for characterizing materials as scientific instruments and measurement techniques, such as dispersive-phonon imaging [2], picosecond laser ultrasonic [3], infrared spectroscopy [4], neutron spectroscopy [5], Raman spectroscopy [6], and Brillouin spectroscopy [7], were developed leading to a significant improvement in accurately measuring the phonon modes propagating in solids. These achievements resulted in revealing unbounded information regarding materials’ properties, such as mechanical [8], [9], thermal and thermoelectric [10], [11], electrical and flexo-electrical [12], [13], and chemical [14], suitable for engineering the materials’ functionalities. Therefore, it is crucial to employ phonon dispersions to understand materials’ properties to optimize their performance or tune their properties.

One of the most interesting applications of phonon dispersions is to understand the mechanical behaviors of materials on nanoscales [15]. Indeed, it was shown that to precisely demonstrate the phonon dispersions at the center of the first Brillouin zone based on the solid mechanics theories, very complicated solid mechanics formulas including several material constants, also called nonclassical properties, are required [9], [15], [16]. These complicated solid mechanics formulas had not been explained based on the classical continuum mechanics definition until the nonclassical solid mechanics theories, such as couple-stress theory (CST) [17] and SGT [18], were developed. Depending on the dimensions of the structures, the nonclassical mechanical properties have significant influences on mechanical behaviors of structures, especially when dimensions are nanometers and submicrometers, or negligible influences on the mechanical behaviors of the structures when the dimensions are millimeters and larger [8]. Equivalently, the classical continuum mechanics, which first proposed for bulk dimensions, fails to illustrate the mechanical behaviors of nanostructures because it ignores the nonclassical mechanical properties (also called length scales or nonclassical properties) at nanoscales [19], [20], [21].

Moreover, it has been shown over the years that classical continuum mechanics fails to precisely demonstrate the effects of surface energy [22], discrete nature of the materials [23], presence of micro- and nanocracks [24], [25], internal strains, and quantum confinements in the nanostructures [26]. These effects can be only investigated and explained using the nonclassical theories [15]. Therefore, it is crucial to implement nonclassical solid mechanics theories including nonclassical mechanical properties, such as gradient elastic theory [27], [28], [29], to illustrate the mechanical performance of nanostructures. However, to be theoretically and experimentally consistent, the first step is to establish a semianalytical approach, to combine the experimental and theoretical observations.

Phonon dispersions can be calculated using well-established quantum mechanics approaches, such as density-functional perturbation theory [30] and density-functional theory [31], or lattice-dynamic approaches, such as deformable-ion model, polarizable bond charges model, and rigid valence shell model [32], [33], [34], [35]. To the knowledge of the authors, the authors of [9] considered dispersive elastic constants for the first time to illustrate the deviation between the experimental results, the ballistic phonon imaging experiment, and the classical continuum mechanics as well as the lattice-dynamics simulation. Later [8], used the second SGT without acceleration gradients in conjugation with lattice-dynamics simulations and molecular-dynamics simulations to connect the dispersive elastic constants to the length scales. In contrast to [8], [36] only focused on the experimental data to find the nonclassical mechanical properties with the same approach. It should be mentioned the difference between the results of [8] and [36] could be due to lack of a connection between the theoretical and experimental approaches because each approach has some extended approximations and errors. A more complicated analysis was proposed by [37], where they described the wave dispersion using a gradient-enriched shell theory and gradient-enriched Timoshenko theory. Their approach resulted in an imaginary group velocity besides unusual cut-off frequencies that are inconsistent with observations [16]. However [38], represented a revision of [37] model leading to an improvement in the simulation results. Furthermore, in conjugation with the SGT, there is always a contradiction regarding the sign of gradient terms along with concerns of uniqueness and stability versus the ability to describe dispersive acoustic waves’ propagation. This controversial discussion was professionally reviewed in [16] while it was not evaluated using other nonclassical theories, such as CST [17] or modified couple-stress theory (MCST) [39], including the acceleration gradients and being combined with both atomic simulations and experimental data to resolve this contradiction fully.

Consequently, here we predict the nonclassical mechanical constants of several materials based on a semianalytical approach. To do so, we rigorously derive the most general forms of the phonon dispersions in the framework of modified couples stress theory including the acceleration gradients (MCST-AG). Then, the phonon-dispersions curves of the materials were calculated using atomic simulations. To fit the atomic simulation results with the derived phonon dispersions, we use the experimental values for the classical constants and then we find the nonclassical constants through fitting the atomic simulation results. To compare the present results of this approach with the earlier works, which were determined based on SGT without considering acceleration gradients, we calculated the phase velocity and group velocity in addition to the natural frequency of a simply–simply supported beam (SSSB) with different static length scale to radius ratios.

2 Governing equations

According to the MCST of [19] the strain energy density is a function of both strain (conjugated with stress) and curvature (conjugated with couple stress). It then follows that the strain energy U of a deformed linear elastic material occupying region Ω is given by

where the stress tensor, σ_ij, strain tensor, ε_ij, deviatoric part of the couple-stress tensor, m_ij, and symmetric curvature tensor, χ_ij, are, respectively, defined by

The C_ijkl and β_ijkl are the material constants expressing the classical constants and nonclassical constants, respectively. Furthermore, u is the displacement vector; and θ is the rotation vector. To consider the effects of the anisotropic properties of the materials, the materials constants are considered in the general form. However, to be consistent with the original form of the MCST and for simplicity, only one static length scale is considered, such that βijkl=8ls2Cijkl, where l_s is the static length scale. The total kinetic energy is given by

This form of the kinetic energy, which is a simplified version of [18] definition, has attracted special attention especially in the dynamical analysis of nanostructures. That is why this definition of the kinetic energy has an extra length scale, named the dynamical length scale, so it eliminates the dynamical instabilities [16]. It should be mentioned, because the focus is on studying the phonon-dispersion relations, the external loadings are ignored. Therefore, the total energy stored in an elastic body is

To derive the governing equations, the minimum energy principle is applied. Therefore,

Applying the variational theorem to Eq. (1) and using the symmetry property of the C_ijkl and β_ijkl, one obtains

Using the definitions of the strain tensor, symmetric curvature tensor, and microrotations vector provided in Eqs. (4)–(6), their variational forms are, respectively, determined by

Substituting Eqs. (11)–(13) into Eq. (10), it can be rearranged as

Applying integration-by-parts on the first and third terms in Eq. (14) gives

The last step to simplify Eq. (15) is to apply the divergence theorem to the second bracket. The divergence theorem converts the integral over the volume to an integral over the surface as follows

To find the variational form of the kinetic energy, the same process can be repeated for Eq. (7). Therefore,

Applying the integration-by-parts technique with respect to time, Eq. (17) can be rearranged as

Again, the integration-by-parts technique must be applied to the second term of the first integral of Eq. (18) with respect to the spatial derivative. Thus,

The last step to derive the variational form of the kinetic energy is to apply the divergence theorem to the last terms under the integrals of Eq. (19). Eventually, one gets

Inserting Eq. (16) and Eq. (20) into Eq. (9), the equations of motion will be achieved as follows,

Equations (21) are the most general form of the equations of motion in the scope of the MCST including the acceleration gradients. It should be emphasized that, because here the focus is to derive the mechanical constants using an atomic simulation, the equilibrium equations are free from external traction terms. However, the external tractions can be readily added into Eqs. (21) based on the domain that they are exerting on. Furthermore, for an isotropic linear elastic material, the material constants will be replaced with

Substituting Eqs. (22) into Eqs. (21), the following equations of motion for an isotropic linear elastic material can be obtained

Ignoring the dynamic length scale in Eqs. (23), the MCST’s governing equations of motions can be achieved as proposed by several authors [18], [39], [40].

3 Dispersion relations

In this section, the phonon-dispersion relations are derived for both longitudinal and transverse phonon modes. To do so, a longitudinal plane wave (LPW) (u₁≠0, u₂=u₃=0) and transverse plane waves (TPWs) (u₁=u₃=0, u₂≠0 and u₁=u₂=0, u₃≠0) are considered and implemented into Eq. (21a). Consequently, the dispersion relations for the longitudinal and transverse modes can be determined as

It should be mentioned that Eqs. (24) carry much interesting information. First, the static length scale does not contribute to the longitudinal dispersion relation, Eq. (24a). This fact can be explained from both solid mechanics and atomic structure point of view. From the solid mechanics view, that is why the curvature tensor is deviatoric, its diagonal terms associated with the longitudinal terms are zero. From the atomic structure view, the longitudinal phonons have higher energy than the transverse phonons, therefore, they should be less dispersed. However, this fact does not mean the longitudinal terms are immune to dispersion because the dynamic length scale appears in the longitudinal phonon-dispersion relation, which highlights the critical role of the dynamic length scale in demonstrating the phonon-dispersion relation. This conclusion is compatible with the prediction of the all previous works using the SGT that they calculated the static length scale for the longitudinal modes to be much smaller than the static length scale for the transverse modes. Second, Eqs. (24b) and (24c) explain the effects of the unit cell asymmetry on the transverse dispersion relations with fewer materials constants comparing the SGT. For example, for materials with the orthorhombic crystal structure, such as agrinierite and krennerite, because all the lattice constants are different, the transverse phonon-dispersion relations in different directions are different. Lastly, all dispersion relations, Eqs. (24), are dynamically stable and they suggest the phonon frequencies to be real and finite regardless of the l_ik (i=s, d) values.

Equations (25) represent the acoustic waves’ propagation in k-(110) and k-(111) directions, respectively. In Eqs. (25), for simplicity, the j-index is summation from 1 to 3. It should be mentioned that the remaining acoustic modes can be readily determined using Eqs. (25).

4 Results and discussion

4.1 Mechanical properties of ZnO, Si, SiC, InSb, and diamond

In this section, mechanical properties of several materials, ZnO, Si, SiC, InSb, and diamond, are studied using a semianalytical approach. First, the phonon-dispersion relations of the materials were determined using an atomic simulation. Then, the results of the simulations were analyzed by the definitions of the phonon dispersions derived in Eqs. (24). Each phonon-dispersion relation consists of at least two material constants, one classical and one nonclassical. According to Eqs. (24), the longitudinal mode includes a classical constant and the dynamic length scale while the transverse mode includes both static and dynamic length scales in addition to the classical constants. Therefore, using the available experimental results for the classical constants, the longitudinal phonon-dispersion mode first was fitted with the atomic simulation results to find the dynamic length scale. Then, using the classical constants of the experiments in conjunction with this dynamic length scale, the transverse phonon-dispersion modes were fitted with the atomic simulation results to achieve the static length scale as well. Because the procedure is the same for all materials, here we explain the procedure only for silicon (Si). Figure 1 shows the phonon dispersion of Si.

Figure 1:

Phonon dispersion of Si calculated using the Tersoff potential provided by [41].

Table 1 shows the estimated values for the length scales and a comparison with the literature. The experimental data used for the classical constants of Si are provided from [42]. It is noteworthy that, in addition to the completely different approaches used in the present work and the earlier works, here the phonon-dispersion relations were fitted for k-wavevectors up to 0.1 while they were fitted for much smaller values of k-wavevectors, for example, |k|<0.05 in [8].

Table 1:

Estimated length scales for Si.

	Transverse				Longitudinal
	This work				This work
Static length scale (nm)	0.21	0.225a	0.07b	0.12c	0	0.13a	0.234b	0.12c
Dynamic length scale (nm)	0.11	–	–	–	0.11	–	–	–

1. ^aReference [8].

2. ^bReference [36].

3. ^cReference [3].

Table 2 provides the predicted values for the dynamic and static length scales of Si, SiC, ZnO, InSb, and diamond. The experimental data used for the classical constants of SiC, ZnO, InSb, and diamond are provided from [43], [42], [44], and [45], respectively.

Table 2:

The static and dynamic length scales of Si, SiC, ZnO, InSb, and diamond.

	Materials
	Si	SiC	ZnO	InSb	InSba	Diamond	Diamonda
Static length scale (nm)	0.21	0.30	0.14	0.34	0.241c, 0.372c	0.60	0.0251b, 0.1922b
Dynamic length scale (nm)	0.11	0.29	0.34	0.21	–	0.38	–

1. ^aSuperscripts 1 and 2 represent the values for longitudinal and transverse directions, respectively.

2. ^bReference [8].

3. ^cReference [36].

4.2 Phase velocity and group velocity

In this section, the phase and group velocities of Si were studied as a criterion for comparing the results of this work and the previous results. That is why the wave propagation analysis verifies the stability of the solid mechanics theories [46], [47]. The phase of an acoustic wave in a nanostructure is the rate at which the phase of the wave propagates. Therefore, the phase velocity is defined as

(26)vpi=ωik, i=L, T,

where ω_i is the angular frequency, L and T stand for longitudinal and transverse modes, respectively. Implementing Eqs. (24) and the estimated length scales in Table 1 into Eq. (26), the normalized longitudinal and transverse phase velocities of the Si were calculated and shown in Figure 2. Furthermore, the normalized longitudinal and transverse phase velocities of the Si were calculated using the strain-gradient theory.

Figure 2:

The normalized longitudinal phase velocity (A), and the normalized transverse phase velocity (B), of Si. a, b, and c curves were plotted according to the SGT and predicted mechanical constants in [8], [36], and [3], respectively.

As can be seen from Figure 2, the previous works based on the SGT, which do not consider acceleration gradients, show unstable behaviors [16]. That is why these models predict the square of the longitudinal and transverse phase velocities to grow unboundedly toward negative infinity as the wavevector increases, resulting in imaginary phase velocities. This unrealistic behavior is due to the sign paradox of the strain gradient terms as well as neglecting the effects of acceleration gradients and considering very small wavevectors for predicting the mechanical constants. In the efforts of eliminating this inconsistency [48], have reformulated the governing equation of SGT to justify its stabilities at high wave numbers. General speaking, the acceleration gradients introduce a propagation characteristic into the governing equations determining the propagations of waves with higher wave numbers [47], which play a critical role in predicting the dispersion relations. Considering very small wavevectors for fitting the atomic simulations or experimental data caused those works to correctly predict the longitudinal and transverse phase velocities for extremely small wavevectors, but failed for larger wavevectors [16].

In addition to the longitudinal and transverse phase velocities, we calculated the longitudinal and transverse group velocities of Si by the same process explained for the longitudinal and transverse phase velocities. According to the definition of the group velocity, which is the velocity with which the overall shape of the wave amplitude propagates, the longitudinal and transverse group velocities can be calculated as

(27)vgi=∂ωi∂k, i=L, T.

Figure 3 shows the normalized longitudinal and transverse group velocities inside the Si. The explanation for the group velocities’ behavior is the same as for the phase velocities.

Figure 3:

The normalized longitudinal group velocity (A), and the normalized transverse group velocity, (B), of Si. a, b, and c curves were plotted according to the SGT and predicted mechanical constants in [8], [36], and [3], respectively.

4.3 Natural frequency

Lastly, in this section to evaluate the MCST-AG relative to SGT for demonstrating the dynamical behavior of nanostructures, we calculate the natural frequencies of an simply-simply supported beam (SSSB). The general forms of the governing equation of motions and the boundary conditions as well as initial conditions of a beam can be determined using Hamilton’s principle as mentioned in the first section, Eq. (9). For a one-dimensional beam, the displacement components can be considered as follows

(28)ux=−zdw(x)dx, uy=0, uz=w(x),

where u_x, u_y, and u_z are the displacement components in the x, y, and z directions, respectively. Here, it is assumed that the beam is along the x-axis in the x–z plane. Substituting the displacement components into Eq. (9), for free vibration of a beam, one gets

(29)∂(U−T) =∂∫t∫0L[12(C1111I+4C1212ls2A)(d2wdx2)2− 12(ρA(dwdt)2+ρld2A(ddx(dwdt))2)]dxdt=∫t∫0L[(C1111I+4C1212ls2A)d4wdx4+ρAd2wdt2− ρld2Ad2dx2(d2wdt2)]∂wdxdt+ ∫t[(C1111I+4C1212ls2A)d3wdx3∂w]0Ldt− ∫0L[(ρAdwdt−ρld2Ad2dx2(dwdt))∂w]tdx+ ∫t[(C1111I+4C1212ls2A)d2wdx2∂(dwdx)]0Ldt− ∫t[ρld2Addx(dwdt)∂(dwdt)]0Ldt=0.

Consequently, the governing equation of motion, initial conditions, and boundary conditions are given by

Governing equation of motion:

(30a)(C1111I+4C1212ls2A)d4wdx4+ρAd2wdt2−ρld2Ad2dx2(d2wdt2)=0.

Initial conditions:

(30b)(ρAdwdt−ρld2Ad2dx2(dwdt))∂w|t=0,ddx(dwdt)∂(dwdt)|t=0.

Boundary conditions:

(30c)w=0 or d3wdx3=0dwdx=0   or d2wdx2=0

Equations (30) provide the most general forms of the governing equation of motion, initial conditions, and boundary conditions of a beam under the free vibration condition. However, because here our goal is to compare the MCST-AG with the SGT, we only consider an SSSB, and then we calculate its first natural frequency. For free vibration of an SSSB, it can be shown that the solution is in the form of ω(x, t)=X(x)exp(iωt). Thus, the governing equation of motion will be

(31)(C1111I+4C1212ls2A)d4Xdx4+ρld2Aω2d2Xdx2−ρAω2X=0.

Because Eq. (31) is a fourth-order differential equation with constant coefficients, the general form of its solution is X(x)∝exp(λt). Therefore, its characteristic equation is given by

(32)Bλ4+Dλ2−G=0,B=C1111I+4C1212ls2A, D=ρld2Aω2, G=ρAω2.

The solutions of the characteristic equation are

(33)λ1, λ2=±jD+D2+4BG2Bλ3, λ4=±−D+D2+4BG2B

Therefore, the X(x) is

(34)X(x)=a1sin(D+D2+4BG2B)+a2cos(D+D2+4BG2B)  +a3sinh(−D+D2+4BG2B)+a4cosh(−D+D2+4BG2B).

For an SSSB, the four boundary conditions are

(35)w(0)=w(L)=0,d2wdx2(0)=d2wdx2(L)=0.

Implementing Eq. (34) into Eqs. (35) results in a system of four equations that can be written as follows

(36)[H(ω)]4×4{a1a2a3a4}=0.

To have a nonzero solution for a_i (i=1, 2, 3, 4) of Eq. (36), the following conditions must be satisfied

(37)det([H(ω)])=0.

Solving Eq. (37) provides all the natural frequencies of the SSSB. The natural frequencies of an SSSB based on the SGT can be determined using the same procedure. Here, we avoided deriving the relations for natural frequencies of an SSSB based on SGT, the interested readers are referred to [49] and [50] for more details.

Figure 4 shows the results of normalized values for the first natural frequency of an SSSB in terms of radius to the static length scale ratio. The normalized natural frequencies based on all methods approach each other as the radius to the static length scale ratio increases. For smaller ratios, the present model predicts lower values for the natural frequency because of the existence of the acceleration gradients that act as an extra term in the inertia matrix.

Figure 4:

Normalized first natural frequency of an SSSB based on the strain-gradient theory and the MCST-AG, l_i is the static length scale in each theory. The curves are normalized with respect to the first natural frequency according to the classical continuum mechanics. a, b, and c are from [8], [35], [3], respectively.

5 Conclusion

Here, we represented the equation of motion and dispersion relations based on the MCST-AG. Using MCST-AG suggests that the static length scale does not contribute to the longitudinal phonon dispersions; however, the dynamic length scale contributes to the phonon-dispersion relations, while both the static and dynamic length scales affect the transverse phonon dispersions. This result is compatible with the SGT, where it predicts the static length scale in the longitudinal direction to be very small compared with the static length scale in transverse directions. Combining these formulas with experimental results and atomic simulations, a semianalytical approach, the nonclassical length scales, static and dynamic length scales, were predicted for ZnO, Si, SiC, InSb, and diamond. To evaluate the accuracy of these formulas and technique, the phase velocity, group velocity, and natural frequency of Si SSSB were examined. This approach suggests the phase velocity and group velocity to be finite and real, which proves this approach provides dynamically stable behavior for nanostructures even for large values of wavevector. Furthermore, because of the appearance of the dynamical length scale, this approach predicts the natural frequencies to be somewhat lower than predicted natural frequencies with the strain-gradient approach.