Time-harmonic loading over a piezoelectric layered half-space

Abstract

Mathematical modeling of multilayered piezoelectric (PE) ceramic substantially acquires attention due to its distinctive advantages of fast response time, positioning, optical systems, vibration feedback, and sensors, such as deformation and vibration control. As such, fundamental solution of a PE structure is essential. This paper presents three-dimensional (3D) static and dynamic solutions (i.e. Green’s functions) in a multilayered transversally isotropic (TI) PE layered half-space. The uniform vertical mechanical load, vertical electrical displacement, and horizontal mechanical load are applied on the surface of the structure. The novel Fourier-Bessel series (FBS) system of vector functions (which is computationally more powerful and streamlined) and the dual-variable and position (DVP) method are employed to solve the related boundary-value problem. Two systems of first-order ordinary differential equations (i.e. the LM- and N-types) are obtained in terms of the FBS system of vector functions, with these expansion coefficients being the Love numbers. A recursive relation for the expansion coefficients is established by using DVP method that facilitates the combination of two neighboring layers into a new one and minimizes the computational effort to a great extent. The corresponding physical-domain solutions are acquired by applying the appropriate boundary/interface conditions. Several numerical examples pertaining to static and dynamic response are solved, and the efficiency and accuracy of the proposed solutions are validated with the existing results for the reduced cases. The solutions provided could be beneficial to better developments of PE materials, configurations, fabrication, and applications in the future.

1. Introduction

Due to the intrinsic mechanical-electrical coupling effect, piezoelectric (PE) materials stimulate an abundance of interdisciplinary research (Azizi et al., 2016). Developing PE materials with structural systems like advanced composite beams, plates, and shells have enormous applications in emerging fields such as aerospace industries (Gasbarri et al., 2014), military operation (Barrett-Gonzalez and Wolf, 2022), and energy harvesting (Elahi et al., 2019). The designed smart structures include the shape control of large space antennas, active or passive control of vibration (Rao and Sunar, 1994; Yoon and Washington, 1998), miniature positioning devices, for example, micro-robots, medical apparatus, micro-pumps (Muralt et al., 1986), highly efficient micro-machined acoustic devices like microphones and micro speakers (Williams et al., 2012), hearing aids (Zargarpour and Zarifi, 2015), and chemical and biomedical ones (Huang and Lee, 2008; Manbachi and Cobbold, 2011). Recently, extensive studies are conducted on improving PE-based devices by utilizing the coupled electromechanical properties of PE materials (Cao and Yang, 2012; Gao and Zhang, 2020; Li et al., 2017; Yang et al., 2013).

Since the electric and mechanical fields are coupled together, there are many challenges associated with piezo-based devices where precise control of the electric field, mechanical stress, and the potential is needed. As such, various rigorous mathematical models are developed. For instance, Wang et al. (2002) acquired a state-vector method for the axisymmetric deformation in multilayered PE media. Three-dimensional Green’s functions in anisotropic tri-materials were studied by Yang and Pan (2002). Li et al. (2002) studied the stress singularities in PE-bonded structures. Pan (2003) presented a general method for deriving the Green’s functions in three-dimensional (3D) anisotropic PE and multilayered media. With the help of the extended Stroh formalism and two-dimensional (2D) Fourier transforms, Yang et al. (2004) investigated the 3D response in a multilayered anisotropic PE structure under the action of a point force and a point charge. Simplified models under different loading conditions and medium variable parameters were developed extensively to simulate the behavior of the PE structure in the past few decades (Chanda et al., 2023; Liew et al., 2003; Lu et al., 2018; Xiang and Shi, 2009; Zhang and Ma, 2015; Zhang et al., 2019; Zhou et al., 2017).

Various methods in modeling and simulation of the multilayered structure were proposed. These include the exact spectral elements method (Abad and Rouzegar, 2017), the finite element method (Hegendörfer et al., 2022), the space-time backward substitution method (Zhang et al., 2022), the dual reciprocity method (DRM) (Zhan et al., 2023), the meshless boundary element (Power et al., 2017), and the Stroh formalism method (Pan et al., 2023). In dealing with layering, the propagator matrix method (PMM) (Pan and Han, 2004), stiffness matrix method (Wang et al., 2022), precise integration matrix method (PIM) (Zhang et al., 2016), generalized reflection–transmission matrices method (Liu et al., 2020), the wave refraction method (Rafiee-Dehkharghani et al., 2021), the analytic layer-element method (Ai et al., 2022), and the reverberation-ray matrix method (Pao et al., 2007) were proposed.

Inspired by the PIM method, a new method named the dual-variable and position method (DVP) was recently proposed (Liu and Pan, 2018; Liu et al., 2018). This new method is much more advanced to address the time-dependent problem in various isotropic and transversally isotropic (TI) materials.

In this paper, we apply the DVP method along with the Fourier-Bessel series (FBS) system (Pan et al., 2021) to solve the time-harmonic response of TI layered PE half-spaces induced by a surface load or surface electric displacement. While the DVP method is unconditionally stable, the FBS system is much efficient than the traditional integral transform. The reason is that the FBS system only needs simple summation over the discrete zero points of the Bessel functions, whilst, in the integral transform method, numerical integration has to be carried out over each interval between these zero points (i.e. Liu et al., 2016; Moshtagh et al., 2017). Furthermore, since the expansion coefficients are discrete in terms of the FBS system, they can be served as Love numbers which can be saved and repeatedly used later for the calculation of the corresponding Green’s functions (i.e. Zhou et al., 2023). These Love numbers are derived and used to find the response/Green’s function in three typical layered structures proposed in this paper. The accuracy is validated by comparing with existing solutions. Numerical examples indicate the strong dependence of the time-harmonic responses on the input frequency, material layering, and loading types.

2. Description of the boundary-value problem

A multilayered transversally isotropic piezoelectric (TI-PE) structure composed of n different layers over a TI-PE half-space is considered in this paper as shown in Figure 1. A global cylindrical coordinate system (r, θ, z) is introduced in such a way that the (r, θ)-plane is on the surface and the layered half-space is in the positive z-direction. Layer j is bounded by the interfaces z = z_j and z = z_j−1 with a thickness of h_j = z_j−z_j−1. Without losing the generalization, we assume that a uniform time-harmonic load is applied within a circular region of radius $a$ on the top surface of the structure (i.e. at z = z₀). Notice that since the load is proportional to the complex exponential function e^iωt, where ω represents the angular frequency, the time factor e^iωt in the response is suppressed for brevity.

Figure 1.

A TI-PE layered half-space under time-harmonic loads (in both vertical direction along z-axis and horizontal direction in the (x, y)-plane) within a circle of radius a on the surface.

2.1. Fundamental equations

In terms of the cylindrical coordinates, the time-harmonic wave equation for the PE material without body force and body charge (i.e. proportional to e^iωt) is given by

$$\begin{array}{c}\frac{\partial {\sigma }_{\mathrm{rr}}}{\partial r}+\frac{1}{r}\frac{\partial {\sigma }_{r\theta }}{\partial \theta }+\frac{\partial {\sigma }_{\mathrm{rz}}}{\partial z}+\frac{{\sigma }_r-{\sigma }_{\theta }}{2}+\rho {\omega }^2{u}_r=0\hfill \\ \frac{\partial {\sigma }_{r\theta }}{\partial r}+\frac{1}{r}\frac{\partial {\sigma }_{\theta \theta }}{\partial \theta }+\frac{\partial {\sigma }_{\theta z}}{\partial z}+\frac{2{\sigma }_{r\theta }}{r}+\rho {\omega }^2{u}_{\theta }=0\hfill \\ \frac{\partial {\sigma }_{\mathrm{rz}}}{\partial r}+\frac{1}{r}\frac{\partial {\sigma }_{\theta z}}{\partial \theta }+\frac{\partial {\sigma }_{\mathrm{zz}}}{\partial z}+\frac{{\sigma }_{\mathrm{rz}}}{r}+\rho {\omega }^2{u}_z=0\hfill \\ \frac{\partial D_r}{\partial r}+\frac{1}{r}\frac{\partial D_{\theta }}{\partial \theta }+\frac{\partial D_z}{\partial z}+\frac{D_r}{r}=0\hfill \end{array}$$

(1)

where σ _ij (N/m²) (i, j = r, θ, z), D_j (C/m²), ρ (kg/m³), u_i (m), and ω = 2πf (f (Hz = 1/s)) are the stress, electric displacement, mass density, elastic displacement, and angular frequency.

The linear constitutive relation for the TI-PE material (with its symmetry axis along z-axis) in the cylindrical coordinate system is given by

$$\begin{array}{c}{\sigma }_{\mathrm{rr}}=c_{11}{\gamma }_{\mathrm{rr}}+c_{12}{\gamma }_{\theta \theta }+c_{13}{\gamma }_{\mathrm{zz}}-e_{31}{E}_z\hfill \\ {\sigma }_{\theta \theta }=c_{12}{\gamma }_{\mathrm{rr}}+c_{11}{\gamma }_{\theta \theta }+c_{13}{\gamma }_{\mathrm{zz}}-e_{31}{E}_z\hfill \\ {\sigma }_{\mathrm{zz}}=c_{13}{\gamma }_{\mathrm{rr}}+c_{13}{\gamma }_{\theta \theta }+c_{33}{\gamma }_{\mathrm{zz}}-e_{33}{E}_z\hfill \\ {\sigma }_{\mathrm{rz}}=2c_{44}{\gamma }_{\mathrm{rz}}-e_{15}{E}_r\hfill \\ {\sigma }_{\theta z}=2c_{44}{\gamma }_{\theta z}-e_{15}{E}_{\theta }\hfill \\ {\sigma }_{r\theta }=2c_{66}{\gamma }_{r\theta }\hfill \\ D_r=e_{15}{\gamma }_{\mathrm{rz}}+{\eta }_{11}{E}_r\hfill \\ D_{\theta }=e_{15}{\gamma }_{\theta z}+{\eta }_{11}{E}_{\theta }\hfill \\ D_z=e_{31}{\gamma }_{\mathrm{rr}}+e_{31}{\gamma }_{\theta \theta }+e_{33}{\gamma }_{\mathrm{zz}}+{\eta }_{33}{E}_z\hfill \end{array}$$

