Modelling and Investigation of the Dynamic Behavior of a Penny-Shaped Interface Crack in Piezoelectric Bimaterials

Abstract

In this section, the dynamic propagation behavior of a penny-shaped interface crack in piezoelectric bimaterials is analyzed. The objective of this paper is to use the boundary conditions of the penny-shaped interface crack to study the dynamic propagation of the crack under the action of load, so as to provide some valuable implications for the fracture mechanics of the piezoelectric bimaterials and simulate the interface crack between piezoelectric bimaterials, it is necessary to establish a suitable model and give appropriate boundary conditions according to the actual situation. The elastic displacement and potential equations are constructed according to the structural characteristics of the circular crack. In the case of a given displacement or stress, the Laplace transform and Hankel transform are used to simplify the problem into an integral equation with unknown functions. According to the boundary conditions, the corresponding unknowns are obtained, and the closed solution is derived. The results show that the fracture toughness of a penny-shaped interface crack in piezoelectric bimaterials is related to the thickness of the material, the impact time, the material characteristics, and the electric field. At the same time, it can be found that different materials have different roles in the crack propagation, so it is very important to study the crack opening displacement (COD) intensity factor of the crack for safety design.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

1.1. Background

Piezoelectric materials are widely used to make all kinds of transducers and sensors because of their piezoelectric effects. They are commonly used in electronic, excited light, ultrasonic, hydroacoustic, microacoustic, red external, navigation, biological, and other technical fields. However, the brittleness of piezoelectric ceramic material in mechanical properties will lead to the generation and expansion of cracks from the stress concentrations, leading to the failure of components. It is necessary to analyze the crack propagation process of piezoelectric materials to improve the operational performance of high-voltage electric ceramic elements and predict their active life.

In order to meet the practical needs, experts have used experimental, theoretical, and simulation methods to study the influence factors of piezoelectric materials in the different situations [1, 2]. The forced vibration of a piezoelectric plate with initial stress and the factors affecting the dynamic stability of a prestressed piezoelectric plate are studied [3, 4].

Based on the above situation, research on the fracture mechanics of piezoelectric ceramics has received attention from the experts [5–13]. Influence of incomplete bonding on the dynamic response of prestressed sandwich plate–strip with elastic layers and a piezoelectric core was studied by Daşdemir [14]. The influence of poling direction and imperfection defects of the prestressed system with a piezoelectric core bonded to elastic faces was analyzed [15]. In most cases, the crack is simulated by a penny shape [16, 17]. A 3D problem of a half-space with penny-shaped cracks has been studied [18]. Zikun [19] solved a penny-shaped crack problem and obtained analytical expressions of the stress field and electric displacement field near the crack tip. Using the Dugdale hypothesis and Hankel transform theory, Danyluk et al. [20] analyzed the penny-shaped crack in a thick transversely isotropic elastic layer. Kogan et al. [21] obtained the stress intensity factors of a penny-shaped crack. The plasticity of a penny-shaped Dugdale crack tip in an infinite elastic medium has also been estimated [22].

A penny-shaped crack has been considered in a piezoelectric medium to analyze the coupling behavior [23]. Penny-shaped cracks in 3D piezoelectric media are analyzed by the extended displacement discontinuity method [24]. The Dugdale plastic zone of a penny-shaped crack under axisymmetric loading has been analyzed [25]. Wu at al. [26] studied a penny-shaped crack in a piezoelectric layer sandwiched between two elastic layers with the electrical saturation and mechanical yielding zones. The crack propagation behavior of an infinite 1D hexagonal piezoelectric quasicrystal plate with a penny-shaped dielectric crack has been analyzed [27].

The crack propagation will change into a dynamic behavior with time under an explosion or impact load. Wang et al. [28] investigated the penny-shaped interface crack configuration in orthotropic multilayers under dynamic torsional loading by utilizing Laplace transform and the Hankel transforms technique. The dynamic behavior of a penny-shaped crack in a magnetoelectroelastic materials has been analyzed in detail [29]. A point force method was proposed for obtaining the transient response of dynamic penny-shaped cracks in multilayer sandwich composites [30]. The dynamic behavior of a magnetoelectroelastic material with a penny-shaped dielectric crack under impact loading has been determined [31]. A mathematical expression for the dynamic behavior (open function circular) of a permeable penny-shaped crack in an infinitely porous elastic solid has been presented [32].

In order to overcome the weak characteristics of the piezoelectric effect of a single piezoelectric material and better expand the piezoelectric effect, the piezoelectric double material has emerged, which is better applied in intelligent structures such as transducers, sensors, and drivers, and it has also produced great application and economic value. At the same time, because of its superiority, it has also been studied and analyzed by many experts and scholars [33–39]. Under loading conditions, the interface of piezoelectric bimaterials is prone to fracture due to surface offset. In addition, crack propagation becomes a dynamic behavior when piezoelectric bimaterials are impacted in practice [40]. There are few studies on the dynamic change of a penny-shaped interface crack in piezoelectric bimaterials. Therefore, the dynamic analysis of a penny-shaped interface crack in piezoelectric bimaterials is very important in designing practical engineering applications.

It is well-known that defects such as interface cracks, holes, and dislocations seriously affect the mechanical behavior and change their strength. In addition, the interface crack boundary conditions are of great significance for selecting fracture criteria and predicting crack propagation. In addition, the penny-shaped interface crack, which is the research target of this paper, is a common type of crack, and the dynamic propagation to be studied is also the most practical and research value. However, the dynamic propagation of a penny-shaped interface crack is difficult in experiments and simulations. Therefore, in order to meet the needs of development and safety, the objective of this paper is to use the boundary conditions of the penny-shaped interface crack to study the dynamic propagation of the crack under the action of load, so as to provide some valuable implications for the fracture mechanics of the piezoelectric bimaterials.

1.2. Outline

After the introduction, the model and constitutive equations are described. Section 3 discusses the calculation methods under the Laplace transform and Hankel transform. In Section 4, the dynamic field intensity factors of crack-tip propagation are derived. In Section 5, numerical results are obtained based on the model and calculation. The final section presents some conclusions drawn from this study.

2. Problem Statement and Formulation

In Figure 1, a penny-shaped interface crack made in piezoelectric materials one and two is considered. The thickness of piezoelectric Materials 1 and 2 are $h_1$ and $h_2$, respectively. The radius of the penny-shaped interface crack is $a$. Let us analyze the case of the poling direction along the $z$-axis in polar the coordinates $r,\theta ,z$. The center of the crack is the origin of coordinates, and its region is $r\le a,z=0$. The names of the parameters can be seen in nomenclature.

[figure(s) omitted; refer to PDF]

According to Tiersten [41], thess balance equation (divergence equation) of the dynamic behavior of a penny-shaped interface crack in piezoelectric bimaterials can be expressed as follows:$$\begin{array}{}\text{(1)}& \begin{array}{c}\frac{\partial {\tau }_{r{r}_k}}{\partial r}+\frac{\partial {\tau }_{r{z}_k}}{\partial z}+\frac{{\tau }_{r{r}_k}-{\tau }_{\theta {\theta }_k}}{r}={\rho }_k\frac{{\partial }^2{u}_{r_k}}{\partial t^2}\\ \frac{\partial {\tau }_{r{z}_k}}{\partial r}+\frac{\partial {\tau }_{z{z}_k}}{\partial z}+\frac{{\tau }_{r{z}_k}}{r}={\rho }_k\frac{{\partial }^2{u}_{z_k}}{\partial t^2}\hfill \\ \frac{\partial D_{r_k}}{\partial r}+\frac{\partial D_{z_k}}{\partial z}+\frac{D_{r_k}}{r}=0\hfill \end{array},k=1,2.\end{array}$$

The constitutive equations (equations of state) [41] can be written as follows:$$\begin{array}{}\text{(2)}& \begin{array}{c}{\tau }_{r{r}_k}=c_{11}^k\frac{\partial u_{r_k}}{\partial r}+c_{12}^k\frac{u_{r_k}}{r}+c_{13}^k\frac{\partial u_{z_k}}{\partial z}+e_{31}^k\frac{\partial {\phi }_k}{\partial z},\\ {\tau }_{r{z}_k}=c_{44}^k\frac{\partial u_{z_k}}{\partial r}+c_{44}^k\frac{\partial u_{r_k}}{\partial z}+e_{15}^k\frac{\partial {\phi }_k}{\partial r},\\ {\tau }_{\theta {\theta }_k}=c_{12}^k\frac{\partial u_{r_k}}{\partial r}+c_{11}^k\frac{u_{r_k}}{r}+c_{13}^k\frac{\partial u_{z_k}}{\partial z}+e_{31}^k\frac{\partial {\phi }_k}{\partial z},\\ {\tau }_{z{z}_k}=c_{13}^k\frac{\partial u_{r_k}}{\partial r}+c_{13}^k\frac{u_{r_k}}{r}+c_{33}^k\frac{\partial u_{z_k}}{\partial z}+e_{33}^k\frac{\partial {\phi }_k}{\partial z},\\ D_{r_k}=e_{15}^k\frac{\partial u_{z_k}}{\partial r}+e_{15}^k\frac{\partial u_{r_k}}{\partial z}-{\epsilon }_{11}^k\frac{\partial {\phi }_k}{\partial r},\\ D_{z_k}=e_{31}^k\frac{\partial u_{r_k}}{\partial r}+e_{31}^k\frac{u_{r_k}}{r}+e_{33}^k\frac{\partial u_{z_k}}{\partial z}-{\epsilon }_{33}^k\frac{\partial {\phi }_k}{\partial z}.\end{array}\end{array}$$

