Second gradient continuum model for anisotropic elastic and piezoelectric structures calibrated based on phonon dispersion relations

Abstract

In this paper, we present a variant of the second gradient continuum model that allows description of the spatial dispersion effects of high-frequency waves in anisotropic crystals. The considered model takes into account the strain and inertia gradient effects and the classical electro-mechanical coupling. The simplified constitutive equations of the model contain seven and four length scale parameters for piezoelectric hexagonal crystals and elastic cubic crystals, respectively. It is shown that all except one length scale parameter can be identified based on the fitting of the continuum model to the lattice dynamics calculations for the dispersion relations of the bulk acoustic phonons. Examples of the identification are presented for hexagonal piezoelectric ZnO and AlN and for cubic diamond and Si crystals. Using a calibrated model, we then obtain the dispersion relations for the high-frequency Bleustein–Gulyaev wave in piezoelectric crystals and for the shear horizontal surface wave in elastic crystals. The last one is predicted by the gradient models only for the high frequencies, while at low frequencies this wave degenerates to the bulk shear wave. Examples of calculations are also provided for the Rayleigh wave propagating through the homogeneous piezoelectric half-space and layered ZnO(AlN)/diamond/Si structures. It is shown that this kind of surface wave can be used to identify the last unknown length scale parameter of the model.

Keywords

Electroelasticity gradient effects length scale parameters anisotropy bulk wave surface wave dispersion relations spatial dispersion piezoelectric films ZnO AlN

1. Introduction

Fitting of the continuum dispersion relations to the lattice dynamics calculations is one of the effective methods to determine additional material constants of generalized continuum theories. This method has been proposed within the simplified strain gradient elasticity [1, 2] and then used in different theories, including the general Mindlin–Toupin second gradient elasticity [3 –7] and flexoelectricity [8]. A similar approach has been used to analyze different lattice metamaterials at the level of structural mechanics [9, 10]. Quasi-static atomistic calculations have also been used to calibrate static second gradient models [11, 12].

In the present paper, we extend the lattice-dynamics-based calibration method for the piezoelectric materials and structures, which are described within the second gradient electro-elasticity theory [13, 14]. To the best of the author’s knowledge, this was not done previously and only elastic materials have been considered in such approaches [1, 4, 5]. In the present study, the variant of non-classical theory with strain and inertia gradient effects is considered. The electric field gradient effects [13, 15] and all the coupling effects apart from the standard piezoelectricity are out of consideration. Our choice is dictated by the need for model simplification, such that all non-classical phenomena can be described using a minimal number of additional parameters. In the present case, we use the model that allows correct description of the spatial dispersion effects for relatively high-frequency waves in the elastic and piezoelectric crystals. For example, the models with electric field gradients [15] would not allow us to describe such effects in purely elastic materials.

Note that the inertia gradient effects should be obligatorily included within the strain gradient theories to correctly describe the normal spatial dispersion [14, 16]. An additional fourth-order tensor of micro-inertia moduli arises in such theories [14, 16 –18]. In the present paper, it is shown that for crystals with hexagonal and cubic symmetry, the model contains at least seven and four different length scale parameters, respectively. Simplifications in the constitutive equations are assumed for the sixth-order tensor of gradient moduli according to previously proposed approaches [19]. The comparison between dispersion relations for the bulk waves in the considered continuum model and corresponding lattice dynamics data validates these results. Namely, it is found that the models with reduced number of length scale parameters cannot describe with appropriate accuracy the lattice dynamics data in a wide frequency range.

Examples of calculations are provided in this paper for the piezoelectric hexagonal crystals ZnO and AlN that are widely used in different applications in sensors, filters, resonators, and so on [20, 21]. The lattice dynamics data for the phonon dispersion relations in these crystals is taken from literature [22, 23]. An example of the application of the calibrated model is presented for the layered ZnO(AlN)/diamond/Si structures, for which the length scale parameters of cubic diamond and silicon crystals are also found. All calculations are performed under quasi-electrostatic assumption.

2. Second gradient electro-elasticity theory

We consider a piezoelectric body that occupies the domain $Ω$ with smooth boundary $\partial Ω$ . The field equations and boundary conditions of the considered model can be obtained by applying the least action principle $\partial S = 0$ on the action functional of the following form:

\begin{matrix} S = \int_{0}^{t} \int_{Ω} L d v d t, L = T - H, \\ T = \frac{1}{2} ρ {\overset{\cdot}{u}}_{i} {\overset{\cdot}{u}}_{i} + \frac{1}{2} κ_{ijkl} {\overset{\cdot}{ε}}_{ij} {\overset{\cdot}{ε}}_{kl}, \\ H = \frac{1}{2} C_{ijkl} ε_{ij} ε_{kl} + \frac{1}{2} A_{ijklmn} ε_{ij, k} ε_{lm, n} - e_{kij} ε_{ij} E_{k} - \frac{1}{2} ϵ_{ij} E_{i} E_{j}, \end{matrix}

(1)

where $L$ is the Lagrangian density function; $T$ is the kinetic energy density; $H$ is the electric enthalpy density; $u_{i}$ is the displacement vector; $ε_{ij} = (u_{i, j} + u_{j, i}) / 2$ is an infinitesimal strain tensor, $u_{i}$ is the displacement vector; $E_{i} = - ϕ_{, i}$ is the electric field vector; $ϕ$ is the electric potential; $ρ$ is the mass density; $κ_{ijkl}$ is the micro-inertia tensor; $C_{ijkl}$ and $A_{ijklmn}$ are the fourth- and sixth-order tensors of the elastic moduli, respectively; $e_{ijk}$ is the third-order tensor of the piezoelectric constants; $ϵ_{ij}$ is the tensor of dielectric permittivity constants; the comma denotes the differentiation with respect to the spatial variables; the superposed dot denotes the differentiation with respect to time; and the repeated indices imply summation.

The model considered is the classical electroelasticity [15] with additional strain gradient and inertia gradient effects according to Mindlin Form II [24]. In the static case, this model is the special variant of a more general theory that was considered previously in literature [25, 26] In the present study, we do not include some other effects such as flexoelectricity, high-grade coupling, and gradient electric field effects. We will try to use such reduced formulation of the model (1) to describe the known lattice dynamics data on the spatial dispersion of high-frequency waves in crystals.

Using a standard variational approach from equation (1), one can obtain the following equations of motion and charge equation in $Ω$ :

\begin{matrix} σ_{ij, j} - μ_{ijk, jk} = ρ {\overset{\cdot\cdot}{u}}_{i} - κ_{ijkl} {\overset{\cdot\cdot}{u}}_{k, lj}, D_{i, i} = 0, \end{matrix}

(2)

where the stress $σ_{ij} = σ_{ji}$ , double stress $μ_{ijk} = μ_{jik}$ , and electric displacement $D_{i}$ are defined by