(2)

where c_ij (N/m²), e_ij (C/m²), and η _ij (C²/N m²) are the elastic constant tensor, PE constant tensor, and dielectric coefficients tensor, respectively; γ _ij (dimensionless) and E_i (V/m) are the strain and electric field; c₆₆ = 0.5 (c₁₁−c₁₂) for the TI material.

The strain components γ _ij and electrical field E_j can be expressed in terms of mechanical displacements u_i and electric potential φ (V) as

$$\begin{array}{c}{\gamma }_{\mathrm{rr}}=\frac{\partial u_r}{\partial r},{\gamma }_{\theta \theta }=\frac{\partial u_{\theta }}{r\partial r}+\frac{u_r}{r},{\gamma }_{\mathrm{zz}}=\frac{\partial u_z}{\partial z}\hfill \\ {\gamma }_{\theta z}=\frac{1}{2}(\frac{\partial u_{\theta }}{\partial z}+\frac{\partial u_z}{r\partial z}),{\gamma }_{\mathrm{rz}}=\frac{1}{2}(\frac{\partial u_z}{\partial r}+\frac{\partial u_r}{\partial z}),\hfill \\ {\gamma }_{r\theta }=\frac{1}{2}(\frac{\partial u_r}{r\partial \theta }+\frac{\partial u_{\theta }}{\partial r}-\frac{u_{\theta }}{r})\hfill \\ E_r=-\frac{\partial \varphi }{\partial r},E_{\theta }=-\frac{1}{r}\frac{\partial \varphi }{\partial \theta },E_z=-\frac{\partial \varphi }{\partial z},\hfill \end{array}$$

(3)

2.2. Boundary and interface conditions

2.2.1. Condition on the surface

It is assumed that, on the surface z = z₀, a uniform load of time-harmonic is applied within the circle of radius a as shown in Figure 1. Then, the boundary condition at z = z₀ (=0) (in term of Cartesian coordinate system, and proportional to e^iωt) is given by

$$\begin{array}{c}{\sigma }_{\mathrm{xz}}=\left\{\begin{array}{c}{p}_{x}r\le a\hfill \\ \multicolumn{1}{c}{0r>0}\end{array}\right\},{\sigma }_{\mathrm{yz}}=\left\{\begin{array}{c}{p}_{y}r\le a\hfill \\ \multicolumn{1}{c}{0r>0}\end{array}\right\},\hfill \\ {\sigma }_{\mathrm{zz}}=\left\{\begin{array}{c}{p}_{z}r\le a\hfill \\ \multicolumn{1}{c}{0r>0}\end{array}\right\},D_z=\left\{\begin{array}{c}{d}_{z}r\le a\hfill \\ \multicolumn{1}{c}{0r>0}\end{array}\right\}\hfill \end{array}$$

(4)

2.2.2. Continuity condition at any interface z = z_k (except for the half-space)

It is assumed that the interface is perfectly bonded. Namely, the mechanical displacements, tractions, electrical displacements, and electric potential are continuous across it, for example, at z = z_k,

$$\begin{array}{c}{u}_i(r,\theta ,z_{k+})=u_i(r,\theta ,z_{k-}),i=r,\theta ,z\hfill \\ {\sigma }_{\mathrm{iz}}(r,\theta ,z_{k+})={\sigma }_{\mathrm{iz}}(r,\theta ,z_{k-}),i=r,\theta ,z\hfill \\ D_z(r,\theta ,z_{k+})=D_z(r,\theta ,z_{k-})\hfill \\ \varphi (r,\theta ,z_{k+})=\varphi (r,\theta ,z_{k-})\hfill \end{array}$$

(5)

In the last homogeneous TI-PE half-space, we require that the solution there is finite. In other words, the elastic displacements and electric potential are zero when z approaches infinity. This condition will be restated later.

3. Fundamental solution in term of the new vector system and the DVP method

3.1. Fourier-Bessel series system of vector functions

To solve the problem, the powerful FBS system of vector functions is introduced (Pan et al., 2021)

$$\begin{array}{c}\mathbf{L}(r,\theta ;{\lambda }_{\mathrm{mk}})={\mathit{e}}_{z}S(r,\theta ;{\lambda }_{\mathrm{mk}})\hfill \\ \mathbf{M}(r,\theta ;{\lambda }_{\mathrm{mk}})=\nabla S=\frac{\partial S}{\partial r}{\mathit{e}}_r+\frac{\partial S}{r\partial \theta }{\mathit{e}}_{\theta }\hfill \\ \mathbf{N}(r,\theta ;{\lambda }_{\mathrm{mk}})=\nabla \times \left(S{\mathit{e}}_z\right)=\frac{\partial S}{r\partial \theta }{\mathit{e}}_r-\frac{\partial S}{\partial r}{\mathit{e}}_{\theta }\hfill \end{array}$$

(6)

where e _r , e _θ , and e _z represent the unit vectors along the r-, θ-, and z-directions of the cylindrical system; λ _mk is the k-th zero of the Bessel function of order m which is scaled by a large value R, that is, $J_m({\lambda }_{\mathrm{mk}}R)=0$ so that the deformation vanishes for r > R. Notice that the vector system L, M, and N is orthogonal.

The scalar function in equation (6) is given by

$$S(r,\theta ;{\lambda }_{\mathrm{mk}})=\frac{1}{\sqrt{2\pi }}{J}_m({\lambda }_{\mathrm{mk}}r){\mathrm{e}}^{\mathrm{i}m\theta },m=0,\pm 1,\pm 2......$$

(7)

which is a solution of the Helmholtz equation

$$\frac{{\partial }^{2}S}{\partial r^2}+\frac{\partial S}{r\partial r}+\frac{{\partial }^{2}S}{r^2\partial {\theta }^2}+{\lambda }_{\mathrm{mk}}^{2}S=0$$

(8)

It is worthy to mention that the FBS system of vector function is motivated by the spectral element method where the approximation is made in the physical domain by truncating the radial distance at a large R. This novel approach is fundamentally different from the traditional Hankel integral transform where approximation is made in the transformed wavenumber domain (Pan et al., 2021). Due to the orthogonality of FBS system (equation (6)), the displacement vector (u), traction vector (t) at z = constant plane, and electric displacement vector (D) can be expanded in terms of the FBS system of vector functions as

$$\begin{array}{l}\mathbf{u}\left(r,\theta ,z\right)=u_r{\mathbf{e}}_r+u_{\theta }{\mathbf{e}}_{\theta }+u_z{\mathbf{e}}_z\\ =\sum _m\sum _k[U_L\left(m,k;z\right)\mathbf{L}\left(r,\theta ;{\lambda }_{mk}\right)+U_M\left(m,k;z\right)\mathbf{M}\left(r,\theta ;{\lambda }_{mk}\right)\\ +U_N\left(m,k;z\right)\mathbf{N}\left(r,\theta ;{\lambda }_{mk}\right)]\end{array}$$

(9)

$$\begin{array}{l}\mathbf{t}\left(r,\theta ,z\right)={\sigma }_{rz}{\mathbf{e}}_r+{\sigma }_{\theta z}{\mathbf{e}}_{\theta }+{\sigma }_{zz}{\mathbf{e}}_z\\ =\sum _m\sum _k[T_L\left(m,k;z\right)\mathbf{L}\left(r,\theta ;{\lambda }_{mk}\right)+T_M\left(m,k;z\right)\mathbf{M}\left(r,\theta ;{\lambda }_{mk}\right)+\\ T_N\left(m,k;z\right)\mathbf{N}\left(r,\theta ;{\lambda }_{mk}\right)]\end{array}$$

(10)

$$\begin{array}{l}\mathbf{D}\left(r,\theta ,z\right)=D_r{\mathbf{e}}_r+D_{\theta }{\mathbf{e}}_{\theta }+D_z{\mathbf{e}}_z\\ =\sum _m\sum _k[D_L\left(m,k;z\right)\mathbf{L}\left(r,\theta ;{\lambda }_{mk}\right)+D_M\left(m,k;z\right)\mathbf{M}\left(r,\theta ;{\lambda }_{mk}\right)\\ +D_N\left(m,k;z\right)\mathbf{N}\left(r,\theta ;{\lambda }_{mk}\right)]\end{array}$$

(11)

whilst the electric potential φ function can be expanded as

$$\varphi (r,\theta ,z)=\sum _m\sum _k\mathrm{\Phi }\left(z\right)S(r,\theta ;{\lambda }_{\mathrm{mk}})$$

(12)