3. Solution to the Problem

Inserting Equation (2) into Equation (1), we have$$\begin{array}{}\text{(3)}& c_{11}^k\frac{{\partial }^2{u}_{r_k}}{\partial r^2}+\frac{1}{r}\frac{\partial u_{r_k}}{\partial r}+c_{44}^k\frac{{\partial }^2{u}_{r_k}}{\partial z^2}-c_{11}^k\frac{u_{r_k}}{r^2}+c_{13}^k+c_{44}^k\frac{{\partial }^2{u}_{z_k}}{\partial r\partial z}+e_{31}^k+e_{15}^k\frac{{\partial }^2{\phi }_k}{\partial r\partial z}={\rho }_k\frac{{\partial }^2{u}_{r_k}}{\partial t^2},\end{array}$$$$\begin{array}{}\text{(4)}& \begin{array}{c}{c}_{44}^k\frac{{\partial }^2{u}_{z_k}}{\partial r^2}+\frac{1}{r}\frac{\partial u_{z_k}}{\partial r}+c_{33}^k\frac{{\partial }^2{u}_{z_k}}{\partial z^2}+c_{13}^k+c_{44}^k\frac{{\partial }^2{u}_{r_k}}{\partial r\partial z}+\frac{1}{r}\frac{\partial u_{r_k}}{\partial z}\\ +e_{15}^k\frac{1}{r}\frac{\partial {\phi }_k}{\partial r}+e_{15}^k\frac{{\partial }^2{\phi }_k}{\partial r^2}+e_{33}^k\frac{{\partial }^2{\phi }_k}{\partial z^2}={\rho }_k\frac{{\partial }^2{u}_{z_k}}{\partial t^2},\end{array}\end{array}$$$$\begin{array}{}\text{(5)}& e_{15}^k\frac{{\partial }^2{u}_{z_k}}{\partial r^2}+\frac{1}{r}\frac{\partial u_{z_k}}{\partial r}+e_{15}^k+e_{31}^k\frac{{\partial }^2{u}_{r_k}}{\partial r\partial z}+\frac{1}{r}\frac{\partial u_{r_k}}{\partial z}-{\epsilon }_{11}^k\frac{{\partial }^2{\phi }_k}{\partial r^2}+\frac{1}{r}\frac{\partial {\phi }_k}{\partial r}+e_{33}^k\frac{{\partial }^2{u}_{z_k}}{\partial z^2}-{\epsilon }_{33}^k\frac{{\partial }^2{\phi }_k}{\partial z^2}=0.\end{array}$$

In order to further solve the problem, the piezoelectric bimaterials are in a mechanical and electrical static state when $t=0$. The initial conditions for piezoelectric bimaterials are zero at $t=0$, which can be expressed by the formulas as follows:$$\begin{array}{}\text{(6)}& {u_{r_k}r,z,0=0,\frac{\partial u_{r_k}r,z,t}{\partial t}}_{t=0}=0,u_{z_k}r,z,0=0,{\frac{\partial u_{z_k}r,z,t}{\partial t}}_{t=0}=0.\end{array}$$

For piezoelectric media, the boundary conditions of the crack surface are well known. In this paper, we consider the boundary conditions in the presence of a circular coin-type interface crack. For piezoelectric bimaterials, we assume that there are impact stresses and electric fields above and below, and it is sufficient to solve the corresponding problems presented in the following boundary conditions according to symmetry. The boundary conditions under impacts of elastic stress and electric field at $t>0$ are as follows:$$\begin{array}{}\text{(7)}& {\tau }_{z{z}_1}r,h_1,t={\tau }_{z{z}_2}r,-h_2,t={\tau }_{0}Ht,r<\mathrm{\infty },\end{array}$$$$\begin{array}{}\text{(8)}& E_{z_1}r,h_1,t=E_{z_2}r,-h_2,t=E_{0}Ht,r<\mathrm{\infty },\end{array}$$$$\begin{array}{}\text{(9)}& u_{r_1}r,h_1,t=u_{r_2}r,-h_2,t,r<\mathrm{\infty }.\end{array}$$

In this paper, the defect problem of piezoelectric bimaterials with finite thickness is studied. According to symmetry, the upper and lower materials are continuous in the $z=0$ plane (crack plane). On the other hand, applying symmetry to the $z=0$ plane at $t>0$, we get$$\begin{array}{}\text{(10)}& {\tau }_{r{z}_1}r,0^{+},t={\tau }_{r{z}_2}r,0^{-},t=0,r<a,\end{array}$$$$\begin{array}{}\text{(11)}& D_{z_1}r,0^{+},t=D_{z_2}r,0^{-},t=d_{0}t,r<a,\end{array}$$$$\begin{array}{}\text{(12)}& {\tau }_{z{z}_1}r,0^{+},t={\tau }_{z{z}_2}r,0^{-},t=0,r<a,\end{array}$$$$\begin{array}{}\text{(13)}& u_{z_1}r,0^{+},t=u_{z_2}r,0^{-},t=0,r\ge a,\end{array}$$$$\begin{array}{}\text{(14)}& {\phi }_{1}r,0^{+},t={\phi }_{2}r,0^{-},t=0,r\ge a,\end{array}$$$$\begin{array}{}\text{(15)}& {\tau }_{r{z}_1}r,0^{+},t={\tau }_{r{z}_2}r,0^{-},t=0,r\ge a.\end{array}$$

To solve, the following is the Laplace transform concerning the time $t$$$\begin{array}{}\text{(16)}& f^{\ast }x,y,p={\int }_0^{+\mathrm{\infty }}fx,y,t{e}^{-pt}dt,fx,y,t=\frac{1}{2\pi i}{\int }_{Br}{f}^{\ast }x,y,p{e}^{pt}dp,\end{array}$$where $Br$ denotes the Bromwich path of integration.

Using the Laplace transform, Equations (7)–(15) can be transformed into$$\begin{array}{}\text{(17)}& {\tau }_{z{z}_1}^{\ast }r,h_1,p={\tau }_{z{z}_2}^{\ast }r,-h_2,p=\frac{{\tau }_0}{p},r<\mathrm{\infty },\end{array}$$$$\begin{array}{}\text{(18)}& E_{z_1}^{\ast }r,h_1,p=E_{z_2}^{\ast }r,-h_2,p=\frac{E_0}{p},r<\mathrm{\infty },\end{array}$$$$\begin{array}{}\text{(19)}& u_{r_1}^{\ast }r,h_1,p=u_{r_2}^{\ast }r,-h_2,p,r<\mathrm{\infty },\end{array}$$$$\begin{array}{}\text{(20)}& {\tau }_{r{z}_1}^{\ast }r,0^{+},p={\tau }_{r{z}_2}^{\ast }r,0^{-},p=0,r<a,\end{array}$$$$\begin{array}{}\text{(21)}& D_{z_1}^{\ast }r,0^{+},p=D_{z_2}^{\ast }r,0^{-},p=d_0^{\ast }p,r<a,\end{array}$$$$\begin{array}{}\text{(22)}& {\tau }_{z{z}_1}^{\ast }r,0^{+},p={\tau }_{z{z}_2}^{\ast }r,0^{-},p=0,r<a,\end{array}$$$$\begin{array}{}\text{(23)}& u_{z_1}^{\ast }r,0^{+},p=u_{z_2}^{\ast }r,0^{-},p=0,r\ge a,\end{array}$$$$\begin{array}{}\text{(24)}& {\phi }_1^{\ast }r,0^{+},p={\phi }_2^{\ast }r,0^{-},p=0,r\ge a,\end{array}$$$$\begin{array}{}\text{(25)}& {\tau }_{r{z}_1}^{\ast }r,0^{+},p={\tau }_{r{z}_2}^{\ast }r,0^{-},p=0,r\ge a.\end{array}$$