σ_{ij} = \frac{\partial H}{\partial ε_{ij}} = C_{ijkl} ε_{kl} - e_{kij} E_{k},

(3)

μ_{ijk} = \frac{\partial H}{\partial ε_{ij, k}} = A_{ijklmn} ε_{lm, n},

(4)

D_{i} = - \frac{\partial H}{\partial E_{i}} = e_{ijk} ε_{jk} + ϵ_{ij} E_{j} .

(5)

Natural and essential boundary conditions on $\partial Ω$ also follow from the variational approach and have the form

\begin{matrix} t_{i} = {\bar{t}}_{i} & or u_{i} = {\bar{u}}_{i}, \\ m_{i} = {\bar{m}}_{i} & or \partial_{n} u_{i} = {\bar{ε}}_{i}, \\ q = \bar{q} & or ϕ = \bar{ϕ}, \end{matrix}

(6)

where the superposed bar denotes the prescribed value of the corresponding quantity on $\partial Ω$ , $\partial_{n} u_{i} = u_{i, j} n_{j}$ is the normal gradient of displacements, $n_{i}$ are the components of the external unit normal on $\partial Ω$ , and the traction $t_{i}$ , double traction $m_{i}$ , and surface charge $q$ are defined by [26, 27]:

\begin{matrix} t_{i} = (σ_{ij} - μ_{ijk, k} + κ_{ijkl} {\overset{\cdot\cdot}{u}}_{k, l}) n_{j} - μ_{ijk, l} (δ_{jl} - n_{j} n_{l}) n_{k} - 2 H m_{i}, \\ m_{i} = μ_{ijk} n_{j} n_{k}, \\ q = D_{i} n_{i}, \end{matrix}

(7)

where $δ_{ij}$ is the Kronecker delta; $H = - D_{i} n_{i} / 2$ is the mean curvature of $\partial Ω$ and $D_{i} = (δ_{ij} - n_{i} n_{j}) \partial_{j}$ is the surface gradient.

Continuity conditions on the contact surface between two different materials must be prescribed with respect to the same quantities that arose in equation (6), as follows:

[t_{i}] = [u_{i}] = [m_{i}] = [\partial_{n} u_{i}] = [q] = [ϕ] = 0,

(8)

where the square brackets denote the difference between the values of the enclosed quantity on the two sides of the contact surface.

In the following, we will consider hexagonal and cubic single-crystal materials. For hexagonal crystals, the constitutive tensors of the model can be represented using the Voigt notation as follows:

Elastic moduli:

C_{ijkl} : (\begin{matrix} C_{11} & C_{12} & C_{13} & 0 & 0 & 0 \\ C_{11} & C_{13} & 0 & 0 & 0 \\ C_{33} & 0 & 0 & 0 \\ C_{44} & 0 & 0 \\ sym & C_{44} & 0 \\ C_{66} \end{matrix}),

(9)

where $C_{66} = (C_{11} - C_{12}) / 2$ .

Piezoelectric constants and dielectric permittivities:

e_{ijk} : (\begin{matrix} 0 & 0 & 0 & 0 & e_{15} & 0 \\ 0 & 0 & 0 & e_{15} & 0 & 0 \\ e_{31} & e_{31} & e_{33} & 0 & 0 & 0 \end{matrix}), ϵ_{ij} : (\begin{matrix} ϵ_{11} & 0 & 0 \\ 0 & ϵ_{11} & 0 \\ 0 & 0 & ϵ_{33} \end{matrix}) .

(10)

For the sixth-order tensor of gradient moduli, we will imply the following simplified assumption about its structure that is usually introduced for anisotropic materials [19]:

G_{ijklmn} = L_{kn} C_{ijlm},

(11)

where $L_{kn}$ is tensor of squared length scale parameters, which structure for the transversely isotropic materials and is as follows:

L_{ij} : (\begin{matrix} l_{1}^{2} & 0 & 0 \\ 0 & l_{1}^{2} & 0 \\ 0 & 0 & l_{3}^{2} \end{matrix}) .

(12)

Tensor of micro-inertia coefficients $κ_{ijkl}$ has the fourth order. General symmetry conditions for this tensor within the considered Mindlin Form II are standard:

κ_{ijkl} = κ_{klij} = κ_{jikl} = κ_{ijlk} .

(13)

Using these conditions (13) and also accounting for the hexagonal symmetry of considered crystals, we can state that the micro-inertia tensor has a similar structure to the classical elasticity tensor, i.e., in Voigt notation, it can be represented as follows:

κ_{ijkl} : (\begin{matrix} κ_{11} & κ_{12} & κ_{13} & 0 & 0 & 0 \\ κ_{11} & κ_{13} & 0 & 0 & 0 \\ κ_{33} & 0 & 0 & 0 \\ κ_{55} & 0 & 0 \\ sym & κ_{55} & 0 \\ κ_{66} \end{matrix}),

(14)

where we have five independent constants and it is valid that $κ_{66} = (κ_{11} - κ_{12}) / 2$ .

Thus, the hexagonal crystal within the considered model should be characterized by the following set of material constants:

5 elastic constants: $C_{11}$ , $C_{12}$ , $C_{13}$ , $C_{33}$ , $C_{44}$

3 piezoelectric constants: $e_{31}$ , $e_{33}$ , $e_{15}$

2 dielectric constants: $ϵ_{11}$ , $ϵ_{33}$

2 length scale parameters: $l_{1}$ , $l_{3}$

5 micro-inertia coefficients: $κ_{11}$ , $κ_{12}$ , $κ_{13}$ , $κ_{33}$ , $κ_{55}$

among which the first 10 constants (elastic, piezoelectric, and dielectric) are standard and can be directly evaluated based on ab initio methods [22, 23]. Additional seven parameters are related to the non-classical gradient effects, and they can be found using several approaches developed within static and dynamic formulations of strain gradient elasticity [1, 4, 11, 12]. In the present case, we will use the method of fitting between the dispersion relations found based on the continuum gradient model and the corresponding relations found using ab initio lattice dynamics for the bulk acoustic phonons in crystals.

For the following analysis, it is useful to introduce the length scale parameters for the micro-inertia effects. Accounting for the dimensions of coefficients $κ_{ijkl}$ and classical mass density $ρ$ (see equation (1)), we can define:

κ_{11} = g_{1}^{2} ρ, κ_{12} = g_{2}^{2} ρ, κ_{33} = g_{3}^{2} ρ, κ_{13} = g_{4}^{2} ρ, κ_{55} = g_{5}^{2} ρ, κ_{66} = g_{6}^{2} ρ,

where $g_{i} (i = 1, \dots, 6) [m]$ are the introduced additional length scale parameters for micro-inertia coefficients, and it is valid that $g_{6}^{2} = (g_{1}^{2} - g_{2}^{2}) / 2$ .