where U_L, U_M, U_N, T_L, T_M, T_N, D_L, D_M, D_N, and $\mathrm{\Phi }$ are the expansion coefficients. Notice that for a given vector or scalar function, its expansion coefficients can be found. For instance, the expansion coefficients of the displacement vector in equation (9) can be expressed as

$$\begin{array}{c}{U}_L(m,k;z)=\frac{1}{N_1(m,k)}\underset{0}{\overset{2\pi }{\int }}\underset{0}{\overset{R}{\int }}\mathbf{u}(r,\theta ,z).\mathbf{L}*(r,\theta ;{\lambda }_{\mathrm{mk}})\mathrm{rdrd}\theta \hfill \\ U_M(m,k;z)=\frac{1}{{\lambda }_{\mathrm{mk}}^2{N}_1(m,k)}\underset{0}{\overset{2\pi }{\int }}\underset{0}{\overset{R}{\int }}\mathbf{u}(r,\theta ,z).\mathbf{M}*(r,\theta ;{\lambda }_{\mathrm{mk}})\mathrm{rdrd}\theta \hfill \\ U_N(m,k;z)=\frac{1}{{\lambda }_{\mathrm{mk}}^2{N}_1(m,k)}\underset{0}{\overset{2\pi }{\int }}\underset{0}{\overset{R}{\int }}\mathbf{u}(r,\theta ,z).\mathbf{N}*(r,\theta ;{\lambda }_{\mathrm{mk}})\mathrm{rdrd}\theta \hfill \end{array}$$

(13)

and the expansion coefficient of the electric potential in equation (12) can be expressed as

$$\varphi (m,k;z)=\frac{1}{N_1(m,k)}\underset{0}{\overset{2\pi }{\int }}\underset{0}{\overset{R}{\int }}\mathrm{\Phi }\left(z\right).\mathbf{S}*(r,\theta ;{\lambda }_{\mathrm{mk}})\mathrm{rdrd}\theta$$

(14)

where * means the complex conjugate, and

$$N_1(m,k)=\frac{R^2}{2}{\left[J_{m+1}\left({\lambda }_{\mathrm{mk}}R\right)\right]}^2=\frac{R^2}{2}{\left[J_{m-1}\left({\lambda }_{\mathrm{mk}}R\right)\right]}^2$$

It is noteworthy, for the axis-symmetric deformation, that is, m = 0, only LM-type vector functions will be involved. However, under a horizontal load, we have m = ±1 and thus both LM-type and N-type are involved.

From the constitutive relations, we find the following relations among the expansion coefficients of the traction, electric displacement, elastic displacement, and electric potential

$$\begin{array}{c}{T}_L=-{\lambda }^2{c}_{13}{U}_M+c_{33}{U\prime }_L+e_{33}\mathrm{\Phi }\prime \hfill \\ T_M=c_{44}({U\prime }_M+U_L)+e_{15}\mathrm{\Phi }\hfill \\ T_N=c_{44}{U\prime }_N\hfill \end{array}$$

(15)

and

$$D_L=-{\lambda }^2{e}_{31}{U}_M+e_{33}U{\prime }_L-{\eta }_{33}\mathrm{\Phi }\prime$$

(16a)

$$\begin{array}{c}{D}_M=e_{15}({U\prime }_M+U_L)-{\eta }_{11}\mathrm{\Phi }\hfill \\ D_N=e_{15}{U\prime }_N\hfill \end{array}$$

(16b)

where prime denotes derivative with respect to the coordinate variable z.

Equations (15) and (16a) contains eight expansion coefficients and needs to be combined with the following three relations (from equations (1) to (3)) to solve them

$$\begin{array}{c}(c_{13}+c_{44})U_L\prime -{\lambda }^2{c}_{11}{U}_M+c_{44}{U}_M″+(e_{15}+e_{31})\mathrm{\Phi }\prime +\rho {\omega }^2{U}_M=0\hfill \\ -{\lambda }^2{c}_{66}{U}_N+c_{44}{U}_N″+\rho {\omega }^2{U}_N=0\hfill \\ -{\lambda }^2{c}_{44}(U_M\prime +U_L)U_L\prime -{\lambda }^2{c}_{13}{U}_M\prime +c_{33}{U}_L″-{\lambda }^2{e}_{15}\mathrm{\Phi }\hfill \\ +e_{33}\mathrm{\Phi }″+\rho {\omega }^2{U}_L=0\hfill \\ -{\lambda }^2(e_{15}+e_{31})U_M\prime +e_{33}{U}_L″-{\lambda }^2{e}_{15}{U}_L-{\eta }_{33}\mathrm{\Phi }″+{\lambda }^2{\eta }_{11}\mathrm{\Phi }=0\hfill \end{array}$$

(17)

Now equations (15), (16a), and (17) contain a total of eight equations for the eight expansion coefficients that is, U_L, U_M, U_N, $\mathrm{\Phi }$ T_L, T_M, T_N, and D_L. They can be converted into a linear system of first-order differential equations. Furthermore, both the LM-type and N-type systems can be decoupled.

For the N-type, we have the following standard linear system of ordinary differential equations

$$\frac{d}{\mathrm{dz}}{\mathbf{W}}_N\left(z\right)={\lambda }_{\mathrm{mk}}{\left[{\mathbf{A}}_N\right]}_{2\times 2}{\mathbf{W}}_N\left(z\right)$$

(18)

For the LM-type, we have

$$\frac{d}{\mathrm{dz}}\mathbf{W}\left(z\right)={\lambda }_{\mathrm{mk}}{\left[\mathbf{A}\right]}_{6\times 6}\mathbf{W}\left(z\right)$$

(19)

where coefficient matrices $\mathbf{A}$ and ${\mathbf{A}}_N$ are listed in Appendix A, and

$$\begin{array}{c}{\mathbf{W}}_N\left(z\right)={\left[\begin{array}{cc}{\lambda }_{\mathrm{mk}}{U}_N\hfill & T_N\hfill \end{array}\right]}^t;\hfill \\ \mathbf{W}\left(z\right)={\left[\begin{array}{cccccc}{U}_L\hfill & {\lambda }_{\mathrm{mk}}{U}_M\hfill & \mathrm{\Phi }\hfill & T_L/{\lambda }_{\mathrm{mk}}\hfill & T_M\hfill & D_L/{\lambda }_{\mathrm{mk}}\hfill \end{array}\right]}^t\hfill \end{array}$$

It is noted here that the LM-type system corresponds to the P-, SV-, and Rayleigh-waves whilst the N-type corresponds to the SH- and Love-waves. If the applied load is axis-symmetric, the N-type solution is zero.

3.2. General solution in each layer

The solutions of equations (18) and (19) can be simply expressed in terms of their eigenvalues and eigenvectors. The general solutions in any layer, say in Layer j bounded by z_j−1 and z_j, that is, $z_{j-1}\le z\le z_j(1\le j\le n)$ can be found as listed below.

For N-type,

$$\begin{array}{c}\left[\begin{array}{c}{\lambda }_{\mathrm{mk}}{U}_N\left(z\right)\hfill \\ \multicolumn{1}{c}{T_N\left(z\right)}\end{array}\right]=\left[\begin{array}{cc}1\hfill & 1\hfill \\ \multicolumn{1}{c}{c_{44}{s}_N}& -c_{44}{s}_N\hfill \end{array}\right]\hfill \\ \left[\begin{array}{cc}{e}^{{\lambda }_{\mathrm{mk}}{s}_N(z-z_j)}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& e^{-{\lambda }_{\mathrm{mk}}{s}_N(z-z_{j-1})}\hfill \end{array}\right]\left[\begin{array}{c}{c}_{N+}\hfill \\ \multicolumn{1}{c}{c_{N-}}\end{array}\right]\hfill \end{array}$$

(20)

where

$$s_N=\sqrt{(c_{66}-\rho {\omega }^2/{{\lambda }_{\mathrm{mk}}}^2)/c_{44}}$$

For LM-type,

$$\begin{array}{c}\left[\begin{array}{c}{\mathbf{U}}_{\mathrm{LM}}\left(z\right)\hfill \\ \multicolumn{1}{c}{{\mathbf{T}}_{\mathrm{LM}}\left(z\right)}\end{array}\right]=\left[\begin{array}{cc}{\mathbf{E}}_{11}\hfill & {\mathbf{E}}_{12}\hfill \\ \multicolumn{1}{c}{{\mathbf{E}}_{21}}& {\mathbf{E}}_{22}\hfill \end{array}\right]\hfill \\ \left[\begin{array}{cc}〈e^{{\lambda }_{\mathrm{mk}}{s}_1(z-z_j)}〉\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 〈e^{{\lambda }_{\mathrm{mk}}{s}_2(z-z_{j-1})}〉\hfill \end{array}\right]\left[\begin{array}{c}{\mathbf{c}}_{\mathbf{*}}\hfill \\ \multicolumn{1}{c}{{\mathbf{c}}_{\mathbf{-}}}\end{array}\right]\hfill \end{array}$$

(21)

where