From Equations (17)–(19), the elastic displacement and electric potential by unknown functions $A_{j_k}$ and $B_{j_k}$ ($j=1,2,3$ and $k=1,2$) in the transformed domain can be obtained as follows:$$\begin{array}{}\text{(26)}& u_{r_k}^{\ast }r,z,p=\sum _{j=1}^3{\int }_0^{\mathrm{\infty }}{A}_{j_k}\xi ,p{e}^{3-2k{\beta }_{j_k}\xi z}+B_{j_k}\xi ,p{e}^{-3-2k{\beta }_{j_k}\xi z}{J}_1\xi rd\xi ,\end{array}$$$$\begin{array}{}\text{(27)}& u_{z_k}^{\ast }r,z,p=-\sum _{j=1}^3{\int }_0^{\mathrm{\infty }}{\gamma }_{1j_k}{\beta }_{j_k}{A}_{j_k}\xi ,p{e}^{3-2k{\beta }_{j_k}\xi z}-B_{j_k}\xi ,p{e}^{-3-2k{\beta }_{j_k}\xi z}{J}_0\xi rd\xi +\frac{A_0^{k}z}{p},\end{array}$$$$\begin{array}{}\text{(28)}& {\phi }_k^{\ast }r,z,p=-\sum _{j=1}^3{\int }_0^{\mathrm{\infty }}{\gamma }_{2j_k}{\beta }_{j_k}{A}_{j_k}\xi ,p{e}^{3-2k{\beta }_{j_k}\xi z}-B_{j_k}\xi ,p{e}^{-3-2k{\beta }_{j_k}\xi z}{J}_0\xi rd\xi -\frac{B_{0}z}{p},\end{array}$$where $A_0^k=\frac{{\tau }_0+e_{33}^k{E}_0}{c_{33}^k}$, $B_0=E_0$, $J_0\xi r$, and $J_1\xi r$ are zeroth order and first order Bessel functions of the first kind, respectively, known constants ${\gamma }_{1j_k},{\gamma }_{2j_k}$ and ${\beta }_{j_k}$ are determined by the Hankel and Laplace transform.

By substituting Equations (26)–(28) to Equations (3)–(5), we find ${\beta }_{j_k}$ satisfy the following$$\begin{array}{}\text{(29)}& \text{det}{R}_k=0,\end{array}$$where$$\begin{array}{}\text{(30)}& R_k=\begin{array}{ccc}{c}_{11}^k-c_{44}^k{\beta }_{j_k}^2{3-2k}^2+{\rho }^k\frac{p^2}{{\xi }^2}& c_{13}^k+c_{44}^k{\beta }_{jk}3-2k& e_{15}^k+e_{31}^k{\beta }_{jk}3-2k\\ c_{13}^k+c_{44}^k{\beta }_{jk}3-2k& c_{33}^k{\beta }_{j_k}^2{3-2k}^2-c_{44}^k-{\rho }^k\frac{p^2}{{\xi }^2}& e_{33}^k{\beta }_{j_k}^2{3-2k}^2-e_{15}^k\\ e_{15}^k+e_{31}^k{\beta }_{jk}3-2k& e_{33}^k{\beta }_{j_k}^2{3-2k}^2-e_{15}^k& {\epsilon }_{11}^k-{\epsilon }_{33}^k{\beta }_{j_k}^2{3-2k}^2\end{array},\end{array}$$and ${\gamma }_{i{j}_k}i=1,2$ meet this condition$$\begin{array}{}\text{(31)}& R_k\begin{array}{c}1\\ -{\gamma }_{1j_k}{\beta }_{j_k}\\ -{\gamma }_{2j_k}{\beta }_{j_k}\end{array}=0.\end{array}$$

From Equation (2) and Equations (26)–(28), the solutions can be obtained as follows:$$\begin{array}{}\text{(32)}& \begin{array}{c}{\tau }_{z{z}_k}^{\ast }r,z,p=c_{13}^k\frac{\partial u_{r_k}^{\ast }}{\partial r}+c_{13}^k\frac{u_{r_k}^{\ast }}{r}+c_{33}^k\frac{\partial u_{z_k}^{\ast }}{\partial z}+e_{33}^k\frac{\partial {\phi }_k^{\ast }}{\partial z}\\ =-\sum _{j=1}^3{\int }_0^{\mathrm{\infty }}{c}_{33}^k{\gamma }_{1j_k}+e_{33}^k{\gamma }_{2j_k}3-2k{\beta }_{j_k}^2-c_{13}^k\xi A_{j_k}\xi ,p{e}^{3-2k{\beta }_{j_k}\xi z}+B_{j_k}\xi ,p{e}^{-3-2k{\beta }_{j_k}\xi z}{J}_0\xi rd\xi +\frac{{\tau }_0}{p}\\ =-\sum _{j=1}^3{\int }_0^{\mathrm{\infty }}{\eta }_{1j_k}\xi A_{j_k}\xi ,p{e}^{3-2k{\beta }_{j_k}\xi z}+B_{j_k}\xi ,p{e}^{-3-2k{\beta }_{j_k}\xi z}{J}_0\xi rd\xi +\frac{{\tau }_0}{p},\end{array}\end{array}$$$$\begin{array}{}\text{(33)}& \begin{array}{c}{\tau }_{r{z}_k}^{\ast }r,z,p=c_{44}^k\frac{\partial u_{r_k}^{\ast }}{\partial z}+c_{44}^k\frac{\partial u_{z_k}^{\ast }}{\partial r}+e_{15}^k\frac{\partial {\phi }_k^{\ast }}{\partial r}\\ =\sum _{j=1}^3{\int }_0^{\mathrm{\infty }}{c}_{44}^{k}3-2k+{\gamma }_{1j_k}+e_{15}^k{\gamma }_{2j_k}{\beta }_{j_k}\xi A_{j_k}\xi ,p{e}^{3-2k{\beta }_{j_k}\xi z}-B_{j_k}\xi ,p{e}^{-3-2k{\beta }_{j_k}\xi z}{J}_1\xi rd\xi \\ =\sum _{j=1}^3{\int }_0^{\mathrm{\infty }}{\eta }_{2j_k}{A}_{j_k}\xi ,p{e}^{3-2k{\beta }_{j_k}\xi z}-B_{j_k}\xi ,p{e}^{-3-2k{\beta }_{j_k}\xi z}{J}_1\xi rd\xi ,\end{array}\end{array}$$$$\begin{array}{}\text{(34)}& \begin{array}{cc}{D}_{z_k}^{\ast }r,z,p& =e_{31}^k\frac{\partial u_{r_k}^{\ast }}{\partial r}+e_{31}^k\frac{u_{r_k}^{\ast }}{r}+e_{33}^k\frac{\partial u_{z_k}^{\ast }}{\partial z}-{\epsilon }_{33}^k\frac{\partial {\phi }_k^{\ast }}{\partial z}\hfill \\ & =-\sum _{j=1}^3{\int }_0^{\mathrm{\infty }}{\eta }_{3j_k}\xi A_{j_k}\xi ,p{e}^{3-2k{\beta }_{j_k}\xi z}+B_{j_k}\xi ,p{e}^{-3-2k{\beta }_{j_k}\xi z}{J}_0\xi rd\xi +\frac{{\bar{D}}_0^k}{p},\hfill \end{array}\end{array}$$where$$\begin{array}{}\text{(35)}& {\bar{D}}_0^{\mathrm{k}}=\frac{e_{33}^k}{c_{33}^k}{\tau }_0+e_{33}^k{E}_0+{\epsilon }_{33}^k{E}_0,{\eta }_{1j_k}=c_{33}^k{\gamma }_{1j_k}+e_{33}^k{\gamma }_{2j_k}3-2k{\beta }_{j_k}^2-c_{13}^k,\end{array}$$$$\begin{array}{}\text{(36)}& {\eta }_{2j_k}=c_{44}^{k}3-2k+{\gamma }_{1j_k}+e_{15}^k{\gamma }_{2j_k}{\beta }_{j_k}\xi ,{\eta }_{3j_k}=e_{33}^k{\gamma }_{1j_k}-{\epsilon }_{33}^k{\gamma }_{2j_k}3-2k{\beta }_{j_k}^2-e_{31}^k.\end{array}$$

According to Li and Lee [42], the relationship between ${\tau }_0,E_0$ and $D_0$ is as follows:$$\begin{array}{}\text{(37)}& D_0=\frac{c_{11}+c_{12}{e}_{33}-2c_{13}{e}_{31}}{c_{11}+c_{12}{c}_{33}-2c_{13}^2}{\tau }_0+\frac{c_{11}+c_{12}{e}_{33}^2+2c_{33}{e}_{31}^2-4c_{13}{e}_{31}{e}_{33}}{c_{11}+c_{12}{c}_{33}-2c_{13}^2}+{\epsilon }_{33}{E}_0.\end{array}$$

According to the boundary conditions, we have$$\begin{array}{}\text{(38)}& A_{j_k}=-e^{-2{\beta }_{j_k}\xi h_k}{B}_{j_k},B_{j_1}=-\frac{1+e^{-2{\beta }_{j_2}\xi h_2}}{1+e^{-2{\beta }_{j_1}\xi h_1}}{B}_{j_2}.\end{array}$$