As a result, in the presented model, we have seven length scale parameters. Two of them $(l_{1}, l_{2})$ characterize the non-local effects in the quasi-static behavior of material, while the other five $(g_{1}, . . ., g_{5})$ are related to its dynamic response. Note that we considered the simplified gradient theory due to assumption (11). In the general strain gradient theory, there arise five additional length scale parameters even for isotropic materials [24].

In the following, we will also consider elastic crystals of cubic symmetry. The constitutive relations for the elastic cubic materials can be defined in a similar form to the presented one for the piezoelectric hexagonal structures (3)–(5) and (9)–(14), while the following sets of material constants should be used:

3 elastic constants: $C_{11} = C_{33}$ , $C_{12} = C_{13}$ , $C_{44}$

1 dielectric constant: $ϵ_{11} = ϵ_{33}$

1 length scale parameter: $l_{1} = l_{3}$

3 micro-inertia coefficients: $κ_{11} = κ_{33} = g_{1}^{2} ρ$ , $κ_{12} = κ_{13} = g_{4}^{2} ρ$ , $κ_{55} = κ_{66} = g_{6}^{2} ρ$

and all piezoelectric constants should be set to zero for non-piezoelectric structures.

3. Bulk waves: calibration of the model

For the model calibration, we consider the following five uncoupled types of bulk waves that are possible in the hexagonal crystals (corresponding classical solutions can be found in Simionescu-Panait [28]):

1. Longitudinal wave in the plane of isotropy.

This case corresponds to the one-dimensional propagation of pressure wave along, e.g., the $x_{1}$ -axis direction with a non-zero displacement component $u_{1}$ . In this case, the piezoelectric effects do not arise. Motion and charge equations (2) become uncoupled. Substituting the constitutive relations (3)–(5) and (9)–(14) into equation (2) and taking into account that $u_{2} = u_{3} = 0$ , we found the following form of the motion equations:

(1 - l_{1}^{2} \partial_{11}) C_{11} u_{1, 11} = (1 - g_{1}^{2} \partial_{11}) ρ {\overset{\cdot\cdot}{u}}_{1},

(15)

where symbol $\partial_{11} = \partial^{2} / \partial x_{1}^{2}$ denotes the second-order derivative.

Assuming harmonic law for the displacement $u_{1} = U e^{i (k x_{1} - ω t)}$ ( $k$ is the wavenumber, $ω$ is the angular frequency, and $U$ is amplitude), we found the dispersion relations for the considered wave:

ω = v_{1} k, v_{1} = \sqrt{\frac{C_{11} (1 + l_{1}^{2} k^{2})}{ρ (1 + g_{1}^{2} k^{2})}} .

(16)

2. Shear horizontal wave in the plane of isotropy.

Consider the propagation of shear wave along, e.g., axis $x_{1}$ with a non-zero displacement component $u_{2}$ . Uncoupled motion equations can be derived from equation (2), assuming that $u_{1} = u_{3} = 0$ in the following form:

(1 - l_{1}^{2} \partial_{11}) C_{66} u_{2, 11} = (1 - g_{6}^{2} \partial_{11}) ρ {\overset{\cdot\cdot}{u}}_{2} .

(17)

Dispersion relations for the harmonic waves $u_{2} = U e^{i (k x_{1} - ω t)}$ take the following form:

ω = v_{2} k, v_{2} = \sqrt{\frac{C_{66} (1 + l_{1}^{2} k^{2})}{ρ (1 + g_{6}^{2} k^{2})}} .

(18)

3. Longitudinal wave along axis of symmetry.

This case corresponds to the one-dimensional propagation of the pressure wave along the $x_{3}$ -axis direction with a non-zero displacement component $u_{3}$ . The motion and charge equations (2) in terms of displacement and electric potential will take the coupled form:

\begin{matrix} (1 - l_{3}^{2} \partial_{33}) C_{33} u_{3, 33} + e_{33} ϕ_{, 33} = (1 - g_{3}^{2} \partial_{33}) ρ {\overset{\cdot\cdot}{u}}_{3} \\ e_{33} u_{, 33} - ϵ_{33} ϕ_{, 33} = 0 \end{matrix} .

(19)

Considering harmonic electro-elastic wave $u_{3} = U e^{i (k x_{3} - ω t)}$ , $ϕ = Φ e^{i (k x_{3} - ω t)}$ , we obtain the following dispersion relations:

ω = v_{3} k, v_{3} = \sqrt{\frac{C_{33} (1 + l_{3}^{2} k^{2}) + \frac{e_{33}^{2}}{ϵ_{33}}}{ρ (1 + g_{3}^{2} k^{2})}} .

(20)

4. Shear wave along axis of symmetry.

Consider the propagation of shear wave along the $x_{3}$ axis with a non-zero displacement component $u_{1}$ . The governing equations (2) become uncoupled for this mode:

\begin{matrix} (1 - l_{3}^{2} \partial_{33}) C_{44} u_{1, 33} = (1 - g_{5}^{2} \partial_{33}) ρ {\overset{\cdot\cdot}{u}}_{1} \end{matrix} .

(21)

The dispersion relations for the harmonic wave $u_{1} = U e^{i (k x_{3} - ω t)}$ are the following:

ω = v_{4} k, v_{4} = \sqrt{\frac{C_{44} (1 + l_{3}^{2} k^{2})}{ρ (1 + g_{5}^{2} k^{2})}} .

(22)

5. Transverse shear wave propagated in the plane of isotropy.

Finally, consider the propagation of the shear wave along the $x_{1}$ -axis direction with a non-zero displacement $u_{3}$ . The coupled governing equations (2) for this mode take the following form:

\begin{matrix} (1 - l_{1}^{2} \partial_{11}) C_{44} u_{3, 11} + e_{15} ϕ_{, 11} = (1 - g_{5}^{2} \partial_{11}) ρ {\overset{\cdot\cdot}{u}}_{3} \\ e_{15} u_{, 11} - ϵ_{11} ϕ_{, 11} = 0 \end{matrix} .

(23)

The dispersion relations for the harmonic electro-elastic wave $u_{3} = U e^{i (k x_{1} - ω t)}$ , $ϕ = Φ e^{i (k x_{1} - ω t)}$ are the following:

ω = v_{5} k, v_{5} = \sqrt{\frac{C_{44} (1 + l_{1}^{2} k^{2}) + \frac{e_{15}^{2}}{ϵ_{11}}}{ρ (1 + g_{5}^{2} k^{2})}} .

(24)

Thus, for the piezoelectric hexagonal crystal, we have five relations (16), (18), (20), (22), and (24) containing nine standard material constants $(C_{11}, C_{33}, C_{44}, C_{66}, e_{15}, e_{33}, ϵ_{11}, ϵ_{33}, ρ)$ that are known for certain kinds of crystalline materials. These relations also contain six length scale parameters $(l_{1}, l_{3}, g_{1}, g_{3}, g_{5}, g_{6})$ that should be found based on the fitting between the continuum and atomistic models.