$$\begin{array}{c}{\mathbf{U}}_{\mathrm{LM}}={\left[\begin{array}{ccc}{U}_L\hfill & {\lambda }_{\mathrm{mk}}{U}_M\hfill & \mathrm{\Phi }\hfill \end{array}\right]}^t;{\mathbf{T}}_{\mathrm{LM}}={\left[\begin{array}{ccc}{T}_L/{\lambda }_{\mathrm{mk}}\hfill & T_M\hfill & D_L/{\lambda }_{\mathrm{mk}}\hfill \end{array}\right]}^t\hfill \\ 〈e^{{\lambda }_{\mathrm{mk}}{s}_{1}z}〉=\mathrm{diag}\left[\begin{array}{ccc}{e}^{{\lambda }_{\mathrm{mk}}{s}_{1}z}\hfill & e^{{\lambda }_{\mathrm{mk}}{s}_{2}z}\hfill & e^{{\lambda }_{\mathrm{mk}}{s}_{3}z}\hfill \end{array}\right];\hfill \\ 〈e^{{\lambda }_{\mathrm{mk}}{s}_{2}z}〉=\mathrm{diag}\left[\begin{array}{ccc}{e}^{{\lambda }_{\mathrm{mk}}{s}_{4}z}\hfill & e^{{\lambda }_{\mathrm{mk}}{s}_{5}z}\hfill & e^{{\lambda }_{\mathrm{mk}}{s}_{6}z}\hfill \end{array}\right]\hfill \end{array}$$

and $\left[\mathbf{E}\right]$ is the eigenmatrix.

Notice that for the N-type in equation (20), s_N has a positive real part. For the LM type in equation (21), the first three eigenvalues all have positive real parts and the other three have negative real parts. Thus, we have Re(s₁)≥…≥Re(s₄).

3.3. Solution in the layered half–space

3.3.1. DVP matrix and recursive relations

In order to obtain the time-harmonic solution in a multi-layered TI-PE half-space, the DVP method is adopted. This new DVP method is unconditionally stable, flexible, and powerful for problems having high frequency, material anisotropy, imperfect interface, and many more. First, equations (20) and (21) can be written on the interfaces of Layer j

$$\begin{array}{c}\left[\begin{array}{c}{\lambda }_{\mathrm{mk}}{U}_N\left(z_j\right)\hfill \\ \multicolumn{1}{c}{T_N\left(z_j\right)}\end{array}\right]=\left[\begin{array}{cc}1\hfill & e^{-{\lambda }_{\mathrm{mk}}{s}_N{h}_j}\hfill \\ \multicolumn{1}{c}{c_{44}{s}_N}& -c_{44}{s}_N{e}^{-{\lambda }_{\mathrm{mk}}{s}_N{h}_j}\hfill \end{array}\right]\left[\begin{array}{c}{c}_{N+}\hfill \\ \multicolumn{1}{c}{c_{N-}}\end{array}\right]\hfill \\ \left[\begin{array}{c}{\lambda }_{\mathrm{mk}}{U}_N\left(z_{j-1}\right)\hfill \\ \multicolumn{1}{c}{T_N\left(z_{j-1}\right)}\end{array}\right]=\left[\begin{array}{cc}{e}^{-{\lambda }_{\mathrm{mk}}{s}_N{h}_j}\hfill & 1\hfill \\ \multicolumn{1}{c}{c_{44}{s}_N{e}^{-{\lambda }_{\mathrm{mk}}{s}_N{h}_j}}& -c_{44}{s}_N\hfill \end{array}\right]\left[\begin{array}{c}{c}_{N+}\hfill \\ \multicolumn{1}{c}{c_{N-}}\end{array}\right]\hfill \end{array}$$

(22)

$$\begin{array}{c}\left[\begin{array}{c}{\mathbf{U}}_{\mathrm{LM}}\left(z_j\right)\hfill \\ \multicolumn{1}{c}{{\mathbf{T}}_{\mathrm{LM}}\left(z_j\right)}\end{array}\right]=\left[\begin{array}{cc}{\mathbf{E}}_{11}\hfill & {\mathbf{E}}_{12}〈e^{{\lambda }_{\mathrm{mk}}{s}_2{h}_j}〉\hfill \\ \multicolumn{1}{c}{{\mathbf{E}}_{21}}& {\mathbf{E}}_{22}〈e^{{\lambda }_{\mathrm{mk}}{s}_2{h}_j}〉\hfill \end{array}\right]\left[\begin{array}{c}{\mathbf{c}}_{\mathbf{*}}\hfill \\ \multicolumn{1}{c}{{\mathbf{c}}_{-}}\end{array}\right]\hfill \\ \left[\begin{array}{c}{\mathbf{U}}_{\mathrm{LM}}\left(z_{j-1}\right)\hfill \\ \multicolumn{1}{c}{{\mathbf{T}}_{\mathrm{LM}}\left(z_{j-1}\right)}\end{array}\right]=\left[\begin{array}{cc}{\mathbf{E}}_{11}〈e^{-{\lambda }_{\mathrm{mk}}{s}_1{h}_j}〉\hfill & {\mathbf{E}}_{12}\hfill \\ \multicolumn{1}{c}{{\mathbf{E}}_{21}〈e^{-{\lambda }_{\mathrm{mk}}{s}_1{h}_j}〉}& {\mathbf{E}}_{22}\hfill \end{array}\right]\left[\begin{array}{c}{\mathbf{c}}_{\mathbf{*}}\hfill \\ \multicolumn{1}{c}{{\mathbf{c}}_{-}}\end{array}\right]\hfill \end{array}$$

(23)

where h_j = z_j−z_j−1. Then, by eliminating the constants $c_{N+},c_{N-},{\mathbf{c}}_{+}$ , and ${\mathbf{c}}_{-}$ in equations (22) and (23), the relations between the expansion coefficients for Layer j can be derived as

$$\begin{array}{c}\left[\begin{array}{c}{\lambda }_{\mathrm{mk}}{U}_N\left(z_{j-1}\right)\hfill \\ \multicolumn{1}{c}{T_N\left(z_j\right)}\end{array}\right]=\left[\begin{array}{cc}{N}_{11}^j\hfill & N_{12}^j\hfill \\ \multicolumn{1}{c}{N_{21}^j}& N_{22}^j\hfill \end{array}\right]\left[\begin{array}{c}{\lambda }_{\mathrm{mk}}{U}_N\left(z_j\right)\hfill \\ \multicolumn{1}{c}{T_N\left(z_{j-1}\right)}\end{array}\right]\hfill \\ \left[\begin{array}{c}{\mathbf{U}}_{\mathrm{LM}}\left(z_{j-1}\right)\hfill \\ \multicolumn{1}{c}{{\mathbf{T}}_{\mathrm{LM}}\left(z_j\right)}\end{array}\right]=\left[\begin{array}{cc}{\mathbf{S}}_{11}^j\hfill & {\mathbf{S}}_{12}^j\hfill \\ \multicolumn{1}{c}{{\mathbf{S}}_{21}^j}& {\mathbf{S}}_{22}^j\hfill \end{array}\right]\left[\begin{array}{c}{\mathbf{U}}_{\mathrm{LM}}\left({\mathrm{z}}_{\mathrm{j}}\right)\hfill \\ \multicolumn{1}{c}{{\mathbf{T}}_{\mathrm{LM}}\left(z_{j-1}\right)}\end{array}\right]\hfill \end{array}$$

(24)

where the submatrices are given by

$$\begin{array}{c}{\mathbf{N}}_{\mathrm{ij}}^j\mathbf{=}\left[\begin{array}{cc}{\mathbf{N}}_{11}^j\hfill & {\mathbf{N}}_{12}^j\hfill \\ \multicolumn{1}{c}{{\mathbf{N}}_{13}^j}& {\mathbf{N}}_{14}^j\hfill \end{array}\right]=\left[\begin{array}{cc}{e}^{-{\lambda }_{\mathrm{mk}}{s}_N{h}_j}\hfill & 1\hfill \\ \multicolumn{1}{c}{s_N{c}_{44}}& -s_N{c}_{44}{e}^{-{\lambda }_{\mathrm{mk}}{s}_N{h}_j}\hfill \end{array}\right]\hfill \\ {\left[\begin{array}{cc}1\hfill & e^{-{\lambda }_{\mathrm{mk}}{s}_N{h}_j}\hfill \\ \multicolumn{1}{c}{s_N{c}_{44}{e}^{-{\lambda }_{\mathrm{mk}}{s}_N{h}_j}}& -s_N{c}_{44}\hfill \end{array}\right]}^{-1}\hfill \\ {\mathbf{S}}_{\mathrm{ij}}^j\mathbf{=}\left[\begin{array}{cc}{\mathbf{S}}_{11}^j\hfill & {\mathbf{S}}_{12}^j\hfill \\ \multicolumn{1}{c}{{\mathbf{S}}_{13}^j}& {\mathbf{S}}_{14}^j\hfill \end{array}\right]=\left[\begin{array}{cc}{\mathbf{E}}_{11}〈e^{-{\lambda }_{\mathrm{mk}}{s}_1{h}_j}〉\hfill & {\mathbf{E}}_{12}\hfill \\ \multicolumn{1}{c}{{\mathbf{E}}_{21}}& {\mathbf{E}}_{22}〈e^{{\lambda }_{\mathrm{mk}}{s}_2{h}_j}〉\hfill \end{array}\right]\hfill \\ {\left[\begin{array}{cc}{\mathbf{E}}_{11}\hfill & {\mathbf{E}}_{12}〈e^{{\lambda }_{\mathrm{mk}}{s}_2{h}_j}〉\hfill \\ \multicolumn{1}{c}{{\mathbf{E}}_{21}〈e^{-{\lambda }_{\mathrm{mk}}{s}_1{h}_j}〉}& {\mathbf{E}}_{22}\hfill \end{array}\right]}^{-1}\hfill \end{array}$$

Similarly, the expansion coefficients for Layer j + 1 are