Inserting boundary conditions Equations (23) and (24) into Equations (27) and (28), respectively, they are$$\begin{array}{}\text{(39)}& -\sum _{j=1}^3{\int }_0^{\mathrm{\infty }}{\gamma }_{1j_k}{\beta }_{j_k}{A}_{j_k}\xi ,p-B_{j_k}\xi ,p{J}_0\xi rd\xi =0,\end{array}$$$$\begin{array}{}\text{(40)}& -\sum _{j=1}^3{\int }_0^{\mathrm{\infty }}{\gamma }_{2j_k}{\beta }_{j_k}{A}_{j_k}\xi ,p-B_{j_k}\xi ,p{J}_0\xi rd\xi =0,\end{array}$$$$\begin{array}{}\text{(41)}& \sum _{j=1}^3{\int }_0^{\mathrm{\infty }}{\eta }_{2j_k}{A}_{j_k}\xi ,p-B_{j_k}\xi ,p{J}_1\xi rd\xi =0.\end{array}$$

Introducing new functions $f_{i_k}s,pi=1,2$ and $q_{i_k}s,pi=1,2$, we have$$\begin{array}{}\text{(42)}& \sum _{j=1}^3{\gamma }_{1j_1}{\beta }_{j_1}{e}^{-2{\beta }_{j_1}\xi h_1}+1B_{j_1}-\sum _{j=1}^3{\gamma }_{1j_2}{\beta }_{j_2}{e}^{-2{\beta }_{j_2}\xi h_2}+1B_{j_2}=-{\int }_0^a{f}_{1}s,p\sin \xi sds=-q_1,\end{array}$$$$\begin{array}{}\text{(43)}& \sum _{j=1}^3{\gamma }_{2j_1}{\beta }_{j_1}{e}^{-2{\beta }_{j_1}\xi h_1}+1B_{j_1}-\sum _{j=1}^3{\gamma }_{2j_2}{\beta }_{j_2}{e}^{-2{\beta }_{j_2}\xi h_2}+1B_{j_2}=-{\int }_0^a{f}_{2}s,p\sin \xi sds=-q_2,\end{array}$$$$\begin{array}{}\text{(44)}& \sum _{j=1}^3{\eta }_{2j_2}{e}^{-2{\beta }_{j_2}\xi h_2}+1B_{j_2}=0.\end{array}$$

From Equations (42)–(44), one can see that the unknown functions $B_{j_2}\xi ,pj=1,2,3$ can be obtained as follows:$$\begin{array}{}\text{(45)}& \begin{array}{c}{B}_{1_2}\xi ,p\\ B_{2_2}\xi ,p\\ B_{3_2}\xi ,p\end{array}=-{C_{ji}}_{3\times 3}{\int }_0^a\begin{array}{c}{f}_1\\ f_2\\ 0\end{array}\sin \xi sds.\end{array}$$

From the study of Abramowitz and Stegun [43], we have$$\begin{array}{}\text{(46)}& {\int }_0^{\mathrm{\infty }}{J}_0\xi r\sin \xi sd\xi =\begin{array}{cc}0,& r<s\\ \frac{Hs-r}{\sqrt{s^2-r^2}},& r>s\end{array}.\end{array}$$

By applying Equation (46) and substituting Equations (27) and (28) into Equations (23) and (24), respectively, the values of elastic displacement $u_{z_k}^{\ast }$ and electric potential ${\phi }_k^{\ast }$ on the $z=0$ can be obtained as follows:$$\begin{array}{}\text{(47)}& u_{z_k}^{\ast }r,0,p=-{\int }_0^{\mathrm{\infty }}{\int }_r^a{f}_{1}s,p\sin \xi s{J}_0\xi rdsd\xi =-{\int }_r^a{f}_{1}s,p\frac{1}{\sqrt{s^2-r^2}}ds,\end{array}$$$$\begin{array}{}\text{(48)}& {\phi }_k^{\ast }r,0,p=-{\int }_0^{\mathrm{\infty }}{\int }_r^a{f}_{2}s,p\sin \xi s{J}_0\xi rdsd\xi =-{\int }_r^a{f}_{2}s,p\frac{1}{\sqrt{s^2-r^2}}ds.\end{array}$$

It is moreover, inserting Equations (32) and (34) into Equations (21) and (22), respectively, we arrive at$$\begin{array}{}\text{(49)}& {\tau }_{z{z}_2}^{\ast }r,0,p=-{\int }_0^{\mathrm{\infty }}\sum _{i=1}^2\sum _{j=1}^3{\eta }_{{}_{1j_2}}{C}_{ji}1-e^{-2{\beta }_{j_2}\xi h_2}{q}_i\xi J_0\xi rd\xi +\frac{{\tau }_0}{p},\end{array}$$$$\begin{array}{}\text{(50)}& D_{z_2}^{\ast }r,0,p=-{\int }_0^{\mathrm{\infty }}\sum _{i=1}^2\sum _{j=1}^3{\eta }_{{}_{3j_2}}{C}_{ji}1-e^{-2{\beta }_{j_2}\xi h_2}{q}_i\xi J_0\xi rd\xi +\frac{{\bar{D}}_0}{p}.\end{array}$$

Then we have$$\begin{array}{}\text{(51)}& {\int }_0^{\mathrm{\infty }}\sum _{i=1}^2\sum _{j=1}^3{\eta }_{{}_{1j_2}}{C}_{ji}1-e^{-2{\beta }_{j_2}\xi h_2}{q}_i\xi J_0\xi rd\xi =\frac{{\tau }_0}{p},\end{array}$$$$\begin{array}{}\text{(52)}& {\int }_0^{\mathrm{\infty }}\sum _{i=1}^2\sum _{j=1}^3{\eta }_{{}_{3j_2}}{C}_{ji}1-e^{-2{\beta }_{j_2}\xi h_2}{q}_i\xi J_0\xi rd\xi =\frac{{\bar{D}}_0}{p}.\end{array}$$

Using what we already know [42]$$\begin{array}{}\text{(53)}& {\int }_0^x\frac{r{J}_0\xi r}{\sqrt{x^2-r^2}}dr=\frac{\sin \xi x}{\xi },\xi J_0\xi r=\sin \xi x.\end{array}$$

We finally obtain$$\begin{array}{}\text{(54)}& {\int }_0^{\mathrm{\infty }}\sum _{i=1}^2\sum _{j=1}^3{\eta }_{{}_{1j_2}}{C}_{ji}1-e^{-2{\beta }_{j_2}\xi h_2}{q}_i\sin \xi xd\xi =\frac{{\tau }_0}{p},\end{array}$$$$\begin{array}{}\text{(55)}& {\int }_0^{\mathrm{\infty }}\sum _{i=1}^2\sum _{j=1}^3{\eta }_{{}_{3j2}}{C}_{ji}1-e^{-2{\beta }_{j2}\xi h_2}{q}_i\sin \xi xd\xi =\frac{{\bar{D}}_0}{p}.\end{array}$$

By substituting Equations (43) and (44) into Equations (54) and (55), respectively, one can obtain the following results$$\begin{array}{}\text{(56)}& {\int }_0^{\mathrm{\infty }}{\int }_0^a\sum _{i=1}^2\sum _{j=1}^3{\eta }_{{}_{1j_2}}{C}_{ji}1-e^{-2{\beta }_{j_2}\xi h_2}{f}_{i}s,p\sin \xi s\sin \xi xd\xi ds=\frac{{\tau }_0}{p},\end{array}$$$$\begin{array}{}\text{(57)}& {\int }_0^{\mathrm{\infty }}{\int }_0^a\sum _{i=1}^2\sum _{j=1}^3{\eta }_{{}_{3j_2}}{C}_{ji}1-e^{-2{\beta }_{j_2}\xi h_2}{f}_{i}s,p\sin \xi s\sin \xi xd\xi ds=\frac{{\bar{D}}_0}{p}.\end{array}$$

Using the known relation [42]$$\begin{array}{}\text{(58)}& {\int }_0^{\mathrm{\infty }}\sin \xi s\sin \xi xd\xi =\frac{\pi }{2}\delta s-x,\end{array}$$where $\delta \ast$ is the Dirac delta function.

By applying Equation (58), Equations (56) and (57) can be further read as follows:$$\begin{array}{}\text{(59)}& {\int }_0^a\sum _{i=1}^2{d}_{1i}{f}_{i}s,p\frac{\pi }{2}\delta s-xds=\frac{{\tau }_0}{p},\end{array}$$$$\begin{array}{}\text{(60)}& {\int }_0^a\sum _{i=1}^2{d}_{2i}{f}_{i}s,p\frac{\pi }{2}\delta s-xds=\frac{{\bar{D}}_0}{p},\end{array}$$where$$\begin{array}{}\text{(61)}& \sum _{j=1}^3{\eta }_{{}_{1j_2}}{C}_{ji}1-e^{-2{\beta }_{j_2}\xi h_2}=d_{1i},\end{array}$$$$\begin{array}{}\text{(62)}& \sum _{j=1}^3{\eta }_{{}_{3j_2}}{C}_{ji}1-e^{-2{\beta }_{j_2}\xi h_2}=d_{2i}.\end{array}$$