For the calibration of the strain gradient elasticity theory for the elastic cubic crystals, we will consider the longitudinal and shear waves, whose dispersion relations are similar to the solutions (16) and (18). For cubic crystal in these two relations, we have three standard parameters $(C_{11}, C_{66}, ρ)$ and three additional length scales $(l_{1}, g_{1}, g_{6})$ .

The classical solutions for the constant phase velocities of the form $v = \sqrt{C_{ij} / ρ}$ follow from the obtained solutions for the case of small values of $l_{i} k$ and $g_{i} k$ that correspond to the relatively low-frequency range of elastic waves. The influence of gradient effects results in the non-constant values of the wave phase velocity in the high-frequency ranges (see Figure 1). This effect is known as spatial dispersion [15]. For natural single-crystal materials, the phase velocities of short elastic waves become lower, and strong dispersion arises in the THz range of frequencies. Therefore, using ab initio or experimental data on the specific spatial dispersion in crystalline materials, we can choose the values of the length scale parameters of the second gradient model to adjust the continuum dispersion relations to the atomistic one. In such a way, we can calibrate the continuum model and identify its length scale parameters. One of the most challenging issues here is to use the rigorous model with minimal number of additional length scale parameters that can be uniquely identified based on the lattice dynamics data for different kinds of waves.

Figure 1.

Illustration for the influence of the length scale parameters on the dispersion relations for the longitudinal wave along the axis of symmetry in hexagonal piezoelectric crystal (solution (20), $v_{0} = \sqrt{C_{33} / ρ}$ , $k_{33} = e_{33} / \sqrt{ϵ_{33} C_{33}} = 0.5$ ).

Up to date, a lot of data exist on the dispersion relations for different crystalline materials obtained based on the ab initio lattice dynamics [22, 23] and validated experimentally [29]. Although the full phonon dispersion for the crystal contains several acoustic branches and a lot of optical branches, for calibration of the second gradient continuum theories, we can use only the former. Optical phonons are out of consideration here, and their continuum description could be performed using only the models with additional internal degrees of freedom [10, 30]. Thus, in the rest of this section, we will use the known lattice dynamics data on the dispersion relations of acoustic phonons in several examples of anisotropic piezoelectric and elastic crystals.

Fitting of the continuum model will be performed using the derived solutions (16), (18), (20), (22), (24) and the lattice dynamics dispersion relations for the elastic waves in crystals. The correspondence between the derived continuum solutions and standard notations for the phonon band structure is presented in Table 1. Here, we show the crystallographic directions of wave propagation and corresponding notation in terms of the high symmetry points of the Brillouin zone in hexagonal crystals. For cubic crystals, we fitted solutions (16) and (18) to the lattice dynamics data for $[100]$ the $(Γ - X)$ direction of propagation of LA and TA waves.

Table 1.

Notations in continuum model and in lattice dynamics for hexagonal crystals.

Phase velocity	$v_{1}$	$v_{2}$	$v_{3}$	$v_{4}$	$v_{5}$
Continuum solution	(16)	(18)	(20)	(22)	(24)
Type of wave	LA	TA	LA	TA	TA
Non-zero displacements	$u_{1}$	$u_{2}$	$u_{3}$	$u_{1}$	$u_{3}$
Propagation direction	$x_{1}$	$x_{1}$	$x_{3}$	$x_{3}$	$x_{1}$
Crystallographic direction	$[100]$	$[100]$	$[001]$	$[001]$	$[100]$
Symmetry points	$Γ - M$	$Γ - M$	$Γ - A$	$Γ - A$	$Γ - M$

LA: longitudinal acoustic; TA: transverse acoustic.

The length scale parameters are identified for the hexagonal piezoelectric crystals ZnO and AlN. The phonon dispersion relations obtained based on lattice dynamics for these materials are taken from literature [22, 23]. The elastic, piezoelectric, and dielectric constants have also been calculated in these references, and they were used in the present calculations (Table 2). Wavevector scaling was performed accounting for the Brillouin zone dimensions. Namely, in hexagonal crystals, it is $\bar{Γ M} = 2 π / (\sqrt{3} a)$ , $\bar{Γ A} = π / c$ and in cubic crystal $\bar{Γ X} = 2 π / a$ , where $a$ and $c$ are the lattice parameters listed in Table 2.

Table 2.

Material constants and lattice parameters used in the calculations [22, 23].

Parameter	Dimensions	ZnO	AlN	Diamond	Si
$ρ$	$g / c m^{3}$	11.44	9.74	5.83	13
$C_{11}$		188	376	1054	144
$C_{13}$		109	98	126	53
$C_{33}$	GPa	205	353	1054	144
$C_{44}$		37	113	562	75
$C_{66}$		39	124	562	75
$e_{31}$		−0.537	−0.58	0	0
$e_{33}$	$C / m^{2}$	1.037	1.46	0	0
$e_{15}$		−0.385	−0.289	0	0
$ϵ_{11}$		10.53	8.23	5.83	13
$ϵ_{33}$	$ϵ_{0}$	11.44	9.74	5.83	13
$a$		3.289	3.129	2.527	3.867
$c$	Å	5.307	5.017	2.527	3.867