$$\begin{array}{c}\left[\begin{array}{c}{\lambda }_{\mathrm{mk}}{U}_N\left(z_j\right)\hfill \\ \multicolumn{1}{c}{T_N\left(z_{j+1}\right)}\end{array}\right]=\left[\begin{array}{cc}{N}_{11}^{j+1}\hfill & N_{12}^{j+1}\hfill \\ \multicolumn{1}{c}{N_{21}^{j+1}}& N_{22}^{j+1}\hfill \end{array}\right]\left[\begin{array}{c}{U}_N\left(z_{j+1}\right)\hfill \\ \multicolumn{1}{c}{{\lambda }_{\mathrm{mk}}{T}_N\left(z_j\right)}\end{array}\right]\hfill \\ \left[\begin{array}{c}{\mathbf{U}}_{\mathrm{LM}}\left(z_j\right)\hfill \\ \multicolumn{1}{c}{{\mathbf{T}}_{\mathrm{LM}}\left(z_{j+1}\right)}\end{array}\right]=\left[\begin{array}{cc}{\mathbf{S}}_{11}^{j+1}\hfill & {\mathbf{S}}_{12}^{j+1}\hfill \\ \multicolumn{1}{c}{{\mathbf{S}}_{21}^{j+1}}& {\mathbf{S}}_{22}^{j+1}\hfill \end{array}\right]\left[\begin{array}{c}{\mathbf{U}}_{\mathrm{LM}}\left(z_{j+1}\right)\hfill \\ \multicolumn{1}{c}{{\mathbf{T}}_{\mathrm{LM}}\left(z_j\right)}\end{array}\right]\hfill \end{array}$$

(25)

Assuming that the interface between the two layers is well bonded (i.e. the elastic displacements, tractions, electric displacement in z-direction, and electric potential are continuous), we then find the following important recursive relation from Layer j to Layer j + 1 as

$$\begin{array}{c}\left[\begin{array}{c}{\lambda }_{\mathrm{mk}}{U}_N\left(z_{j-1}\right)\hfill \\ \multicolumn{1}{c}{T_N\left(z_{j+1}\right)}\end{array}\right]=\left[\begin{array}{cc}{N}_{11}^{j:j+1}\hfill & N_{12}^{j:j+1}\hfill \\ \multicolumn{1}{c}{N_{21}^{j:j+1}}& N_{22}^{j:j+1}\hfill \end{array}\right]\left[\begin{array}{c}{\lambda }_{\mathrm{mk}}{U}_N\left(z_{j+1}\right)\hfill \\ \multicolumn{1}{c}{T_N\left(z_{j-1}\right)}\end{array}\right]\hfill \\ \left[\begin{array}{c}{\mathbf{U}}_{\mathrm{LM}}\left(z_{j-1}\right)\hfill \\ \multicolumn{1}{c}{{\mathbf{T}}_{\mathrm{LM}}\left(z_{j+1}\right)}\end{array}\right]=\left[\begin{array}{cc}{\mathbf{S}}_{11}^{j:j+1}\hfill & {\mathbf{S}}_{12}^{j:j+1}\hfill \\ \multicolumn{1}{c}{{\mathbf{S}}_{21}^{j:j+1}}& {\mathbf{S}}_{22}^{j:j+1}\hfill \end{array}\right]\left[\begin{array}{c}{\mathbf{U}}_{\mathrm{LM}}\left(z_{j+1}\right)\hfill \\ \multicolumn{1}{c}{{\mathbf{T}}_{\mathrm{LM}}\left(z_{j-1}\right)}\end{array}\right]\hfill \end{array}$$

(26)

where the superscripts j:j+1 denote the propagated submatrix from Layer j to Layer j+1 given by

$\begin{array}{c}{N}_{11}^{j:j+1}=N_{11}^j{N}_{11}^{j+1}+N_{11}^j{N}_{12}^{j+1}{[1-N_{21}^j{N}_{12}^{j+1}]}^{-1}{N}_{21}^j{N}_{11}^{j+1}\hfill \\ N_{12}^{j:j+1}=N_{12}^j+N_{11}^j{N}_{12}^{j+1}{[1-N_{21}^j{N}_{12}^{j+1}]}^{-1}{N}_{22}^j\hfill \\ N_{21}^{j:j+1}=N_{21}^j+N_{22}^{j+1}{[1-N_{21}^j{N}_{12}^{j+1}]}^{-1}{N}_{21}^j{N}_{11}^{j+1}\hfill \\ N_{22}^{j:j+1}=N_{22}^{j+1}{[1-N_{21}^j{N}_{12}^{j+1}]}^{-1}{N}_{22}^j\hfill \\ {\mathbf{S}}_{11}^{j:j+1}={\mathbf{S}}_{11}^j{\mathbf{S}}_{11}^{j+1}+{\mathbf{S}}_{11}^j{\mathbf{S}}_{12}^{j+1}{[\mathit{I}-{\mathbf{S}}_{21}^j{\mathbf{S}}_{12}^{j+1}]}^{-1}{\mathbf{S}}_{21}^j{\mathbf{S}}_{11}^{j+1}\hfill \\ {\mathbf{S}}_{12}^{j:j+1}={\mathbf{S}}_{12}^j+{\mathbf{S}}_{11}^j{\mathbf{S}}_{12}^{j+1}{[\mathit{I}-{\mathbf{S}}_{21}^j{\mathbf{S}}_{12}^{j+1}]}^{-1}{\mathbf{S}}_{22}^j\hfill \\ {\mathbf{S}}_{21}^{j:j+1}={\mathbf{S}}_{21}^j+{\mathbf{S}}_{22}^{j+1}{[\mathit{I}-{\mathbf{S}}_{21}^j{\mathbf{S}}_{12}^{j+1}]}^{-1}{\mathbf{S}}_{21}^j{\mathbf{S}}_{11}^{j+1}\hfill \\ {\mathbf{S}}_{22}^{j:j+1}={\mathbf{S}}_{22}^{j+1}{[\mathit{I}-{\mathbf{S}}_{21}^j{\mathbf{S}}_{12}^{j+1}]}^{-1}{\mathbf{S}}_{22}^j\hfill \end{array}$

The new recursive relation (i.e. equation (26)) can be propagated in the layered structure, as long as the interfaces are perfect, or well-bonded. By applying it repeatedly from the bottom layer to the top or vice versa, the unknown expansion coefficients at the top and bottom interfaces can be connected. By utilizing the boundary conditions given, the involved known coefficients can be solved. After that, equation (26) can be utilized to find the expansion coefficients at any depth.

4. Boundary conditions

4.1. Boundary conditions on the surface in terms of expansion coefficients

For the given time-harmonic load within the circle of radius a on the surface (Figure 1), the following expansion coefficients in term of FBS system are found as (Pan et al., 2021)

$$\begin{array}{c}{T}_L\left(z_0\right)=\frac{\sqrt{2\pi }{p}_z{J}_1\left({\lambda }_{0k}a\right)a}{{\lambda }_{0k}{N}_1},m=0\hfill \\ D_L\left(z_0\right)=\frac{\sqrt{2\pi }{d}_z{J}_1\left({\lambda }_{0k}a\right)a}{{\lambda }_{0k}{N}_1},m=0\hfill \\ T_M\left(z_0\right)=\frac{\sqrt{\pi }}{2}\frac{a{J}_1\left({\lambda }_{1k}a\right)}{N_2}(\pm p_x-\mathrm{i}{p}_y),m=\pm 1\hfill \\ T_N\left(z_0\right)=\frac{\sqrt{\pi }}{2}\frac{a{J}_1\left({\lambda }_{1k}a\right)}{N_2}(-\mathrm{i}{p}_x\mp p_y),m=\pm 1\hfill \end{array}$$

(27)

The following three separate cases are considered in this paper.

• Case 1: Under uniform vertical load p_z, we have

$${\mathbf{T}}_{\mathrm{LM}}\left(z_0\right)={\left[\begin{array}{ccc}\frac{\sqrt{2\pi }{p}_z{J}_1\left({\lambda }_{0k}a\right)a}{{\lambda }_{0k}{N}_1}\hfill & 0\hfill & 0\hfill \end{array}\right]}^t$$

• Case 2: Under uniform vertical electric displacement d_z, we have

$${\mathbf{T}}_{\mathrm{LM}}\left(z_0\right)={\left[\begin{array}{ccc}0\hfill & 0\hfill & \frac{\sqrt{2\pi }{d}_z{J}_1\left({\lambda }_{0k}a\right)a}{{\lambda }_{0k}{N}_1}\hfill \end{array}\right]}^t$$

• Case 3: Under uniform horizontal load (p_x, p_y) in both x- and y-directions, we have

$${\mathbf{T}}_{\mathrm{LM}}\left(z_0\right)={\left[\begin{array}{ccc}0\hfill & \frac{\sqrt{\pi }}{2}\frac{a{J}_1\left({\lambda }_{1k}a\right)}{N_2}(\pm p_x-\mathrm{i}{p}_y)\hfill & 0\hfill \end{array}\right]}^t,\mathrm{for}m=1$$

4.2. Expansion coefficients of displacements and electric potential on the surface

With the given general surface traction expansion coefficients and the derived recursive relation, we can propagate the expansion coefficients from the surface z₀ to the last layer interface z_n to obtain the following relations