Let us consider the following limiting form$$\begin{array}{}\text{(63)}& \underset{\xi \to \mathrm{\infty }}{\text{lim}}{d}_{1i}\xi ,p=l_{1i},i=1,2,\end{array}$$$$\begin{array}{}\text{(64)}& \underset{\xi \to \mathrm{\infty }}{\text{lim}}{d}_{2i}\xi ,p=l_{2i},i=1,2.\end{array}$$

Equations (59) and (60) can be further written as follows:$$\begin{array}{}\text{(65)}& \sum _{i=1}^2{l}_{1i}{f}_{i}x,p+{\int }_0^a\sum _{i=1}^2{L}_{1i}s,x,p{f}_{i}s,pds=\frac{2{\tau }_0}{p\pi }x,x<a,\end{array}$$$$\begin{array}{}\text{(66)}& \sum _{i=1}^2{l}_{2i}{f}_{i}x,p+{\int }_0^a\sum _{i=1}^2{L}_{2i}s,x,p{f}_{i}s,pds=\frac{2{\bar{D}}_0}{p\pi }x,x<a,\end{array}$$where$$\begin{array}{}\text{(67)}& L_{1i}s,x,p=\frac{2}{\pi }{\int }_0^{\mathrm{\infty }}{l}_{1i}-d_{1i}\sin \xi t\sin \xi xd\xi ,i=1,2,\end{array}$$$$\begin{array}{}\text{(68)}& L_{2i}s,x,p=\frac{2}{\pi }{\int }_0^{\mathrm{\infty }}{l}_{2i}-d_{2i}\sin \xi t\sin \xi xd\xi ,i=1,2.\end{array}$$

For the convenience of calculation, the following transformations are introduced as follows:$$\begin{array}{}\text{(69)}& \stackrel{⌢}{s}=\frac{s}{a},\stackrel{⌢}{x}=\frac{x}{a},g_i\stackrel{⌢}{x},p=\frac{p{f}_{i}x,p}{a}.\end{array}$$

Hence, Equations (65) and (66) can be expressed as follows:$$\begin{array}{}\text{(70)}& \sum _{i=1}^2{l}_{1i}{f}_i\stackrel{⌢}{x},p+{\int }_0^1\sum _{i=1}^2{\stackrel{⌢}{L}}_{1i}\stackrel{⌢}{s},\stackrel{⌢}{x},p{f}_i\stackrel{⌢}{s},pd\stackrel{⌢}{s}=\frac{2{\tau }_0}{\pi }\stackrel{⌢}{x},\end{array}$$$$\begin{array}{}\text{(71)}& \sum _{i=1}^2{l}_{2i}{f}_i\stackrel{⌢}{x},p+{\int }_0^1\sum _{i=1}^2{L}_{2i}\stackrel{⌢}{s},\stackrel{⌢}{x},p{f}_i\stackrel{⌢}{s},pd\stackrel{⌢}{s}=\frac{2{\bar{D}}_0}{\pi }\stackrel{⌢}{x}.\end{array}$$

4. Dynamic Field Intensity Factors

From the perspective of fracture mechanics, the dynamic field intensity factor is an essential factor in characterize the penny-shaped crack. From Equation (32), we have$$\begin{array}{}\text{(72)}& {\tau }_{zz}^{\ast }r,p=-{\int }_0^{\mathrm{\infty }}{\int }_0^a\sum _{i=1}^2\sum _{j=1}^3{\eta }_{{}_{1j2}}{C}_{ji}1-e^{-2{\beta }_{j2}\xi h_2}{f}_{i}s,p\xi \sin \xi s{J}_{0}r\xi d\xi ds,\end{array}$$$$\begin{array}{}\text{(73)}& D_z^{\ast }r,p=-{\int }_0^{\mathrm{\infty }}{\int }_0^a\sum _{i=1}^2\sum _{j=1}^3{\eta }_{{}_{3j2}}{C}_{j{i}_2}1-e^{-2{\beta }_{j_2}\xi h_2}{f}_{i}s,p\xi \sin \xi s{J}_{0}r\xi d\xi ds.\end{array}$$

After some calculation, the expressions for ${\tau }_{zz}^{\ast }r,p$ and $D_z^{\ast }r,p$ can be obtained as follows:$$\begin{array}{}\text{(74)}& {\tau }_{zz}^{\ast }r,p=\sum _{i=1}^2{l}_{1i}{f}_{i}a,p{\int }_0^{\mathrm{\infty }}\cos \xi a{J}_0\xi rd\xi +o1,\end{array}$$$$\begin{array}{}\text{(75)}& D_z^{\ast }r,p=\sum _{i=1}^2{l}_{2i}{f}_{i}a,p{\int }_0^{\mathrm{\infty }}\cos \xi a{J}_0\xi rd\xi +o1.\end{array}$$

We note that$$\begin{array}{}\text{(76)}& {\int }_0^{\mathrm{\infty }}\cos \xi a{J}_0\xi rd\xi =\frac{1}{\sqrt{r^2-a^2}},r>a.\end{array}$$

The dynamic field intensity factor in the Laplace transform domain can be defined as follows:$$\begin{array}{}\text{(77)}& K_{\tau }^{\ast }p=\underset{r\to a^{+}}{\text{lim}}\sqrt{2\pi r-a}{\tau }_{zz}^{\ast }r,p=\frac{\sum _{i=1}^2{l}_{1i}{f}_{i}1,p}{p}\sqrt{\pi a},\end{array}$$$$\begin{array}{}\text{(78)}& K_D^{\ast }p=\underset{r\to a^{+}}{\text{lim}}\sqrt{2\pi r-a}{D}_z^{\ast }r,p=\frac{\sum _{i=1}^2{l}_{2i}{f}_{i}1,p}{p}\sqrt{\pi a}.\end{array}$$

Similarly, the field intensity factor is related to the propagation displacement of the crack tip and the electric potential, so it can be defined as follows:$$\begin{array}{}\text{(79)}& K_{\text{COD}}^{\ast }p=\underset{r\to a^{-}}{\text{lim}}\sqrt{\frac{\pi }{2a-r}}{u}_z^{\ast }r,p=\frac{f_{1}a,p}{a}\sqrt{\pi a},\end{array}$$$$\begin{array}{}\text{(80)}& K_{\phi }^{\ast }p=\underset{r\to a^{-}}{\text{lim}}\sqrt{\frac{\pi }{2a-r}}{\phi }^{\ast }r,p=\frac{f_{2}a,p}{a}\sqrt{\pi a}.\end{array}$$

According to Equations (69)–(71), one can obtain$$\begin{array}{}\text{(81)}& K_{\text{COD}}=\frac{f_{1}1,p}{a}\sqrt{\pi a},\end{array}$$$$\begin{array}{}\text{(82)}& K_{\phi }=\frac{f_{2}1,p}{a}\sqrt{\pi a}.\end{array}$$

5. Numerical Examples

When the dynamic field intensity factors exceed the corresponding critical value, the crack expands. Some numerical calculations are shown to describe the dynamic behavior of a penny-shaped crack. The material coefficients of the numerical examples are shown in Table 1.

Table 1

The properties of piezoelectric material.

Material	PZT-6B	PZT-5H	BaTiO3
$c_{11}$ (1010 N/m2)	16.8	12.6	22.6
$c_{33}$ (1010 N/m2)	16.3	11.7	21.6
$c_{44}$ (1010 N/m2)	2.71	3.53	4.4
$c_{12}$ (1010 N/m2)	6	5.3	12.5
$c_{13}$ (1010 N/m2)	6	5.5	12.4
$e_{31}$ (C/m2)	−0.9	−6.5	−2.2
$e_{33}$ (C/m2)	7.1	23.3	9.3
$e_{15}$ (C/m2)	4.6	17	5.8
${\epsilon }_{11}$ (10−10 F/m)	36	15.052	56.4
${\epsilon }_{33}$ (10−10 F/m)	34	13	63.5
$\rho$ (103 kg/m3)	7.55	7.5	5.7

Let us take the inverse Laplace transform as Stehfest [44] proposed. In our example, the stress impact load ${\tau }_0$ is selected to be 4.2 MPa. The parameters relating the electrical and mechanical loadings are expressed by $I_{\ast }=e_{33}{E}_0/{\tau }_0$ and $c_v=vt$. Let us define the dimensional function as $F/F_0$, where $F_0$ is the static COD intensity factor ($F$ is $K_{\text{COD}}$ in Equation (81)).

When the double material degenerates into a single material, we can calculate the corresponding value (Appendix B). Figure 2 is drawn to describe the trend of $F/F_0$ relative to $h/a$ ($E_{0_{\ast }}=1$). The two piezoelectric single materials are PZT-6B and BaTiO₃, respectively. The $F/F_0$ of both materials decreases monotonically and tends to change smoothly with increasing $h/a$. On the other hand the changes in the two curves are shown in Figure 2. We can see that the change in the value of the material BaTiO₃ is greater than that of material PZT-6B, which in a particular sense implies the correctness of our results.

[figure(s) omitted; refer to PDF]