The results of fitting between continuum solutions and phonon dispersion relations are presented in Figure 2. The dependence of frequency $f = ω / (2 π)$ on the wavenumber $k$ is shown. Lattice dynamics data is shown by dots, and adjusted continuum solutions are shown by lines. The values of the identified length scale parameters are presented in Table 3. The fitting was performed using the least square method in the limited range of frequencies, in which an appropriate accuracy can be achieved. Namely, for piezoelectric crystals and for silicon, the best fit of the continuum model can be obtained up to wavenumber values $k = 6, . . ., 8 n m^{- 1}$ . For the diamond crystal, the good correspondence was achieved up to $k = 12 n m^{- 1}$ (Figure 2(c)). The identified values of the length scales are rather small and they are even smaller than the lattice parameters (see parameters $a$ and $c$ in Table 2). This situation is typical for the strain gradient elasticity of single crystal materials [1, 4, 11, 12]. The standard hypotheses of continuum mechanics (of the “first gradient" theories) become invalid at this level and second gradient effects should be considered to provide the correspondence between the phenomenological continuum modeling and atomistic data.

Figure 2.

The adjusted dispersion relations of the continuum model to the lattice dynamics data [22, 23] for ZnO (a), AlN (b), diamond (c), and Si (d). The dashed lines show the classical continuum solutions with constant phase velocities.

Table 3.

Identified length scale parameters (Å).

Parameter	ZnO	AlN	Diamond	Si
$l_{1}$	0.99	0.86	0.36	0.55
$l_{3}$	0.53	0.51	0.36	0.55
$g_{1}$	1.29	1.14	0.63	0.72
$g_{3}$	0.78	0.58	0.63	0.72
$g_{5}$	1.21	0.96	0.73	1.88
$g_{6}$	1.93	1.46	0.73	1.88

The important result of identification is that all the introduced length scales have different values (see Table 3). Therefore, we cannot use a more simple model with reduced number of additional parameters. Such an example of fitting the model with a single micro-inertia parameter $g = g_{1} = g_{3} = g_{5} = g_{6}$ for ZnO is shown in Figure 3(a). It is seen that for the shear waves, this simplified model cannot provide appropriate accuracy of description for the frequencies higher than $f = 1.5 THz (k > 3 n m^{- 1})$ . The corresponding branches of continuum dispersion relations (yellow, green, and purple lines in Figure 3(a)) are overlapped on the presented plot. Therefore, it is important to use a more general dynamic model with at least five additional micro-inertia length scales introduced in the present study.

Figure 3.

Example of adjusted dispersion relations of simplified continuum models to the lattice dynamics data for ZnO: (a) model with single micro-inertia coefficient $g = g_{1} =$ $g_{3} = g_{5} = g_{6} = 0.85$ Å, $l_{1} = 0.38$ Å, $l_{3} = 0.33$ Å and (b) model with zero elastic length scale parameters: $l_{1} = l_{3} = 0$ , $g_{1} = 0.71$ Å, $g_{3} = 0.7$ Å, $g_{5} = 0.82$ Å, $g_{5} = 1.5$ Å.

On the contrary, it can be assumed that the dispersion relations could be fitted by considering the micro-inertia terms only. Such an example, for which we use zero values of the elastic length scale parameters $l_{1} = l_{3} = 0$ in solutions (16), (18), (20), (22), and (24) is presented in Figure 3(b) for ZnO. It is seen that the accuracy of the continuum model with inertia gradient effects only is also reduced for the shear waves in comparison with the full model presented in Figure 2(a). Thus, in the present case, it is better to use the full formulation of the model with seven length scale parameters. However, other approaches may exist; for example, one can use the general definition for the gradient moduli tensor without hypothesis (11).

4. Surface waves: predictions of the model

In this section, we present the results of the calculations for the surface waves obtained by considering the electro-elasticity model with identified length scale parameters. At first, we obtain a prediction for the dispersion relations of anti-plane shear horizontal waves propagated over the surface of piezoelectric and elastic crystals. These predictions do not contain any unknown parameters, and they are totally defined by classical and gradient constants presented in Tables 2 and 3. For the Rayleigh waves, we obtain the solutions, which contain the last length scale parameter $(g_{4})$ that does not arise in the solutions for the bulk waves. It is shown that it is possible to identify this parameter based on lattice dynamics or the experimental data for the surface phonon dispersion relations.

4.1. Bleustein–Gulyaev wave

Consider a special type of anti-plane surface wave propagated along the $x_{1}$ -axis direction over the half-space polarized along the $x_{3}$ axis. The outer normal of the half-space is oriented along the $x_{2}$ axis (Figure 4(a)). Such an uncoupled mode of the shear wave is allowable for hexagonal crystals [15]. The displacement and the electric potential fields are assumed to be:

u_{1} = u_{2} = 0, u_{3} = w (x_{1}, x_{2}, t), ϕ = ϕ (x_{1}, x_{2}, t) .

(25)

Figure 4.

Coordinate systems and direction of motions in the problems with surface waves: (a) anti-plane shear horizontal wave, (b) Rayleigh wave on a half-space, and (c) Rayleigh wave in a three-layered structure.

It can be shown that the non-trivial equation of motion and the charge equation for the Bleustein–Gulyaev wave in second gradient electro-elasticity are the following [14]:

\begin{matrix} C_{44} (1 - ℓ_{1}^{2} Δ) Δ w + e_{15} Δ ϕ = ρ (1 - g_{5}^{2} Δ) \overset{\cdot\cdot}{w}, \\ e_{15} Δ w - ϵ_{11} Δ ϕ = 0 . \end{matrix}

(26)

Introducing a new field variable $ψ = ϕ - (e_{15} / ϵ_{11}) w$ , we can obtain system of uncoupled equations:

\begin{matrix} v^{2} (1 + k_{15}^{2} - ℓ_{1}^{2} Δ) Δ w = (1 - g_{5}^{2} Δ) \overset{\cdot\cdot}{w}, \\ Δ ψ = 0 . \end{matrix}

(27)

where $k_{15}^{2} = e_{15}^{2} / (ϵ_{11} C_{44})$ is the piezoelectric coupling factor and $v^{2} = C_{44} / ρ$ is the standard bulk shear wave velocity.

Stress-free boundary conditions on the surface $x_{2} = 0$ should be prescribed with respect to traction and double traction that are defined according to equation (7):

\begin{matrix} t_{3} = σ_{32} - μ_{321, 1} - μ_{322, 2} - μ_{312, 1} + g_{5}^{2} ρ {\overset{\cdot\cdot}{w}}_{, 2} \\ = C_{44} ((1 + k_{15}^{2}) w_{, 2} - ℓ_{1}^{2} (2 w_{, 112} + w_{, 222})) + e_{15} ψ_{, 2} + g_{5}^{2} ρ {\overset{\cdot\cdot}{w}}_{, 2} = 0 \\ m_{3} = μ_{322} = C_{44} ℓ_{1}^{2} w_{, 22} = 0 . \end{matrix}

(28)

Electrostatic boundary conditions at $x_{2} = 0$ are prescribed assuming a contact between the piezoelectric half-space and surrounding free space (air or vacuum) such that

\begin{matrix} ϕ = \hat{ϕ}, D_{2} = {\hat{D}}_{2} \end{matrix},

(29)

where $D_{2} = - ϵ_{11} ψ_{, 2}$ is the electric displacement in a half-space; $\hat{ϕ}$ and ${\hat{D}}_{2} = - ϵ_{0} {\hat{ϕ}}_{, 2}$ are the electric potential and electric displacement in free space, respectively; and $ϵ_{0} = 8.85 \cdot 10^{- 12} F / m$ is the dielectric constant.

The solution of equations (27) for the harmonic surface wave with decaying behavior from the surface $(x_{2} \to - \infty)$ and solution for the electric potential in a free space should be presented in the following form [14]:

\begin{matrix} w = \sum_{n = 1}^{2} A_{n} e^{β_{n} k x_{2}} e^{i (k x_{1} - ω t)}, \\ ψ = B e^{k x_{2}} e^{i (k x_{1} - ω t)}, \\ \hat{ϕ} = Φ e^{- k x_{2}} e^{i (k x_{1} - ω t)}, \end{matrix}

(30)

where $β_{n} > 0$ are the attenuation coefficients.

Values of $β_{n}$ and description of the following standard procedure for numerical calculations of dispersion relations $ω (k)$ are provided in Appendix 1. The examples of evaluated dispersion relations for the high-frequency Bleustein–Gulyaev waves in ZnO and AlN crystals are shown in Figure 5 (blue and yellow solid lines). Material constants and length scale parameters for these calculations are taken from Tables 2 and 3. It is found that the notable spatial dispersion effects arise at the frequencies range $f > 5 THz$ when the phase velocities deviate more than 5% from the classical solution:

v_{BG} = \sqrt{\frac{C_{44} (1 + k_{15}^{2})}{ρ} (1 - \frac{k_{15}^{4}}{{(1 - ϵ_{11} / ϵ_{0})}^{2}})} .

Figure 5.

Predicted dispersion relations of high-frequency shear horizontal surface waves in ZnO and AlN crystals (Bleustein–Gulyaev waves) and in cubic diamond and Si crystals (strain gradient effects): (a) frequency $f (k) = \frac{1}{2 π} ω (k)$ and (b) phase velocity $v = \frac{1}{k} ω (k)$ .

The corresponding classical values of Bleustein–Gulyaev phase velocities (2654 m/s for ZnO and 5953 m/s for AlN) are shown by dashed blue and yellow lines in Figure 5(b).

In Figure 5, we also show the calculated dispersion relations for the shear horizontal surface waves in elastic cubic diamond and Si crystals. Such waves are restricted within classical elastodynamics. Nevertheless, it was shown recently that these waves are allowed within the generalized continuum theories with surface, strain gradient, or flexoelectric effects [27, 31]. The presented results correspond to the predictions of strain gradient elasticity theory, which have been obtained for isotropic materials in Eremeyev et al. [27]. In the present study, we used the considered model (25)–(30) and reduced it to the solution for the elastic crystal assuming the absence of piezoelectric effect $(e_{15} \to 0)$ . Notably, in the limiting case of small frequencies, these shear horizontal waves in elastic crystals degenerate to the bulk shear wave (its attenuation coefficient tends to zero [27]. The corresponding values of bulk shear wave velocities (11,514 m/s in diamond and 4467 m/s in Si) are shown by red and green dashed lines in Figure 5(b). Notably, the dispersive nature of these non-classical surface shear waves in elastic crystals is predicted to be much stronger in comparison to the Bleustein–Gulyaev waves in piezoelectric crystals.

4.2. Rayleigh wave

Now consider the Rayleigh waves, i.e., the surface shear waves with transverse and longitudinal motions of particles. These types of waves in piezoelectric materials are widely used in non-destructive testing and in microelectronics [21]. In the present section, we derive the high-frequency dispersion relations within considered electro-elasticity theory accounting for the strain gradient and inertia gradient effects.

We consider the hexagonal piezoelectric half-space with its axis of symmetry $(x_{3})$ oriented along the outer surface normal. The wave propagates along $x_{1}$ -axis direction (Figure 4(b)). The non-zero components of displacements are $u_{1}$ (longitudinal motions) and $u_{3}$ (transverse motions). It can be shown that in this case the two non-trivial motion equations and charge equation (2) can be reduced to the following form:

\begin{matrix} (1 - Δ_{L}) (C_{11} u_{1, 11} + (C_{13} + C_{44}) u_{3, 31} + C_{44} u_{1, 33}) + (e_{31} + e_{15}) ϕ_{, 13} \\ = ρ ({\overset{\cdot\cdot}{u}}_{1} - g_{1}^{2} {\overset{\cdot\cdot}{u}}_{1, 11} - (g_{4}^{2} + g_{5}^{2}) {\overset{\cdot\cdot}{u}}_{3, 31} - g_{5}^{2} {\overset{\cdot\cdot}{u}}_{1, 33}) \\ (1 - Δ_{L}) (C_{33} u_{3, 33} + (C_{13} + C_{44}) u_{1, 13} + C_{44} u_{3, 11}) + e_{33} ϕ_{, 33} + e_{15} ϕ_{, 11} \\ = ρ ({\overset{\cdot\cdot}{u}}_{3} - g_{3}^{2} {\overset{\cdot\cdot}{u}}_{3, 33} - (g_{4}^{2} + g_{5}^{2}) {\overset{\cdot\cdot}{u}}_{1, 13} - g_{5}^{2} {\overset{\cdot\cdot}{u}}_{3, 11}) \\ (e_{15} + e_{31}) u_{1, 13} + e_{15} u_{3, 11} + e_{33} u_{3, 33} - ϵ_{11} ϕ_{, 11} - ϵ_{33} ϕ_{, 33} = 0, \end{matrix}

(31)

where $Δ_{L} = ℓ_{1}^{2} \frac{\partial^{2}}{\partial x_{1}^{2}} + ℓ_{3}^{2} \frac{\partial^{2}}{\partial x_{3}^{2}}$ is a modified 2D Laplace operator.

Mechanical boundary conditions should be specified on the free surface of half-space $(x_{3} = 0)$ with respect to the normal and tangent components of tractions and double tractions:

\begin{matrix} t_{1} = σ_{13} - μ_{133, 3} - μ_{113, 1} - μ_{131, 1} + g_{5}^{2} ρ ({\overset{\cdot\cdot}{u}}_{1, 3} + {\overset{\cdot\cdot}{u}}_{3, 1}) = 0, \\ t_{3} = σ_{33} - μ_{333, 3} - μ_{331, 1} + ρ (g_{4}^{2} {\overset{\cdot\cdot}{u}}_{1, 1} + g_{3}^{2} {\overset{\cdot\cdot}{u}}_{3, 3}) = 0, \\ m_{1} = μ_{133} = 0, m_{3} = μ_{333} = 0 . \end{matrix}

(32)

Representation of surface tractions in terms of displacements and electric potential can be obtained by substituting the constitutive relations (3)–(5) into equation (32). We do not present these relations here for brevity. The electric boundary conditions at $x_{3} = 0$ imply the contact with external free space:

\begin{matrix} ϕ = \hat{ϕ}, D_{3} = {\hat{D}}_{3}, \end{matrix}

(33)

where $D_{3} = e_{31} ε_{11} + e_{33} ε_{33} - ϵ_{33} ϕ_{, 3}$ and ${\hat{D}}_{3} = - ϵ_{0} {\hat{ϕ}}_{, 3}$ .

The displacement components and electric potential in the solid material and in the free space are assumed to be the following functions:

\begin{matrix} u_{1} = \sum_{n = 1}^{5} A_{n} e^{β_{n} k x_{2}} e^{i (k x_{1} - ω t)}, u_{3} = \sum_{n = 1}^{5} A_{n} U_{n} e^{β_{n} k x_{2}} e^{i (k x_{1} - ω t)}, \\ ϕ = \sum_{n = 1}^{5} A_{n} Φ_{n} e^{β_{n} k x_{2}} e^{i (k x_{1} - ω t)}, \hat{ϕ} = A_{6} e^{- k x_{3}} e^{i (k x_{1} - ω t)}, \end{matrix}

(34)

where attenuation coefficients $β_{n} > 0$ and constants $U_{n}$ , $Φ_{n}$ $(n = 1, \dots, 5)$ are one set of elementary solutions of the eigenvalue problem defined by governing equations (31); $A_{n}$ $(n = 1, \dots, 6)$ are the unknown constants defined by the boundary conditions (32) and (33).

The dispersion relations $ω (k)$ can be found using a standard procedure similar to those described in Appendix 1 for the shear horizontal waves.

These calculations were performed numerically using material constants from Tables 2 and 3. The single unknown parameter in these calculations was the length scale $g_{4}$ which does not arise in the solutions for the bulk waves and was not identified previously. For calculations, we assumed that this parameter has the values $g_{4} = 0, 0.5 g_{1}, g_{1}$ , or $2 g_{1}$ . The evaluated dispersion relations for the Rayleigh waves in ZnO and AlN crystals are presented in Figure 6(a). It is seen that the value of the length scale parameter $g_{4}$ influences the intensity of spatial dispersion of the high-frequency waves, which becomes stronger for higher values of $g_{4}$ . For the low frequencies, the phase velocities tend to corresponding classical solutions: $v_{R} = 1963 m / s$ in ZnO and $v_{R} = 4841 m / s$ in AlN. Therefore, we can propose an identification approach for the $g_{4}$ parameter based on the fitting of the Rayleigh wave dispersion relations obtained from the considered continuum model to those of lattice dynamics or experiments.

Figure 6.

Predicted dispersion relations $v = \frac{1}{k} ω (k)$ of high-frequency Rayleigh waves in piezoelectric crystals (a) and in three-layered structures (b) for different values of the micro-inertia length scale parameter $g_{4}$ .

Considering the Rayleigh waves in the plane $x_{1} x_{2}$ (instead of $x_{1} x_{3}$ , see Figure 4b), one can obtain additional data for independent validation of material constants in the plane of isotropy of transversely isotropic crystals. Namely, the isotropy condition for the identified length scale parameters $g_{6}^{2} = (g_{1}^{2} - g_{2}^{2}) / 2$ could be checked. The corresponding analytical model for such a problem can be easily derived as the particular case of equations (31)–(34). The check of isotropy condition could be an important issue for the validation of the presented model for the high-frequency ranges and, particularly, for validation of the used simplified hypothesis (11). If a breaking of symmetry occurs, it would mean that more complicated forms of a sixth-order elasticity tensor $G_{ijklmn}$ should be considered.

Finally, we consider the three-layered structures consisting of thin piezoelectric film (ZnO or AlN), diamond sub-layer, and silicon substrate (Figure 4(c)). Such structures are used in the high-frequency surface acoustic wave filters [20]. The axis of symmetry of the film materials is oriented along the surface normal (consider the so-called highly oriented piezoelectric films). The Rayleigh-type wave is propagated along the $x_{1}$ -axis direction with non-zero longitudinal $u_{1} (x_{1}, x_{3})$ and transverse $u_{3} (x_{1}, x_{3})$ components. Governing equations, definitions for tractions, and surface charge for the piezoelectric films have the form (31)–(33). Similar relations for the elastic cubic crystals (diamond and Si) can be obtained from equations (31) to (33) assuming the absence of electro-mechanical coupling and the equivalence between corresponding material constants $(C_{33} = C_{11}, l_{3} = l_{1}, etc .)$ . The solution for the displacements and electric potential in the layers is given by the following:

\begin{matrix} u_{1}^{(s)} = \sum_{n = 1}^{N (s)} A_{n}^{(s)} e^{β_{n}^{(s)} k x_{2}} e^{i (k x_{1} - ω t)}, \\ u_{3}^{(s)} = \sum_{n = 1}^{N (s)} A_{n}^{(s)} U_{n}^{(s)} e^{β_{n}^{(s)} k x_{2}} e^{i (k x_{1} - ω t)}, \\ ϕ^{(s)} = \sum_{n = 1}^{N (s)} A_{n}^{(s)} Φ_{n}^{(s)} e^{β_{n}^{(s)} k x_{2}} e^{i (k x_{1} - ω t)}, \end{matrix}

(35)

where $s = 1, 2, 3$ corresponds to the piezoelectric film, elastic sub-layer, and substrate, respectively; according to the order of governing equations (31), it is valid that $N (1) = N (2) = 10$ , and in the substrate accounting for the radiation conditions, we have $N (3) = 5$ ; constants $A_{n}^{(s)} (s = 1, \dots, 3, n = 1, . . ., 10 (5))$ should be found from the boundary and continuity conditions; attenuation coefficients $β_{n}^{(s)}$ and constants $U_{n}^{(s)}$ , $Φ_{n}^{(s)}$ are the sets of solutions for the eigenvalue problems defined by the governing equations (31) for a piezoelectric film and by the corresponding motion equations for elastic cubic crystals, where we have $Φ_{n}^{(1, 3)} = 1$ due to the absence of electromechanical coupling.

Boundary conditions at the free surface of the piezoelectric film $(z = h_{1} + h_{2})$ are similar to those used in the previous problem with homogeneous half-space (equations (32) and (33)). Continuity conditions are prescribed according to equation (8) at the contact surfaces between piezoelectric film and sub-layer $(z = h_{S})$ and between sub-layer and substrate $(z = 0)$ . In total, we have a system of 26 boundary and continuity conditions with a corresponding number of unknown constants (including the amplitude of the electric potential in the external free space). Equating to zero the determinant of this system, we can find the dispersion relations of the considered problem. Similar to classical solutions, there arise several branches that correspond to different modes of surface waves. In the present study, we evaluate only the first root that corresponds to the fundamental mode. The obtained dispersion relations for the ZnO(AlN)/diamond/Si structures are presented in Figure 6(b). Here, it is seen that the spatial dispersion in three-layered structures may become stronger in comparison with pure piezoelectric materials due to the influence of the substrate properties. In the present case, we consider the silicon substrate, whose length scale parameters are relatively high and provide significant gradient effects in the solution.

5. Conclusion

In this paper, we present a variant of the second gradient electro-elastic model that allows the correct description of the spatial dispersion effects in piezoelectric and elastic anisotropic crystals. We identify the length scale parameters of the model based on fitting its dispersion relations in bulk to the lattice dynamics calculations for the single-crystal materials. The examples of calculations within the calibrated model are presented for the high-frequency surface shear waves in piezoelectric and elastic structures. The presented predictions for the shear horizontal waves do not contain any unknown parameters. These relations can be used to directly predict experimental data (if some other effects such as surface elasticity do not arise [32], etc.). The solutions for the Rayleigh waves contain a single unknown length scale related to the micro-inertia effects. This parameter can be identified based on the lattice dynamics data for the velocities of the surface waves at THz range. A calibrated second gradient model can then be used for the refined analysis of the high-frequency processes in the advanced small-scale piezoelectric and elastic structures within the continuum approaches [33, 34]. An approach similar to the presented one can also be developed for wave propagation analysis in inhomogeneous materials with periodic internal structure (composites and meta-materials), in which the length scale parameters will have the order of the unit cell size and the spatial dispersion effects will arise even for a relatively low-frequency range [35].

Footnotes

Appendix 1 Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the Ministry of Science and Higher Education of the Russian Federation (agreement No. 075-15-2022-1023 dated May 17, 2022).

ORCID iD

Yury Solyaev

References

Maranganti

Sharma

A novel atomistic approach to determine strain-gradient elasticity constants: tabulation and comparison for various metals, semiconductors, silica, polymers and the (ir)relevance for nanotechnologies. J Mech Phys Solids 2007; 55(9): 1823–1852.

Askes

Aifantis

EC.

Gradient elasticity in statics and dynamics: an overview of formulations, length scale identification procedures, finite element implementations and new results. Int J Solids Struct 2011; 48(13): 1962–1990.

Kotoul

Skalka

Profant

, et al. Ab initio aided strain gradient elasticity theory in prediction of nanocomponent fracture. Mech Mater 2019; 136: 103074.

Shodja

Tehranchi

A formulation for the characteristic lengths of FCC materials in first strain gradient elasticity via the Sutton–Chen potential. Philos Mag 2010; 90(14): 1893–1913.

Danescu

Grenet

Continuum strain-gradient elasticity from discrete valence force field model for diamond-like crystals. Int J Fracture 2012; 174(1): 95–102.

Every

Jakata

Perspectives on spatial dispersion in cubic crystals provided by neutron scattering, phonon imaging, picosecond laser ultrasound and lattice dynamics models. Proc Mtgs Acoust 2008; 4(1): 045018.

Seppecher

Alibert

Lekszycki

, et al. Pantographic metamaterials: an example of mathematically driven design and of its technological challenges. Continuum Mech Therm 2019; 31(4): 851–884.

Stengel

Unified ab initio formulation of flexoelectricity and strain-gradient elasticity. Phys Rev B 2016; 93(24): 245107.

Madeo

. Wave propagation in relaxed micromorphic continua: modeling meta-materials exhibiting frequency band-gaps. In: 2015-sustainable industrial processing summit, vol. 1. Mont-Royal, QC, Canada: Flogen Star Outreach, pp. 177–178.

10.

Rosi

Placidi

Auffray

On the validity range of strain-gradient elasticity: a mixed static-dynamic identification procedure. Eur J Mech A-Solids 2018; 69: 179–191.

11.

Gusev

Lurie

SA.

Strain-gradient elasticity for bridging continuum and atomistic estimates of stiffness of binary Lennard–Jones crystals. Adv Eng Mater 2010; 12(6): 529–533.

12.

Lazar

Agiasofitou

Böhlke

Mathematical modeling of the elastic properties of cubic crystals at small scales based on the Toupin–Mindlin anisotropic first strain gradient elasticity. Continuum Mech Therm 2022; 34(1): 107–136.

13.

Arvanitakis

Kalpakides

Hadjigeorgiou

Electric field gradients and spontaneous quadrupoles in elastic ferroelectrics. Acta Mech 2011; 218(3): 269–294.

14.

Solyaev

Lurie

Electric field, strain and inertia gradient effects on anti-plane wave propagation in piezoelectric materials. J Sound Vib 2021; 494: 115898.

15.

Yang

Special topics in the theory of piezoelectricity. Cham: Springer Science & Business Media, 2010.

16.

Metrikine

On causality of the gradient elasticity models. J Sound Vib 2006; 297(3–5): 727–742.

17.

Giorgio

Galantucci

Hamdan

, et al. Wave reflection at a free interface in an anisotropic pyroelectric medium with nonclassical thermoelasticity. Continuum Mech Therm 2016; 28(1–2): 67–84.

18.

Placidi

Rosi

Giorgio

, et al. Reflection and transmission of plane waves at surfaces carrying material properties and embedded in second-gradient materials. Math Mech Solids 2014; 19(5): 555–578.

19.

Polizzotto

Anisotropy in strain gradient elasticity: simplified models with different forms of internal length and moduli tensors. Eur J Mech A-Solids 2018; 71: 51–63.

20.

Assouar

Elmazria

Rioboo

, et al. Modelling of SAW filter based on ZnO/diamond/Si layered structure including velocity dispersion. Appl Surf Sci 2000; 164(1–4): 200–204.

21.

Luo

Nguyen

, et al. Advances in piezoelectric thin films for acoustic biosensors, acoustofluidics and lab-on-chip applications. Prog Mater Sci 2017; 89: 31–91.

22.

De Jong

Chen

Angsten

, et al. Charting the complete elastic properties of inorganic crystalline compounds. Sci Data 2015; 2(1): 1–13.

23.

De Jong

Chen

Geerlings

, et al. A database to enable discovery and design of piezoelectric materials. Sci Data 2015; 2(1): 1–13.

24.

Mindlin

RD.

Microstructure in linear elasticity. Technical report, Department of Civil Engineering and Engineering Mechanics, Columbia University, New York, 1963.

25.

Kalpakides

Agiasofitou

On material equations in second gradient electroelasticity. J Elasticity 2002; 67(3): 205–227.

26.

Iesan

A theory of thermopiezoelectricity with strain gradient and electric field gradient effects. Eur J Mech A-Solids 2018; 67: 280–290.

27.

Eremeyev

Rosi

Naili

Comparison of anti-plane surface waves in strain-gradient materials and materials with surface stresses. Math Mech Solids 2019; 24(8): 2526–2535.

28.

Simionescu-Panait

Waves in strained/polarized media. In: Yang

(ed.) Special topics in the theory of piezoelectricity. Cham: Springer, 2009, pp. 169–246.

29.

Calzolari

Nardelli

MB.

Dielectric properties and Raman spectra of ZnO from a first principles finite-differences/finite-fields approach. Sci Rep 2013; 3(1): 1–6.

30.

Chaki

Eremeyev

Singh

AK.

Surface and interfacial anti-plane waves in micropolar solids with surface energy. Math Mech Solids 2021; 26(5): 708–721.

31.

Eliseev

Morozovska

Glinchuk

, et al. Lost surface waves in nonpiezoelectric solids. Phys Rev B 2017; 96(4): 045411.

32.

Eremeyev

Rosi

Naili

Surface/interfacial anti-plane waves in solids with surface energy. Mech Res Commun 2016; 74: 8–13.

33.

Shen

Qiu

, et al. A low-power circuit for piezoelectric vibration control by synchronized switching on voltage sources. Sensor Actuat A: Phys 2010; 161(1–2): 245–255.

34.

Darleux

Lossouarn

Giorgio

, et al. Electrical analogs of curved beams and application to piezoelectric network damping. Math Mech Solids 2022; 27(4): 578–601.

35.

Demore

Rizzi

Collet

, et al. Unfolding engineering metamaterials design: relaxed micromorphic modeling of large-scale acoustic meta-structures. J Mech Phys Solids 2022; 168: 104995.