$$\begin{array}{c}\left[\begin{array}{c}{\lambda }_{\mathrm{mk}}{U}_N\left(z_0\right)\hfill \\ \multicolumn{1}{c}{T_N\left(z_n\right)}\end{array}\right]=\left[\begin{array}{cc}{N}_{11}^{j:n}\hfill & N_{12}^{j:n}\hfill \\ \multicolumn{1}{c}{N_{21}^{j:n}}& N_{22}^{j:n}\hfill \end{array}\right]\left[\begin{array}{c}{\lambda }_{\mathrm{mk}}{U}_N\left(z_n\right)\hfill \\ \multicolumn{1}{c}{T_N\left(z_0\right)}\end{array}\right]\hfill \\ \left[\begin{array}{c}{\mathbf{U}}_{\mathrm{LM}}\left(z_0\right)\hfill \\ \multicolumn{1}{c}{{\mathbf{T}}_{\mathrm{LM}}\left(z_n\right)}\end{array}\right]=\left[\begin{array}{cc}{\mathbf{S}}_{11}^{j:n}\hfill & {\mathbf{S}}_{12}^{j:n}\hfill \\ \multicolumn{1}{c}{{\mathbf{S}}_{21}^{j:n}}& {\mathbf{S}}_{22}^{j:n}\hfill \end{array}\right]\left[\begin{array}{c}{\mathbf{U}}_{\mathrm{LM}}\left(z_n\right)\hfill \\ \multicolumn{1}{c}{{\mathbf{T}}_{\mathrm{LM}}\left(z_0\right)}\end{array}\right]\hfill \end{array}$$

(28)

On the last interface z = z_n, the unknown coefficients in equation (28) are related as, for the LM-type (similar expression for the N-type)

$${\mathbf{T}}_{\mathrm{LM}}\left(z_n\right)=\left[{\mathbf{E}}_{22}{\left({\mathbf{E}}_{12}\right)}^{-1}\right]{\mathbf{U}}_{\mathrm{LM}}\left(z_n\right)$$

(29)

By substituting equation (29) into equation (28), we have the following relations for different loading cases.

• (1) Under the vertical symmetric load (i.e. for m = 0) we have only the LM-type solution, which is

$$\left[\begin{array}{c}{\mathbf{U}}_{\mathrm{LM}}\left(z_0\right)\hfill \\ \multicolumn{1}{c}{{\mathbf{E}}_{22}{\left({\mathbf{E}}_{12}\right)}^{-1}{\mathbf{U}}_{\mathrm{LM}}\left(z_n\right)}\end{array}\right]=\left[\begin{array}{cc}{\mathbf{S}}_{11}^{j:n}\hfill & {\mathbf{S}}_{12}^{j:n}\hfill \\ \multicolumn{1}{c}{{\mathbf{S}}_{21}^{j:n}}& {\mathbf{S}}_{22}^{j:n}\hfill \end{array}\right]\left[\begin{array}{c}{\mathbf{U}}_{\mathrm{LM}}\left(z_n\right)\hfill \\ \multicolumn{1}{c}{{\mathbf{T}}_{\mathrm{LM}\left(0\right)}\left(z_0\right)}\end{array}\right]$$

for which the solution on the surface is given by

$${\mathbf{U}}_{\mathrm{LM}}\left(z_0\right)=[{\mathbf{S}}_{11}^{j:n}{({\mathbf{E}}_{22}{\left({\mathbf{E}}_{12}\right)}^{-1}-{\mathbf{S}}_{21}^{j:n})}^{-1}{\mathbf{S}}_{22}^{j:n}+{\mathbf{S}}_{12}^{j:n}]{\mathbf{T}}_{\mathrm{LM}\left(0\right)}\left(z_0\right)$$

(30)

where

$$\begin{array}{c}{\mathbf{T}}_{\mathrm{LM}\left(0\right)}\left(z_0\right)={\left[\begin{array}{ccc}\frac{\sqrt{2\pi }{p}_z{J}_1\left({\lambda }_{0k}a\right)a}{{\lambda }_{0k}^2{N}_1}\hfill & 0\hfill & 0\hfill \end{array}\right]}^t;\hfill \\ {\mathbf{T}}_{\mathrm{LM}\left(0\right)}\left(z_0\right)={\left[\begin{array}{ccc}0\hfill & 0\hfill & \frac{\sqrt{2\pi }{d}_z{J}_1\left({\lambda }_{0k}a\right)a}{{\lambda }_{0k}^2{N}_1}\hfill \end{array}\right]}^t\hfill \end{array}$$

• (2) Under the horizontal load (i.e. for m = ±1), we have both N-type and LM-type solutions, which are given below

For m = 1, we have

$$\begin{array}{c}\left[\begin{array}{c}\left\{{\lambda }_{1k}{U}_N\left(z_0\right)\right\}\hfill \\ \multicolumn{1}{c}{E_{22}{\left(E_{12}\right)}^{-1}\left\{{\lambda }_{1k}{U}_N\left(z_n\right)\right\}}\end{array}\right]=\left[\begin{array}{cc}{N}_{11}^{j:n}\hfill & N_{12}^{j:n}\hfill \\ \multicolumn{1}{c}{N_{21}^{j:n}}& N_{22}^{j:n}\hfill \end{array}\right]\left[\begin{array}{c}\left\{{\lambda }_{1k}{U}_N\left(z_n\right)\right\}\hfill \\ \multicolumn{1}{c}{T_{N(1)}\left(z_0\right)}\end{array}\right]\hfill \\ \left[\begin{array}{c}{\mathbf{U}}_{\mathrm{LM}}\left(z_0\right)\hfill \\ \multicolumn{1}{c}{{\mathbf{E}}_{22}{\left({\mathbf{E}}_{12}\right)}^{-1}{\mathbf{U}}_{\mathrm{LM}}\left(z_n\right)}\end{array}\right]=\left[\begin{array}{cc}{\mathbf{S}}_{11}^{j:n}\hfill & {\mathbf{S}}_{12}^{j:n}\hfill \\ \multicolumn{1}{c}{{\mathbf{S}}_{21}^{j:n}}& {\mathbf{S}}_{22}^{j:n}\hfill \end{array}\right]\left[\begin{array}{c}{\mathbf{U}}_{\mathrm{LM}}\left(z_n\right)\hfill \\ \multicolumn{1}{c}{{\mathbf{T}}_{\mathrm{LM}\left(1\right)}\left(z_0\right)}\end{array}\right]\hfill \end{array}$$

with the solution on the surface being

$$\begin{array}{c}\left\{{\lambda }_{1k}{U}_N\left(z_0\right)\right\}=[N_{11}^{j:n}{(E_{22}{\left(E_{12}\right)}^{-1}-N_{21}^{j:n})}^{-1}{N}_{22}^{j:n}+N_{12}^{j:n}]T_{N\left(1\right)}\left(z_0\right)\hfill \\ {\mathbf{U}}_{\mathrm{LM}}\left(z_0\right)=[{\mathbf{S}}_{11}^{j:n}{({\mathbf{E}}_{22}{\left({\mathbf{E}}_{12}\right)}^{-1}-{\mathbf{S}}_{21}^{j:n})}^{-1}{\mathbf{S}}_{22}^{j:n}+{\mathbf{S}}_{12}^{j:n}]{\mathbf{T}}_{\mathrm{LM}\left(1\right)}\left(z_0\right)\hfill \end{array}$$

(31)

where

$${\mathbf{T}}_{\mathrm{LM}\left(1\right)}\left(z_s\right)={\left[\begin{array}{ccc}0\hfill & \sqrt{\frac{\pi }{2}}\frac{a{J}_1\left({\lambda }_{1k}a\right)}{N_2}(p_x-\mathrm{i}{p}_y)\hfill & 0\hfill \end{array}\right]}^t$$

Similarly, for m = −1, we have

$$\begin{array}{c}\left[\begin{array}{c}\left\{{\lambda }_{1k}{U}_N\left(z_0\right)\right\}\hfill \\ \multicolumn{1}{c}{E_{22}{\left(E_{12}\right)}^{-1}\left\{{\lambda }_{1k}{U}_N\left(z_n\right)\right\}}\end{array}\right]=\left[\begin{array}{cc}{N}_{11}^{j:n}\hfill & N_{12}^{j:n}\hfill \\ \multicolumn{1}{c}{N_{21}^{j:n}}& N_{22}^{j:n}\hfill \end{array}\right]\left[\begin{array}{c}\left\{{\lambda }_{1k}{U}_N\left(z_n\right)\right\}\hfill \\ \multicolumn{1}{c}{T_{N(1)}\left(z_0\right)}\end{array}\right]\hfill \\ \left[\begin{array}{c}{\mathbf{U}}_{\mathrm{LM}}\left(z_0\right)\hfill \\ \multicolumn{1}{c}{{\mathbf{E}}_{22}{\left({\mathbf{E}}_{12}\right)}^{-1}{\mathbf{U}}_{\mathrm{LM}}\left(z_n\right)}\end{array}\right]=\left[\begin{array}{cc}{\mathbf{S}}_{11}^{j:n}\hfill & {\mathbf{S}}_{12}^{j:n}\hfill \\ \multicolumn{1}{c}{{\mathbf{S}}_{21}^{j:n}}& {\mathbf{S}}_{22}^{j:n}\hfill \end{array}\right]\left[\begin{array}{c}{\mathbf{U}}_{\mathrm{LM}}\left(z_n\right)\hfill \\ \multicolumn{1}{c}{{\mathbf{T}}_{\mathrm{LM}\left(1\right)}\left(z_0\right)}\end{array}\right]\hfill \end{array}$$