For the $F/F_0$ plot with $E_0$ change in Figure 3, a thinner piezoelectric layer can be found to cause an increase in $F/F_0$, which means that the contribution of the electric field to the fracture toughness is evident for the thin piezoelectric layer in Figure 3. When $h_1/a>2,F/F_0$ is also insensitive.

[figure(s) omitted; refer to PDF]

Figure 4(a)–4(c) discuss the influence of the applied electric fields $E_0$ and $vt/a$ on $F/F_0$. Under the impact problem, this time normalization is denoted by $vt/a$, where $v$ is the shear wave velocity. In Figure 4, for a given pressure ${\tau }_0,F/F_0$ increases with increasing of electric fields $E_0$ and $vt/a$. At a fixed time, it can be found that an increase in $E_0$ will lead to an increase in $F/F_0$, which is consistent with the result shown in Figure 3. For a fixed electric field $E_0$, an increase in $vt/a$ leads to an increase in $F/F_0$. However, we can see that the increase in $F/F_0$ flattens out in the later period. From Figures 2–4 show that different materials play different roles in the crack propagation.

[figure(s) omitted; refer to PDF]

Figure 5 shows that when the substrate is constant, the thicker the material 1, the less influence it has. In other words, when the thickness of Material 2 is constant, the thickness of Material 1 is small, and the COD intensity factor changes greatly, indicating that the thinner of Material 1 when appropriate, the more conducive it is to safe design.

[figure(s) omitted; refer to PDF]

It can be seen from Figures 3 and 4 that both positive and negative electric fields will affect the propagation of cracks. At the same time, it can be found that the COD intensity factor value is large under a positive electric field, while the COD intensity factor value is small under a negative electric field, which indicates that the positive electric field will cause the expansion of the piezoelectric material, while the negative electric field will cause the contraction of the piezoelectric material. This conclusion is consistent with that of the study by Li and Lee [42]. In other words, crack propagation can be controlled by the regulation of electric field.

For the $F/F_0$ plot with $I_{\ast }$ in Figure 6, a thinner piezoelectric layer can be found to cause an increase in $F/F_0$. In addition, observations of Figures 4–6 suggest that the electromechanical coupling coefficient plays an important role in crack propagation: the stronger the electromechanical coupling effect, the smaller $F/F_0$ is.

[figure(s) omitted; refer to PDF]

All the c figures in Figures 4–6 studied the crack growth of different piezoelectric bimaterials. By comparing the images, we can find that when Materials 1 and 2 are different, the variation trend of COD intensity of crack growth is basically the same as that when Materials 1 and 2 are the same. That is, different material parameters will affect the change of COD intensity factor, so in practical applications, the appropriate material should be selected according to the actual situation.

It can be seen from Figures 5 and 6 that when the thickness of Material 2 is constant, the thickness of Material 1 is small, and the COD intensity factor changes greatly, indicating that the thinner the material when appropriate, the more conducive it is to safe design.

6. Conclusion

Based on the piezoelectric theory, the dynamic penny-shaped interface crack propagation of piezoelectric bimaterials is analyzed. The boundary conditions are transformed into a nonlinear Fredholm integral equation by using the Hankel transformation technique. Numerical solutions are given. According to the given surface displacement and stress of the layer, the corresponding models are constructed, and the displacement functions that meet the conditions are established. Solving the model, the influences of the electric field, impact time and layer thickness on the dynamic COD are analyzed accordingly. The results show that the COD increases with decreasing $h_1$. At the same time, when the speed $v$ and crack size $a$ are fixed, the COD increases with increasing impact time, and tends to be flat after reaching the peak value. The stronger the electromechanical coupling effect is, the smaller $F/F_0$ is. At the same time, different materials have different roles in crack propagation, so it is very important to study the COD factor of cracks for safety design.

Acknowledgments

This study was supported by the National Natural Science Foundation of China (11972019 and 52275567), Shanxi Provincial Key Research and Development Project (202102090301027), and Shanxi Postgraduate Innovation Project (2022Y668).

Appendix

A. Elements of Matrix [C]

$$\begin{array}{}\text{(A.1)}& A_1={\gamma }_{{11}_1}{\beta }_{1_1}1+e^{-2{\beta }_{1_1}\xi h_1}1-e^{-2{\beta }_{1_2}\xi h_2}/1-e^{-2{\beta }_{1_1}\xi h_1}-{\gamma }_{{11}_2}{\beta }_{1_2}1+e^{-2{\beta }_{1_2}\xi h_2},\end{array}$$$$\begin{array}{}\text{(A.2)}& B_1={\gamma }_{{12}_1}{\beta }_{2_1}1+e^{-2{\beta }_{2_1}\xi h_1}1-e^{-2{\beta }_{2_2}\xi h_2}/1-e^{-2{\beta }_{2_1}\xi h_1}-{\gamma }_{{12}_2}{\beta }_{2_2}1+e^{-2{\beta }_{2_2}\xi h_2},\end{array}$$$$\begin{array}{}\text{(A.3)}& C_1={\gamma }_{{13}_1}{\beta }_{3_1}1+e^{-2{\beta }_{3_1}\xi h_1}1-e^{-2{\beta }_{3_2}\xi h_2}/1-e^{-2{\beta }_{3_1}\xi h_1}-{\gamma }_{{13}_2}{\beta }_{3_2}1+e^{-2{\beta }_{3_2}\xi h_2},\end{array}$$$$\begin{array}{}\text{(A.4)}& A_2={\gamma }_{{21}_1}{\beta }_{1_1}1+e^{-2{\beta }_{1_1}\xi h_1}1-e^{-2{\beta }_{1_2}\xi h_2}/1-e^{-2{\beta }_{1_1}\xi h_1}-{\gamma }_{{21}_2}{\beta }_{1_2}1+e^{-2{\beta }_{1_2}\xi h_2},\end{array}$$$$\begin{array}{}\text{(A.5)}& B_2={\gamma }_{{22}_1}{\beta }_{2_1}1+e^{-2{\beta }_{2_1}\xi h_1}1-e^{-2{\beta }_{2_2}\xi h_2}/1-e^{-2{\beta }_{2_1}\xi h_1}-{\gamma }_{{22}_2}{\beta }_{2_2}1+e^{-2{\beta }_{2_2}\xi h_2},\end{array}$$$$\begin{array}{}\text{(A.6)}& C_2={\gamma }_{{23}_1}{\beta }_{3_1}1+e^{-2{\beta }_{3_1}\xi h_1}1-e^{-2{\beta }_{3_2}\xi h_2}/1-e^{-2{\beta }_{3_1}\xi h_1}-{\gamma }_{{23}_2}{\beta }_{3_2}1+e^{-2{\beta }_{3_2}\xi h_2},\end{array}$$$$\begin{array}{}\text{(A.7)}& A_3={\eta }_{{21}_2}1+e^{-2{\beta }_{1_2}\xi h_2},\end{array}$$$$\begin{array}{}\text{(A.8)}& B_3={\eta }_{{22}_2}1+e^{-2{\beta }_{2_2}\xi h_2},\end{array}$$$$\begin{array}{}\text{(A.9)}& C_3={\eta }_{{23}_2}1+e^{-2{\beta }_{3_2}\xi h_2},\end{array}$$$$\begin{array}{}\text{(A.10)}& C_0=A_1{B}_2{C}_3-C_2{B}_3-A_2{B}_1{C}_3-C_1{B}_3+A_3{B}_1{C}_2-B_2{C}_1,\end{array}$$$$\begin{array}{}\text{(A.11)}& C_{11}=B_2{C}_3-B_3{C}_2/C_0,\end{array}$$$$\begin{array}{}\text{(A.12)}& C_{12}=B_3{C}_1-B_1{C}_3/C_0,\end{array}$$$$\begin{array}{}\text{(A.13)}& C_{13}=B_1{C}_2-B_2{C}_1/C_0,\end{array}$$$$\begin{array}{}\text{(A.14)}& C_{21}=C_2{A}_3-C_3{A}_2/C_0,\end{array}$$$$\begin{array}{}\text{(A.15)}& C_{22}=C_3{A}_1-C_1{A}_3/C_0,\end{array}$$$$\begin{array}{}\text{(A.16)}& C_{23}=C_1{A}_2-C_2{A}_1/C_0,\end{array}$$$$\begin{array}{}\text{(A.17)}& C_{31}=A_2{B}_3-A_3{B}_2/C_0,\end{array}$$$$\begin{array}{}\text{(A.18)}& C_{32}=A_3{B}_1-A_1{B}_3/C_0,\end{array}$$$$\begin{array}{}\text{(A.19)}& C_{33}=A_1{B}_2-A_2{B}_1/C_0.\end{array}$$

B. The Element Value of Matrix [C] in Single Material

$$\begin{array}{}\text{(B.1)}& a=\frac{{\eta }_{23}{\beta }_{1}1+e^{2{\beta }_1\xi h}{\gamma }_{11}-{\eta }_{21}{\beta }_{3}1+e^{2{\beta }_3\xi h}{\gamma }_{13}{\beta }_{2}1+e^{2{\beta }_2\xi h}}{{\eta }_{22}{\beta }_{1}1+e^{2{\beta }_1\xi h}{\gamma }_{11}-{\eta }_{21}{\beta }_{2}1+e^{2{\beta }_2\xi h}{\gamma }_{12}{\beta }_{3}1+e^{2{\beta }_3\xi h}}-\frac{{\gamma }_{23}{\gamma }_{11}-{\gamma }_{21}{\gamma }_{13}}{{\gamma }_{22}{\gamma }_{11}-{\gamma }_{21}{\gamma }_{12}},\end{array}$$$$\begin{array}{}\text{(B.2)}& b_1=\frac{{\eta }_{21}{\beta }_{2}1+e^{2{\beta }_2\xi h}{\gamma }_{23}{\gamma }_{11}-{\gamma }_{21}{\gamma }_{13}}{{\eta }_{21}{\beta }_{2}1+e^{2{\beta }_2\xi h}{\gamma }_{12}-{\eta }_{22}{\beta }_{1}1+e^{2{\beta }_1\xi h}{\gamma }_{11}{\gamma }_{22}{\gamma }_{11}-{\gamma }_{21}{\gamma }_{12}},\end{array}$$$$\begin{array}{}\text{(B.3)}& b_2=\frac{{\eta }_{23}{\beta }_{1}1+e^{2{\beta }_1\xi h}{\gamma }_{11}-{\eta }_{21}{\beta }_{3}1+e^{2{\beta }_3\xi h}{\gamma }_{13}{\beta }_{2}1+e^{2{\beta }_2\xi h}{\gamma }_{21}}{{\eta }_{22}{\beta }_{1}1+e^{2{\beta }_1\xi h}{\gamma }_{11}-{\eta }_{21}{\beta }_{2}1+e^{2{\beta }_2\xi h}{\gamma }_{12}{\beta }_{3}1+e^{2{\beta }_3\xi h}{\gamma }_{21}{\gamma }_{12}-{\gamma }_{22}{\gamma }_{11}},\end{array}$$$$\begin{array}{}\text{(B.4)}& b_3=-\frac{{\eta }_{21}{\beta }_{2}1+e^{2{\beta }_2\xi h}}{{\eta }_{21}{\beta }_{2}1+e^{2{\beta }_2\xi h}{\gamma }_{12}-{\eta }_{22}{\beta }_{1}1+e^{2{\beta }_1\xi h}{\gamma }_{11}}+\frac{{\gamma }_{21}}{{\gamma }_{22}{\gamma }_{11}-{\gamma }_{21}{\gamma }_{12}},\end{array}$$$$\begin{array}{}\text{(B.5)}& b_4=\frac{{\eta }_{23}{\beta }_{1}1+e^{2{\beta }_1\xi h}{\gamma }_{11}-{\eta }_{21}{\beta }_{3}1+e^{2{\beta }_3\xi h}{\gamma }_{13}{\beta }_{2}1+e^{2{\beta }_2\xi h}{\gamma }_{11}}{{\eta }_{22}{\beta }_{1}1+e^{2{\beta }_1\xi h}{\gamma }_{11}-{\eta }_{21}{\beta }_{2}1+e^{2{\beta }_2\xi h}{\gamma }_{12}{\beta }_{3}1+e^{2{\beta }_3\xi h}{\gamma }_{22}{\gamma }_{11}-{\gamma }_{21}{\gamma }_{12}},\end{array}$$$$\begin{array}{}\text{(B.6)}& b_5=\frac{{\beta }_{1}1+e^{2{\beta }_1\xi h}{\beta }_{2}1+e^{2{\beta }_2\xi h}{\gamma }_{11}}{{\eta }_{22}{\beta }_{1}1+e^{2{\beta }_1\xi h}{\gamma }_{11}-{\eta }_{21}{\beta }_{2}1+e^{2{\beta }_2\xi h}{\gamma }_{12}},\end{array}$$$$\begin{array}{}\text{(B.7)}& b_6=\frac{{\beta }_{1}1+e^{2{\beta }_1\xi h}{\beta }_{2}1+e^{2{\beta }_2\xi h}{\gamma }_{11}{\gamma }_{11}{\gamma }_{23}-{\gamma }_{21}{\gamma }_{13}}{{\eta }_{22}{\beta }_{1}1+e^{2{\beta }_1\xi h}{\gamma }_{11}-{\eta }_{21}{\beta }_{2}1+e^{2{\beta }_2\xi h}{\gamma }_{12}{\gamma }_{11}{\gamma }_{22}-{\gamma }_{21}{\gamma }_{12}},\end{array}$$$$\begin{array}{}\text{(B.8)}& b_7=\frac{{\gamma }_{12}}{{\gamma }_{11}},b_8=\frac{{\gamma }_{13}}{{\gamma }_{11}},b_9=\frac{1}{{\gamma }_{11}},b_{10}=\frac{{\gamma }_{11}}{{\gamma }_{22}{\gamma }_{11}-{\gamma }_{21}{\gamma }_{12}},b_{11}=\frac{1}{{\beta }_{1}1+e^{2{\beta }_1\xi h}},\end{array}$$$$\begin{array}{}\text{(B.9)}& b_{12}=\frac{1}{{\beta }_{2}1+e^{2{\beta }_2\xi h}},b_{13}=\frac{1}{{\beta }_{3}1+e^{2{\beta }_3\xi h}},\end{array}$$$$\begin{array}{}\text{(B.10)}& c_{11}=\frac{b_1{b}_7-b_2{b}_7+b_{9}a-b_8{b}_3}{a}{b}_{11},c_{12}=\frac{b_8{b}_{10}-b_4{b}_7}{a}{b}_{11},c_{13}=\frac{b_6{b}_7-b_5{b}_8}{a}{b}_{11},c_{21}=\frac{b_2-b_1}{a}{b}_{12},c_{22}=\frac{b_4{b}_{12}}{a},c_{23}=-\frac{b_6}{a}{b}_{12},\end{array}$$$$\begin{array}{}\text{(B.11)}& c_{31}=\frac{b_3}{a}{b}_{13},c_{32}=-\frac{b_{10}}{a}{b}_{13},c_{33}=\frac{b_5{b}_{13}}{a}.\end{array}$$

The value of $c_{ji}$ corresponds to the value of $C_{ji}$ for the matrix in Appendix A.

References

[1] V. Y. Topolov, C. R. Bowen, P. Bisegna, Piezo-Active Composites, 2013.

[2] A. Daşdemir, "A finite element model for a bi-layered piezoelectric plate–strip with initial stresses under a time-harmonic force," Journal of the Brazilian Society of Mechanical Sciences and Engineering, vol. 44,DOI: 10.1007/s40430-022-03659-9, 2022.

[3] A. Daşdemir, "A modal analysis of forced vibration of a piezoelectric plate with initial stress by the finite-element simulation," Mechanics of Composite Materials, vol. 58 no. 1, pp. 69-80, DOI: 10.1007/s11029-022-10012-7, 2022.

[4] A. Daşdemir, "Effect of interaction between polarization direction and inclined force on the dynamic stability of a pre-stressed piezoelectric plate," Mechanics Research Communications, vol. 131,DOI: 10.1016/j.mechrescom.2023.104150, 2023.

[5] W. F. Deeg, The analysis of dislocation, crack, and inclusion problems in piezoelectric solids, 1980. Ph. D. Thesis

[6] S. B. Park, C. T. Sun, "Effect of electric field on fracture of piezoelectric ceramics," International Journal of Fracture, vol. 70, pp. 203-216, DOI: 10.1007/BF00012935, 1995.

[7] Y. E. Pak, "Crack extension force in a piezoelectric material," Journal of Applied Mechanics, vol. 57 no. 3, pp. 647-653, DOI: 10.1115/1.2897071, 1990.

[8] M. Nazemizadeh, F. Bakhtiari-Nejad, A. Assadi, B. Shahriari, "Nonlinear vibration of piezoelectric laminated nanobeams at higher modes based on nonlocal piezoelectric theory," Acta Mechanica, vol. 231, pp. 4259-4274, DOI: 10.1007/s00707-020-02736-1, 2020.

[9] S. Hu, J. Liu, J. Li, "Fracture analysis of Griffith interface crack in fine-grained piezoelectric coating/substrate under thermal loading," Advances in Mathematical Physics, vol. 2020,DOI: 10.1155/2020/4201591, 2020.

[10] Y. Zhang, J. Li, X. Xie, "Dynamic analysis of interfacial multiple cracks in piezoelectric thin film/substrate," Acta Mechanica, vol. 234, pp. 705-727, DOI: 10.1007/s00707-022-03390-5, 2023.

[11] Y. E. Pak, "Linear electro-elastic fracture mechanics of piezoelectric materials," International Journal of Fracture, vol. 54 no. 1, pp. 79-100, DOI: 10.1007/BF00040857, 1992.

[12] J. Claeyssen, R. Copetti, L. Tonetto, J. Carvalho, "Nanobeams and AFM subject to piezoelectric and surface scale effects," Advances in Mathematical Physics, vol. 2018,DOI: 10.1155/2018/9268973, 2018.