with the solution on the surface being

$$\begin{array}{c}\left\{{\lambda }_{1k}{U}_N\left(z_0\right)\right\}=[N_{11}^{j:n}{(E_{22}{\left(E_{12}\right)}^{-1}-N_{21}^{j:n})}^{-1}{N}_{22}^{j:n}+N_{12}^{j:n}]T_{N(-1)}\left(z_0\right)\hfill \\ {\mathbf{U}}_{\mathrm{LM}}\left(z_0\right)=[{\mathbf{S}}_{11}^{j:n}{({\mathbf{E}}_{22}{\left({\mathbf{E}}_{12}\right)}^{-1}-{\mathbf{S}}_{21}^{j:n})}^{-1}{\mathbf{S}}_{22}^{j:n}+{\mathbf{S}}_{12}^{j:n}]{\mathbf{T}}_{\mathrm{LM}(-1)}\left(z_0\right)\hfill \end{array}$$

(32)

where

$${\mathbf{T}}_{\mathrm{LM}(-1)}\left(z_s\right)={\left[\begin{array}{ccc}0\hfill & \sqrt{\frac{\pi }{2}}\frac{a{J}_1\left({\lambda }_{1k}a\right)}{N_2}(-p_x-\mathrm{i}{p}_y)\hfill & 0\hfill \end{array}\right]}^t$$

Once the expansion coefficients on the surface and last interface are solved, we can just simply carry out the summation (i.e. equations (9)–(12)) to find the frequency-domain solution at any physical point on the surface. After that, we only need to replace z₀ by z in equation (28) in order to find the solutions at different z-levels. We point out again that since these coefficients are discrete numbers, they are the exact Love numbers which can be pre-calculated and saved for later use in the corresponding Green’s function calculation (i.e. Zhou et al., 2023).

5. Numerical studies

Before presenting numerical examples, we first remark that our newly derived solution has been verified with existing solutions in the literature. These results are consistent with those reported in Han et al. (2006), Pan et al. (2021), and Qi et al. (2023). Appendix B contains a few graphs which demonstrate the accuracy of the present new method as compared with those in Pan et al. (2021) and Qi et al. (2023). Then, the present solution is applied to the time-harmonic loads (mechanical and electrical), which are uniformly applied over a circular patch of radius a on the surface of a layered TI-PE half-space structure for an electrically open circuit case. An FBS MATLAB program is written, which is very powerful and can handle easily multiple loads along with material anisotropy. To avoid the singularities, the elastic properties are perturbed away from their real values by multiplying a small material damping factor β. Namely, c_kj, e_kj, and ε _kj are replaced by c_kj (1 + 2iβ), e_kj (1 + 2iβ), and ε _kj (1 + 2iβ) in the numerical calculation. Also in the paper, all numerical results are presented in the dimensionless form, defined as

$\overline{a}$  = a/a; ${\overline{p}}_x$ =p_x/C_max; $\overline{p_z}$  = p_z/C_max; $\overline{d_z}$  = d_z/E_max; $\overline{x}$  = x/a; $\overline{z}$  = z/a; ${\overline{u}}_x$  = u_x/a; ${\overline{u}}_z$  = u_z/a; $\overline{\phi }$  = φ/φ _max

${\overline{\sigma }}_{zz}=$ σ _zz /C_max; $\overline{\omega }$  = ωa√(ρ _max /C_max).

Here φ_max = E_max/a. The value of C_max and E_max for Model 1 is 210.9 × 10⁹ N/m² and 17 C/m², respectively. For Model 2, they are C_max = 174.1 × 10⁹ N/m² and E_max = 20.95 C/m².

To study how the static and dynamic loads affect the piezoelectric response, two models are constructed, along with the three loading cases, that is, vertical electrical load $\overline{d_z}$ , the vertical mechanical loading $\overline{p_z}$ , and the horizontal mechanical loading ${\overline{p}}_x$ . The detailed models are provided in Table 1 and shown by Figure 2, respectively.

Table 1.

Three different loading cases and two different material models where the material properties are listed in Appendix C.

Loading cases ( $\overline{a}$  = 1)	Case 1	Uniform $\overline{d_z}$ is applied
	Case 2	Uniform $\overline{p_z}$ is applied
	Case 3	Uniform ${\overline{p}}_x$ is applied
Material models	Model 1	Four layered model with stacking sequence BaTiO3/PZT-5H/PZT-4/ZnO
	Model 2	Two layered model with stacking sequence (BaTiO3 (50%) + PZT-5H (50%))/(PZT-4 (50%) + ZnO (50%))

Figure 2.

Geometry of the four-layered structure (a) and two-layered structure (b).

5.1. Four-layered model

The time-harmonic response subjected to the circular vertical electrical load $\overline{d_z}$ (=d_z/E_max = 1/E_max), vertical mechanical load $\overline{p_z}$ (=p_z/C_max = 1/C_max) and horizontal mechanical load ${\overline{p}}_x$ (=p_x/C_max = 1/C_max) on the surface of a four layered structure are shown in Figures 3 to 5, respectively. The model used in this study is Model 1, which comprises three layers with thickness h₁/a, h₂/a, and h₃/a (=0.1, 0.2, and 0.2, respectively). These results highlight the importance of understanding the relationship between the input frequency and the resulting behavior of the layered system. For instance, neglecting the frequency-dependent behavior can lead to flawed or insufficient analyses, and thus, can potentially result in errors in system design and performance. Therefore, it is necessary to incorporate frequency-dependent behavior while studying the piezoelectric functionality and reliability when under time-dependent loads.

Figure 3.

Variation of dimensionless displacement ${\overline{u}}_x$ (a) and electric potential $\overline{\phi }$ (b), from (0, 0, 0) to (0.7, 0, 0), and displacement ${\overline{u}}_z$ (c) and normal stress ${\overline{\sigma }}_{\mathrm{zz}}$ (d), from (0, 0, 0) to (0, 0, 0.7), induced by the vertical electrical load $\overline{d_z}$ (=d_z/E_max = 1/E_max) Here C_max = 210.9 × 10⁹ N/m² and E_max = 17 C/m².

Figure 4.

Variation of the dimensionless displacement ${\overline{u}}_x$ (a) and electric potential $\overline{\phi }$ (b), from (0, 0, 0) to (0.7, 0, 0), and displacement ${\overline{u}}_z$ (c), normal stress ${\overline{\sigma }}_{\mathrm{zz}}$ (d), from (0, 0, 0) to (0, 0, 0.7), due to the vertical mechanical load $\overline{p_z}$ (=p_z/C_max = 1/C_max) Here C_max = 210.9 × 10⁹ N/m² and E_max = 17 C/m².

Figure 5.

Variation of the dimensionless displacement ${\overline{u}}_x$ (a) and electric potential $\overline{\phi }$ (b), from (0, 0, 0) to (0.7, 0, 0), and displacement ${\overline{u}}_z$ (c) and normal stress ${\overline{\sigma }}_{\mathrm{zz}}$ (d), from (0, 0, 0) to (0, 0, 0.7), due to the horizontal mechanical load ${\overline{p}}_x$ (=p_x/C_max = 1/C_max). Here C_max = 210.9 × 10⁹ N/m² and E_max = 17 C/m².

5.1.1. Under vertical electric displacement

Figure 3(a) and (b) depicts the spatial distribution of the dimensionless mechanical displacement ${\overline{u}}_x$ and electric potential function $\overline{\phi }$ along the surface line (varying from $\overline{x}$  = 0 to 0.7), caused by the vertical electric displacement $\overline{d_z}$ It is observed that, with increasing frequency or increasing $\overline{x}$ , the magnitude of ${\overline{u}}_x$ increases, and this trend becomes clearer when frequency $\overline{\omega }$ is near 0.9 and $\overline{x}$ near 0.7. As for $\overline{\phi }$ , its variation is opposite, namely, with increasing frequency or increasing $\overline{x}$ , its value decreases. Figure 3(c) and (d) exhibits the variation of displacement ${\overline{u}}_z$ and normal stress ${\overline{\sigma }}_{\mathrm{zz}}$ along z-axis. It is observed that while the magnitude of ${\overline{u}}_z$ decreases with increasing depth, the variation of $\overline{\phi }$ is peculiar along $\overline{z}$ -axis when frequency is relatively high. It should be also noted that when frequency $\overline{\omega }$  = 0.6 or higher, there is a maximum stress value in this depth range with its location depending on the given frequency (Figure 3(d)). At $\overline{\omega }$  = 0.6, its maximum is 1.0119 × 10^–3 at $\overline{z}$  = 0.62, whist at $\overline{\omega }$  = 0.9, its maximum is 0.83826 × 10^–3 at $\overline{z}$  = 0.29. This indicates that under a time-harmonic electric displacement, the maximum stress value and its location are sensitive to both the frequency and the layered structure.

5.1.2. Under vertical mechanical load

Figure 4(a) shows that the magnitude of the dimensionless displacement ${\overline{u}}_x$ increases linearly with increasing distance $\overline{x}$ , except for $\overline{\omega }$  = 0.9, where an opposite trend is observed. This indicates that the elastic displacement is sensitive to the input frequency, particularly when it is high. Comparing to Figure 4(a), it appears that the frequency has a reverse effect on the electrical potential as shown in Figure 4(b). For a given frequency, its maximum value is at the loading center and it then decreases with increasing distance $\overline{x}$ . Figure 4(c) shows the variation of ${\overline{u}}_z$ along $\overline{z}$ -axis, which are all piecewise in linear. It is interesting to note that the slope of the curves is discontinuous and such a discontinuity is more apparent when frequency is high (i.e. for $\overline{\omega }$  = 0.9). The normal stress ${\overline{\sigma }}_{\mathrm{zz}}$ variation along $\overline{z}$ -axis is presented in Figure 4(d) where its boundary value is unit one. This partially verifies that our formulation and the corresponding code are correct. Comparing to Figure 4(c) for the displacement variation, we notice that the normal stress is not piecewise linear but in a higher-order polynomial. We further observe that for a fixed $\overline{z}$ , ${\overline{\sigma }}_{\mathrm{zz}}$ increases with increasing frequency, except for $\overline{\omega }$  = 0.9 where its value is significantly reduced compared to the lower-frequency results. This implies that there exists certain higher-frequency range over which the stress values can be substantially reduced, for frequency-mediated stress control.

5.1.3. Under horizontal mechanical load

Figure 5(a) and (b) show the variation of the dimensionless displacement ${\overline{u}}_x$ and electric potential $\overline{\phi }$ along $\overline{x}$ -axis for different frequencies, induced by the horizontal mechanical load ${\overline{p}}_x$ . While the results for cases of static and $\overline{\omega }$  = 0.3 are very close to each other, those at higher-frequencies are substantially different. This can be observed in Figure 5(a) where the magnitude of ${\overline{u}}_x$ at $\overline{\omega }$  = 0.9 is only about one-fourth to one-third of that at the static case.

Figure 5(c) shows the variation of ${\overline{u}}_z$ along $\overline{z}$ -axis. Compared with Figure 4(c), we note that ${\overline{u}}_z$ is more sensitive to the horizontal load where the slope discontinuity is clear, even at lower frequency and static case. This indicates that ${\overline{u}}_z$ under the horizontal load could be more correlated to the material layering. It is further noticed that, when frequency is high, the displacement changes from compressive to tensile. As in Figure 5(d) for the normal stress, it is interesting that in the first layer, the response is nearly the same under different frequencies and that there is a minimum within this layer. After that, its value increases with increasing depth $\overline{z}$ (for a given frequency).

5.2. Two-layered model

We now study the time-harmonic response in the “equivalent two-layered half-space,” that is, Model 2. Notice that the material properties are from those in the four-layered half-space by taking the two-layer average. The applied loads are the same, namely, a uniform vertical electric load $\overline{d_z}$ (=d_z/E_max = 1/E_max), a vertical mechanical load $\overline{p_z}$ (=p_z/C_max = 1/C_max), and a horizontal mechanical load ${\overline{p}}_x$ (=p_x/C_max = 1/C_max) on the surface of the two-layered half-space, as shown in Figure 2(b). We discuss the results below one by one.

5.2.1. Under vertical electric displacement

Figure 6(a) to (d) depict the variation of the dimensionless elastic displacement ${\overline{u}}_x$ and electric potential $\overline{\phi }$ along $\overline{x}$ -axis, and of the mechanical displacement ${\overline{u}}_z$ and normal stress ${\overline{\sigma }}_{\mathrm{zz}}$ along z-axis.

Figure 6.

Variation of the dimensionless displacement ${\overline{u}}_x$ (a) and electric potential $\overline{\phi }$ (b), from (0, 0, 0) to (0.5, 0, 0), and displacement ${\overline{u}}_z$ (c), normal stress ${\overline{\sigma }}_{\mathrm{zz}}$ (d), from (0, 0, 0) to (0, 0, 0.5), due to the vertical electrical load $\overline{d_z}$ (=d_z/E_max = 1/E_max). Here C_max = 174.1 × 10⁹ N/m² and E_max = 20.95 C/m².

Compared to those in the four-layered half-space in Figure 3(a) to (d), we observe that, while some patterns are similar, significant differences exist as clearly shown in Figure 6(b) and (c). Notably, Figure 6(b) reveals that, for fixed $\overline{x}$ , increasing frequency leads to a rise in the electric potential $\overline{\phi }$ which is opposite to the trend in Figure 3(b). This indicates that a higher frequency would excite a larger electric voltage in the equivalent two-layered half-space, potentially for a larger electric energy. In Figure 6(c), we observe that, as frequency increases, the magnitude of displacement ${\overline{u}}_z$ decreases, which is not evident in Figure 3(c). Also, the distribution of the normal stress is more complicated in Figure 6(d) than that in Figure 3(d).

5.2.2. Under vertical mechanical load

Figure 7(a) to (d) show the distribution of dimensionless mechanical displacement ${\overline{u}}_x$ and electric potential $\overline{\phi }$ along $\overline{x}$ -axis, and the mechanical displacement ${\overline{u}}_z$ and normal stress ${\overline{\sigma }}_{\mathrm{zz}}$ along $\overline{z}$ -axis, due to the uniform vertical mechanical load $\overline{p_z}$ . It can be observed that the variation of the displacements ( ${\overline{u}}_x$ and ${\overline{u}}_z$ ) and ${\overline{\sigma }}_{\mathrm{zz}}$ is very similarly to that in Figure 4(a), (c), and (d), with the exception of the electric potential as discussed below.

Figure 7.

Variation of the dimensionless displacement ${\overline{u}}_x$ (a) and electric potential $\overline{\phi }$ (b), from (0, 0, 0) to (0.5, 0, 0), and displacement ${\overline{u}}_z$ (c) and normal stress ${\overline{\sigma }}_{\mathrm{zz}}$ (d), from (0, 0, 0) to (0, 0, 0.5), due to the vertical mechanical load $\overline{p_z}$ (=p_z/C_max = 1/C_max). Here C_max = 174.1 × 10⁹ N/m² and E_max = 20.95 C/m².

Figure 7(b) shows a significant impact of frequency on the magnitude of the electric potential. Namely, an increase in frequency could cause the electric potential $\overline{\phi }$ to rise throughout the entire $\overline{x}$ -range, which is much different from the trend observed in Figure 4(b). It has the remarkable aptitude to store a significant amount of energy, which in turn causes an increase in the electric potential function. These results evince that the two-layered structure could be better equipped to handle and boost the electric potential $\overline{\phi }$ as compared to the corresponding four-layered structure. This remarkable characteristic of the two-layered material structure opens up new possibilities for its application in advanced technologies that solicit high-performance electrical properties.

5.2.3. Under horizontal mechanical load

Figure 8(a) to (d) depict the distribution of dimensionless mechanical displacement ${\overline{u}}_x$ and electric potential $\overline{\phi }$ along $\overline{x}$ -axis, and the mechanical displacement ${\overline{u}}_z$ and normal stress ${\overline{\sigma }}_{\mathrm{zz}}$ along $\overline{z}$ -axis, due to the horizontal load of magnitude ${\overline{p}}_x$ . Here, it can be observed that, the results are very similarly to those in Figure 5(a) to (d).

Figure 8.

Variation of the dimensionless displacement ${\overline{u}}_x$ (a) and electric potential $\overline{\phi }$ (b), from (0, 0, 0) to (0.5, 0, 0), and the displacement ${\overline{u}}_z$ (c) and normal stress ${\overline{\sigma }}_{\mathrm{zz}}$ (d), from (0, 0, 0) to (0, 0, 0.5), due to the horizontal mechanical load ${\overline{p}}_x$ (=p_x/C_max = 1/C_max). Here C_max = 174.1 × 10⁹ N/m² and E_max = 20.95 C/m².

6. Conclusion

The time-harmonic response of a transversely isotropic piezoelectric layered half-space induced by surface loading is solved. Two different material models, four-layered and two-layered, are considered (each constituted by a different class of PE materials) to understand the complexity of the responses in the coupled structure. The novel and computationally powerful Fourier-Bessel series (FBS) system of vector functions and the dual-variable and position (DVP) method are employed to solve the boundary-value problem accurately and efficiently. Via numerical examples, the following conclusions can be drawn:

• The effect of input frequency on the horizontal displacement ${\overline{u}}_x$ (in $\overline{x}$ -direction) and normal stress ${\overline{\sigma }}_{\mathrm{zz}}$ (in $\overline{z}$ -direction) are similar in the two models under the vertical electrical displacement. However, the induced electric potential function $\overline{\phi }$ (in $\overline{x}$ -direction) and vertical displacement ${\overline{u}}_z$ are significantly different in these two models. Furthermore, the reduced equivalent two-layer model tends to be more prominent in amplifying the electric potential and vertical displacement, as compared to the corresponding four-layer model. This indicates an alternative way for achieving high voltage output for possible better PE device design.

• Under a vertical mechanical load, the effect of frequency on the displacement ${\overline{u}}_x$ (in $\overline{x}$ -direction), displacement ${\overline{u}}_z$ (in $\overline{z}$ -direction), and normal stress ${\overline{\sigma }}_{\mathrm{zz}}$ (in $\overline{z}$ -direction) are similar in these two models. However, in the four-layered structure, the magnitude of the electric potential $\overline{\phi }$ can be both increased and decreased with varying frequency.

• Under a horizontal mechanical load, the response behavior is similar to those under a vertical mechanical load. However, their variations along $\overline{z}$ -direction are more sensitive to the layered structure, indicating that a horizontal surface load could be more suited in the corresponding inverse problem (i.e. inverting the layered profile).