[13] Z. Wang, X. M. Wang, Y. H. He, Y. M. Hu, H. S. Gu, Y. Wang, "Piezoelectric nanowires in energy harvesting applications," Advances in Materials Science and Engineering, vol. 2015,DOI: 10.1155/2015/165631, 2015.

[14] A. Daşdemir, "Effect of imperfect bonding on the dynamic response of a pre-stressed sandwich plate–strip with elastic layers and a piezoelectric core," Acta Mechanica Solida Sinica, vol. 30 no. 6, pp. 658-667, DOI: 10.1016/j.camss.2017.10.005, 2017.

[15] A. Daşdemir, "Effects of poling direction and imperfection of interfaces on a pre-stressed system with a piezoelectric core bonded to elastic faces," Proceedings of the Institution of Mechanical Engineers Part C: Journal of Mechanical Engineering Science, vol. 236 no. 9, pp. 4671-4688, DOI: 10.1177/09544062211045997, 2022.

[16] K. N. Srivastava, K. Singh, "The effect of penny-shaped crack on the distribution of stress in a semi-infinite solid," International Journal of Engineering Science, vol. 7 no. 5, pp. 469-490, DOI: 10.1016/0020-7225(69)90034-2, 1969.

[17] Y. M. Tsai, "Initiation and propagation of a penny-shaped crack in a finitely deformed incompressible elastic medium," Engineering Fracture Mechanics, vol. 14 no. 3, pp. 627-636, DOI: 10.1016/0013-7944(81)90049-7, 1981.

[18] A. N. Guz, V. M. Nazarenko, "Symmetric failure of the halfspace with penny-shaped cracks in compression," Theoretical and Applied Fracture Mechanics, vol. 3 no. 3, pp. 233-245, DOI: 10.1016/0167-8442(85)90033-3, 1985.

[19] W. Zikun, "Penny-shaped crack in transversely isotropic piezoeleetric materials," Acta Mechanica Sinica, vol. 10, pp. 49-60, DOI: 10.1007/BF02487657, 1994.

[20] H. T. Danyluk, B. M. Singh, J. Vrbik, "A Dugdale-type estimation of the plastic zone for a penny-shaped crack in a thick transversely isotropic layer due to radial shear," Engineering Fracture Mechanics, vol. 51 no. 5, pp. 735-740, DOI: 10.1016/0013-7944(94)00312-6, 1995.

[21] L. Kogan, C.-Y. Hui, V. Molkov, "Stress and induction field of a spheroidal inclusion or a penny-shaped crack in a transversely isotropic piezo-electric material," International Journal of Solids and Structures, vol. 33 no. 19, pp. 2719-2737, DOI: 10.1016/0020-7683(95)00182-4, 1996.

[22] X. Li, W. Chen, H. Wang, G. Wang, "Crack tip plasticity of a penny-shaped Dugdale crack in a power-law graded elastic infinite medium," Engineering Fracture Mechanics, vol. 88,DOI: 10.1016/j.engfracmech.2012.03.006, 2012.

[23] J. H. Huang, "A fracture criterion of a penny-shaped crack in transversely isotropic piezoelectric media," International Journal of Solids and Structures, vol. 34 no. 20, pp. 2631-2644, DOI: 10.1016/S0020-7683(96)00160-6, 1997.

[24] C. Y. Fan, Z. H. Guo, H. Y. Dang, M. H. Zhao, "Extended displacement discontinuity method for nonlinear analysis of penny-shaped cracks in three-dimensional piezoelectric media," Engineering Analysis with Boundary Elements, vol. 38,DOI: 10.1016/j.enganabound.2013.09.014, 2014.

[25] K. Hu, Z. Chen, "Dugdale plastic zone of a penny-shaped crack in a piezoelectric material under axisymmetric loading," Acta Mechanica, vol. 227 no. 3, pp. 899-912, DOI: 10.1007/s00707-015-1501-5, 2016.

[26] Q. Wu, Y. Wan, Z. Zhong, "A penny-shaped crack in a piezoelectric layer sandwiched between two elastic layers with the electrical saturation and mechanical yielding zones," Engineering Fracture Mechanics, vol. 260,DOI: 10.1016/j.engfracmech.2021.108173, 2022.

[27] Y.-B. Zhou, X.-F. Li, "Fracture analysis of an infinite 1D hexagonal piezoelectric quasicrystal plate with a penny-shaped dielectric crack," European Journal of Mechanics-A/Solids, vol. 76, pp. 224-234, DOI: 10.1016/j.euromechsol.2019.04.011, 2019.

[28] B. L. Wang, J. C. Han, S. Y. Du, "Fracture mechanics for multilayers with penny-shaped cracks subjected to dynamic torsional loading," International Journal of Engineering Science, vol. 38 no. 8, pp. 893-901, DOI: 10.1016/S0020-7225(99)00068-3, 2000.

[29] W. J. Feng, E. Pan, X. Wang, "Dynamic fracture analysis of a penny-shaped crack in a magnetoelectroelastic layer," International Journal of Solids and Structures, vol. 44 no. 24, pp. 7955-7974, DOI: 10.1016/j.ijsolstr.2007.05.020, 2007.

[30] H. Y. Yu, K. P. Cooper, "Dynamic penny-shaped cracks in multilayer sandwich composites," Theoretical and Applied Fracture Mechanics, vol. 51 no. 3, pp. 181-188, DOI: 10.1016/j.tafmec.2009.05.003, 2009.

[31] X.-C. Zhong, K.-S. Zhang, "Dynamic analysis of a penny-shaped dielectric crack in a magnetoelectroelastic solid under impacts," European Journal of Mechanics-A/Solids, vol. 29 no. 2, pp. 242-252, DOI: 10.1016/j.euromechsol.2009.10.002, 2010.

[32] Y. Song, H. Hu, J. W. Rudnicki, "Dynamic stress intensity factor (Mode I) of a permeable penny-shaped crack in a fluid-saturated poroelastic solid," International Journal of Solids and Structures, vol. 110-111, pp. 127-136, DOI: 10.1016/j.ijsolstr.2017.01.034, 2017.

[33] B. L. Wang, J. C. Han, "An analytical model for electrode-ceramic interaction in multilayer piezoelectric actuators," Acta Mechanica Sinica, vol. 23, pp. 199-208, DOI: 10.1007/s10409-007-0064-7, 2007.

[34] H. G. Beom, K. M. Jeong, Y. H. Kim, "Intensity factors for subinterface cracks in dissimilar anisotropic piezoelectric media," Archive of Applied Mechanics, vol. 73, pp. 184-198, DOI: 10.1007/s00419-003-0280-y, 2003.

[35] Z. C. Ou, Y. H. Chen, "Near-tip stress fields and intensity factors for an interface crack in metal/piezoelectric bimaterials," International Journal of Engineering Science, vol. 42 no. 13-14, pp. 1407-1438, DOI: 10.1016/j.ijengsci.2004.01.008, 2004.

[36] P. S. Yang, J. Y. Liou, J. C. Sung, "Subinterface crack in an anisotropic piezoelectric bimaterial," International Journal of Solids and Structures, vol. 45 no. 18-19, pp. 4990-5014, DOI: 10.1016/j.ijsolstr.2008.05.001, 2008.

[37] Y. Zhang, J. Li, X. Xie, "Dynamic propagation characteristics of a mode-III interfacial crack in piezoelectric bimaterials," Advances in Materials Science and Engineering, vol. 2022,DOI: 10.1155/2022/1733011, 2022.

[38] C.-F. Gao, M.-Z. Wang, "Collinear permeable cracks between dissimilar piezoelectric materials," International Journal of Solids and Structures, vol. 37 no. 36, pp. 4969-4986, DOI: 10.1016/S0020-7683(99)00184-5, 2000.

[39] S.-H. Chi, Y.-L. Chung, "Cracking in coating–substrate composites with multi-layered and FGM coatings," Engineering Fracture Mechanics, vol. 70 no. 10, pp. 1227-1243, DOI: 10.1016/S0013-7944(02)00114-5, 2003.

[40] X.-F. Li, "Transient response of a piezoelectric material with a semi-infinite mode-III crack under impact loads," International Journal of Fracture, vol. 111 no. 2, pp. 119-130, DOI: 10.1023/A:1012208524059, 2001.

[41] H. F. Tiersten, Linear Piezoelectric Plate Vibrations, 1969.

[42] X.-F. Li, K. Y. Lee, "Effects of electric field on crack growth for a penny-shaped dielectric crack in a piezoelectric layer," Journal of the Mechanics and Physics of Solids, vol. 52 no. 9, pp. 2079-2100, DOI: 10.1016/j.jmps.2004.02.012, 2004.

[43] M. Abramowitz, I. A. Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, 1965.

[44] H. Stehfest, "Algorithm 368: numerical inversion of laplace transforms [D5]," Communications of the ACM, vol. 13 no. 1, pp. 47-49, DOI: 10.1145/361953.361969, 1970.